Why finite differences are not the last word
The four guides before this one built a complete, trustworthy machine: lay a grid, swap each derivative for a finite-difference stencil, check consistency with truncation error, check stability with von Neumann analysis and the CFL condition, and let the Lax theorem promise convergence. For a smooth solution of the heat or wave equation, that machine is hard to beat. But two cracks appear the moment the world gets less tame, and this guide is about the families of methods built to fill them.
The first crack is conservation. Many of the most important PDEs are conservation laws — they say a quantity (mass, momentum, energy) is neither created nor destroyed, only moved around, exactly the picture from deriving a PDE from a conservation law. A plain finite-difference scheme approximates the differential form u_t + f(u)_x = 0 pointwise, and there is no reason its little errors should respect the global budget. After a long run the total mass on the grid can quietly drift, which for a physical simulation is a quiet disaster.
The second crack is shocks. From the nonlinear rung you know that a smooth solution of a scalar conservation law like Burgers' equation u_t + u u_x = 0 can steepen until a shock forms in finite time, after which the honest object is a weak solution with a genuine jump. Centred finite differences placed across such a jump produce violent oscillations — the discrete cousin of the Gibbs phenomenon — because a chord across a cliff is a wild lie. We need methods that expect cliffs, not ones that assume smoothness everywhere.
Finite volumes: bookkeeping by the cell
The finite-volume method fixes the conservation crack with one decisive change of attitude. Instead of tracking the pointwise value U_j at a grid node, it tracks the cell average of u over a little interval (cell) C_j = [x_{j-1/2}, x_{j+1/2}]. Then it does not approximate the differential equation at all — it integrates the conservation law exactly over each cell. The integral of u_t over the cell is the rate of change of the cell's total amount; the integral of f(u)_x collapses, by the fundamental theorem of calculus, to the flux at the right face minus the flux at the left face. That is just honest bookkeeping: the stuff in a box changes only by what flows in one side and out the other.
cell average: Ubar_j ~ (1/h) * integral over C_j of u dx
update rule: Ubar_j^{n+1} = Ubar_j^n - (k/h) * ( F_{j+1/2} - F_{j-1/2} )
F_{j-1/2} cell C_j F_{j+1/2}
| Ubar_j (average) |
----+-------------------------------+----
x_{j-1/2} x_{j+1/2}
flux in --> --> flux outStare at that update rule and the magic of conservation jumps out. The numerical flux F_{j+1/2} is the flow across the face between cells j and j+1, and it appears with a plus sign for cell j and a minus sign for cell j+1. So whatever leaves one cell enters its neighbour exactly — when you sum over all cells, the interior fluxes cancel in pairs (a telescoping sum) and only the boundary fluxes remain. The total quantity on the grid is conserved to machine precision, by construction, no matter how coarse the mesh. This is the structural reason finite volumes dominate computational fluid dynamics.
Shock capturing: which way does the wind blow?
How should we choose the face flux? The most principled answer is the Godunov scheme. At each face, treat the two neighbouring cell averages as a tiny initial-value problem with a single jump — exactly a Riemann problem, the elementary building block from the nonlinear rung. Solve that little problem (it has an exact, known answer for simple fluxes), read off the value at the face, and evaluate the physical flux f there. Godunov's insight was that the local Riemann solution already knows which way information travels, so the scheme automatically takes its data from the upwind side.
That upwind instinct matters because it is what makes the scheme pick out the physically correct weak solution. Recall the honest warning from the nonlinear rung: once a shock forms, a weak solution is not unique, and only the one satisfying the entropy condition is real — the others are mathematical ghosts like spontaneous rarefaction shocks that never occur in nature. A well-built Godunov-type scheme has just enough built-in dissipation to reject those ghosts and converge to the entropy solution. It is the discrete echo of vanishing viscosity: a whisker of numerical smearing, applied in the right place, that selects the admissible jump.
But there is no free lunch, and here is the central tension of shock capturing. A first-order scheme like plain Godunov is rock-solid and oscillation-free, yet its built-in numerical diffusion smears a sharp shock over several cells into a gentle ramp — the cliff survives but loses its edge. A naive high-order scheme keeps the edge sharp but resurrects the oscillations near the jump. We want high accuracy in the smooth regions and a clean, monotone profile at the shock, and a famous theorem (Godunov's order barrier) says no fixed linear scheme can give both at once.
Flux limiters: the best of both worlds
The way out of Godunov's barrier is to be nonlinear — to let the scheme change its mind depending on what the solution looks like nearby. This is the job of a flux limiter. The idea is to blend two fluxes: a safe, low-order flux that never oscillates, and a sharp, high-order flux that resolves smooth regions beautifully. A limiter function inspects a local smoothness sensor — usually the ratio of consecutive slopes — and decides, cell by cell, how much of the high-order flux to trust.
- At each face, measure a smoothness ratio r from neighbouring slopes — large and positive where u is smooth, near zero or negative where a slope reverses (a sign of an extremum or a shock).
- Feed r into a limiter function (minmod, van Leer, superbee, ...) that returns a weight between 0 and a small cap.
- Where the solution is smooth, the weight is near 1: use almost the full high-order flux and enjoy second-order accuracy.
- Near a shock or an extremum, the weight collapses toward 0: fall back to the safe low-order flux, killing the oscillation.
The payoff is a scheme that is high-resolution: second-order accurate where the flow is smooth, yet free of spurious wiggles at the discontinuity. The technical goal a good limiter enforces is that the scheme be total-variation diminishing (TVD) — the total up-and-down movement of the discrete solution never grows, which is precisely the property that forbids new oscillations from being born. This is the genuine state of the art for everyday hyperbolic conservation laws, and it is why the cell-average viewpoint of finite volumes was worth setting up so carefully.
Spectral methods: spending smoothness on accuracy
Now swing to the opposite extreme. Suppose your solution is genuinely smooth — no shocks, no corners, perhaps a periodic flow or a wave with gentle features. Finite differences and finite volumes spend their points locally, each stencil seeing only a few neighbours, and that locality caps their accuracy. The spectral method makes the opposite bet: represent the unknown not by local grid values but by a sum of global smooth basis functions — a truncated Fourier series for periodic problems, or Chebyshev polynomials on an interval — and let every coefficient feel the whole domain at once.
The reward is breathtaking when it applies. A derivative becomes trivial: differentiating sin(n x) just multiplies by a factor in n, so applying u_xx is a single multiplication on each Fourier coefficient — and the Fast Fourier Transform moves between grid values and coefficients in nearly linear time. Better still, for a smooth periodic function the error of the truncated series decays faster than any power of the grid spacing — so-called spectral accuracy or exponential convergence. Where a second-order finite-difference scheme needs four times the points to cut the error by four, a spectral method can cut it by orders of magnitude. Doubling the resolution can buy you ten more correct digits.
Step back and the whole rung becomes a map of trade-offs rather than a single algorithm. Finite differences are simple and transparent; the finite-element method from the previous guide bends gracefully to complicated geometries; finite volumes guarantee conservation and tame shocks; spectral methods deliver extraordinary accuracy on smooth, simple domains. None is universally best. The mature numerical analyst's skill is not memorising one scheme but reading the structure of the problem — conservation? shocks? smoothness? geometry? — and choosing the method whose strengths line up with it.