54 Numerical ODEs and PDEs

Stepping through time and space

54.1 The setup: replacing dy/dx with Δy/Δx

The initial-value problem has the form

\[ y' = f(x, y), \quad y(x_0) = y_0. \]

Read this as: “the slope of \(y\) at any point \((x, y)\) is given by the function \(f\), and we know where we start.” The question is: what is \(y\) at some later \(x\)?

If \(f\) is simple enough, calculus gives you \(y(x)\) directly. When it isn’t, you fall back on the most basic approximation: the tangent-line step.

The step size \(h\) — a small positive number — is your control parameter. Smaller \(h\) means more accuracy and more work. The methods below differ in how many function evaluations they use per step, and therefore in how quickly the error shrinks as \(h\) decreases.

54.2 Euler’s method

54.2.1 Forward Euler

Approximate \(y'(x_n) \approx \dfrac{y_{n+1} - y_n}{h}\), rearrange, and you get

\[ y_{n+1} = y_n + h \, f(x_n, y_n). \]

\(y_n\) means “the approximate value of \(y\) at step \(n\).” The update rule says: take the current value, add the step size times the slope at the current point.

Local truncation error. Taylor-expand \(y(x_{n+1})\) about \(x_n\):

\[ y(x_{n+1}) = y(x_n) + h y'(x_n) + \tfrac{h^2}{2} y''(x_n) + \cdots \]

Forward Euler keeps only the first-derivative term, so the error in one step is \(O(h^2)\). Over \(N = (b - a)/h\) steps, the errors accumulate: global error is \(O(h)\). Halve the step size, halve the global error.

54.2.2 Backward (implicit) Euler

\[ y_{n+1} = y_n + h \, f(x_{n+1}, y_{n+1}). \]

The slope is evaluated at the new point, which makes this an equation you must solve for \(y_{n+1}\) at each step — more expensive per step. The payoff is unconditional stability for stiff equations (equations where the true solution decays quickly and explicit methods require an extremely small \(h\) to stay bounded). For non-stiff problems, forward Euler is fine.

What “stiff” means. A stiff ODE has solution components that decay at very different rates. The fast-decaying component forces \(h\) to be tiny for forward Euler to remain stable, even though that component itself becomes negligible quickly. Implicit methods remain stable at larger \(h\) because they incorporate information from the next step.

54.2.3 Worked example: forward Euler on \(y' = -2y\)

Exact solution: \(y = e^{-2x}\). Let \(h = 0.5\), \(y_0 = 1\).

\(n\)	\(x_n\)	\(y_n\)	\(f = -2y_n\)	\(y_{n+1} = y_n + h f\)
0	0.0	1.0000	−2.0000	1 + 0.5(−2) = 0.0000

That’s the stability boundary. For an equation of the form \(y' = \lambda y\) — where \(\lambda\) is the coefficient multiplying \(y\), here \(\lambda = -2\) — forward Euler multiplies the solution by \((1 + h\lambda)\) at each step. Stability requires \(|1 + h\lambda| \leq 1\), which for \(\lambda = -2\) gives \(h \leq 1\). With \(h = 0.5\) we land exactly at the edge: the solution collapses to zero and stays there. For \(h > 1\) the factor \((1 + h\lambda)\) exceeds 1 in magnitude, and the computed solution oscillates with growing amplitude — blowing up — even though the true solution decays smoothly to zero. Use \(h = 0.2\) to see the method tracking the decay:

\(n\)	\(x_n\)	\(y_n\)	\(f = -2y_n\)	\(y_{n+1}\)	Exact \(e^{-2x}\)
0	0.0	1.0000	−2.0000	0.6000	—
1	0.2	0.6000	−1.2000	0.3600	0.6703
2	0.4	0.3600	−0.7200	0.2160	0.4493

The error at \(x = 0.4\) is about \(0.089\), or roughly 20%. Three steps with \(h = 0.2\) on a decaying exponential is not a recipe for accuracy — but it illustrates the mechanics.

54.3 Runge-Kutta methods

The key idea: use several evaluations of \(f\) within a single step — at different intermediate points — then combine them into a weighted average slope. More work per step, but the error drops much faster with decreasing \(h\).

54.3.1 RK2: two-slope methods

Midpoint method. Take a half-step to estimate the midpoint, then use the slope there:

\[ k_1 = f(x_n, y_n), \qquad k_2 = f\!\left(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2} k_1\right) \]

\[ y_{n+1} = y_n + h k_2. \]

Heun’s method (modified Euler). Use the slope at the beginning and the predicted end, average them:

\[ k_1 = f(x_n, y_n), \qquad k_2 = f(x_n + h,\; y_n + h k_1) \]

\[ y_{n+1} = y_n + \tfrac{h}{2}(k_1 + k_2). \]

Both are \(O(h^2)\) per step, giving global error \(O(h^2)\). Halving \(h\) divides the error by 4.

54.3.2 RK4: the standard workhorse

Evaluate four slopes — \(k_1\), \(k_2\), \(k_3\), \(k_4\) — and combine with weights \(1, 2, 2, 1\):

\[ k_1 = f(x_n,\; y_n) \]

\[ k_2 = f\!\left(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2} k_1\right) \]

\[ k_3 = f\!\left(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2} k_2\right) \]

\[ k_4 = f(x_n + h,\; y_n + h k_3) \]

\[ y_{n+1} = y_n + \tfrac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4). \]

\(k_1\) is the slope at the left end. \(k_2\) and \(k_3\) are two estimates of the midpoint slope (the second corrects the first). \(k_4\) is the slope at the right end. The weights \(1, 2, 2, 1\) come from Simpson’s rule: the interval has three evaluation points — left end, midpoint, right end — with Simpson weights \(1, 4, 1\). RK4 has two independent midpoint estimates rather than one, so their weight is split, giving \(1, 2, 2, 1\) after normalising to sum to \(6\).

Error. Local error is \(O(h^5)\); global error is \(O(h^4)\). Halving \(h\) divides the global error by 16. This is why RK4 is the default choice for non-stiff problems: it is usually accurate enough with step sizes that are not punishingly small.

Adaptive step size. RK4 with a fixed \(h\) wastes work in regions where the solution is smooth and underestimates error in regions where it changes rapidly. Adaptive methods (e.g. RK45, the Dormand-Prince method) compute two approximations of different orders at each step, compare them to estimate the local error, and adjust \(h\) to keep the error within a specified tolerance. Embedded pairs like RK23 and RK45 reuse slope evaluations between the two orders, keeping the extra cost small. The implementation is routine; the idea is simply: use \(h\) proportional to what the solution needs.

54.3.3 Worked example: RK4 on \(y' = y\), one step

\(y(0) = 1\), \(h = 1\). Exact: \(y(1) = e \approx 2.71828\).

\[ k_1 = f(0, 1) = 1 \] \[ k_2 = f(0.5,\; 1 + 0.5) = f(0.5, 1.5) = 1.5 \] \[ k_3 = f(0.5,\; 1 + 0.5 \times 1.5) = f(0.5, 1.75) = 1.75 \] \[ k_4 = f(1,\; 1 + 1 \times 1.75) = f(1, 2.75) = 2.75 \] \[ y_1 = 1 + \tfrac{1}{6}(1 + 3 + 3.5 + 2.75) = 1 + \tfrac{10.25}{6} \approx 2.7083. \]

Error: \(|2.7083 - 2.71828| \approx 0.00998\). With a single step of \(h = 1\), RK4 gets the answer to 3 significant figures. Forward Euler with \(h = 1\) gives \(y_1 = 2.0\) — an error of 26%.

54.4 Systems of ODEs

Many physical problems involve second- (or higher-) order ODEs. A second-order ODE can always be rewritten as a system of two first-order ODEs, and then Euler or RK4 applies to the system directly.

Reduction procedure. Given

\[ y'' + p(x)y' + q(x)y = r(x), \]

introduce \(u_1 = y\) and \(u_2 = y'\). Then:

\[ u_1' = u_2 \]

\[ u_2' = r(x) - p(x) u_2 - q(x) u_1. \]

Write this as \(\mathbf{u}' = \mathbf{F}(x, \mathbf{u})\) where \(\mathbf{u} = [u_1, u_2]^T\). Apply Euler or RK4 to the vector equation: each \(k_i\) becomes a two-component vector.

54.4.1 Worked example: damped oscillator

\(y'' + 0.5y' + 4y = 0\), \(y(0) = 1\), \(y'(0) = 0\), \(h = 0.1\).

Let \(u_1 = y\), \(u_2 = y'\). Then:

\[ u_1' = u_2, \qquad u_2' = -4 u_1 - 0.5 u_2. \]

One Euler step:

\[ k_1^{(1)} = u_2^{(0)} = 0, \qquad k_1^{(2)} = -4(1) - 0.5(0) = -4. \]

\[ u_1^{(1)} = 1 + 0.1 \times 0 = 1.0, \qquad u_2^{(1)} = 0 + 0.1 \times (-4) = -0.4. \]

So after one step: \(y(0.1) \approx 1.0\), \(y'(0.1) \approx -0.4\). The method is mechanical once the system is written down; the reduction step is the only place judgment is required.

54.5 Finite differences for PDEs

The same idea — replace a derivative with a finite ratio — extends to PDEs, where you discretise both time and space.

54.5.1 The 1D heat equation

\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad 0 \le x \le L, \quad t \ge 0, \]

with given initial condition \(u(x, 0)\) and boundary conditions. \(\alpha\) is thermal diffusivity (a positive constant).

Lay a grid: \(x_j = j \, \Delta x\) for \(j = 0, 1, \ldots, M\) and \(t^n = n \, \Delta t\) for \(n = 0, 1, \ldots\). Write \(u_j^n\) for the approximation to \(u(x_j, t^n)\).

FTCS scheme (forward-time, central-space). Use a forward difference in time and a central difference in space:

\[ \frac{u_j^{n+1} - u_j^n}{\Delta t} = \alpha \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{(\Delta x)^2}. \]

Define the Fourier number (also called the mesh ratio)

\[ r = \frac{\alpha \, \Delta t}{(\Delta x)^2}. \]

Rearrange:

\[ u_j^{n+1} = r \, u_{j-1}^n + (1 - 2r) \, u_j^n + r \, u_{j+1}^n. \]

This is explicit: each new value depends only on known values at the previous time level.

Stability condition. Look at the FTCS update coefficients: \(r\), \((1-2r)\), \(r\). When \(r \leq 1/2\), the coefficient \((1-2r)\) is non-negative, so the new value at each interior node is a weighted average of the neighbouring values — physically reasonable, bounded, and stable. When \(r > 1/2\), the coefficient \((1-2r)\) turns negative. The centre node’s own value now enters with a negative weight: the update is no longer an average. A small perturbation at one node gets amplified rather than smoothed; each time step multiplies the perturbation, and the solution oscillates with growing amplitude until it overflows. The physical heat equation diffuses and decays; the numerical scheme with an oversized time step does the opposite — it invents energy out of rounding errors. Von Neumann stability analysis makes this precise: \(r \le \tfrac{1}{2}\) is the condition. In practice, for \(\alpha = 1\), \(\Delta x = 0.1\): the maximum stable \(\Delta t\) is \(\tfrac{1}{2} (0.1)^2 = 0.005\).

54.5.2 Crank-Nicolson: unconditional stability

Average the spatial derivative between time levels \(n\) and \(n+1\):

\[ \frac{u_j^{n+1} - u_j^n}{\Delta t} = \frac{\alpha}{2} \left[ \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{(\Delta x)^2} + \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} \right]. \]

The unknowns at time level \(n+1\) appear on both sides. Collecting them gives a tridiagonal linear system to solve at each time step — implicit, but only marginally more expensive than FTCS because tridiagonal systems have a special structure: each equation involves only its immediate neighbours, so a forward sweep followed by back-substitution solves the whole system in \(O(N)\) operations. This procedure is called the Thomas algorithm — Gaussian elimination stripped down to exploit the tridiagonal sparsity pattern. The method is unconditionally stable: any \(r\) is admissible, so you can take large time steps without blowup. The global error is \(O((\Delta t)^2 + (\Delta x)^2)\) — second order in both.

54.5.3 Boundary-value problems: tridiagonal systems from central differences

For a 1D BVP:

\[ y'' = f(x), \quad y(a) = \alpha, \quad y(b) = \beta, \]

discretise the interior points \(x_1, \ldots, x_{N-1}\) with spacing \(h = (b-a)/N\). The central-difference approximation \(y'' \approx (y_{i-1} - 2y_i + y_{i+1})/h^2\) at each interior point gives:

\[ y_{i-1} - 2y_i + y_{i+1} = h^2 f(x_i), \quad i = 1, \ldots, N-1. \]

Collect these equations (incorporating \(y_0 = \alpha\) and \(y_N = \beta\) into the first and last equations) and you get a tridiagonal linear system. The structure is

\[ \begin{pmatrix} -2 & 1 & & \\ 1 & -2 & 1 & \\ & \ddots & \ddots & 1 \\ & & 1 & -2 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{pmatrix} = h^2 \begin{pmatrix} f(x_1) - \alpha/h^2 \\ f(x_2) \\ \vdots \\ f(x_{N-1}) - \beta/h^2 \end{pmatrix}. \]

Solve with the Thomas algorithm (forward sweep + back substitution) in \(O(N)\) operations.

Why the tridiagonal structure matters. Gaussian elimination on a dense \(N \times N\) system costs \(O(N^3)\). Exploiting the tridiagonal structure reduces this to \(O(N)\). For \(N = 1000\) spatial nodes, that is a factor of \(10^9\) in computational work. Recognising sparsity is not a nicety; it is what makes large-scale simulation possible.

54.6 Worked example: FTCS one step

Heat equation on \([0, 1]\) with 5 spatial nodes (\(\Delta x = 0.25\)), \(\alpha = 1\), \(\Delta t = 0.02\).

First check stability: \(r = 1 \times 0.02 / (0.25)^2 = 0.02 / 0.0625 = 0.32 < 0.5\). Stable.

Initial condition: \(u(x, 0) = \sin(\pi x)\). Interior nodes: \(x_1 = 0.25\), \(x_2 = 0.50\), \(x_3 = 0.75\).

Values at \(t = 0\): \(u_0^0 = 0\), \(u_1^0 = \sin(0.25\pi) \approx 0.7071\), \(u_2^0 = \sin(0.5\pi) = 1.0\), \(u_3^0 = \sin(0.75\pi) \approx 0.7071\), \(u_4^0 = 0\).

Apply \(u_j^1 = 0.32 u_{j-1}^0 + 0.36 u_j^0 + 0.32 u_{j+1}^0\) (coefficients: \(r = 0.32\), \(1 - 2r = 0.36\)):

\[ u_1^1 = 0.32(0) + 0.36(0.7071) + 0.32(1.0) = 0.2546 + 0.3200 = 0.5746 \]

\[ u_2^1 = 0.32(0.7071) + 0.36(1.0) + 0.32(0.7071) = 0.2263 + 0.3600 + 0.2263 = 0.8126 \]

\[ u_3^1 = 0.5746 \quad \text{(by symmetry)} \]

Exact solution: \(u(x, t) = e^{-\pi^2 t} \sin(\pi x)\). At \(t = 0.02\): \(e^{-0.197} \approx 0.821\). So \(u(0.5, 0.02)^{\text{exact}} = 0.821\). Numerical: \(0.8126\). Error \(\approx 0.008\).

54.7 Where this goes

The methods here are the direct input to adaptive ODE solvers (e.g. scipy.integrate.solve_ivp, MATLAB’s ode45). Those solvers use embedded Runge-Kutta pairs — the same structure as above — with automatic step-size control. Understanding the underlying stepping scheme lets you diagnose solver failures: a stiff ODE switched to an explicit method, or a step-size tolerance too loose for the problem’s dynamics.

Finite differences for PDEs extend naturally to finite element methods (FEM), where the domain is discretised into irregular elements rather than a regular grid. FEM is the foundation of structural analysis, fluid simulation, and electromagnetic modelling. The tridiagonal system from 1D BVPs generalises to sparse banded systems in higher dimensions — still tractable, still solved by exploiting the sparsity pattern.

54.8 Applications

Where numerical ODE/PDE methods appear

Structural engineering. Finite element codes (ANSYS, Abaqus) step through nonlinear load-displacement ODEs and solve elliptic PDEs for stress fields.
Aerospace. Flight simulators use RK4 (or RK45) to integrate the 6-DOF equations of motion. Autopilot design uses frequency-domain tools (Laplace, Bode analysis), but the plant model is validated by time-domain simulation using exactly these stepping methods.
Climate science. Atmospheric models discretise the primitive equations on a spherical grid with time steps of minutes; finite-difference or spectral schemes advance the state.
Biomedical. Hodgkin-Huxley neuron models are stiff ODEs solved with implicit methods; cardiac electrophysiology models are stiff PDE systems.
Computer graphics. Game physics engines use symplectic Euler (a variant of forward Euler that conserves energy better) to integrate rigid-body dynamics at 60–120 Hz.
Finance. The Black-Scholes PDE for option pricing is discretised with Crank-Nicolson on a log-price grid.

54.9 Exercises

Exercise 1. A bacterial culture’s population decreases due to an antibiotic. The concentration \(y\) (mg/L) satisfies \(dy/dx = -2y\), \(y(0) = 1\). Apply forward Euler with \(h = 0.1\) for three steps. Compare with the exact values \(e^{-0.2}\), \(e^{-0.4}\), \(e^{-0.6}\).

Exercise 2. Apply RK4 to \(dy/dx = -2y\), \(y(0) = 1\) with a single step of \(h = 0.3\). Compute all four slopes \(k_1, k_2, k_3, k_4\) and the resulting \(y_1\). Compare with the exact \(e^{-0.6}\).

Exercise 3. A second-order system satisfies \(y'' - 3y' + 2y = 0\), \(y(0) = 0\), \(y'(0) = 1\), \(h = 0.1\).

Reduce to a first-order system.
Apply one step of forward Euler to get \(y(0.1)\) and \(y'(0.1)\).

Exercise 4. A 1D diffusion problem is modelled with diffusivity \(\alpha = 0.5\) (dimensionless or in consistent units scaled to the rod length). The rod lies on \([0, 1]\) with 5 spatial nodes (\(\Delta x = 0.25\)). Initial condition: \(u(x,0) = \sin(\pi x)\). Boundary conditions: \(u(0,t) = u(1,t) = 0\). Take \(\Delta t = 0.01\).

Apply one FTCS time step. Give the value of \(u\) at each interior node.

Exercise 5. Using the parameters from Exercise 4, verify analytically that the scheme is stable. Now choose \(\Delta t = 0.04\) (with the same \(\Delta x = 0.25\), \(\alpha = 0.5\)). Compute \(r\) and explain what will happen if you run the scheme forward in time.

Exercise 6. Solve the boundary-value problem \(y'' = -\pi^2 \sin(\pi x)\), \(y(0) = 0\), \(y(1) = 0\) numerically using central differences with 4 interior nodes (\(h = 0.2\)). Set up the \(4 \times 4\) linear system and solve it. The exact solution is \(y(x) = \sin(\pi x)\).