Numerical Stability

The explicit scheme is wonderfully simple — but it hides a trap that has nothing to do with ordinary inaccuracy. Every computer stores numbers with only finitely many digits, so a tiny rounding error is baked into every value on the grid, at every single step. Normally that is harmless: the error just sits there, far too small to matter. But if the scheme is unstable, each time step multiplies that tiny error by more than one in magnitude — and an error that doubles (or worse) every step becomes enormous astonishingly fast. After only a few dozen steps a perfectly sensible-looking simulation can be pure numerical garbage, even though every line of the code is correct.

This is a genuinely different failure mode from simply using too coarse a grid. A method is stable only if errors stay bounded as it marches forward; an unstable method guarantees they explode, no matter how carefully everything else was set up.

The stability condition

Von Neumann's test asks how the scheme treats a single Fourier mode u_i^n = G^n e^{i k x_i}. Substituting into the update gives an amplification factor

G = 1 - 4r\sin^2\!\Big(\frac{k\,\Delta x}{2}\Big).

Stability needs |G| \le 1 for every mode. The danger is the highest frequency, where \sin^2 = 1 and G = 1 - 4r. Requiring |1 - 4r| \le 1 gives the famous bound — the one-dimensional heat-equation version of the Courant–Friedrichs–Lewy (CFL) condition:

\boxed{\,r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \tfrac12.\,}

Cross r = \tfrac12 and the highest mode has |G| > 1: it grows by a factor each step, flipping sign as it goes, and swamps the true solution. This is a conditional stability — and an expensive one, since halving \Delta x forces \Delta t down by a factor of four. Violate it, and it does not matter how carefully everything else about the simulation was set up: the numbers will explode.

The explicit heat scheme amplifies mode k by G = 1 - 4r\sin^2(k\Delta x/2).
Stable iff |G| \le 1 for all modes, i.e. r = \alpha\Delta t/\Delta x^2 \le \tfrac12.
Beyond r = \tfrac12 the highest frequency blows up — a sign-flipping, doubling instability.

Worked example: watch the amplification factor

Take the worst-case, highest-frequency mode, where G = 1 - 4r exactly. Compare two mesh ratios:

r = 0.4: \quad G = 1 - 4(0.4) = -0.6, \qquad |G| = 0.6 < 1 \;\; \text{(stable)}.

r = 0.6: \quad G = 1 - 4(0.6) = -1.4, \qquad |G| = 1.4 > 1 \;\; \text{(unstable)}.

An error of size \varepsilon in that mode grows to G^n \varepsilon after n steps. With r = 0.4 it shrinks steadily: after 5 steps it is only (0.6)^5 \varepsilon \approx 0.078\,\varepsilon — under a tenth of its starting size. With r = 0.6, though, it is (-1.4)^5 \varepsilon \approx 5.38\,\varepsilon — more than five times larger, and it has flipped sign at every single step along the way. Keep going another 20 steps and that factor is well past a million. That is the entire mechanism of numerical blow-up, in two lines of arithmetic.

Push it past one-half

The explicit scheme run on the initial profile \sin x (which should simply decay). The bold curve is the numerical solution after a fixed time; the faint curve is the exact answer e^{-t}\sin x. Keep the mesh ratio r \le \tfrac12 and they agree; nudge r past 0.5 and the numerical curve erupts into a saw-toothed blow-up — a rapidly growing, sign-flipping oscillation from grid point to grid point — while the true solution sits quietly near zero. The visual signature of instability is unmistakable once you know to look for it: a jagged zig-zag whose height keeps doubling, not a smooth curve that merely drifts off course.

The practical fix

Once a scheme is diagnosed as unstable, there are exactly two honest ways out:

Shrink the time step. Reduce \Delta t until r = \alpha\Delta t/\Delta x^2 \le \tfrac12 holds again. This always works for the explicit scheme, but it can be painfully slow: a fine spatial grid (small \Delta x) forces an extremely small \Delta t, since the bound depends on \Delta x^2.
Switch to an implicit scheme. Methods such as Crank–Nicolson write the spatial difference at the new time level instead of the old one, which requires solving a small linear system at every step — more work per step — but in exchange the scheme is unconditionally stable: it never blows up, whatever \Delta t you choose. You trade cheap-but-restricted steps for expensive-but-unrestricted ones.

Which is better depends on the problem: if the physics itself changes slowly, a big stable implicit step can win outright; if it changes fast, you may need a small step anyway, and the cheap explicit scheme is hard to beat.

Concretely: shrinking \Delta x by a factor of ten, to resolve some fine detail, shrinks the explicit stability limit on \Delta t by a factor of one hundred (since r \propto \Delta t/\Delta x^2) — so reaching the same final time now costs a hundred times as many steps. An implicit scheme sidesteps that entirely: it might comfortably take a step fifty times larger than the explicit limit allows, needing only a tiny fraction of the steps, even though each of those steps individually costs more (solving a linear system rather than one line of arithmetic). For a very fine grid or a long time span, that trade tips decisively in the implicit scheme's favour.

Instability is not the same thing as inaccuracy, and it does not look the same either. An inaccurate-but-stable scheme gives you a wrong answer that stays wrong by roughly the same amount, step after step — you could, in principle, spot the error by comparing with a known case. An unstable scheme is far sneakier: for the first handful of steps it can look completely reasonable, tracking the true solution closely, before the smallest sliver of rounding error finally gets amplified enough to dominate — and then it diverges catastrophically, often within just a few more steps.

That is exactly why you check the stability condition before running a long simulation, using pure algebra (as above), rather than trusting a quick glance at the first few outputs. "It looked fine for the first ten steps" is not evidence of stability — it may just mean you haven't run it for the eleventh step yet.

The CFL condition is named after three mathematicians — Richard Courant, Kurt Friedrichs, and Hans Lewy — who derived it in 1928, to study whether certain finite-difference schemes could even in principle converge to the true solution. At the time there was no electronic computer on Earth capable of actually running one of these schemes; the entire analysis was pencil-and-paper pure mathematics, decades ahead of any machine that would need it.

Then digital computers arrived in the 1940s and 50s, numerical simulation of PDEs became routine, and the CFL condition turned out to be one of the single most load-bearing results in the whole field — every weather model, every fluid simulation, every finite-difference method you will ever meet depends on getting it right. A purely theoretical curiosity became, twenty years later, an indispensable engineering rule.