Simulation and numerical experimentation

ODE Integration & Euler's Method

When analytical solutions fail, we use numerical methods to approximate physics via step-by-step calculation.

The Euler method is the simplest ODE integrator: $y_{n+1} = y_n + h f(t_n, y_n)$. Here, $y_n$ is the current state, $h$ is the step size, $f(t_n, y_n)$ is the slope, and $y_{n+1}$ is the next state.

Common pitfall: A simulation that runs without errors is not necessarily correct. Always validate against a known analytic case, like exponential decay $e^{-t}$, before trusting new output.

Check convergence by running your simulation with step size $h$ and $h/2$. If results match your desired precision, the solution has converged.

Monte Carlo Methods

Monte Carlo methods use random sampling to estimate numerical quantities when analytical integration is intractable.

To estimate $\pi$, place random points in a unit square. The fraction falling inside the inscribed quarter-circle ($x^2 + y^2 \leq 1$) converges to $\pi/4$:

$$\pi \approx 4 \cdot \frac{\text{points inside}}{\text{total points}}$$

inside = sum(1 for _ in range(N)
             if random.random()**2 
             + random.random()**2 <= 1)
pi_est = 4 * inside / N

Accuracy scales as $O(1/\sqrt{N})$. This means quadrupling points only halves the error, requiring substantial samples for high precision.

Tip: Always validate your simulation against a known analytical case or conservation law before applying it to uncharted problems.