Past exam of the mathematics course of the University of Cambridge / 2002 / iii / Paper 60 / 7 / Solution

Stability of a numerical method controls how perturbations in initial data, roundoff and local residuals propagate. For a linear evolution scheme $U^{n+1}=S_h{U}^n+k{F}^n$, a suitable finite-time requirement iswith $C_T$ 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: $C_T=e^{CT}$ 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 $e^{n+1}=S_h{e}^n+k{\tau }^n$, thenThus 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 ${\hat{U}}^{n+1}=g(\theta ){\hat{U}}^n$, the Parseval identity converts a uniform multiplier-power bound into an L2 norm bound. The contraction condition $|g(\theta )|\le 1$ is sufficient; a more general $|g(\theta )|\le 1+Ck$ also yields a finite-time bound. For the centered-space explicit heat method,so $|g|\le 1$ exactly for $0\le \mu \le 1/2$. Larger fixed Courant numbers amplify high frequencies. The Backward Euler method gives instead $g=[1+4\mu {\sin }^2(\theta /2){]}^{-1}$ and is contractive for every $\mu \ge 0$.

For $u_t=u_x$, forward time and centered space give $g=1+i\nu \sin \theta$, with $\nu =k/h$. Hence $|g{|}^2=1+{\nu }^2{\sin }^2\theta >1$ for nonzero frequencies, causing instability at a fixed nonzero Courant number. The directionally correct upwind finite difference scheme uses the neighbor at $m+1$ and hasIt is contractive exactly for $0\le \nu \le 1$. 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 $J=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has $J^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$. 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 $0<\mu <1$ 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 $(U^{n+1}-U^n)/k=D_{xx}{U}^{n+1}$ with $U^{n+1}$. Discrete integration by parts givesThe 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 $L_h$ is symmetric with nonpositive eigenvalues and the time update is $R(k{L}_h)$, an orthonormal eigenbasis givesThis 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.