Codex @codex  0

The domain of dependence of a solution value consists of the initial and boundary locations whose data can influence it. For $u_t+a{u}_x=0$ on the whole line, $u(x,t)=u(x-at,0)$, so its initial-data domain of dependence is the foot of the characteristic. For a finite-difference scheme, repeated stencil dependencies determine a numerical domain of dependence. The Courant–Friedrichs–Lewy condition requires the numerical domain to cover the physical one for convergence; it is necessary under its usual hyperbolic hypotheses, but not in general sufficient for stability.

The Engquist-Osher numerical flux is $F(a,b)=f(0)+{\int }_0^a\max (f\prime (s),0)ds+{\int }_0^b\min (f\prime (s),0)ds$. It sends positive-speed contributions from the left state and negative-speed contributions from the right state. Its consistency $F(u,u)=f(u)$ and opposite monotonicities in the two arguments yield a monotone conservative scheme under the appropriate Courant–Friedrichs–Lewy condition. For the Inviscid Burgers equation, it is $\frac{1}{2}\max (a,0{)}^2+\frac{1}{2}\min (b,0{)}^2$.

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.

The Engquist-Osher method is an upwind finite volume method. For cell values $U_j^n$, mesh width $h$ and time step $k$, put $r=k/h$ and use the conservative updatewhere the numerical flux isConsistency of a numerical method follows from $F(u,u)=f(u)$. Positive characteristic speeds use the left state, and negative speeds use the right state. The flux difference telescopes when summed over cells, giving a discrete conservation law. The method is first order in space and time for smooth solutions; its conservation law form is also appropriate when shock waves develop.

Let $s_{\ast }$ be the unique sonic point of a scalar flux, $f^{\prime }(s_{\ast })=0$. Convexity implies $f^{\prime }<0$ to its left and $f^{\prime }>0$ to its right. Put $f_{\ast }=f(s_{\ast })$ and splitThen $f=f_{\ast }+f^{+}+f^{-}$, $(f^{+}{)}^{\prime }\ge 0$, $(f^{-}{)}^{\prime }\le 0$, andThis is the sonic-point splitting of a convex Engquist-Osher flux. If both states are above $s_{\ast }$, $F=f(a)$; if both are below it, $F=f(b)$. Across a transsonic rarefaction wave, $a\le s_{\ast }\le b$, the flux is the minimum $f_{\ast }$. In the opposite crossing it is $f(a)+f(b)-f_{\ast }$, which in general differs from the exact Godunov numerical flux. This distinction does not affect the monotonicity proof.

Here stability requires the usual Courant–Friedrichs–Lewy condition. Let the initial values lie in $[m,M]$, and chooseA differentiable convex function has a continuous derivative on this interval, so the maximum is finite. Write the three-input update as $H(a,b,c)$. Differentiating givesThus the update $T$ preserves componentwise order. Constants are fixed, since $H(c,c,c)=c$, so comparison with the constants $m,M$ proves $m\le U_j^n\le M$ for every time level. In particular the same CFL condition remains applicable. This is a mesh-independent maximum-principle stability bound.

A stronger stability statement is contraction in discrete $L^1$. Consider two solutions $U,V$ in the same invariant interval, on a periodic mesh or on the infinite lattice with summable differences. Let $Q=U\vee V$ be their componentwise maximum. Monotonicity gives $TQ\ge TU,TV$, henceThe discrete conservation law for the difference gives the equality: on a finite periodic mesh the flux sum cancels exactly, and on the infinite lattice it follows by truncation and the fact that the flux is a Lipschitz function on $[m,M]$. Interchanging $U,V$ and adding provesThis is L1 contraction of a monotone conservative scheme. Since $T$ commutes with the cell shift, applying this contraction to $U$ and its shift gives ${\sum }_j|U_{j+1}^{n+1}-U_j^{n+1}|\le {\sum }_j|U_{j+1}^n-U_j^n|$. The scheme is therefore a total variation diminishing scheme as well.

The same order argument proves a discrete entropy inequality, explaining why this stable discretization selects entropy solutions rather than nonphysical expansion shock waves. For a constant $c\in [m,M]$, setBecause $T(U\vee c)\ge TU,c$ and $T(U\wedge c)\le TU,c$, their componentwise difference is at least $|TU-c|$. Subtracting their updates givesFor constants outside $[m,M]$, the entropy inequality is an equality by the invariant bound and the conservative update. At equal states $Q_c(u,u)=\mathrm{sgn}(u-c)[f(u)-f(c)]$, the Kruzhkov entropy flux. The flux is a Lipschitz function on the invariant interval: writing $L=\max |f^{\prime }|$, the conservative update gives $h{\sum }_j|U_j^{n+1}-U_j^n|\le 2kL{\sum }_j|U_{j+1}^n-U_j^n|$. This time-translation bound, the invariant bound and the total variation estimate give local space-time compactness for bounded-variation initial conditions as the mesh is refined; the conservative update passes to the weak conservation law, and the displayed inequality passes to its Kruzhkov entropy inequality. This identifies the limit as the entropy solution.

Finally the time-step qualification is essential. For the convex Inviscid Burgers equation flux $f(u)=u^2/2$, whose unique sonic point of a scalar flux is zero, linearize about a positive constant $U$. The perturbation update becomes $v_j^{n+1}=(1-rU)v_j^n+rU{v}_{j-1}^n$. At Fourier frequency $\pi$ its multiplier is $1-2rU$, with modulus greater than one if $rU>1$. Wave packets near this phase give linear instability even for summable Cauchy perturbations. Thus the assumptions on $f$ imply the stability results above under the CFL condition, not for arbitrary $k/h$. This is CFL necessity for explicit Engquist-Osher stability.

Fourier stability analysis exploits translation invariance of a linear constant-coefficient finite difference method. On an infinite uniform grid of spacing $h$, a spatial shift acts on $e^{ij\theta }$ by multiplication by $e^{i\theta }$. Hence a scalar two-level scheme with new- and old-time symbols $L(\theta ),N(\theta )$ propagates each Fourier mode by $G(\theta )=N(\theta )/L(\theta )$, provided the new-time operator is invertible. An explicit stencil $U_j^{n+1}={\sum }_{\ell }{a}_{\ell }{U}_{j+\ell }^n$ simply has $G(\theta )={\sum }_{\ell }{a}_{\ell }{e}^{i\ell \theta }$. This converts a spatial evolution into scalar algebra at each frequency.

The Plancherel theorem gives the full norm estimate, rather than just a heuristic about waves. For the discrete Fourier transform $\hat{U}(\theta )={\sum }_j{U}_j{e}^{-ij\theta }$, interpreted in the $L^2$ sense if necessary,Thus a scalar mode bound gives an exact discrete L2 norm bound. Stability on a fixed physical interval requires constants independent of the refining spatial and time meshes:The contraction condition $|G|\le 1$ is sufficient and customary for equations whose exact solutions do not grow. More generally $|G|\le 1+Ck$ suffices, since $(1+Ck{)}^n\le e^{CT}$. At a fixed Courant number, when $G$ is independent of $k$, a mode with modulus strictly greater than one violates the estimate as $n\sim T/k\to \infty$. One must keep this scaling: physical growth such as $e^{\kappa t}$ does not by itself mean finite-time numerical instability. This is finite-time stability versus power boundedness.

For systems, the symbol is an amplification matrix. The requirement is a uniform bound on all its powers, not merely its eigenvalues. A unit-modulus Jordan block $\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has powers $\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$ and is unstable under refinement. A uniformly conditioned eigenbasis or an energy estimate can justify a modal conclusion; eigenvalues alone cannot control a general nonnormal matrix. Multilevel schemes likewise produce an amplification polynomial of a multilevel finite difference scheme, with a uniform root condition and control of starting errors. Repeated unit-modulus roots can give polynomial growth, and roots approaching each other require attention to mesh-uniform constants.

Several elementary examples show how the calculation works. For $u_t+c{u}_x=0$ with $c>0$, forward Euler with centered spatial differences gives $G=1-i\mu \sin \theta$, where $\mu =ck/h$. Since $|G{|}^2=1+{\mu }^2{\sin }^2\theta$, it is unstable at every fixed nonzero Courant number. The upwind finite difference scheme instead gives $G=1-\mu +\mu e^{-i\theta }$ andso its exact contraction range is $0\le \mu \le 1$. The Lax-Wendroff advection scheme has $G=1-i\mu \sin \theta +{\mu }^2(\cos \theta -1)$; expansion gives $|G{|}^2=1-4{\mu }^2(1-{\mu }^2){\sin }^4(\theta /2)$, hence $|\mu |\le 1$. The Courant–Friedrichs–Lewy condition relates these restrictions to the numerical domain of dependence, but is a necessary geometric check rather than a complete stability proof. Implicit schemes can have a global spatial dependence and different Courant ranges, as Question 3 illustrates.

For the heat equation $u_t=\nu u_{xx}$, put $r=\nu k/h^2$. Central spatial differences with forward Euler give $G=1-4r{\sin }^2(\theta /2)$, so $0\le r\le 1/2$. With backward Euler the factor is $[1+4r{\sin }^2(\theta /2){]}^{-1}$, a contraction for every $r\ge 0$. For the Crank-Nicolson diffusion scheme it is $[1-2r{\sin }^2(\theta /2)]/[1+2r{\sin }^2(\theta /2)]$, again a contraction for all $r\ge 0$. At large $r$, however, high-frequency factors approach minus one rather than zero: unconditional stability does not imply strong damping or accurate resolution of every mode.

The connection with convergence is the Lax equivalence theorem for a consistent approximation to a well-posed linear evolution problem. Its forward implication follows directly from the error recursion. If $e^{n+1}=G_h{e}^n+k{\tau }^n$ and $\Vert G_h^j\Vert \le C_T$ for $jk\le T$, thenA convergent initialization and a vanishing normalized local truncation error therefore give convergence. Smooth-data consistency can be extended by density when the approximation and initialization are uniformly bounded. For rough $L^2$ data, cell averages or another L2-compatible initialization of grid data are preferable to undefined point samples. A high formal order does not rescue an unstable amplification operator.

Boundary conditions determine whether Fourier diagonalization is exact. On a periodic grid the update is a circulant matrix, and the finite discrete Fourier transform diagonalizes it at the frequencies $2\pi \ell /N$. Here an interior symbol controls the entire scheme, including the wrapping boundary. On the whole line, the continuous frequency transform serves the same purpose and there is no boundary closure to check.

For a finite interval with homogeneous Dirichlet boundary conditions, centered diffusion is a symmetric tridiagonal problem diagonalized by discrete sine modes. One tests the same diffusion symbol at the admitted sine frequencies, not at arbitrary periodic modes. For homogeneous Neumann boundary conditions, a suitable reflected closure uses discrete cosine modes; the constant mode is present and is neutral for pure diffusion. With the standard reflected endpoint stencil, the natural norm uses half endpoint weights, in which the matrix is self-adjoint. Consequently stable diffusion need not make every solution tend to zero under Neumann conditions. These special orthogonal bases supply an exact substitute for Fourier modes when the operator and closure have the required symmetry.

For transport, prescribe only the incoming physical boundary data. With $u_t+c{u}_x=0$, $c>0$, the inflow is the left endpoint; for $u_t=u_x$ it is the right endpoint. Imposing unrelated data at both ends generally overdetermines the first-order PDE. Outflow rows or ghost points must be chosen compatibly with the interior method. A general boundary closure of a difference scheme destroys the circulant structure, so whole-line von Neumann stability analysis is necessary in many settings but not sufficient for the full initial-boundary problem. A defective boundary row that doubles one boundary unknown each step can create growth even beside a stable interior stencil.

A more refined boundary calculation seeks half-line modes $U_j^n=w^n{\zeta }^j$ with $|\zeta |<1$; these are localized near the boundary. Insert them into both the interior recurrence and boundary equations. A compatible mode with $|w|>1$ proves boundary instability. Absence of one isolated growing mode is not the whole theorem: weakly controlled limiting modes can still destroy a uniform estimate. The Uniform Kreiss--Lopatinskii condition requires a uniformly bounded inverse boundary system on the incoming/decaying mode space. Direct energy estimates are often a clearer alternative to solving that transformed boundary problem.

For example, suppose a first-derivative matrix obeys the discrete summation by parts identity $HD+D^{T}H=e_N{e}_N^T-e_0{e}_0^T$, with positive norm matrix $H$. The semidiscrete advection equation $U_t=-cDU$ then givesWith homogeneous inflow $U_0=0$ and $c>0$, only a nonpositive outflow term remains; prescribed inflow contributes a controlled forcing term. This is the inflow advection energy estimate from summation by parts, and it tests the boundary as part of the operator. For a dissipative semidiscrete matrix $A$, backward Euler also preserves the estimate, because its update givesThis distinguishes spatial energy stability from stability of an arbitrary chosen time integrator.

Finally, nonhomogeneous boundaries can be homogenized by subtracting a suitable lifting, leaving a forced equation whose solution needs a bound in terms of both initial and boundary data. Variable coefficients and nonuniform grids remove translation invariance, so freezing coefficients gives only local diagnostic information. A full energy method, a mesh-uniform operator estimate, or a dedicated boundary-symbol analysis is then needed. Fourier symbols give complete stability tests when an orthogonal mode decomposition diagonalizes the full operator; otherwise the boundaries, matrix powers and norm estimates must be checked as well.

Choose a spatial interval with periodic boundary conditions and a uniform grid, so no separate inflow boundary closure of a difference scheme is required. The discrete Fourier transform is a unitary operator in the discrete L2 norm. Every circulant stencil is diagonal in that transform, and a one-step Fourier multiplier is a contraction exactly when its modulus is at most one for every grid frequency. For a mesh-uniform statement, use all $\theta \in [-\pi ,\pi ]$, the limiting frequency range.

A Fourier mode $e^{im\theta }$ giveswith the coefficients from (a). Squaring moduli and putting $v=\cos \theta$ yieldsThis has the required nonnegative sign for every frequency exactly whenOne must also check that the implicit equation is solvable, rather than divide by a vanishing symbol. HereFor $\mu \ne 1/2$, a zero would have to occur at zero or $\pi$; the latter occurs only at $(1\pm \sqrt{3})/2$, both outside the displayed stable intervals. At $\mu =1/2$, $D=3/4+(1/4)\cos \theta \ge 1/2$. Thus there is no denominator zero in the stable range. For fixed such $\mu$, continuity on the compact frequency interval gives a uniformly bounded inverse, and $|G|\le 1$ proves $\Vert U^n{\Vert }_{2,h}\le \Vert U^0{\Vert }_{2,h}$ for all steps.

Outside these intervals the displayed difference is negative at every nonzero frequency where $D\ne 0$, so $|G|>1$. On periodic refining meshes there are frequencies approaching any selected such value. Their amplification is bounded away from one and is repeated $O(h^{-1})$ times over a fixed physical interval when $k=\mu h$, violating a mesh-uniform stability bound. At the two singular parameters the implicit operator itself fails on meshes containing the Nyquist mode. These observations prove necessity as well as sufficiency.

For the physical convention $k>0$, $h>0$, the result is $0<\mu \le 1$ or $\mu \ge 2$, with zero admitted as the identity limit. The negative intervals describe the algebraic contraction result for signed Courant numbers. In particular, the additional large positive interval must not be discarded by imposing the usual explicit-scheme Courant–Friedrichs–Lewy condition: this scheme is implicit. Stability at large fixed Courant number is not a guarantee of a small error constant.

Eigenvalue analysis proves mesh-uniform stability when it is accompanied by control of eigenvectors or an energy norm; the signs of eigenvalues alone are insufficient. For an evolution equation, the method of lines first produces $u_h^{\prime }=L_h{u}_h$. A time integrator then gives an amplification matrix $S_{h,\delta }$, for example $S_{h,\delta }=R(\delta L_h)$ for a Runge-Kutta method. The useful stability requirement on a fixed physical time interval iswhere $C_T$ is independent of the refining spatial and temporal meshes. The semidiscrete analogue is $\Vert e^{t{L}_h}{\Vert }_h\le C_T$ for $0\le t\le T$. These are bounds on propagation of initial errors and residuals, not merely statements that each finite matrix has a bounded solution over its own finite time interval.

If $L_h$ is a normal matrix in the chosen discrete inner product, it has a unitary eigenbasis. ThenFor dissipative problems it is sufficient that every $\delta {\lambda }_j$ lie in the integrator's linear stability domain. A diagonalizable family also suffices if its eigenvector matrices have uniformly bounded condition numbers: the same maximum is multiplied by $\Vert V_h\Vert \Vert V_h^{-1}\Vert$. Without that uniform bound the argument does not establish mesh stability. For finite-time estimates, an amplification modulus $1+O(\delta )$ can also be acceptable; strict contraction at every step is a stronger property.

For a constant-coefficient stencil on the full lattice or a periodic grid, von Neumann stability analysis provides the eigenvectors explicitly. A Fourier mode diagonalizes each translation-invariant spatial difference, giving a scalar symbol ${\lambda }_h(\theta )$ or a small matrix symbol for a system. The discrete Fourier transform and Parseval identity turn a uniform modal bound into an actual discrete L2 norm bound. The two-dimensional stencil in question 3 has $\mathrm{Re}{\lambda }_h=-(1-\cos \eta {)}^2/h\le 0$, so its continuous-time amplification is contractive; a later time integrator must still be checked. With boundaries one must analyze the boundary rows as well: a whole-line Fourier symbol cannot detect an unstable boundary closure.

For the heat equation $u_t=\kappa u_{xx}$ on a unit interval with zero Dirichlet data, the centred second difference has a sine eigenbasis and eigenvaluesExplicit Euler method multiplies mode $j$ by $1+\delta {\lambda }_j$. Requiring all factors to have modulus at most one yields $\kappa \delta /h^2\le 1/[2{\sin }^2((N-1)\pi /(2N))]$ on a given grid, and the sharp grid-independent sufficient limit isThese parabolic time-step restrictions express the growth of the largest-magnitude diffusion eigenvalue like $h^{-2}$. Backward Euler method instead has factor $(1-\delta {\lambda }_j{)}^{-1}$ and the Crank-Nicolson method factor $(1+\delta {\lambda }_j/2)/(1-\delta {\lambda }_j/2)$, both bounded by one for every $\delta \ge 0$. Their unconditional stability follows from A-stability, with much stronger damping of the shortest scales for backward Euler.

For $u_t=-a{u}_x$, $a>0$, explicit time stepping with a backward spatial difference givesThus the upwind finite difference scheme is contractive exactly for $0\le \nu \le 1$, a Courant–Friedrichs–Lewy condition matching the direction and speed of information propagation. A centred spatial difference has purely imaginary eigenvalues, and forward Euler gives $|S{|}^2=1+{\nu }^2{\sin }^2\theta >1$. Under a fixed positive hyperbolic Courant number, high-frequency data grow by a fixed factor in each of $O(1/h)$ steps, proving instability as the mesh is refined. This is not a claim that no exceptionally small coupled step can ever give a finite-time bound: if $\delta =O(h^2)$, the excess per step is only $O(\delta )$, but the usual hyperbolic scaling is unstable. Higher-order explicit integrators can possess nontrivial imaginary-axis stability intervals and are then suitable for skew-adjoint advection or the energy-form first-order wave equation.

Variable coefficients usually prevent Fourier diagonalization, but may preserve a self-adjoint spatial operator. The positive-face diffusion discretization in question 5 is an example. Its discrete energy identity makes $A+B$ symmetric negative definite, so a unitary eigenbasis still exists even without a closed eigenvalue formula. More generally, a dissipative operator $L_h$ in a discrete energy inner product obeyswhich directly proves a uniform semigroup estimate. The dissipative Cayley-transform contraction also follows without normality: for $\alpha \ge 0$,Taking $v=(I-\alpha L_h{)}^{-1}u$ proves contraction of the Crank–Nicolson factor. The contractivity of split Crank-Nicolson diffusion then follows from multiplying two contraction bounds; simultaneous eigenvectors of the two directional matrices are unnecessary.

The negative spectra do not imply uniform semidiscrete stability example shows what can go wrong. LetBoth eigenvalues are $-1$ for every mesh, but applying the propagator to the second unit vector at $t=1$ gives norm at least $1/(eh)$. The evolution is therefore not mesh-uniformly stable. For a discrete matrix, a unit-modulus eigenvalue with a nontrivial Jordan block gives powers growing polynomially in the step number. Even strictly interior eigenvalues can allow large transient growth when eigenvectors are badly conditioned. A pseudospectrum or numerical range of an operator, and especially a direct energy estimate, can reveal behavior missed by eigenvalues.

A resolvent estimate provides another way to see the needed uniformity. If $\Vert S_h^n\Vert \le C$ for every $n$, then for $|z|>1$,This shows directly that a uniform power estimate imposes more than the location of the spectrum. Conversely, resolvent-to-power estimates in varying dimension must themselves be uniform; a bound with constants depending on the number of grid points is not enough.

For a linear multistep method, each spatial eigenvalue gives an amplification polynomial, and all roots must satisfy the root condition for a multistep method, including simplicity of any unit-modulus root. Merely following the root approximating $e^{\delta \lambda }$ misses parasitic modes. The temporal and spatial analyses must therefore be combined with the appropriate zero-stability and uniform modal bounds.

Finally stability connects to convergence. If the exact grid samples have normalized residual ${\tau }^n$ and the error obeys $e^{n+1}=S_h{e}^n+\delta {\tau }^n$, iteration givesA residual tending uniformly to zero then gives convergence on $[0,T]$. For a well-posed linear initial-value problem and a consistent linear finite-difference approximation, this is the forward implication of the Lax equivalence theorem. The practical conclusion is to analyze the full operator including boundaries, select a time method whose stability region contains the relevant scaled spectrum, and verify that eigenvectors or an energy estimate make the resulting bounds uniform as the grid is refined.

Fourier stability analysis applies most directly to linear, constant-coefficient evolution schemes on uniform infinite or periodic grids. Translation invariance makes spatial shifts diagonal in discrete Fourier modes. If a one-step scheme has $U^{n+1}={\sum }_j{A}_j{U}_{m+j}^n$, substituting $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$ gives the Fourier amplification symbolThe discrete Parseval identity identifies the discrete L2 norm with the squared norm of Fourier coefficients. Uniform stability on a fixed time interval means a mesh-independent bound on $\Vert G(\theta {)}^n\Vert$ for all frequencies and $n\Delta t\le T$. In the scalar case $|G|\le 1$ yields contraction. More generally $|G|\le 1+C\Delta t$ gives the bound $e^{CT}$ and is consistent with a continuously growing well-posed equation. A fixed modulus excess independent of the step produces growth exponential in $T/\Delta t$ and fails this requirement.

For systems, eigenvalue moduli alone are insufficient. The uniform power bound for matrix Fourier symbols also controls eigenvector conditioning and Jordan growth. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has only unit eigenvalues but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$, unbounded as $\Delta t\to 0$ at fixed $T$. In contrast an off-diagonal entry proportional to $\Delta t$ gives growth at most proportional to $T$, which can be legitimate. A uniformly conditioned diagonalization is a sufficient reduction to scalar modal bounds; nonnormal matrices require direct power or energy estimates. Multilevel schemes similarly require the root condition for their frequency-dependent amplification polynomial and a uniform bound on the associated companion evolution, including limiting repeated-root cases.

The simplest parabolic example is the Forward Euler method with centered spatial diffusion. For $\mu =\Delta t/h^2$,is the sharp mesh-uniform contraction range. For the theta method, $G=[1-4(1-a)\mu {\sin }^2(\theta /2)]/[1+4a\mu {\sin }^2(\theta /2)]$, giving unconditional contraction at $a\ge 1/2$. Implicit stability permits larger steps, but accuracy still constrains them, and the Crank-Nicolson method does not strongly damp high-frequency rough data.

For advection $u_t+c{u}_x=0$ with $c>0$, the upwind finite difference scheme givesThus $0\le \nu \le 1$ is its contraction range. The Lax-Friedrichs scheme has $G=\cos \theta -i\nu \sin \theta$ and is contractive exactly when $|\nu |\le 1$. These express a Courant–Friedrichs–Lewy condition: the numerical domain of dependence must encompass the physical one for convergence. That necessary condition is not sufficient by itself. Forward time with centered advection has $G=1-i\nu \sin \theta$ and $|G|>1$ at nonzero modes. At a fixed nonzero Courant number it is unstable even if $|\nu |\le 1$. Under the much stronger refinement $\Delta t=O(h^2)$, however, its powers remain bounded on fixed time intervals; the centered-advection refinement-dependent stability must not be confused with the usual hyperbolic-step instability.

Symbols also quantify dispersion and dissipation. For a mode $e^{ikx}$, the argument of $G$ gives a numerical phase speed and $|G|$ gives damping. Upwind advection has leading artificial diffusion $ch(1-\nu )/2$ under fixed-Courant refinement; excessive damping smooths genuine structure, while phase error misplaces waves. Stability controls boundedness, not either error by itself. Smooth low-frequency Taylor expansions establish consistency; high-frequency modal bounds establish stability, so analyzing only fixed physical wavenumbers can miss grid-scale instabilities.

For a well-posed linear initial-value problem, the Lax equivalence theorem makes consistency plus a uniform stability estimate equivalent to convergence in the chosen norm. Boundaries and variable coefficients need extra care: Fourier diagonalization no longer handles the full update, and an interior or frozen-coefficient symbol does not prove boundary or global stability. Boundary energy estimates, a suitable matrix analysis or a separate boundary-stability argument must supplement it. Nonlinear schemes also require controls beyond this linear symbol test. These limitations delimit what the method proves rather than reducing its usefulness: in its proper setting, it turns an arbitrarily large grid evolution into exact and interpretable modal algebra.

For fixed finite matrices, the matrix exponential and a Taylor expansion showThis agrees through degree two with $e^{h(A+B)}$, even when $A$ and $B$ do not commute. On a fixed time interval, $\Vert S(h)\Vert \le e^{h(\Vert A\Vert +\Vert B\Vert )}$ and the telescoping identityturns the $O(h^3)$ one-step defect into $O(N{h}^3)=O(h^2)$ global error for $Nh\le T$. Hence Strang splitting is second order in general; for commuting matrices it is exact.

For the spatial discretization, use $d=2/(M+1)$, $x_m=-1+md$ and the interior vector $U=(u_1,\dots ,u_M{)}^T$, with $u_0=u_{M+1}=0$. The central finite differences giveHere $\mathrm{tridiag}$ lists lower diagonal, main diagonal and upper diagonal, in that order. Thus $A^T=A$ and $B^T=-B$. Discrete summation by parts givesThe matrix exponential $e^{hA/2}$ is consequently a contraction in the Euclidean norm, while $e^{hB}$ is an orthogonal matrix for real $\alpha$. ThereforeRepeated split steps are contractive in the mesh-weighted norm $d\sum |u_m{|}^2$, uniformly in the spatial mesh and without a Courant–Friedrichs–Lewy condition. The split semidiscretization is unconditionally stable. No commutativity or shared eigenvectors are needed for this Strang splitting contraction for symmetric diffusion and skew advection. Its time-order proof is for fixed spatial matrices; as $d\to 0$, commutators can grow, so this stability result alone is not a uniform-in-mesh second-order error estimate.

The Engquist-Osher method is a conservative finite volume method built by separating positive and negative characteristic speeds. Take $f\in C^1$, a uniform cell width $d$, a time step $k$, and cell averages $U_j^n$. Write $\lambda =k/d$ and define the numerical fluxThe Engquist-Osher flux uses the left state for right-going transport and the right state for left-going transport. Its conservative update isIt is consistent because $F(u,u)=f(u)$. The integrals have their ordinary oriented meaning even for negative states. For the Inviscid Burgers equation, $f(u)=u^2/2$ gives the useful exampleIn a region where all characteristic speeds have one sign, the flux reduces to the corresponding upwind flux. In smooth regions, this basic piecewise-constant, forward-time version is first order in space and time; its main virtue is robust nonlinear stability across shocks.

Here is a precise stability proof. Assume the data lie in a bounded interval $J=[m,M]$ and chooseWork on a periodic grid, or on the infinite grid with summable differences; on a finite interval the inflow boundary data and boundary fluxes must be treated monotonically as well. Define the three-point update mapWriting $f_{+}^{\prime }=\max (f^{\prime },0)$ and $f_{-}^{\prime }=\min (f^{\prime },0)$, its three partial derivatives areThus it is a monotone conservative scheme. Since $H(s,s,s)=s$, comparison with the constant states provesThe interval is invariant, so the same Courant–Friedrichs–Lewy condition remains valid at every step. This is an $L^{\infty }$ bound relative to constant states, not a claim that arbitrary pairs of solutions are contractive in $L^{\infty }$.

For a stronger perturbation estimate, let $\mathcal{T}$ denote the global update and use componentwise maxima. Monotonicity gives $\mathcal{T}(U\vee V)\ge \mathcal{T}U,\mathcal{T}V$. ThereforeSum over the grid. The shared interface flux cancels telescopically, so conservation givesInterchanging $U,V$ and adding yields the L1 contraction of a monotone conservative scheme:This is a mesh-independent nonlinear stability estimate. On an infinite grid the telescoping argument is justified by summability of the differences and the bounded characteristic-speed range, or by truncation followed by a limit. For unequal prescribed boundary data, additional boundary flux terms enter the estimate.

The method is also a total variation diminishing scheme. Let $\mathcal{S}$ shift a grid sequence by one cell. Translation invariance gives $\mathcal{T}\mathcal{S}=\mathcal{S}\mathcal{T}$. Applying the preceding $L^1$ contraction to $U$ and $\mathcal{S}U$ showsHence initially finite discrete total variation cannot grow.

Finally, monotonicity supplies a discrete entropy inequality, explaining why the stability is compatible with the physically admissible entropy solution. For a constant $c\in J$, defineComparison of $H$ with the clipped triples above and below $c$ givesIndeed $H(U\vee c)\ge \max (H(U),c)$ and $H(U\wedge c)\le \min (H(U),c)$, and subtracting these two update formulas yields exactly the right side. Also $Q(u,u;c)=\mathrm{sgn}(u-c)(f(u)-f(c))$, the Kruzhkov entropy flux. Under the displayed CFL condition the method therefore has an invariant state range, $L^1$ contraction, decreasing discrete total variation and entropy dissipation. These estimates remain meaningful at discontinuities, where linearized Fourier analysis alone cannot establish the corresponding nonlinear stability.

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.