Past exam of the mathematics course of the University of Cambridge / 2009 / iii / Paper 72 / 6

The eigenvalue stability analysis of a finite difference method reduces a linear evolution discretization to its modal behavior, but it must retain the norm and the dependence on the mesh. For a homogeneous system obtained by the method of lines, $U^{\prime }=A_{h}U$, the solution is $e^{t{A}_h}U(0)$. A one-step time discretization gives $U^{n+1}=G_{h,k}{U}^n$; for a Runge-Kutta method, where its stage matrices are invertible, $G_{h,k}=R(k{A}_h)$. With a multistep method one instead uses an amplification matrix on the augmented vector of several time levels.

The relevant finite-time stability estimate iswith $C_T$ independent of the refining spatial mesh and allowed time steps. Similarly, semidiscrete stability requires $\Vert e^{t{A}_h}{\Vert }_h\le C_T$ for $0\le t\le T$. The norm should represent the continuum problem, such as the discrete L2 norm $\Vert U{\Vert }_{2,h}^2=h{\sum }_j|U_j{|}^2$ in one dimension. A fixed finite matrix having decaying solutions as $t\to \infty$ is not by itself a mesh-uniform result for the partial differential equation.

If $G=X\Lambda X^{-1}$, thenThus eigenvalues of modulus at most one give a uniform bound when the diagonalizing bases have uniformly bounded condition numbers. For a normal matrix in the chosen inner product, the spectral theorem for normal operators supplies an orthonormal basis, so the condition number is one and the norm of $G^n$ is exactly ${\max }_j|{\lambda }_j{|}^n$. For a self-adjoint or skew-adjoint spatial matrix, applying a scalar stability function preserves this favorable modal structure. This is the setting where eigenvalue calculations provide especially clean, reliable step restrictions.

A precise fixed-matrix statement is also useful: a matrix is power bounded for all nonnegative integers $n$ if and only if all its eigenvalues lie in the closed unit disk and every Jordan block at a unit-modulus eigenvalue has size one. Indeed, powers of an interior Jordan block contain a polynomial in $n$ times a decaying geometric factor and are bounded; a larger block on the unit circle produces unbounded polynomial growth. For a family of matrices, the resulting bounds must still be uniform. Moreover, finite-time stability versus power boundedness distinguishes this all-time condition from the physical-time estimate above. For example $G_k=I+kN$, $N^2=0$, has $G_k^n=I+nkN$, bounded for $nk\le T$ despite a nontrivial Jordan block at one. Even the scalar update $1+k$ is stable on fixed intervals, since $(1+k{)}^n\le e^T$, although it represents a growing equation rather than a contraction.

For a constant-coefficient periodic difference operator, Fourier modes diagonalize the spatial matrix. The discrete Fourier transform is a unitary operator, so Parseval's identity identifies its modal maximum with the operator norm. The eigenvalue stability analysis of a finite difference method then becomes von Neumann stability analysis. For the Dirichlet discrete Laplacian, a discrete sine transform plays the same role. These facts explain its advantages: a large matrix calculation becomes a scalar symbol or a known spectral interval, and the extremal eigenvalue directly determines a safe time step.

As a first example, centered space discretization of the heat equation $u_t=u_{xx}$ has eigenvaluesForward Euler method gives $R(k\lambda )=1-4r{\sin }^2(\theta /2)$ with $r=k/h^2$. Requiring all factors in $[-1,1]$ yields $0\le r\le 1/2$ for all periodic meshes. In two dimensions the corresponding uniform bound is $r\le 1/4$. Because these are normal matrices, these scalar bounds prove contraction rather than merely suggest it.

Backward Euler method has $R(z)=1/(1-z)$, so the centered discretization of the heat equation is contractive for every $k>0$. The Crank-Nicolson method has $R(z)=(1+z/2)/(1-z/2)$ and is likewise contractive for these negative real modes. However, its factor tends to $-1$ for very stiff modes, leaving high-frequency oscillations weakly damped; Backward Euler method tends to zero. The eigenvalues therefore expose the presence of stiff differential equations and damping as well as binary stability. Unconditional stability does not remove accuracy requirements.

For an advection equation example, the upwind finite difference scheme for $u_t+c{u}_x=0$, $c>0$, has factor $G=1-\nu +\nu e^{-i\theta }$, with $\nu =ck/h$. Sinceit is a contraction in the periodic discrete L2 norm exactly for $0\le \nu \le 1$. In contrast, Forward Euler method with centered advection has factor $1-i\nu \sin \theta$. Its modulus exceeds one at nonzero active frequencies, and repeated amplification makes it unstable under fixed nonzero Courant refinement. A much smaller scaling $k=O(h^2)$ can bound its finite-time growth, demonstrating again that exact contraction and the most general stability estimate are different claims.

The major limitation is the possible presence of a nonnormal matrix. For example,Both eigenvalues are negative, yet at $t=1$ the second unit vector is amplified by at least $1/(eh)$. This proves that negative spectra do not imply uniform semidiscrete stability. Ill-conditioned eigenvectors, Jordan normal form and transient amplification can invalidate a spectral-only argument. Direct energy estimates or estimates for the resolvent can supply the missing norm control.

Boundary conditions are another limitation. A periodic symbol proves a periodic or whole-line result, not an arbitrary initial-boundary-value problem. Boundary closure of a difference scheme can introduce growing modes or mesh-dependent amplification. A spatial operator with variable coefficients can still be studied through its assembled matrix, but then whether it is a normal matrix, its condition number and the physical inner product must be checked rather than borrowed from a Fourier stability analysis. For the finite element method the natural mass matrix inner product often turns a generalized eigenvalue problem into a self-adjoint one.

For time-dependent or split updates, separate spectra are insufficient to control products. For instance, $G_1=\left(\begin{array}{cc}0& K\\ 0& 0\end{array}\right)$ and $G_2=\left(\begin{array}{cc}0& 0\\ K& 0\end{array}\right)$ each have only zero eigenvalues, while $G_2{G}_1$ has eigenvalue $K^2$. For $K>1$, alternating the two is unstable. A common contractive norm would control the products, but their individual eigenvalues do not.

Finally, stability and consistency of a numerical method have different roles. The Lax equivalence theorem states that, for a well-posed linear initial-value problem and a consistent finite difference method on compatible grid spaces, stability on every fixed time interval is equivalent to convergence of a numerical method as the admissible meshes refine. Bounded grid representation and the correct initial and boundary setup are part of that framework. An eigenvalue estimate can establish its stability hypothesis in the favorable settings described above; it neither proves consistency of a numerical method nor bypasses its uniformity requirements. The method is therefore a powerful modal tool, provided its spectral calculation is connected to a genuine norm estimate.