Global Error Accumulation and the CFL Condition

Two error analyses in numerical differential equations look superficially similar and answer different questions. The first asks how fast the computed solution converges to the exact one as the grid is refined; it is a statement about accuracy. The second asks whether the computed solution stays bounded at all; it is a statement about stability. A scheme can be consistent, with a truncation error that vanishes as the grid refines, and still produce divergent output at every step size, because consistency and stability are independent properties. The classic example is the explicit central-difference scheme for the advection equation, and walking through it exposes why the Courant number is the quantity that matters.

A discrete grid and a discrete equation

Consider the linear advection equation with constant speed $c > 0$:

$$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 $$

Discretize space on a uniform grid with spacing $\Delta x$ and time with spacing $\Delta t$, writing $u_j^n$ for the value at $x_j = j\Delta x$, $t^n = n\Delta t$. Replace the time derivative by a forward difference and the space derivative by a centered difference:

$$ u_j^{n+1} = u_j^n - \frac{c,\Delta t}{2,\Delta x}\left(u_{j+1}^n - u_{j-1}^n\right) $$

Insert the exact solution into the same expression and expand in Taylor series. The forward time difference is $u_t + (\Delta t/2)u_{tt} + \cdots$ and the centered space difference is $u_x + (\Delta x^2/6)u_{xxx} + \cdots$. The first derivatives cancel against the PDE, leaving a residual of order $\Delta t$ from the time term and $\Delta x^2$ from the space term. The scheme is consistent, and its local truncation error vanishes under refinement. On those grounds alone it looks acceptable.

Running it reveals otherwise. It is unconditionally unstable for this version of the update, and giving the update a forward space difference instead of a centered one restores stability only under a restriction on $\Delta t / \Delta x$.

Exponential solutions and the amplification factor

Stability is tested by asking what happens to a single Fourier mode. Since the equation is linear and the grid is uniform, any solution can be written as a sum of modes $u_j^n = \hat u^n e^{i k x_j}$, where $k$ is the wavenumber. Substituting this into the update gives the amplification factor $g(k)$, defined by $\hat u^{n+1} = g,\hat u^n$. For the centered scheme:

$$ g(\theta) = 1 - i,\nu \sin\theta, \qquad \nu = \frac{c,\Delta t}{\Delta x},\quad \theta = k\Delta x $$

and therefore

$$ |g|^2 = 1 + \nu^2 \sin^2\theta $$

The magnitude is at least 1 for every wavenumber, and strictly greater than 1 for every non-trivial one. Each step multiplies the amplitude by more than one, so after $n$ steps the amplitude grows as $(1 + \nu^2 \sin^2\theta)^{n/2}$. The scheme is unstable for every choice of $\Delta t$ and $\Delta x$. Refinement makes the growth factor closer to 1 but never below it. The number of steps also increases as the grid refines, and the two effects conspire to leave the computed solution unusable. This is the sharp version of the point that consistency does not imply convergence.

The von Neumann analysis extends to systems and to multi-level schemes by treating the vector of values at each grid point and looking for the spectral radius of the amplification matrix. The condition for stability is $\rho(g(\theta)) \le 1$ for all real $\theta$, uniformly in the parameters. For discretizations of problems with variable coefficients, the analysis is treated as a local statement frozen at each point, which is what the standard Lax equivalence theorem assumes when it states that consistency plus stability implies convergence.

The upwind fix and where the Courant number appears

Take the same PDE with a one-sided difference in the direction the information travels. For $c > 0$ the characteristic moves rightward, so the value at $x_j$ should be influenced by $x_{j-1}$:

$$ u_j^{n+1} = u_j^n - \nu\left(u_j^n - u_{j-1}^n\right) $$

whose amplification factor is $g(\theta) = 1 - \nu(1 - e^{-i\theta})$, and therefore

$$ |g|^2 = 1 - 2\nu(1 - \nu)(1 - \cos\theta) $$

The factor $|g|^2 \le 1$ for all $\theta$ exactly when $0 \le \nu \le 1$. That is the Courant-Friedrichs- Lewy condition for this scheme: the time step cannot exceed the time for information to cross one cell. The condition has a geometric reading that predates the algebra. The numerical domain of dependence of $u_j^{n+1}$ is the interval between $x_{j-1}$ and $x_j$; the physical domain of dependence is the single point $x_j - c\Delta t$. For the scheme to have any chance of representing the solution, the physical domain must lie inside the numerical one, which for this stencil means $0 \le c\Delta t \le \Delta x$.

The same constraint generalizes in the multidimensional case to a sum over directions. For advection in $d$ dimensions with speeds $c_d$ on a Cartesian grid, the bound is $\sum_d |c_d|\Delta t/\Delta x_d \le 1$ for the simplest explicit upwind scheme, so adding dimensions tightens the step. For parabolic problems, where the appropriate quantity involves $\Delta t / \Delta x^2$ rather than $\Delta t/\Delta x$, the consequence is more severe and is the subject of explicit diffusion and the parabolic step limit.

The upwind scheme’s exact stability boundary at $\nu = 1$ has an interpretation worth stating: at that value the update becomes $u_j^{n+1} = u_{j-1}^n$, which is exact translation by one cell. The scheme reproduces the exact solution at grid points when the departure point of a characteristic lands on a grid node.

The upwind scheme also exposes a cost that the centered scheme hides. Its first-order truncation error is a diffusive term with coefficient $c,\Delta x(1 - \nu)/2$, so the scheme solves a modified equation with an added viscosity proportional to the mesh spacing. On a smooth solution the effect is to smear features over a width that shrinks with $\Delta x$, and the smear is why a first-order scheme’s error is $O(\Delta x)$ in smooth regions rather than the $O(\Delta x^2)$ that the formal truncation analysis of the centered stencil suggests. The remedy, a higher-order reconstruction with a limiter that keeps the profile from oscillating; is a separate design with its own constraints, covered in limiters, TVD conditions, and shock capture.

Numerical dispersion, which stability does not prevent

Stability says the amplitude stays bounded; it says nothing about the phase. For the upwind scheme with $\nu$ between 0 and 1, the amplification factor’s modulus is below 1 for most wavenumbers, so short wavelengths decay even though the exact advection equation does not dissipate anything. Separately, the phase of $g$ differs from the exact phase $e^{-i\nu\theta}$, so different wavenumbers travel at different speeds. The error is largest for the shortest resolved wavelengths, at $\theta$ near $\pi$, where the grid has only two points per wavelength.

This is the origin of a familiar artefact: a sharp profile with a jump is advected as a wave train, with oscillations trailing or leading the discontinuity depending on the scheme. Refining the grid moves the artefact to shorter wavelengths and reduces its amplitude, but a first-order upwind scheme adds numerical viscosity of order $c,\Delta x(1 - \nu)/2$, which smooths a step over several cells no matter how fine the grid. The remedy is a higher-order reconstruction, traded against the stencil width and the stability bound, a trade examined in limiters, TVD conditions, and shock capture.

Implementing the bound as a guard

The CFL condition is a hard constraint on the step, and code that computes the step from a target accuracy must clamp it. The clamp belongs at the point where the step is chosen, not in a comment:

/* Advect with first-order upwind, c > 0. Returns steps taken. */
size_t advance(double *u, size_t n, double dx, double c, double t_end,
               double dt_requested) {
    const double nu_max = 0.9;                 /* safety factor below 1 */
    double dt_max = nu_max * dx / c;
    double dt = dt_requested < dt_max ? dt_requested : dt_max;
    size_t steps = 0;
    double *u_new = malloc(n * sizeof *u_new);
    for (double t = 0.0; t + 0.5 * dt < t_end; t += dt, steps++) {
        u_new[0] = 0.0;                         /* inflow boundary */
        for (size_t j = 1; j < n; j++)
            u_new[j] = u[j] - (c * dt / dx) * (u[j] - u[j - 1]);
        memcpy(u, u_new, n * sizeof *u);
    }
    free(u_new);
    return steps;
}

The safety factor below 1 is not decoration. At exactly $\nu = 1$ the scheme is stable but any rounding error that pushes an effective coefficient above 1 flips the growth rate positive, and nonlinear problems localize such effects at steep gradients. A margin also leaves room for the rounding in dt itself: accumulating t by repeated addition drifts from the nominal grid, and the inequality must hold for the actual coefficients used, not the intended ones.

Boundary conditions decide stability independently of the interior

The von Neumann analysis assumes a uniform infinite grid. Real domains have boundaries, and the interior bound is necessary but not sufficient. An outflow boundary treated with a naive condition that reflects can make an interiorly stable scheme unstable, and the reflection appears as a mode growing at the boundary rather than in the interior. The standard practice is to test the boundary closure separately, either by a localized normal-mode analysis or by checking that the discrete energy is non-increasing when the boundary terms are included.

For systems, the count of boundary conditions comes from the characteristic decomposition: the number of incoming characteristics at each boundary decides how many conditions the continuous problem admits, and a discrete scheme must match that count. Supplying one condition too few leaves the boundary underspecified and produces a scheme with a null mode that drifts; one too many overconstrains the solution and generates a spurious reflected wave. Neither failure shows up in a consistency check, and both show up immediately as the computed solution departs from the exact one after the wave reaches the edge.

Taken together, the two error concepts give the practical rule for choosing a scheme: truncation error sets the accuracy at a fixed grid, stability sets the region of grid parameters in which the truncation analysis has any meaning, and the two constraints meet at the boundary of the stable region where accuracy is usually best per unit of work.