Stability needs a test, not a feeling
The previous guide left us with a sharp division of labour: by the Lax equivalence theorem, for a consistent scheme on a well-posed linear problem, convergence holds if and only if the scheme is stable. Consistency is the easy half — you just Taylor-expand the stencil and confirm the truncation error vanishes as the grid shrinks. Stability is the genuinely hard half: it asks whether the tiny errors born at every grid point each step stay bounded as you take more and more steps, or whether they snowball. We need a concrete, mechanical test for it. That test is von Neumann analysis.
Here is the intuition before any formula. A linear update is just a machine that takes the numbers on the grid at time level n and produces the numbers at level n+1. Because it is linear, an error riding on top of the solution obeys the very same update rule the solution does. So the real question is: when this machine is applied over and over, does the error get amplified or damped? If we could break any error into simple pieces, watch what the machine does to each piece, and find that no piece ever grows, we would be done. Fourier modes are exactly those simple pieces.
Feed the scheme a single wave
The trick is to test the scheme on one pure Fourier mode and let linearity handle the rest. Write the grid value at point j and time level n as a complex wave U_j^n = G^n e^{i k x_j}, where k is the wavenumber (how many wiggles per unit length) and G is an unknown complex number called the amplification factor. The wave's spatial shape e^{i k x_j} is fixed; the only thing that changes from step to step is that it gets multiplied by G. After one step the amplitude is |G| times bigger; after n steps it is |G|^n times bigger. So the whole fate of this mode is decided by one number, |G|.
Now the demand becomes obvious. If |G| > 1 for even one wavenumber k, that mode multiplies itself bigger every step, and after enough steps it explodes — that is exactly numerical instability, a solution that blows up to garbage or NaN. If |G| <= 1 for ALL k, every mode stays bounded and the scheme is stable. So the von Neumann stability condition is simply: |G(k)| <= 1 for every wavenumber the grid can represent. (Strictly we allow |G| <= 1 + C*dt for a constant C, which still keeps growth bounded over a finite time; for the plain constant-coefficient schemes here, |G| <= 1 is the clean statement.)
- Substitute the trial mode U_j^n = G^n e^{i k x_j} into the difference scheme, replacing every U at every grid point and time level.
- Cancel the common factor G^n e^{i k x_j} that appears in every term; you are left with one algebraic equation for G in terms of k.
- Use the identity e^{i theta} + e^{-i theta} = 2 cos(theta) to collapse the neighbour terms into cosines; the bookkeeping becomes much cleaner.
- Solve for G(k), then find the WORST mode — the wavenumber that makes |G| largest — and demand even that worst |G| stays <= 1.
- The inequality that survives is your stability condition: usually a ceiling on the time step dt relative to the grid spacing dx.
A worked example: heat with FTCS
Take the heat equation u_t = alpha u_xx and the explicit FTCS scheme from the second guide: a forward difference in time and a centred second difference in space. The update is U_j^{n+1} = U_j^n + r (U_{j+1}^n - 2 U_j^n + U_{j-1}^n), where r = alpha dt / dx^2 is the single dimensionless number that controls everything. Substitute the trial mode and cancel the common factor, and the cosine identity turns the three neighbour terms into a tidy expression for G.
FTCS for u_t = alpha u_xx, with r = alpha dt / dx^2 :
U_j^{n+1} = U_j^n + r ( U_{j+1}^n - 2 U_j^n + U_{j-1}^n )
Insert U_j^n = G^n e^{i k x_j} , let theta = k dx , divide out G^n e^{i k x_j} :
G = 1 + r ( e^{i theta} - 2 + e^{-i theta} )
= 1 + r ( 2 cos(theta) - 2 )
= 1 - 4 r sin^2(theta/2) # since 1 - cos = 2 sin^2(half)
Worst mode theta = pi (the +,-,+ sawtooth) gives G = 1 - 4 r .
Need |G| <= 1 for ALL theta ==> -1 <= 1 - 4 r ==> r <= 1/2 .
Stability: alpha dt / dx^2 <= 1/2 .Read what r <= 1/2 actually means, because it is brutal. It says alpha dt <= dx^2 / 2, so the time step must shrink like dx^2. Halve the grid spacing to double your spatial resolution, and you must take FOUR times as many time steps to reach the same final time. Want 10x finer space? Pay 100x more time steps. This is the explicit scheme's curse on parabolic problems: refinement in space is punished quadratically in time. It is the concrete reason the implicit schemes BTCS and Crank-Nicolson exist — their von Neumann analysis gives |G| <= 1 for EVERY r, so they are unconditionally stable and you choose dt for accuracy alone, not survival.
The CFL condition: information cannot outrun the grid
For wave-like (hyperbolic) problems the same machinery produces a different and famous bound. Take the advection equation u_t + a u_x = 0, which simply transports a profile to the right at speed a. The natural explicit choice, the upwind scheme, gives an amplification factor whose modulus stays <= 1 exactly when the Courant number C = |a| dt / dx is at most 1. That single inequality, |a| dt / dx <= 1, is the CFL condition (after Courant, Friedrichs, and Lewy, 1928).
There is a beautiful physical picture behind the algebra. The true solution at a point depends only on its analytic domain of dependence: for advection at speed a, the answer at (x, t) is fixed entirely by the point that distance a*dt upstream. The scheme, meanwhile, has its own numerical domain of dependence — the set of grid points that actually feed into the update, reaching one cell of width dx per step. CFL is simply the demand that the numerical domain of dependence must CONTAIN the true one. If in one time step the physical signal travels farther than one grid cell (C > 1), the scheme is trying to compute an answer from data it never looked at — information has outrun the grid — and no clever weighting can rescue it.
When G is complex: dissipation and dispersion
Stability only asks about the SIZE of G, |G| <= 1. But G is a complex number, so it also has a PHASE, and the phase carries the rest of the story. The exact advection mode should keep its amplitude exactly (a wave neither grows nor shrinks as it travels) and shift its phase by precisely the amount that moves it at speed a. A real scheme misses on both counts, and the two kinds of miss have names. When |G| < 1 for the modes that should be preserved, the scheme is quietly draining their amplitude — that is numerical dissipation, and it shows up as a sharp front being smeared into a soft blur.
The phase error has its own name: numerical dispersion. If different wavenumbers travel at slightly different numerical speeds — when they should all move at the single speed a — a sharp pulse made of many modes drifts apart, and you see spurious wiggles, usually trailing a steep front like ripples behind a boat. Upwind schemes are heavily dissipative (fronts blur but stay smooth); the centred leapfrog scheme is nearly dissipation-free but strongly dispersive (fronts stay sharp but grow oscillatory wiggles). There is rarely a free lunch: damp the wiggles and you blur the front; preserve the front and you risk the wiggles.
Hold all of this honestly. Von Neumann analysis on its periodic, constant-coefficient model is the right first question to ask, and getting |G| <= 1 is non-negotiable — an unstable scheme produces literal garbage. But passing the stability test does NOT mean your answer is good: a stable scheme can dissipate or disperse a clean wave into mush, and on top of that everything runs in floating-point arithmetic, so the computed G is itself rounded and a marginally stable scheme sitting right at |G| = 1 can be nudged over the edge by round-off. The trustworthy verdict comes from convergence (the next rung's manufactured-solution tests), not from stability alone.