Virtually every physics simulation boils down to integrating a system of ordinary differential equations. The choice of numerical method determines accuracy, stability, and computational cost — and the wrong choice on a stiff system (one where some quantities change much faster than others, forcing tiny steps to keep the fast part stable) can mean results that are orders of magnitude off or outright diverging. This article surveys the key methods: from the textbook Euler method through the classical RK4, to adaptive Dormand–Prince, multistep Adams–Bashforth, and geometry-preserving symplectic integrators.

TL;DR: This guide compares ODE solvers by accuracy, stability, and cost: Euler is simple but unstable; RK4 is the accurate workhorse for smooth non-stiff problems; Dormand–Prince adapts its step size automatically; Adams–Bashforth reuses past evaluations for cheap high order; implicit BDF methods handle stiff systems; and symplectic integrators like Verlet conserve energy exactly in long-run orbital and molecular simulations.

1. The Initial Value Problem

We seek to solve:

dy/dt = f(t, y), y(t₀) = y₀, t ∈ [t₀, T]

where y ∈ ℝⁿ and f: ℝ × ℝⁿ → ℝⁿ. All classical physical systems (Newton's second law, Lorenz equations, chemical kinetics) fit this form after reducing higher-order equations to first-order systems.

A numerical integrator approximates y(t) at a sequence of discrete time points tₖ = t₀ + k·h (fixed step) or at adaptively chosen tₖ. The core challenge is balancing three competing demands: accuracy (local truncation error per step), stability (error doesn't grow unboundedly), and cost (function evaluations per step).

2. Euler Method (Order 1)

yₙ₊₁ = yₙ + h · f(tₙ, yₙ)

Local truncation error (LTE): O(h²). Global error: O(h). The simplest method and almost never used in practice — it requires extremely small h for acceptable accuracy. Explicit Euler is also A-unstable: for the test equation dy/dt = λy with Re(λ) < 0, the method requires h < 2/|λ| to avoid growing oscillations, which is very restrictive for stiff problems.

Symplectic Euler — a simple improvement

For Hamiltonian systems (q position, p momentum), the symplectic Euler method updates p first, then uses the new p to update q:
pₙ₊₁ = pₙ + h·f(qₙ)
qₙ₊₁ = qₙ + h·pₙ₊₁/m
Despite being only first-order accurate, it preserves the symplectic structure (area in phase space), so energy error remains bounded over exponentially long times — unlike the Runge-Kutta methods which slowly accumulate energy drift.

3. Classical Runge–Kutta (RK4, Order 4)

The "gold standard" workhorse of physics simulation:

k₁ = f(tₙ, yₙ) k₂ = f(tₙ + h/2, yₙ + h/2·k₁) k₃ = f(tₙ + h/2, yₙ + h/2·k₂) k₄ = f(tₙ + h, yₙ + h·k₃) yₙ₊₁ = yₙ + h/6·(k₁ + 2k₂ + 2k₃ + k₄)

LTE: O(h⁵). Global error: O(h⁴). Four function evaluations per step. The error constant is very small: for the harmonic oscillator the error per unit time is roughly 1/180 · h⁴ · max|y⁽⁵⁾|. For most non-stiff problems with smooth solutions, RK4 with h ≈ 0.01–0.1 is highly accurate. Stability region in the complex h·λ plane is a disc-like region covering most of the left half-plane up to |hλ| ≈ 2.8. 

// RK4 in JavaScript
function rk4(f, t, y, h) {
  const k1 = f(t,       y);
  const k2 = f(t + h/2, y.map((v, i) => v + h/2 * k1[i]));
  const k3 = f(t + h/2, y.map((v, i) => v + h/2 * k2[i]));
  const k4 = f(t + h,   y.map((v, i) => v + h   * k3[i]));
  return y.map((v, i) => v + h/6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]));
}

4. Adaptive Step-Size Control: Dormand–Prince RK45

Fixed-step RK4 wastes function evaluations in smooth regions and steps too coarsely in rapidly varying ones. Adaptive methods control the step size to keep estimated local error within a tolerance (rtol, atol).

The Dormand–Prince method (Matlab's ode45, SciPy's RK45) uses a pair of embedded RK formulas sharing the same 6 stages (+ one FSAL "first same as last" reuse): one 5th-order and one 4th-order. Their difference estimates the local error without extra cost.

err ≈ yₙ₊₁⁽⁵⁾ − yₙ₊₁⁽⁴⁾ (error estimate) h_new = h · (tol / err)^{1/5} · safety_factor (safety = 0.9)

The FSAL trick reuses k₇ from step n as k₁ of step n+1, reducing cost from 6 to 5 evaluations per accepted step. Most practical simulations use rtol = 1e-6, atol = 1e-9 as defaults for scientific accuracy.

Method Order Evals/step Adaptive? Best use
Euler 1 1 No Teaching only
Heun (RK2) 2 2 No Very coarse approximations
RK4 (classical) 4 4 No Non-stiff, known smooth ODE
Dormand–Prince RK45 4(5) 6 (5 FSAL) Yes General non-stiff default
Cash–Karp RK45 4(5) 6 Yes Smooth non-stiff with tight tol
DOP853 8(5,3) 13 Yes High-precision non-stiff
Verlet / leapfrog 2 1 No Hamiltonian, N-body, MD
Störmer–Verlet 2 1 No Symplectic, energy conservation

5. Multistep Methods: Adams–Bashforth & Adams–Moulton

Single-step methods (Runge–Kutta) discard all previously computed function values after each step. Multistep methods reuse the last k function evaluations, achieving high order with just one or two new evaluations per step.

The explicit Adams–Bashforth p-step formula (AB-p):

yₙ₊₁ = yₙ + h · Σⱼ₌₀^{p−1} βⱼ · f(tₙ₋ⱼ, yₙ₋ⱼ)

AB-4 (4-step) is order 4 with only 1 new evaluation per step (vs 4 for RK4). Implicit Adams–Moulton (AM-p) achieves order p+1 with p steps. The combination explicit predictor + implicit corrector is the PECE (Predict–Evaluate–Correct–Evaluate) scheme used in Matlab's ode113.

Limitation: multistep methods require a startup procedure (use RK4 for the first k steps) and cannot easily change step size (require restarting or interpolation).

6. Stiff Systems and Implicit Methods

A system is stiff when the largest eigenvalue of the Jacobian ∂f/∂y has magnitude much larger than the inverse of the time scale of interest: stiffness ratio S = max|Re(λᵢ)| · T ≫ 1.

Examples: chemical kinetics with fast/slow reactions (S can reach 10¹⁵), RC circuit with widely separated time constants, diffusion equations after spatial discretisation (eigenvalues ∝ −1/h²).

A-stability and L-stability

An integrator is A-stable if for all λ with Re(λ)≤0, |R(hλ)| ≤ 1 for all h>0 — the amplification factor never exceeds 1. Explicit RK methods are not A-stable (their stability regions in the left half-plane are bounded). Implicit methods can be A-stable.

L-stability (stronger): additionally R(∞)=0, so stiff components decay immediately. Backward-Euler (order 1), Trapezoidal (order 2, A-stable but not L-stable), SDIRK methods, and Rosenbrock methods provide A/L-stable integration.

For stiff systems, the implicit Backward Euler scheme requires solving a nonlinear system at each step (Newton iteration), but allows arbitrarily large h without instability. The BDF-k (Backward Differentiation Formula) methods (Gear methods), used in MATLAB's ode15s, achieve order up to 6 while remaining stiffly stable.

7. Symplectic Integrators for Hamiltonian Systems

For conservative mechanical systems with Hamiltonian H = T(p) + V(q), Runge–Kutta methods introduce a slow secular energy drift even at high order. A symplectic integrator preserves the symplectic 2-form ω = Σ dqᵢ ∧ dpᵢ exactly, which in practice means:

The Verlet/leapfrog integrator is 2nd-order symplectic:

pₙ₊½ = pₙ − h/2 · ∇V(qₙ) qₙ₊₁ = qₙ + h · pₙ₊½ / m pₙ₊₁ = pₙ₊½ − h/2 · ∇V(qₙ₊₁)

Higher-order symplectic methods (Yoshida 4th-order, Ruth–Forest 4th-order) apply the Verlet map with cleverly scaled sub-steps. The price: no adaptive step size — the step must stay fixed for symplecticity to hold.

8. Interactive: Compare Methods on the Damped Oscillator

The damped harmonic oscillator ÿ + 2γẏ + ω₀²y = 0 has exact solution y(t) = e^{−γt}·cos(ωd·t) where ωd = √(ω₀²−γ²). Adjust h (step size) and γ (damping) to see how each method tracks the true solution. Large h reveals how Euler diverges while RK4 stays accurate. The absolute energy error over time is shown in the lower panel.

Exact   Numerical   Energy error (scaled)
Verlet with γ=0 conserves exact energy; all explicit RK methods accumulate energy drift proportional to h^order.

← Fractal Geometry Linear Algebra Simulations →