Finite differences for hyperbolic equations: upwind, Lax-Wendroff and the CFL condition
The advection equation, instability of the centred scheme by von Neumann analysis, upwind and numerical diffusion, Lax-Friedrichs, Lax-Wendroff, Godunov's theorem and transmission lines.
This in-depth page completes the lesson on finite differences with the class of equations that is hardest to solve numerically: hyperbolic equations, which describe transport and propagation at finite speed, such as a pollutant carried by a river, a signal on a transmission line, a shock wave. We start from the simplest model, the linear advection equation, and see why the most “natural” scheme is unstable, how it is fixed and what the Courant condition says.
The advection equation
![\[ u_t + a\,u_x = 0,\qquad a \gt 0 \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-0b1cafa3d21d825c9d6084ae4248f071_l3.png "Rendered by QuickLaTeX.com")
has the exact solution 
: the initial profile moves right at speed 
without changing shape. Along the characteristics 
the solution is constant: information travels from left to right. It is an ideal test case, because the exact solution is known and any distortion of the profile is an error of the method.
The centred scheme is unstable
With a forward difference in time and a centred one in space, the most accurate:
![\[ u_i^{n+1} = u_i^n - \frac C2\big(u_{i+1}^n - u_{i-1}^n\big),\qquad C = \frac{a\,\Delta t}{\Delta x} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-6c51b0359e7562d3198ff6d7664d9f0b_l3.png "Rendered by QuickLaTeX.com")
It looks reasonable, but it is unstable for every 
. To see this, use von Neumann analysis: substitute 
(a wave of wavenumber 
) and compute the amplification factor 
:
![\[ g = 1 - iC\sin\theta,\qquad |g|^2 = 1 + C^2\sin^2\theta \gt 1 \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-24dacf7b66686d2ebd78e05d410b4980_l3.png "Rendered by QuickLaTeX.com")
Every component grows at every step. The scheme is consistent but not stable, and by the Lax-Richtmyer theorem it does not converge.
The upwind scheme
The fix is to look where the information comes from. With 
it comes from the left, so use the backward difference:
![\[ u_i^{n+1} = u_i^n - C\big(u_i^n - u_{i-1}^n\big) = (1 - C)\,u_i^n + C\,u_{i-1}^n \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-58d8e97b50e3025472814b905612b191_l3.png "Rendered by QuickLaTeX.com")
If 
the new value is a weighted average with positive weights: the scheme is stable and creates no oscillations. The amplification factor has 
exactly when 
: the Courant-Friedrichs-Lewy condition. With 
the scheme is exact (
, the profile moves one cell per step).
The price is accuracy: upwind is first order and introduces numerical diffusion. A step is carried along but its corners are gradually rounded, as if there were a small diffusion term the true equation does not have. Writing the modified equation shows it: the scheme actually solves 
.
Lax-Friedrichs and Lax-Wendroff
- Lax-Friedrichs replaces
in the centred scheme with the average of the neighbours, 
. It becomes stable for , but it is even more diffusive than upwind. - Lax-Wendroff uses the second-order Taylor expansion in time and the equation itself to substitute
:
![\[ u_i^{n+1} = u_i^n - \frac C2\big(u_{i+1}^n - u_{i-1}^n\big) + \frac{C^2}{2}\big(u_{i+1}^n - 2u_i^n + u_{i-1}^n\big) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-db31c1e18c1e453aed867334eb979ac5_l3.png "Rendered by QuickLaTeX.com")
It is second order and stable for , with little diffusion. But near jumps it produces spurious oscillations (overshoots before and after the step), similar to the Gibbs phenomenon of Fourier series.
Godunov’s theorem and modern schemes
In 1959 Godunov proved that no linear scheme of order higher than one can carry a jump without creating oscillations. Either you accept diffusion (first order) or oscillations (higher order). Modern schemes get around the theorem by being nonlinear: they use a second-order method where the solution is smooth and switch automatically to a first-order one near jumps, with so-called flux limiters (TVD, MUSCL, WENO schemes). These are the schemes used in computational fluid dynamics programs to simulate shock waves and supersonic flows.
Hyperbolic systems and variable coefficients
The telegrapher’s equations for a transmission line,
![\[ \frac{\partial v}{\partial t} = -\frac1C\frac{\partial i}{\partial x},\qquad \frac{\partial i}{\partial t} = -\frac1L\frac{\partial v}{\partial x} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-6a83aa1233d9e1c834788dc16acda589_l3.png "Rendered by QuickLaTeX.com")
(with 
and 
the inductance and capacitance per unit length), form a hyperbolic system. Diagonalising the coefficient matrix (the eigenvalues are 
) splits the system into two advection equations, one wave going forward and one coming back: the same forward and reflected waves of transmission line analysis. Upwind is applied to each in the right direction, and the Courant condition is imposed with the largest speed. With variable coefficients (a line whose impedance changes along its length) the local speed changes from point to point, and the time step must be chosen from the maximum speed over the whole grid.
