Numerical ODE Methods: RK4, Adaptive Steps & Stiff Systems

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 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.

1. The Initial Value Problem

We seek to solve:

dy/dt = f(t, y), y(t₀) = y₀, t \in [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 \cdot 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\cdotk₁) k₃ = f(tₙ + h/2, yₙ + h/2\cdotk₂) k₄ = f(tₙ + h, yₙ + h\cdotk₃) yₙ₊₁ = yₙ + h/6\cdot(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:

Energy error is bounded (oscillatory) not drifting
The modified Hamiltonian H̃ = H + O(hᵖ) is exactly conserved
Ideal for N-body gravity, molecular dynamics, and long-term orbital mechanics

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

pₙ₊½ = pₙ - h/2 \cdot \nablaV(qₙ) qₙ₊₁ = qₙ + h \cdot pₙ₊½ / m pₙ₊₁ = pₙ₊½ - h/2 \cdot \nablaV(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.

Method: h: 0.20 γ: 0.10

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 →