Von Neumann stability analysis

In numerical analysis, von Neumann stability analysis (also known as Fourier stability analysis) is a procedure used to check the stability of finite difference schemes as applied to linear partial differential equations. The analysis is based on the Fourier decomposition of numerical error and was developed at Los Alamos National Laboratory after having been briefly described in a 1947 article by British researchers John Crank and Phyllis Nicolson. This method is an example of explicit time integration where the function that defines governing equation is evaluated at the current time. Later, the method was given a more rigorous treatment in an article co-authored by John von Neumann.

The stability of numerical schemes is closely associated with numerical error. A finite difference scheme is stable if the errors made at one time step of the calculation do not cause the errors to be magnified as the computations are continued. A neutrally stable scheme is one in which errors remain constant as the computations are carried forward. If the errors decay and eventually damp out, the numerical scheme is said to be stable. If, on the contrary, the errors grow with time the numerical scheme is said to be unstable. The stability of numerical schemes can be investigated by performing von Neumann stability analysis. For time-dependent problems, stability guarantees that the numerical method produces a bounded solution whenever the solution of the exact differential equation is bounded. Stability, in general, can be difficult to investigate, especially when the equation under consideration is nonlinear.

In certain cases, von Neumann stability is necessary and sufficient for stability in the sense of Lax–Richtmyer (as used in the Lax equivalence theorem): The PDE and the finite difference scheme models are linear; the PDE is constant-coefficient with periodic boundary conditions and has only two independent variables; and the scheme uses no more than two time levels. Von Neumann stability is necessary in a much wider variety of cases. It is often used in place of a more detailed stability analysis to provide a good guess at the restrictions (if any) on the step sizes used in the scheme because of its relative simplicity.

The von Neumann method is based on the decomposition of the errors into Fourier series. To illustrate the procedure, consider the one-dimensional heat equation $${\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}$$ defined on the spatial interval ${\displaystyle L}$, with the notation $${\displaystyle u_{j}^{n}=u(x_{j},t^{n})}$$ where ${\displaystyle x_{j}}$ are the specific x values, and ${\displaystyle t^{n}}$ are the sequence of t values.

We can discretize the heat equation as

where $${\displaystyle r={\frac {\alpha \,\Delta t}{\left(\Delta x\right)^{2}}}}$$

Then the solution ${\displaystyle u_{j}^{n}}$ of the discrete equation approximates the analytical solution ${\displaystyle u(x,t)}$ of the PDE on the grid.

Define the round-off error ${\displaystyle \epsilon _{j}^{n}}$ as $${\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}$$ where ${\displaystyle u_{j}^{n}}$ is the solution of the discretized equation (1) that would be computed in the absence of round-off error, and ${\displaystyle N_{j}^{n}}$ is the numerical solution obtained in finite precision arithmetic. Since the exact solution ${\displaystyle u_{j}^{n}}$ must satisfy the discretized equation exactly, the error ${\displaystyle \epsilon _{j}^{n}}$ must also satisfy the discretized equation. Here we assumed that ${\displaystyle N_{j}^{n}}$ satisfies the equation, too (this is only true in machine precision). Thus