Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
Hyperbolic Equations: Advection, Waves and the CFL Condition
60m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

Hyperbolic Equations: Advection, Waves and the CFL Condition

Advection and wave equations, upwind, Lax–Friedrichs and Lax–Wendroff, the CFL condition, numerical diffusion and dispersion

1922: a forecast six weeks in the making
Lewis Fry Richardson — whom we met in Unit 4 relaxing stresses in a masonry dam — believed that weather could be forecast by computing the equations of fluid motion forward from today's observations. In 1922 he published Weather Prediction by Numerical Process and, as proof, a single forecast: starting from observations over central Europe at 7 a.m. on 20 May 1910, the conditions six hours later, computed by hand in the intervals of his service with an ambulance unit on the Western Front. His forecast predicted that surface pressure would rise by 145 hPa in six hours — a change larger than any hurricane produces. In reality the pressure was almost steady. Richardson published the failure honestly. Decades later the meteorologist Peter Lynch showed that the method was sound and the trouble lay in unbalanced initial data, which excited fast gravity waves that the computation could not control. In 1928 Richard Courant, Kurt Friedrichs and Hans Lewy in Göttingen proved the theorem that tells every wave computation how fast it is allowed to march. In 1950, respecting it, Jule Charney, Ragnar Fjørtoft and John von Neumann ran the first successful numerical forecast on ENIAC.
Hyperbolic equations describe things that travel at finite speed. The simplest is advection, ut+a ux=0u_t + a\,u_x = 0: any initial profile moves at speed aa without changing shape, u(x,t)=u0(x−at)u(x, t) = u_0(x - at). A pollutant carried by a river, a temperature front in a pipe, a density wave in traffic. The wave equation utt=c2uxxu_{tt} = c^2u_{xx} — a vibrating string, sound, a pressure surge in a penstock — consists of two advections in opposite directions. Hyperbolic solutions keep their sharp features: a shock or a step does not smooth out. That is exactly what makes them hard to compute.
💡The CFL condition
In time Δt\Delta t the true solution at xix_i depends on data a distance aΔta\Delta t upstream. An explicit scheme computes uin+1u_i^{n+1} from neighbours within a few Δx\Delta x. If the physical dependence reaches further than the numerical stencil, the scheme cannot possibly be right, and in fact it is unstable. For a one-point-each-side stencil this gives the Courant–Friedrichs–Lewy condition: C=∣a∣ ΔtΔx≤1.C = \dfrac{|a|\,\Delta t}{\Delta x} \le 1. CC is the Courant number. Note the difference from diffusion: Δt∝Δx\Delta t \propto \Delta x, not Δx2\Delta x^2. The condition is necessary for every explicit scheme; it is not sufficient — a scheme can satisfy CFL and still be unstable (centred differences in space with forward Euler in time are unstable for any C>0C > 0).
Derivation
The exact solution, and why upwind is stable but centred is not
Translation solves it. For u=u0(x−at)u = u_0(x - at) the chain rule gives ut=−a u0′(x−at)u_t = -a\,u_0'(x - at) and ux=u0′(x−at)u_x = u_0'(x - at), so ut+aux=0u_t + au_x = 0 for any profile u0u_0.
Upwind. Insert the Fourier mode ujn=Gneiφju_j^n = G^ne^{\mathrm i\varphi j} (lesson 10-2) into ujn+1=ujn−C(ujn−uj−1n)u_j^{n+1} = u_j^n - C(u_j^n - u_{j-1}^n):
G=1−C(1−e−iφ)=(1−C+Ccos⁡φ)−i Csin⁡φ.G = 1 - C\big(1 - e^{-\mathrm i\varphi}\big) = (1 - C + C\cos\varphi) - \mathrm i\,C\sin\varphi.
Its squared size, after expanding and using cos⁡2+sin⁡2=1\cos^2 + \sin^2 = 1:
∣G∣2=(1−C+Ccos⁡φ)2+C2sin⁡2φ=1−2C(1−C)(1−cos⁡φ).|G|^2 = (1 - C + C\cos\varphi)^2 + C^2\sin^2\varphi = 1 - 2C(1 - C)(1 - \cos\varphi).
Since 1−cos⁡φ≥01 - \cos\varphi \ge 0, ∣G∣≤1|G| \le 1 for every mode exactly when C(1−C)≥0C(1 - C) \ge 0: 0≤C≤10 \le C \le 1, the CFL condition.
Centred differences with forward Euler, ujn+1=ujn−C2(uj+1n−uj−1n)u_j^{n+1} = u_j^n - \tfrac C2(u_{j+1}^n - u_{j-1}^n), give
G=1−C2(eiφ−e−iφ)=1−iCsin⁡φ,∣G∣2=1+C2sin⁡2φ>1:G = 1 - \tfrac C2\big(e^{\mathrm i\varphi} - e^{-\mathrm i\varphi}\big) = 1 - \mathrm iC\sin\varphi, \qquad |G|^2 = 1 + C^2\sin^2\varphi > 1:
every mode grows, at any Courant number. Satisfying CFL is not enough.
upwind: uin+1=uin−C (uin−ui−1n),Lax–Wendroff: uin+1=uin−C2(ui+1n−ui−1n)+C22(ui+1n−2uin+ui−1n)\text{upwind: } u_i^{n+1} = u_i^n - C\,(u_i^n - u_{i-1}^n),\qquad \text{Lax–Wendroff: } u_i^{n+1} = u_i^n - \tfrac C2(u_{i+1}^n - u_{i-1}^n) + \tfrac{C^2}{2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)

Two classic explicit schemes for u_t + a u_x = 0 (a > 0). Upwind takes its difference from the side the information comes from; first order. Lax–Wendroff adds a second-derivative correction from the Taylor series in time; second order. Both are stable for C ≤ 1.

Derivation
Lax–Wendroff from a Taylor series in time, and what upwind really solves
Lax–Wendroff. Step forward in time with two Taylor terms, u(t+Δt)=u+Δt ut+12Δt2uttu(t + \Delta t) = u + \Delta t\,u_t + \tfrac12\Delta t^2u_{tt}, and use the equation itself to turn time derivatives into space derivatives: ut=−auxu_t = -au_x, and differentiating again, utt=−auxt=a2uxxu_{tt} = -au_{xt} = a^2u_{xx}. So
u(t+Δt)≈u−aΔt ux+12a2Δt2 uxx.u(t + \Delta t) \approx u - a\Delta t\,u_x + \tfrac12a^2\Delta t^2\,u_{xx}.
Replace uxu_x by the central difference (uj+1−uj−1)/2Δx(u_{j+1} - u_{j-1})/2\Delta x and uxxu_{xx} by the second difference (uj+1−2uj+uj−1)/Δx2(u_{j+1} - 2u_j + u_{j-1})/\Delta x^2. With C=aΔt/ΔxC = a\Delta t/\Delta x the coefficients become C/2C/2 and C2/2C^2/2 — the formula above.
The modified equation of upwind. Ask which equation the upwind scheme solves exactly to the next order. Expand both differences by Taylor: (ujn+1−ujn)/Δt=ut+Δt2utt+…(u_j^{n+1} - u_j^n)/\Delta t = u_t + \tfrac{\Delta t}{2}u_{tt} + \dots and (uj−uj−1)/Δx=ux−Δx2uxx+…(u_j - u_{j-1})/\Delta x = u_x - \tfrac{\Delta x}{2}u_{xx} + \dots. The scheme says their weighted sum is zero:
ut+aux=aΔx2uxx−Δt2utt=aΔx2(1−C)uxx,u_t + au_x = \frac{a\Delta x}{2}u_{xx} - \frac{\Delta t}{2}u_{tt} = \frac{a\Delta x}{2}\big(1 - C\big)u_{xx},
using utt≈a2uxxu_{tt} \approx a^2u_{xx} once more. The right side is a diffusion term with coefficient aΔx2(1−C)\tfrac{a\Delta x}{2}(1 - C): upwind computes advection plus artificial diffusion, which is the smearing. It vanishes at C=1C = 1, where upwind shifts the data exactly one cell per step.
💡Numerical diffusion and dispersion
Modified-equation analysis shows what each scheme actually solves. Upwind solves advection plus a diffusion term aΔx2(1−C)uxx\tfrac{a\Delta x}{2}(1 - C)u_{xx}: it smears. Carrying a square pulse once around a periodic domain at C=0.5C = 0.5, each sharp edge is spread over 0.18 of the domain; at C=0.8C = 0.8 over 0.10; at exactly C=1C = 1 upwind is exact (it shifts by one cell per step). Lax–Wendroff's leading error is a third derivative — dispersion: different wavelengths travel at slightly different speeds, so a smooth Gaussian is carried with error 0.1 but a square pulse grows ripples that overshoot to 1.23 and undershoot to −0.23. A negative concentration of pollutant is not a small error; it is nonsense. Modern high-resolution schemes (flux limiters, Godunov, WENO) switch between the two behaviours to keep fronts sharp without oscillations.
Young man, in mathematics you don't understand things. You just get used to them.
– John von Neumann, reply to a physicist who did not understand the method of characteristics, as quoted in The Dancing Wu Li Masters (1979)
📐A plucked string
A guitar string 0.65 m long with wave speed 120 m/s (fundamental 92.3 Hz), plucked 3 mm at a quarter of its length, stepped with the centred leapfrog scheme. At C=1C = 1 the scheme is exact for the wave equation: after one period the shape returns to within 10−1410^{-14} mm. At C=0.9C = 0.9 it returns with a 0.03 mm error from dispersion. At C=1.01C = 1.01 — one percent over the limit — it blows up to 102510^{25} mm in a single period.
Derivation
The leapfrog scheme for the wave equation, and why C=1C = 1 is exact
The scheme. Replace both second derivatives in utt=c2uxxu_{tt} = c^2u_{xx} by second differences (lesson 3-3), one in time and one in space:
ujn+1−2ujn+ujn−1Δt2=c2 uj+1n−2ujn+uj−1nΔx2.\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{\Delta t^2} = c^2\,\frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{\Delta x^2}.
Solve for the new value, with C=cΔt/ΔxC = c\Delta t/\Delta x:
ujn+1=2ujn−ujn−1+C2(uj+1n−2ujn+uj−1n).u_j^{n+1} = 2u_j^n - u_j^{n-1} + C^2\big(u_{j+1}^n - 2u_j^n + u_{j-1}^n\big).
It needs two time levels, like a two-step method (lesson 9-3). To start from rest, the fictitious u−1u^{-1} is set equal to u1u^{1}, which gives a first step with half the second difference — the lab's hint.
Why C=1C = 1 is exact. Then the ujnu_j^n terms cancel and the scheme is ujn+1=uj+1n+uj−1n−ujn−1u_j^{n+1} = u_{j+1}^n + u_{j-1}^n - u_j^{n-1}. Every solution of the wave equation is a sum of two travelling waves, u=F(x−ct)+G(x+ct)u = F(x - ct) + G(x + ct) (d'Alembert). With Δx=cΔt\Delta x = c\Delta t, a wave moving right moves exactly one cell per step, so FF at (j,n+1)(j, n+1) equals FF at (j−1,n)(j - 1, n), and also FF at (j,n−1)(j, n-1) equals FF at (j+1,n)(j + 1, n); together these make the formula hold exactly for FF, and likewise for GG. Any other CC makes waves of different lengths travel at slightly different speeds — the dispersion seen at C=0.9C = 0.9 — and C>1C > 1 violates the CFL condition.
💡Where you will use this
CFD and Aerodynamics: the Euler and Navier–Stokes equations are hyperbolic in their convective part; every explicit solver chooses its time step from a CFL number. Rocket Propulsion: unsteady nozzle flow, combustion instability, water hammer in feed lines. Hydraulics: river flood waves, dam-break flows and pressure surges in pipelines (method of characteristics). Structures: impact and blast analysis with explicit dynamics, where the time step is set by the fastest wave in the smallest element.
📚Further reading
  • Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 10, 'Advection Equations and Hyperbolic Systems' — upwind, Lax–Wendroff, modified equations and the CFL condition.
  • Randall LeVeque, Finite Volume Methods for Hyperbolic Problems (Cambridge University Press, 2002) — high-resolution methods that keep fronts sharp without ripples.
  • Peter Lynch, The Emergence of Numerical Weather Prediction: Richardson's Dream (Cambridge University Press, 2006) — Richardson's forecast recomputed, and why it failed.
↓But there is a problem
Diffusion and waves both evolve in time. A great deal of engineering asks about the state that remains after everything has settled: the steady temperature in a heat sink, the stress in a loaded bracket, the pressure in groundwater flowing steadily under a dam. Time disappears, and every point balances its neighbours at once. On 5 June 1976 the Teton Dam in Idaho, newly completed, was filling for the first time. In the morning muddy water appeared at the downstream face. By noon a hole had eaten through the embankment, and the dam collapsed, sending a flood down the valley that killed eleven people. The water had found a path through the soil — a path governed by Laplace's equation.