Stability of a numerical method controls how perturbations in initial data, roundoff and local residuals propagate. For a linear evolution scheme , a suitable finite-time requirement is
with independent of the refining mesh. The norm should approximate the norm of the continuous problem, such as an area-weighted L2 norm. The bound need not imply decay: is appropriate when the equation itself permits finite-time growth. For a multilevel scheme the state includes all required time levels.
The error recurrence shows why stability matters. If , then
Thus a consistent residual tending to zero and a convergent start produce convergence when the stability constant is uniform. Tiny defects need not remain tiny without that bound. The Lax equivalence theorem makes this precise: for a well-posed linear initial value problem and a consistent linear approximation, stability is equivalent to convergence in the compatible solution norms. It is not a theorem for arbitrary nonlinear schemes, and a high formal order does not substitute for stability.
For constant coefficients on a periodic or infinite grid, von Neumann stability analysis diagonalizes spatial translations using Fourier modes. If , the Parseval identity converts a uniform multiplier-power bound into an L2 norm bound. The contraction condition is sufficient; a more general also yields a finite-time bound. For the centered-space explicit heat method,
so exactly for . Larger fixed Courant numbers amplify high frequencies. The Backward Euler method gives instead and is contractive for every .
For , forward time and centered space give , with . Hence for nonzero frequencies, causing instability at a fixed nonzero Courant number. The directionally correct upwind finite difference scheme uses the neighbor at and has
It is contractive exactly for . This also has a direct maximum-norm proof: the update is a convex combination of two old values. The continuum domain of dependence must fit inside the numerical one for a local explicit hyperbolic method to converge, giving the Courant–Friedrichs–Lewy condition. That necessary geometric restriction is not by itself a sufficient stability proof.
For systems and multilevel methods, each Fourier multiplier becomes a matrix. All amplification roots must obey the appropriate root bound, but their moduli alone are insufficient. A repeated unit root in a companion matrix gives a Jordan block; for instance has . Even distinct roots require control of the diagonalizing matrices uniformly in frequency and mesh. The uniform power bound from separated amplification roots supplies such control for a bounded two-root family with a positive gap. Question5 illustrates both sides: strict gives a uniform gap, whereas the endpoints have a double unit root and are unstable despite unit moduli.
Energy estimates are especially useful when coefficients vary or boundaries prevent Fourier diagonalization. As a simple example, for implicit heat evolution with homogeneous Dirichlet boundary conditions, take the grid inner product of with . Discrete integration by parts gives
The boundary term vanishes because the boundary values are zero, proving contraction. For variable positive diffusion coefficients the corresponding flux-weighted squared differences remain nonnegative, so the argument extends without a constant-coefficient Fourier symbol. With forcing or lower-order terms, a discrete energy inequality and Gronwall inequality yield a mesh-independent finite-time bound.
A related technique is eigenvalue stability analysis of a finite difference method. If the spatial matrix is symmetric with nonpositive eigenvalues and the time update is , an orthonormal eigenbasis gives
This proves stability when all scaled spatial eigenvalues lie in the time method's stability domain. For a nonnormal spatial or amplification matrix, the eigenvalues alone do not control powers: ill-conditioned eigenvectors or Jordan terms can cause transient growth. Direct operator-norm estimates, symmetrizers or energy estimates then provide the missing information. In particular an A-stable scalar time formula does not automatically give mesh-uniform bounds for every ill-conditioned family of spatial matrices.
Boundary conditions and starting procedures are part of the discretization. Interior Fourier analysis may fail to detect a growing boundary mode, and a multilevel method can excite a parasitic mode through a poor start. One must therefore analyze the actual finite-domain operator or include its boundary terms in an energy proof. Stability is a uniform error-propagation bound; Fourier symbols, energy identities, matrix norms and spectral analysis are complementary ways of establishing that bound.
The stability convention matters. The Cambridge numerical analysis notes define a strict linear stability domain through decay . For a multistep recurrence, robust decay requires every amplification root to lie strictly inside the unit disk. With that convention we can prove the stronger result that this method's strict linear stability domain is empty, and therefore bounded. Mere boundedness of oscillatory sequences gives a different answer, described below.
On the Dahlquist test equation, . With , the amplification polynomial of a multistep method is
Exclude , where the intended two-step update cannot be solved. If both polynomial roots have modulus below one, their product requires . But
so necessarily . The complex quadratic Schur criterion also requires
Its left-hand side equals , while for its right-hand side is only . The inequality is impossible. Thus
For comparison, the bounded-root stability set of the symmetric two-derivative formula is unbounded. On , rotate the amplification polynomial into a real reciprocal quadratic. Its two roots are distinct and on the unit circle precisely when
equivalently or . In particular, arbitrarily large imaginary parameters give bounded, undamped numerical oscillations. At the two equality thresholds, repeated unit-circle roots produce unbounded solutions. Hence the question's bounded-domain claim holds for the strict decay convention; it does not hold if the phrase is interpreted using the non-growing root condition for a multistep method.
The strict decay convention for a linear stability domain demands that every numerical mode tends to zero on the Dahlquist test equation. For a one-step stability function it requires ; for a multistep method it requires every amplification root to have modulus below one. This differs from a boundedness convention that allows simple unit-circle roots, and from mesh-uniform finite-time numerical stability. Cambridge numerical analysis notes use the decay convention. Explicitly specifying the convention is essential for methods with undamped modes.
For , write an amplification root as . Its equation is . The roots have modulus one and a positive frequency-uniform separation exactly for , giving L2 norm stability through a uniformly bounded diagonalization. At either endpoint the Nyquist symbol has a repeated unit root and a Jordan block, producing linear growth. The normalized local truncation error is generally second order; at and exactly sampled characteristics solve the recurrence, although the latter parameter is unstable for arbitrary perturbations.