Notes on Numerical Integration of ODEs

Most simulations of physical systems reduce to the same task: advance an ordinary differential equation $\dot{y} = f(t, y)$ in small time steps of size $h$. The choice of method decides how much accuracy a given step size buys, and whether the result stays bounded at all. These notes compare two standard methods on the simplest useful test problem.

The test equation

Consider $\dot{y} = \lambda y$ with $y(0) = 1$, whose exact solution is $y(t) = e^{\lambda t}$. Setting $\lambda = -1$ gives a decaying signal that is easy to check against.

Explicit Euler takes one slope evaluation per step:

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

On the test equation this becomes $y_{n+1} = (1 + h\lambda)\,y_n$. The error made in one step is of order $h^2$, and accumulated over a fixed time interval the global error is of order $h$.

Classical Runge-Kutta (RK4) evaluates the slope four times per step:

$$\begin{aligned} k_1 &= f(t_n, y_n), \\ k_2 &= f\bigl(t_n + \tfrac{h}{2},\, y_n + \tfrac{h}{2}k_1\bigr), \\ k_3 &= f\bigl(t_n + \tfrac{h}{2},\, y_n + \tfrac{h}{2}k_2\bigr), \\ k_4 &= f\bigl(t_n + h,\, y_n + h\,k_3\bigr), \end{aligned}$$

and combines them into the next value:

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

Its global error is of order $h^4$.

Measuring the order

import math

def euler(f, t, y, h, n):
    for _ in range(n):
        y += h * f(t, y)
        t += h
    return y

def rk4(f, t, y, h, n):
    for _ in range(n):
        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 += h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        t += h
    return y

f = lambda t, y: -y
exact = math.exp(-1.0)
for h in (0.1, 0.05, 0.025):
    n = round(1.0 / h)
    print(h, abs(euler(f, 0.0, 1.0, h, n) - exact), abs(rk4(f, 0.0, 1.0, h, n) - exact))

Absolute error at $t = 1$, from running this script:

Step size $h$	Euler	RK4
0.1	$1.92 \times 10^{-2}$	$3.33 \times 10^{-7}$
0.05	$9.39 \times 10^{-3}$	$2.00 \times 10^{-8}$
0.025	$4.65 \times 10^{-3}$	$1.22 \times 10^{-9}$

Halving $h$ roughly halves the Euler error and divides the RK4 error by about sixteen, consistent with first-order and fourth-order convergence.

Accuracy is not stability

For Euler on the test equation with real $\lambda < 0$, the numerical solution decays only if $|1 + h\lambda| \le 1$, which means $h \le 2/|\lambda|$. With $\lambda = -50$ and $h = 0.05$ the factor is $1 - 2.5 = -1.5$, so the computed solution grows in magnitude and flips sign at every step, even though the true solution decays quickly.

Choosing a step size by asking how accurate the answer should be is only half the job. For fast decaying dynamics the step size may be limited by stability instead.

RK4 has a larger stability region than Euler, but it is also bounded. On the negative real axis it reaches about $-2.79$, so the step $h\lambda = -2.5$ that makes Euler blow up is stable for RK4.

Stiff problems, where fast and slow dynamics coexist, usually call for implicit methods, which are outside the scope of these notes.