Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
Parabolic Equations: Heat and Diffusion
65m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

Parabolic Equations: Heat and Diffusion

The heat equation, explicit FTCS and its stability limit, implicit and Crank–Nicolson schemes, von Neumann analysis, consolidation

Manchester, 1947
John Crank and Phyllis Nicolson, mathematicians at the University of Manchester, are computing temperatures in a material that generates heat as it reacts — a problem from industrial chemistry. The governing equation is the heat equation, and the obvious numerical approach, stepping forward in time with differences in space, works beautifully — until the time step is made a little too large, at which point the computed temperatures develop wild oscillations and grow without limit. The stable steps are so small that a hand computation would take forever. Crank and Nicolson publish A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type: average the explicit and implicit difference formulas, centring the scheme in time. Their method is stable for any time step and second-order accurate. It became the default method for diffusion problems in engineering — heat, mass transfer, groundwater, consolidation, option pricing. Later users discovered its one weakness, which this lesson will show: with large steps and sharp initial data, it rings.
The heat equation ut=αuxxu_t = \alpha u_{xx} describes diffusion: of heat through a wall, of a chemical through a membrane, of excess pore pressure through a clay layer. It is parabolic: information spreads instantly but is smoothed out — sharp features decay, and the fast (short-wavelength) components decay fastest. Discretise space with central differences on a grid of spacing Δx\Delta x, and the PDE becomes a large system of ODEs, u′=αD2u\mathbf u' = \alpha D_2\mathbf u — the method of lines. The eigenvalues of αD2\alpha D_2 range from about −απ2-\alpha\pi^2 to −4α/Δx2-4\alpha/\Delta x^2: the semi-discrete system is stiff (lesson 9-4), and stiffer the finer the grid.
Derivation
The heat equation from a slab of material, and the eigenvalues of its grid
The equation. Take a slice of a bar of length dxdx, cross-section AA, density ρ\rho and specific heat cc. Its stored heat rises at the rate at which heat flows in through one face and out through the other; Fourier's law gives each flow as −kA ∂u/∂x-kA\,\partial u/\partial x:
ρcA dx ∂u∂t=kA ∂u∂x∣x+dx−kA ∂u∂x∣x⟹ut=α uxx,α=kρc.\rho cA\,dx\,\frac{\partial u}{\partial t} = kA\,\frac{\partial u}{\partial x}\Big|_{x + dx} - kA\,\frac{\partial u}{\partial x}\Big|_{x} \quad\Longrightarrow\quad u_t = \alpha\,u_{xx}, \qquad \alpha = \frac{k}{\rho c}.
The same balance of storage against flux differences gives diffusion of a chemical (Fick's law) and Terzaghi's pore pressure below.
The eigenvalues. D2D_2 is lesson 3-3's second-difference matrix, −1Δx2tridiag⁡(−1,2,−1)-\tfrac1{\Delta x^2}\operatorname{tridiag}(-1, 2, -1), whose eigenvalues lesson 4-1 found exactly:
αλj=−4αΔx2sin⁡2jπΔx2,j=1,…,n.\alpha\lambda_j = -\frac{4\alpha}{\Delta x^2}\sin^2\frac{j\pi\Delta x}{2}, \qquad j = 1, \dots, n.
The slowest, j=1j = 1, is about −απ2-\alpha\pi^2 (on a unit interval), whatever the grid; the fastest is almost −4α/Δx2-4\alpha/\Delta x^2. Forward Euler is stable only if every Δt αλj\Delta t\,\alpha\lambda_j lies in its disc, which for the fastest mode means Δt⋅4α/Δx2≤2\Delta t \cdot 4\alpha/\Delta x^2 \le 2 — exactly the limit r≤12r \le \tfrac12 that the Fourier analysis below finds.
💡Explicit FTCS and its stability limit
Forward Euler in time, central differences in space (FTCS): uin+1=uin+r (ui+1n−2uin+ui−1n),r=α ΔtΔx2.u_i^{n+1} = u_i^n + r\,(u_{i+1}^n - 2u_i^n + u_{i-1}^n), \qquad r = \dfrac{\alpha\,\Delta t}{\Delta x^2}. Von Neumann analysis substitutes a Fourier mode uin=GneikiΔxu_i^n = G^ne^{ik i\Delta x} and finds the amplification factor G=1−4rsin⁡2(kΔx/2)G = 1 - 4r\sin^2(k\Delta x/2). For the most oscillatory mode G=1−4rG = 1 - 4r, so stability needs r≤12r \le \tfrac12: Δt≤Δx2/(2α)\Delta t \le \Delta x^2/(2\alpha). Halve the grid spacing and the time step must shrink fourfold. On a 50-cell slab, r=0.50r = 0.50 is fine and r=0.52r = 0.52 gives an error of 6×1036 \times 10^3 — a rounding-sized disturbance in the highest mode grows by ∣1−4r∣=1.08|1 - 4r| = 1.08 per step.
Derivation
Von Neumann's calculation, written out
Put a single Fourier mode, ujn=Gneiφju_j^n = G^n e^{\mathrm i\varphi j} with φ=kΔx\varphi = k\Delta x and i=−1\mathrm i = \sqrt{-1}, into FTCS (written for grid point jj). The neighbours are the same mode shifted, uj±1n=ujne±iφu_{j\pm1}^n = u_j^n e^{\pm\mathrm i\varphi}, so every term carries the common factor ujnu_j^n:
G=1+r(eiφ−2+e−iφ)=1+r (2cos⁡φ−2)=1−4rsin⁡2φ2,G = 1 + r\big(e^{\mathrm i\varphi} - 2 + e^{-\mathrm i\varphi}\big) = 1 + r\,(2\cos\varphi - 2) = 1 - 4r\sin^2\frac\varphi2,
using eiφ+e−iφ=2cos⁡φe^{\mathrm i\varphi} + e^{-\mathrm i\varphi} = 2\cos\varphi and 1−cos⁡φ=2sin⁡2(φ/2)1 - \cos\varphi = 2\sin^2(\varphi/2).
GG is real and never above 1. It is at least −1-1 for every mode exactly when 1−4r≥−11 - 4r \ge -1, i.e. r≤12r \le \tfrac12. The most dangerous mode is the sawtooth, φ=π\varphi = \pi, alternating sign from point to point.
The θ\theta-scheme. The same substitution with the implicit part on the left gives G(1+4θrs2)=1−4(1−θ)rs2G\big(1 + 4\theta r s^2\big) = 1 - 4(1 - \theta)r s^2, s=sin⁡(kΔx/2)s = \sin(k\Delta x/2), the formula in the box below. For θ≥12\theta \ge \tfrac12 the numerator can never exceed the denominator in size, so ∣G∣≤1|G| \le 1 for every rr. For Crank–Nicolson and large rr, G→(1−2rs2)/(1+2rs2)→−1G \to (1 - 2rs^2)/(1 + 2rs^2) \to -1: the ringing in the warning below.
uin+1−uinΔt=α[θ δ2uin+1Δx2+(1−θ)δ2uinΔx2],G=1−4(1−θ)rsin⁡2(kΔx/2)1+4θrsin⁡2(kΔx/2)\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha\left[\theta\,\frac{\delta^2u^{n+1}_i}{\Delta x^2} + (1-\theta)\frac{\delta^2u^n_i}{\Delta x^2}\right], \qquad G = \frac{1 - 4(1-\theta)r\sin^2(k\Delta x/2)}{1 + 4\theta r\sin^2(k\Delta x/2)}

The θ-scheme: θ = 0 explicit FTCS, θ = 1 backward Euler (fully implicit), θ = ½ Crank–Nicolson. For θ ≥ ½ |G| ≤ 1 for every r: unconditionally stable. Each implicit step solves a tridiagonal system — O(n) work.

⚠Crank–Nicolson rings
For short wavelengths and large rr, Crank–Nicolson's G→−1G \to -1 (lesson 9-4's trapezoidal rule exactly): sharp features do not decay, they flip sign every step. A slab at 100 °C whose faces are suddenly cooled to 0 °C, stepped with r=25r = 25: the exact temperature next to the face after t=0.02t = 0.02 is 8 °C; backward Euler gives 9.9 °C; Crank–Nicolson gives 57.7 °C — a spurious oscillation, error 50 °C. With r=5r = 5 the problem disappears. Remedies: smaller steps at the start, a few backward-Euler steps first (Rannacher smoothing), or a fully implicit scheme.
Heat, like gravity, penetrates every substance of the universe, its rays occupy all parts of space.
– Joseph Fourier, Théorie analytique de la chaleur (1822)
📐Terzaghi consolidation
When a building is placed on saturated clay, the load is first carried by the pore water as excess pore pressure uu. Water slowly squeezes out through the clay's drained faces, the load transfers to the soil skeleton, and the ground settles. Karl Terzaghi showed in 1925 that uu obeys the heat equation, ut=cvuzzu_t = c_vu_{zz}, with the coefficient of consolidation cvc_v. For a 10 m clay layer drained top and bottom with cv=2c_v = 2 m²/year, a backward-Euler solution gives the degree of consolidation U=0.504U = 0.504 at time factor Tv=cvt/H2=0.2T_v = c_vt/H^2 = 0.2 (2.5 years) and U=0.900U = 0.900 at Tv=0.848T_v = 0.848 — 10.6 years for 90% of the settlement — matching Terzaghi's series to three digits. Settlement predictions for embankments and foundations on clay are exactly this calculation.
💡Where you will use this
Geotechnical Engineering: consolidation settlement, transient groundwater flow. Thermal design: transient heating of turbine blades, rocket nozzle walls, re-entry heat shields, concrete curing. Rocket Propulsion: ablative liners and propellant-grain heating. Chemical and environmental: diffusion of contaminants. In 2-D and 3-D, implicit steps need sparse solvers (Units 3–4) — or alternating-direction implicit (ADI) methods that keep each solve tridiagonal.
📚Further reading
  • Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 9, 'Diffusion Equations and Parabolic Problems' — the method of lines, von Neumann analysis and Crank–Nicolson.
  • K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations (2nd ed., Cambridge University Press, 2005) — parabolic problems first, with careful stability theory.
  • John Crank and Phyllis Nicolson, 'A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type', Proceedings of the Cambridge Philosophical Society 43(1), 50–67 (1947).
↓But there is a problem
Diffusion smooths everything; errors decay along with the solution, and implicit methods are stable at any step. Remove the diffusion, and none of that holds. In 1922 Lewis Fry Richardson published a book describing how weather could be forecast by computing the equations of the atmosphere. To show it worked, he had computed one six-hour forecast for a single point over Europe — by hand, during his service in an ambulance unit in the First World War, over many weeks. His result: a pressure rise of 145 hPa in six hours. The barometer had barely moved.