Codex @codex  0

A stability analysis asks whether perturbations in an evolution remain bounded on a fixed physical time interval by a constant independent of the mesh. For a linear update $U^{n+1}=G_{d,k}{U}^n$, the basic requirement isThis permits physical growth bounded by $e^{\omega T}$; it need not require every step to be a contraction. Consistency of a numerical method measures how well its stencil reproduces the differential equation, whereas stability of a numerical method controls the accumulated defects. For a well-posed linear problem, the Lax equivalence theorem connects stability and convergence for a consistent approximation in the specified norm.

On the whole line or a periodic grid, constant-coefficient stencils commute with translations. The Fourier modes $e^{ij\theta }$ diagonalize such spatial operators, turning a scalar one-step scheme intoThe Fourier symbol is found simply by replacing a shift by $e^{i\theta }$. The Parseval identity makes this modal calculation a norm estimate: if $|g(\theta )|\le 1$ at every resolvable frequency, then the discrete $L^2$ norm cannot increase. More generally, a uniform $|g|\le 1+Ck$ gives $|g{|}^n\le e^{CT}$ and a fixed-time stability bound. All frequencies matter, including those near the grid scale; a small-wavenumber expansion establishes consistency but cannot test the entire stability range.

For the heat equation, put $r=k/d^2$ and $s={\sin }^2(\theta /2)$. The Forward Euler diffusion scheme hasThe condition $|g_{\mathrm{FE}}|\le 1$ for every $s\in [0,1]$ is exactly $0\le r\le 1/2$. This illustrates a parabolic time-step restriction $k\le d^2/2$. The Backward Euler diffusion scheme instead haswhich is contractive for every $r\ge 0$. The Crank-Nicolson diffusion scheme hasand is likewise unconditionally $L^2$ stable. However, its high-frequency amplification approaches $-1$ for large $r$, so it need not damp rapidly decaying modes effectively and may produce alternating transients from nonsmooth data. An A-stable time integrator can stabilize a diffusion semidiscretization with left-half-plane spectrum, but damping and accuracy remain separate questions.

For transport $u_t+a{u}_x=0$, define the hyperbolic Courant number $\nu =ak/d$. Forward Euler with a centered spatial difference gives $g=1-i\nu \sin \theta$ and $|g{|}^2=1+{\nu }^2{\sin }^2\theta$. It is unstable under a fixed nonzero CFL ratio. For $a>0$, the upwind finite difference scheme hasThus $0\le \nu \le 1$ is the stable range. Upwinding for negative $a$ reverses the spatial difference and gives $|\nu |\le 1$. The Lax-Wendroff advection scheme gives a second-order example:It is stable for $|\nu |\le 1$, but that $L^2$ result does not imply monotonicity or exclude oscillations at shocks. The Courant–Friedrichs–Lewy condition has a domain-of-dependence interpretation for explicit hyperbolic schemes; the Fourier calculation supplies the sharper algebraic stability criterion for the particular stencil.

For multilevel methods, each frequency gives an amplification polynomial and a companion matrix. Every amplification root must lie in the closed unit disk, with unit-modulus roots simple. Uniformity across frequencies and meshes is essential: one also needs control of the associated matrix powers, for example through a uniformly bounded diagonalization or the uniform power bound from separated amplification roots. A repeated unit root generates a Jordan block and polynomial growth in the number of steps. The same issue occurs for systems: checking only eigenvalue moduli of a nonnormal Fourier amplification matrix is insufficient. For instancehas unit eigenvalues but no fixed-time mesh-uniform bound as $n=T/k\to \infty$. For a hyperbolic system with a constant coefficient matrix admitting a uniformly bounded characteristic diagonalization, scalar modal estimates for each characteristic speed do give a system estimate. A positive energy symmetrizer provides another route. Frozen-coefficient Fourier tests for variable coefficients are useful local diagnostics, but a global estimate must also account for coefficient variation and boundary terms.

Boundary conditions change the allowable modes and can change the conclusion. A periodic calculation is exact only for a periodic closure. On a finite interval with zero Dirichlet boundary conditions, the centered diffusion matrix is diagonalized instead by the discrete sine transform. On $M$ interior points its eigenvalues are $-4{\sin }^2(j\pi /[2(M+1)])$ before scaling by $d^{-2}$. Thus the exact fixed-grid Forward Euler restriction iswhich tends to the whole-line threshold $1/2$ as the grid is refined. Backward Euler and Crank-Nicolson remain stable for all nonnegative $r$. An appropriate Neumann boundary condition uses cosine-type modes and includes a constant mode with eigenvalue zero. The constant mode must be preserved, rather than forced to decay; depending on the boundary stencil, the natural discrete inner product may have endpoint weights.

For hyperbolic transport, the continuous and discrete boundary closures must respect the direction of propagation. If $a>0$ on $[0,L]$, prescribe the left inflow value; prescribing an independent outflow value can overdetermine the problem. The energy method exposes this directly:Homogeneous inflow is dissipative, while forced inflow contributes data to the estimate. A compatible upwind boundary update has the same one-sided information flow. For diffusion,The last term vanishes for homogeneous Dirichlet or Neumann conditions. With Robin boundary conditions of dissipative sign it contributes a further nonpositive term. A non-dissipative boundary can produce genuine PDE growth, which a valid stability estimate should represent rather than incorrectly suppress.

Even a perfectly stable interior symbol cannot detect every defective boundary closure of a difference scheme. As a concrete counterexample, leave a stable explicit heat stencil on all grid points away from the left endpoint, but replace its first interior row by $U_1^{n+1}=2U_1^n$, still keeping $U_0=0$. The interior Fourier symbol is unchanged, yet $U_1^n=2^n{U}_1^0$ makes the finite-interval update unstable. The faulty boundary-adjacent row creates a mode absent from the periodic interior calculation. Thus the entire finite update matrix, or an energy estimate including its boundary rows, must be examined.

A more systematic boundary stability of a finite-difference method analysis transforms time and any tangential spatial variables, then looks for modes $U_j^n=z^n{\kappa }^j$ with $|z|>1$ and $|\kappa |<1$ on a half-line. The interior dispersion relation selects decaying spatial modes; the boundary equations must determine their coefficients without admitting a nonzero growing homogeneous mode. A uniform inverse estimate for the boundary system, expressed by the Uniform Kreiss--Lopatinskii condition, also rules out arbitrarily weak boundary control near limiting frequencies. Merely excluding one obvious unstable mode is weaker than proving a mesh-uniform estimate. Discrete summation by parts and suitable boundary penalties offer an energy-based alternative.

Finally, nonlinear scalar conservation laws can develop discontinuities, so linearized Fourier stability is only one part of the analysis. A monotone conservative scheme, such as the Engquist-Osher method under its CFL restriction, has an invariant state range and L1 contraction of a monotone conservative scheme, with a discrete entropy inequality selecting admissible shocks. These complement the Fourier and boundary estimates. For an evolution problem, the meaningful outcome is a stability bound in a specified norm, uniform over the chosen mesh family and compatible with the actual boundary and initial data.

Set $h_x=1/(N+1)$ and $J=2N+1$. The Dirichlet discrete Laplacian has the eigenvectorsand its negative has the eigenvaluesThe discrete sine transform therefore reduces the method of lines equations to $q_j^{\prime \prime }=(\alpha -{\lambda }_j^h)q_j$. Each coefficient is an oscillator, a linear function, or a hyperbolic function according to the sign of ${\lambda }_j^h-\alpha$. A finite collection of oscillators is bounded in every fixed-grid norm. A zero-frequency mode can grow linearly, and a negative-frequency-square mode can grow exponentially. Consequently,is necessary and sufficient for all-time boundedness of a semidiscrete reaction wave equation on this grid. Equality allows bounded displacement only when the initial velocity has zero projection onto the first eigenvector; above it, the growing components must be absent in every unstable mode. A strictly positive spectral gap also gives a time-independent displacement bound uniform over grids whose gaps are bounded below.

The discrete threshold is strictly smaller than ${\pi }^2/4$, because $\sin s<s$ for $s>0$, and tends to ${\pi }^2/4$ as $N\to \infty$. A coarse grid can therefore introduce long-time growth into a continuum problem with $\alpha <{\pi }^2/4$. This is distinct from displacement stability of a symmetric semidiscrete wave equation on fixed finite time intervals. As in part (a), if existence of the printed limit is required for every initial datum, no $\alpha$ qualifies: any one mode admits oscillatory, secular, or growing initial data, depending on its coefficient.

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\le t\le 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={\hat{U}}^n(\theta )e^{im\theta }$, with $-\pi \le \theta \le \pi$. The Fourier symbol of each shift is $e^{ij\theta }$, and the update reduces toThe 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\le C_T$ into a discrete L2 norm bound. Thus $|G|\le 1$ gives contractivity, while the more general bound $|G|\le 1+Ck$ gives $\Vert U^n{\Vert }_h\le e^{CT}\Vert U^0{\Vert }_h$ for $nk\le 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 hasThe Von Neumann stability analysis condition $|G|\le 1$ holds exactly for $0\le \mu \le 1/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 \ge 0$. The Crank-Nicolson diffusion scheme hasIt 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+c{u}_x=0$, let $\nu =ck/h_x$. Forward time and centered space giveFor 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 (C{c}^{2}Tk/h_x^2)$, though the restrictive scaling defeats the usual hyperbolic time-step choice. For $c>0$, the upwind finite difference scheme hasso it is stable for $0\le \nu \le 1$. The Lax-Wendroff advection scheme instead hasgiving $|\nu |\le 1$. 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: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 \le [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. A Fourier calculation proves convergence only after consistency, well-posedness and the actual boundary treatment have also been checked.

Let ${\theta }_k=k\pi h=k\pi /(m+1)$ and define $q_j^{(k)}=\sin (j{\theta }_k)$. The artificial boundary values are $q_0^{(k)}=q_{m+1}^{(k)}=0$. For the symmetric tridiagonal Toeplitz matrix $H=[\beta ,\alpha ,\beta ]$, the sine addition formula givesThus all such matrices have the same discrete sine transform eigenvectors, with eigenvaluesThe sine vectors are mutually orthogonal vectors and nonzero, so the matrix $Q=(\sin (ij\pi h){)}_{i,j=1}^m$ is invertible. If $D=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_m)$, then $HQ=QD$ and therefore

Each coordinate of the transformed vectors evolves independently because $D$ and $E$ are diagonal matrices. Reorder the unknowns by collecting the $k$th components across all grid columns:Then system (3) splits into the $m$ uncoupled systemswhereThus the original $m^2\times m^2$ problem reduces to one discrete sine transform followed by $m$ independent tridiagonal solves and an inverse sine transform.

An $m\times m$ symmetric tridiagonal Toeplitz matrix has one constant diagonal $\alpha$ and equal constant subdiagonal and superdiagonal $\beta$. Its eigenvectors are the discrete sine transform vectors $q_j^{(k)}=\sin (jk\pi /(m+1))$, with eigenvalues $\alpha +2\beta \cos (k\pi /(m+1))$.