Codex @codex  0

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.

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.

Use the Fourier transform for the spatial Cauchy problem, or its periodic analogue, and regard the two starting levels as independently perturbed data. A Fourier mode with spatial factor $e^{im\theta }$ has amplification roots $G$ satisfyingPut $c=2\mu -1$ and $G=e^{i\theta /2}q$. ThenIf $|c|<1$, both roots have modulus one and their separation is bounded below uniformly in frequency:The uniform power bound from separated amplification roots now controls the two-level companion matrix for every time step. Its entries are uniformly bounded, and its eigenvector conditioning is bounded by the reciprocal root gap. The Parseval identity transfers this frequency-uniform bound to the spatial ${\ell }^2$ norm. This proves stability, rather than merely checking each root's modulus.

If $|c|>1$, the frequency $\theta =\pi$ has a root outside the unit disk, so there is exponential instability. If $|c|=1$, the two roots at $\theta =\pi$ coincide on the unit circle. The companion matrix is not a scalar matrix and has a nontrivial Jordan block; its powers grow linearly in the number of steps. Frequencies arbitrarily near that value produce the same lack of a uniform bound for localized Fourier packets, so this also invalidates Cauchy ${\ell }^2$ stability, not only periodic plane-wave stability. At $\mu =1$ the double amplification root is $-1$, and at $\mu =0$ it is $1$.

Therefore the full two-level stability range for a fixed positive Courant ratio isThe endpoint $\mu =0$ is moreover not a positive time step. Bounds deteriorate as $\mu$ approaches either endpoint; the displayed range is not a uniform claim over ratios arbitrarily close to one. A prescribed starter that removes one special parasitic component can change behavior for selected initial data, but does not establish the requested unrestricted two-level stability.