Solution (source code)

= Solution

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'=A_hU$, the solution is $e^{tA_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(kA_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 is
$$
\boxed{\|G_{h,k}^n\|_h\leq C_T\quad\text{whenever }0\leq nk\leq T,}
$$
with $C_T$ independent of the refining spatial mesh and allowed time steps. Similarly, semidiscrete stability requires $\|e^{tA_h}\|_h\leq C_T$ for $0\leq t\leq T$. The <norm> should represent the continuum problem, such as the <discrete L2 norm> $\|U\|_{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}$, then
$$
\|G^n\|\leq\|X\|\|X^{-1}\|\max_j|\lambda_j|^n.
$$
Thus <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\leq T$ despite a nontrivial <Jordan block> at one. Even the scalar update $1+k$ is stable on fixed intervals, since $(1+k)^n\leq 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 <eigenvalues>
$$
\lambda(\theta)=-\frac4{h^2}\sin^2(\theta/2).
$$
<Forward 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\leq r\leq1/2$ for all periodic meshes. In two dimensions the corresponding uniform bound is $r\leq1/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$. Since
$$
|G|^2=1-2\nu(1-\nu)(1-\cos\theta),
$$
it is a contraction in the periodic <discrete L2 norm> exactly for $0\leq\nu\leq1$. 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,
$$
A_h=\begin{pmatrix}-1&h^{-1}\\0&-1\end{pmatrix},\qquad
e^{tA_h}=e^{-t}\begin{pmatrix}1&t/h\\0&1\end{pmatrix}.
$$
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=\begin{pmatrix}0&K\\0&0\end{pmatrix}$ and $G_2=\begin{pmatrix}0&0\\K&0\end{pmatrix}$ each have only zero <eigenvalues>, while $G_2G_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.