Compendium

ODE Numerical Methods — Math Reference

22 min read

ODE Numerical Methods — Math Reference

Ordinary differential equations (ODEs) relate a function and its derivatives with respect to a single independent variable (usually time). They are the lingua franca of dynamical modelling: every continuous-time system from a falling apple to a feedforward residual network can be expressed as an ODE. This reference covers analytical techniques for the small set of ODEs that admit closed-form solutions, then the numerical machinery — explicit and implicit Runge-Kutta, multistep methods, adaptive step control, stiff solvers, boundary-value problems, delay and stochastic extensions, and symplectic integrators — that handles the rest, along with the production software ecosystem and engineering / ML applications.

1. At a Glance

What ODEs model

DomainEquation formVariables
Mechanical motionm·ẍ + c·ẋ + k·x = F(t)position x(t)
Thermal transientsρ c_p · dT/dt = q̇ − h(T − T_∞)temperature T(t)
Chemical kineticsdC_i/dt = Σ ν_{ij} r_j(C)concentration C_i(t)
Electrical circuitsL · dI/dt + R·I + (1/C)·∫I dt = V(t)current I(t)
Population dynamicsdN/dt = r·N·(1 − N/K) (logistic)population N(t)
FinancedS = µ·S dt + σ·S dW (geometric Brownian)asset price S(t)
Neural ODEsdh/dt = f(h, t; θ) (Chen 2018)hidden state h(t)
Epidemics (SIR)Ṡ = −β SI/N, İ = β SI/N − γI, Ṙ = γIcompartment populations
Orbital mechanicsr̈ = −µ r / ‖r‖³ (Kepler)position vector r(t)

Classification

  • Linear vs nonlinear — linear if both y and its derivatives appear to the first power and not as arguments of nonlinear functions. The principle of superposition holds only for linear ODEs.
  • Autonomous vs non-autonomous — autonomous when f(t, y) does not depend explicitly on t (e.g. ẏ = −y vs ẏ = sin(t)·y).
  • Single equation vs system — n-th order ODEs are converted to a system of n first-order ODEs by stacking derivatives.
  • Initial-value problem (IVP) vs boundary-value problem (BVP) — IVP specifies all state at one t = t₀; BVP specifies values at two or more points (typical for spatial problems on a finite domain).
  • Stiff vs non-stiff — stiff when widely separated time scales coexist; explicit methods are forced into impractically small steps to maintain stability.
  • Ordinary vs delay vs stochastic — DDE y’(t) = f(t, y(t), y(t − τ)) introduces history; SDE adds a Wiener increment.

The same physical system can sit in different categories depending on which question is asked: simulating a rocket trajectory is an IVP, finding a trajectory that hits a target is a BVP, modelling sensor noise on it is a SDE.

2. Analytical Methods for Simple ODEs

Closed-form solutions exist for a small but important set.

2.1 Separable equations

If dy/dx = f(x) g(y), separate and integrate:

∫ dy/g(y) = ∫ f(x) dx + C

Example. Newton’s cooling dT/dt = −k(T − T_∞). With u = T − T_∞, du/u = −k dt, ln|u| = −k t + C, T(t) = T_∞ + (T₀ − T_∞) e^{−k t}.

2.2 First-order linear

dy/dx + p(x) y = q(x).

Multiply by integrating factor µ(x) = exp(∫ p(x) dx):

d/dx [µ y] = µ q ⇒ y(x) = (1/µ) [∫ µ q dx + C].

Example. RC circuit RC · dV/dt + V = V_s. Here p = 1/(RC), µ = e^{t/RC}, V(t) = V_s + (V₀ − V_s) e^{−t/RC}.

2.3 Exact equations

M(x, y) dx + N(x, y) dy = 0 is exact when ∂M/∂y = ∂N/∂x. A potential F with F_x = M, F_y = N exists; the solution is F(x, y) = C. Non-exact equations can sometimes be made exact with an integrating factor µ(x) or µ(y).

2.4 Linear homogeneous with constant coefficients

a_n y^{(n)} + a_{n−1} y^{(n−1)} + … + a₁ y’ + a₀ y = 0.

Substitute y = e^{λ t} → characteristic polynomial p(λ) = a_n λ^n + … + a₀ = 0. For each root λ_k:

  • Real distinct: contributes C_k e^{λ_k t}.
  • Repeated of multiplicity m: contributes (C₀ + C₁ t + … + C_{m−1} t^{m−1}) e^{λ t}.
  • Complex conjugate α ± iβ: contributes e^{α t}(A cos β t + B sin β t).

For a system ẋ = A x, eigenvalues of A play the role of λ; eigenvectors give the mode shapes. This is the bridge between ODE theory and linear algebra (see [[Math/linear-algebra-essentials]]).

2.5 Second-order linear — the canonical resonator

y” + 2 ζ ω_n y’ + ω_n² y = f(t)

with ω_n the undamped natural frequency and ζ the damping ratio. The free response (f = 0) splits into three regimes:

RegimeζRootsResponse
Overdampedζ > 1both real negativesum of two decaying exponentials
Critically dampedζ = 1double real −ω_n(A + B t) e^{−ω_n t}
Underdamped0 < ζ < 1−ζ ω_n ± i ω_d, ω_d = ω_n √(1 − ζ²)e^{−ζ ω_n t}(A cos ω_d t + B sin ω_d t)

This single equation captures mass-spring-damper systems, RLC circuits, pendulums under small-angle approximation, and the dynamics of second-order PID-controlled plants.

2.6 Laplace transform method

For linear ODEs with constant coefficients and right-hand side f(t):

L{y^{(n)}} = s^n Y(s) − s^{n−1} y(0) − … − y^{(n−1)}(0)

Solve algebraically in the s-domain, then invert. Particularly natural for piecewise-defined forcing (step, impulse, ramp), feedback control transfer functions, and convolution. See [[Engineering/classical-control]] for control-engineering usage.

3. Initial-Value Problem Numerical Methods

For ẏ = f(t, y), y(t₀) = y₀, step from t_n to t_{n+1} = t_n + h.

3.1 Forward (explicit) Euler

y_{n+1} = y_n + h · f(t_n, y_n)

Local truncation error O(h²), global O(h). Stable only when |1 + h λ| ≤ 1 — a disk of radius 1 in the complex plane. Useful as a building block; rarely the right choice in production.

Worked example. ẏ = −2 y, y(0) = 1, h = 0.1 (true solution y(t) = e^{−2t}):

nt_ny_nTrue y(t_n)Error
00.01.0001.0000.000
10.10.8000.8190.019
20.20.6400.6700.030
50.50.3280.3680.040

Step y_1 = 1.0 + 0.1 · (−2 · 1.0) = 0.8.

3.2 Backward (implicit) Euler

y_{n+1} = y_n + h · f(t_{n+1}, y_{n+1})

Requires solving (typically via Newton iteration on the residual r(y_{n+1}) = y_{n+1} − y_n − h f(t_{n+1}, y_{n+1})). A-stable — entire left half plane is in its stability region — so good for stiff problems.

3.3 Trapezoidal (Crank-Nicolson for ODEs)

y_{n+1} = y_n + [h/2](f(t_n, y_n) + f(t_{n+1}, y_{n+1}))

Second order, A-stable, but not L-stable: stiff transients ring without damping.

3.4 RK2 — midpoint and Heun

Midpoint:

k₁ = f(t_n, y_n) y_{n+1} = y_n + h · f(t_n + h/2, y_n + (h/2) k₁)

Heun (improved Euler):

k₁ = f(t_n, y_n) k₂ = f(t_n + h, y_n + h k₁) y_{n+1} = y_n + (h/2)(k₁ + k₂)

Both 2nd-order; explicit; useful when accuracy demands of RK4 are unnecessary.

3.5 Classical RK4

The most widely taught method (Kutta 1901):

k₁ = f(t_n, y_n) k₂ = f(t_n + h/2, y_n + (h/2) k₁) k₃ = f(t_n + h/2, y_n + (h/2) k₂) k₄ = f(t_n + h, y_n + h k₃) y_{n+1} = y_n + (h/6)(k₁ + 2 k₂ + 2 k₃ + k₄)

4th-order accurate, four function evaluations per step. Worked example on ẏ = −2y, y(0) = 1, h = 0.1:

  • k₁ = −2·1.0 = −2.0
  • k₂ = −2·(1.0 + 0.05·(−2.0)) = −2·0.9 = −1.8
  • k₃ = −2·(1.0 + 0.05·(−1.8)) = −2·0.91 = −1.82
  • k₄ = −2·(1.0 + 0.1·(−1.82)) = −2·0.818 = −1.636
  • y_1 = 1.0 + (0.1/6)·(−2.0 + 2·(−1.8) + 2·(−1.82) + (−1.636)) = 1.0 + 0.01667·(−10.876) = 0.81873

vs. true 0.81873 — agreement to five places at the first step. Compare Euler’s 0.019 error at the same h.

3.6 Dormand-Prince RK4(5) — the standard adaptive choice

Dormand and Prince (J. Comput. Appl. Math. 1980) gave a seven-stage embedded pair: the fifth-order solution is propagated, the fourth-order solution provides an error estimate by difference. Implemented as MATLAB’s ode45, SciPy solve_ivp(method='RK45'), Julia’s Tsit5 / DP5. FSAL (“first same as last”) property — the last stage of step n is the first stage of step n+1 — so seven stages cost six evaluations per accepted step.

3.7 Cash-Karp RK4(5)

Cash and Karp (ACM TOMS 1990) — alternative six-stage embedded RK4(5); slightly different error coefficients, behaves better on some problems.

3.8 Verner methods

Verner constructed high-order embedded pairs — Verner 6(5), 8(7), 9(8). Common in DifferentialEquations.jl (Vern7, Vern9); preferred when tight tolerances (1e−10 or below) are required because per-step accuracy dominates over per-step cost.

3.9 DOP853

Dormand-Prince 8(5,3) — Hairer’s high-order workhorse. Excellent for smooth non-stiff problems demanding small error.

4. Stiff Solvers

A problem is stiff when its Jacobian has eigenvalues with widely separated real parts, the stiff modes being strongly stable. Explicit RK is then bound to taking ∼1/|λ_max| step even when the slow dynamics evolve on 1/|λ_min| timescale; ratios of 10⁶ – 10¹² are routine in chemical kinetics and circuit simulation.

4.1 BDF — Backward Differentiation Formulae

Gear (Numerical Initial Value Problems in ODEs, 1971). Multistep implicit:

BDF1 (= backward Euler): y_{n+1} − y_n = h f_{n+1} BDF2: (3/2) y_{n+1} − 2 y_n + (1/2) y_{n−1} = h f_{n+1} BDF3: (11/6) y_{n+1} − 3 y_n + (3/2) y_{n−1} − (1/3) y_{n−2} = h f_{n+1}

Orders 1 – 6; orders 1, 2 are A-stable, 3 – 6 are only A(α)-stable (some sliver near the imaginary axis is unstable, fine for dissipative systems). The implementation in SUNDIALS CVODE — variable-step, variable-order BDF with Newton + sparse linear solver — is the industrial gold standard.

4.2 Radau IIA

Implicit RK based on Radau quadrature; 3-stage = 5th order, A-stable and L-stable. Hairer’s RADAU5 reference Fortran code; SciPy solve_ivp(method='Radau'); Julia RadauIIA5. Better than BDF for highly oscillatory stiff problems because it damps fast modes hard.

4.3 Rosenbrock methods

Linearly implicit RK — the Jacobian is evaluated once per step, the stages are linear solves rather than nonlinear. Avoids the inner Newton iteration; useful when the Jacobian is expensive but well-behaved. MATLAB ode23s is a Rosenbrock-Wanner method. Trade-off: order-2 or order-4 only; less robust on hard stiff problems than BDF/Radau.

4.4 TR-BDF2

Trapezoidal-rule + 2nd-order BDF composite (Bank, Coughran et al. 1985). MATLAB ode23t (TR alone, non-stiff) and ode23tb (TR-BDF2, stiff). Designed for moderately stiff problems with discontinuities, common in circuit simulation.

4.5 LSODA

Petzold + Hindmarsh: automatically switches between Adams (non-stiff) and BDF (stiff) based on online stiffness detection. SciPy solve_ivp(method='LSODA'); widely used because the user does not have to classify the problem ahead of time.

4.6 Identifying stiffness

  • Run with ode45 / RK45; if the solver takes thousands of tiny steps over a smooth region, suspect stiffness.
  • Estimate eigenvalues of J = ∂f/∂y. If max|λ| / min|λ| is large (and max Re λ < 0), the problem is stiff.
  • Chemical kinetics example: H₂ + O₂ combustion — radical species evolve at 10⁹ s⁻¹, bulk species at 10⁰ s⁻¹.

5. Multistep Methods

Use prior values y_{n−k} to extrapolate. Cheap (one f-evaluation per step) but require startup with a one-step method.

5.1 Adams-Bashforth (explicit)

AB2: y_{n+1} = y_n + (h/2)(3 f_n − f_{n−1}) AB3: y_{n+1} = y_n + (h/12)(23 f_n − 16 f_{n−1} + 5 f_{n−2}) AB4: y_{n+1} = y_n + (h/24)(55 f_n − 59 f_{n−1} + 37 f_{n−2} − 9 f_{n−3})

5.2 Adams-Moulton (implicit)

AM2 (trapezoidal): y_{n+1} = y_n + (h/2)(f_{n+1} + f_n) AM3: y_{n+1} = y_n + (h/12)(5 f_{n+1} + 8 f_n − f_{n−1}) AM4: y_{n+1} = y_n + (h/24)(9 f_{n+1} + 19 f_n − 5 f_{n−1} + f_{n−2})

5.3 Predictor-corrector (ABMx)

Predict with AB, correct with AM evaluated at the predicted value. Two function evaluations per step; gives an estimate of local error as the predictor-corrector difference — used by CVODE’s Adams mode.

5.4 Gear’s method

Synonym for BDF in older literature; Gear’s 1971 textbook normalised the variable-order, variable-step implicit-multistep approach now standard in stiff codes.

6. Adaptive Step Size + Error Control

Embedded pairs produce y_{n+1}^{(p)} (lower order) and y_{n+1}^{(p+1)} (higher order). Local error estimate:

err = ‖y_{n+1}^{(p+1)} − y_{n+1}^{(p)}‖

Per-component tolerance:

sc_i = atol + rtol · max(|y_n,i|, |y_{n+1,i}|) err_norm = sqrt((1/N) Σ (err_i / sc_i)²) (RMS form used by Hairer)

Step controller (PI controller, Gustafsson 1994):

h_new = h · safety · (1/err_norm)^{1/(p+1)} · (err_norm_prev / err_norm)^{α}

Typical values: safety = 0.9, α ≈ 0.1, p = 4 for RK45. Step is rejected if err_norm > 1 and retried at h_new ≥ h/10; accepted otherwise, and h_new ≤ 10 h.

Typical user-facing tolerances:

  • Smooth engineering: rtol = 1e−3, atol = 1e−6.
  • Astrodynamics: rtol = 1e−10, atol scaled per body.
  • ML neural-ODE training: rtol = 1e−5, atol = 1e−5 for stability of the adjoint.

7. System ODEs

Higher-order ODEs become first-order systems. For the resonator m·y” + c·y’ + k·y = F(t), define state x = [y, y’]ᵀ:

ẋ = [ x₂, (F(t) − c·x₂ − k·x₁) / m ]ᵀ

For controlled LTI plants: ẋ = A x + B u, with A ∈ ℝ^{n×n}, B ∈ ℝ^{n×m}. Nonlinear plants: ẋ = f(x, u, t). Software solvers all accept f as a vector-valued function; users implement f and pass it in.

For very large systems (n > 10⁴) — e.g. spatial discretisations of PDEs — sparse Jacobians and Krylov linear solvers (GMRES inside Newton) become essential. SUNDIALS CVODE supports preconditioned Krylov via SPGMR.

8. Stability

For a linear test problem ẏ = λ y, λ ∈ ℂ, a method with stability function R(z), z = h λ, is stable when |R(z)| ≤ 1. The stability region is { z : |R(z)| ≤ 1 }.

MethodStability regionProperties
Forward Eulerdisk:1 + z
Backward Euler1 − z
Trapezoidalleft half plane exactlyA-stable; R(∞) = 1 → not L-stable
RK4irregular region, bounded ≈ circle radius 2.78 along negative real axisexplicit; not A-stable
BDF2exterior of a cardioid in right half planeA-stable
BDF3-6not A-stable; A(α)-stable with α decreasing in orderusable on dissipative stiff problems
Radau IIA (3-stage)entire left half plane; R(∞) = 0A-stable and L-stable

A-stable — stability region ⊇ {Re z ≤ 0}. Necessary for unconditionally stable stiff integration.

L-stable — A-stable AND R(∞) = 0. Required when fast modes should be damped to zero rather than oscillate; canonical example is the immediate-post-pulse response in a chemical reactor.

Dahlquist’s second barrier theorem (1963): no explicit linear multistep method can be A-stable; among A-stable LMMs the order is at most 2 (trapezoidal is the optimal one). This is why stiff solvers are all either implicit RK or implicit multistep — no escape via cleverer explicit construction.

9. Boundary Value Problems

BVPs specify y at multiple points: y(a), y(b), or mixed conditions like y’(a) − k y(a) = 0. Examples: steady-state heat conduction in a rod, beam deflection, optimal-control two-point boundary problems via Pontryagin.

9.1 Shooting method

Treat the BVP as a parameterised IVP: guess the missing initial slope s, integrate y(b; s), iterate s via Newton or secant on the residual y(b; s) − y_target. Cheap when the dynamics are non-stiff and weakly nonlinear; explodes when the IVP is unstable.

9.2 Multiple shooting

Partition [a, b] into sub-intervals, shoot from each, add continuity constraints. Each sub-shot is bounded so the IVP instability does not explode; equality constraints solved as a block-tridiagonal Newton system. Bock + Plitt 1984 (in optimal control / Maximum Principle).

9.3 Finite difference + collocation

Discretise y on a mesh, replace derivatives by FD stencils or collocation polynomials, solve the resulting (usually sparse, banded) nonlinear system by Newton. MATLAB bvp4c / bvp5c use Lobatto collocation of order 4 and 5. SciPy solve_bvp is a collocation method (Kierzenka + Shampine 2008). Robust on stiff BVPs.

9.4 Sparse linear solves

The Jacobian of the discretised BVP is block-tridiagonal (for collocation on adjacent intervals) or banded (for FD). Sparse direct solvers (SuiteSparse, UMFPACK, MUMPS) make this O(n) instead of O(n³).

10. Delay Differential Equations

y’(t) = f(t, y(t), y(t − τ₁), …, y(t − τ_k)), with history y(t) = φ(t) for t ≤ t₀.

The right-hand side depends on past values. Naively interpolate y on the integration history (typically using a Hermite cubic on the step interval); the order of the underlying RK becomes the order of the DDE solver, provided the interpolant order matches.

Software:

  • MATLAB dde23 (RK2(3) with cubic Hermite interpolation, Shampine 2001).
  • DifferentialEquations.jl MethodOfSteps(Tsit5()) — wraps the user-chosen ODE method.
  • Discontinuities at multiples of τ propagate forward and must be handled by tracking and either stepping exactly to them or restarting.

Applications: control with transport delay (chemical mixing pipes, network round-trip), epidemiology with incubation period, population dynamics with gestation.

11. Stochastic Differential Equations

dX_t = µ(X_t, t) dt + σ(X_t, t) dW_t

with W_t a Wiener process. Two interpretations of the noise integral:

  • Itô — non-anticipating; ⟨W_t² ⟩ = t; natural for finance.
  • Stratonovich — symmetric; chain rule from ordinary calculus holds; natural for physics where σ is a smoothed-out fast process.

Conversion: Itô SDE dX = µ dt + σ dW ⇔ Stratonovich dX = (µ − ½ σ σ’) dt + σ ∘ dW.

11.1 Euler-Maruyama

X_{n+1} = X_n + µ(X_n, t_n) · Δt + σ(X_n, t_n) · ΔW_n, ΔW_n ∼ N(0, Δt)

Strong order 0.5, weak order 1.

11.2 Milstein

Adds the Itô-correction term:

X_{n+1} = X_n + µ Δt + σ ΔW + ½ σ σ’ (ΔW² − Δt)

Strong order 1.0; needs ∂σ/∂x.

11.3 Higher-order

Stochastic RK (Rößler 2010, Burrage + Burrage). Tau-leaping for jump-diffusion. DifferentialEquations.jl SDE solvers (SOSRI, SRA) are the most actively developed.

11.4 Applications

  • Black-Scholes geometric Brownian motion dS = µ S dt + σ S dW, closed-form solution for European options; numerical for path-dependent.
  • Langevin dynamics in molecular simulation: m ẍ = −∇U − γ ẋ + √(2 γ k_B T) ξ(t).
  • Neural SDEs (Tzen + Raginsky 2019, Li et al. 2020) — drift and diffusion parameterised by neural nets; trained by stochastic adjoint.
  • Score-based diffusion models for image generation (Song et al. 2021) — reverse-time SDE / probability-flow ODE sampling.

12. Symplectic Integrators

For Hamiltonian systems H(q, p) with q̇ = ∂H/∂p, ṗ = −∂H/∂q, standard methods drift in energy on long integrations. Symplectic integrators preserve the canonical 2-form dq ∧ dp; the numerical energy oscillates near the true value forever instead of drifting.

12.1 Leap-frog / Störmer-Verlet

For separable H = T(p) + V(q):

p_{n+1/2} = p_n − (h/2) ∇V(q_n) q_{n+1} = q_n + h · ∇T(p_{n+1/2}) [= h p_{n+1/2}/m for T = p²/2m] p_{n+1} = p_{n+1/2} − (h/2) ∇V(q_{n+1})

2nd order, time-reversible, symplectic. Backbone of MD codes for fifty years.

12.2 Velocity Verlet

Algebraically equivalent rearrangement; the form used in LAMMPS, GROMACS, AMBER, NAMD.

12.3 Yoshida 4th-order

Composition of three 2nd-order Verlet sub-steps with carefully chosen coefficients (Yoshida 1990):

w₁ = 1 / (2 − 2^{1/3}), w₀ = 1 − 2 w₁ step = Verlet(w₁ h) ∘ Verlet(w₀ h) ∘ Verlet(w₁ h)

12.4 Implicit symplectic — Gauss-Legendre RK

s-stage Gauss-Legendre is symplectic and of order 2s, A-stable, suitable for stiff Hamiltonian problems but expensive per step.

12.5 Applications

  • Molecular dynamics (LAMMPS, GROMACS, OpenMM) — protein folding, materials, drug discovery.
  • N-body astrophysics (Barnes-Hut tree codes, REBOUND).
  • Hamiltonian Monte Carlo in Bayesian inference (Neal 2010, Stan, NumPyro) — leap-frog is the proposal step.
  • Plasma particle-in-cell simulation (Boris pusher is a symplectic-like scheme for charged-particle motion).

13. Software Ecosystem

StackMethodsNotes
SciPy solve_ivpRK23, RK45, DOP853, Radau, BDF, LSODAPython default; events, dense output, vectorised f.
MATLABode45 (DP5), ode23 (BS23), ode113 (Adams), ode15s (NDF/BDF), ode23s (Rosenbrock), ode23t (TR), ode23tb (TR-BDF2)Shampine + Reichelt 1997 suite.
SUNDIALS (LLNL)CVODE (BDF + Adams), CVODES (sensitivity), IDA (DAE), IDAS, ARKode (ARK explicit + implicit + IMEX), KINSOLC library; the gold standard in production. Fortran/Python bindings (scikits.odes, sundialsTB).
DifferentialEquations.jlDozens of RK, multistep, IMEX, exponential, symplectic, SDE, DDE, DAE, jump processesRackauckas + Nie 2017; the most comprehensive ecosystem; GPU via DiffEqGPU; adjoint via SciMLSensitivity.
torchdiffeqRK4, DOPRI5, adaptive Adams, implicit Adams; full adjointChen + Rubanova et al. 2018 NeurIPS — Neural ODEs paper.
torchdynHigher-level neural-ODE / neural-CDE / neural-SDE on top of torchdiffeq
diffrax (JAX)Tsit5, DOPRI5, DOPRI8, Kvaerno3/4/5 (stiff), Heun, Euler; SDE; adjointKidger 2022; JIT + autodiff + GPU/TPU; the JAX-native counterpart.
Boost.odeint (C++)RK4, RK54-CK, RK78-Fehlberg, controlled adaptive, symplecticHeader-only; embedded / robotics simulation.
GSL gsl_odeiv2RK45, RK45-CK, RK8PD, BSIMPOlder C reference; well-tested.

For embedded / real-time systems (motion control, MPC at kHz rates), hand-coded RK4 in C is the norm — variable-step controllers add too much jitter.

14. Applications

14.1 Robotics dynamics

Manipulator dynamics M(q) q̈ + C(q, q̇) q̇ + G(q) = τ is a stiff-ish second-order ODE. Forward-dynamics simulation for training and offline planning uses RK4 or implicit Runge-Kutta when contact is involved (MuJoCo, Bullet, Drake’s RK3 + Radau options). Model-predictive-control rollouts must run faster than the control rate; RK4 with a fixed step is typical.

14.2 Vehicle dynamics + flight simulation

6-DoF rigid-body equations plus aerodynamics; tens of states. DOP853 or DP5; tolerances 1e−6. Real-time flight sim runs at fixed 60 – 240 Hz with RK4 — simulation fidelity traded for determinism.

14.3 Chemical reactor models

PFR / CSTR with multiple species and Arrhenius kinetics — stiffness ratios of 10⁸ – 10¹². Cantera + SUNDIALS CVODE BDF is the standard pipeline.

14.4 Electrical circuits + power systems

SPICE-class simulators (ngspice, Xyce, LTspice) use TR-BDF2 or BDF for transient analysis. Power-grid transient stability simulations use trapezoidal rule with corrector iterations on networks of millions of buses.

14.5 Population biology + epidemiology

SIR/SEIR ODE systems are small and non-stiff — RK45 suffices. Spatial epidemiology adds reaction-diffusion PDEs; method-of-lines + stiff solver.

14.6 Neural ODEs

Chen, Rubanova, Bettencourt, Duvenaud (NeurIPS 2018): treat a residual block as one Euler step h_{n+1} = h_n + f(h_n; θ_n) and take the continuous limit dh/dt = f(h, t; θ). Training uses the adjoint sensitivity method — solve a second ODE backward in time for the gradient — giving O(1) memory in depth. Applications: continuous normalising flows, irregular-time-series modelling (latent-ODE).

14.7 Diffusion models

Score-based generative models (Song et al. 2021): forward SDE adds Gaussian noise over t ∈ [0, T], reverse SDE / probability-flow ODE removes it. Sampling is ODE integration of dx/dt = −½ β(t) x − ½ β(t) ∇_x log p_t(x). DPM-Solver (Lu et al. 2022), DEIS, UniPC — specialised exponential-integrator solvers cut sampling steps from 1000 to 10 – 25.

14.8 Optimal control + reachability

Pontryagin’s maximum principle produces a coupled state/adjoint BVP. Direct shooting + sequential quadratic programming (CasADi, ACADOS) and collocation methods (GPOPS-II) are both used.

15. Common Pitfalls

  • Stiff problem with non-stiff solver — RK45 takes microscopic steps or blows up. Symptom: solver step count is huge despite smooth-looking solution. Switch to BDF / Radau or use LSODA’s autoselect.
  • Tolerance too tight — wastes compute, sometimes amplifies round-off in finite-precision arithmetic. Too loose — silently incorrect solutions, often by orders of magnitude. Always sanity-check by halving the tolerance and confirming the answer is stable.
  • Discontinuous f — RHS jumps at known events (impact, switch toggle, valve open) destroy a method’s order assumption. Use event detection (root-finding on a switching function); SciPy events=, MATLAB Events. Integrate up to the event, restart the solver after.
  • Jacobian not provided — implicit solvers fall back to finite-difference Jacobian: slow, inaccurate, sometimes non-convergent. Provide analytic Jacobian (or AD-generated Jacobian via JAX / Zygote / PyTorch). Sparsity patterns matter for large systems — declare them.
  • Conservation laws drift — energy, momentum, mass conservation violated by standard methods. Use symplectic for Hamiltonians; use projection methods (project onto the manifold after each step) for general invariants.
  • Initial-condition mismatch in DAEs — index-1 DAEs need consistent initial conditions. IDA’s IDACalcIC; Sundials emits a clear error otherwise.
  • Mass matrix not symmetric/positive definite — implicit methods may fail to converge. Reformulate (e.g. constrained Lagrangian → quasi-coordinate form) before solving.
  • Random seed not fixed in SDE — irreproducible runs in research; always set the seed and document it.
  • Sensitivity / adjoint instability — the backward-in-time adjoint ODE can be stiff even when the forward problem is not, particularly near a fixed point. Use a stiff solver in the adjoint or use checkpointing.

16. Cross-References

  • [[Math/_index]] — discipline overview.
  • [[Math/linear-algebra-essentials]] — eigenvalues / eigenvectors for linear ODEs and stability analysis.
  • [[Math/calculus-multivariate]] — vector calculus, Jacobians, divergence/curl underpin ODE systems and conservation laws (TBD).
  • [[Math/convex-optimization]] — inner Newton iterations in implicit solvers; sensitivity-based optimisation.
  • [[Math/probability-fundamentals]] — stochastic integration, Itô vs Stratonovich, Brownian motion.
  • [[Engineering/classical-control]] — Laplace-domain solutions, transfer functions, second-order resonator.
  • [[Engineering/vibration-dynamics]] — mechanical second-order systems, modal analysis.
  • [[Engineering/cfd-deep]] — PDE temporal advance uses ODE integrators after spatial discretisation (method of lines).

17. Citations

  • Hairer, E., Nørsett, S. P., Wanner, G. — Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer 1993.
  • Hairer, E., Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer 1996. (The canonical pair.)
  • Ascher, U. M., Petzold, L. R. — Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM 1998.
  • Butcher, J. C. — Numerical Methods for Ordinary Differential Equations, 3rd ed., Wiley 2016.
  • Dormand, J. R., Prince, P. J. — “A family of embedded Runge-Kutta formulae”, J. Comput. Appl. Math. 6(1):19-26, 1980.
  • Cash, J. R., Karp, A. H. — “A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides”, ACM TOMS 16(3):201-222, 1990.
  • Gear, C. W. — Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall 1971.
  • Shampine, L. F., Reichelt, M. W. — “The MATLAB ODE Suite”, SIAM J. Sci. Comput. 18(1):1-22, 1997.
  • Shampine, L. F. — “Solving DDEs in MATLAB”, Appl. Numer. Math. 37:441-458, 2001.
  • Kierzenka, J., Shampine, L. F. — “A BVP Solver that Controls Residual and Error”, JNAIAM 3:27-41, 2008.
  • Yoshida, H. — “Construction of higher order symplectic integrators”, Phys. Lett. A 150(5-7):262-268, 1990.
  • Hindmarsh, A. C. et al. — “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers”, ACM TOMS 31(3):363-396, 2005.
  • Rackauckas, C., Nie, Q. — “DifferentialEquations.jl — A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia”, J. Open Research Software 5:15, 2017.
  • Chen, R. T. Q., Rubanova, Y., Bettencourt, J., Duvenaud, D. — “Neural Ordinary Differential Equations”, NeurIPS 2018.
  • Kidger, P. — “On Neural Differential Equations”, PhD thesis Oxford / diffrax library, 2022.
  • Song, Y. et al. — “Score-Based Generative Modeling through Stochastic Differential Equations”, ICLR 2021.
  • Lu, C. et al. — “DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps”, NeurIPS 2022.
  • Bank, R. E. et al. — “Transient simulation of silicon devices and circuits”, IEEE TCAD 4(4):436-451, 1985 (TR-BDF2).
  • Gustafsson, K. — “Control-theoretic techniques for stepsize selection in implicit Runge-Kutta methods”, ACM TOMS 20(4):496-517, 1994.
  • Dahlquist, G. — “A special stability problem for linear multistep methods”, BIT 3:27-43, 1963 (the second Dahlquist barrier).

Graph View