📐 Euler vs Runge-Kutta: Which Integration Method for Your Simulation?

The Integration Problem

Most physics simulations are governed by ordinary differential equations (ODEs). A pendulum's angle θ satisfies d²θ/dt² = −(g/L)·sin(θ). A gravitational n-body system satisfies Newton's second law for each body. A fluid is governed by the Navier-Stokes equations. To simulate these systems, we must integrate the ODEs forward in time.

The core challenge: computers can only represent discrete states. We pick a timestep Δt and ask "given the state at time t, what is the state at time t + Δt?" The answer is always an approximation — and the quality of that approximation depends on the integration method. Poor choices lead to energy drift, instability, or outright numerical explosion.

Euler's Method: Simple but Dangerous

Leonhard Euler's method, dating to 1768, is the most straightforward approach. For a first-order ODE dy/dt = f(y, t), it steps forward as:

y(t + Δt) ≈ y(t) + Δt · f(y(t), t)

This is a first-order Taylor expansion — it uses only the derivative at the current point to extrapolate. The local truncation error (error per step) is proportional to Δt², and the global error accumulated over a fixed time interval is O(Δt) — hence "first-order method."

The practical problem is energy drift. For conservative systems like pendulums or orbits, Euler's method systematically adds energy to the system. A pendulum simulated with Euler integration will slowly spiral outward, gaining amplitude indefinitely. Reduce Δt and the drift slows — but never stops. This is a fundamental property of the method, not a bug to be fixed by smaller timesteps.

Runge-Kutta 4: The Workhorse

The fourth-order Runge-Kutta method (RK4) evaluates the derivative at four points within each timestep and combines them with a weighted average:

k1 = f(y,       t)
k2 = f(y + k1·Δt/2,  t + Δt/2)
k3 = f(y + k2·Δt/2,  t + Δt/2)
k4 = f(y + k3·Δt,    t + Δt)

y(t+Δt) = y(t) + (Δt/6)·(k1 + 2k2 + 2k3 + k4)

The local truncation error is O(Δt⁵) and the global error is O(Δt⁴). Doubling the timestep multiplies the error by 16. Halving it divides the error by 16. This dramatic sensitivity to step size is why RK4 is so powerful — modest timesteps achieve high accuracy.

RK4 requires four function evaluations per step (versus one for Euler), so it costs about 4× more per step. But because it achieves accuracy with larger timesteps, the total cost to reach a given accuracy level is dramatically lower. For most smooth physics problems, RK4 is the default choice.

The Energy Drift Problem: Symplectic Integrators

RK4 is dramatically better than Euler, but it still has a fundamental flaw for conservative mechanical systems: it is not symplectic. A symplectic integrator preserves the geometric structure of Hamiltonian mechanics — specifically, it conserves the phase-space volume (Liouville's theorem). Non-symplectic methods violate this and produce secular energy drift.

The simplest symplectic integrator is the Leapfrog (or Stormer-Verlet) method. For a system with position x and velocity v:

// Half-step velocity update
v_half = v(t) + (Δt/2) · a(x(t))

// Full-step position update
x(t+Δt) = x(t) + Δt · v_half

// Compute acceleration at new position
a_new = f(x(t+Δt))

// Complete velocity step
v(t+Δt) = v_half + (Δt/2) · a_new

Leapfrog is only second-order accurate — worse than RK4 per step. But because it is symplectic, it conserves energy on average over long timescales. The energy oscillates around the true value rather than drifting monotonically. For long-running simulations of conservative systems — planetary orbits, molecular dynamics, particle accelerators — symplectic integrators are essential.

Higher-Order Symplectic Methods

Forest and Ruth (1990) and Yoshida (1990) showed how to construct higher-order symplectic integrators by combining Leapfrog steps with carefully chosen coefficients. A fourth-order symplectic integrator requires roughly 4–7 Leapfrog sub-steps per full step, achieving RK4-level accuracy while preserving symplectic structure. These methods are standard in long-duration molecular dynamics and N-body astronomy codes.

The trade-off matrix:

• Euler: Fast, first-order, energy-drifting. Use only for quick prototypes or non-conservative systems.
• RK4: Four evaluations per step, fourth-order, slightly energy-drifting. Best for accurate short-to-medium simulations of smooth systems.
• Leapfrog/Verlet: Two evaluations per step, second-order, energy-preserving on average. Best for long conservative-system simulations.
• 4th-order symplectic: 5–7 evaluations per step, fourth-order, symplectic. Best of both worlds for long, demanding simulations.

Adaptive Step Size

Fixed-timestep methods struggle when the solution changes rapidly in some regions and slowly in others. Adaptive methods estimate the local error by comparing solutions of different orders (e.g., RK4 vs RK5) and adjust Δt dynamically. The Dormand-Prince method (used in MATLAB's ode45 and SciPy's solve_ivp) achieves this with just six function evaluations. Step size is increased when error is small (to gain speed) and decreased when error is large (to maintain accuracy).