von Neumann stability analysis

Von Neumann analysis inserts the Fourier mode $u_m^n=G^n{e}^{im\theta }$ into a constant-coefficient difference scheme. A one-step method is stable in the discrete $2$-norm when its amplification factor satisfies $|G(\theta )|\le 1$ for every resolvable wavenumber.

For a system or a multilevel method, each Fourier mode evolves through an amplification matrix $G(\theta )$. Stability requires bounds on $G(\theta {)}^n$ uniform in frequencies and mesh parameters. Eigenvalue moduli alone do not suffice: $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has unit eigenvalues but $G^n$ has off-diagonal entry $n$. Uniformly bounded eigenvector matrices, or a uniformly positive energy symmetrizer, can supply the missing power bound.

For a uniformly bounded two-level Fourier companion matrix $C(\theta )$ with distinct roots $g_1,g_2$ of modulus at most one, its powers satisfyA positive frequency-uniform gap therefore gives a uniform bound on every power. The Parseval identity then supplies spatial L2 norm stability. When the gap closes at a repeated unit root, a Jordan block can produce growth proportional to the number of steps despite both roots having modulus one.

For the periodic-grid recurrencethe exact von Neumann stability analysis range at fixed $\mu >0$ is $\mu \ge 3$. The Fourier symbol gives $A(\theta ){\xi }^2-4\xi +1=0$, $A=3e^{2i\theta }+8\mu {\sin }^2(\theta /2)$. The root near one obeys $|\xi {|}^2=1+(6-2\mu ){\theta }^2+O({\theta }^4)$, proving instability for $\mu <3$.

For $\mu \ge 3$, $\mathrm{Re}A-3=6(\cos \theta -1{)}^2+4(\mu -3)(1-\cos \theta )\ge 0$. Writing the modal recurrence as $3v^{n+2}-4v^{n+1}+v^n=-(A-3)v^{n+2}$, the BDF2 discrete energy identity proves nonincrease of $|v^{n+1}{|}^2+|2v^{n+1}-v^n{|}^2$. Summing modes proves stability uniformly even when $\mu$ varies within $[3,\infty )$.

This shifted scheme is not a consistent standard heat discretization under $k=\mu h^2$: dividing its extra spatial-shift defect by $2k/3$ produces $3h{u}_x/k+O(h^2/k)$. On a finite interval with Dirichlet boundary conditions it also requires an extra boundary closure. Thus the periodic stability theorem is not a theorem of convergence to the heat equation or of stability for an unspecified boundary treatment.

For a translation-invariant difference operator $(Lu{)}_m={\sum }_j{a}_j{u}_{m+j}$, its Fourier symbol is $\ell (\theta )={\sum }_j{a}_j{e}^{ij\theta }$. Applying $L$ to a Fourier mode $e^{im\theta }$ multiplies that mode by $\ell (\theta )$. Symbols turn constant-coefficient stencil equations into scalar algebraic equations and expose their von Neumann stability analysis.

For a two-level Fourier recurrence ${\hat{u}}^{n+1}=a(\theta ){\hat{u}}^n+b(\theta ){\hat{u}}^{n-1}$, the amplification matrix isStability requires its powers to be uniformly bounded in the time index and the mesh frequencies. If its two amplification roots have modulus at most one and their separation has a positive mesh-independent lower bound, the eigenvectors $({\xi }_j,1{)}^T$ give a uniformly bounded diagonalization of a matrix, proving stability. A repeated unit-modulus root of this companion matrix instead produces a Jordan block and linear growth in time. Thus checking only the moduli of the roots is insufficient.

For $U_m^{n+1}=(1-2\mu )(U_m^n-U_{m+1}^n)+U_{m+1}^{n-1}$, write an amplification root as $G=e^{i\theta /2}\lambda$. Its equation is ${\lambda }^2+2i(1-2\mu )\sin (\theta /2)\lambda -1=0$. The roots have modulus one and a positive frequency-uniform separation exactly for $0<\mu <1$, giving L2 norm stability through a uniformly bounded diagonalization. At either endpoint the Nyquist symbol has a repeated unit root and a Jordan block, producing linear growth. The normalized local truncation error is generally second order; at $\mu =1/2$ and $1$ exactly sampled characteristics solve the recurrence, although the latter parameter is unstable for arbitrary perturbations.

For $u_t=u_{xx}+\alpha u_x$ on the whole line, apply centered differences in space and the Forward Euler method in time. With $r=\Delta t/(\Delta x{)}^2$ and $c=\alpha \Delta t/\Delta x$, the amplification factor isPutting $s={\sin }^2(\theta /2)$ givesFor $\Delta t>0$, the affine expression in braces is nonpositive on $[0,1]$ exactly when $r\le 1/2$ and $c^2\le 2r$. Hence von Neumann stability analysis givesThese conditions are weaker than requiring every stencil weight to be nonnegative; stability in the discrete L2 norm need not imply monotonicity.

A Laurent operator on the whole integer lattice is a translation-invariant bi-infinite matrix. Under the discrete Fourier transform it becomes multiplication by its symbol, so its $2$-operator norm is the essential supremum of the symbol's modulus.

A Toeplitz operator is a constant-diagonal operator on a one-sided sequence space. It models the restriction of a translation-invariant stencil to a half-line, where the boundary prevents direct diagonalization by the full-lattice Fourier transform.

Split a smooth symbol into its strictly negative and strictly positive Fourier modes. Their Toeplitz operators $A,B$ exponentiate exactly to the Toeplitz operators of the corresponding exponentials, and $T(e^f)T(e^{-f})=e^A{e}^B{e}^{-A}{e}^{-B}$. The trace of $[A,B]$ is ${\sum }_{n>0}n{f}_n{f}_{-n}$ because the shift relation gives $\mathrm{Tr}[(S^{\ast }{)}^m,S^n]=n{\delta }_{mn}$. Smoothness makes the series converge in trace norm, so the Fredholm determinant of an exponential commutator applies.

On the Hardy space of the circle, a continuous symbol gives a Fredholm Toeplitz operator exactly when it never vanishes. Nonvanishing gives a parametrix $T_{1/f}$ modulo compact semicommutators. Write $f(z)=z^n{e}^{h(z)}$, where $n$ is its counterclockwise winding number and $h$ is continuous. Homotopy through nonvanishing symbols and local constancy of the Fredholm index reduce the index to that of the shift $T_{z^n}$, namely $-n$. If $f$ vanishes at a point, normalized Hardy reproducing kernels concentrating there are weakly null unit vectors on which $T_f$ tends to zero, contradicting a Fredholm parametrix.

Boundary stability is a mesh-uniform estimate for a finite-difference initial-boundary problem with forced boundary data. Stability of the whole-line Cauchy scheme is necessary but may fail to control boundary-supported numerical modes.

A boundary closure supplies the boundary rows or ghost values needed by a spatial finite difference method. It must reflect the PDE boundary data and preserve a mesh-uniform stability estimate. Interior von Neumann stability analysis cannot detect an arbitrary unstable boundary row; for example a row $U_1^{n+1}=2U_1^n$ produces growth even beside an otherwise stable heat stencil.

After transforming time and tangential variables, the uniform Kreiss--Lopatinskii condition requires the boundary equations to determine every decaying normal mode with an inverse bounded uniformly over transformed frequencies. It rules out growing or weakly controlled boundary modes.

For a normal amplification matrix $Q$, the bound $\rho (Q)\le 1$ with semisimple unit-modulus eigenvalues controls every power $Q^n$. For a nonnormal family, eigenvalues alone do not control transient growth or give a mesh-uniform bound; eigenvector conditioning, pseudospectra, or a direct energy estimate is then needed.

Stability of the partial differential equation requires a bound uniform over refining meshes and time steps on every fixed physical time interval. This is distinct from demanding an all-time bound for one fixed matrix. A fixed matrix is power bounded exactly when all eigenvalues are in the closed unit disk and its unit-modulus Jordan blocks have size one. A time-step family can nevertheless have harmless polynomial physical-time growth: $G_{\tau }=I+\tau N$ with $N^2=0$ gives $G_{\tau }^n=I+n\tau N$, uniformly bounded for $n\tau \le T$ despite its nontrivial unit-modulus Jordan block. Mesh dependence and the scaling of the step must therefore be kept in an eigenvalue argument.

Every eigenvalue of this nonnormal matrix is $-1$, but $e^{t{L}_h}=e^{-t}\left(\begin{array}{cc}1& t/h\\ 0& 1\end{array}\right)$. Applied to the second coordinate unit vector, its norm at $t=1$ is at least $1/(eh)$. Thus negative eigenvalues give no mesh-uniform semigroup bound by themselves. Normal matrices, uniformly conditioned eigenvector bases, or a direct dissipative energy estimate supply the missing control.

On a finite inflow grid, explicit upwinding has a triangular amplification matrix whose only eigenvalue is $1-\mu$. For $1<\mu <2$ this eigenvalue lies inside the unit disk, although interior high-frequency data are amplified by $|1-2\mu |>1$ before reaching the boundary. This shows why eigenvalue analysis without uniform normality can give the wrong stability range.

For the centered spatial second difference, backward Euler has amplification factor $G(\theta )=[1+4\mu {\sin }^2(\theta /2){]}^{-1}$ and is stable for every physical Courant number $\mu \ge 0$.

For a smooth solution of the heat equation, backward Euler with the centered second difference has residual divided by $k$ equal to $-(k/2+d^2/12)u_{xxxx}+O(k^2+k{d}^2+d^4)$. Both leading terms have the same sign. Under parabolic mesh refinement their coefficient is $-d^2(r/2+1/12)$, which cannot vanish for a positive diffusion Courant number. Consistency of an update $U^{n+1}-U^n=\alpha (r){\delta }^2{U}^{n+1}$ forces $\alpha (r)=r$ for every fixed positive $r$.

On $J$ interior Dirichlet grid points, backward Euler diffusion is stable for $\mu \ge 0$ and also for the nonphysical branch $\mu \le -[2{\sin }^2(\pi /(2(J+1))){]}^{-1}$. The latter makes every amplification denominator at most $-1$.

For the centered spatial second difference, forward Euler has the updateIt is stable for $0\le \mu \le 1/2$. In this range the update is a convex combination, which gives a direct maximum-norm proof as well as the Fourier proof.

Let $D_h$ be the scaled Dirichlet discrete Laplacian and $V_h$ a bounded diagonal reaction matrix. If $0\le k/h_x^2\le 1/2$, the heat update $H_h=I+k{D}_h$ is contractive in both the maximum norm and the mesh-weighted L2 norm. For $\Vert V_h\Vert \le A$,This proves finite-time mesh-uniform stability without requiring the full update to have nonnegative entries or to be contractive. Negative reaction terms can destroy nonnegativity at the endpoint $k/h_x^2=1/2$, so it is the heat part alone that is treated as a nonnegative contraction.

For the schemewhere $L=[1,-2,1]$ is the Dirichlet second-difference matrix and $\mu >0$, the modal amplification factor isMesh-uniform stability holds exactly whenFor consistency with $u_t=u_{xx}$ under the displayed normalization one additionally chooses $a=1/2$; stability then requires $c\le 1/2$.

For a one-step translation-invariant recurrence, the amplification factor is the multiplier satisfyingThe Parseval identity converts the pointwise bound $|H|\le 1$ into non-growth of the discrete $2$-norm.

Forthe Fourier convention $\hat{u}(\theta )={\sum }_m{e}^{-im\theta }{u}_m$ givesprovided the denominator does not vanish.

For a schemethe Fourier ansatz gives the amplification polynomialIts roots are the amplification factors of that Fourier mode.

Centered space and leapfrog time discretization of $u_t+c{u}_x=0$ givesIts amplification polynomial is $G^2+2i\nu \sin \theta G-1$. The two roots have unit modulus for $|\nu \sin \theta |\le 1$, but a repeated unit root at equality violates the uniform 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.

For the centered two-level advection recurrence, the Fourier roots at Courant number one coalesce at $i$ for phase $\pi /2$. The resulting nontrivial Jordan block admits linearly growing solutions despite both eigenvalues having modulus one. Uniform two-level stability for arbitrary starting perturbations therefore requires the strict inequality $0<\mu <1$. A specially consistent startup can select the exact translating branch at the boundary, which is a restriction on the data rather than a uniform bound on the full update.

Applying a centered leapfrog step in time and the centered second difference in space to $u_t=u_{xx}$ givesFor every $\mu >0$, each nonconstant Fourier mode has an amplification root of modulus greater than one, so the scheme is unconditionally unstable.

For the centered spatial second difference, the Crank-Nicolson diffusion scheme has amplification factorIt is stable for every $\mu \ge 0$.

For the two-dimensional heat equation with zero Dirichlet boundary conditions, the Crank-Nicolson diffusion scheme preserves the interval between the minimum and maximum of the initial nodal values and zero under the displayed restriction. The old-time coefficients form a convex combination, while the new-time matrix has positive diagonal, nonpositive off-diagonal entries and row diagonal-dominance margin one. A hypothetical new value above the old maximum makes the new-time row at least that value while its right side is at most the old maximum, a contradiction; apply the reversed inequality at a new minimum. In particular initial nonnegative or nonpositive values preserve their sign. The step restriction guarantees these value bounds even though scalar A-stability holds without it.

Let $D$ be the real skew-symmetric tridiagonal matrix with superdiagonal $1$ and subdiagonal $-1$. The centered-space Crank-Nicolson discretization of $u_t=u_x$ has amplification matrixThe eigenvalues of $D$ are $2i\cos (j\pi /(M+1))$, so those of $Q$ are Cayley transformsThey all have modulus one, and their common orthonormal eigenbasis makes $Q$ normal. Hence the method is stable for every $\mu >0$.

Forthe amplification polynomial isAll modes have unit-modulus roots only when $\rho =1$ and $0\le \mu \le 1$.