Finite-Difference Method

The finite-difference method is a family of techniques in numerical analysis for approximating derivatives and solving differential equations. It replaces continuous differential operators with weighted differences between values at discrete points, usually arranged on a computational grid. The resulting equations determine numerical approximations to the unknown function at those points. Finite differences are used for both ordinary and partial differential equations, including problems involving diffusion, transport, waves, and equilibrium fields. (damtp.cam.ac.uk)

Basic construction

Let $u(x)$ be a sufficiently smooth function sampled on a uniform grid,

x_i=x_0+ih,

where $h>0$ is the grid spacing. Write $u_i$ for either the exact sample $u(x_i)$ or, when solving an equation, its numerical approximation. The simplest approximations to the first derivative are

\begin{aligned} D_+u_i&=\frac{u_{i+1}-u_i}{h} &&\text{(forward difference)},\\ D_-u_i&=\frac{u_i-u_{i-1}}{h} &&\text{(backward difference)},\\ D_0u_i&=\frac{u_{i+1}-u_{i-1}}{2h} &&\text{(central difference)}. \end{aligned}

Forward and backward differences have first-order truncation error, $O(h)$ ; the central difference has second-order error, $O(h^2)$ . Here big-O notation describes how the error scales as $h$ tends to zero, under suitable smoothness assumptions. (damtp.cam.ac.uk)

These formulas follow from Taylor expansions. For example,

u(x_i\pm h) =u(x_i)\pm hu'(x_i) +\frac{h^2}{2}u''(x_i) \pm\frac{h^3}{6}u'''(x_i)+\cdots.

Subtracting the expansions cancels the even-power terms and yields the central first difference. Adding them gives the standard second-derivative approximation,

u''(x_i) =\frac{u_{i-1}-2u_i+u_{i+1}}{h^2}+O(h^2).

Thus, symmetry can eliminate leading error terms without requiring additional function evaluations. (tsapps.nist.gov)

The set of points and coefficients used in a difference formula is called a stencil. More general formulas take the form

u^{(m)}(x_*)\approx\sum_{j=0}^{s}w_j u(x_j),

where the weights depend on the derivative order, evaluation point, and node locations. Formulas need not use equally spaced nodes. Bengt Fornberg’s 1988 algorithm provides recursive procedures for generating weights for arbitrary one-dimensional node distributions and derivative orders. (ams.org)

Discretizing boundary-value problems

For a boundary-value problem, finite differences replace the differential equation at interior grid points, while separate equations impose the boundary data. Consider

-u''(x)=f(x),\qquad 0<x<L,

with prescribed endpoint values. On a uniform grid, the interior equations become

\frac{-u_{i-1}+2u_i-u_{i+1}}{h^2}=f(x_i).

After incorporating the endpoint values into the right-hand side, these equations form a linear system $A\mathbf u=\mathbf b$ . Each row couples only neighboring unknowns, so $A$ is a tridiagonal sparse matrix. (ocw.mit.edu)

In two dimensions, the standard approximation to the Laplacian on a square grid is

\Delta_hu_{i,j} =\frac{ u_{i+1,j}+u_{i-1,j} +u_{i,j+1}+u_{i,j-1} -4u_{i,j}}{h^2}.

This five-point stencil is second-order accurate for smooth functions. Applied to Poisson’s equation, $-\Delta u=f$ , it produces a sparse system whose coefficients reflect the grid connectivity and boundary conditions. (dlmf.nist.gov)

Boundary treatment is part of the numerical method, not merely an implementation detail. Dirichlet conditions prescribe function values; Neumann conditions prescribe normal derivatives; and Robin conditions prescribe combinations of the two. One-sided differences or auxiliary points outside the physical domain, called ghost points, can supply equations near boundaries. Curved boundaries and interfaces may require modified stencils or extensions of the solution; these constructions must preserve the intended accuracy and stability. (arxiv.org)

Time-dependent equations

A time-dependent problem can be discretized in space and time simultaneously, or first in space alone. The latter approach, known as the method of lines, turns a partial differential equation into a system of ordinary differential equations, to which a time-integration method is applied. (damtp.cam.ac.uk)

For the one-dimensional heat equation,

u_t=\alpha u_{xx},\qquad \alpha>0,

let $u_i^n$ approximate $u(x_i,t_n)$ , with $t_n=n\Delta t$ , and define

r=\frac{\alpha\Delta t}{h^2}.

A forward time difference and centered spatial difference give

u_i^{n+1} =u_i^n+r\left(u_{i-1}^n-2u_i^n+u_{i+1}^n\right).

This explicit scheme computes the next time level directly from known values. For the standard uniform-grid problem, it is stable when $0\le r\le\tfrac12$ . Its truncation error is $O(\Delta t+h^2)$ . Consequently, refining the spatial grid requires a time step proportional to $h^2$ or smaller. (damtp.cam.ac.uk)

A backward time difference instead produces the implicit scheme

u_i^{n+1} -r\left(u_{i-1}^{n+1}-2u_i^{n+1}+u_{i+1}^{n+1}\right) =u_i^n.

It requires solving a coupled system at each time step. For this linear diffusion problem, backward Euler is unconditionally stable: stability does not impose an upper bound on $r$ . Accuracy still restricts useful time-step sizes. (damtp.cam.ac.uk)

The Crank–Nicolson method averages the spatial operator at the old and new time levels:

u_i^{n+1}-u_i^n =\frac r2\left[ \delta^2u_i^{n+1}+\delta^2u_i^n \right], \qquad \delta^2u_i=u_{i-1}-2u_i+u_{i+1}.

For smooth solutions, it is second-order accurate in both time and space and unconditionally stable for the standard linear heat equation. Nevertheless, stability does not ensure strong damping: its amplification factor approaches $-1$ for sufficiently stiff spatial modes, so large time steps can produce slowly damped, alternating numerical components. (math.mit.edu)

Consistency, stability, and convergence

Three properties organize the analysis of a difference scheme:

Consistency: substituting a smooth exact solution into the discrete equations leaves a residual that tends to zero under refinement.
Stability: perturbations in the discrete data or intermediate computation remain controlled, with bounds uniform over the relevant refinements.
Convergence: the numerical solution approaches the exact solution, in a specified norm, as the discretization is refined.

The Lax–Richtmyer equivalence theorem connects these properties: for a consistent linear approximation to a well-posed linear initial-value problem, stability is equivalent to convergence, under the theorem’s assumptions. It is not a blanket guarantee for nonlinear equations or arbitrary boundary treatments. (epubs.siam.org)

Von Neumann analysis studies how a scheme amplifies individual Fourier modes. For the explicit heat scheme, substitution of

u_i^n=G^n e^{\mathrm{i}i\theta}

gives

G(\theta)=1-4r\sin^2(\theta/2).

Requiring $|G(\theta)|\le1$ for every frequency yields $r\le\tfrac12$ . This analysis is especially useful for constant-coefficient problems on infinite or periodic grids; physical boundaries require additional attention. (math.mit.edu)

For propagation problems, the Courant–Friedrichs–Lewy condition relates the numerical domain of dependence to that of the differential equation. For example, the standard centered scheme for $u_{tt}=c^2u_{xx}$ has the stability restriction

\frac{|c|\Delta t}{h}\le1.

A CFL restriction is scheme-dependent, and satisfying a domain-of-dependence condition alone does not establish every other requirement for stability or convergence. (damtp.cam.ac.uk)

Accuracy and numerical limitations

The formal order of a stencil concerns smooth functions and local approximation error. The observed error in a computed solution also depends on boundary discretization, regularity, time integration, and the accuracy of the algebraic solver. A high-order interior stencil therefore does not, by itself, establish high-order convergence of the complete method. (epubs.siam.org)

Refinement also encounters limits from floating-point arithmetic. Numerical differentiation subtracts nearby function values and divides by a small spacing, which can magnify cancellation and evaluation errors. For a first-derivative formula of order $p$ , a simplified error model is

E(h)\approx C_{\mathrm{tr}}h^p+ C_{\mathrm{round}}\frac{\varepsilon}{h},

where $\varepsilon$ represents a floating-point error scale. The model explains why decreasing $h$ indefinitely need not improve the derivative estimate. (tsapps.nist.gov)

Transport equations introduce additional difficulties. Upwind differences incorporate the direction of propagation but can introduce artificial diffusion; other schemes can introduce phase errors or oscillations. When solutions contain shocks or steep fronts, smooth-solution error estimates are insufficient, and specialized flux formulations, limiters, or nonlinear reconstructions may be required. (damtp.cam.ac.uk)

Applications and relation to other methods

Finite differences are used to approximate equations describing diffusion, wave propagation, fluid flow, and electromagnetic fields. Their local stencils are particularly convenient on structured grids. In image processing, differences approximate image gradients, while discrete Poisson equations reconstruct images from prescribed or modified gradient fields. (math.mit.edu)

Finite differences are distinguished from the finite-element method, which typically approximates a weak or variational formulation using functions defined over mesh elements, and the finite-volume method, which balances fluxes over control volumes. These distinctions concern the derivation and interpretation of the discretization: on simple grids, different approaches can produce identical or closely related algebraic equations. Finite-volume formulations make local conservation explicit, whereas a finite-difference scheme must be constructed appropriately to retain that property. (math.mit.edu)

Historical development

A major foundation of finite-difference theory was the 1928 paper by Richard Courant, Kurt Friedrichs, and Hans Lewy on partial difference equations in mathematical physics. It investigated discrete approximations to differential equations and established the domain-of-dependence principle associated with the CFL condition. Later developments linked consistency, stability, and convergence through the Lax–Richtmyer framework and extended difference formulas to higher orders and nonuniform grids; Fornberg’s 1988 work provided a systematic algorithm for generating differentiation weights. (web.stanford.edu)

References

Course 18.336: Numerical Methods for Partial Differential Equationsmath.mit.edu
DLMF: §3.4 Differentiationdlmf.nist.gov
Numerical Differentiationtsapps.nist.gov
Generation of Finite Difference Formulas on Arbitrarily Spaced Gridsams.org
Finite Differences and Fast Poisson Solversocw.mit.edu
A Finite Difference Ghost-Cell Multigrid Approach for Poisson Equation with Mixed Boundary Conditions in Arbitrary Domainarxiv.org
High-Order Finite-Difference Discretization for Elliptic Problems on Complex Domainsmath.mit.edu
The Heat Equation and Convection-Diffusionmath.mit.edu