Numerical ODE Solvers

Prerequisites

• Module 0.1: Calculus (Taylor series, integration)
• Module 1.1: First-Order ODEs (existence, uniqueness, analytical solutions for comparison)

Why This Matters

Most differential equations that arise in science and engineering have no closed-form solution. The Lorenz equations governing atmospheric convection, the Navier-Stokes equations for fluid flow, and the N-body equations of celestial mechanics all require numerical methods. Numerical ODE solvers are the workhorses of computational science: they power weather forecasting, spacecraft trajectory planning, circuit simulation, and molecular dynamics. Understanding how these solvers work—their accuracy, stability, and computational cost—is essential for anyone who uses or builds scientific software.

This module builds solvers from scratch so you understand exactly what happens inside scipy.integrate.solve_ivp and similar library functions. You will implement forward Euler, the classical fourth-order Runge-Kutta method (RK4), and an adaptive step-size controller, then compare their accuracy and efficiency on concrete test problems.

Learning Objectives

After completing this module you will be able to:

1. Implement forward Euler’s method from a Taylor-series derivation.
2. Implement the classical fourth-order Runge-Kutta method (RK4) and explain the role of each stage.
3. Define local truncation error and global truncation error, and relate them through convergence theory.
4. State the order of accuracy of a method and verify it numerically on a log-log error plot.
5. Explain A-stability and identify which methods are A-stable.
6. Implement adaptive step-size control using embedded Runge-Kutta pairs (RK4/5).
7. Compare solver accuracy versus computational cost for smooth, oscillatory, and near-stiff problems.
8. Use SciPy’s solve_ivp with appropriate method selection and tolerance settings.
9. Frame the numerical solver as an agent that chooses step sizes to balance accuracy and cost.

Core Concepts

Common Misconceptions

Worked Examples

Interactive Code Labs

import numpy as np
import matplotlib.pyplot as plt

# ── Right-hand side ──────────────────────────────────────────
def f(t, y):
    return -2.0 * y

# ── Exact solution ───────────────────────────────────────────
def exact(t):
    return np.exp(-2.0 * t)

# ── Forward Euler ────────────────────────────────────────────
def euler(f, t0, y0, T, h):
    t_vals = [t0]
    y_vals = [y0]
    t, y = t0, y0
    while t < T - 1e-12:
        y = y + h * f(t, y)
        t = t + h
        t_vals.append(t)
        y_vals.append(y)
    return np.array(t_vals), np.array(y_vals)

# ── Classical RK4 ────────────────────────────────────────────
def rk4(f, t0, y0, T, h):
    t_vals = [t0]
    y_vals = [y0]
    t, y = t0, y0
    while t < T - 1e-12:
        k1 = f(t, y)
        k2 = f(t + h/2, y + h/2 * k1)
        k3 = f(t + h/2, y + h/2 * k2)
        k4 = f(t + h, y + h * k3)
        y = y + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        t_vals.append(t)
        y_vals.append(y)
    return np.array(t_vals), np.array(y_vals)

# ── Solve with both methods ──────────────────────────────────
h = 0.1
T = 2.0
t_eu, y_eu = euler(f, 0, 1.0, T, h)
t_rk, y_rk = rk4(f, 0, 1.0, T, h)
t_ex = np.linspace(0, T, 200)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Left: solutions
axes[0].plot(t_ex, exact(t_ex), 'k-', lw=2, label='Exact')
axes[0].plot(t_eu, y_eu, 'ro--', ms=4, label=f'Euler h={h}')
axes[0].plot(t_rk, y_rk, 'bs--', ms=4, label=f'RK4 h={h}')
axes[0].set_xlabel('t'); axes[0].set_ylabel('y')
axes[0].legend(); axes[0].set_title('Solutions')
axes[0].grid(True, alpha=0.3)

# Right: error vs step size (log-log)
step_sizes = [0.5, 0.25, 0.1, 0.05, 0.025, 0.01]
euler_errors = []
rk4_errors = []

for hs in step_sizes:
    _, ye = euler(f, 0, 1.0, T, hs)
    _, yr = rk4(f, 0, 1.0, T, hs)
    euler_errors.append(abs(ye[-1] - exact(T)))
    rk4_errors.append(abs(yr[-1] - exact(T)))

axes[1].loglog(step_sizes, euler_errors, 'ro-', label='Euler (order 1)')
axes[1].loglog(step_sizes, rk4_errors, 'bs-', label='RK4 (order 4)')
# Reference slopes
h_ref = np.array(step_sizes)
axes[1].loglog(h_ref, 0.5*h_ref**1, 'r:', alpha=0.5, label='$O(h)$')
axes[1].loglog(h_ref, 0.5*h_ref**4, 'b:', alpha=0.5, label='$O(h^4)$')
axes[1].set_xlabel('Step size h'); axes[1].set_ylabel('|Error at t=2|')
axes[1].legend(); axes[1].set_title('Convergence Order')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print numerical convergence rates
print("Euler convergence rates:")
for i in range(1, len(step_sizes)):
    rate = np.log(euler_errors[i-1]/euler_errors[i]) / np.log(step_sizes[i-1]/step_sizes[i])
    print(f"  h={step_sizes[i]:.3f}: rate = {rate:.2f}")
print("RK4 convergence rates:")
for i in range(1, len(step_sizes)):
    rate = np.log(rk4_errors[i-1]/rk4_errors[i]) / np.log(step_sizes[i-1]/step_sizes[i])
    print(f"  h={step_sizes[i]:.3f}: rate = {rate:.2f}")

import numpy as np
import matplotlib.pyplot as plt

# ── Stiff-ish problem: y' = -1000(y - sin(t)) + cos(t) ──────
def f_stiff(t, y):
    return -1000.0 * (y - np.sin(t)) + np.cos(t)

def exact_stiff(t):
    # Exact: y(t) = sin(t) + C*exp(-1000t), y(0)=0 => C = -sin(0) = 0
    # Actually y(0)=0, sin(0)=0, so C=0 and y(t)=sin(t) for t>>0.001
    # More precisely: y(t) = sin(t) - sin(0)*exp(-1000t) = sin(t)
    return np.sin(t)

# ── Adaptive RK4/5 (simplified Dormand-Prince) ──────────────
def rk45_adaptive(f, t0, y0, T, tol=1e-6, h_init=0.001,
                  h_min=1e-8, h_max=0.5):
    t_vals = [t0]
    y_vals = [y0]
    h_vals = []
    t, y, h = t0, y0, h_init
    n_reject = 0

    while t < T - 1e-12:
        if t + h > T:
            h = T - t

        # RK4 step
        k1 = f(t, y)
        k2 = f(t + h/2, y + h/2 * k1)
        k3 = f(t + h/2, y + h/2 * k2)
        k4 = f(t + h, y + h * k3)
        y4 = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)

        # RK5 step (extra stage for error estimate)
        k5 = f(t + h/4, y + h/4 * k1)
        k6 = f(t + 3*h/4, y + h*(3*k1/8 + 9*k3/8 - 3*k2/4 + 3*k4/8))
        # Simplified 5th-order estimate
        y5 = y + (h/90) * (7*k1 + 32*k2 + 12*k3 + 32*k5 + 7*k4)

        # Error estimate
        err = abs(y5 - y4)
        if err < 1e-16:
            err = 1e-16

        # Step-size control
        if err <= tol:
            # Accept step
            t = t + h
            y = y4
            t_vals.append(t)
            y_vals.append(y)
            h_vals.append(h)
        else:
            n_reject += 1

        # Compute new step size
        safety = 0.9
        h_new = h * safety * (tol / err) ** 0.2
        h_new = max(h_min, min(h_max, h_new))
        h = h_new

    print(f"Adaptive: {len(t_vals)-1} steps, {n_reject} rejected")
    return np.array(t_vals), np.array(y_vals), np.array(h_vals)

# ── Fixed-step RK4 for comparison ────────────────────────────
def rk4_fixed(f, t0, y0, T, h):
    t_vals = [t0]
    y_vals = [y0]
    t, y = t0, y0
    while t < T - 1e-12:
        k1 = f(t, y)
        k2 = f(t + h/2, y + h/2 * k1)
        k3 = f(t + h/2, y + h/2 * k2)
        k4 = f(t + h, y + h * k3)
        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        t_vals.append(t)
        y_vals.append(y)
    print(f"Fixed RK4: {len(t_vals)-1} steps (h={h})")
    return np.array(t_vals), np.array(y_vals)

# ── Run both solvers ─────────────────────────────────────────
T = 1.0
ta, ya, ha = rk45_adaptive(f_stiff, 0, 0.0, T, tol=1e-6)
tf, yf = rk4_fixed(f_stiff, 0, 0.0, T, h=0.001)
t_ex = np.linspace(0, T, 1000)

fig, axes = plt.subplots(2, 1, figsize=(10, 7))

# Top: solutions
axes[0].plot(t_ex, exact_stiff(t_ex), 'k-', lw=2, label='Exact')
axes[0].plot(ta, ya, 'r.-', ms=3, lw=0.8, label=f'Adaptive ({len(ta)-1} steps)')
axes[0].plot(tf, yf, 'b-', lw=0.5, alpha=0.7, label=f'Fixed ({len(tf)-1} steps)')
axes[0].set_xlabel('t'); axes[0].set_ylabel('y')
axes[0].legend(); axes[0].set_title('Solution Comparison')
axes[0].grid(True, alpha=0.3)

# Bottom: step sizes chosen by adaptive solver
axes[1].semilogy(ta[:-1], ha, 'r.-', ms=3)
axes[1].axhline(y=0.001, color='b', ls='--', alpha=0.5, label='Fixed h=0.001')
axes[1].set_xlabel('t'); axes[1].set_ylabel('Step size h')
axes[1].set_title('Adaptive Step Sizes')
axes[1].legend(); axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Error comparison
err_adapt = abs(ya[-1] - exact_stiff(T))
err_fixed = abs(yf[-1] - exact_stiff(T))
print(f"Adaptive error at T={T}: {err_adapt:.2e}")
print(f"Fixed RK4 error at T={T}: {err_fixed:.2e}")

Agent Lens: The Solver as an Agent

An adaptive ODE solver can be viewed as a decision-making agent operating in a sequential environment:

Classical controllers as policies. The standard step-size formula is a proportional (P) controller on the log-error. More sophisticated PI and PID step-size controllers (Gustafsson 1991, Soderlind 2002) use the error history from the previous one or two steps, exactly analogous to PID control in engineering. These hand-crafted policies work well for most problems but are not optimal in any formal sense.

Could an RL agent do better? An RL agent could learn a step-size policy that accounts for problem-specific structure: for example, learning that periodic solutions allow aggressive step-size increases between peaks, or that boundary-layer regions require caution. Research (Dellnitz et al. 2021) has explored neural network–based step-size controllers that outperform classical heuristics on specific problem classes.

The exploration-exploitation trade-off. A conservative policy (small $h$) guarantees accuracy but wastes computation. An aggressive policy (large $h$) risks rejected steps. The optimal policy balances these, much like the exploration-exploitation dilemma in reinforcement learning.

Exercises

Analytical

Exercise 1. Show that the local truncation error of the forward Euler method applied to $y' = f(t,y)$ is $\tau_{n+1} = \frac{h^2}{2}y''(\xi_n)$ for some $\xi_n \in [t_n, t_{n+1}]$.

Exercise 2. For the test equation $y' = \lambda y$ with $\lambda \in \mathbb{C}$, find the stability function $R(z)$ where $z = \lambda h$ for (a) forward Euler, (b) backward Euler, (c) the trapezoidal rule.

Exercise 3. The Butcher tableau for the midpoint method is: $c = [0, 1/2]^T$, $A = [[0,0],[1/2,0]]$, $b = [0, 1]^T$. Verify that it satisfies the order-2 conditions: $\sum b_i = 1$ and $\sum b_i c_i = 1/2$.

Exercise 4. Prove that for a method of order $p$, the global error satisfies $|e_N| = O(h^p)$ where $N = (T-t_0)/h$, given that the LTE is $O(h^{p+1})$ and $f$ is Lipschitz in $y$.

Computational

Exercise 5. Implement the midpoint method ($y_{n+1} = y_n + hf(t_n + h/2, y_n + \frac{h}{2}f(t_n,y_n))$) and verify that it has order 2 by computing the convergence rate on $y' = -2y$.

Exercise 6. Implement backward Euler for $y' = -1000y$, $y(0) = 1$ using Newton’s method to solve the implicit equation at each step. Compare with forward Euler for $h = 0.01$ (is forward Euler stable? is backward Euler stable?).

Exercise 7. Solve the Lotka-Volterra system $x' = 1.5x - xy$, $y' = -3y + xy$ with both Euler and RK4 for $t \in [0, 15]$. Compare the phase portraits. Does Euler conserve the invariant?

Exercise 8. Use SciPy’s solve_ivp with method 'RK45' and method 'Radau' to solve $y' = -1000(y - \cos t) - \sin t$ on $[0, 10]$. Compare the number of function evaluations reported by sol.nfev.

Agentic

Exercise 9. Design a reward function for a step-size agent that penalizes (a) rejected steps, (b) excessive function evaluations, and (c) errors above tolerance. Run your adaptive solver on three test problems and compute the total reward. How does your reward function compare to the standard PI controller?

Exercise 10. Implement a “learning” step-size controller: keep a lookup table indexed by $\log_{10}(\text{err}/\text{tol})$ that stores the best step-size multiplier. After each step, update the table entry using exponential moving average. Does this “agent” learn a better policy than the standard formula over repeated solves of the same problem class?

Assessment

Quiz

Mini-Project: Solver Benchmark Dashboard

Build a benchmark comparing three solvers on three test problems.

Solvers: Forward Euler, RK4 (fixed step), adaptive RK4/5.

Test problems:

1. Smooth: $y' = -y$, $y(0) = 1$, $t \in [0, 5]$
2. Oscillatory: $y'' + 25y = 0$ (rewritten as a 2D system), $y(0) = 1, y'(0) = 0$, $t \in [0, 4\pi]$
3. Near-stiff: $y' = -500(y - \cos t) - \sin t$, $y(0) = 1$, $t \in [0, 5]$

Deliverables: For each solver–problem pair, produce (a) a solution plot, (b) the error at $t = T$, (c) the number of function evaluations. Summarize in a $3 \times 3$ table.

References & Next Steps

References

1. Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations, 3rd ed. Wiley.
2. Hairer, E., Nørsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems. Springer.
3. Hairer, E. & Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer.
4. Dormand, J. R. & Prince, P. J. (1980). “A family of embedded Runge-Kutta formulae.” Journal of Computational and Applied Mathematics, 6(1), 19–26.
5. Söderlind, G. (2002). “Automatic control and adaptive time-stepping.” Numerical Algorithms, 31, 281–310.
6. SciPy solve_ivp documentation: scipy.integrate.solve_ivp

Next Steps

• You now have the tools to solve any ODE numerically. In Module 2.1, we extend these ideas to partial differential equations.
• The adaptive step-size agent concept reappears in Module 3.4, where we replace hand-crafted policies with learned RL policies.