Three questions hiding in one
In the last two guides you laid a grid, replaced derivatives by stencils, and built schemes — FTCS, BTCS, Crank-Nicolson — that march a heat profile forward in time. Each produces a table of numbers. But a nagging question remains: as you refine the grid, do those numbers actually approach the true solution of the partial differential equation, or are you just generating confident-looking nonsense? That single goal — the numerical answer approaching the truth as the grid shrinks — is called convergence, and it is the only thing you ultimately care about.
The trouble is that convergence is brutally hard to check head-on: to compare your answer with the true solution, you would need the true solution — which is exactly what you do not have. The whole reason you are computing is that the PDE cannot be solved by hand. So the genius of this subject is to split that one impossible-looking question into two checkable ones. The first is consistency: does the discrete equation, as h and k shrink, really turn back into the PDE you meant to solve? The second is stability: do small errors stay small as you take step after step, or do they snowball?
Consistency: are you even solving the right equation?
Consistency is the minimum sanity check, and it has a precise meaning. Take the EXACT solution of the PDE and substitute it into your difference scheme. The true solution will not satisfy the discrete equation exactly — there is a leftover, called the local truncation error, that measures by how much the true solution fails to obey the stencil. The scheme is consistent if that leftover tends to zero as the grid spacings h and k tend to zero. A scheme that is NOT consistent is solving some other equation entirely, and refining the grid will not save it.
You compute the truncation error with Taylor series — the same tool that built the stencils in the first place. Expand every grid value about a central point, and the PDE's own terms cancel out, leaving the error as a series in powers of h and k. The lowest surviving powers give the scheme's order. For FTCS on the heat equation u_t = alpha*u_xx, the leftover works out to a term of size O(k) plus a term of size O(h^2): first order in time, second order in space, written O(k + h^2). Crank-Nicolson, by centring at the half time-step, lifts the time part to O(k^2 + h^2). Both vanish as the grid shrinks, so both are consistent.
Plug the TRUE solution u into FTCS for u_t = alpha*u_xx :
(u_j^{n+1} - u_j^n)/k - alpha*(u_{j-1}^n - 2 u_j^n + u_{j+1}^n)/h^2
= [u_t - alpha*u_xx] <- the PDE: equals 0
+ (k/2) u_tt <- O(k) time error
- (alpha h^2 / 12) u_xxxx <- O(h^2) space error
+ (higher-order terms)
local truncation error = O(k) + O(h^2) -> 0 as k,h -> 0 (CONSISTENT)Stability: do the mistakes snowball?
Consistency alone is almost worthless, and here is the vivid reason. A consistent scheme has a tiny truncation error at each step — but a numerical solution is thousands of steps long, and at every cell there is also a tiny round-off blip from floating-point arithmetic. Stability asks: do those blips stay bounded as you march, or does each step amplify them so they swell into a tidal wave that drowns the real answer? You saw this concretely with FTCS: push the ratio r = alpha*k/h^2 past 1/2 and within a handful of steps the numbers detonate into wild oscillations.
The everyday tool for checking stability is von Neumann analysis — the subject of the very next guide — so here just hold the idea: write a single error wave, see what factor each step multiplies it by (the amplification factor g), and demand |g| <= 1 for every wavenumber so no mode grows without bound. The crucial honest point is that stability is the genuinely SUBTLE leg. Consistency you check almost mechanically with Taylor series; stability is where schemes that look perfectly reasonable quietly fail. FTCS is consistent yet only conditionally stable; BTCS and Crank-Nicolson are unconditionally stable. That difference between explicit and implicit schemes is, at bottom, a difference in stability.
The Lax equivalence theorem: the bridge
Now the payoff. The Lax equivalence theorem (also called Lax-Richtmyer) states it with startling cleanliness: for a well-posed linear initial-value problem and a consistent finite-difference scheme, the scheme CONVERGES if and only if it is STABLE. In slogan form, the one every numerical analyst keeps in their pocket: consistency + stability = convergence. The hard, unobservable question (does my answer reach the truth?) is replaced by two checkable ones, and once consistency is established, stability becomes EQUIVALENT to convergence — you get the prize you wanted by proving the part you can actually prove.
Watch why neither leg can stand alone. A consistent-but-unstable scheme has the right target — it tends to the PDE — but its errors blow up, so the numbers fly off to infinity and never converge: think FTCS run with r = 0.6. A stable-but-inconsistent scheme keeps its errors bounded, but it is faithfully solving the WRONG equation, so the numbers settle calmly onto the wrong answer. Only when both hold do the numbers settle onto the RIGHT answer. The theorem is the reason every textbook spends pages on stability and almost none on convergence proofs: for a consistent scheme, proving stability IS proving convergence.
- Confirm the problem is well-posed and linear (the heat and wave equations qualify; this is where the theorem lives).
- Check consistency: Taylor-expand the scheme, watch the PDE cancel, confirm the truncation error -> 0 as h, k -> 0, and read off the order.
- Check stability: run von Neumann analysis to find the amplification factor g and the condition |g| <= 1 for all wavenumbers.
- Invoke Lax: consistent + stable => convergent — no separate, unobservable convergence proof needed.
What it tells you in practice — and where it stops
Run the loop on FTCS for the heat equation and the whole story snaps together. Consistency: order O(k + h^2), so the truncation error vanishes. Stability: von Neumann gives g = 1 - 4r*sin^2(theta/2), which satisfies |g| <= 1 exactly when r = alpha*k/h^2 <= 1/2. By Lax, FTCS therefore CONVERGES if and only if you keep r <= 1/2 — and when it does, the global error inherits the consistency order, falling like O(k) + O(h^2). That is the practical meaning of the order of convergence: halve h (and quarter k to stay under the stability bound) and the error drops by about four. The implicit schemes dodge the r <= 1/2 leash entirely, which is exactly why you reach for them on stiff problems.
Be honest about the fine print, because the slogan is often over-applied. The Lax theorem is stated for LINEAR, well-posed problems, and assumes consistency holds in a suitable norm. For nonlinear PDEs — shocks in compressible flow, turbulence — there is no such clean equivalence: consistency plus linear stability does NOT guarantee convergence to the physically correct (entropy-satisfying) solution, and you need stronger machinery (the Lax-Wendroff theorem, total-variation-diminishing bounds, entropy conditions). And remember the round-off floor from earlier rungs: every computed answer is approximate, floating-point is not exact real arithmetic, so shrinking h forever does not keep helping — past a point, round-off creeps back in. Convergence is an as-h-tends-to-zero promise, not a guarantee at the one finite grid you actually run.