ODE Solvers Compared: Euler, RK4, Verlet, Adaptive Methods

1. The ODE Problem & Error Concepts

An ordinary differential equation (ODE) expresses the time derivative of a state vector y as a function of that same state:

dy/dt = f(t, y) y(t₀) = y₀ Examples: Pendulum: y = [θ, ω], dy/dt = [ω, -(g/L)sin θ] Orbit: y = [x, y, vx, vy], dy/dt = [vx, vy, -GMx/r³, -GMy/r³] Logistic: y = N, dy/dt = rN(1 - N/K) A numerical integrator advances y from t to t+Δt using only f evaluations — it never sees a closed-form solution.

Two types of error matter:

Local truncation error (LTE): the error made in a single step, ignoring prior error. A method of order p has LTE ∝ Δt^(p+1).
Global error: the accumulated error after integrating from t₀ to T. For a p-th order method, global error ∝ Δt^p (one power less than LTE, because there are T/Δt steps).

A 2nd-order global error method halves its error when Δt is halved — a 4th-order method reduces error by 16× — meaning higher-order methods earn back their extra computation cost very quickly when accuracy is needed.

2. Forward Euler: Simple but Dangerous

The simplest possible integrator: take the derivative at the current point and step in that direction.

// Forward (explicit) Euler y(t+Δt) = y(t) + Δt \cdot f(t, y(t)) Cost: 1 function evaluation per step Order: 1 (global error \propto Δt) Example — harmonic oscillator: y'' = -y \to y = [x, v] x_new = x + Δt \cdot v v_new = v + Δt \cdot (-x) Exact solution: x(t) = A cos(t), energy = constant Euler result: energy grows without bound for any Δt > 0 (orbit spirals outward — integrator is unstable for oscillations)

Forward Euler is unconditionally energy-increasing for conservative oscillations — the numerical orbit gains energy at every step, no matter how small Δt. With Δt = 0.01 for a unit harmonic oscillator, amplitude grows by ~0.5% per orbit period. After 1000 periods the amplitude has doubled. This makes it unsuitable for any simulation requiring long-term energy conservation.

When Forward Euler is fine: dissipative systems (where energy should decrease) and transient simulations where only short-time accuracy matters — e.g., one step of a game physics update, or a rough first-pass trajectory. It is also the foundation for understanding all higher-order methods.

Stability Region

For the test equation dy/dt = λy, Euler is stable when |1 + λΔt| ≤ 1. For purely imaginary λ (oscillations), this circle never contains the imaginary axis — confirming that Euler is always unstable for conservative oscillators regardless of timestep size.

3. Symplectic Euler & Why It Matters

A one-character change from Forward Euler produces a dramatically better integrator for Hamiltonian (energy-conserving) systems: use the updated velocity when computing the new position.

// Symplectic (Semi-implicit) Euler v_new = v + Δt \cdot a(x) // update velocity first x_new = x + Δt \cdot v_new // use NEW velocity for position vs. Forward Euler: v_new = v + Δt \cdot a(x) // same x_new = x + Δt \cdot v // uses OLD velocity Cost: 1 force evaluation per step (same as Forward Euler) Order: 1 (same local accuracy) BUT: exactly preserves a modified Hamiltonian \to no energy drift!

The key difference is geometric: Symplectic Euler preserves the symplectic structure of phase space — areas of phase-space regions are conserved under the map. This is a discrete version of Liouville's theorem. The consequence is that a symplectic integrator's numerical orbit lies exactly on a slightly perturbed energy surface, oscillating around the true energy rather than drifting away from it.

Energy comparison over 10 000 orbital periods: Forward Euler (Δt=0.01): energy drifts +47% (orbit escapes) Symplectic Euler (Δt=0.01): energy oscillates \pm0.01% (bounded) RK4 (Δt=0.01): energy oscillates \pm10⁻⁷ (very small drift) For a game physics engine running at 60 fps (Δt\approx0.016 s), Symplectic Euler keeps objects from spontaneously gaining energy — Forward Euler would cause visible "popping" in spring systems.

This is why virtually every game physics engine (Box2D, Bullet, Unity's PhysX integration) uses Symplectic Euler, not Forward Euler. It costs the same but behaves fundamentally better.

4. Runge-Kutta 4: The Workhorse

RK4 evaluates the derivative at four carefully chosen points within each timestep and combines them in a weighted average that cancels error terms through 4th order.

// Classical 4th-order Runge-Kutta (RK4) k₁ = f(t, y) k₂ = f(t + Δt/2, y + Δt/2 \cdot k₁) k₃ = f(t + Δt/2, y + Δt/2 \cdot k₂) k₄ = f(t + Δt, y + Δt \cdot k₃) y(t+Δt) = y(t) + (Δt/6)(k₁ + 2k₂ + 2k₃ + k₄) Cost: 4 function evaluations per step Order: 4 (global error \propto Δt⁴) Error comparison vs Forward Euler at same Δt: Euler (order 1): error \propto Δt RK4 (order 4): error \propto Δt⁴ \to ~10⁴\times smaller at Δt=0.1

The coefficients (1/6, 1/3, 1/3, 1/6) match Simpson's rule — RK4 is essentially applying Simpson's rule to the integral form of the ODE. The two midpoint evaluations k₂ and k₃ refine each other, cancelling the error term that a single midpoint evaluation (the 2nd-order midpoint method) would leave.

Practical Accuracy

RK4 with Δt = 0.1 is typically more accurate than Forward Euler with Δt = 0.001 — using the same total number of function evaluations. This is the core argument for higher-order methods: they buy accuracy far more cheaply than simply shrinking the timestep.

RK4 Drawbacks

Not symplectic: For long-time Hamiltonian dynamics, RK4 still exhibits slow energy drift (though much slower than Forward Euler). Leapfrog with the same cost outperforms RK4 on energy conservation over millions of steps.
Four evaluations: If f is expensive (e.g., gravity from N bodies — O(N²) per evaluation), the 4× cost matters.
Fixed step size: Classic RK4 doesn't adapt Δt. Stiff problems require very small Δt throughout even if most of the solution varies slowly.

5. Verlet & Leapfrog: Long-Time Accuracy

The Verlet integrator was developed in 1967 for molecular dynamics simulation — contexts where trillions of timesteps must be computed while preserving total energy. It exploits the time-reversibility and symplectic structure of Newtonian mechanics.

// Störmer-Verlet (original form) x(t+Δt) = 2\cdotx(t) - x(t-Δt) + Δt² \cdot a(t) Requires TWO previous positions — awkward to start and to track velocities. // Velocity Verlet (equivalent, standard form) x(t+Δt) = x(t) + v(t)\cdotΔt + ½\cdota(t)\cdotΔt² a(t+Δt) = f(x(t+Δt)) / m // recompute force at new position v(t+Δt) = v(t) + ½ \cdot (a(t) + a(t+Δt)) \cdot Δt Cost: 1 force evaluation per step (a(t+Δt) is reused next step) Order: 2 for positions and velocities BUT: symplectic \to energy error is bounded, not drifting

Leapfrog Form

Leapfrog is algebraically equivalent to Velocity Verlet but stores velocities half a timestep offset ("leaping over" each other):

// Leapfrog (half-step velocities) v(t + Δt/2) = v(t - Δt/2) + Δt \cdot a(t) x(t + Δt) = x(t) + Δt \cdot v(t + Δt/2) Positions and velocities alternate: x at integer steps, v at half steps. Exactly time-reversible: run backwards \to exactly retraces the orbit. Symplectic: phase-space volume preserved exactly. Energy error oscillates around 0 — bounded by O(Δt²) — forever. No secular drift even after 10⁹ steps.

Molecular dynamics benchmark: Integrating a system of 1000 Lennard-Jones particles for 10⁶ steps: Leapfrog RMS energy error = 0.003%. RK4 (same cost) = 0.0002% but growing — after 10⁸ steps RK4 error exceeds Leapfrog. For astronomy (N-body orbital dynamics), leapfrog is the standard for the same reason.

Higher-Order Symplectic Methods

Forest-Ruth (4th order) and Yoshida (6th, 8th order) methods compose multiple leapfrog steps with specific positive and negative sub-steps to cancel higher-order error while preserving symplecticity. They achieve the accuracy of RK4/RK8 with the long-time energy bounds of leapfrog — at the cost of more sub-steps per full step (6–15 evaluations).

6. Adaptive Step-Size Control

Fixed-step integrators waste computation during smooth parts of a solution and may miss rapid changes. Adaptive methods continuously estimate the local error and adjust Δt to maintain a user-specified tolerance:

// Embedded Runge-Kutta pair (e.g. Dormand-Prince RK45) Two solutions computed using the SAME k evaluations: y4 — 4th-order solution y5 — 5th-order solution (1 extra evaluation) Error estimate: err = |y5 - y4| (difference between orders) Step size controller (PI controller): If err > tol: reject step, reduce Δt: Δt_new = Δt \cdot (tol/err)^(1/5) If err < tol: accept step, possibly grow Δt: Δt_new = Δt \cdot (tol/err)^(1/5) \cdot safety

The Dormand-Prince RK45 pair (the basis of MATLAB's ode45 and SciPy's RK45) uses 6 evaluations per step with 5th-order accuracy. The 4th-order estimate comes "for free" from the same k values — the clever choice of Butcher tableau coefficients ensures this.

Cash-Karp and Bogacki-Shampine

Bogacki-Shampine (RK23) is a cheaper 3rd/2nd-order pair using only 4 evaluations — ideal for real-time applications or when the solution has limited smoothness. Cash-Karp (RK45) was the previous MATLAB standard; Dormand-Prince is slightly more accurate for the same cost and has become the modern default.

Tolerance setting: Adaptive solvers take a relative tolerance rtol and absolute tolerance atol. Use rtol = 1e-6, atol = 1e-9 for accurate scientific work. rtol = 1e-3 for interactive real-time use. Very tight tolerances (< 1e-12) are usually not meaningful because floating-point round-off dominates.

7. Stiff ODEs & Implicit Methods

A stiff ODE has components with vastly different timescales. The classic example: a chemical reaction network where some species react in microseconds while others vary over hours. An explicit solver must use the smallest relevant Δt (microseconds) throughout the full simulation (hours) — 10¹⁰ steps. An implicit solver can handle the fast component analytically and take large steps for the slow dynamics.

// Implicit (Backward) Euler — simplest implicit method y(t+Δt) = y(t) + Δt \cdot f(t+Δt, y(t+Δt)) // uses FUTURE state This is a nonlinear equation in y(t+Δt) — must solve iteratively (Newton's method per step). Each iteration requires a Jacobian \partialf/\partialy. Stability: A-stable — stable for ALL λΔt with Re(λ) < 0. Explicit RK4 is only stable for |λΔt| in a bounded region. Cost per step: much higher than explicit (Newton iterations) But: Δt can be 10³-10⁶\times larger on stiff problems \to net WIN

BDF and Rosenbrock Methods

Gear's Backward Differentiation Formula (BDF) methods — used in MATLAB's ode15s — are the industry standard for stiff chemical kinetics, circuit simulation, and structural dynamics with widely varying spring stiffnesses. LSODA (used in SciPy) automatically switches between non-stiff (Adams-Bashforth) and stiff (BDF) integrators based on detected stiffness.

Rosenbrock methods linearise the implicit equation analytically rather than using Newton iteration — cheaper per step but limited to moderate stiffness ratios.

Is your problem stiff? Compute the Jacobian eigenvalues λᵢ. If max|Re(λᵢ)| / min|Re(λᵢ)| exceeds ~10³, it is stiff. Practical indicators: explicit solvers require tiny Δt to stay stable even while the solution appears smooth; the integrator keeps rejecting and shrinking steps on what looks like a boring part of the trajectory.

8. Which Solver Should You Use?

Use Case	Recommended Solver	Reason
Game physics (rigid bodies, springs)	Symplectic Euler	1 eval, no energy drift, cheap, ~60 fps
Orbital mechanics, N-body	Leapfrog / Velocity Verlet	Symplectic, 1 eval, millions of steps stable
Molecular dynamics	Velocity Verlet	Industry standard, time-reversible, NVE ensemble
Smooth trajectory, moderate accuracy	RK4	4th order, simple to implement, widely understood
High accuracy, unknown smoothness	Dormand-Prince RK45	Adaptive, automatic error control, `ode45`
Stiff chemistry / circuits	BDF (ode15s / LSODA)	A-stable, handles 10⁶:1 timescale ratio
Long-time Hamiltonian (high order)	Yoshida 6th order	Symplectic + high accuracy, no secular drift
Quick prototype / teaching	Forward Euler	5 lines of code, visualise errors directly

Summary: The Three Key Properties

Order (accuracy per step): Higher order → fewer steps needed → less total work. RK4 (order 4) beats Euler (order 1) by 10⁴× at Δt = 0.1.
Symplecticity (long-time energy conservation): Leapfrog and Symplectic Euler preserve a modified Hamiltonian. Non-symplectic methods (including RK4) drift on conserved quantities over very long integrations.
A-stability (stiffness handling): Implicit methods are A-stable — they do not have a step-size restriction from stability. Explicit methods require Δt < O(1/|λ_max|) regardless of how slowly the interesting dynamics actually evolve.

No single integrator wins on all three. The best choice depends on the physics (conservative vs dissipative vs stiff), the required accuracy, and the cost of one force/derivative evaluation.