Finite Difference Methods

Pdes

Finite difference methods discretize derivatives by replacing them with differences between nearby grid values. The simplest and oldest approach to PDEs: lay down a uniform mesh, write down a stencil that approximates each derivative to some order in the spacing, and assemble the resulting linear (or nonlinear) algebraic system. Three questions dominate the analysis: ACCURACY (how fast does error go to zero as the mesh refines?), STABILITY (do small perturbations grow without bound, or stay bounded?), and CONSISTENCY (does the discrete equation actually approximate the continuous one?). The famous Lax equivalence theorem says: for well-posed linear PDEs, CONSISTENCY + STABILITY ↔ CONVERGENCE. Get those two right and the discrete solution converges to the true one.

Discretizing derivatives

The CENTRAL second-derivative stencil:

u^{''} (x_{j}) \approx \frac{u _{j - 1} - 2 u _{j} + u _{j + 1}}{h ^{2}} .

Taylor expansion of $u_{j \pm 1}$ around $x_{j}$ :

u_{j \pm 1} = u_{j} \pm h u_{j}^{'} + \frac{h ^{2}}{2} u_{j}^{''} \pm \frac{h ^{3}}{6} u_{j}^{'''} + \frac{h ^{4}}{24} u_{j}^{(4)} + \dots

Summing gives $u_{j + 1} + u_{j - 1} - 2 u_{j} = h^{2} u_{j}^{''} + \frac{h ^{4}}{12} u_{j}^{(4)} + O (h^{6})$ , so the central second-derivative stencil is SECOND-ORDER ACCURATE: the leading error is $\frac{h ^{2}}{12} u^{(4)}$ . The first-derivative analogue $(u_{j + 1} - u_{j - 1}) / (2 h)$ is also second-order. ONE-SIDED differences $(u_{j} - u_{j - 1}) / h$ are FIRST-order — useful for upwind discretizations of advection but pay a one-order accuracy penalty.

Worked example: 1D Poisson

For the elliptic problem $- u^{''} (x) = f (x)$ on $[0, 1]$ with Dirichlet boundary conditions $u (0) = u (1) = 0$ , the central second-derivative stencil yields a tridiagonal linear system:

\frac{1}{h ^{2}} 2 - 1 - 1 2 ⋱ - 1 ⋱ - 1 ⋱ 2 u_{1} u_{2} ⋮ u_{N - 2} = f (x_{1}) f (x_{2}) ⋮ f (x_{N - 2}) .

A symmetric positive definite tridiagonal system. The condition number of this matrix scales as $h^{- 2}$ — so refining the mesh by 10× makes the linear solve $1 0^{4} \times$ harder (in flop count for an iterative method). For 1D this is irrelevant; for 2D and 3D the $h^{- 2}$ conditioning is what motivates preconditioning and multigrid. See condition number for the analysis.

Time-stepping for parabolic equations

For the heat equation $u_{t} = α u_{xx}$ , the simplest scheme is FTCS (forward in time, central in space):

u_{j}^{n + 1} = u_{j}^{n} + r (u_{j - 1}^{n} - 2 u_{j}^{n} + u_{j + 1}^{n}), r = α Δ t / h^{2} .

First-order in time, second-order in space, EXPLICIT (each new value computed in one pass, no linear solve). The von Neumann stability analysis: substitute $u_{j}^{n} = G^{n} e^{ikjh}$ (a single Fourier mode of growth factor $G$ ). Plug in:

G (k) = 1 - 4 r sin^{2} (kh /2) .

Stability requires $∣ G (k) ∣ \leq 1$ for ALL wave numbers $k$ . The most-unstable mode is $kh = π$ (the highest representable frequency), giving $G = 1 - 4 r$ . The condition $∣1 - 4 r ∣ \leq 1$ is $0 \leq r \leq 1/2$ . EXACTLY at the boundary, $r = 1/2$ gives $G = - 1$ : marginal stability. Above the boundary, the $kh = π$ mode grows exponentially.

The $r \leq 1/2$ constraint is restrictive: it forces $Δ t \leq h^{2} / (2 α)$ . Refining the mesh by 10× makes the time step shrink by 100×, so total work scales as $h^{- 3}$ in 1D — punishing at fine resolutions. The standard fix is implicit time stepping: BACKWARD EULER or CRANK-NICOLSON. Both are unconditionally stable (any $Δ t$ works) but require solving a linear system at every step. Crank-Nicolson is second-order in time as well; it's the standard default for parabolic problems.

The CFL condition for hyperbolic equations

For the linear transport equation $u_{t} + c u_{x} = 0$ with $c > 0$ , the FIRST-ORDER UPWIND scheme is

u_{j}^{n + 1} = u_{j}^{n} - σ (u_{j}^{n} - u_{j - 1}^{n}), σ = c Δ t / h (the CFL number) .

Von Neumann growth factor: $G (k) = (1 - σ) + σ e^{- ikh}$ . The magnitude squared simplifies to $∣ G ∣^{2} = 1 - 4 σ (1 - σ) sin^{2} (kh /2)$ . Stability $∣ G ∣ \leq 1$ requires $0 \leq σ \leq 1$ — the CFL CONDITION named after Courant, Friedrichs, and Lewy (1928).

Geometrically: the discrete scheme has a finite domain of dependence ( $u_{j}^{n + 1}$ depends only on $u_{j - 1}^{n}, u_{j}^{n}$ ). The CFL condition says that this discrete domain of dependence must CONTAIN the continuous one (the characteristic line back to the previous time level). If $σ > 1$ , the foot of the characteristic falls outside the upwind stencil, the scheme is asking values "from the future" relative to its information, and any numerical error in the unobserved direction blows up.

Code: three demonstrations

# Finite difference methods for PDEs — three demonstrations:
#   (a) 1D Poisson convergence at O(h²)
#   (b) FTCS heat-equation stability boundary at r = 1/2
#   (c) Upwind advection CFL boundary at CFL = 1

import numpy as np

# ─── (a) Convergence of central FD on 1D Poisson ────────────────────────
# -u''(x) = pi² sin(pi x),  u(0) = u(1) = 0,  exact u(x) = sin(pi x).
def poisson_convergence():
    print("  (a) Poisson -u'' = pi² sin(pi x) — central FD, second order:")
    print(f"  {'N':>5s}  {'h':>10s}  {'||err||_inf':>14s}  {'ratio':>8s}")
    errs = []
    for N in [11, 21, 41, 81, 161, 321]:
        x = np.linspace(0, 1, N)
        h = x[1] - x[0]
        A = (np.diag(-2*np.ones(N-2)) + np.diag(np.ones(N-3), 1)
             + np.diag(np.ones(N-3), -1)) / h**2
        rhs = -np.pi**2 * np.sin(np.pi * x[1:-1])
        u = np.zeros(N)
        u[1:-1] = np.linalg.solve(A, rhs)
        err = np.max(np.abs(u - np.sin(np.pi * x)))
        ratio = errs[-1] / err if errs else float('nan')
        errs.append(err)
        print(f"  {N:>5d}  {h:>10.4f}  {err:>14.4e}  {ratio:>8.2f}")
    print("  (h halves → error /≈4: O(h²) confirmed)")

# ─── (b) Stability boundary for FTCS heat equation ──────────────────────
# Use random initial data so high-frequency modes are present.
def heat_stability():
    print("\n  (b) FTCS heat eqn:  u^{n+1}_j = u^n_j + r (u^n_{j+1} - 2 u^n_j + u^n_{j-1})")
    print("      Stability bound r = alpha dt/dx² <= 1/2.")
    print(f"    {'r':>6s}  {'final max|u|':>16s}  {'status':>10s}")
    N = 51
    h = 1.0 / (N - 1)
    np.random.seed(1)
    u0 = np.random.randn(N); u0[0] = u0[-1] = 0.0
    for r in [0.40, 0.50, 0.51, 0.55, 0.60]:
        dt = r * h**2
        nsteps = int(0.05 / dt)
        u = u0.copy()
        for _ in range(nsteps):
            u_new = u.copy()
            u_new[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
            u_new[0] = u_new[-1] = 0.0
            u = u_new
        peak = np.max(np.abs(u))
        status = "STABLE" if peak < 100 else "BLOWUP"
        print(f"    {r:>6.2f}  {peak:>16.3e}  {status:>10s}")

# ─── (c) CFL boundary for upwind advection ──────────────────────────────
def upwind_cfl():
    print("\n  (c) Upwind advection u_t + u_x = 0:  CFL = c dt/dx must be <= 1.")
    print(f"    {'CFL':>6s}  {'final max|u|':>16s}  {'status':>10s}")
    N = 401
    h = 1.0 / N
    np.random.seed(2)
    u0 = np.random.randn(N)
    for cfl in [0.5, 0.9, 1.0, 1.1, 1.5]:
        dt = cfl * h
        nsteps = int(0.4 / dt)
        u = u0.copy()
        for _ in range(nsteps):
            u_new = u.copy()
            u_new[1:] = u[1:] - cfl * (u[1:] - u[:-1])
            u_new[0]  = u[0]  - cfl * (u[0]  - u[-1])
            u = u_new
        peak = np.max(np.abs(u))
        status = "STABLE" if peak < 100 else "BLOWUP"
        print(f"    {cfl:>6.2f}  {peak:>16.3e}  {status:>10s}")

poisson_convergence()
heat_stability()
upwind_cfl()

Output:

  (a) Poisson -u'' = pi² sin(pi x) — central FD, second order:
      N           h     ||err||_inf     ratio
     11      0.1000      8.2654e-03       nan
     21      0.0500      2.0587e-03      4.01
     41      0.0250      5.1420e-04      4.00
     81      0.0125      1.2852e-04      4.00
    161      0.0063      3.2128e-05      4.00
    321      0.0031      8.0319e-06      4.00
  (h halves → error /≈4: O(h²) confirmed)

  (b) FTCS heat eqn:  u^{n+1}_j = u^n_j + r (u^n_{j+1} - 2 u^n_j + u^n_{j-1})
      Stability bound r = alpha dt/dx² <= 1/2.
         r      final max|u|      status
      0.40         7.152e-02      STABLE
      0.50         2.166e-01      STABLE
      0.51         2.343e+03      BLOWUP
      0.55         1.731e+17      BLOWUP
      0.60         5.348e+29      BLOWUP

  (c) Upwind advection u_t + u_x = 0:  CFL = c dt/dx must be <= 1.
       CFL      final max|u|      status
      0.50         4.310e-01      STABLE
      0.90         7.524e-01      STABLE
      1.00         4.109e+00      STABLE
      1.10         2.118e+11      BLOWUP
      1.50         4.313e+31      BLOWUP

Three things land cleanly. (1) O(h²) Poisson convergence: error ratio is 4.0 every time the mesh is doubled — exactly the predicted second-order behavior. By $N = 321$ we have $\sim 1 0^{- 5}$ error from a 6th-grade-arithmetic discretization. (2) The r = 1/2 stability boundary is sharp: r=0.50 gives $O (1 0^{- 1})$ output, r=0.51 has already blown up to 2000, r=0.60 reaches $1 0^{29}$ — the high-frequency mode amplifying by $∣1 - 4 \cdot 0.6 ∣^{\sim 50}$ per step. (3) The CFL boundary is equally sharp: CFL ≤ 1 stays bounded, CFL = 1.1 already explodes to $1 0^{11}$ , CFL = 1.5 reaches $1 0^{31}$ . These are not soft transitions; the stability criterion is a hard threshold.

Higher-order schemes and trade-offs

A few rule-of-thumb extensions worth knowing.

5-point stencil for $u^{''}$ : $(- u_{j - 2} + 16 u_{j - 1} - 30 u_{j} + 16 u_{j + 1} - u_{j + 2}) / (12 h^{2})$ is fourth-order accurate. Doubles the stencil width but raises convergence from $O (h^{2})$ to $O (h^{4})$ : at fixed accuracy, the mesh can be coarser by a factor of $N_{2nd} /1$ .
Compact / Padé differences give high-order accuracy with the same stencil width as the basic scheme, by treating the derivatives implicitly: $α u_{j - 1}^{'} + u_{j}^{'} + α u_{j + 1}^{'} = (β / h) (u_{j + 1} - u_{j - 1})$ . Lele's compact schemes hit 6th and 10th order.
Crank-Nicolson: $(u^{n + 1} - u^{n}) /Δ t = (1/2) (L u^{n + 1} + L u^{n})$ . Second-order in time, unconditionally stable, second-order in space. The default time integrator for parabolic problems.
Runge-Kutta in time: discretize space first, then integrate the resulting ODE system with RK4 — fourth-order in time. See method of lines for the framework.
Higher-order upwind / WENO: for hyperbolic problems with shocks, basic upwind is too diffusive; weighted essentially-non-oscillatory schemes give shock-capturing high-order accuracy.

When FD is the wrong tool

Complex geometry: FD on a uniform grid struggles with curved boundaries. Use finite elements on an unstructured mesh, or boundary-fitted curvilinear coordinates.
High-order accuracy on smooth solutions: spectral methods have super-exponential convergence and reach machine precision with vastly fewer degrees of freedom.
Operator splitting required: when the PDE has multiple physical effects (advection + diffusion + reaction), operator splitting integrates each part with the method best suited to it.

The three classical equations — Poisson, heat, wave; the canonical FD targets.
Spectral methods — the high-accuracy alternative for smooth solutions.
Method of lines — discretize space with FD, then integrate the resulting ODE.
Condition number — the discrete Laplacian has $κ \sim h^{- 2}$ conditioning, motivating preconditioning at fine meshes.
Finite element method — the alternative discretization strategy, especially for unstructured meshes.