The Parabolic Step Limit in Explicit Diffusion Schemes

The heat equation $u_t = \alpha u_{xx}$ is the simplest parabolic problem, and every explicit discretization of it is subject to a stability constraint whose form is unlike the CFL condition for advection. The step must scale as the square of the spatial spacing. That single dependence explains most of the practical differences between explicit and implicit approaches to diffusion, because halving the mesh costs a factor of four in time steps rather than a factor of two.

The scheme and its amplification factor

Discretize with forward time and centered space:

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

Writing $r = \alpha \Delta t/\Delta x^2$ and substituting a Fourier mode $e^{ikx_j}$, the second difference produces the factor $2(\cos\theta - 1) = -4\sin^2(\theta/2)$ with $\theta = k\Delta x$, so

$$ g(\theta) = 1 - 4r\sin^2(\theta/2) $$

Stability requires $|g| \le 1$ for all $\theta$. The term $4r\sin^2(\theta/2)$ ranges over $[0, 4r]$, so the condition is $-1 \le 1 - 4r \le 1$, which reduces to $r \le 1/2$. The bound is familiar and its form follows from the largest possible value of the second difference: at $\theta = \pi$, adjacent values alternate, and the difference is the extreme $4$.

The choice $r = 1/2$ is special. It makes the amplification factor $g(\theta) = \cos\theta$, which is exactly what the mode’s exact evolution gives in the long-wavelength limit but also has a specific property: $g(\pi) = -1$, so the highest-frequency mode is not damped, merely flipped in sign each step. Any initial data with a grid-scale component keeps a parasitic oscillation of constant amplitude. Choosing $r$ slightly below $1/2$ damps that mode, and $r = 1/6$ makes $g \ge 0$ for all $\theta$, eliminating sign alternation entirely. Schemes whose amplification factor is nonnegative are sometimes preferred for this reason, at the cost of a step smaller than the stability limit.

Reading the bound as a physical constraint

The constraint $r \le 1/2$ states that the diffusion length over one step cannot exceed the mesh spacing by much. A random walk of a diffusing quantity spreads a distance of order $\sqrt{2\alpha t}$, so the condition is that the rms spread per step stay comparable to $\Delta x$. The exact statement is $\sqrt{\alpha \Delta t} \lesssim \Delta x/\sqrt{2}$, and the bound is the discrete expression of the fact that explicit time stepping samples the diffusion process at intervals short enough to see the quantum of spreading that a mesh cell represents.

Refining the mesh is therefore quadratically expensive in time. Halving $\Delta x$ in one dimension quarters the step, so the work grows by a factor of eight: the point count doubles while the step count quadruples. Refining in three dimensions multiplies the point count by eight and quarters the step, a factor of 32. This is the arithmetic that makes explicit integration of a diffusion-dominated problem unattractive at fine resolution, and it is the reason the parabolic case differs so sharply from CFL, where the cost of refining is linear in the cell count.

What the same analysis looks like for the eigenvalues

The von Neumann result has a matching statement in terms of the discretization matrix. The centered second-difference operator on $N$ interior points with Dirichlet boundaries has eigenvalues

$$ \lambda_k = -\frac{4\alpha}{\Delta x^2}\sin^2\left(\frac{k\pi}{2(N+1)}\right), \quad k = 1, \ldots, N $$

whose extremes are slow, of order $-\alpha \pi^2/L^2$ for a domain of length $L$, and fast, of order $-4\alpha/\Delta x^2$. The ratio of extremes is the stiffness ratio of the semi-discrete system, and it scales as $1/\Delta x^2$. An explicit method’s stability region must contain all of them, and the fastest sets the step. This is stiffness generated by discretization rather than by the model, and the general treatment of stiffness in stiffness: when explicit steps collapse applies directly.

The eigenvalue view adds information the Fourier analysis does not give for finite domains: the slowest mode sets the time scale over which the solution actually changes, and it is a factor $1/\Delta x^2$ away from the fastest. For a grid fine enough to resolve a boundary layer, that gap can be enormous, and an explicit method integrates the interesting dynamics with a step set by the least interesting mode in the system.

Implicit treatment and its structure

Backward Euler applied to the same discretization gives

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

a linear system with a tridiagonal matrix in one dimension. The amplification factor is $g = 1/(1 + 4r\sin^2(\theta/2))$, which is positive, bounded by 1, and tends to 0 as $r \to \infty$: the method is L-stable and damps the grid-scale mode completely. The step is limited only by accuracy.

The cost structure is what decides the trade. In one dimension the tridiagonal solve is $O(N)$ per step, so the implicit step is a small multiple of the explicit one and the elimination of the stability bound is almost free: any factor by which the step grows is a factor of work saved. In two dimensions the five-point matrix admits banded factorization whose cost grows like $N^{3/2}$ for an $N$-point grid versus $N$ for the explicit sweep, so implicit costs more per step and wins only when the step grows by more than the ratio. In three dimensions the gap widens further, with a dense or sparse direct factorization growing like $N^2$ or $N^{7/3}$ depending on ordering, and iterative methods or operator-splitting schemes take over. The scaling of implicit solves with problem structure is the subject of Newton iterations for implicit operators, which covers the nonlinear case where the linear solve is nested inside an outer iteration.

The Crank-Nicolson scheme, which averages the explicit and implicit operators, $g = (1 - 2r\sin^2(\theta/2))/(1 + 2r\sin^2(\theta/2))$, is A-stable and second order but not L-stable: as $r \to \infty$ the factor tends to $\cos\theta$, so the stiffest mode oscillates rather than decays. On smooth initial data this is invisible. On data with a discontinuity or after a change in the boundary condition it produces a lasting oscillation. This is the distinction explained in A-stability and the Dahlquist barrier, and it is a common source of a solution that looks stable but has a decaying ripple the user did not expect.

A middle path is the so-called IMEX or semi-implicit treatment: integrate the diffusion implicitly and the reaction or advection explicitly, computing a single linear solve per step with the same matrix coefficients. For constant diffusion the matrix is fixed, so its factorization is computed once and reused for every step, and the per-step cost becomes a triangular solve.

Anisotropy and the directional sum

For an anisotropic or multidimensional problem the bound becomes a sum over directions. On a Cartesian grid with spacings $\Delta x_1, \ldots, \Delta x_d$ and diffusivities $\alpha_1, \ldots, \alpha_d$, the explicit scheme with the standard stencil requires

$$ \sum_{m=1}^{d} \frac{2\alpha_m \Delta t}{\Delta x_m^2} \le 1 $$

A single thin cell in one direction therefore bounds the step for the whole problem, even though the solution may be essentially one-dimensional elsewhere. This is a discretization artifact with a concrete fix: treating a direction implicitly, or using an alternating-direction scheme that applies one-dimensional implicit solves in sequence, removes that direction’s constraint from the sum.

The bound also justifies grid design. Since the step scales with the square of the smallest spacing, a mesh refined by a factor of three near a boundary costs nine times the steps for the same physical time, whereas a graded mesh that coarsens away from the layer limits the penalty. The interaction between mesh design and step selection is a joint optimization, not two independent decisions.

Failure signatures in practice

An explicit diffusion code that violates $r \le 1/2$ does not merely become inaccurate. The grid-scale mode’s amplification factor passes $-1$ and then exceeds 1 in magnitude, and the computed solution develops a checkerboard pattern whose amplitude doubles roughly every step. The pattern appears first in the highest-frequency content of the initial data or at a boundary where the condition changes, and spreads as the instability feeds on itself. Because the growth is per step and per mode, the failure can be invisible for dozens of steps if the initial data is smooth and the odd-even modes are weakly excited, and then appear suddenly.

A useful diagnostic is to compute, each step, the discrete energy $\sum_j u_j^2 \Delta x$ and compare its decay to the exact rate. For the continuous problem the energy decays at a rate given by the Dirichlet form; a discrete energy that grows, or that decays more slowly than the physical rate by a large factor, indicates that the grid-scale modes are being amplified rather than damped. The check costs one pass over the array and catches the failure while the solution is still recognizably a solution.

The other class of failure is more subtle. A scheme with $r$ just below the limit has amplification factors near $-1$ for the grid-scale modes, so short-wavelength noise persists instead of being damped. If the initial data or the boundary conditions contain such noise, the result is a solution whose smooth part is correct while a high-frequency haze rides on top. Tightening the tolerance and refining the grid both reduce the haze’s amplitude but do not remove the mechanism; only a scheme that damps the grid-scale modes removes it, which is the same requirement the L-stability discussion isolates.