The Finite Difference Method

Finite differences replace derivatives with difference quotients on a grid, turning a differential equation into a system of algebraic equations.

Derivatives on a grid

The finite difference method covers the domain with a grid and replaces each derivative in a PDE by a difference quotient of neighboring grid values. A second derivative becomes (u_{i+1} - 2u_i + u_{i-1})/h^2. Substituting these stencils everywhere converts the differential equation into a large system of algebraic equations for the grid values.

A worked example

For the one-dimensional Poisson equation -u'' = f on [0,1] with zero boundary values, the standard three-point stencil gives a tridiagonal linear system A u = h^2 f, where A has 2 on the diagonal and -1 on the off-diagonals. Solving it yields the discrete solution.

import numpy as np
def poisson_1d(f, n):
    h = 1.0/(n+1)
    x = np.linspace(h, 1-h, n)
    A = (np.diag(2*np.ones(n)) - np.diag(np.ones(n-1),1)
         - np.diag(np.ones(n-1),-1))
    b = h**2 * f(x)
    return x, np.linalg.solve(A, b)

Accuracy and convergence

The order of a scheme is the power of h in its truncation error. Central differences are second order; wider stencils give higher order at the cost of larger bandwidth and trickier boundaries. Convergence requires both consistency (the stencil approximates the derivative) and stability, tied together by the Lax equivalence theorem.

Strengths and limits

Finite differences are simple to derive and implement on regular grids and are the natural choice for structured domains. Their weakness is complex geometry: irregular boundaries are awkward compared with finite elements. On rectangular or mapped grids, though, they are efficient and easy to reason about.

Structured finite-difference grids are used for many field and transport solves in the breeder Hyperion models where the geometry maps cleanly to a regular mesh.