Finite Differences

The heat equation and the wave equation happen to have tidy exact formulas — a Fourier series, d'Alembert's travelling waves. That is the exception, not the rule. Give a PDE an irregular domain (a car's dashboard, not a straight rod), a nonlinear term (fluid drag, not plain diffusion), or a messy source term (a real weather system, not a clean sine wave), and no formula on Earth writes down the exact solution. Calculus hands you an unsolvable equation.

Finite differences sidestep the impasse by giving up on an exact formula entirely. Instead of solving for u(x,t) everywhere, chop space and time into a grid of discrete points and track only the values there. Every partial derivative in the PDE gets replaced by ordinary subtraction and division between neighbouring grid values. An unsolvable calculus problem turns into a solvable — if large — arithmetic problem, exactly the kind a computer chews through in milliseconds.

From a limit to a subtraction

Recall how the derivative is defined in the first place: u'(x) = \lim_{h \to 0} \dfrac{u(x+h) - u(x)}{h}. A finite difference is that exact same fraction — with the limit simply never taken. Pick a small but fixed grid spacing \Delta x instead of letting h \to 0, and you get a number you can actually compute:

u_x(x_i) \approx \frac{u(x_i + \Delta x) - u(x_i)}{\Delta x} = \frac{u_{i+1} - u_i}{\Delta x}.

This is the same idea that gave us the derivative — a ratio of a small change in output to a small change in input — except now the "small" is a genuine, non-zero number sitting on a grid, not an idealised limit. That is the whole trick of finite differences: keep every piece of the calculus definition except the limit.

For the heat equation u_t = \alpha u_{xx} on a grid with spacing \Delta x and time step \Delta t, write u_i^n \approx u(x_i, t_n) and approximate both derivatives this way:

u_t \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}, \qquad u_{xx} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}.

The second formula is just the first difference applied twice — a difference of differences — which is why it needs three grid points instead of two.

The explicit update rule

Substitute and solve for the one unknown at the new time level, u_i^{n+1}. Writing the mesh ratio r = \alpha\,\Delta t/\Delta x^2, the explicit scheme is

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

Each new value is computed directly from three old neighbours — no equations to solve, just arithmetic marched forward in time. The pattern of grid points it touches is the stencil: three points at the current level feeding one at the next. It is the discrete echo of the heat equation's meaning — each point edges toward the average of its neighbours.

Replace u_t by a forward time difference and u_{xx} by a centred space difference.
Update: u_i^{n+1} = u_i^n + r(u_{i+1}^n - 2u_i^n + u_{i-1}^n), with r = \alpha\Delta t/\Delta x^2.
The stencil is three old points → one new point; no linear system to solve (it is "explicit").

The stencil

The three filled points at the lower level are the known values u_{i-1}^n, u_i^n, u_{i+1}^n; the arrows feed them into the single new value u_i^{n+1} directly above. Sweep this stencil across the grid and step forward to advance the whole solution one time level.

Worked example: a five-point grid, by hand

Take a thin rod with \alpha = 1, held at zero temperature at both ends, and split it into just five grid points, i = 0, 1, 2, 3, 4, with spacing \Delta x = 0.25. Suppose the initial "bump" of heat is

u_0^0 = 0, \quad u_1^0 = 1, \quad u_2^0 = 2, \quad u_3^0 = 1, \quad u_4^0 = 0.

Choose the mesh ratio r = 0.4 (we'll see exactly why that particular number matters on the next page). Because the ends are held at zero, only the three interior points update; each uses the explicit rule with its two neighbours:

u_2^{1} = u_2^0 + r\big(u_3^0 - 2u_2^0 + u_1^0\big) = 2 + 0.4\big(1 - 4 + 1\big) = 2 + 0.4(-2) = 1.2.

The same subtraction-and-multiply, applied at every interior point, gives the whole new row in one pass:

u_1^{1} = 1 + 0.4(2 - 2 + 0) = 1, \qquad u_2^{1} = 1.2, \qquad u_3^{1} = 1 + 0.4(0 - 2 + 2) = 1.

So after a single time step the profile has gone from [0, 1, 2, 1, 0] to [0, 1, 1.2, 1, 0] — the tall centre point has sagged towards its neighbours, exactly the smoothing-out behaviour the heat equation predicts. Repeat this same arithmetic, row after row, and the whole bump gradually flattens: a PDE "solved" by nothing more exotic than multiplication and subtraction, applied over and over.

Halving \Delta x on a rod (a 1-D problem) roughly doubles the number of grid points — cheap enough. But halve \Delta x on a flat metal plate (2-D) and the grid points quadruple, since you double in both directions at once. On a 3-D block, they go up by a factor of eight.

And that's before accounting for time: as the next page shows, the explicit scheme also demands that \Delta t shrink roughly as \Delta x^2 — so halving the spatial grid quarters the allowed time step too, meaning four times as many time levels on top of the extra spatial points. Accuracy and cost are genuinely opposite ends of the same dial: turn it towards precision and the computation bill rises fast, especially in three dimensions.

A finite-difference answer always looks like a real function — a smooth-ish curve of numbers — whether or not it is actually close to the truth. The approximation is only as good as its grid resolution. A grid that is too coarse can miss real features of the solution entirely: a sharp spike gets flattened away, a fast-changing region gets smeared into nothing, simply because there weren't enough points nearby to notice it.

The dangerous part is that a too-coarse simulation rarely announces itself — it just quietly produces a plausible-looking wrong answer. And coarseness isn't even the only way a scheme can fail: choose too large a time step for a given grid spacing, and the method can fail far more dramatically, as the next page shows.

Tomorrow's weather forecast, the crash simulation that certifies a new car design, the airflow model that shapes an aircraft wing, the ripple tank predicting how a tsunami will hit a coastline — nearly every one of these is, underneath the graphics, a grid of numbers being updated by exactly this kind of subtraction-based rule, run billions of times over on a supercomputer.

None of those problems has a tidy closed-form solution — the atmosphere is nonlinear, the car is not a perfect cylinder, the coastline is not a straight line. Finite differences (and their close cousins, finite elements and finite volumes) are quite simply how the real, messy world of PDEs gets solved. The elegant formulas you meet in a calculus course are the tip of a very large iceberg whose bulk is arithmetic on a grid.