Chapter 22 — Key Takeaways (Numerical Integration and Differentiation)
A one-page reference for computing derivatives and integrals from function samples — and for knowing how wrong the answer is.
The finite differences (approximating $f'$ and $f''$)
| Formula | Approximation | Leading error | Order |
|---|---|---|---|
| Forward | $\dfrac{f(x+h) - f(x)}{h}$ | $\tfrac{h}{2}f''(x)$ | $O(h)$ |
| Backward | $\dfrac{f(x) - f(x-h)}{h}$ | $-\tfrac{h}{2}f''(x)$ | $O(h)$ |
| Central | $\dfrac{f(x+h) - f(x-h)}{2h}$ | $\tfrac{h^2}{6}f'''(x)$ | $O(h^2)$ |
| Second difference (the stencil) | $\dfrac{f(x+h) - 2f(x) + f(x-h)}{h^2}$ | $\tfrac{h^2}{12}f^{(4)}(x)$ | $O(h^2)$ |
Rule to memorize: symmetry buys an order. The central difference is $O(h^2)$; the one-sided ones are $O(h)$. The $(1, -2, 1)/h^2$ second difference is the 1-D Laplacian stencil (→ Ch. 24's five-point stencil).
The quadrature rules (approximating $\int_a^b f$)
| Rule | Weights (× $h$ or × factor) | Error | Exact for | Nodes |
|---|---|---|---|---|
| Trapezoid | $\tfrac12, 1, 1, \ldots, 1, \tfrac12$ (× $h$) | $-\tfrac{(b-a)h^2}{12}f''$ | degree 1 | equally spaced |
| Simpson ($n$ even) | $1,4,2,4,\ldots,4,1$ (× $\tfrac{h}{3}$) | $-\tfrac{(b-a)h^4}{180}f^{(4)}$ | degree 3 | equally spaced |
| Gauss-$n$ | tabulated $w_i$ | — | degree $2n-1$ | roots of Legendre $P_n$ |
Gauss-2 nodes on $[-1,1]$: $\pm\tfrac{1}{\sqrt3}\approx\pm0.57735$, weights $1,1$. Gauss-3: $0,\pm\sqrt{3/5} \approx\pm0.7746$, weights $\tfrac89,\tfrac59,\tfrac59$. Map $[-1,1]\to[a,b]$ with $x=\tfrac{b-a}{2}t + \tfrac{a+b}{2}$, Jacobian $\tfrac{b-a}{2}$.
Order of accuracy — the halving test
$$ E(h) \approx C h^p \quad\Longrightarrow\quad \frac{E(h)}{E(h/2)} \to 2^p. $$
| Observed error ratio per halving | Order $p$ | It's a… |
|---|---|---|
| $\approx 2$ | 1 | forward/backward difference |
| $\approx 4$ | 2 | central difference, trapezoid, the stencil |
| $\approx 16$ | 4 | Simpson's rule |
The universal correctness check: refine, take the error ratio, confirm it is $2^p$. Wrong ratio → bug.
Error estimates and extrapolation
- Adaptive error estimate (Simpson): $\text{error}(S_{2n}) \approx \dfrac{S_{2n} - S_n}{15}$ — compare two resolutions to estimate the error without knowing the answer.
- Richardson extrapolation (order $p$): $\dfrac{2^p A(h/2) - A(h)}{2^p - 1}$ cancels the $h^p$ term. Applied to the trapezoid at halved $h$ = Romberg integration.
⚠️ The round-off floor (differentiation only)
A finite-difference derivative subtracts near-equal values → catastrophic cancellation (Ch. 20). Total error is U-shaped: truncation $\sim h^p$ falls, round-off $\sim \varepsilon/h$ rises.
| Method | Optimal step $h^*$ | Best error |
|---|---|---|
| Forward difference | $\sim \sqrt{\varepsilon} \approx 10^{-8}$ | $\sim \sqrt{\varepsilon} \approx 10^{-8}$ |
| Central difference | $\sim \varepsilon^{1/3} \approx 10^{-5}$ | $\sim \varepsilon^{2/3} \approx 10^{-11}$ |
Smaller $h$ is NOT always better. Below $h^*$ the derivative gets worse. Integration is immune (it adds, never subtracts near-equals).
Multidimensional (the note)
- Product rule: $\sum_i\sum_j w_i w_j f(x_i,y_j)$ — $O(h^p)$ per axis, but $n^d$ nodes = curse of dimensionality.
- Monte Carlo: $\int \approx \tfrac{V}{N}\sum_k f(\mathbf{x}_k)$, error $O(1/\sqrt N)$ — slow but dimension-independent; wins above $d\approx4$–$8$.
Fortran idioms introduced
! Pass an integrand as a procedure argument (needs a pure abstract interface)
abstract interface
pure function scalar_fn(x) result(y)
import :: dp
real(dp), intent(in) :: x
real(dp) :: y
end function scalar_fn
end interface
pure function simpson(f, a, b, n) result(s)
procedure(scalar_fn) :: f ! f inherits `pure` from the interface
...
end function
! Central difference as a whole-array section (interior of u(1:n)) — vectorizes
d(2:n-1) = (u(3:n) - u(1:n-2)) / (2.0_dp*h)
| Feature | What it does |
|---|---|
abstract interface + procedure(scalar_fn) :: f |
Pass any conforming function as an argument (Ch. 6). |
pure function interface |
Required so a pure integrator may call f; frees the optimizer. |
| Array-section derivative | (u(3:n) - u(1:n-2))/(2h) — whole-array central difference, no loop. |
Decision aid — which rule?
| Situation | Reach for |
|---|---|
| Derivative in a hot loop / on a grid | Central difference (array section); it's the stencil |
| Need a derivative accurately | Central difference, $h \approx \varepsilon^{1/3}\max(|x|,1)$; or analytic if available |
| Smooth analytic integrand, evals cheap | Simpson (simple) |
| Smooth analytic integrand, evals expensive | Gauss-$n$ (most accuracy per evaluation) |
| Sampled data (nodes fixed) | Trapezoidal rule — you can't move the nodes |
| Integrand spiky / uneven | Adaptive (compare-and-subdivide) |
| $> 4$–$8$ dimensions | Monte Carlo |
Numbers worth memorizing
- Central difference & trapezoid & stencil: $O(h^2)$ (ratio 4). Simpson: $O(h^4)$ (ratio 16).
- Gauss-$n$ is exact to degree $2n-1$; Simpson to degree 3.
- Forward-difference optimal step: $h^* \sim \sqrt\varepsilon \approx 10^{-8}$ (double precision).
- Product-rule cost in $d$ dimensions: $n^d$; Monte Carlo error: $1/\sqrt N$.
The heat-solver piece added
No new file — a test. heat-solver/tests/test_stencil_order.f90 verifies the spatial stencil
$(u_{i-1} - 2u_i + u_{i+1})/h^2$ is $O(h^2)$ by refining the grid and confirming the error ratio is $4$. It
is the first regression test of the suite (→ Ch. 37) and guards the accuracy of $\nabla^2 u$ before Ch. 24
builds the real PDE core.