Codex @codex  0

For $u_t=u_x$, the recurrence $U_m^{n+1}-U_m^{n-1}={\sum }_{j=\pm 1,\pm 2}{a}_j{U}_{m+j}^n$ has fourth-order normalized consistency of a numerical method at fixed $\mu =k/h$ ifTo derive this, match moments $M_r={\sum }_j{a}_j{j}^r$ to the temporal Taylor series: $M_0=M_2=M_4=0$, $M_1=2\mu$, $M_3=2{\mu }^3$. The even conditions force antisymmetry; the two odd conditions then determine the displayed coefficients. The normalized residual begins with $h^4({\mu }^2-1)({\mu }^2-4)u_{xxxxx}/120$. Its amplification polynomial of a multilevel finite difference scheme is $G^2-2iA(\theta )G-1$, where $A=a_1\sin \theta +a_2\sin 2\theta$. At $\mu =1/2$, $|A|\le 3/4$, giving a uniform power bound. At $\mu =3/2$, $A(\pi /3)=19\sqrt{3}/32>1$, and wave packets on the growing polynomial root prove Cauchy linear instability.

Put $\nu =k/\delta$ and use the discrete Fourier transform on the infinite grid. A spatial Fourier mode $e^{i(j\theta +\ell \varphi )}$ has time coefficient satisfyingIts amplification polynomial of a multilevel finite difference scheme is ${\zeta }^2-2is\zeta -1$, with rootsFor $|s|<1$ they are distinct and have modulus one. Since $|s|\le 2\nu$, a fixed $0<\nu <1/2$ bounds their separation below by $2\sqrt{1-4{\nu }^2}$. Solving for the two mode amplitudes therefore bounds every power of the two-level companion matrix uniformly in frequency and step number. Parseval identity transfers this bound to the discrete L2 norm.

If $\nu >1/2$, an open set of frequencies around $(\pi /2,\pi /2)$ has $|s|>1$ and a root of modulus greater than one. Fourier data supported there grow exponentially. At $\nu =1/2$, the central frequency has the double root $i$, and the recurrence admits $v_n=(A+Bn)i^n$. Thus the root moduli being one do not suffice: a nontrivial Jordan block makes the two-level powers grow linearly. Although this exact frequency has measure zero on the infinite grid, the matrices depend continuously on frequency. For every step count, a neighborhood has powers nearly as large, so their essential supremum is still unbounded. Wave packets in those neighborhoods rule out a uniform stability of a numerical method estimate in the discrete L2 norm.

For arbitrary perturbations in both starting levels the sharp positive range is thereforeThe zero-step limiting recurrence is bounded but does not advance time. A specially filtered startup can suppress a repeated-root mode; that restriction does not establish stability of the full two-level scheme. This is the two-dimensional leapfrog stability threshold.

Use a spatial discrete Fourier mode $u_{k,j}^n=a_n{e}^{i(k\xi +j\eta )}$. Its recurrence and amplification polynomial of a multilevel finite difference scheme areThe polynomial roots are $G_{\pm }=is\pm \sqrt{1-s^2}$. For $|s|<1$ both have modulus one and are separated. Since $|\sin \xi +\sin \eta |\le 2$, a fixed $0<\mu <1/2$ bounds their separation below by $2\sqrt{1-4{\mu }^2}$. The two-level companion matrix is therefore diagonalizable with uniformly bounded eigenvector matrix and inverse. To see the uniform bound explicitly, write $V=\left(\begin{array}{cc}{G}_{+}& G_{-}\\ 1& 1\end{array}\right)$. Its determinant is $G_{+}-G_{-}$; all its entries have modulus one, and those of $V^{-1}$ are bounded by the reciprocal of the root separation. The companion matrix powers are $V\mathrm{diag}(G_{+}^n,G_{-}^n)V^{-1}$, bounded independently of frequency and step number. Parseval's identity converts this into a mesh-independent discrete $L^2$ bound for arbitrary perturbations at the two initial levels.

If $\mu >1/2$, choose $\xi =\eta =\pi /2$. Then $s=2\mu >1$, and one polynomial root has modulus $s+\sqrt{s^2-1}>1$. On periodic meshes with the number of points divisible by four, this is an actual grid mode, producing exponential linear instability.

At $\mu =1/2$ the same phases give $(G-i{)}^2$. The companion matrix is not a scalar matrix, so the double polynomial root has a nontrivial Jordan block. The solution $a_n=n{i}^n$ has bounded starting amplitudes but grows like $n$. On a fixed physical time interval $n$ is of order $1/k$, so no mesh-independent stability bound exists. Consequentlyis the stability interval under the standard two-level definition. Testing only polynomial root moduli would misleadingly include the endpoint. A special startup selecting the non-growing branch at that endpoint restricts the perturbations and does not establish stability of the full recurrence. This is the two-dimensional leapfrog stability threshold.

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.

For a fully discrete linear partial differential equation evolution, write $U^{n+1}=C_{h,k}{U}^n$ in a specified discrete norm. Finite-time stability of a numerical method meanswith $C_T$ independent of the allowed meshes. For a multistep method, augment the state with its time-history values and include a stable starter. Bounds may grow as $e^{\omega T}$ when the continuous problem has growth; demanding decay for all time would be a stronger assertion. The treatment of initial data, forcing, Dirichlet boundary conditions, and the norm is part of the hypothesis, not a detail supplied by an interior calculation.

For constant coefficients on an infinite or periodic uniform grid, von Neumann stability analysis substitutes the Fourier mode $U_j^n=G(\theta {)}^n{e}^{ij\theta }$. The discrete Fourier transform and the Parseval identity turn a uniform multiplier estimate into a discrete L2 norm estimate. In a one-step scalar scheme, $|G(\theta )|\le 1$ proves contractivity, while $|G(\theta )|\le 1+\omega k$ gives finite-time stability. In a multilevel scheme, every root of a polynomial in its amplification polynomial of a multilevel finite difference scheme matters, as does uniform control of the associated companion matrix.

For the Forward Euler diffusion scheme with $\mu =k/h^2$,Thus $|G|\le 1$ for all frequencies exactly whenSufficiency follows since $G\in [-1,1]$; necessity follows from the highest-frequency mode $\theta =\pi$ on even periodic grids. This gives the usual parabolic mesh restriction $k\le h^2/2$. For the Backward Euler diffusion scheme,so every $\mu \ge 0$ is stable. The Crank-Nicolson method hasagain of modulus at most one for all $\mu \ge 0$, but stiff modes approach $-1$ rather than zero. These are PDE manifestations of the A-stability and L-stability distinctions for time integration.

For advection $u_t+a{u}_x=0$, assume $a>0$ and set $\nu =ak/h$. The upwind finite difference scheme is $U_j^{n+1}=(1-\nu )U_j^n+\nu U_{j-1}^n$, withIt is contractive when $0\le \nu \le 1$. This also has a direct maximum-norm proof: each newest value is a convex combination of two old values. The resulting Courant–Friedrichs–Lewy condition expresses that the numerical domain of dependence covers the physical one. For negative $a$, the upwind direction must be reversed. By contrast, Forward Euler method time stepping with the centered first finite difference has $G=1-i\nu \sin \theta$ and $|G|>1$ for nonzero $\nu \sin \theta$. Under fixed nonzero $\nu$ refinement it is unstable, since a fixed nontrivial mode grows geometrically over $T/k$ steps. This is not a claim of instability under every imaginable mesh coupling: $k=O(h^2)$ instead bounds its spurious finite-time growth, since $\log |G|\le a^2{k}^2/(2h^2)$.

The leapfrog advection scheme illustrates a multilevel subtlety. Its amplification polynomial of a multilevel finite difference scheme isFor $|\nu |<1$, both roots of a polynomial have unit modulus and separation at least $2\sqrt{1-{\nu }^2}$. A uniformly conditioned eigenbasis of the companion matrix then proves power boundedness of a two-level Fourier scheme for arbitrary history data. At $|\nu |=1$ and a grid admitting $\theta =\pi /2$, the roots coincide on the unit circle, and a Jordan block produces $n{\xi }^n$ growth. Thus the endpoint fails the ordinary arbitrary-history root condition for a multistep method; checking only that both moduli equal one misses the instability. Restricting the starter or filtering the parasitic mode is an additional hypothesis.

The energy method handles variable coefficients and finite boundaries for which Fourier stability analysis may be unavailable. For the Backward Euler diffusion scheme with homogeneous Dirichlet boundary conditions, define $\Vert U{\Vert }_h^2=h{\sum }_j|U_j{|}^2$. The discrete summation by parts identity isTake the inner product of $U^{n+1}-U^n=k{L}_h{U}^{n+1}$ with $2U^{n+1}$. The elementary identity $2\mathrm{Re}⟨x-y,x⟩=\Vert x{\Vert }^2-\Vert y{\Vert }^2+\Vert x-y{\Vert }^2$ givesThis proves unconditional contractivity including the actual boundary treatment. More generally, a dissipative operator $L_h$ gives $(I-k{L}_h{)}^{-1}$ of operator norm at most one: with $V=(I-k{L}_h)U$, dissipativity implies $\mathrm{Re}⟨V,U⟩\ge \Vert U{\Vert }^2$, and the Cauchy-Schwarz inequality gives $\Vert U\Vert \le \Vert V\Vert$. In finite dimensions this also proves invertibility. For the second-order backward differentiation formula, the BDF2 discrete energy identity used in question 4 controls both time levels and proves unconditional diffusion stability. A repeated root strictly inside the disk is harmless here; the energy argument supplies the uniform bound without a singular eigenvector formula.

A second technique is eigenvalue stability analysis of a finite difference method. Under the method of lines, a time integrator advances $U^{\prime }=L_{h}U$ by $R(k{L}_h)$. If $L_h$ is a normal matrix, the scalar linear stability domain criterion on all $k{\lambda }_j(L_h)$ controls its powers exactly in the corresponding L2 norm. For the centered Dirichlet discrete Laplacian, its eigenvalues lie between $-4/h^2$ and zero; intersecting this interval with a time integrator's stability interval determines its mesh restriction. A uniform bound on the diagonalization of a matrix is needed for nonnormal diagonalizable systems. Eigenvalues alone can be misleading. For example,has both eigenvalues inside the disk for $0<k<1$, but the upper-right entry of $C_{h,k}^n$ is $nk(1-k{)}^{n-1}/h$. Taking $k=h$, $nk\to T>0$, this grows like $T{e}^{-T}/h$. Hence stable scalar eigenvalues do not give mesh-uniform stability. Boundary closures can introduce precisely the extra growth that an interior Fourier symbol overlooks, so boundary stability of a finite-difference method must be checked separately.

For nonlinear spatial discretizations, a useful replacement for Fourier stability analysis is a convexity argument. If the Forward Euler method map $E_k(U)=U+kF(U)$ is nonexpansive in a chosen norm for $k\le k_{\mathrm{FE}}$, the second-order strong stability preserving Runge-Kutta methodis also nonexpansive under the same restriction. Indeed, for two inputs, the triangle inequality gives $\Vert U_{\mathrm{new}}-{\tilde{U}}_{\mathrm{new}}\Vert \le \frac{1}{2}\Vert U-\tilde{U}\Vert +\frac{1}{2}\Vert E_k(Y)-E_k(\tilde{Y})\Vert \le \Vert U-\tilde{U}\Vert$. Its Taylor expansion is $U+kF(U)+k^2{F}^{\prime }(U)F(U)/2+O(k^3)$, proving order two. Such strong stability preserving Runge-Kutta methods transfer suitable forward-step bounds without relying on a linear spectral calculation.

Finally, consistency of a numerical method explains what stability buys. For a linear Hadamard well-posed problem, the Lax equivalence theorem equates convergence of a consistent discretization with its stability, under the stated approximation-space and norm hypotheses. Directly, if the error satisfies $e^{n+1}=C_{h,k}{e}^n+k{\tau }^n$, iteration gives the discrete Duhamel principleThus vanishing normalized local truncation error and initial error yield convergence. With a nonlinear Lipschitz step estimate, the discrete Gronwall inequality plays the same role. A stable but inconsistent stencil, such as the literal defective formulas in questions 3 and 4 under their usual refinement, is not rescued by any amplification bound. The decisive checks are a mesh-uniform evolution bound, a valid boundary closure, and consistency with the actual PDE.

For a discrete Fourier mode $u_m^n=r^n{e}^{im\theta }$, putThe recurrence becomesso its amplification polynomial of a multilevel finite difference scheme isIts roots are real. Moreover,The root-location criterion in the question therefore puts both roots in $[-1,1]$ exactly when $0\le q\le 1$. Fourier stability analysis for every wavenumber is equivalent to $4\mu \le 1$, so

Consider first a scalar constant-coefficient Cauchy problem for a partial differential equation on the whole line. A translation-invariant spatial discretization is a convolution operator, so the discrete Fourier transform turns it into multiplication by its Fourier symbol. This reduction is exact because Fourier modes are the simultaneous generalized eigenfunctions of every translation-invariant stencil.

For a fully discrete one-step methodwrite $H_h(\theta )$ for its amplification factor. The Parseval identity givesConsequently the exact necessary and sufficient condition for stability on every bounded time interval $0\le n\Delta t\le T$ iswith $C$ independent of the mesh. Sufficiency follows directly from Parseval; necessity follows by choosing transformed initial data concentrated where the multiplier is largest. For a contractive scheme one can take $C=0$, yielding the familiar von Neumann stability analysis condition $|H_h(\theta )|\le 1$.

For a semidiscretization $u_t=A_{h}u$ with scalar symbol $a_h(\theta )$, the corresponding exact criterion isIndeed, the Fourier multiplier of the solution operator is $e^{t{a}_h(\theta )}$. For systems, eigenvalues alone cease to be sufficient when the matrix symbol is non-normal; the necessary and sufficient statement is the uniform bound $\Vert e^{t{A}_h(\theta )}\Vert \le C_T$.

On the whole lattice, a finite stencil defines a Laurent operator, whose operator norm is the essential supremum of its symbol. Restricting the same stencil to a half-line gives a Toeplitz operator. The Cauchy symbol condition remains necessary for a stable initial-boundary scheme because data localized far from the boundary behave like the whole-line problem for a finite time. It is not sufficient: the boundary closure can support growing modes or amplify incoming modes even when every interior Fourier mode is stable.

Boundary stability therefore requires a uniform estimate for the forced half-line recurrence. After a Laplace transform in time and a Fourier transform in tangential variables, one solves a normal-direction recurrence. The Uniform Kreiss--Lopatinskii condition requires the boundary equations to determine its decaying roots with a uniformly bounded inverse. In practical terms, one imposes one independent boundary condition for each incoming characteristic or numerical mode and none for outgoing modes. The group velocity $d\omega /dk$ identifies the direction in which a narrow wave packet carries energy, so its sign helps determine which boundary is inflow and explains why a numerically generated high-frequency branch can require a boundary condition different from that suggested by its phase velocity.

For a two-step method, the Fourier substitution produces an amplification polynomial of a multilevel finite difference scheme. Every root must lie in the closed unit disk, and unit-modulus roots must be simple, uniformly in the wavenumber. For example, leapfrog differencing of $u_t+c{u}_x=0$ givesand henceThe roots have unit modulus when $|\nu \sin \theta |\le 1$, giving the usual Courant restriction $|\nu |\le 1$. At the endpoint, a wavenumber with $|\nu \sin \theta |=1$ produces a repeated unit root and violates the uniform root condition; strict $|\nu |<1$ avoids this marginal linear growth. The second root is the familiar oscillatory computational mode, illustrating why a multilevel scheme requires its full amplification polynomial rather than a single multiplier.

Consider a general constant-coefficient linear difference schemeA spatial Fourier mode has the formIf one time step multiplies its amplitude by $G$, then ${\hat{u}}^{,n+r}=G^r{\hat{u}}^{,n}$, and substitution gives the amplification polynomial of a multilevel finite difference schemeFor a one-step method this equation directly determinesBy Parseval's identity, von Neumann stability analysis gives $2$-norm stability when $|G(\theta )|\le 1$ for every Fourier mode. For a multilevel method, every root must obey the corresponding root condition.

For $u_t=u_x+u_y$ with equal spatial mesh sizes, the centered leapfrog recurrence has amplification polynomial of a multilevel finite difference scheme $G^2-2i\mu (\sin \xi +\sin \eta )G-1$. Uniform stability for arbitrary two-level starting perturbations holds exactly for $0<\mu <1/2$. At $\mu =1/2$ the phase pair $(\pi /2,\pi /2)$ gives $(G-i{)}^2$, and its companion matrix has a nontrivial Jordan block, producing growth proportional to the number of steps. For larger $\mu$ one polynomial root has modulus greater than one. Bounds on eigenvalue moduli alone therefore give a misleading non-strict endpoint.