Finite difference

A finite difference compares the values of a function at finitely separated arguments. For example, the forward finite difference with step $h$ isIt is the discrete analogue of a derivative and lowers the degree of a nonconstant polynomial by one.

A Newton series expands a function in the basis of binomial coefficients: $f(x)={\sum }_{k\ge 0}\left(\genfrac{}{}{0ex}{}{x}{k}\right){\Delta }^{k}f(0)$. For a polynomial, the sum terminates.

The finite increment ${\Delta }_{h}f(x)=f(x+h)-f(x)$ for a specified step $h$. Higher forward differences are obtained by iterating this operation. For a polynomial of degree $d$, its $(d+1)$st forward difference is zero.

The first difference of a sequence records the change between consecutive terms. Its partial sums telescope: $u_n=u_0+{\sum }_{j=0}^{n-1}\Delta u_j$. Thus a compact description of the differences can give an efficient way to reconstruct distant terms of the original sequence.

Matching the monomials $1,x,x^2$ gives $f^{\prime \prime }(2)\approx (2/3)f(0)-f(1)+(1/3)f(3)$. With error defined as the true derivative minus this approximation, its Peano kernel is $K(\theta )=21_{\{\theta <2\}}+(1-\theta {)}_{+}^2-(3-\theta {)}^2/3$ on $[0,3]$. The derivative evaluation causes a jump at two. Applying the Peano kernel theorem gives $\lambda (f)=\frac{1}{2}{\int }_0^{3}K{f}^{(3)}$.

Given noisy evaluations at $x_0+h{Z}_i$, where the $Z_i$ are independent Rademacher random variables, the estimatoraverages one-sided finite differences from both directions. If $|g^{\prime \prime }|\le M$ and the noise variance is ${\sigma }^2$, its mean squared error is at most $h^2{M}^2/4+{\sigma }^2/(N{h}^2)$.

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

For the heat equation with homogeneous Neumann boundary conditions at both ends of a uniform grid, even reflection supplies ghost values for a centered finite difference method. A fourth-order centered second derivative needs two ghost values at each end, obtained with the displayed rule for $j=1,2$. For a sufficiently smooth compatible solution, differentiating the Neumann boundary conditions in time and using the heat equation makes every odd spatial derivative vanish at the endpoints. Its even extension therefore retains the smoothness required by the fourth-order stencil; arbitrary incompatible initial data need care at the initial boundary corners.

The discrete cosine modes diagonalize this reflected spatial operator, with eigenvalues $-4h^{-2}{\sin }^2(\theta /2)[1+{\sin }^2(\theta /2)/3]$ and $\theta =j\pi /N$. Its most negative eigenvalue is $-16/(3h^2)$, so the forward Euler method is stable in the corresponding trapezoid-weighted norm when $\Delta t/h^2\le 3/8$. The Crank-Nicolson diffusion scheme is stable for every time step for this closure, but its amplification factor approaches $-1$ on highly damped spatial modes, allowing sign-alternating numerical transients. The ghost values are defined here as auxiliary samples outside the physical interval, rather than additional unknown initial data.

For $u_t=(a{u}_x{)}_x$, this conservative finite difference method uses coefficient samples at edge midpoints. For smooth $a,u$, its spatial residual is $h^2(a{u}_{xxxx}/12+a^{\prime }{u}_{xxx}/6+a^{\prime \prime }{u}_{xx}/8+a^{\prime \prime \prime }{u}_x/24)+O(h^4)$. With zero Dirichlet boundary conditions, $-v^T{L}_{h}v=h^{-2}{\sum }_j{a}_{j+1/2}(v_{j+1}-v_j{)}^2$. If $0<a\le a_{+}$, its eigenvalues lie in $[-4a_{+}/h^2,0]$, giving the mesh-independent forward Euler method bound $\Delta t/h^2\le 1/(2a_{+})$. Constant $a=a_{+}$ proves sharpness as the grid is refined.

Finite difference coefficients are the weights in a linear combination of sampled function values approximating a derivative. Matching polynomial moments determines the weights and approximation order. Forward, backward and centered stencils have different coefficients.

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.

A Mehrstellen method increases the accuracy of a compact finite difference method discretization by using the differential equation to express its leading truncation terms in terms of the prescribed source. For $\Delta u=f$, the square-grid nine-point finite-difference stencil for the Laplacian satisfies $D_{9}u=\Delta u+h^2{\Delta }^{2}u/12+O(h^4)$. Replacing ${\Delta }^{2}u$ by $\Delta f$ and approximating this source derivative gives $D_{9}U=(I+h^2{D}_5/12)f$, with fourth-order normalized consistency of a numerical method while retaining a compact nine-point matrix. The uncorrected nine-point finite-difference stencil operator is only second-order accurate on general functions.

On a square grid, let $C_h$ sum the four axial neighbors and $K_h$ the four diagonal neighbors. With $M_h=2I/3+C_h/12$ and $L_h=(-10I/3+2C_h/3+K_h/6)/h^2$, the method of lines has fourth-order normalized spatial numerical consistency. Its mass symbol is $(4+\cos \xi +\cos \eta )/6\ge 1/3$, and its diffusion symbol is $[-10+4(\cos \xi +\cos \eta )+2\cos \xi \cos \eta ]/(3h^2)\le 0$. Every continuous-time mode therefore contracts, and Parseval identity gives mesh-independent discrete L2 stability. No time-step restriction is involved until a time integrator is selected.

Consider a one-step stencil with new-level coefficients $a=\mu (1+\mu )/2$, $b=(1+\mu )(2-\mu )$ and $c=(1-\mu )(2-\mu )/2$ at offsets $-1,0,1$, and old-level coefficients $d=2-\mu$, $e=1+\mu$ at offsets $0,1$. For $u_t=u_x$, its exact-solution residual starts at fourth degree with coefficient $\mu (\mu -2)(\mu -1)(\mu +1)/24$. Thus its normalized local truncation error is third order at fixed nonzero Courant number, except at the exact-shift values $-1,1,2$; zero gives the zero-step identity. For its Fourier symbol, $|D{|}^2-|N{|}^2=\mu (\mu -2)(\mu -1)(\mu +1)(1-\cos \theta {)}^2$. Uniformly invertible stable steps occur for $\mu \le -1$, $0\le \mu \le 1$ or $\mu \ge 2$.

With zero endpoints, the centered first-difference matrix has opposite off-diagonal entries and is a skew-adjoint operator in the mesh-weighted Euclidean inner product. Adding a real multiple of it to the negative symmetric second-difference matrix does not change the real energy form. The semidiscrete convection-diffusion solution therefore obeys a mesh-independent energy estimate, even when its drift is large compared with its diffusion coefficient.

On a rectangular grid with zero Dirichlet boundary conditions, the nine-point Poisson equation matrix has diagonal $10/3$, axial off-diagonal entries $-2/3$, and diagonal-neighbor entries $-1/6$. Its quadratic form is the sum of squared differences over internal links, weighted by $2/3$ or $1/6$, plus the weighted squares of values linked to the fixed zero boundary. The connected axial grid forces a zero quadratic form to have all entries zero. This proves positive-definite matrix in every ordering. The Jacobi method converges because its nonnegative iteration matrix has row sums at most one, with a strict deficit at boundary-adjacent vertices; a hypothetical unit-modulus eigenvector would propagate its maximum modulus to such a deficit row, which is impossible.

A three-point implicit advection stencil can match translated data through cubic Taylor expansion while retaining a rational Fourier amplification symbol. For the coefficient family $a=\mu (1+\mu )/6$, $b=(2-\mu )(1+\mu )/3$, $c=(2-\mu )(1-\mu )/6$, $d=(2-\mu )/3$, $e=(1+\mu )/3$, its factor is $G=(d+e{e}^{i\theta })/(a{e}^{-i\theta }+b+c{e}^{i\theta })$. Direct subtraction gives $|D{|}^2-|N{|}^2=\mu (\mu -2)(\mu -1)(\mu +1)(1-\cos \theta {)}^2/9$. The displayed ranges are the contraction ranges in the discrete L2 norm, with a nonvanishing denominator there. For forward time and positive speed parameter, retain the nonnegative ranges. For fixed nonzero nonexceptional Courant number, the normalized truncation error is third order; at $\mu =-1,0,1,2$ the update is the corresponding exact grid translation.

For $u_t+c{u}_x=0$ on a uniform grid, the update $U_m^{n+1}=(U_{m-1}^n+U_{m+1}^n)/2-\nu (U_{m+1}^n-U_{m-1}^n)/2$ has symbol $\cos \theta -i\nu \sin \theta$, where $\nu =c\Delta t/h$. Its squared modulus is $1+({\nu }^2-1){\sin }^2\theta$, giving contraction exactly for $|\nu |\le 1$. Stability is separate from the dissipative and dispersive approximation errors.

An arbitrary $L^2$ initial function has no well-defined nodal point samples. Cell averages or a bounded projection to the discrete approximation space supply grid data independently of its pointwise representative. Jensen's inequality bounds the cell-average discrete L2 norm by the continuous L2 norm. Combining consistent bounded initialization, stability, and density of smooth data can extend convergence to rough initial functions, without claiming their nonexistent classical derivatives.

The cell-average projection replaces a function by its average on each grid cell. It is independent of the representative of an L2 function, unlike point sampling. The Cauchy-Schwarz inequality gives $h|U_m{|}^2\le {\int }_{\text{cell }m}|u{|}^{2}dx$, hence $h{\sum }_m|U_m{|}^2\le \Vert u{\Vert }_{L^2}^2$. Viewed as a piecewise constant function, the projection is the orthogonal projection onto the grid's piecewise constant linear subspace. Density of continuous compactly supported functions and the contraction estimate give convergence in the L2 norm as the mesh tends to zero.

On a periodic uniform grid, shifting an index in the discrete inner product gives the displayed skew-adjoint identity for the centered first difference. It also gives $⟨U,D_{xx}V{⟩}_h=-⟨D_{+}U,D_{+}V{⟩}_h$. With finite boundaries there are additional boundary terms unless the prescribed endpoint traces eliminate them. These identities are discrete counterparts of integration by parts and prove energy estimates for finite difference methods.

For $U_t=-cDU$ with $c>0$, a positive discrete norm matrix $H$ satisfying the displayed identity gives $d(U^{T}HU)/dt=-c(U_N^2-U_0^2)$. Homogeneous inflow data $U_0=0$ therefore leave only the nonpositive outflow flux. Inhomogeneous inflow contributes a controlled boundary forcing term. This estimate includes the numerical boundary closure, unlike an interior-only Fourier stability analysis. For any dissipative semidiscrete operator in the $H$ norm, Backward Euler method steps preserve the contraction by the identity $\Vert U^{n+1}{\Vert }_H^2-\Vert U^n{\Vert }_H^2=2k⟨U^{n+1},A{U}^{n+1}{⟩}_H-\Vert U^{n+1}-U^n{\Vert }_H^2$.

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 $0<a(x)\le \beta$, the explicit conservative diffusion update has weights $k{a}_{m-1/2}/h^2$, $1-k(a_{m-1/2}+a_{m+1/2})/h^2$, and $k{a}_{m+1/2}/h^2$. Under the displayed restriction they are nonnegative and sum to one. With zero Dirichlet boundary conditions, each new value is a convex combination of old interior and boundary values, proving contraction in the discrete maximum norm. The spatial matrix is also symmetric with $v^T{L}_{h}v=-h^{-2}\sum a_{m+1/2}(v_{m+1}-v_m{)}^2$. This puts its eigenvalues in $[-4\beta /h^2,0]$, so the same restriction gives L2 norm contraction. A stable recurrence need not be consistent with a different proposed differential equation.

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.

For the equation $u_t=u_x$, the second-order forward difference has Fourier symbol $[-(1-\cos \theta {)}^2+i\sin \theta (2-\cos \theta )]/h$. The nonpositive real part gives a contractive semidiscrete evolution in the discrete L2 norm. Equivalently, for the unitary lattice shift $S$, $\mathrm{Re}⟨u,D_{+,2}u⟩=-\Vert (S-I{)}^{2}u{\Vert }^2/(4h)$. A centred difference in another direction contributes only a skew-adjoint term, preserving this energy estimate.

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.

For a smooth solution with bounded fourth derivatives on the unit square, the five-point local residual is $O(h^2)$. The smallest eigenvalue of the unscaled Dirichlet matrix is $8{\sin }^2(\pi h/2)\ge 8h^2$. With $M^2$ interior points and $M\asymp h^{-1}$, the residual has unweighted Euclidean norm $O(h)$. The equation $Ae=h^2\tau$ therefore gives unweighted error $O(h)$, equivalent to second-order error in the area-weighted grid norm.

On a unit-square uniform mesh with $h=1/(M+1)$ and zero Dirichlet boundary conditions, the semidiscrete five-point Laplacian has eigenvalues $-4h^{-2}[{\sin }^2(j\pi h/2)+{\sin }^2(k\pi h/2)]$. Adding reaction $\kappa I$ gives a contractive continuous-time grid evolution exactly under the displayed threshold. The threshold is $2{\pi }^2-{\pi }^4{h}^2/6+O(h^4)$, slightly below the continuum threshold. This is an all-time non-growth criterion; finite-time mesh-uniform bounds also hold for fixed larger reaction rates.

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 the Dirichlet discrete Laplacian on a square, let $T=h^{-2}\mathrm{tridiag}(1,-2,1)$. The standard axial five-point operator is $T\otimes I+I\otimes T$, and the axial-weight $2/3$, corner-weight $1/6$ nine-point operator additionally contains $h^{2}T\otimes T/6$. Its eigenvalues are therefore strictly larger than the five-point values for the same sine modes. Both lie above the corresponding negative continuum Laplacian eigenvalues, so the five-point operator has smaller spectral error. For fixed mode numbers the leading errors are respectively ${\pi }^4{h}^2(k^4+l^4)/12$ and ${\pi }^4{h}^2(k^2+l^2{)}^2/12$.

For the square-grid nine-point finite-difference stencil for the Laplacian,If $\Delta u=0$, both displayed correction terms vanish: the mixed sixth derivatives sum to $\Delta u_{xxyy}=0$, and the pure sixth derivatives sum to the negative of that sum. Expanding two orders further yields $D_{9}u=h^6{u}_{xxxxxxxx}/3024+O(h^8)$ for a smooth harmonic function. Thus the compact harmonic scheme has sixth-order normalized consistency of a numerical method and, with the discrete maximum principle inverse bound, sixth-order nodal convergence on a square with exact smooth Dirichlet boundary data. For example $\mathrm{Re}(x+iy{)}^8$ shows that the sixth-order term need not vanish.

Let ${\Gamma }_9$ have center weight $-10/3$, axial weights $2/3$ and diagonal weights $1/6$, and let $M_h$ have center weight $2/3$ and axial weights $1/12$. Taylor expansion gives $h^{-2}{\Gamma }_{9}u=\Delta u+h^2{\Delta }^{2}u/12+O(h^4)$ and $M_{h}u=u+h^2\Delta u/12+O(h^4)$. On a smooth solution of the constant-coefficient Helmholtz equation, the second-order term becomes $h^2\Delta (\Delta u+\lambda u)/12=0$. Thus the normalized PDE residual is fourth order, while the original unscaled stencil residual is sixth order. This local cancellation alone does not prove uniform global accuracy at or near a resonance.

The Jacobi method has axial neighbour weights $1/5$ and diagonal weights $1/20$. Interior row sums are one, boundary-adjacent sums are smaller, and the positive-weight graph is connected. If an eigenvalue had modulus one, a maximal-modulus eigenvector component would force all its neighbours to have the same modulus. Propagation reaches a deficient row, a contradiction. This proves strict spectral-radius control even though the infinity norm of the iteration matrix can equal one.

Extend the interior grid vector by zero. With axial weights $2/(3h^2)$ and diagonal weights $1/(6h^2)$, the negative central coefficient equals the total neighbour weight. Expanding the edge-square sum gives the displayed quadratic form. Equality forces all axial differences to vanish, and connection to the zero exterior forces the vector to be zero. Simultaneous row and column permutations preserve symmetry and negative definiteness.

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.

The nine-point finite-difference stencil for the Laplacian has expansion $D_{9}u=\Delta u+h^2{\Delta }^{2}u/12+h^4(u_{xxxxxx}+u_{yyyyyy}+5u_{xxxxyy}+5u_{xxyyyy})/360+O(h^6)$. If $\Delta u=f$, differentiating this equation gives $u_{xxxxyy}+u_{xxyyyy}=f_{xxyy}$ and $u_{xxxxxx}+u_{yyyyyy}={\Delta }^{2}f-3f_{xxyy}$. Hence the corrected compact equationhas sixth-order normalized consistency of a numerical method. Source derivatives can be supplied analytically; numerical approximation of $\Delta f$ must be fourth order and of the fourth derivatives second order. With exact square-grid Dirichlet boundary data and uniformly smooth solution, the discrete maximum principle inverse bound converts this residual into sixth-order nodal convergence. The unknown stencil stays compact even if more distant known source values are used.

On the unit square with exact Dirichlet boundary data, let $D_{9}U=B_{h}f$ and suppose the normalized residual on the exact smooth solution has supremum norm at most $\epsilon$. The error satisfies $D_{9}e=-\tau$ and vanishes on the boundary. The barrier $q=[x(1-x)+y(1-y)]/4$ has $D_{9}q=-1$, $q\ge 0$ on the boundary, and maximum $1/8$. Apply the discrete maximum principle to $e-\epsilon q$ and $-e-\epsilon q$, obtaining $|e|\le \epsilon q$. Thus $\Vert e{\Vert }_{\infty }\le \epsilon /8$, and fourth-order residual gives fourth-order supremum norm nodal convergence. The proof needs a uniform smoothness bound and a compatible boundary discretization.

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.

A constant-coefficient translation-invariant grid update is diagonalized spatially by discrete Fourier modes. Its scalar or matrix multiplier is the displayed symbol. In a semidiscrete evolution, the analogous symbol is a continuous-time generator and its exponential propagates each mode. Discrete Parseval identity converts uniform modal bounds to discrete L2 norm bounds.

For a scalar constant-coefficient grid update, Fourier transformation turns $S$ into multiplication by its symbol $G$. Parseval identity proves the displayed norm equality; necessity is obtained by concentrating Fourier data near a maximizing frequency. Contractivity for every step is equivalent to $|G|\le 1$. Mesh-uniform fixed-time stability, on $nk\le T$, permits $\sup |G|\le 1+Ck$: powers are bounded by $e^{CT}$. Conversely a uniform bound at $n=\lfloor T/k\rfloor$ gives this $1+O(k)$ restriction as $k\to 0$. For matrix symbols, use the existing uniform power bound for matrix Fourier symbols rather than only eigenvalue moduli.

This uniform bound is the actual fixed-time stability requirement for a vector-valued periodic grid evolution. Eigenvalue moduli alone are insufficient for nonnormal matrices. A unit Jordan block with an order-one off-diagonal entry grows like $n$, while an entry proportional to $\Delta t$ can give bounded growth on $n\Delta t\le T$. Uniform eigenvector conditioning or a suitable discrete energy estimate can establish the bound. Mesh-dependent diagonalizations cannot be treated as automatically harmless.

The centered unit-step approximation to a second derivative is

The weights $1,-2,1$ are finite difference coefficients for this centered second-derivative approximation.

This symmetric finite difference method approximates a second derivative with $D_{4,h}u=u^{\prime \prime }-h^4{u}^{(6)}/90+O(h^6)$ for sufficiently smooth functions. Matching the constant, second-derivative and fourth-derivative coefficients of a symmetric five-node Taylor expansion determines the displayed weights uniquely; its sixth-derivative coefficient is nonzero. The Fourier symbol is $-4h^{-2}{\sin }^2(\theta /2)[1+{\sin }^2(\theta /2)/3]\le 0$. Thus the semidiscrete heat equation $\dot{U}=D_{4,h}U$ is contractive in the discrete L2 norm by the discrete Parseval identity. This says nothing by itself about stability of a subsequently chosen time update.

For $f\in C^3[x-h,x+h]$, the Peano kernel for the error of the central secant slope has one sign. Integrating its absolute value gives the displayed constant $h^2/6$, and a cubic polynomial attains it. At $x=1$, $h=1$, the unnormalized order-two kernel is $-{\theta }^2/2$ on $[0,1]$ and $-(2-\theta {)}^2/2$ on $[1,2]$; the integral representation includes the factor $1/2!$.

Forward time with centered advection has amplification magnitude greater than one at nonzero modes. Under a fixed nonzero Courant number, its powers blow up as the time step vanishes. If instead $\Delta t/h^2$ stays bounded, then $\log |G|\le c^2\Delta t^2/(2h^2)$, so $|G{|}^n\le \exp (c^{2}T\Delta t/(2h^2))$ is uniformly bounded on fixed time intervals. This severe parabolic refinement supplies stability without contraction; it does not validate the usual hyperbolic-step choice.

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.