Finite difference method

A finite difference method replaces derivatives by algebraic combinations of values on a discrete space-time grid.

Averaging a smooth function at the four points $x\pm h{e}_1,x\pm h{e}_2$ cancels odd Taylor expansion terms and leaves $h^2\Delta u/4$. The fourth-order correction is $h^4(u_{1111}+u_{2222})/48$. Equality to the center value at every sufficiently small radius therefore forces $\Delta u=0$ at that point, but does not establish harmonicity throughout a neighbourhood.

A parabolic refinement keeps the time step proportional to the square of the spatial step: $k=r{d}^2$ for a fixed positive diffusion Courant number $r$. A consistency bound $O(k+d^2)$ then becomes $O(d^2)$, though its separate time order remains one.

For the advection equation $u_t+c{u}_x=0$ on a uniform grid, the Lax-Wendroff update isIts amplification factor is $G=1-i\nu \sin \theta +{\nu }^2(\cos \theta -1)$, with $|G{|}^2=1-4{\nu }^2(1-{\nu }^2){\sin }^4(\theta /2)$. Thus Von Neumann stability analysis on an infinite or periodic grid gives stability for $|\nu |\le 1$. The scheme has second order of a numerical method for sufficiently smooth solutions; finite-interval boundary conditions still need separate treatment.

Evaluating a variable diffusion coefficient at both adjacent half-grid points makes odd spatial Taylor terms cancel. The conservative spatial bracket is $h^2(a{u}_x{)}_x+O(h^4)$, yielding a fourth-order one-step residual when the forward time step is proportional to $h^2$.

For constant-speed advection $u_t+a{u}_x=0$ with $a>0$, the explicit upwind scheme uses $U_j^{n+1}=(1-\nu )U_j^n+\nu U_{j-1}^n$, $\nu =ak/h$. For $a<0$ the neighbour must be on the other side. When $0\le \nu \le 1$, the update is a convex combination and is contractive in the maximum norm on a periodic grid, or with appropriate controlled inflow data.

On a periodic grid, Von Neumann stability analysis gives $G=1-\nu +\nu e^{-i\theta }$ and $|G{|}^2=1-2\nu (1-\nu )(1-\cos \theta )$, proving the same exact contraction range in the discrete L2 norm. Boundary estimates are still required on a finite interval; stable scalar eigenvalues of a triangular boundary update alone need not give a mesh-uniform bound.

On a square mesh, the dimensionless five-point Laplacian is ${\Gamma }_{5}u=u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}$. The scaled operator $h^{-2}{\Gamma }_5$ approximates the Laplacian with error $O(h^2)$ for sufficiently smooth functions.

A nine-point finite-difference stencil on a square grid couples a central value to its four axial neighbours and four diagonal neighbours. Ordering the grid by columns produces a block tridiagonal matrix whose diagonal and off-diagonal blocks are symmetric tridiagonal Toeplitz matrices.

For a smooth solution of $\Delta u=f$, the corrected formula ${\Gamma }_{9}u=h^{2}f+(h^2/12){\Gamma }_{5}f$ has unscaled residual $O(h^6)$ and normalized local truncation error $O(h^4)$. The inverse interior matrix has Euclidean operator norm $O(h^{-2})$; on $O(h^{-2})$ grid points this gives unweighted Euclidean nodal error $O(h^3)$. Uniform derivatives up to the boundary are required for a grid-independent remainder bound.

Fourier stability analysis transforms a translation-invariant finite-difference scheme into scalar or matrix recurrences indexed by wavenumber. The Parseval identity then converts uniform multiplier bounds into discrete $2$-norm bounds.

The centered unit-step approximation to a second derivative is

On a uniform grid of spacing $h$, the second-order central differenceapproximates $f^{\prime \prime }(x)$ with error $O(h^2)$ for a sufficiently smooth function.

The forward difference of a sequence or function on the integers is

On an $N$-pixel line with last forward difference zero, $D^{\ast }p$ has first entry $-p_1$, interior entries $p_{i-1}-p_i$, and last entry $p_{N-1}$, with $D^{\ast }=0$ for $N=1$. On a square grid add these expressions along rows and columns. Unused last-edge components contribute nothing. The identity $⟨Du,p⟩=⟨u,D^{\ast }p⟩$ fixes every boundary sign. On the unscaled square grid, $\Vert D{\Vert }^2=8{\cos }^2(\pi /(2N))\le 8$ for $N>1$.

A discrete antiderivative of $f$ is a function $F$ satisfying $\Delta F=f$. Pascal's identity gives

If ${\Delta }^{k+1}f=0$, then $f$ is an integer linear combination of $1,\left(\genfrac{}{}{0ex}{}{n}{1}\right),\dots ,\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ whenever $f$ is integer-valued. Successive differences recover the coefficients.

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_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.

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.

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.

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$.

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$.

The Courant number is the dimensionless ratio of physical propagation during one time step to one spatial grid spacing, commonly $\mu =c\Delta t/\Delta x$.

The diffusion Courant number is $r=\nu k/d^2$ for diffusivity $\nu$, time step $k$ and spatial step $d$. Keeping $r$ fixed is a parabolic mesh refinement. For the Forward Euler diffusion scheme it is bounded by $1/2$ on a uniform one-dimensional grid; the Backward Euler diffusion scheme has no positive upper stability restriction.

A Courant–Friedrichs–Lewy condition requires the numerical domain of dependence to contain the differential equation's physical domain of dependence. It commonly bounds a ratio such as $c\Delta t/\Delta x$ for an explicit hyperbolic scheme.

On $J$ interior points, the one-dimensional centered second-difference matrix with homogeneous Dirichlet boundaries has eigenvalues

If $L$ is the negative semidefinite Dirichlet discrete Laplacian, the Crank-Nicolson amplification matrixhas eigenvalues $(1+\mu {\lambda }_j/2)/(1-\mu {\lambda }_j/2)$ of modulus at most one for every $\mu \ge 0$.

If $T$ is the one-dimensional Dirichlet second-difference matrix on $m$ interior points, the two directional parts of the two-dimensional five-point Laplacian areThey commute, and $A_x+A_y$ is the Kronecker-sum discretization of the Laplacian.

The split stephas amplification matrixOn a common directional eigenvector its eigenvalue is $(1+\mu {\lambda }_p)/(1-\mu {\lambda }_q)$.

On the $m$ by $m$ Dirichlet grid with $h=1/(m+1)$, the exact discrete stability condition for $m>1$ isThe mesh-independent sufficient condition is $0<\mu \le 1/2$, and the exact bound tends to $1/2$ as $h\to 0$.

Adding the correctiongivesIts common-basis eigenvalues areSince ${\lambda }_p,{\lambda }_q<0$, one has $0<d_{pq}\le 1$ for every $\mu >0$.

On a three-dimensional Cartesian grid, the seven-point Dirichlet Laplacian adds the six nearest-neighbour values and subtracts six times the central value. Extending a grid vector by zero on the boundary gives the discrete energy identityso its matrix is symmetric negative definite.

With homogeneous Dirichlet boundary conditions, centered differences turn one-dimensional diffusion into a symmetric negative-definite matrix and constant advection into a skew-symmetric matrix. Their sum dissipates the discrete Euclidean norm, giving mesh-uniform energy stability for every fixed advection coefficient.

A discrete elliptic operator satisfies a discrete maximum principle when a grid function with nonnegative discrete Laplacian cannot have a positive interior maximum unless it is constant. It yields uniqueness and max-norm stability estimates for Dirichlet finite-difference problems.

The method of lines discretizes space in a time-dependent partial differential equation while leaving time continuous, producing a finite system of ordinary differential equations.

Let $K_h$ be a symmetric nonnegative matrix and consider $U^{\prime \prime }=(\alpha I-K_h)U$ with fixed real $\alpha$. In any inner product for which $K_h$ is self-adjoint, diagonalization of a matrix givesIndeed each modal coefficient is either a cosine and sine divided by its frequency, a linear function at zero frequency, or a hyperbolic cosine and sine divided by its growth rate. The bounds $|\cos (\omega t)|\le 1$, $|\sin (\omega t)/\omega |\le t$, and $\sinh (\beta t)/\beta \le t{e}^{\beta t}$ prove the estimate. It is uniform in the spectrum of $K_h$ and therefore gives finite-time displacement stability. Uniform velocity estimates need additional control of the initial spatial energy. All-time displacement bounds require a strictly positive lower bound for $K_h-\alpha I$, a different condition.

For $U^{\prime \prime }=(\alpha I-K_h)U$ with symmetric positive definite $K_h$, all-data displacement is bounded for all time on a fixed grid exactly when $\alpha <{\lambda }_{\min }(K_h)$. Diagonalization of a matrix reduces this assertion to oscillators. At equality, an initial velocity in a zero-frequency mode produces linear growth; above it a growing hyperbolic mode is available. A lower bound on ${\lambda }_{\min }(K_h)-\alpha$ uniform in the mesh gives a uniform displacement estimate. This all-time statement is stronger than displacement stability of a symmetric semidiscrete wave equation on fixed finite time intervals.

If a spatial method of lines has a Hermitian matrix $H$, its generator $-iH$ is a skew-Hermitian matrix. The matrix exponential $e^{-itH}$ is a unitary matrix, so it preserves the Euclidean norm and any constant-volume discrete L2 norm. Equivalently, differentiating the squared norm gives $2\mathrm{Re}(-i{U}^{\ast }HU)=0$. A real sampled potential preserves the Hermitian property of a symmetric discrete Laplacian. Subsequent time discretization must be assessed separately.

For a well-posed linear initial-value problem and a consistent finite-difference approximation, stability is equivalent to convergence.