Chapter 24 — Key Takeaways (PDEs and Finite Differences)

A one-page reference for the heat/wave equations, the five-point stencil, explicit stepping, the CFL stability limit, and boundary conditions — the chapter where the solver became real.

The equations

The five-point stencil (discrete Laplacian, $\Delta x = \Delta y = h$)

$$ \nabla^2 u\big|_{i,j} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4\,u_{i,j}}{h^2} $$

• The $1/h^2$ is not optional — it turns a neighbour comparison into a real second derivative.
• Second-order accurate: error $O(h^2)$; halving $h$ quarters the error. Exact for quadratics.
• General $\Delta x \neq \Delta y$: $\dfrac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} + \dfrac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2}$.

lap(2:n-1,2:n-1) = ( u(1:n-2,2:n-1) + u(3:n,2:n-1)      &   ! N + S
                   + u(2:n-1,1:n-2) + u(2:n-1,3:n)      &   ! E + W
                   - 4.0_dp*u(2:n-1,2:n-1) ) / h**2         ! -4*centre, scaled

The explicit FTCS update (the solver's core)

$$ u^{n+1}_{i,j} = u^n_{i,j} + r\left(u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4u^n_{i,j}\right), \qquad r = \frac{\alpha\,\Delta t}{h^2} $$

• Explicit = new value from known (old) values only; no system to solve.
• Update from a snapshot — evaluate the Laplacian on $u^n$ before writing, or use two buffers. In-place overwrite silently becomes Gauss–Seidel.

The CFL / stability limit — memorize this

• Cross the limit and it detonates: the worst mode's amplitude multiplies by $|1-8r|>1$ every step (2D) — errors double and overflow to Inf/NaN.
• $\Delta t \sim h^2$: halve $h$ ⇒ quarter $\Delta t$ ⇒ $4\times$ more steps. Refining is expensive.
• Never hard-code $\Delta t$. Derive it: dt = safety * h**2 / (4*alpha) (2D), safety ≈ 0.9.
• Wave CFL (the original): $C = c\,\Delta t/h \le 1$ — a wave crosses at most one cell per step; scales as $\Delta t \sim h$.

Boundary conditions

• Interior sweep is do i = 2, n-1 (and j) — writing 1, n reads off-grid.
• Use modulo (not mod) for periodic wrap — it returns a non-negative index at both ends.

The project's heat_solver (canonical interfaces)

pure function laplacian(u, dx, dy) result(lap)   ! scaled 5-point Laplacian, interior filled
subroutine    step(field, alpha, dt)             ! FTCS step; field is a field_t (has %dx)
pure function stable_dt(alpha, dx, dy, safety)    ! CFL-safe dt from the grid

• step's signature is frozen since Chapter 6; only the body became real, and field is now a field_t so it carries field%dx for the scaling.
• The stencil is a memory-bound kernel — loop order (inner over the first index) is the whole performance story (Ch. 27).

Intrinsics / features used

Pitfalls

• Dropping the $1/h^2$ — answer changes when you refine the mesh (a broken discretisation).
• Hard-coded $\Delta t$ — stable on a coarse grid, explodes on a fine one ($r \propto 1/h^2$).
• In-place update — mixes time levels, silently changing the scheme.
• do i = 1, n at the boundary — off-grid stencil read; crash with -fcheck=all, corruption without.
• Wrong loop order (j inner) — several-fold slowdown on a column-major array.

Numbers worth memorizing

• 2D explicit heat is stable iff $r = \alpha\Delta t/h^2 \le 1/4$. This is the whole chapter in one line.
• Refine $h \to h/2$: work grows $16\times$ in 2D ($4\times$ cells $\times\ 4\times$ steps).
• The five-point stencil is $O(h^2)$; the nine-point is $O(h^4)$.