58  Numerical ODEs and PDEs

Stepping through time and space

58.1 What this chapter helps you do

Symbols to keep handy

These are the bits of notation you'll see a lot. If a line of symbols feels like a fence, read it out loud once, then keep going.

  • r = t / (x)^2: the Fourier number — the ratio that controls stability of the heat equation scheme

  • k_1, k_2, k_3, k_4: the four slope estimates in one RK4 step

Definitions to keep handy

These are the words we keep coming back to. If one feels slippery, come back here and steady it before you push on.

  • time stepping: A way to solve an ODE by moving forward in small steps instead of finding a closed-form formula.

  • explicit vs implicit: Explicit uses only known information; implicit solves an equation each step but can be more stable.

  • local vs global error: Local is the error made in one step; global is the accumulated error after many steps.

  • stiffness: When a system has very fast and very slow scales together, forcing tiny steps unless you use a stable method.

  • finite difference: Replacing derivatives with difference quotients on a grid so a PDE becomes algebra.

Here is the main move this chapter is making, in plain terms. You do not need to be fast. You just need to keep the thread.

  • Coming in: You can write down the differential equation. Closed-form solutions exist only for special cases; real engineering problems are rarely special.

  • Leaving with: Numerical time-stepping and finite-difference discretisation turn a continuous problem into an algebraic one a computer can solve exactly.

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

When you see an initial-value problem, you are being asked for a whole function, not a single number. You know how the system changes (the slope), and you know where it starts. The problem is: if we start here, where do we end up after we move forward?

The standard form is:

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

Read this as

  • y' = f(x,y): “the slope of the curve at the point (x,y) is given by the rule f”
  • y(x_0)=y_0: “we start at the point (x_0,y_0)”

If f is simple enough, calculus gives you y(x) directly. When it isn’t, we do the most human thing you can do: take a small step using the slope you have right now, then repeat.

The step size h is how far you move each time. Smaller h means you follow the curve more closely (better accuracy) but you have to take more steps (more computation). Different methods spend more work per step in order to reduce the error faster.

58.3 Euler’s method

58.3.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.

58.3.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.

58.3.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.

Use the step-size slider to compare how Euler and RK4 accumulate error on the same simple growth problem.

58.4 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.

58.4.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.

58.4.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.

58.4.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%.

58.5 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.

58.5.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.

58.6 Finite differences for PDEs

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

58.6.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.

Use the controls to see the FTCS heat solver behave like a smoother when the stability condition holds, and like a source of fake oscillations when it does not.

58.6.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.

58.6.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.

58.7 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.

58.8 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.

58.9 Applications

TipWhere 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.

58.10 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.

  1. Reduce to a first-order system.
  2. 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 you expect to happen if you run the scheme forward in time. Then, as a separate check, choose a larger \Delta t that makes r > 1/2 (for these numbers, \Delta t = 0.07 works) and describe the instability you would expect.

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).