By the end of this week you should be able to:
write a system of ODEs and a higher-order ODE as a first-order system in vector form;
derive Euler’s method from Taylor expansion and from numerical integration;
implement Euler’s method and apply it to scalar and vector IVPs;
compare the accuracy of Euler’s and Heun’s methods and explain why Heun’s method converges faster.
In earlier weeks we solved ODEs analytically. Many real-world models, however, cannot be solved in closed form, and we must turn to numerical methods.
We focus on initial value problems (IVPs): given a function \(f\), an initial time \(t_0\) and an initial value \(y_0\), find \(y(t)\) for \(t > t_0\) satisfying \[\begin{equation} y'(t) = f(t,\, y(t)), \qquad y(t_0) = y_0. \label{eq:scalar-ivp} \end{equation}\]
Example: Exponential growth/decay
The IVP \(y' = \lambda y\), \(y(t_0) = y_0\) has the exact solution \(y(t) = y_0\,e^{\lambda(t - t_0)}\). For \(\lambda > 0\) this models population growth with unlimited resources; for \(\lambda < 0\) it models radioactive decay. We will use this as a test problem throughout, since the exact solution is known.
Example: Time-dependent coefficient
The IVP \(y' = -2t\,y\), \(y(0) = y_0\) has the exact solution \(y(t) = y_0\,e^{-t^2}\). Here the right-hand side \(f(t,y) = -2ty\) depends explicitly on \(t\).
Many models involve more than one unknown function evolving simultaneously. A system of \(m\) first-order ODEs takes the form \[\begin{equation} \mathbf{y}'(t) = \mathbf{f}(t,\, \mathbf{y}(t)), \qquad \mathbf{y}(t_0) = \mathbf{y}_0, \label{eq:system-ivp} \end{equation}\] where \(\mathbf{y}(t) = (y_1(t),\, y_2(t),\, \ldots,\, y_m(t))^\top\) is a vector in \(\mathbb{R}^m\) and \(\mathbf{f}\) maps \((t, \mathbf{y})\) to \(\mathbb{R}^m\). Written out component by component: \[\begin{aligned} y_1' &= f_1(t,\, y_1,\, y_2,\, \ldots,\, y_m), &\qquad y_1(t_0) &= y_{1,0}, \\ y_2' &= f_2(t,\, y_1,\, y_2,\, \ldots,\, y_m), &\qquad y_2(t_0) &= y_{2,0}, \\ &\;\;\vdots & &\;\;\vdots \\ y_m' &= f_m(t,\, y_1,\, y_2,\, \ldots,\, y_m), &\qquad y_m(t_0) &= y_{m,0}. \end{aligned}\] The numerical methods we develop for the scalar case [eq:scalar-ivp] apply directly to systems — we simply replace scalars with vectors.
Example: Lotka–Volterra equations
The interaction between prey (\(y\)) and predators (\(z\)) can be modelled by \[\begin{aligned} y'(t) &= \alpha\,y - \beta\,y\,z, \\ z'(t) &= \delta\,y\,z - \gamma\,z, \end{aligned}\] where \(\alpha\), \(\beta\), \(\delta\), \(\gamma > 0\) are parameters. In vector form, \(\mathbf{y} = (y, z)^\top\) and \[\mathbf{f}(t, \mathbf{y}) = \begin{pmatrix} \alpha y - \beta y z \\ \delta y z - \gamma z \end{pmatrix}.\] This system is nonlinear (the terms \(yz\) couple the two equations) and has no closed-form solution in general.
Example: SIR model
A simple model for the spread of an infectious disease divides the population into three classes: susceptible (\(S\)), infected (\(I\)) and removed (\(R\)). The ODE system is \[S' = -\beta S I, \qquad I' = \beta S I - \gamma I, \qquad R' = \gamma I,\] where \(\beta\) is the infection rate and \(\gamma\) is the removal rate. Since \(S + I + R\) is constant (check by adding the three equations), the system is effectively two-dimensional.
An \(m\)-th order ODE \[u^{(m)} = f\!\left(t,\, u,\, u',\, \ldots,\, u^{(m-1)}\right), \qquad u(t_0) = u_0,\; u'(t_0) = u_0',\; \ldots,\; u^{(m-1)}(t_0) = u_0^{(m-1)}\] can always be rewritten as a system of \(m\) first-order ODEs by introducing \[y_1 = u, \quad y_2 = u', \quad y_3 = u'', \quad \ldots, \quad y_m = u^{(m-1)}.\] Then \[\begin{aligned} y_1' &= y_2, &\quad y_1(t_0) &= u_0, \\ y_2' &= y_3, &\quad y_2(t_0) &= u_0', \\ &\;\;\vdots & &\;\;\vdots \\ y_{m-1}' &= y_m, &\quad y_{m-1}(t_0) &= u_0^{(m-2)}, \\ y_m' &= f(t,\, y_1,\, y_2,\, \ldots,\, y_m), &\quad y_m(t_0) &= u_0^{(m-1)}. \end{aligned}\] This is now a first-order system of the form [eq:system-ivp] and can be solved numerically using any of the methods in this chapter.
Example: Van der Pol oscillator
The Van der Pol equation is the second-order ODE \[u'' = \mu(1 - u^2)\,u' - u, \qquad u(0) = u_0, \quad u'(0) = u_0'.\] Setting \(y_1 = u\) and \(y_2 = u'\), this becomes the first-order system \[\begin{aligned} y_1' &= y_2, &\quad y_1(0) &= u_0, \\ y_2' &= \mu(1 - y_1^2)\,y_2 - y_1, &\quad y_2(0) &= u_0'. \end{aligned}\] For \(\mu > 0\) the solutions exhibit self-sustaining oscillations (a limit cycle). Common test values are \(\mu = 2\), \(u_0 = 2\), \(u_0' = 0\).
Key Idea
Any ODE or system of ODEs, of any order, can be reduced to a first-order system \(\mathbf{y}' = \mathbf{f}(t, \mathbf{y})\). We therefore only need to develop numerical methods for first-order systems as they cover all cases.
We now turn to our first numerical method. Given the IVP [eq:scalar-ivp] and a final time \(T\), we want to compute an approximation to \(y(t)\) on \([t_0, T]\).
Starting from \(t_0\), choose a step size \(\tau > 0\) and set \(t_1 = t_0 + \tau\). Taylor-expanding the exact solution about \(t_0\): \[y(t_0 + \tau) = y(t_0) + \tau\,y'(t_0) + \tfrac{1}{2}\tau^2\,y''(t_0) + \cdots\] For small \(\tau\), drop everything beyond the first two terms and use the ODE to replace \(y'(t_0) = f(t_0, y_0)\): \[y(t_0 + \tau) \approx y_0 + \tau\,f(t_0, y_0).\] This gives the numerical approximation \(y_1 \approx y(t_1)\). Repeating from each new point produces a sequence of approximations.
Key Idea
Given \(f(t,y)\), initial value \((t_0, y_0)\), final time \(T\) and step size \(\tau_k = t_{k+1} - t_k\): \[\boxed{y_{k+1} = y_k + \tau_k\,f(t_k,\, y_k), \qquad t_{k+1} = t_k + \tau_k.}\] Repeat until \(t_k \ge T\). For a constant step size with \(N\) steps, \(\tau = (T - t_0)/N\).
Euler’s method is an example of a one-step method (OSM): to advance from \(t_k\) to \(t_{k+1}\), it only needs the current value \(y_k\), not any earlier values.
Between the nodes, the numerical solution can be extended by connecting each pair \((t_k, y_k)\) and \((t_{k+1}, y_{k+1})\) with a straight line — this is consistent with the linear Taylor approximation underlying the method.
Rearranging Euler’s formula: \[\frac{y_{k+1} - y_k}{\tau} = f(t_k, y_k).\] This replaces the derivative \(y'(t_k)\) with a forward difference quotient, the simplest finite-difference approximation. We will discuss finite difference approaches in more detail in year 2 of your degree.
In Python, Euler’s method can be implemented in a few lines:
def explicit_euler(y0, t0, T, f, Nmax):
ys = [y0]
ts = [t0]
dt = (T - t0) / Nmax
while ts[-1] < T:
t, y = ts[-1], ys[-1]
ys.append(y + dt * f(t, y))
ts.append(t + dt)
return (np.array(ts), np.array(ys))This works for both scalar and vector problems, if y0 is
a NumPy array, the same code handles systems automatically.
Example: Testing Euler on exponential growth
For \(y' = \lambda y\), \(y(0) = 1\) with \(\lambda = 1\) on \([0, 1]\), using \(N = 4\) steps gives \(\tau = 0.25\). The numerical solution is computed step by step: \[\begin{aligned} y_0 &= 1, \\ y_1 &= y_0 + 0.25 \times 1 \times y_0 = 1.25, \\ y_2 &= y_1 + 0.25 \times 1 \times y_1 = 1.5625, \\ y_3 &= y_2 + 0.25 \times 1 \times y_2 \approx 1.9531, \\ y_4 &= y_3 + 0.25 \times 1 \times y_3 \approx 2.4414. \end{aligned}\] The exact solution at \(t = 1\) is \(y(1) = e \approx 2.7183\), so the error is about \(0.277\). As we increase \(N\), the approximation improves.
How does the error behave as we refine the step size? Define the global error at step \(k\) as \(e_k = |y(t_k) - y_k|\), and the maximum error over all steps as \[E(\tau) = \max_{0 \le k \le N} |y(t_k) - y_k|.\]
For Euler’s method with constant step size \(\tau\), one can show that \(E(\tau) = \mathcal{O}(\tau)\) — the method is first-order accurate. In practice, this means:
Key Idea: Convergence of Euler’s method
If you halve the step size (double \(N\)), the error is approximately halved. The error reduction factor is 2.
We can verify this numerically by computing \(E(\tau)\) for \(N = 4, 8, 16, 32, 64, 128\) and checking the ratio of successive errors.
Euler’s method approximates the solution by taking a single step using the slope at the start of the interval. But the slope changes as we move from \(t_k\) to \(t_{k+1}\) — using only the starting slope systematically misses this change, which is why the method is only first-order accurate.
Heun’s method improves on this with a simple idea: use the average of two slopes.
Step 1 — predict.Take a trial Euler step to estimate where we’ll end up: \[\tilde{y}_{k+1} = y_k + \tau_k\,f(t_k,\, y_k).\]
Step 2 — evaluate.Compute the slope at this predicted point: \[k_1 = f(t_k,\, y_k), \qquad k_2 = f(t_k + \tau_k,\, \tilde{y}_{k+1}).\] Now \(k_1\) is the slope at the start and \(k_2\) is an estimate of the slope at the end.
Step 3 — correct.Take the actual step using the average of the two slopes: \[\boxed{y_{k+1} = y_k + \frac{\tau_k}{2}(k_1 + k_2).}\]
This “predictor-corrector” strategy costs one extra evaluation of \(f\) per step, but as we shall see pays off significantly in terms of the overall accuracy of the method.
Key Idea
Heun’s method: Given \(f(t,y)\), initial value \((t_0, y_0)\), final time \(T\) and step size \(\tau_k\): \[\boxed{ \begin{aligned} k_1 &= f(t_k,\, y_k), \\ k_2 &= f(t_k + \tau_k,\, y_k + \tau_k\,k_1), \\ y_{k+1} &= y_k + \frac{\tau_k}{2}(k_1 + k_2). \end{aligned} }\]
Heun’s method is also called the improved Euler method, it requires two evaluations of \(f\) per step (compared to one for Euler), but gains significantly in accuracy.
def heun(y0, t0, T, f, Nmax):
ys = [y0]
ts = [t0]
dt = (T - t0) / Nmax
while ts[-1] < T:
t, y = ts[-1], ys[-1]
k1 = f(t, y)
k2 = f(t + dt, y + dt * k1)
ys.append(y + 0.5 * dt * (k1 + k2))
ts.append(t + dt)
return (np.array(ts), np.array(ys))For Heun’s method with constant step size \(\tau\), the global error satisfies \(E(\tau) = \mathcal{O}(\tau^2)\) — the method is second-order accurate.
Key Idea
If you halve the step size (double \(N\)), the error is approximately reduced by a factor of 4. Heun’s method converges much faster than Euler’s method, at the cost of one extra function evaluation per step.
The following table summarises the comparison:
| Method | Order | Error reduction when \(\tau \to \tau/2\) |
|---|---|---|
| Euler | 1 | \(\times\, 2\) |
| Heun | 2 | \(\times\, 4\) |
In practice, for a given accuracy target, Heun’s method typically requires far fewer steps than Euler’s method, even accounting for the extra function evaluation.
Example: Predator–prey dynamics
Consider the Lotka–Volterra system with parameters \(\alpha = 2\), \(\beta = 1\), \(\delta = 0.5\), \(\gamma = 1\) and initial conditions \(y(0) = 2\), \(z(0) = 0.5\), solved over \([0, 20]\).
The right-hand side in Python:
def lotka_volterra(t, y):
alpha, beta, delta, gamma = 2, 1, 0.5, 1
dy = np.array([alpha*y[0] - beta*y[0]*y[1],
delta*y[0]*y[1] - gamma*y[1]])
return dyThe exact solution is not known in closed form, but the solutions are known to be periodic and positive. This provides a qualitative check: if the numerical solution drifts or goes negative, the step size is too large.
With Euler’s method and \(\tau = 0.1\), the solution visibly drifts (the amplitude grows over time). Reducing to \(\tau = 0.02\) or \(\tau = 0.002\) improves matters. Heun’s method gives much better results for the same step size.
Example: Limit cycle oscillation
The Van der Pol equation with \(\mu = 2\), \(u(0) = 2\), \(u'(0) = 0\) is converted to the first-order system \[y_1' = y_2, \qquad y_2' = 2(1 - y_1^2)\,y_2 - y_1.\] In Python:
def van_der_pol(t, y):
mu = 2
dy = np.array([y[1],
mu*(1 - y[0]**2)*y[1] - y[0]])
return dySolving over \([0, 20]\) with Euler’s method and \(\tau = 0.1\) produces visible inaccuracies, particularly in the steep parts of the oscillation. Heun’s method with the same step size is significantly more accurate. Experimenting with different values of \(\mu\) reveals that larger \(\mu\) creates sharper transitions, requiring smaller step sizes.
Exercise 6.1
Consider the IVP \(y' = \lambda y\), \(y(0) = 1\) with \(\lambda = 1\) on \([0, 1]\).
Compute the Euler approximation by hand for \(N = 4\) steps and compare with \(y(1) = e\).
Write code to compute the maximum error \(E(\tau) = \max_k |y(t_k) - y_k|\) for \(N = 4, 8, 16, 32, 64, 128\). Verify that the error halves each time \(N\) is doubled.
Repeat part (b) using Heun’s method. Verify that the error is reduced by a factor of 4 each time \(N\) is doubled.
Exercise 6.2
Solve the Lotka–Volterra system from the worked example using both Euler’s and Heun’s methods.
Plot the prey and predator populations against time for \(\tau = 0.1\), \(0.02\) and \(0.002\) using Euler’s method. For which step sizes does the solution remain periodic and positive?
Repeat with Heun’s method. How does the required step size compare?
Plot \(z\) against \(y\) (a phase portrait). The exact orbits are closed loops. Which method preserves this structure better?
Exercise 6.3
The Van der Pol oscillator with \(\mu = 2\), \(u(0) = 2\), \(u'(0) = 0\) is solved over \([0, 20]\).
Write down the conversion to a first-order system explicitly.
Solve using Euler’s method with \(\tau = 0.1\). Plot \(y_1(t)\) and \(y_2(t)\).
Solve using Heun’s method with \(\tau = 0.1\) and compare. How many Euler steps would you need to obtain a similar quality solution?
Experiment with \(\mu = 0.1\), \(1\), \(5\) and \(10\). How does the character of the solution change? Which values of \(\mu\) require smaller step sizes, and why?