ODE Solvers

Ordinary differential equation solvers compute numerical approximations to y' = f(t, y), y(t₀) = y₀. They are essential for simulating dynamic systems — physics, biology, chemistry, and control.

Euler's Method

Explicit (Forward) Euler

The simplest ODE solver. Approximates the derivative with a forward difference:

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

Derivation: Taylor expansion: y(t+h) = y(t) + hy'(t) + O(h²). Drop the O(h²) term.

Local truncation error: O(h²) per step. Global error: O(h) — first-order method.

Stability: Conditionally stable. For y' = λy (λ < 0), stable iff |1 + hλ| < 1, i.e., h < 2/|λ|.

FUNCTION EULER(f, t0, y0, h, n)
    result ← [(t0, y0)]
    t ← t0
    y ← y0
    FOR i ← 1 TO n DO
        y ← y + h * f(t, y)
        t ← t + h
        APPEND (t, y) TO result
    RETURN result

Implicit (Backward) Euler

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

The right side depends on yₙ₊₁ — requires solving a (possibly nonlinear) equation at each step.

Advantage: Unconditionally stable for y' = λy (λ < 0). Essential for stiff equations.

Cost: More expensive per step (need Newton iteration for nonlinear f).

Runge-Kutta Methods

Higher-order methods that evaluate f at multiple points within each step.

Classical RK4

The workhorse fourth-order method:

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

Local error: O(h⁵). Global error: O(h⁴). Four function evaluations per step.

Stability: Larger stability region than Euler, but still conditionally stable.

FUNCTION RK4(f, t0, y0, h, n)
    result ← [(t0, y0)]
    t ← t0
    y ← y0
    FOR i ← 1 TO n DO
        k1 ← f(t, y)
        k2 ← f(t + h/2, y + h * k1 / 2)
        k3 ← f(t + h/2, y + h * k2 / 2)
        k4 ← f(t + h, y + h * k3)
        y ← y + h/6 * (k1 + 2*k2 + 2*k3 + k4)
        t ← t + h
        APPEND (t, y) TO result
    RETURN result

Adaptive Runge-Kutta

Automatically adjust step size h based on error estimation.

Embedded pairs: Use two RK methods of different orders (computed from the same function evaluations) to estimate error.

Dormand-Prince (DOPRI5/ode45): A 4(5) pair — 4th-order solution with 5th-order error estimate. 6 function evaluations per step (7 with FSAL — First Same As Last).

Error estimate: ‖ŷ₅ - ŷ₄‖
If error < tolerance: accept step, maybe increase h
If error > tolerance: reject step, decrease h

Step size control: h_new = h · (tolerance/error)^(1/(p+1)) · safety_factor.

Fehlberg (RKF45): Another popular 4(5) pair. MATLAB's ode45 originally used this.

Butcher Tableau

A compact notation for RK methods:

c₁ |
c₂ | a₂₁
c₃ | a₃₁  a₃₂
 ⋮ |  ⋮
cₛ | aₛ₁  aₛ₂  ... aₛ,ₛ₋₁
---|------------------------
   | b₁    b₂   ... bₛ        (weights)
   | b̂₁    b̂₂   ... b̂ₛ        (embedded weights, for error estimate)

Explicit RK: aᵢⱼ = 0 for j ≥ i (lower triangular). Implicit RK (DIRK, SDIRK, fully implicit): Non-zero entries above diagonal.

Multi-Step Methods

Use information from multiple previous steps (not just yₙ).

Adams-Bashforth (Explicit)

AB2: yₙ₊₁ = yₙ + h/2 · (3fₙ - fₙ₋₁)
AB3: yₙ₊₁ = yₙ + h/12 · (23fₙ - 16fₙ₋₁ + 5fₙ₋₂)
AB4: yₙ₊₁ = yₙ + h/24 · (55fₙ - 59fₙ₋₁ + 37fₙ₋₂ - 9fₙ₋₃)

Only one new f evaluation per step (efficient!). Need a starting method (e.g., RK) for the first few steps.

Adams-Moulton (Implicit)

AM2: yₙ₊₁ = yₙ + h/12 · (5fₙ₊₁ + 8fₙ - fₙ₋₁)

More stable than Adams-Bashforth. Often used as corrector in predictor-corrector schemes.

Predictor-Corrector

Predict with Adams-Bashforth (explicit).
Evaluate f at the predicted point.
Correct with Adams-Moulton (implicit, using the predicted value — no iteration needed).

BDF (Backward Differentiation Formulas)

BDF1 (backward Euler): yₙ₊₁ = yₙ + h·fₙ₊₁
BDF2: 3yₙ₊₁ - 4yₙ + yₙ₋₁ = 2h·fₙ₊₁

Key property: A-stable for BDF1 and BDF2. The method of choice for stiff equations.

BDF methods up to order 6 are stable. BDF7 and higher are unstable (Dahlquist's barrier).

Stiff Equations

An ODE is stiff if it has widely separated time scales — some components change rapidly, others slowly.

Example: y' = -1000y + 1000cos(t). The homogeneous part decays on time scale 1/1000, but the solution follows cos(t) on time scale 1. Explicit methods must use h < 2/1000 = 0.002, even though the solution is smooth.

Stiffness Ratio

For y' = Ay, the stiffness ratio is |λ_max|/|λ_min| where λᵢ are eigenvalues. Large ratio → stiff.

Implicit Methods for Stiff Systems

Implicit methods (backward Euler, BDF, implicit RK) are A-stable or L-stable — they can take large steps without instability, even for stiff problems.

Cost tradeoff: Each step requires solving a nonlinear system (Newton iteration with Jacobian). But far fewer steps are needed than with explicit methods.

A-stability: The numerical solution decays for all y' = λy with Re(λ) < 0, regardless of step size.

L-stability: A-stable + numerical solution → 0 as h → ∞. Backward Euler is L-stable; trapezoidal rule is A-stable but not L-stable.

Stability Regions

The stability region of a method is the set of hλ ∈ ℂ where the method is stable for y' = λy.

Forward Euler: Circle of radius 1 centered at -1. Small region.
RK4: Larger region, roughly [-2.8, 0] on real axis.
Backward Euler: Everything outside a circle of radius 1 centered at +1. Includes entire left half-plane.
Trapezoidal: Entire left half-plane (A-stable).

Error Control and Adaptive Step-Sizing

Local Error Estimation

Embedded RK methods provide two solutions of different orders. The difference estimates the local error.

Step Size Selection

h_new = h_old · min(max_factor, max(min_factor, safety · (tol/error)^(1/(p+1))))

Typical: safety = 0.9, max_factor = 5, min_factor = 0.2.

Dense Output (Continuous Extension)

Interpolate the solution between mesh points. Enables event detection (zero-crossings), output at arbitrary times, and coupling with other systems.

Symplectic Integrators

For Hamiltonian systems (energy-conserving):

dq/dt = ∂H/∂p,   dp/dt = -∂H/∂q

Standard methods (RK4) gradually drift in energy. Symplectic integrators preserve the symplectic structure, giving much better long-term energy conservation.

Störmer-Verlet (leapfrog):

pₙ₊₁/₂ = pₙ - h/2 · ∂H/∂q(qₙ)
qₙ₊₁ = qₙ + h · ∂H/∂p(pₙ₊₁/₂)
pₙ₊₁ = pₙ₊₁/₂ - h/2 · ∂H/∂q(qₙ₊₁)

Second-order, symplectic, time-reversible. Used in molecular dynamics, celestial mechanics, and HMC (Hamiltonian Monte Carlo).

Summary: Choosing a Method

| Problem Type | Recommended Method | |---|---| | Non-stiff, moderate accuracy | RK4 (fixed step) or DOPRI5 (adaptive) | | Non-stiff, high accuracy | High-order RK (DOP853) or Adams methods | | Stiff | BDF (LSODA, CVODE) or implicit RK (Radau) | | Hamiltonian / long-term | Symplectic integrator (Verlet) | | Very high accuracy | Extrapolation methods (Bulirsch-Stoer) |

Applications in CS

Game physics: Rigid body dynamics. Verlet integration for position-based dynamics. RK4 for smooth motion.
Robotics: Forward dynamics simulation. Trajectory optimization.
Machine learning: Neural ODEs (continuous-depth networks). Adjoint method for gradients through ODE solvers.
Molecular dynamics: Symplectic integrators for atomic simulations (LAMMPS, GROMACS).
Circuit simulation: SPICE uses BDF methods for stiff circuit equations.
Climate modeling: Coupled ODE/PDE systems with multiple time scales.
Population dynamics: SIR epidemic models, predator-prey (Lotka-Volterra).
Control systems: Simulating system response. Real-time control requires fast solvers.
Financial modeling: Stochastic differential equations (SDEs) for option pricing. Euler-Maruyama method.