10.3 Finite differences for the wave and heat equations

The PDE simulations in Foundations 6 — the wave/heat side-by-side, the heat diffusion, the Laplace solver — all run on the same underlying technique: replace the spatial and temporal derivatives in the PDE by finite differences on a $(x, t)$ grid. The resulting recurrence advances the field one time step at a time. The technique works, but with a sharp constraint: the time step and the grid spacing have to satisfy a stability condition that depends on the equation. Crossing that boundary makes the simulation blow up in a handful of steps.

This lesson sets up the schemes, derives the CFL stability condition for the wave equation, and explains the analogous (and tighter) constraint for the heat equation.

Discretising in space and time

Place a uniform grid on the spatial axis: $x_j = j \Delta x$ for $j = 0, 1, \ldots, N$ . Place a uniform grid on the time axis: $t_n = n \Delta t$ for $n = 0, 1, 2, \ldots$ . We approximate the unknown field $u(x, t)$ by its values at the grid points,

u^n_j \;\approx\; u(x_j, t_n).

The wave equation is $\partial_t^2 u = c^2 \partial_x^2 u$ . Apply the second-order centred difference (from 10.1) to both derivatives:

\frac{u^{n+1}_j - 2 u^n_j + u^{n-1}_j}{\Delta t^2} \;=\; c^2 \, \frac{u^n_{j+1} - 2 u^n_j + u^n_{j-1}}{\Delta x^2}.

Solve for the only unknown at the new time level, $u^{n+1}_j$ :

\boxed{\;u^{n+1}_j \;=\; 2 u^n_j - u^{n-1}_j + C^2 \bigl(u^n_{j+1} - 2 u^n_j + u^n_{j-1}\bigr),\;}

where $C \equiv c\, \Delta t / \Delta x$ is the Courant number. This is the leapfrog scheme for the wave equation. It needs the field at two preceding time levels to compute the next one — natural for a second-order-in-time equation, which needs two pieces of initial data.

The heat equation $\partial_t u = D \partial_x^2 u$ is first-order in time. Apply forward Euler in time (from 10.2) and centred difference in space:

\boxed{\;u^{n+1}_j \;=\; u^n_j + \frac{D \Delta t}{\Delta x^2}\, \bigl(u^n_{j+1} - 2 u^n_j + u^n_{j-1}\bigr).\;}

This is forward-time, centred-space (FTCS) for the heat equation. It needs only the previous time level — natural for first-order time.

Both schemes are explicit: the update at each grid point uses only known values from previous time levels. No linear system needs to be solved per step.

Stability: the CFL condition

The wave-equation leapfrog scheme is conditionally stable: it works only if the Courant number satisfies

\boxed{\;C \;=\; \frac{c\, \Delta t}{\Delta x} \;\leq\; 1.\;}

This is the Courant–Friedrichs–Lewy (CFL) condition. Cross the threshold and the simulation blows up exponentially within a few steps.

Courant number C = c·Δt/Δx = 0.900

The leapfrog finite-difference scheme for the 1-D wave equation is conditionally stable: it works only when the Courant number C = c·Δt/Δx ≤ 1. The threshold is sharp. At C = 0.9 the simulation is stable and shape is preserved; at C = 1.01 numerical errors grow exponentially and the simulation blows up within a few timesteps. The bound is fundamental — it comes from von Neumann stability analysis, which lifts the algebra into Fourier space and asks under what conditions each Fourier mode is preserved in amplitude rather than amplified.

Slide the Courant number across $C = 1$ . The boundary is sharp: at $C = 0.9$ the simulation runs fine; at $C = 1.01$ the numerical wave grows by a factor of $\sim 1.2$ per step and is unrecognisable by step 20. Try it.

▶ Deriving the CFL bound by von Neumann analysis

The standard way to analyse stability of a linear finite-difference scheme is to look at each Fourier mode separately and ask under what conditions the scheme amplifies or attenuates that mode. This is von Neumann stability analysis.

Step 1 — Try a Fourier mode. Substitute the ansatz $u^n_j = G^n e^{i k j \Delta x}$ into the leapfrog scheme. Here $G$ is the amplification factor — a complex number that tells you how much one timestep multiplies the mode’s amplitude by. Stability requires $|G| \leq 1$ for every spatial wavenumber $k$ .

Step 2 — Substitute. The scheme is

u^{n+1}_j \;=\; 2 u^n_j - u^{n-1}_j + C^2 (u^n_{j+1} - 2 u^n_j + u^n_{j-1}).

Plug in the ansatz and divide through by $G^{n-1} e^{i k j \Delta x}$ :

G^2 \;=\; 2 G - 1 + C^2 G (e^{i k \Delta x} - 2 + e^{-i k \Delta x}) \;=\; 2 G - 1 + 2 C^2 G (\cos(k \Delta x) - 1).

Use the identity $\cos\theta - 1 = -2 \sin^2(\theta/2)$ :

G^2 - 2 G + 1 + 4 C^2 G \sin^2(k \Delta x / 2) \;=\; 0.

Step 3 — Solve the quadratic in G. Let $\sigma \equiv 2 C \sin(k \Delta x / 2)$ . Then

G^2 - (2 - \sigma^2)\, G + 1 \;=\; 0

with roots

G \;=\; \frac{2 - \sigma^2 \pm \sqrt{(2 - \sigma^2)^2 - 4}}{2}.

Step 4 — Check $|G| \leq 1$ . The product of the two roots is $G_+ G_- = 1$ , so either both roots have $|G| = 1$ (and the scheme is on the stability boundary) or one has $|G| > 1$ (and the scheme is unstable).

The discriminant is $(2 - \sigma^2)^2 - 4$ . It is negative iff $|2 - \sigma^2| < 2$ , iff $0 < \sigma^2 < 4$ . When the discriminant is negative the roots are complex conjugates with magnitude exactly 1 — the marginal-stable case.

So the scheme is stable iff $\sigma^2 \leq 4$ , i.e. $|C \sin(k \Delta x / 2)| \leq 1$ , i.e. (taking the worst-case $\sin = 1$ ):

\boxed{\;C \;=\; \frac{c \Delta t}{\Delta x} \;\leq\; 1.\;}

The CFL condition is a statement that the numerical scheme must not propagate information faster than the underlying physics does. The numerical “signal speed” in the leapfrog scheme is $\Delta x / \Delta t$ — one cell per timestep. The physical signal speed in the wave equation is $c$ . CFL says $c \leq \Delta x / \Delta t$ , equivalently $C \leq 1$ .

Stability for the heat equation

The same analysis on FTCS for the heat equation produces a tighter condition:

\boxed{\;\frac{D\, \Delta t}{\Delta x^2} \;\leq\; \frac{1}{2}.\;}

Notice the asymmetry: it’s $\Delta x^2$ , not $\Delta x$ . Halving the grid spacing forces a quartering of the time step. This is brutal for fine spatial resolution; doubling the spatial resolution costs four times the run time. The standard fix is to use implicit schemes for parabolic equations, where the stability constraint disappears entirely — Crank–Nicolson is the workhorse. We won’t develop the implicit case here; the bookshelf’s HeatDiffusion interactive uses FTCS with a stability cap on the time step.

The reason heat has $\Delta x^2$ rather than $\Delta x$ : the heat equation has one time derivative against two spatial derivatives, so dimensionally the stability bound has to balance $\Delta t$ against $\Delta x^2$ . The wave equation has two time derivatives, balancing $\Delta t^2$ against $\Delta x^2$ , hence $\Delta t \propto \Delta x$ .

Two big lessons

Stability ≠ accuracy. A stable scheme can still be wildly inaccurate (large truncation error). And a high-order accurate scheme can be unconditionally unstable. The two analyses are orthogonal.

Implicit schemes trade work per step for freedom in step size. Explicit methods need just an evaluation per grid point per step; implicit methods solve a linear system, which is more expensive per step but lets you use any $\Delta t$ you like. For stiff parabolic problems where the stability bound on explicit schemes forces $\Delta t$ down to ridiculous values, implicit methods pay off enormously.

What we use this for

Every PDE simulation on this site uses one of these schemes:

WaveVsHeat (MATH-31) and HeatDiffusion (MATH-12) — leapfrog for the wave equation, FTCS for the heat equation, each with the appropriate CFL/parabolic stability cap baked in. The visible diffusion-then-shape-loss is the consequence of forward-Euler-in-time.
TravelingWave (MATH-11) — leapfrog with appropriate boundary handling.

The next lesson, 10.4, tackles elliptic PDEs — Laplace’s equation — where time is absent and the algorithm is a fundamentally different beast: iterative averaging rather than time-stepping.