Simulation and numerical experimentation

When analytical solutions are impossible, we use numerical methods to approximate the physics.

Euler method — The simplest ODE integrator:

$$y_{n+1} = y_n + h\,f(t_n,\,y_n)$$

def euler(f, y0, t0, t_end, h):
    t, y = t0, y0
    trajectory = [(t, y)]
    while t < t_end:
        y = y + h * f(t, y)
        t = t + h
        trajectory.append((t, y))
    return trajectory

Better integrators

• Midpoint (RK2): Evaluate slope at the midpoint for better accuracy.
• RK4 (Runge-Kutta 4): Four slope evaluations per step — the workhorse of physics simulations.
• Error scales as $O(h^{p})$ where $p$ is the method's order (Euler: $p=1$, RK4: $p=4$).

Monte Carlo methods — Use random sampling to estimate quantities:

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

Convergence check — Run the simulation with step size $h$ and $h/2$. If results agree to the desired precision, the solution is converged.

Tip: Always validate your simulation against a known analytical case before applying it to a new problem. If Euler's method can't match $e^{-t}$ for exponential decay, it won't be reliable for anything more complex.

Common pitfall: A simulation that runs without errors is not a simulation that is right. Validate against a case with a known answer (analytic solution, conservation law) before trusting any new output.

Monte Carlo Simulation

Monte Carlo methods use random sampling to estimate numerical quantities.

To estimate $\pi$: place random points in a unit square. The fraction that falls 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}}$$

Accuracy improves as $O(1/\sqrt{N})$ — quadrupling points halves the error.

Monte Carlo methods are essential when analytical integration is intractable.