Numerical Methods for PDEs
A partial differential equation ties the value of an unknown function to its rates of
change in space and time: how heat spreads through a bar, how a drumhead vibrates, how a stock
option is priced. Almost none of the PDEs that matter in the real world have a tidy closed-form
solution. So we do what numerical analysts always do — we stop asking for a formula and start asking for
numbers on a grid.
The oldest and most transparent recipe is the finite difference method. The idea is a
single, disarmingly simple move: chop the domain into a mesh of points, and everywhere a derivative
appears in the PDE, replace it with a difference quotient — a ratio of nearby grid
values. A differential equation, which lives on a continuum, becomes a big pile of algebraic
equations relating neighbouring points. Those we can solve. This page teaches that one idea end to end:
turning derivatives into differences, and reading the two great model problems — a steady
elliptic problem and a marching parabolic one — straight off the grid.
(For the same method viewed from the PDE side, see
finite
differences for PDEs.)
From derivatives to differences
Everything rests on the definition of the derivative, with the limit
\Delta x \to 0 simply not taken. Keep a small but finite spacing
\Delta x and three natural approximations to
u'(x) fall out, differing only in which neighbours they lean on:
u'(x) \approx \underbrace{\frac{u_{i+1}-u_i}{\Delta x}}_{\text{forward}}
\quad\text{or}\quad \underbrace{\frac{u_i-u_{i-1}}{\Delta x}}_{\text{backward}}
\quad\text{or}\quad \underbrace{\frac{u_{i+1}-u_{i-1}}{2\,\Delta x}}_{\text{central}}.
Here u_i = u(x_i) is the value at grid point
x_i = i\,\Delta x. Taylor's theorem grades their accuracy. The one-sided
forward and backward quotients are first order — their error shrinks like
O(\Delta x). The central difference is the star: the
first-order error terms cancel by symmetry, leaving an error of
O(\Delta x^2). Halve the spacing and its error drops by a factor of
four, not two — a free upgrade in accuracy for the same stencil width.
For the second derivative we combine the forward and backward first differences into
the workhorse of the whole subject — the central second difference:
-
The second derivative is approximated by the three-point stencil
u''(x_i) \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}.
-
It is second-order accurate: the truncation error is
-\tfrac{1}{12}u''''(x_i)\,\Delta x^2 + O(\Delta x^4), i.e.
O(\Delta x^2).
-
The weights (1,\,-2,\,1)/\Delta x^2 are the discrete Laplacian in one
dimension — the building block of every problem below.
Let us not take the accuracy claim on faith. Take a function whose second derivative we know exactly —
u(x)=\sin x, so u''(x)=-\sin x — and watch the
stencil's error fall by four each time we halve \Delta x:
// Second-order check: (u_{i+1} - 2u_i + u_{i-1})/Δx² vs the true u''.
const u = (x: number): number => Math.sin(x);
const uSecondExact = (x: number): number => -Math.sin(x); // d²/dx² sin = -sin
const x0: number = 1.0;
let prevErr: number = 0;
for (const h of [0.2, 0.1, 0.05, 0.025]) {
const approx: number = (u(x0 + h) - 2 * u(x0) + u(x0 - h)) / (h * h);
const err: number = Math.abs(approx - uSecondExact(x0));
const ratio: string = prevErr === 0 ? "—" : (prevErr / err).toFixed(2);
console.log(`Δx=${h.toFixed(3)} approx=${approx.toFixed(6)} error=${err.toExponential(3)} ratio=${ratio}`);
prevErr = err;
}
console.log("Error ratio ≈ 4 each time Δx halves ⇒ O(Δx²).");
The ratio parks itself at almost exactly 4: numerical proof that the stencil
is second order. That single, humble three-point formula is about to build both model problems.
Laying down a grid
Before any differences make sense we need a mesh. For a one-dimensional problem on
[0, L] we scatter N+1 equally spaced points
x_i = i\,\Delta x with \Delta x = L/N. If time is
involved too, we add a second axis of steps t^n = n\,\Delta t, giving a
space–time lattice of nodes u_i^{\,n} \approx u(x_i, t^n) — subscript for
space, superscript for time. The boundary nodes carry the values the problem hands us;
the interior nodes are the unknowns we must find.
Two shapes of problem now split apart, and they demand completely different numerical treatment:
-
Elliptic — a steady state. No time at all; the solution just sits there in
equilibrium. Laplace's equation -u''=0 and Poisson's
-u''=f are the models. Every interior node is coupled to its neighbours at
once, so we get one big linear system to solve.
-
Parabolic — diffusion in time. The heat equation
u_t = \alpha\,u_{xx} smears an initial profile out as time runs. We
march: knowing every node at time level n, we compute time level
n+1, and repeat.
Elliptic: Poisson's equation as a tridiagonal system
Take the boundary-value problem -u''(x) = f(x) on
(0,1) with fixed ends u(0)=a,
u(1)=b. Substitute the central second difference at every interior node
i = 1,\dots,N-1 and multiply through by
\Delta x^2:
-u_{i-1} + 2u_i - u_{i+1} = \Delta x^2\, f_i.
Read the Laplace case f=0 and something beautiful pops out:
u_i = \tfrac12(u_{i-1} + u_{i+1}) — every interior value is the
average of its two neighbours. A discrete harmonic function has no bumps of its own; it just
interpolates smoothly between the boundaries. In two dimensions the same stencil makes each node the
average of its four neighbours, which is exactly why relaxation methods (and, much later, image blurring)
work.
Stacking one equation per interior node gives a linear system
A\mathbf{u} = \mathbf{d} whose matrix A has
2 on the diagonal and -1 on the two off-diagonals —
it is tridiagonal. That structure is a gift: a tridiagonal system of
N unknowns solves in O(N) work by the Thomas
algorithm (a streamlined Gaussian elimination), rather than the
O(N^3) of a dense solve. Let us build and solve one, choosing
f = 1 with zero boundaries so the exact answer is the parabola
u(x) = \tfrac12 x(1-x):
// Solve -u'' = 1 on (0,1), u(0)=u(1)=0. Exact: u(x) = x(1-x)/2.
const N: number = 6; // sub-intervals ⇒ N-1 interior unknowns
const dx: number = 1 / N;
const n: number = N - 1;
// Right-hand side d_i = Δx² · f_i with f = 1.
const d: number[] = [];
for (let i = 0; i < n; i++) d.push(dx * dx * 1);
// Tridiagonal (-1, 2, -1). Thomas algorithm (forward sweep + back-substitution).
const lo: number = -1, di: number = 2, up: number = -1;
const cp: number[] = new Array(n);
const dp: number[] = new Array(n);
cp[0] = up / di;
dp[0] = d[0] / di;
for (let i = 1; i < n; i++) {
const m: number = di - lo * cp[i - 1];
cp[i] = up / m;
dp[i] = (d[i] - lo * dp[i - 1]) / m;
}
const u: number[] = new Array(n);
u[n - 1] = dp[n - 1];
for (let i = n - 2; i >= 0; i--) u[i] = dp[i] - cp[i] * u[i + 1];
console.log(" x u (finite diff) exact x(1-x)/2");
for (let i = 0; i < n; i++) {
const x: number = (i + 1) * dx;
const exact: number = 0.5 * x * (1 - x);
console.log(`${x.toFixed(3)} ${u[i].toFixed(6)} ${exact.toFixed(6)}`);
}
Every computed value lands on the exact parabola. That is no fluke: the second difference is
exact for quadratics (its error involves the fourth derivative, which vanishes here), so this
particular problem is captured perfectly even on a coarse grid. For a wavier
f the answer would carry the usual O(\Delta x^2)
error and refining the mesh would sharpen it. Big meshes can also be relaxed iteratively
(Jacobi, Gauss–Seidel, SOR) by repeatedly resetting each node toward the average of its neighbours — the
iterative
methods for linear systems that shine when the grid is enormous.
Because the stencil is local. Node i's equation mentions only
u_{i-1}, u_i and
u_{i+1} — nobody further away. So row i of
A has non-zeros only in columns i-1, i, i+1, and
everything else is a structural zero. Locality of the physics becomes sparsity of the
matrix, and sparsity is what makes large PDE solves affordable. In 2-D the five-point stencil gives a
banded, block-tridiagonal matrix — still overwhelmingly zeros — which is why sparse linear algebra is
the beating heart of scientific computing.
Parabolic: marching the heat equation with FTCS
Now let time back in. The 1-D heat equation u_t = \alpha\, u_{xx} says the
temperature at a point rises or falls in proportion to how much it differs from the average of its
surroundings — hot spots bleed into cold ones. Discretise the time derivative with a
forward difference and the space derivative with the central
second difference, and you get the explicit FTCS scheme (Forward-Time,
Central-Space):
-
Update rule — each new value is an old value plus a weighted curvature of its neighbours:
u_i^{\,n+1} = u_i^{\,n} + r\left(u_{i+1}^{\,n} - 2u_i^{\,n} + u_{i-1}^{\,n}\right).
-
The dimensionless diffusion number bundles the step sizes:
r = \frac{\alpha\,\Delta t}{\Delta x^2}.
-
It is explicit: the right-hand side uses only known values from time level
n, so each new node is a direct calculation — no linear system to solve.
This looks like a free lunch — no matrix, just a sweep. But explicitness has a price, and it is the
single most important fact about this scheme. The update is only stable — meaning tiny
errors decay rather than grow — when the diffusion number is small enough:
-
The explicit FTCS scheme is stable if and only if
r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \tfrac12.
-
Equivalently, the time step is capped:
\Delta t \le \dfrac{\Delta x^2}{2\alpha}.
-
Cross the line and every rounding error is amplified each step, so the numerical solution erupts
into growing sawtooth oscillations — even though the true solution is smoothly decaying.
Watch it happen. We start the same smooth bump and march it with two diffusion numbers straddling the
limit — r = 0.5 (just legal) and r = 0.6 (just
over). We track the peak magnitude: it should melt away in the stable run and explode in the unstable
one.
// FTCS for u_t = u_xx on [0,1], u(0)=u(1)=0, smooth initial bump.
function march(r: number, steps: number): number[] {
const m: number = 11; // 11 grid points
const u: number[] = [];
for (let i = 0; i < m; i++) u.push(Math.sin(Math.PI * i / (m - 1))); // half-sine bump
for (let s = 0; s < steps; s++) {
const next: number[] = u.slice(); // boundaries stay 0
for (let i = 1; i < m - 1; i++) {
next[i] = u[i] + r * (u[i + 1] - 2 * u[i] + u[i - 1]);
}
for (let i = 0; i < m; i++) u[i] = next[i];
}
return u;
}
function peak(a: number[]): number {
let mx: number = 0;
for (const v of a) mx = Math.max(mx, Math.abs(v));
return mx;
}
console.log("start peak = 1.000");
console.log(`r = 0.50 (stable): peak after 60 steps = ${peak(march(0.50, 60)).toFixed(5)}`);
console.log(`r = 0.60 (unstable): peak after 60 steps = ${peak(march(0.60, 60)).toExponential(3)}`);
console.log("Stable run decays toward 0; unstable run blows up by many orders of magnitude.");
At r=0.5 the bump diffuses away exactly as heat should. At
r=0.6 — a mere 20\% larger time step — the answer
detonates into astronomically large, alternating values. Nothing about the arithmetic looks dangerous
line by line; the instability is a property of the whole iteration, and only the stability
analysis reveals it in advance.
The FTCS update u_i + r(u_{i+1}-2u_i+u_{i-1}) looks utterly benign — a value
plus a small correction. Yet the instant r > \tfrac12 that correction
over-shoots: each step multiplies the highest-frequency grid wobble by a factor bigger than one in
magnitude, and repeated multiplication turns a rounding speck into an explosion. The tell-tale sign is
a sawtooth that flips sign every node and grows every step.
Worse, the limit \Delta t \le \Delta x^2/(2\alpha) punishes you for
refining space. Because r \propto 1/\Delta x^2, halving
\Delta x forces you to quarter
\Delta t just to stay stable — and with twice as many nodes to update, a
finer mesh can cost roughly eight times the work per unit of simulated time. Accuracy and stability
pull in opposite directions, and stability wins the veto.
The escape hatch is to go implicit. Evaluate the space difference at the
new time level instead (backward Euler), or average old and new
(Crank–Nicolson), and the scheme becomes unconditionally stable —
any \Delta t at all. The cost: the unknowns are now coupled, so every step
requires solving a tridiagonal system — the very same elliptic machinery from the Poisson problem
above. You trade a cheap-but-fragile explicit sweep for a costlier-but-robust solve.
Seeing the diffusion
With a stable diffusion number the half-sine bump keeps its shape but shrinks in height — the heat
equation's signature: high spatial frequencies decay fastest, smoothing every profile toward its
boundary values. The chart shows the initial temperature profile and the same profile after many stable
FTCS steps.
The one big idea
Both model problems came from the same three-point stencil
(u_{i+1}-2u_i+u_{i-1})/\Delta x^2. Freeze time and it assembles a tridiagonal
linear system whose solution is the steady elliptic answer. Let time flow and it drives
an explicit march for the parabolic heat equation — fast and matrix-free, but shackled
by the stability limit r \le \tfrac12, unless you pay for an implicit solve to
buy unconditional stability. Replace derivatives with differences on a grid, mind your accuracy order,
and above all respect stability: that is the finite difference method in one breath.