= Solution
For an evolutionary <partial differential equation>, <stability of a numerical method> means that errors in the starting data and forcing stay controlled on each fixed interval $0\leq t\leq T$, with constants independent of the mesh. <Fourier stability analysis> is especially effective for a constant-coefficient <finite difference method> on a uniform infinite or periodic grid: translation invariance makes different <Fourier modes> evolve independently.
For a scalar one-step scheme on the integer lattice, insert $U_m^n=\widehat U^n(\theta)e^{im\theta}$, with $-\pi\leq\theta\leq\pi$. The <Fourier symbol> of each shift is $e^{ij\theta}$, and the update reduces to
$$
\widehat U^{n+1}(\theta)=G(\theta)\widehat U^n(\theta).
$$
The <amplification factor> $G$ must be defined for every relevant frequency; for an implicit scheme this includes checking that its denominator is nonzero. The <Parseval identity> turns a bound $\sup_\theta|G(\theta)|^n\leq C_T$ into a discrete <L2 norm> bound. Thus $|G|\leq1$ gives contractivity, while the more general bound $|G|\leq1+Ck$ gives $\|U^n\|_h\leq e^{CT}\|U^0\|_h$ for $nk\leq T$. Conversely, frequencies with amplification uniformly greater than one produce unstable wave packets or periodic <Fourier modes> under refinement.
For a system, $G(\theta)$ is an amplification <matrix>; for a multilevel scheme, a companion <matrix> evolves the vector of time levels. The actual requirement is a uniform bound on matrix powers, not just their <spectral radius>. A nontrivial <Jordan block> at a unit-modulus <eigenvalue> creates polynomial growth in the time index. Even simple <eigenvalues> inside the unit disk can fail to give a uniform bound if the <eigenvector> matrices become ill-conditioned as the mesh changes. Uniformly controlled <diagonalization of a matrix>, or a suitable quadratic energy estimate, resolves this issue. The <power boundedness of a two-level Fourier scheme> illustrates why root multiplicities and conditioning matter.
For the <heat equation>, put $\mu=k/h_x^2$. Centered space with the <Forward Euler method> has
$$
G(\theta)=1-4\mu\sin^2(\theta/2).
$$
The <Von Neumann stability analysis> condition $|G|\leq1$ holds exactly for $0\leq\mu\leq1/2$ on the full frequency interval. The <Backward Euler diffusion scheme> instead has $G=(1+4\mu\sin^2(\theta/2))^{-1}$ and is stable for every $\mu\geq0$. The <Crank-Nicolson diffusion scheme> has
$$
G=\frac{1-2\mu\sin^2(\theta/2)}{1+2\mu\sin^2(\theta/2)}.
$$
It too is unconditionally stable, but poorly resolved high-frequency modes have $G\approx-1$ for large $\mu$, giving oscillatory numerical transients. <A-stability> therefore does not guarantee strong damping; <L-stability> distinguishes the damping of the <Backward Euler method>.
For the <advection equation> $u_t+cu_x=0$, let $\nu=ck/h_x$. Forward time and centered space give
$$
G=1-i\nu\sin\theta,\qquad |G|^2=1+\nu^2\sin^2\theta.
$$
For nonzero fixed $\nu$, repeated steps amplify some modes by a fixed factor greater than one, so the method is unstable under the usual refinement $k\propto h_x$. This calculation also shows the refinement qualification: if $k=O(h_x^2)$, its finite-time growth can be bounded by $\exp(Cc^2Tk/h_x^2)$, though the restrictive scaling defeats the usual hyperbolic time-step choice. For $c>0$, the <upwind finite difference scheme> has
$$
G=1-\nu+\nu e^{-i\theta},\qquad
|G|^2=1-4\nu(1-\nu)\sin^2(\theta/2),
$$
so it is stable for $0\leq\nu\leq1$. The <Lax-Wendroff advection scheme> instead has
$$
G=1-i\nu\sin\theta+\nu^2(\cos\theta-1),\qquad
|G|^2=1-4\nu^2(1-\nu^2)\sin^4(\theta/2),
$$
giving $|\nu|\leq1$. These examples separate the effects of the spatial stencil, temporal approximation and <Courant number>; a higher <order of a numerical method> by itself does not establish <stability>.
A <periodic boundary condition> is ideal for this analysis. The <discrete Fourier transform> diagonalizes the circulant update, with only the discrete frequencies $2\pi j/J$ needed on a grid of $J$ points. A bound on all $\theta$ is a convenient guarantee over every such grid. On the whole line, the <Fourier transform> uses a continuous frequency interval and the <Parseval identity> proves the corresponding square-summable-data estimate.
Homogeneous <Dirichlet boundary conditions> do not admit arbitrary complex exponential <Fourier modes>. For the standard centered second difference on $M$ interior points, a <discrete sine transform> diagonalizes the actual boundary-value <matrix>:
$$
v_j(m)=\sin\frac{jm\pi}{M+1},\qquad
-D_hv_j=\frac4{h_x^2}\sin^2\frac{j\pi}{2(M+1)}v_j.
$$
The same heat amplification formulas then apply at these sine frequencies. For example, the exact forward-Euler contractivity limit on that fixed grid is $\mu\leq[2\cos^2(\pi/(2(M+1)))]^{-1}$; the mesh-independent sufficient limit $1/2$ follows by including all frequencies. Homogeneous <Neumann boundary conditions> can similarly permit a <discrete cosine transform>, provided the endpoint discretization and its weighted <inner product> are chosen consistently. The constant mode is then present, reflecting conservation of the heat equation's spatial mean. For inhomogeneous <boundary conditions>, subtract a suitable lifting of the prescribed boundary data and estimate the induced source using the homogeneous evolution's bound and a <Duhamel principle> estimate. Bounds for the lifting and source must themselves be uniform in the mesh.
General <boundary conditions> require separate <boundary stability of a finite-difference method>. For an <advection equation>, incoming data are prescribed at the inflow boundary, while an outflow closure must respect the outgoing characteristics. An interior periodic <Fourier symbol> cannot detect a growing mode confined near a boundary. A normal-mode test $U_m^n=z^n\kappa^m$ on a half-line seeks modes with $|z|>1$ and $|\kappa|<1$ that satisfy both the interior recurrence and boundary closure. Their existence proves instability, and uniform control requires more than merely excluding isolated growing roots; the <Uniform Kreiss--Lopatinskii condition> addresses boundary resolvent bounds. Alternatively, a direct <energy method> using <summation by parts> can include the boundary terms and establish the required estimate.
Variable coefficients and nonuniform meshes usually destroy exact <Fourier transform> diagonalization. Frozen-coefficient <Fourier stability analysis> is then a useful diagnostic, but it is not automatically a proof for the full variable-coefficient boundary problem. The mesh-uniform <energy method> in question 5 is an example of a direct proof. Also, an <L2 norm> proof is not automatically a mesh-uniform maximum-<norm> proof: finite-dimensional norm-equivalence constants may grow with the number of grid points.
Finally, <stability> measures error propagation; <consistency of a numerical method> measures the defect of inserting the exact solution. For a well-posed linear initial-value problem in the chosen norm, the <Lax equivalence theorem> says that a consistent approximation is convergent exactly when it is stable. Boundary discretization and starting data must be included in that assertion. \b[A Fourier calculation proves convergence only after consistency, well-posedness and the actual boundary treatment have also been checked.]
Back to article page
