ODE Integration Methods

Marching forward in time: from Euler's simple step to adaptive Runge-Kutta, implicit methods for stiff systems, and symplectic integrators that conserve energy for Hamiltonian dynamics.

The Initial Value Problem

Given \(\frac{dy}{dt} = f(t, y)\) with \(y(t_0) = y_0\), compute \(y(t)\) for\(t > t_0\). This is the core problem of dynamics: planetary orbits, chemical kinetics, neural networks, population models, circuit simulation.

Key questions: How small must \(h = \Delta t\) be? Is the method stable? Does it conserve physical invariants (energy, symplectic structure)?

1. Euler's Method

The simplest ODE integrator: step along the tangent line at the current point.

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

First-order method: local error \(O(h^2)\), global error \(O(h)\)

Error Analysis

Taylor expand: \(y(t+h) = y(t) + h y'(t) + \frac{h^2}{2} y''(t) + \cdots\)Euler keeps only the first two terms, so the local truncation error is\(\tau = \frac{h^2}{2} y''(\xi) = O(h^2)\). Over \(N = T/h\) steps, errors accumulate to give global error \(O(h)\).

Euler Method: Simple but Instructive

Python

script.py36 lines

import numpy as np

def euler(f, y0, t_span, h):
    """Forward Euler method."""
    t0, tf = t_span
    t = np.arange(t0, tf + h/2, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for i in range(len(t) - 1):
        y[i+1] = y[i] + h * np.atleast_1d(f(t[i], y[i]))
    return t, y

# Test: dy/dt = -y, y(0) = 1 => exact: y = exp(-t)
f = lambda t, y: -y
y0 = np.array([1.0])

print("Euler Method: dy/dt = -y, y(0) = 1")
print(f"Exact solution: y(t) = exp(-t)\n")

print(f"{'h':>8}  {'y(1) Euler':>14}  {'y(1) exact':>14}  {'|error|':>12}  {'order':>6}")
print("-" * 60)

prev_err = None
for h in [0.2, 0.1, 0.05, 0.025, 0.0125]:
    t, y = euler(f, y0, (0, 1), h)
    y_euler = y[-1, 0]
    y_exact = np.exp(-1.0)
    err = abs(y_euler - y_exact)
    if prev_err is not None:
        order = np.log2(prev_err / err)
        print(f"{h:8.4f}  {y_euler:14.10f}  {y_exact:14.10f}  {err:12.2e}  {order:6.2f}")
    else:
        print(f"{h:8.4f}  {y_euler:14.10f}  {y_exact:14.10f}  {err:12.2e}  {'---':>6}")
    prev_err = err

print("\nOrder ~1.0 confirmed: halving h halves the error")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. The Classic Runge-Kutta 4 (RK4)

The workhorse of ODE integration. Four function evaluations per step yield fourth-order accuracy:

\(k_1 = h \, f(t_n, y_n)\)

\(k_2 = h \, f(t_n + h/2, \; y_n + k_1/2)\)

\(k_3 = h \, f(t_n + h/2, \; y_n + k_2/2)\)

\(k_4 = h \, f(t_n + h, \; y_n + k_3)\)

\(y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\)

Local error \(O(h^5)\), global error \(O(h^4)\). The most popular fixed-step ODE solver.

RK4 Implementation and Convergence Test

Python

script.py47 lines

import numpy as np

def rk4(f, y0, t_span, h):
    """Classic 4th-order Runge-Kutta."""
    t0, tf = t_span
    t = np.arange(t0, tf + h/2, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for i in range(len(t) - 1):
        yi = y[i]
        ti = t[i]
        k1 = h * np.atleast_1d(f(ti, yi))
        k2 = h * np.atleast_1d(f(ti + h/2, yi + k1/2))
        k3 = h * np.atleast_1d(f(ti + h/2, yi + k2/2))
        k4 = h * np.atleast_1d(f(ti + h, yi + k3))
        y[i+1] = yi + (k1 + 2*k2 + 2*k3 + k4) / 6
    return t, y

# Test with a harder problem: harmonic oscillator
# y'' + y = 0 => y1' = y2, y2' = -y1
# Exact: y1 = cos(t), y2 = -sin(t)
def harmonic(t, y):
    return np.array([y[1], -y[0]])

y0 = np.array([1.0, 0.0])  # cos(0) = 1, -sin(0) = 0

print("RK4: Harmonic Oscillator y'' + y = 0")
print(f"Initial: y(0) = 1, y'(0) = 0")
print(f"Exact: y(t) = cos(t)\n")

print(f"{'h':>10}  {'y(2pi) RK4':>14}  {'|error|':>12}  {'order':>6}")
print("-" * 48)

prev_err = None
for h in [0.2, 0.1, 0.05, 0.025, 0.0125]:
    t, y = rk4(harmonic, y0, (0, 2*np.pi), h)
    y_rk4 = y[-1, 0]
    y_exact = np.cos(2*np.pi)  # = 1.0
    err = abs(y_rk4 - y_exact)
    if prev_err is not None and err > 0:
        order = np.log2(prev_err / err)
        print(f"{h:10.4f}  {y_rk4:14.10f}  {err:12.2e}  {order:6.2f}")
    else:
        print(f"{h:10.4f}  {y_rk4:14.10f}  {err:12.2e}  {'---':>6}")
    prev_err = err

print("\nOrder ~4.0 confirmed: halving h reduces error by 16x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Adaptive Step Size: RK45

Fixed step sizes waste effort in smooth regions and miss features in rapidly changing regions. Embedded Runge-Kutta pairs (like Dormand-Prince RK45) compute two solutions of different order to estimate the local error, then adjust \(h\) automatically.

Step Size Control

Estimate error \(\hat{e} = |y_5 - y_4|\) (difference between 5th and 4th order solutions). If \(\hat{e} > \text{tol}\), reject step and shrink \(h\). If \(\hat{e} \ll \text{tol}\), accept and grow \(h\). Optimal step factor: \(h_{\text{new}} = h \cdot 0.9 \cdot (\text{tol}/\hat{e})^{1/5}\).

Adaptive ODE Solving with scipy.integrate

Python

script.py43 lines

import numpy as np
from scipy.integrate import solve_ivp

# Stiff-ish problem: Van der Pol oscillator
# y'' - mu*(1-y^2)*y' + y = 0
mu = 10  # moderately stiff

def vanderpol(t, y):
    return [y[1], mu * (1 - y[0]**2) * y[1] - y[0]]

y0 = [2.0, 0.0]
t_span = (0, 3 * mu)  # several periods

# Solve with RK45 (adaptive Dormand-Prince)
sol_rk45 = solve_ivp(vanderpol, t_span, y0, method='RK45',
                      rtol=1e-8, atol=1e-10, dense_output=True)

print(f"Van der Pol oscillator (mu = {mu})")
print(f"\nRK45 (adaptive Dormand-Prince):")
print(f"  Time steps used: {len(sol_rk45.t)}")
print(f"  Function evaluations: {sol_rk45.nfev}")
print(f"  Status: {'success' if sol_rk45.success else 'failed'}")

# Show step sizes vary dramatically
dt = np.diff(sol_rk45.t)
print(f"  Min step size: {dt.min():.6f}")
print(f"  Max step size: {dt.max():.6f}")
print(f"  Ratio max/min: {dt.max()/dt.min():.1f}x")

# Compare with fixed-step RK4 (would need many more steps)
# For same accuracy, fixed step needs min(dt) everywhere
n_fixed = int((t_span[1] - t_span[0]) / dt.min())
print(f"\n  Fixed-step would need ~{n_fixed} steps")
print(f"  Adaptive used only {len(sol_rk45.t)} steps")
print(f"  Speedup: ~{n_fixed / len(sol_rk45.t):.0f}x")

# Sample solution at regular intervals
t_eval = np.linspace(0, t_span[1], 20)
y_eval = sol_rk45.sol(t_eval)
print(f"\n{'t':>8}  {'y':>10}  {'dy/dt':>10}")
print("-" * 32)
for i in range(len(t_eval)):
    print(f"{t_eval[i]:8.2f}  {y_eval[0,i]:10.4f}  {y_eval[1,i]:10.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Stiff Systems and Implicit Methods

A system is stiff when it contains both fast and slow timescales. Explicit methods must use tiny steps (dictated by the fast scale) even when tracking the slow dynamics. Implicit methods like backward Euler and BDF (backward differentiation formulas) remain stable with large steps.

Backward Euler (implicit)

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

A-stable: stable for any \(h > 0\) on \(y' = \lambda y\) with Re(\(\lambda\)) < 0

Forward Euler (explicit)

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

Stable only if \(|1 + h\lambda| < 1\), requires \(h < 2/|\lambda|\)

Stiff System: Explicit vs Implicit

Python

script.py46 lines

import numpy as np
from scipy.integrate import solve_ivp

# Classic stiff problem: chemical kinetics (Robertson)
# dy1/dt = -0.04*y1 + 1e4*y2*y3
# dy2/dt =  0.04*y1 - 1e4*y2*y3 - 3e7*y2^2
# dy3/dt =  3e7*y2^2
# Timescales differ by ~10^7!

def robertson(t, y):
    y1, y2, y3 = y
    dy1 = -0.04*y1 + 1e4*y2*y3
    dy2 =  0.04*y1 - 1e4*y2*y3 - 3e7*y2**2
    dy3 =  3e7*y2**2
    return [dy1, dy2, dy3]

y0 = [1.0, 0.0, 0.0]
t_span = (0, 1e5)

# BDF (implicit, designed for stiff problems)
sol_bdf = solve_ivp(robertson, t_span, y0, method='BDF',
                     rtol=1e-8, atol=1e-10)

# Radau (implicit Runge-Kutta, also good for stiff)
sol_radau = solve_ivp(robertson, t_span, y0, method='Radau',
                       rtol=1e-8, atol=1e-10)

# RK45 would struggle or fail on this problem
sol_rk45 = solve_ivp(robertson, t_span, y0, method='RK45',
                      rtol=1e-8, atol=1e-10, max_step=100)

print("Robertson Chemical Kinetics (stiff system)")
print(f"Timescale ratio: ~10^7\n")
print(f"{'Method':>10}  {'Steps':>8}  {'f-evals':>8}  {'Jac-evals':>10}  {'Status':>10}")
print("-" * 52)
print(f"{'BDF':>10}  {len(sol_bdf.t):8d}  {sol_bdf.nfev:8d}  {sol_bdf.njev:10d}  {'OK' if sol_bdf.success else 'FAIL':>10}")
print(f"{'Radau':>10}  {len(sol_radau.t):8d}  {sol_radau.nfev:8d}  {sol_radau.njev:10d}  {'OK' if sol_radau.success else 'FAIL':>10}")
print(f"{'RK45':>10}  {len(sol_rk45.t):8d}  {sol_rk45.nfev:8d}  {'N/A':>10}  {'OK' if sol_rk45.success else 'FAIL':>10}")

print(f"\nBDF and Radau handle stiff problems efficiently.")
print(f"They use Jacobian information to take large stable steps.")
print(f"\nFinal concentrations (t = {t_span[1]:.0e}):")
print(f"  y1 = {sol_bdf.y[0,-1]:.6f}")
print(f"  y2 = {sol_bdf.y[1,-1]:.6e}")
print(f"  y3 = {sol_bdf.y[2,-1]:.6f}")
print(f"  Conservation: y1+y2+y3 = {sum(sol_bdf.y[:,-1]):.10f} (should be 1.0)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Symplectic Integrators

For Hamiltonian systems \(H(q,p) = T(p) + V(q)\), standard integrators cause energy to drift over long times. Symplectic integrators preserve the symplectic structure of phase space, keeping energy bounded (oscillating, not drifting) forever.

Stormer-Verlet (Leapfrog)

\(p_{n+1/2} = p_n - \frac{h}{2} \nabla V(q_n)\) (half-kick)

\(q_{n+1} = q_n + h \, p_{n+1/2} / m\) (drift)

\(p_{n+1} = p_{n+1/2} - \frac{h}{2} \nabla V(q_{n+1})\) (half-kick)

Second-order, time-reversible, symplectic. Used in molecular dynamics, celestial mechanics, particle physics.

Symplectic vs Non-Symplectic: Energy Conservation

Python

script.py64 lines

import numpy as np

def euler_orbit(q0, p0, h, N):
    """Forward Euler (non-symplectic) for Kepler problem."""
    q, p = q0.copy(), p0.copy()
    energies = []
    for _ in range(N):
        r = np.sqrt(q[0]**2 + q[1]**2)
        E = 0.5*(p[0]**2 + p[1]**2) - 1.0/r
        energies.append(E)
        # Euler step
        a = -q / r**3
        q_new = q + h * p
        p_new = p + h * a
        q, p = q_new, p_new
    return energies

def verlet_orbit(q0, p0, h, N):
    """Stormer-Verlet (symplectic) for Kepler problem."""
    q, p = q0.copy(), p0.copy()
    energies = []
    for _ in range(N):
        r = np.sqrt(q[0]**2 + q[1]**2)
        E = 0.5*(p[0]**2 + p[1]**2) - 1.0/r
        energies.append(E)
        # Half-kick
        a = -q / r**3
        p = p + 0.5 * h * a
        # Drift
        q = q + h * p
        # Half-kick
        r = np.sqrt(q[0]**2 + q[1]**2)
        a = -q / r**3
        p = p + 0.5 * h * a
    return energies

# Kepler problem: circular orbit
q0 = np.array([1.0, 0.0])
p0 = np.array([0.0, 1.0])
h = 0.05
N = 2000  # ~16 orbits

E_euler = euler_orbit(q0, p0, h, N)
E_verlet = verlet_orbit(q0, p0, h, N)

E0 = E_euler[0]
drift_euler = [abs(E - E0) for E in E_euler]
drift_verlet = [abs(E - E0) for E in E_verlet]

print("Kepler Orbit: Symplectic vs Non-Symplectic")
print(f"Step size h = {h}, {N} steps (~{N*h/(2*np.pi):.0f} orbits)\n")

print(f"{'Orbit':>6}  {'Euler dE':>14}  {'Verlet dE':>14}")
print("-" * 38)
steps_per_orbit = int(2 * np.pi / h)
for orbit in range(0, N // steps_per_orbit + 1, 2):
    idx = min(orbit * steps_per_orbit, N - 1)
    print(f"{orbit:6d}  {drift_euler[idx]:14.6e}  {drift_verlet[idx]:14.6e}")

print(f"\nEuler: energy drifts monotonically (WRONG for long times)")
print(f"Verlet: energy oscillates but stays bounded (CORRECT)")
print(f"\nMax energy error - Euler:  {max(drift_euler):.6e}")
print(f"Max energy error - Verlet: {max(drift_verlet):.6e}")
print(f"Ratio: {max(drift_euler)/max(drift_verlet):.0f}x better with symplectic")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Stability Regions

For the test equation \(y' = \lambda y\), a method is stable if the numerical solution does not blow up. The stability region is the set of\(h\lambda \in \mathbb{C}\) for which \(|y_{n+1}/y_n| \leq 1\).

Forward Euler

\(y_{n+1} = (1 + h\lambda) y_n\), stable when \(|1 + h\lambda| \leq 1\): a disk of radius 1 centered at \(-1\).

Backward Euler

\(y_{n+1} = \frac{y_n}{1 - h\lambda}\), stable when \(|1 - h\lambda| \geq 1\): the complement of a disk. A-stable!

Stability Region Computation

Python

script.py52 lines

import numpy as np

# Compute stability regions for test equation y' = lambda*y
# Method is stable when |R(z)| <= 1 where z = h*lambda

z_real = np.linspace(-5, 3, 200)
z_imag = np.linspace(-4, 4, 200)
Z_re, Z_im = np.meshgrid(z_real, z_imag)
Z = Z_re + 1j * Z_im

# Forward Euler: R(z) = 1 + z
R_euler = np.abs(1 + Z)

# RK4: R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24
R_rk4 = np.abs(1 + Z + Z**2/2 + Z**3/6 + Z**4/24)

# Backward Euler: R(z) = 1/(1-z)
R_back = np.abs(1 / (1 - Z))

# Count stable points (|R| <= 1) along negative real axis
print("Stability on negative real axis (z = h*lambda, lambda < 0):")
print(f"{'h*lambda':>10}  {'Euler':>8}  {'RK4':>8}  {'BackEuler':>10}")
print("-" * 40)
for z_val in [-0.5, -1.0, -2.0, -2.5, -3.0, -4.0, -5.0]:
    z = complex(z_val, 0)
    r_e = abs(1 + z)
    r_rk4 = abs(1 + z + z**2/2 + z**3/6 + z**4/24)
    r_be = abs(1 / (1 - z))
    s_e = "STABLE" if r_e <= 1 else f"{r_e:.3f}"
    s_rk4 = "STABLE" if r_rk4 <= 1 else f"{r_rk4:.3f}"
    s_be = "STABLE" if r_be <= 1 else f"{r_be:.3f}"
    print(f"{z_val:10.1f}  {s_e:>8}  {s_rk4:>8}  {s_be:>10}")

# For imaginary axis (oscillatory problems)
print("\nStability on imaginary axis (z = i*omega*h):")
print(f"{'omega*h':>10}  {'Euler':>8}  {'RK4':>8}")
print("-" * 30)
for om in [0.5, 1.0, 2.0, 2.83, 3.0]:
    z = complex(0, om)
    r_e = abs(1 + z)
    r_rk4 = abs(1 + z + z**2/2 + z**3/6 + z**4/24)
    s_e = "STABLE" if r_e <= 1.0001 else f"{r_e:.3f}"
    s_rk4 = "STABLE" if r_rk4 <= 1.0001 else f"{r_rk4:.3f}"
    print(f"{om:10.2f}  {s_e:>8}  {s_rk4:>8}")

print("\nForward Euler is NEVER stable on imaginary axis!")
print("(That's why it can't do undamped oscillators)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server