Stability, Consistency, and Convergence of Numerical Schemes

Every numerical simulation of a partial differential equation makes the same silent bargain: it throws away the continuum and keeps a finite grid, hoping the answer it computes is close to the answer it cannot. Whether that bargain is honest rests on three properties that sound alike but are logically distinct. Consistency asks whether the discrete equation, fed the exact solution, leaves only a small residual — the local truncation error — that vanishes as the mesh shrinks. Stability asks whether the scheme keeps small errors small as it marches forward, rather than amplifying round-off into nonsense. Convergence asks the only question a physicist ultimately cares about: does the computed field approach the true one as the grid is refined? The trouble is that convergence is exactly the thing you cannot check directly — if you knew the exact solution you would not be simulating. The triumph of mid-century numerical analysis, the Lax equivalence theorem of 1956, makes the untestable question testable: for a linear well-posed initial-value problem with a consistent scheme, convergence holds if and only if the scheme is stable. Consistency you verify with a Taylor expansion; stability with a Fourier analysis; the two together certify the one you want. The everyday tool for the stability half is von Neumann analysis — feed the scheme a Fourier mode, read off its per-step amplification factor, demand it not blow up — whose most famous product is the Courant–Friedrichs–Lewy condition: a scheme must not be asked to carry information faster than its stencil allows. This article assembles those pieces into the organizing principle of computational physics: consistency plus stability equals convergence.

Key takeaways

Convergence — the numerical solution approaching the true one as the mesh is refined — is what you want but cannot test directly, because the exact solution is unknown.
Consistency is checkable by Taylor expansion: the local truncation error must satisfy τ → 0 as Δt, Δx → 0, so the discrete operator really does reproduce the PDE in the limit.
Stability is checkable by Fourier or operator-norm analysis: the powers of the update operator stay uniformly bounded, ‖C(Δt)ⁿ‖ ≤ K for nΔt ≤ T, so errors do not explode.
The Lax equivalence theorem ties them together: for a consistent scheme on a well-posed linear problem, stability is necessary and sufficient for convergence — consistency + stability ⇔ convergence.
Von Neumann analysis substitutes a mode uⱼⁿ = gⁿ e^(ikjΔx) and reads off the amplification factor; the condition |g(k)| ≤ 1 + O(Δt) produces bounds like the diffusion limit r ≤ 1/2 and the Courant limit ν ≤ 1. Order of accuracy, by contrast, comes from consistency alone and says nothing about stability.
Key idea: convergence is what you want but cannot test; consistency and stability are what you can test — and the Lax equivalence theorem says, for consistent schemes on well-posed linear problems, those two checkable properties together are exactly equivalent to the one you want.

The core insight

Convergence is what you want but cannot test; consistency and stability are what you can test. The Lax equivalence theorem says that for a consistent scheme approximating a well-posed linear initial-value problem, stability is necessary and sufficient for convergence — so the two properties you can check certify the one you cannot.

Where this lives today

The trio consistency–stability–convergence is the daily grammar of every large simulation code. In computational fluid dynamics, climate models, and magnetohydrodynamics the time step is set automatically each cycle by a CFL number below one, recomputed from the fastest local wave speed — the direct descendant of the 1928 condition, now driving exascale runs on millions of cores. Numerical relativity is the sharpest modern test: to compute the gravitational-wave templates that LIGO matched against GW150914, evolution codes had to march Einstein's equations stably through a binary black-hole merger, and it took the 2005 breakthroughs in well-posed formulations and gauge conditions (moving punctures, the BSSN and generalized-harmonic systems) to make those evolutions stable at all — a reminder that well-posedness is a hypothesis of the Lax theorem, not a given. The live frontier is where the linear theory runs out: for nonlinear conservation laws with shocks the equivalence theorem gives no guarantee, and the design of provably stable schemes — total-variation-diminishing, entropy-stable, and structure-preserving methods — carries Lax's original program forward.

The physical picture

A time-marching scheme is one operator applied over and over: uⁿ⁺¹ = C uⁿ, so uᴺ = Cᴺ u⁰. Consistency is a statement about a single step — one application of C nearly reproduces the true evolution over Δt. Convergence is a statement about many steps — that single nearly-right step raised to the Nth power. The gap is the whole story: a per-step error that looks negligible can be multiplied by an exploding ‖Cᴺ‖ until it swamps the answer, or held in check by a bounded ‖Cᴺ‖. Stability is precisely that bound — which is why it, not the order of accuracy, decides whether a simulation lives or dies. The Courant–Friedrichs–Lewy condition puts a physical face on that bound. In one time step a signal in a hyperbolic PDE travels a distance c Δt, while an explicit scheme, reading only neighbouring cells, passes information a distance Δx per step. If c Δt exceeds Δx, the true domain of dependence reaches outside the points the scheme looked at — it computes from data it never saw, and goes unstable. The Courant number ν = c Δt / Δx measures that ratio, shown in the first figure. Von Neumann analysis supplies the algebra: on a periodic constant-coefficient grid the modes e^(ikx) are eigenvectors of C, each multiplied per step by g(k), so boundedness for all time collapses to sup_k |g(k)| ≤ 1 + C Δt — the second figure.

Figure

Space-time diagram of the CFL condition: from the apex grid point the numerical scheme depends only on the shaded stencil cone, whose edges have slope Δx per step. With Courant number ν = cΔt/Δx = 0.8 the physical characteristic stays inside the cone and the scheme is stable; with ν = 1.4 the characteristic escapes the cone and no explicit scheme can converge

The second figure is the algebra the physicist actually runs. For the explicit forward-time centred-space (FTCS) scheme applied to the heat equation uₜ = α uₓₓ, von Neumann analysis gives the amplification factor of each mode in closed form, g(θ) = 1 − 4r sin²(θ/2), with θ = kΔx and r = αΔt/Δx² the diffusion number. Every mode must satisfy |g| ≤ 1, and the worst is the sawtooth θ = π (wavelength 2Δx), where sin² = 1. That single worst case fixes the whole stability boundary at r ≤ 1/2 — a hard limit on the time step. Below it the curve stays inside the band and the scheme converges; above it the sawtooth mode has |g| > 1 and grows by a fixed factor each step until the run fills with NaNs.

Figure

Von Neumann amplification factor g(θ) = 1 − 4r sin²(θ/2) for the FTCS heat scheme across the mode angle θ = kΔx. With diffusion number r = 0.4 the whole curve stays in the stable band −1 ≤ g ≤ 1; at r = 0.5 the worst mode θ = π just reaches g = −1; at r = 0.6 it dips to g = −1.4, so that sawtooth mode grows and the scheme diverges. The stability limit is r ≤ 1/2

The core equations

A one-step marching scheme rewrites the PDE as an update rule on grid values. Collecting the spatial stencil into a single discrete operator turns the whole method into repeated multiplication by one evolution operator — the object every stability statement is about.

\frac{u_j^{n+1}-u_j^n}{\Delta t} = \mathcal{L}_\Delta\, u_j^n \;\;\Longleftrightarrow\;\; u^{n+1}=C\,u^n,\quad C=I+\Delta t\,\mathcal{L}_\Delta \qquad (1)

where: $u_j^n$ — numerical value at node $x_j=j\,\Delta x$ , time level $t_n=n\,\Delta t$ ; $\Delta t,\,\Delta x$ — time step and mesh spacing; $\mathcal{L}_\Delta$ — discrete spatial operator approximating $\mathcal{L}$ (say $\partial_x$ or $\partial_x^2$ ); $C$ — one-step evolution operator. Reading: the entire history is $u^N=C^N u^0$ , so everything hinges on the powers of $C$ . Consistency is measured by the residual left when the exact solution is dropped into the scheme. If it vanishes as the mesh is refined, the discrete equation really is approximating the PDE and not some other equation.

\tau_j^n \equiv \frac{u(x_j,t_{n+1})-u(x_j,t_n)}{\Delta t} - \mathcal{L}_\Delta\, u(x_j,t_n) \xrightarrow[\Delta t,\Delta x\to 0]{} 0 \qquad (2)

where: $\tau_j^n$ — local truncation error; $u(x_j,t_n)$ — the true PDE solution sampled at the node. Reading: the largest power of $\Delta t$ or $\Delta x$ surviving in $\tau$ is the order of accuracy; consistency is $\tau\to 0$ , order is how fast. Stability is the demand that no error is amplified without bound over the fixed time interval of interest. In operator form it is a uniform bound on every power of $C$ arising before the final time $T$ .

\big\| C(\Delta t)^n \big\| \le K, \qquad 0\le n\,\Delta t \le T \qquad (3)

where: $C(\Delta t)$ — one-step operator, allowed to depend on $\Delta t$ ; $K$ — a constant independent of $n$ and $\Delta t$ (it may depend on $T$ ); $\|\cdot\|$ — an operator norm on the solution space. Reading: stability forbids $\|C^n\|$ from growing with $n$ ; the bound $K$ is the leash on error amplification. For a constant-coefficient scheme on a periodic grid the Fourier modes diagonalize $C$ , so stability reduces to a scalar test on the per-step amplification factor of each mode. This is the von Neumann condition, the most used tool in the subject.

u_j^n = g(k)^n\, e^{i k j \Delta x} \;\Longrightarrow\; |g(k)| \le 1 + C\,\Delta t \;\;\text{for all }k \qquad (4)

where: $g(k)$ — amplification factor of the mode of wavenumber $k$ ; on the grid $k\,\Delta x\in(-\pi,\pi]$ ; $C$ — a constant independent of $k$ and $\Delta t$ . Reading: the slack $+C\Delta t$ tolerates benign growth (as in reaction terms) but forbids growth surviving the continuum limit; for a dissipative or neutral scheme it tightens to $|g|\le 1$ . For hyperbolic problems the mode test collapses to one geometric inequality — the Courant–Friedrichs–Lewy condition — capping the time step by the fastest characteristic speed.

\nu \equiv \frac{c\,\Delta t}{\Delta x} \le 1 \qquad (5)

where: $\nu$ — the Courant number; $c$ — the characteristic speed of $u_t+c\,u_x=0$ ; for a system, replace $c$ by the largest signal speed $c_{\max}$ . Reading: information must not cross more than one cell per step, or the numerical domain of dependence fails to contain the true one. The Lax equivalence theorem is the payoff: it converts the two checkable properties into the one you cannot check, under stated hypotheses.

\text{consistency} + \text{stability} \;\Longleftrightarrow\; \text{convergence} \qquad (6)

where: the equivalence holds for a linear, well-posed initial-value problem approximated by a consistent scheme; $\Longleftrightarrow$ is genuine logical equivalence. Reading: given consistency, stability is both necessary and sufficient for convergence — drop linearity or well-posedness and only fragments survive. Implicit schemes buy stability by solving a linear system each step; Crank–Nicolson for the heat equation is the canonical example, unconditionally stable for any time step.

g_{\text{CN}}(k) = \frac{1 - 2r\,\sin^2(k\Delta x/2)}{1 + 2r\,\sin^2(k\Delta x/2)}, \qquad |g_{\text{CN}}| \le 1 \;\;\text{for all } r>0 \qquad (7)

where: $r=\alpha\Delta t/\Delta x^2$ , now unrestricted; $g_{\text{CN}}$ — amplification factor of the trapezoidal-in-time scheme. Reading: the numerator never exceeds the denominator in magnitude, so every mode is bounded for every $r$ — the price is a tridiagonal solve, the reward freedom from the $r\le 1/2$ leash.

Derivation AConsistency of the FTCS heat scheme by Taylor expansion

The goal is to prove that the explicit forward-time centred-space scheme for $u_t=\alpha u_{xx}$ is consistent and to extract its order of accuracy, by substituting the exact solution into the discrete operator and Taylor-expanding every term about the node $(x_j,t_n)$ . The temporal difference will contribute a first-order error and the spatial difference a second-order one, and the leading bracket cancels only because the exact solution obeys the PDE. This is the archetype for checking consistency of any finite-difference scheme, and it exposes why order of accuracy is a statement about the truncation error and nothing else.

Step 1. Expand the forward time difference of the exact solution: $\dfrac{u(t_{n+1})-u(t_n)}{\Delta t}=u_t+\tfrac12\Delta t\,u_{tt}+O(\Delta t^2)$ . Why: a one-sided difference in time keeps only the first Taylor term cleanly, leaving a remainder proportional to $\Delta t$ — this is the origin of the scheme's first-order-in-time accuracy.

Step 2. Expand the centred second difference in space: $\dfrac{u_{j+1}-2u_j+u_{j-1}}{\Delta x^2}=u_{xx}+\tfrac{1}{12}\Delta x^2\,u_{xxxx}+O(\Delta x^4)$ . Why: symmetry kills the odd-order remainders, so the centred stencil is second-order accurate — the spatial and temporal orders differ, and the worse of the two governs the scheme.

Step 3. Assemble the truncation error $\tau = (u_t-\alpha u_{xx}) + \tfrac12\Delta t\,u_{tt} - \tfrac{\alpha}{12}\Delta x^2\,u_{xxxx}+\cdots$ and invoke the PDE. Why: the exact solution satisfies $u_t-\alpha u_{xx}=0$ , so the leading bracket vanishes, leaving $\tau=O(\Delta t)+O(\Delta x^2)\to 0$ — consistent, first order in time, second in space.

Step 4. Note the moral: nowhere did stability enter. Why: consistency and order come purely from Taylor remainders of the true solution, so the same expansion holds whether the scheme is stable or not — order can never certify convergence on its own. *Verified: $\tau=\tfrac12\Delta t\,u_{tt}-\tfrac{\alpha}{12}\Delta x^2\,u_{xxxx}+\cdots$ , order $O(\Delta t)+O(\Delta x^2)$ ; matches Strikwerda §3.1 and the base-article expansion.*

Go deeper · Derivation B · 2 minProof of the Lax equivalence theorem

Derivation BProof of the Lax equivalence theorem

The goal is the sufficiency half of the theorem — that for a consistent scheme on a well-posed linear problem, stability forces convergence — proved by a telescoping identity that splits the global error into a propagated initial error plus a sum of propagated local truncation errors. Let $U^n$ be the vector of exact values and $u^n$ the numerical solution, with global error $e^n=u^n-U^n$ . The numerical update is $u^{n+1}=Cu^n$ , while by definition of the truncation error the exact values obey $U^{n+1}=CU^n+\Delta t\,\tau^n$ . The argument is the engine of the whole subject: the precise sense in which a bounded $\|C^n\|$ converts small per-step errors into a small final error.

Step 1. Subtract the two update relations to get the error recursion $e^{n+1}=Ce^n-\Delta t\,\tau^n$ . Why: linearity lets the same operator $C$ act on the error that acts on the solution, so the error obeys the scheme's own dynamics driven by the truncation error as a source.

Step 2. Unroll the recursion into $e^N=C^N e^0-\Delta t\sum_{n=0}^{N-1}C^{N-1-n}\,\tau^n$ . Why: telescoping exhibits the global error as the propagated initial error minus the accumulation of every local error, each aged by the appropriate power of $C$ — this is where "one step raised to the Nth power" becomes explicit.

Step 3. Take norms with the triangle inequality and submultiplicativity: $\|e^N\|\le\|C^N\|\,\|e^0\|+\Delta t\sum_{n=0}^{N-1}\|C^{N-1-n}\|\,\|\tau^n\|$ . Why: this lets a single uniform bound on the powers of $C$ control the entire accumulation at once.

Step 5. Let the mesh shrink: consistency sends $\max_n\|\tau^n\|\to 0$ and matched data sends $\|e^0\|\to 0$ , so $\|e^N\|\to 0$ — convergence. Why: with the amplification bounded once and for all, the only thing left driving the error is the truncation error, which consistency kills; the converse (instability breeds a mode with $\|C^n\|\to\infty$ that destroys convergence) closes the equivalence.

\|e^N\| \le K\,\Big( \|e^0\| + T\,\max_n \|\tau^n\| \Big) \;\xrightarrow[\;\Delta t,\Delta x\to0\;]{}\; 0 \qquad (\star)

Verified: telescoping identity and norm inequalities reproduce Lax and Richtmyer (1956); cf. Strikwerda Thm 10.1.1 and Richtmyer & Morton (1967).

Go deeper · Derivation C · 1 minVon Neumann analysis of the upwind scheme and the CFL limit

Derivation CVon Neumann analysis of the upwind scheme and the CFL limit

The goal is to derive the exact stability window of the first-order upwind scheme for the advection equation $u_t+cu_x=0$ with $c>0$ , watching the Courant–Friedrichs–Lewy condition $0\le\nu\le1$ emerge as an algebraic inequality on the amplification factor. Upwinding differences in the direction the information comes from: $u_j^{n+1}=u_j^n-\nu(u_j^n-u_{j-1}^n)$ with $\nu=c\Delta t/\Delta x$ . We substitute the Fourier mode $u_j^n=g^n e^{ikj\Delta x}$ , compute $|g|^2$ in closed form, and read off the condition for every mode to be non-amplifying — the concrete calculation behind the physical picture's stencil cone.

Step 1. Divide the mode out of the update to get $g=1-\nu(1-e^{-i\theta})=(1-\nu)+\nu e^{-i\theta}$ with $\theta\equiv k\Delta x$ . Why: the mode is an eigenvector of the shift, so the whole scheme reduces to one complex number — a convex combination of $1$ and $e^{-i\theta}$ weighted by $1-\nu$ and $\nu$ .

Step 2. Square the modulus using $|a+be^{-i\theta}|^2=a^2+b^2+2ab\cos\theta$ with $a=1-\nu,\ b=\nu$ : $|g|^2=(1-\nu)^2+\nu^2+2\nu(1-\nu)\cos\theta$ . Why: it is $|g|^2$ that controls per-step growth, and the cross term isolates the only $\theta$ -dependence.

Step 3. Simplify with $(1-\nu)^2+\nu^2=1-2\nu(1-\nu)$ to reach $|g|^2=1-2\nu(1-\nu)(1-\cos\theta)$ . Why: the whole stability question now rides on the sign of $\nu(1-\nu)$ , since $1-\cos\theta\ge 0$ always.

Step 4. Demand $|g|^2\le 1$ for all $\theta$ , i.e. $\nu(1-\nu)\ge 0$ , giving $0\le\nu\le 1$ . Why: precisely the CFL condition — both weights $1-\nu$ and $\nu$ must be non-negative, geometrically the statement that the characteristic foot lands inside the cell the stencil reads; upwinding takes its data from the physically correct side.

|g|^2 = 1 - 2\nu(1-\nu)\,(1-\cos\theta) \le 1 \;\Longleftrightarrow\; 0\le\nu\le1,\qquad \theta=k\Delta x \qquad (\star\star)

Verified: reduces to the Courant–Friedrichs–Lewy bound (1928); $\nu=0.8,\theta=\pi$ gives $|g|^2=0.36$ , matching direct evaluation; cf. Strikwerda §1.3.

Limiting cases

The four rows map out where the equivalence bites, where it breaks, and what each failure of the trio looks like.

Case	Result	Meaning
Lax regime	linear, well-posed, consistent	Stability ⇔ convergence. Check consistency by Taylor, stability by Fourier, and convergence is guaranteed.
Consistent but unstable	$\lvert g(k)\rvert > 1 + C\,\Delta t$	e.g. FTCS advection, $\lvert g\rvert^2=1+\nu^2\sin^2\theta$ . Convergence fails; short modes grow like $\lvert g\rvert^n$ regardless of order.
Stable but inconsistent	$\tau \not\to 0$	Converges to the wrong equation (a "modified equation" with spurious diffusion or dispersion) — reproducible, but not the problem you posed.
Nonlinear / variable-coefficient	$\partial_t u + \partial_x f(u)=0$	Lax theorem does not apply. Von Neumann gives only a necessary frozen-coefficient condition; genuine guarantees need energy estimates, the Kreiss theorem, or entropy/TVD bounds.

Common trapOrder of accuracy proves stability

It proves nothing of the kind, and believing it wrecks simulations. Consistency and its order come from a Taylor expansion of the exact solution; stability comes from a Fourier analysis of error growth — the two calculations never touch. The forward-time centred-space scheme for advection is the standing counterexample: its truncation error is $O(\Delta t)+O(\Delta x^2)$ , formally respectable, yet $g=1-i\nu\sin\theta$ gives $|g|^2=1+\nu^2\sin^2\theta>1$ for every non-trivial mode and every time step. It is unconditionally unstable — no choice of $\Delta t$ saves it. Always test $\sup_k|g(k)|$ separately; never read stability off the order.

Worked examples

1 · The FTCS stability leash on a real diffusion problem (fully worked)

A copper bar (thermal diffusivity α = 1.11×10⁻⁴ m²/s) is simulated with the explicit FTCS heat scheme on a grid of spacing Δx = 1.0 mm. What is the largest stable time step, and what does refining the mesh cost?

Step 1. The stability condition from the amplification factor is $r=\alpha\Delta t/\Delta x^2\le\tfrac12$ , so the maximum time step is $\Delta t_{\max}=\Delta x^2/(2\alpha)$ . Insert the numbers: $\Delta t_{\max}=(10^{-3})^2/(2\times1.11\times10^{-4})=10^{-6}/2.22\times10^{-4}$ . Δt_max = 4.5×10⁻³ s (about 4.5 ms).

Step 2. Now refine the space grid to resolve finer structure, halving Δx to 0.5 mm. Because the bound scales as $\Delta x^2$ , the allowed step falls by a factor of four: $\Delta t_{\max}=(0.5\times10^{-3})^2/(2\times1.11\times10^{-4})=2.5\times10^{-7}/2.22\times10^{-4}$ . Δt_max = 1.13×10⁻³ s (about 1.13 ms).

Step 3. Count the cost. Halving Δx doubles the spatial nodes and quarters the time step, so a fixed final time takes four times as many steps over twice as many points — an eightfold increase in work for one mesh halving. This Δt ∝ Δx² coupling is exactly why stiff diffusion problems migrate to implicit Crank–Nicolson, stable for any r at the price of one tridiagonal solve per step. Stability, not accuracy, sets the time step here.

2 · The CFL time step for a tsunami simulation

A trans-oceanic tsunami is modelled by the shallow-water wave speed c = √(gh) over an ocean of depth h = 4000 m (g = 9.81 m/s²), on a grid of spacing Δx = 10 km. With an explicit scheme run at a safety Courant number ν = 0.8, what time step is allowed?

Step 1. Compute the fastest signal speed: $c=\sqrt{gh}=\sqrt{9.81\times4000}=\sqrt{39240}$ . c ≈ 198 m/s (about 710 km/h — jetliner speed).

Step 2. The CFL condition caps the step at $\Delta t\le\Delta x/c$ ; running at ν = 0.8 gives $\Delta t=\nu\,\Delta x/c=0.8\times10\,000/198$ . Δt ≈ 40 s. A full step advances the wave $c\Delta t=198\times40\approx7.9$ km, which is 0.79 of a cell — safely inside one grid spacing, as ν < 1 demands.

Step 3. Interpret it: a wave crossing a 10 000 km ocean at 198 m/s takes about 14 hours, so the run needs roughly $14\times3600/40\approx1300$ time steps. Push ν above 1 to save steps and the physical characteristic escapes the stencil cone — the scheme goes unstable and the forecast fills with noise. By contrast the same FTCS-advection stencil run centred instead of upwind has $|g|^2=1+\nu^2\sin^2\theta>1$ and blows up at any ν: same PDE, same order of consistency, opposite fate — the difference is entirely stability.

Physics of the portrait

The equivalence at the heart of this article carries the name of the mathematician who, with Robert Richtmyer, turned computational physics from an art into a theory.

Origins · 1 minPeter Lax and how this was discovered

Peter Lax

1926–2025 · numerical analysis and the theory of hyperbolic conservation laws

Peter Lax was born in Budapest in 1926 and reached New York in 1941, a prodigy known to the emigre mathematical community before he finished high school. As a teenage soldier he was assigned to the Manhattan Project at Los Alamos, where turning the equations of shock waves into numbers on the first electronic computers gave him a lifelong conviction that analysis and computation are one subject. He spent his career at the Courant Institute of New York University and, more than anyone, built the modern theory of solving time-dependent partial differential equations numerically. With Richtmyer in 1956 he proved the equivalence theorem that bears his name — for consistent schemes on well-posed linear problems, stability and convergence are the same thing — giving practitioners licence to certify a method they could never test directly. With Burton Wendroff he devised the Lax–Wendroff and Lax–Friedrichs schemes that made shock capturing practical; his work on hyperbolic conservation laws introduced the entropy conditions that select the physically correct weak solution when the linear theory fails, and his discovery of Lax pairs revealed the integrable structure behind solitons. He received the Abel Prize in 2005 and remained active nearly to his death in 2025. His signature idea runs through this article: a scheme earns trust not by looking accurate but by provably not amplifying error.

History

1928

— Courant, Friedrichs, and Lewy publish the condition that bounds an explicit time step by the mesh

1947

— Crank and Nicolson introduce their unconditionally stable implicit scheme for diffusion

1950

— von Neumann's Fourier stability method reaches print in the Charney–Fjortoft–von Neumann weather-prediction paper

1956

— Lax and Richtmyer prove that consistency plus stability is equivalent to convergence

1967

— Richtmyer and Morton's monograph systematizes difference methods for initial-value problems

Going further

Texts. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed. (2004), Chs. 1, 2, 10 · Richtmyer & Morton, Difference Methods for Initial-Value Problems, 2nd ed. (1967) · LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (2007) · Stone & Goldbart, Mathematics for Physics, numerical-methods chapter · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed. · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed.

Foundational papers. Courant, Friedrichs & Lewy, Math. Ann. 100, 32 (1928) · Lax & Richtmyer, "Survey of the stability of linear finite difference equations", Comm. Pure Appl. Math. 9, 267 (1956) · Lax & Wendroff, "Systems of conservation laws", Comm. Pure Appl. Math. 13, 217 (1960).

Provenance. Every stability result reproduces by inserting $u_j^n=g^n e^{ikj\Delta x}$ : FTCS heat gives $g=1-4r\sin^2(\theta/2)$ and $r\le1/2$ ; upwind gives $|g|^2=1-2\nu(1-\nu)(1-\cos\theta)$ and $0\le\nu\le1$ ; FTCS advection gives $|g|^2=1+\nu^2\sin^2\theta>1$ ; Crank–Nicolson gives $|g|\le1$ for all $r$ . The convergence bound follows by unrolling $e^{n+1}=Ce^n-\Delta t\,\tau^n$ and applying $\|C^m\|\le K$ . Worked numbers checked independently: copper Δt_max = 4.5 ms (α = 1.11×10⁻⁴ m²/s, Δx = 1 mm); tsunami c = 198 m/s, Δt ≈ 40 s (ν = 0.8, Δx = 10 km). Framing follows Strikwerda, Richtmyer & Morton, and Stone & Goldbart.

Prerequisites. Order Symbols and Gauge Functions · Asymptotic Series and Poincaré Expansions · Fourier analysis and norms · Well-posedness of PDEs. Leads to. WKB Approximation · Optimal Truncation and Superasymptotics · Spectral and pseudospectral methods.

Ask the tutor: Ask me to run von Neumann analysis on any scheme you name — I will insert the Fourier mode, hand you the amplification factor, and turn it into the exact stability condition on your time step.

What is the difference between stability and convergence — aren't they the same idea?

They look alike but are independent, and only the Lax theorem ties them together. Stability is about error growth: repeated application of the update operator does not amplify small perturbations without bound, $\|C^n\|\le K$ . Convergence is about accuracy: the computed solution approaches the true one as the mesh shrinks. A scheme can be stable yet converge to the wrong equation (if inconsistent), or consistent — perfectly accurate per step — yet diverge (if unstable). The equivalence theorem proves that once you have consistency the two become one: stability is exactly what upgrades a consistent scheme to a convergent one, because it keeps the accumulated per-step errors from exploding.

Why can't I just refine the mesh and watch whether the answer converges?

You can — it is called a convergence study — but it is empirical, not a proof, and it can mislead. If a scheme is unstable, refining the mesh usually makes things worse: a finer grid supports shorter, more strongly amplified modes and often forces a smaller time step whose round-off feeds the instability. If a scheme is stable but inconsistent, refinement produces a beautifully convergent sequence heading toward the wrong answer. And a study only samples a few meshes over a finite time, so it cannot certify the limit. The Lax theorem is valuable because it replaces this fragile experiment with two analyses — Taylor for consistency, Fourier for stability — that together guarantee the limit you never actually reach.

Does the Lax equivalence theorem work for nonlinear problems?

Not directly, and this is the most important caveat. The proof leans on linearity — the error obeys the same operator $C$ as the solution, which makes the telescoping identity work — and on well-posedness, a bounded solution operator on a Banach space. For a nonlinear conservation law $u_t+\partial_x f(u)=0$ neither holds globally: shocks form, solutions are non-unique without an entropy condition, and the operator depends on the state. Von Neumann analysis survives only as a necessary local test on the frozen-state linearization; it flags instability but cannot certify convergence. Genuine nonlinear guarantees come from energy methods, the Kreiss matrix theorem, and the total-variation-diminishing and entropy-stability frameworks — much of which grew out of Lax's own later work.

What does the Courant number ν physically mean, and why exactly one?

The Courant number ν = cΔt/Δx compares the physical signal speed c with the grid speed Δx/Δt at which an explicit stencil passes information; equivalently, how far a signal travels in one step versus one cell width, with ν = 1 where they match. Below it the true disturbance stays within the cell the scheme reads each step, so the numerical domain of dependence contains the physical one and the update has the data it needs. Above it the disturbance outruns the stencil: the point's true value depends on data never sampled, and no arithmetic on the wrong inputs gives the right output — the scheme goes unstable. That is the CFL condition, necessary for any explicit scheme regardless of the discretization.

Why is von Neumann analysis restricted to constant coefficients and periodic domains?

Because it is diagonalization in disguise, and only there are the Fourier modes exact eigenvectors of the update operator. On a periodic constant-coefficient grid the stencil shifts $e^{ikx}$ into a pure multiple of itself, so the operator acts mode-by-mode as multiplication by the scalar $g(k)$ and stability collapses to a scalar inequality. Variable coefficients couple the modes; boundaries break the exactness of the discrete Fourier basis. In practice one still applies von Neumann locally, freezing the coefficients at each point — a reliable necessary condition and usually the right time-step restriction, but it can miss instabilities driven by boundaries or coefficient variation. Those need the matrix (Kreiss) theory or an energy estimate.

My simulation suddenly blew up with NaNs. How do I diagnose it?

Look first at the wavelength of what is growing. A classic numerical instability appears as oscillations at the shortest representable wavelength, 2Δx (the sawtooth mode θ = π), growing exponentially — because for most explicit schemes the worst amplification factor sits at high wavenumber. Check your time step against the relevant limit: r = αΔt/Δx² ≤ 1/2 for explicit diffusion, ν = cΔt/Δx ≤ 1 for explicit advection, using the fastest local signal speed. If a wave speed rose during the run — a shock steepened, a temperature spiked — a step stable at first can cross the threshold mid-run. Fixes: shrink Δt, add controlled artificial viscosity, or switch to an implicit or upwind scheme. The unifying diagnosis is always the same — compute sup_k |g(k)| and confirm it exceeds one.