Finite differences

Deep dive

How it works

In: grid spacing Δx, time step Δt, initial and boundary values, scheme (stencil + time integrator) · Out: the field on the grid at each time level (or the solution of the sparse system) · State: the field at the previous time level(s) · Cost: O(N) per explicit step; stability caps Δt (heat: Δt ≤ Δx²/(2α))

1. Grid: x_i = iΔx, t_n = nΔt; boundary values from Dirichlet/Neumann conditions (ghost cells for Neumann: u_{−1} = u_1 for zero flux).
2. Stencils from Taylor series: forward (u_{i+1} − u_i)/Δx is O(Δx); central (u_{i+1} − u_{i−1})/(2Δx) is O(Δx²); second derivative (u_{i−1} − 2u_i + u_{i+1})/Δx² is O(Δx²).
3. Time: explicit Euler (FTCS), RK, or implicit (backward Euler, Crank-Nicolson) on du/dt = L·u.
4. Stability (von Neumann): insert u_j = gⁿ·e^{ikjΔx}; require |g| ≤ 1 for all k. FTCS heat: g = 1 − 4r·sin²(kΔx/2) ⇒ r ≤ 1/2. Upwind advection: |g| ≤ 1 ⇔ ν ≤ 1.
5. Consistency + stability ⇒ convergence (Lax equivalence theorem, linear problems).

heat FTCS (1D):           u_i^{n+1} = u_i + r (u_{i-1} - 2u_i + u_{i+1})
                                 ●
                                 │
                          ●──────●──────●      level n
5-point Laplacian (2D):      [  ·  1  ·  ]
                             [  1 -4  1  ] / h²
                             [  ·  1  ·  ]

Code

/// FTCS step for u_t = alpha u_xx with fixed (Dirichlet) end values. r = alpha dt / dx^2 <= 0.5.
pub fn heat_ftcs(u: &[f64], out: &mut [f64], r: f64) {
    assert!(r <= 0.5, "unstable: r = {r} > 1/2");
    let n = u.len();
    out[0] = u[0];
    out[n - 1] = u[n - 1];
    for i in 1..n - 1 { out[i] = u[i] + r * (u[i - 1] - 2.0 * u[i] + u[i + 1]); }
}

/// Thomas algorithm for tridiagonal systems: sub a, diag b, super c, rhs d (overwritten with x).
/// Stable without pivoting for diagonally dominant matrices.
pub fn thomas(a: &[f64], b: &[f64], c: &[f64], d: &mut [f64], cp: &mut [f64]) {
    let n = d.len();
    cp[0] = c[0] / b[0];
    d[0] /= b[0];
    for i in 1..n {
        let m = b[i] - a[i] * cp[i - 1];
        cp[i] = if i + 1 < n { c[i] / m } else { 0.0 };
        d[i] = (d[i] - a[i] * d[i - 1]) / m;
    }
    for i in (0..n - 1).rev() { d[i] -= cp[i] * d[i + 1]; }
}

/// Crank-Nicolson for the heat equation, Dirichlet ends fixed: (I - r/2 L) u' = (I + r/2 L) u.
/// Any r is stable; accuracy O(dt^2 + dx^2).
pub struct CrankNicolson { a: Vec<f64>, b: Vec<f64>, c: Vec<f64>, rhs: Vec<f64>, cp: Vec<f64> }

impl CrankNicolson {
    pub fn new(n: usize, r: f64) -> Self {
        let m = n - 2; // interior unknowns
        CrankNicolson { a: vec![-r / 2.0; m], b: vec![1.0 + r; m], c: vec![-r / 2.0; m], rhs: vec![0.0; m], cp: vec![0.0; m] }
    }
    pub fn step(&mut self, u: &mut [f64]) {
        let (n, h) = (u.len(), -self.a[0]); // h = r/2
        for i in 1..n - 1 {
            self.rhs[i - 1] = h * u[i - 1] + (1.0 - 2.0 * h) * u[i] + h * u[i + 1];
        }
        self.rhs[0] += h * u[0];      // boundary values move to the right-hand side
        self.rhs[n - 3] += h * u[n - 1];
        thomas(&self.a, &self.b, &self.c, &mut self.rhs, &mut self.cp);
        u[1..n - 1].copy_from_slice(&self.rhs);
    }
}

/// First-order upwind for u_t + c u_x = 0 on a periodic grid; nu = c dt/dx, stable iff |nu| <= 1.
pub fn advect_upwind(u: &[f64], out: &mut [f64], nu: f64) {
    let n = u.len();
    for i in 0..n {
        let (l, r) = (u[(i + n - 1) % n], u[(i + 1) % n]);
        out[i] = if nu >= 0.0 { u[i] - nu * (u[i] - l) } else { u[i] - nu * (r - u[i]) };
    }
}

/// Lax-Wendroff: second order in space and time, |nu| <= 1, dispersive near jumps.
pub fn advect_lax_wendroff(u: &[f64], out: &mut [f64], nu: f64) {
    let n = u.len();
    for i in 0..n {
        let (l, r) = (u[(i + n - 1) % n], u[(i + 1) % n]);
        out[i] = u[i] - 0.5 * nu * (r - l) + 0.5 * nu * nu * (r - 2.0 * u[i] + l);
    }
}

/// 5-point Laplacian on an nx-by-ny row-major grid with spacing h; boundary rows/cols untouched.
pub fn laplacian_5pt(u: &[f64], out: &mut [f64], nx: usize, ny: usize, h: f64) {
    let inv = 1.0 / (h * h);
    for j in 1..ny - 1 {
        let (up, row, dn) = (&u[(j - 1) * nx..j * nx], &u[j * nx..(j + 1) * nx], &u[(j + 1) * nx..(j + 2) * nx]);
        for i in 1..nx - 1 { // contiguous rows: the compiler vectorizes this loop
            out[j * nx + i] = (up[i] + dn[i] + row[i - 1] + row[i + 1] - 4.0 * row[i]) * inv;
        }
    }
}

/// Central-difference derivative with the step that balances truncation (h^2) against rounding (eps/h).
pub fn derivative(f: impl Fn(f64) -> f64, x: f64) -> f64 {
    let h = f64::EPSILON.cbrt() * x.abs().max(1.0); // ~6e-6 relative
    let (xp, xm) = (x + h, x - h);
    (f(xp) - f(xm)) / (xp - xm) // divide by the representable step
}

Complexity and constants

• Explicit stencil step: O(N) work, memory-bound. A 5-point Laplacian with good cache reuse moves about 16 B per point (one f64 read, one write), so order of 5-10·10⁹ points/s per socket at 100 GB/s.
• Stability limits cost: for heat, halving Δx quarters Δt, so explicit cost grows as Δx⁻³ in 1D, Δx⁻⁴ in 2D for a fixed final time. That is the point where implicit schemes (and multigrid or CG for the solves) win.
• Thomas algorithm: 8N flops, two sweeps.
• Numerical differentiation in f64: forward difference best h ≈ √ε ≈ 1.5·10⁻⁸ (error ≈ 10⁻⁸); central best h ≈ ε^(1/3) ≈ 6·10⁻⁶ (error order of 10⁻¹¹; 2.9·10⁻¹² for sin at x = 1), both relative to |x|.
• Measured with the code above (heat, sin(πx) on [0, 1], 51 points): at t = 0.1, FTCS with r = 0.5 has max error ≈ 2.4·10⁻⁴ and Crank-Nicolson with r = 5 (10× larger Δt) ≈ 1.1·10⁻⁴. FTCS at r = 0.51 amplifies a 10⁻¹² perturbation by 1.04× per step: error ≈ 800 at t = 0.2, ≈ 10¹¹ at t = 0.3.

Efficiency

• Cost: an explicit stencil step is O(N) work and memory-bound: about 16 B per point (one f64 read, one write) with good cache reuse, so order of 5-10·10⁹ points/s per socket at 100 GB/s. Stability drives total cost: for heat, explicit work grows as Δx⁻³ (1D) or Δx⁻⁴ (2D) at fixed final time. Thomas solve: 8N flops per implicit 1D step. Bound: memory bandwidth.
• Levers:
  • implicit Crank-Nicolson — r = 5 (10× larger Δt) gives max error 1.1·10⁻⁴ vs 2.4·10⁻⁴ for FTCS at r = 0.5 (measured above) — diffusion, when Δt is stability-limited
  • cache blocking and temporal blocking — several time steps per tile, fewer memory passes — 2D/3D stencils
  • higher-order or compact stencils — same accuracy on a coarser grid (Fornberg 1988) — smooth solutions
  • multigrid-preconditioned CG for the implicit solves — near-O(N) per step — implicit 2D/3D
  • generated stencil code — blocked, vectorized loops from a symbolic spec (Luporini et al. 2020) — seismic and wave workloads
• Naive → standard → state of the art:
  • explicit FTCS at the stability limit · Δt ∝ Δx², cost ∝ Δx⁻⁴ in 2D
  • Crank-Nicolson, tridiagonal or CG solve · Δt set by accuracy, O(N) per step in 1D
  • blocked, vectorized (or GPU) stencils, multigrid-preconditioned implicit solves · near memory-bandwidth speed

Implementation notes

• Upwind direction follows the sign of the velocity; for systems, upwind per characteristic (Riemann solvers in finite-volume codes).
• Numerical diffusion: first-order upwind behaves like adding diffusion c·Δx·(1 − ν)/2; it smears fronts. Higher order needs limiters (TVD, WENO) to avoid oscillations.
• Ghost cells / halos: pad arrays so the inner loop has no boundary branches; in parallel codes exchange halos between domains.
• Blocking: 2D/3D stencils benefit from cache blocking and temporal blocking (several time steps per tile).
• Implicit 2D/3D: the system is a sparse, symmetric (for diffusion) matrix with 5/7 nonzeros per row: CSR + CG preconditioned with multigrid, or ADI splitting into tridiagonal solves.
• Staggered grids (Arakawa C, Yee): store velocities/fields at faces or edges; avoids checkerboard modes in fluid and Maxwell solvers.
• Check order empirically: halve Δx (and Δt per the scheme) and confirm the error ratio (4 for second order).

Pitfalls

• Violating CFL slightly (r = 0.51, ν = 1.01): the error grows exponentially but looks fine for a while.
• Crank-Nicolson with large r on discontinuous data: decaying but sign-flipping oscillations (it is not L-stable); start with a few backward-Euler steps (Rannacher smoothing).
• One-sided differences at boundaries of lower order than the interior: the boundary dominates the global error.
• Numerical derivative with h = 10⁻¹² "for accuracy": rounding error ε/h dominates, you get garbage.
• Centred differences for advection-dominated flows (high Péclet number): spurious oscillations.

Variants and alternatives

In production

Electromagnetics: FDTD solvers (Meep). Atmosphere: WRF uses finite differences on an Arakawa C-grid. Finance: PDE engines solve Black-Scholes with Crank-Nicolson (Duffy). Geophysics: seismic wave propagation codes, many generated with Devito's stencil compiler. Imaging: diffusion filters. HPC: PETSc DMDA structured-grid applications.

Implementations

Sources

• Randall J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007. book
• John C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed., SIAM, 2004. book
• Richard Courant, Kurt Friedrichs, Hans Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen 100, 1928. paper
• John Crank, Phyllis Nicolson, A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-Conduction Type, Proc. Cambridge Philosophical Society 43, 1947. paper
• Bengt Fornberg, Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, Mathematics of Computation 51, 1988. paper
• Ardavan Oskooi et al., Meep: A flexible free-software package for electromagnetic simulations by the FDTD method, Computer Physics Communications 181, 2010. paper · code
• William C. Skamarock et al., A Description of the Advanced Research WRF Model Version 4, NCAR Technical Note, 2019. report
• Daniel J. Duffy, Finite Difference Methods in Financial Engineering, Wiley, 2006. book
• Fabio Luporini et al., Architecture and Performance of Devito, a System for Automated Stencil Computation, ACM TOMS 46(1), 2020. paper · code