Chapter 22: Numerical Integration and Differentiation

"Essentially, all models are wrong, but some are useful." — George E. P. Box, statistician

Overview

The calculus you learned in school is exact and symbolic: the derivative of $\sin x$ is $\cos x$, the integral of $x^2$ is $x^3/3$, and both come from rules you apply on paper. The calculus a computer does is neither exact nor symbolic. Your machine cannot take a limit as $h \to 0$, and it usually does not even have a formula for the function it is differentiating — it has a subroutine that returns $f(x)$ for any $x$ you ask, and nothing more. Numerical differentiation and integration are the art of recovering derivatives and integrals from that one humble capability: sampling the function at a handful of points and combining the samples with cleverly chosen weights. Every physical simulation does this constantly. The heat solver you are building computes a second derivative — the Laplacian $\nabla^2 u$ — at every grid point of the plate, on every timestep, and it does so with exactly the finite-difference formulas of this chapter.

This is the second stop in Part V's numerical tour. Chapter 20 taught you what your numbers are and how rounding limits them; Chapter 21 solved linear systems by handing them to LAPACK. Here we turn from algebra to calculus done numerically, and the flavour is different: the methods are short enough to write yourself in a dozen lines, but each one carries an error you must be able to name, predict, and measure. The through-line of the chapter is that error — where it comes from (the terms a Taylor series drops), how fast it shrinks as you refine the step (the method's order), and the surprising point where shrinking the step stops helping and starts hurting, straight out of Chapter 20's arithmetic. By the end you will be able to compute a derivative without ever writing down its formula, integrate a function your calculus teacher could not, and — most importantly — put a rigorous number on how wrong the answer is.

In this chapter, you will learn to:

• Approximate a first derivative three ways — forward, backward, and central differences — and read each one's accuracy straight off a Taylor expansion.
• Say precisely what "$O(h)$" and "$O(h^2)$" mean for a method, and why the second-order formula is worth a few extra function calls.
• Estimate a definite integral with the trapezoidal rule, Simpson's rule, and Gauss-Legendre quadrature, and verify each against an integral you can do by hand.
• Measure a method's order of accuracy numerically — halve the step, watch the error fall by $2^p$ — and use the same idea to estimate the error and drive adaptive integration.
• Recognize the round-off floor: the moment a smaller step makes a derivative worse, because subtracting nearly equal function values loses digits (Chapter 20's cancellation, returned to bite).
• Know what changes — and what does not — when the integral has more than one dimension.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — §22.1 (the derivatives you will discretize) and §22.4 (how to prove your method converges) are the load-bearing sections; the Project Checkpoint is the skill you will use on your own solver. Read §22.2–22.3 for the integration toolkit and skim the multidimensional note. - 📖 Standard — read straight through; the finite differences here become the heat stencil in Chapter 24 and the ODE steppers in Chapter 23. - 🔧 Legacy — old codes are full of hand-rolled quadrature and one-sided differences; §22.1 and §22.2 are the vocabulary for reading them, and the ⚠️ pitfall in §22.4 explains a class of "it got less accurate when I refined the grid" bugs you will inherit. - ⚡ HPC — the sums in §22.2 are reductions that vectorize and parallelize; §22.3 (Gauss quadrature) is the "more accuracy per function evaluation" lever, and the multidimensional note in §22.5 is where the FLOPS explode. Performance is not accidental: fewer, better-placed samples beat a brute-force grid.

22.1 Finite Differences: Derivatives from Samples

The definition of the derivative you already know is a limit:

$$ f'(x) = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}. $$

A computer cannot take that limit — but it can evaluate the fraction inside it for a small, finite $h$. Stop shrinking $h$ at some concrete value and you have an approximation, and that single act of stopping early is the whole idea of the subject.

Definition (finite difference). A finite difference approximates a derivative by evaluating the function at a few nearby points a finite distance $h$ apart and combining the values, rather than taking the limit $h \to 0$. The step $h$ is the grid spacing (or step size). The name distinguishes it from the exact, infinitesimal difference of calculus: we keep $h$ finite because a machine must. (This is the derivative-approximation sense of the term; the same word names the whole method for solving partial differential equations, which we build in Chapter 24 — the five-point stencil there is made of exactly these differences.)

The most direct approximation just drops the limit from the definition. It is called the forward difference, because it looks forward from $x$ to $x+h$:

$$ f'(x) \approx D_+ f(x) = \frac{f(x+h) - f(x)}{h}. $$

Look backward instead, to $x-h$, and you get the backward difference:

$$ f'(x) \approx D_- f(x) = \frac{f(x) - f(x-h)}{h}. $$

And if you straddle $x$ symmetrically — one step each way — you get the central difference, which will turn out to be markedly better than either one-sided formula for the same $h$:

$$ f'(x) \approx D_0 f(x) = \frac{f(x+h) - f(x-h)}{2h}. $$

Why is central better? Not for any reason you can see by staring at the formulas — you have to expand them. The tool that reveals every finite difference's accuracy is the Taylor series, and it is worth doing the expansion once, by hand, because it is the source of every error estimate in the chapter.

Write $f(x+h)$ and $f(x-h)$ as their Taylor expansions about $x$:

$$ f(x+h) = f(x) + h f'(x) + \tfrac{h^2}{2} f''(x) + \tfrac{h^3}{6} f'''(x) + \cdots $$ $$ f(x-h) = f(x) - h f'(x) + \tfrac{h^2}{2} f''(x) - \tfrac{h^3}{6} f'''(x) + \cdots $$

Substitute the first into the forward difference and cancel $f(x)$:

$$ \frac{f(x+h) - f(x)}{h} = f'(x) + \underbrace{\tfrac{h}{2} f''(x) + \cdots}_{\text{the error}}. $$

The forward difference equals the true derivative plus a leftover that starts at $\tfrac{h}{2} f''(x)$. That leftover is the price of stopping the limit early, and it has a name.

Definition (truncation error). The truncation error (or discretization error) of a finite difference is the difference between the exact derivative and the formula, and it arises from truncating the Taylor series — throwing away the higher-order terms. For the forward difference it is $\tfrac{h}{2}f''(x) + O(h^2)$, so its leading term is proportional to $h$. Truncation error is a property of the formula and the function, entirely separate from the rounding error of Chapter 20; it is the error you would still have on a hypothetical machine with infinite precision.

Now do the same for the central difference. Subtract the two expansions: the $f(x)$ terms cancel, the $f''$ terms cancel (they have the same sign), and — the crucial event — the odd-order terms reinforce:

$$ f(x+h) - f(x-h) = 2h f'(x) + \tfrac{h^3}{3} f'''(x) + \cdots $$

Divide by $2h$:

$$ \frac{f(x+h) - f(x-h)}{2h} = f'(x) + \underbrace{\tfrac{h^2}{6} f'''(x) + \cdots}_{\text{the error}}. $$

The leading error term is now proportional to $h^2$, not $h$. The symmetric formula's leading error cancelled, leaving a smaller one behind. This is the entire reason central differences dominate one-sided ones in practice, and it motivates the vocabulary that organizes the rest of the chapter.

🚪 Threshold Concept — the exponent of $h$ is everything. When we write that the forward difference is $O(h)$ and the central difference is $O(h^2)$, we are naming the power of the step size in the leading error term — and that exponent, not the constant in front of it, decides how the method behaves as you refine. Halve $h$ for the $O(h)$ formula and the error roughly halves; halve $h$ for the $O(h^2)$ formula and the error drops by a factor of four. Refine ten times and the gap is $2^{10} \approx 1000$ versus $4^{10} \approx 10^6$. Once you learn to read a method's error as "$C h^p$" and to care about $p$ above all, the whole of numerical analysis — differences, quadrature, ODE solvers — sorts itself into orders, and you will always reach first for the highest order you can afford.

Here are all three differences in a single program, applied to $f(x) = x^3$ at $x = 2$, where we know the exact answer is $f'(2) = 3x^2 = 12$. A cubic is the ideal first test because every term in the Taylor expansion is something we can compute exactly by hand:

program finite_differences
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), parameter :: x = 2.0_dp, h = 0.1_dp
  real(dp)            :: fwd, bwd, cen, exact

  exact = 3.0_dp * x**2                        ! f'(x) = 3x^2 = 12 for f(x) = x^3
  fwd = (cube(x + h) - cube(x)    ) / h
  bwd = (cube(x)     - cube(x - h)) / h
  cen = (cube(x + h) - cube(x - h)) / (2.0_dp * h)

  print '(a, f12.8)',  'exact  df/dx at 2 = ', exact
  print '(a, f12.8)',  'forward  (O(h))   = ', fwd
  print '(a, f12.8)',  'backward (O(h))   = ', bwd
  print '(a, f12.8)',  'central  (O(h^2)) = ', cen
  print '(a, es12.4)', 'forward  error    = ', fwd - exact
  print '(a, es12.4)', 'central  error    = ', cen - exact
contains
  pure function cube(t) result(y)
    real(dp), intent(in) :: t
    real(dp)             :: y
    y = t**3
  end function cube
end program finite_differences

$ gfortran -std=f2018 -Wall finite_differences.f90 -o fd && ./fd
exact  df/dx at 2 =  12.00000000
forward  (O(h))   =  12.61000000
backward (O(h))   =  11.41000000
central  (O(h^2)) =  12.01000000
forward  error    =   6.1000E-01
central  error    =   1.0000E-02

Work the arithmetic by hand and the orders reveal themselves. The forward difference is $\big((2.1)^3 - 2^3\big)/0.1 = (9.261 - 8)/0.1 = 12.61$, off by $0.61$. Expanding $(2+h)^3 = 8 + 12h + 6h^2 + h^3$ shows the forward estimate is exactly $12 + 6h + h^2$, so the error is $6h + h^2 = 0.61$ — dominated by the $6h$ term, the hallmark of an $O(h)$ method. The central difference is $\big((2.1)^3 - (1.9)^3\big)/0.2 = (9.261 - 6.859)/0.2 = 12.01$, off by only $0.01$: algebra gives the central estimate as exactly $12 + h^2$, so its error is $h^2 = 0.01$, an $O(h^2)$ method. Same step, same function, and the central difference is sixty-one times more accurate — because its leading error term cancelled.

💡 Intuition — the central difference "sees" the curve symmetrically. A forward difference measures the slope of the chord from $x$ to $x+h$; because it only looks to one side, it is fooled by the curvature between $x$ and $x+h$ to first order in $h$. The central difference averages the forward and backward chords (indeed $D_0 = \tfrac12(D_+ + D_-)$), and the two one-sided curvature errors — equal in size, opposite in sign — cancel. Symmetry buys you an order for free, and this is a pattern you will see again and again: symmetric formulas are more accurate than lopsided ones.

The same machinery gives higher derivatives, and one of them is the star of your project. Add the two Taylor expansions instead of subtracting them: the odd terms cancel and the $f''$ terms survive,

$$ f(x+h) - 2f(x) + f(x-h) = h^2 f''(x) + \tfrac{h^4}{12} f^{(4)}(x) + \cdots, $$

which rearranges into the three-point second difference, the one-dimensional heart of the Laplacian:

$$ f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}, \qquad \text{error } \tfrac{h^2}{12} f^{(4)}(x) = O(h^2). $$

Three samples, one subtraction pattern $(1, -2, 1)$, and a second-order-accurate second derivative. When Chapter 24 builds the five-point stencil for $\nabla^2 u$ on the plate, it is precisely this formula applied in $x$ and again in $y$. Your Project Checkpoint below will verify that it really is $O(h^2)$, numerically.

🐍 Python Comparison. NumPy ships numpy.gradient, which applies exactly the central difference in the interior and one-sided differences at the array ends, and SciPy's scipy.misc.derivative (and the newer scipy.differentiate) wrap adaptive finite differences. They are convenient, and for a array of samples they are the right tool. But when the derivative sits inside a hot loop — a stencil swept across a billion-cell grid, a million times — you do not want a Python function call per point; you want the whole-array (u(3:) - u(1:n-2)) / (2*h) written in Fortran, which the compiler turns into one vectorized pass over memory. This is the same story as everywhere in the book: Python for convenience at the boundary, Fortran for the kernel that runs a trillion times.

🔄 Check Your Understanding. 1. Write the forward, backward, and central difference for $f'(x)$, and state the order of each. 2. From the Taylor expansion, why does the central difference's leading error term vanish while the forward difference's does not? 3. What is the three-point formula for $f''(x)$, and what is its order of accuracy?

Answers

(1) Forward $\big(f(x+h)-f(x)\big)/h$ and backward $\big(f(x)-f(x-h)\big)/h$ are both $O(h)$; central $\big(f(x+h)-f(x-h)\big)/(2h)$ is $O(h^2)$. (2) Subtracting $f(x-h)$ from $f(x+h)$ cancels the even-order Taylor terms (including the $\tfrac{h^2}{2}f''$ term that is the forward difference's leading error), leaving the first surviving error at order $h^2$; the forward difference keeps that $f''$ term, so its error starts at order $h$. (3) $f''(x) \approx \big(f(x+h) - 2f(x) + f(x-h)\big)/h^2$, order $O(h^2)$, with leading error $\tfrac{h^2}{12}f^{(4)}(x)$.

22.2 Numerical Integration: The Trapezoidal and Simpson's Rules

Integration is differentiation's friendlier twin. Differentiation subtracts nearby values and divides by a small number — a recipe for amplifying error, as §22.4 will show. Integration adds values and multiplies by a small number, which is numerically gentle: errors tend to average out rather than blow up. The general name for the enterprise is quadrature, and the general shape of every rule is the same — sample the integrand at chosen points and take a weighted sum.

Definition (quadrature). Quadrature is the numerical approximation of a definite integral $\int_a^b f(x)\,dx$ by a weighted sum of function values, $\sum_{i} w_i\, f(x_i)$, at a finite set of nodes $x_i$ with weights $w_i$. The word is inherited from ancient geometry, where "quadrature" meant constructing a square of equal area — literally finding the area under a curve. A quadrature rule is defined entirely by its nodes and weights; the different rules in this chapter are just different choices of where to sample and how much to count each sample.

The simplest useful rule joins adjacent samples with straight lines and sums the areas of the resulting trapezoids. Over a single interval $[a, b]$ the area under the chord is $\tfrac{b-a}{2}\big(f(a) + f(b)\big)$. To make it accurate you cut $[a,b]$ into $n$ equal panels of width $h = (b-a)/n$, apply the trapezoid on each, and add — collecting the shared interior points, each of which is an endpoint of two panels and so appears with full weight:

$$ \int_a^b f(x)\,dx \approx T_n = h\left[\tfrac{1}{2}f_0 + f_1 + f_2 + \cdots + f_{n-1} + \tfrac{1}{2}f_n\right], \qquad f_i = f(a + i h). $$

This is the trapezoidal rule (composite form). Its error, from integrating the Taylor remainder of the chord, is $-\tfrac{(b-a)h^2}{12} f''(\xi)$ for some $\xi \in [a,b]$ — an $O(h^2)$ method, second-order like the central difference and for the same underlying reason (a symmetric, linear fit).

Simpson's rule does better by fitting parabolas instead of straight lines. Take the panels two at a time, pass a parabola through the three points $(x_{i-1}, x_i, x_{i+1})$, and integrate the parabola exactly. The weights that fall out follow the memorable $1\text{–}4\text{–}2\text{–}4\text{–}\cdots\text{–}4\text{–}1$ pattern (endpoints once, odd interior nodes four times, even interior nodes twice), and $n$ must be even:

$$ \int_a^b f(x)\,dx \approx S_n = \tfrac{h}{3}\left[f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 2f_{n-2} + 4f_{n-1} + f_n\right]. $$

Simpson's rule has error $-\tfrac{(b-a)h^4}{180} f^{(4)}(\xi)$ — it is $O(h^4)$, fourth-order, two orders better than the trapezoid for the cost of nothing but rearranged weights. That two-order jump is why Simpson is the default hand-rolled quadrature rule in scientific code.

Here are both, as pure functions that take the integrand as a procedure argument — the modern-Fortran way to write a routine that works for any $f$, using the procedure interface you met in Chapter 6:

module quadrature_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: dp, trapezoid, simpson

  abstract interface
    pure function scalar_fn(x) result(y)   ! any real -> real function
      import :: dp
      real(dp), intent(in) :: x
      real(dp)             :: y
    end function scalar_fn
  end interface

contains

  pure function trapezoid(f, a, b, n) result(s)
    procedure(scalar_fn)  :: f
    real(dp), intent(in)  :: a, b
    integer,  intent(in)  :: n
    real(dp)              :: s, h
    integer               :: i
    h = (b - a) / real(n, dp)
    s = 0.5_dp * (f(a) + f(b))
    do i = 1, n - 1
      s = s + f(a + real(i, dp) * h)
    end do
    s = s * h
  end function trapezoid

  pure function simpson(f, a, b, n) result(s)   ! n must be even
    procedure(scalar_fn)  :: f
    real(dp), intent(in)  :: a, b
    integer,  intent(in)  :: n
    real(dp)              :: s, h
    integer               :: i
    h = (b - a) / real(n, dp)
    s = f(a) + f(b)
    do i = 1, n - 1
      if (mod(i, 2) == 1) then
        s = s + 4.0_dp * f(a + real(i, dp) * h)   ! odd node, weight 4
      else
        s = s + 2.0_dp * f(a + real(i, dp) * h)   ! even node, weight 2
      end if
    end do
    s = s * h / 3.0_dp
  end function simpson

end module quadrature_mod

program integrate_demo
  use quadrature_mod, only: dp, trapezoid, simpson
  implicit none

  print '(a, f12.8)', 'trapezoid  x^2  n=2 : ', trapezoid(sq,      0.0_dp, 1.0_dp, 2)
  print '(a, f12.8)', 'trapezoid  x^2  n=4 : ', trapezoid(sq,      0.0_dp, 1.0_dp, 4)
  print '(a, f12.8)', 'simpson    x^2  n=2 : ', simpson  (sq,      0.0_dp, 1.0_dp, 2)
  print '(a, f12.8)', 'simpson    x^4  n=2 : ', simpson  (quartic, 0.0_dp, 1.0_dp, 2)
  print '(a, f12.8)', 'simpson    x^4  n=4 : ', simpson  (quartic, 0.0_dp, 1.0_dp, 4)

contains
  pure function sq(x) result(y)
    real(dp), intent(in) :: x
    real(dp)             :: y
    y = x * x
  end function sq
  pure function quartic(x) result(y)
    real(dp), intent(in) :: x
    real(dp)             :: y
    y = x**4
  end function quartic
end program integrate_demo

$ gfortran -std=f2018 -Wall quadrature.f90 -o quad && ./quad
trapezoid  x^2  n=2 :   0.37500000
trapezoid  x^2  n=4 :   0.34375000
simpson    x^2  n=2 :   0.33333333
simpson    x^4  n=2 :   0.20833333
simpson    x^4  n=4 :   0.20052083

Every one of those numbers is worth checking by hand, because doing so proves you understand the rules rather than trusting them. The true value $\int_0^1 x^2\,dx = \tfrac13 = 0.3333\ldots$ Trapezoid with two panels ($h = \tfrac12$) is $\tfrac12\big[\tfrac12 f(0) + f(\tfrac12) + \tfrac12 f(1)\big] = \tfrac12\big[0 + 0.25 + 0.5\big] = 0.375$; with four panels it tightens to $0.34375$. Halving $h$ took the error from $0.0417$ to $0.0104$ — a factor of four, the $O(h^2)$ signature. Now the striking line: Simpson on $x^2$ returns $0.33333333$ — the exact answer. That is not luck. Simpson fits parabolas, and a parabola through three points of a parabola is that parabola, so Simpson integrates any quadratic exactly; in fact its error term $f^{(4)}$ vanishes for any cubic too, so Simpson is exact for all polynomials up to degree three.

To see Simpson's finite error you must feed it something it cannot fit exactly, so the last two lines use $\int_0^1 x^4\,dx = \tfrac15 = 0.2$. Simpson with two panels gives $\tfrac{0.5}{3}\big[f(0) + 4f(0.5) + f(1)\big] = \tfrac{0.5}{3}\big[0 + 4(0.0625) + 1\big] = \tfrac{0.5}{3}(1.25) = 0.20833333$, error $0.00833$; with four panels, $0.20052083$, error $0.00052$. The ratio is $0.00833 / 0.00052 = 16$ — halving $h$ cut the error by sixteen, which is $2^4$, the unmistakable fingerprint of a fourth-order method. Hold onto that number 16; §22.4 makes it the basis of measuring order and of adaptive integration.

⚡ Performance Note — quadrature is a reduction, and reductions parallelize. Look at the shape of both rules: evaluate $f$ at many independent points, then sum. The function evaluations have no dependencies on one another, and the sum is a reduction — exactly the pattern that vectorizes and parallelizes cleanly. On a single core the compiler can vectorize the accumulation loop; across cores, an OpenMP reduction(+:s) clause (Chapter 33) splits the panels among threads and combines their partial sums. The one caution is the one from Chapter 20: floating-point addition is not associative, so a parallel sum can differ from a serial one in its last bits. For a well-behaved integrand that is harmless; when it matters, a compensated (Kahan) summation restores reproducibility.

🔗 Connection — the numerical libraries are Fortran, again. You will rarely ship your own Simpson's rule for serious work, for the same reason you did not ship your own dgesv in Chapter 21: a battle-tested library does it better. The standard one is QUADPACK, a Fortran 77 automatic-quadrature library from the 1980s, and it is what sits underneath scipy.integrate.quad in Python — so when a data scientist integrates a function in SciPy, the arithmetic happens in the same decades-old, still-unbeaten Fortran this chapter teaches you to read and write. Fortran is not dead; it is under the notebook.

22.3 Gaussian Quadrature and Adaptive Integration

The trapezoidal and Simpson's rules sample at equally spaced points and ask: given fixed nodes, what are the best weights? Gaussian quadrature asks a bolder question: if we are free to choose the nodes as well, where should we put them? The answer is remarkable. With $n$ points chosen optimally, you can integrate every polynomial up to degree $2n - 1$ exactly — twice the reach of any equally-spaced rule with the same number of samples. Two well-placed points do the work of Simpson's three.

The optimal nodes are not obvious — they are the roots of the Legendre polynomials, and the weights come with them — but you look them up in a table, you do not derive them at run time. The rules are tabulated on the reference interval $[-1, 1]$; the two-point Gauss-Legendre rule is

$$ \int_{-1}^{1} g(t)\,dt \approx g\!\left(-\tfrac{1}{\sqrt{3}}\right) + g\!\left(+\tfrac{1}{\sqrt{3}}\right), \qquad \text{nodes } \pm\tfrac{1}{\sqrt{3}} \approx \pm 0.57735, \ \text{weights } 1, 1. $$

To integrate over a general interval $[a, b]$ you change variables with the linear map $x = \tfrac{b-a}{2}t + \tfrac{a+b}{2}$, which stretches $[-1,1]$ onto $[a,b]$ and contributes a Jacobian factor $\tfrac{b-a}{2}$:

$$ \int_a^b f(x)\,dx = \tfrac{b-a}{2}\int_{-1}^{1} f\!\left(\tfrac{b-a}{2}t + \tfrac{a+b}{2}\right)dt \approx \tfrac{b-a}{2}\sum_{i} w_i\, f\!\left(\tfrac{b-a}{2}t_i + \tfrac{a+b}{2}\right). $$

Here is the two-point rule, with the transformation built in:

module gauss_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: dp, gauss2

  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

contains

  pure function gauss2(f, a, b) result(s)
    procedure(scalar_fn) :: f
    real(dp), intent(in) :: a, b
    real(dp)             :: s, half, mid
    real(dp), parameter  :: node = 0.5773502691896257_dp   ! 1/sqrt(3)
    half = 0.5_dp * (b - a)
    mid  = 0.5_dp * (a + b)
    s = half * ( f(mid - half*node) + f(mid + half*node) )  ! weights are both 1
  end function gauss2

end module gauss_mod

program gauss_demo
  use gauss_mod, only: dp, gauss2
  implicit none
  print '(a, f12.8)', 'Gauss-2  x^2 on [0,1] : ', gauss2(sq,      0.0_dp, 1.0_dp)
  print '(a, f12.8)', 'Gauss-2  x^4 on [0,1] : ', gauss2(quartic, 0.0_dp, 1.0_dp)
contains
  pure function sq(x) result(y)
    real(dp), intent(in) :: x
    real(dp)             :: y
    y = x * x
  end function sq
  pure function quartic(x) result(y)
    real(dp), intent(in) :: x
    real(dp)             :: y
    y = x**4
  end function quartic
end program gauss_demo

$ gfortran -std=f2018 -Wall gauss.f90 -o gauss && ./gauss
Gauss-2  x^2 on [0,1] :   0.33333333
Gauss-2  x^4 on [0,1] :   0.19444444

The first line is the payoff. Two function evaluations — at $x = \tfrac12 \mp \tfrac{1}{2\sqrt3} \approx 0.2113$ and $0.7887$ — reproduce $\int_0^1 x^2\,dx = \tfrac13$ *exactly*, because $x^2$ has degree $2 \le 2n - 1 = 3$. You can verify it without the roots: with the nodes written as $\tfrac12 \pm c$ where $c = \tfrac{1}{2\sqrt3}$, the sum $x_1^2 + x_2^2 = (\tfrac12 - c)^2 + (\tfrac12 + c)^2 = \tfrac12 + 2c^2 = \tfrac12 + \tfrac16 = \tfrac23$, times the Jacobian $\tfrac12$ gives $\tfrac13$. The second line shows the limit: $x^4$ has degree $4 > 3$, so two points can no longer be exact, and Gauss-2 returns $0.19444$, an error of $-0.0056$. But notice it still beats Simpson's three-point result on the same integral ($0.20833$, error $+0.0083$) — fewer, smarter samples. A three-point Gauss rule (nodes $0, \pm\sqrt{3/5}$, weights $\tfrac89, \tfrac59, \tfrac59$) would nail $x^4$ exactly, being good to degree $5$.

🚪 Threshold Concept — accuracy per function evaluation is the real currency. When each call to $f$ is expensive — and in real science one "function evaluation" can be an entire subsimulation costing seconds or minutes — the figure of merit is not error per unit $h$ but error per function call. That reframing is why Gaussian quadrature exists and why it dominates smooth integrands: it extracts the most accuracy from the fewest, most carefully placed samples. The same instinct — spend your evaluations where they buy the most — drives adaptive methods next, and it is the numerical scientist's version of the performance theme that runs through this whole book: the fastest computation is the one you arrange never to do.

Gaussian quadrature is superb for smooth integrands, but it commits to its node count in advance. What if the integrand is placid over most of $[a,b]$ and viciously curved in one small region — a spike, a boundary layer, a near-singularity? A fixed rule either wastes samples on the flat part or starves the hard part. The cure is adaptive integration: put the samples where the function needs them.

Definition (adaptive integration). Adaptive integration is any quadrature scheme that automatically refines its sampling where the integrand is hardest, guided by a run-time error estimate. The standard construction is recursive: apply a rule to an interval, apply it again to the two halves, and compare. If the two results agree to within the requested tolerance, accept the finer one; if not, subdivide each half and recurse. Effort concentrates automatically in the regions of large error, and the user specifies an accuracy, not a node count.

The engine of every adaptive method is a cheap, reliable error estimate, and you already have everything needed to build one. Because Simpson's rule is $O(h^4)$, halving $h$ cuts its error by $16$; so if $S_n$ and $S_{2n}$ are the coarse and fine estimates, their true errors are in the ratio $16 : 1$, and a little algebra gives the error of the finer one almost for free:

$$ \text{error}(S_{2n}) \approx \frac{S_{2n} - S_n}{15}. $$

Test it on the numbers you already computed for $\int_0^1 x^4\,dx$: $S_2 = 0.208333$ and $S_4 = 0.200521$, so the estimate is $(0.200521 - 0.208333)/15 = -0.000521$ — and the true error of $S_4$ is $0.200521 - 0.2 = 0.000521$. The estimate nailed the magnitude without ever knowing the true answer. That is exactly what an adaptive routine leans on: it never knows the true integral, but it can estimate its own error by comparing two resolutions, and refine until that estimate drops below your tolerance.

Push the same idea one step further and you get Richardson extrapolation: instead of merely estimating the error, subtract it. Adding the estimate back, $S_4 + (S_4 - S_2)/15 = 0.200521 - 0.000521 = 0.200000$ — the exact answer, because for a quartic the extrapolation cancels the entire $h^4$ error term. Applying this systematically to the trapezoidal rule at successively halved $h$ is Romberg integration, a compact and classic way to wring high-order accuracy out of the humble trapezoid.

Definition (Richardson extrapolation). Richardson extrapolation combines two approximations of different step size to cancel the leading error term and produce a higher-order estimate. If a method has error $C h^p$ and you have results $A(h)$ and $A(h/2)$, then the combination $\big(2^p A(h/2) - A(h)\big)/(2^p - 1)$ eliminates the $h^p$ term, leaving an error of higher order. It is the algebraic engine behind Romberg integration and behind the error estimates that drive adaptive quadrature — a general-purpose lever you apply whenever you know a method's order.

🔗 Connection — production adaptive quadrature. The adaptive Simpson sketch here is the teaching version; the professional article is QUADPACK's family of routines (qag, qags, qagi for infinite intervals), which pair high-order Gauss-Kronrod rules with exactly this compare-and-subdivide logic. When you call scipy.integrate.quad and it returns both an answer and an error estimate, that second number is the adaptive error control of this section, computed by Fortran you now understand in principle.

🔄 Check Your Understanding. 1. An $n$-point Gauss-Legendre rule is exact for polynomials up to what degree? 2. Simpson estimates $S_2 = 1.4627$ and $S_4 = 1.4637$ for some integral. Estimate the error of $S_4$. 3. In one sentence, how does an adaptive routine decide where to place more points?

Answers

(1) Degree $2n - 1$ — the reason two Gauss points reach as far as Simpson's three equally-spaced ones. (2) $\text{error}(S_4) \approx (S_4 - S_2)/15 = (1.4637 - 1.4627)/15 = 0.001/15 \approx 6.7\times10^{-5}$. (3) It subdivides any interval whose local error estimate (coarse-vs-fine disagreement) exceeds the tolerance, so effort concentrates where the integrand is hardest to fit.

22.4 Error Analysis and Convergence

Every method so far came with an error term — $O(h)$, $O(h^2)$, $O(h^4)$ — pulled from a Taylor expansion. Those are theoretical orders. The mark of a careful computational scientist is to never take them on faith: you measure the order of your own code, on your own machine, and confirm it matches the theory. If it does not, you have a bug — a mis-typed weight, an off-by-one in a loop bound, a boundary handled wrong — and the convergence test is how you find it. This section makes "order of accuracy" precise and shows you the two-line experiment that verifies it.

Definition (order of accuracy). A method has order of accuracy $p$ if its error behaves like $E(h) \approx C h^{p}$ as the step $h \to 0$, for some constant $C$ independent of $h$. Equivalently, the error is $O(h^p)$. The order is the exponent $p$: forward and backward differences are first-order ($p=1$), central differences, the trapezoidal rule, and the three-point second difference are second-order ($p=2$), and Simpson's rule is fourth-order ($p=4$). A higher order means the error falls faster as you refine, so above some resolution a higher-order method is not just better but dramatically better.

Definition (convergence). A numerical method converges if its approximation approaches the exact answer as the discretization is refined — $E(h) \to 0$ as $h \to 0$ (equivalently, as the number of nodes $n \to \infty$). Convergence is the minimum requirement of any method worth using; the order of accuracy then says how fast it converges. A method that does not converge, or converges only until round-off intervenes (below), is telling you something is wrong — with the method, the code, or the problem's conditioning.

The experiment that measures $p$ is disarmingly simple. If $E(h) \approx C h^p$, then halving the step changes the error by a fixed ratio,

$$ \frac{E(h)}{E(h/2)} = \frac{C h^p}{C (h/2)^p} = 2^p, $$

independent of the unknown constant $C$. So you compute the error at $h, h/2, h/4, \ldots$, take successive ratios, and read off the order: a ratio near $2$ means first-order, near $4$ means second-order, near $16$ means fourth-order. Here is the whole method, applied to the forward ($O(h)$) and central ($O(h^2)$) differences of $f(x) = x^3$ at $x = 2$ so you can check every digit:

program convergence_study
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), parameter :: x = 2.0_dp, exact = 12.0_dp   ! d/dx x^3 = 3x^2 = 12 at x=2
  real(dp) :: h, ef, ec
  integer  :: k

  print '(a)', '      h            forward_err         central_err'
  h = 0.1_dp
  do k = 1, 4
    ef = abs((cube(x+h) - cube(x)   ) / h          - exact)   ! forward,  O(h)
    ec = abs((cube(x+h) - cube(x-h) ) / (2.0_dp*h) - exact)   ! central,  O(h^2)
    print '(f9.5, es18.6, es18.6)', h, ef, ec
    h = h / 2.0_dp
  end do
contains
  pure function cube(t) result(y)
    real(dp), intent(in) :: t
    real(dp)             :: y
    y = t**3
  end function cube
end program convergence_study

$ gfortran -std=f2018 -Wall convergence_study.f90 -o conv && ./conv
      h            forward_err         central_err
  0.10000      6.100000E-01      1.000000E-02
  0.05000      3.025000E-01      2.500000E-03
  0.02500      1.506250E-01      6.250000E-04
  0.01250      7.515625E-02      1.562500E-04

Take successive ratios of each column and the orders leap off the page. The forward-error column falls by a factor of about two at each halving — $0.610000 \to 0.302500 \to 0.150625 \to 0.075156$, ratios $2.02, 2.01, 2.00$ — confirming $p = 1$; the small excess above $2$ is the higher-order tail shrinking away, and the ratio marches toward exactly $2$ as $h \to 0$. The central-error column falls by a factor of exactly four — $0.010000 \to 0.002500 \to 0.000625 \to 0.000156$ — confirming $p = 2$, and it is exactly $4$ here because the central difference of a cubic has error precisely $h^2$, with no higher-order tail. This four-line table is the single most useful diagnostic in numerical computing: whenever you write a method whose order you know, refine, take ratios, and confirm. When the ratio comes out wrong, trust the table over your ego and go find the bug.

Now the twist that makes this a Fortran chapter and not just a calculus one. The theory says $E(h) \to 0$ as $h \to 0$, so smaller $h$ is always better. On a real machine, it is not — and the reason is the entire subject of Chapter 20.

⚠️ Common Pitfall — a smaller $h$ is not always better. A finite-difference derivative subtracts two nearly equal function values, $f(x+h)$ and $f(x)$, and divides by a tiny $h$. As $h$ shrinks, $f(x+h)$ and $f(x)$ agree in more and more leading digits, so their difference suffers catastrophic cancellation (Chapter 20) — the subtraction throws away the digits they share and promotes their round-off into the result, which is then amplified by the division by a small $h$. Truncation error falls with $h$, but round-off error grows like $\varepsilon/h$. Their sum is U-shaped: the total error decreases, bottoms out at an optimal step $h^*$, and then gets worse as round-off takes over. Pushing $h$ below $h^*$ makes your derivative less accurate, not more.

You can estimate exactly where the bottom lies. For the forward difference, the truncation error is about $\tfrac{h}{2}|f''|$ and the round-off error about $2\varepsilon|f|/h$ (a rounding of size $\varepsilon|f|$ in the numerator, divided by $h$). Their sum is smallest when the two are balanced, at

$$ h^{*} \approx 2\sqrt{\dfrac{\varepsilon\, |f|}{|f''|}} \sim \sqrt{\varepsilon} \approx 10^{-8} \quad(\text{double precision}), $$

giving a best-possible error of only about $\sqrt{\varepsilon} \approx 10^{-8}$ — you lose half your sixteen digits, no matter how carefully you code, simply because the operation is subtraction of near-equals. The central difference, with truncation error $\tfrac{h^2}{6}|f'''|$, balances round-off at a larger $h^ \sim \varepsilon^{1/3} \approx 10^{-5}$ and reaches a better floor near $\varepsilon^{2/3} \approx 10^{-11}$ — another reason to prefer it. The following table is illustrative* (the exact digits in the round-off-dominated regime are machine-dependent, so we do not print them as certified output), but its shape is real and reproducible:

  forward-difference error for sin'(1),  double precision (illustrative)
    h        total error      regime
  1e-1       ~4e-2            truncation dominates  (error ~ h/2)
  1e-3       ~4e-4            truncation dominates
  1e-6       ~5e-7            approaching the floor
  1e-8       ~6e-9            *** near the optimum h* ***
  1e-10      ~7e-7            round-off dominates   (error ~ eps/h)
  1e-13      ~1e-3            round-off dominates

💡 Intuition — truncation and round-off pull in opposite directions. Think of two error curves on a log–log plot against $h$: truncation error is a line sloping down to the right (smaller $h$, smaller error, slope $p$), round-off error is a line sloping up to the right (smaller $h$, bigger cancellation). The total is their sum, a V (or a smooth valley) whose bottom is the best you can do. Higher-order methods have a steeper truncation line, so they reach a lower valley floor at a larger $h$ — they need less refinement to get more accuracy, and they bottom out before cancellation gets bad. This picture, drawn once, explains why "just use a tiny $h$" is one of the most common and most confident mistakes beginners make.

🔗 Connection — integration does not have this problem. Notice the round-off floor afflicts differentiation, not integration. Quadrature adds function values and multiplies by a small $h$; there is no subtraction of near-equals, so shrinking $h$ (more panels) keeps improving the answer until you simply run out of patience or accumulate round-off in the sum (which grows only like $\sqrt{n}\,\varepsilon$, gently). This asymmetry — differentiation is numerically dangerous, integration is numerically safe — is worth internalizing: it is why we discretize $\tfrac{\partial u}{\partial t}$ and $\nabla^2 u$ with the smallest defensible $h$, but never with an absurdly small one.

22.5 Multidimensional Integration: A Note

Everything so far integrated a function of one variable. Real problems often want the integral of a field over an area or a volume — the total thermal energy stored in your plate, $\int\!\!\int_\Omega u(x,y)\,dx\,dy$, is a two-dimensional integral over the domain $\Omega$. The good news is that the one-dimensional rules extend directly; the bad news is a scaling wall that reshapes how you approach high dimensions.

The direct extension is a product rule: apply a 1-D rule in each direction and multiply the weights. In two dimensions,

$$ \int_c^d\!\!\int_a^b f(x,y)\,dx\,dy \approx \sum_{i}\sum_{j} w_i\, w_j\, f(x_i, y_j), $$

which in Fortran is simply a nested loop over a grid of nodes, accumulating $w_i w_j f(x_i, y_j)$ — a do concurrent or a doubly-nested reduction, and every bit as vectorizable as the 1-D sum. For a rectangular domain with a smooth integrand this is exactly right, and a 2-D Simpson or Gauss product rule is a dozen lines. The order of accuracy carries over: a product of $O(h^p)$ rules is $O(h^p)$ in each direction.

The wall is the curse of dimensionality. A product rule with $n$ nodes per dimension uses $n^d$ nodes in $d$ dimensions. Ten nodes per axis is a reasonable 1-D rule; in three dimensions it is $1{,}000$ function evaluations, in six dimensions $10^6$, and by the time a computational chemist integrates over the coordinates of a handful of particles — dozens of dimensions — a product grid would need more evaluations than there are atoms in the room. Gaussian quadrature's efficiency does not save you: $n^d$ grows explosively no matter how good each axis is.

Above a handful of dimensions the answer is to stop laying down a grid at all and sample randomly. Monte Carlo integration estimates $\int_\Omega f\,dV \approx \tfrac{V}{N}\sum_{k=1}^{N} f(\mathbf{x}_k)$ by averaging the integrand over $N$ points thrown uniformly into the domain. Its error shrinks like $1/\sqrt{N}$ — slowly, and regardless of the dimension $d$. That last clause is the whole point: a grid-based rule's cost to reach a fixed accuracy explodes with $d$, while Monte Carlo's does not, so above some crossover (often around $d \approx 4$–$8$) randomness wins decisively. It is how high-dimensional integrals in finance, particle physics, and statistical mechanics are actually computed, and the reason a "slow" $1/\sqrt{N}$ method is a workhorse rather than a curiosity.

🔗 Connection. Your heat solver lives in two dimensions, comfortably inside product-rule territory — when Chapter 24 computes a diagnostic like the plate's total or average temperature, a 2-D trapezoid over the grid is the honest tool, and it is nothing more than a double sum with half-weight edges. Monte Carlo waits in the wings for the day your model grows a high-dimensional parameter space to integrate over — an uncertainty quantification, say — but for a field on a 2-D grid, the deterministic product rule is both faster and more accurate.

Project Checkpoint

Your heat solver's spatial core is the discrete Laplacian — the second derivative of temperature, computed with the three-point stencil $\big(u_{i-1} - 2u_i + u_{i+1}\big)/h^2$ from §22.1 (applied in $x$ and $y$ to make Chapter 24's five-point stencil). Everything the simulation predicts rests on that stencil being as accurate as we claim: second-order accurate in the grid spacing. This checkpoint proves it, numerically, with the convergence experiment of §22.4 — the same test you should run on every discretization you ever write.

The plan: pick a function whose second derivative you know exactly, apply the stencil on a sequence of refining grids $h, h/2, h/4, \ldots$, and confirm the error falls by a factor of $4$ at each halving — the $O(h^2)$ signature. We use $f(x) = x^4$ at $x = 1$, where $f''(x) = 12x^2 = 12$, because the algebra is exact: $\big((1+h)^4 - 2 + (1-h)^4\big)/h^2 = 12 + 2h^2$, so the stencil's error is exactly $2h^2$, and it must quarter when $h$ halves.

program stencil_order
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), parameter :: x = 1.0_dp, exact = 12.0_dp   ! d2/dx2 x^4 = 12x^2 = 12 at x=1
  real(dp)         :: h, err, err_prev, d2
  character(len=9) :: ratio
  integer          :: k

  print '(a)', '     h         stencil_d2      error        ratio'
  h = 0.1_dp
  err_prev = 0.0_dp
  do k = 1, 4
    d2  = (quartic(x+h) - 2.0_dp*quartic(x) + quartic(x-h)) / h**2   ! three-point 2nd difference
    err = abs(d2 - exact)
    if (k == 1) then
      ratio = '     --  '
    else
      write(ratio, '(f9.2)') err_prev / err
    end if
    print '(f9.5, f15.8, es14.4, a9)', h, d2, err, ratio
    err_prev = err
    h = h / 2.0_dp
  end do
contains
  pure function quartic(t) result(y)
    real(dp), intent(in) :: t
    real(dp)             :: y
    y = t**4
  end function quartic
end program stencil_order

$ gfortran -std=f2018 -Wall stencil_order.f90 -o stencil && ./stencil
     h         stencil_d2      error        ratio
  0.10000    12.02000000    2.0000E-02     --  
  0.05000    12.00500000    5.0000E-03     4.00
  0.02500    12.00125000    1.2500E-03     4.00
  0.01250    12.00031250    3.1250E-04     4.00

The error column is $2h^2 = 0.02, 0.005, 0.00125, 0.0003125$, and the ratio column is a flat $4.00$ — the stencil converges at exactly second order, as advertised. You have now verified, not assumed, the accuracy of the operator at the heart of your solver: if Chapter 24's $\nabla^2 u$ ever fails to show this ratio-4 behaviour on a refined grid, the stencil is coded wrong, and this test will catch it. Save the program as heat-solver/tests/test_stencil_order.f90; it is the first member of the regression suite that Chapter 37 turns into an automated, always-run guardrail. One caution from §22.4: do not push $h$ far below the values here — around $h \approx 10^{-5}$ the numerator's cancellation begins to spoil the ratio, and the "convergence" would reverse. Refine enough to confirm the order, and no further.

Summary

This chapter did calculus the way a computer must — from samples, with a quantified error.

The three things to memorize. First, the central difference is $O(h^2)$ and the one-sided differences are $O(h)$ — symmetry buys an order. Second, to check any method's order, halve $h$ and confirm the error ratio is $2^p$. Third, a finite-difference derivative has an optimal step near $\sqrt\varepsilon$; going smaller makes it worse, because subtraction of near-equals (Chapter 20) takes over.

Spaced Review

Four questions revisiting Chapter 6 (procedures) and Chapter 20 (floating point). Answer from memory before opening the details.

1. (Ch. 6) The trapezoid and simpson functions in §22.2 take the integrand f as a procedure argument and are themselves declared pure. What must be true of the abstract interface scalar_fn for the pure on the integrators to be legal, and why does pure help the compiler here?
   

   Answer

   The interface must itself declare `pure function` — a `pure` procedure may only call other `pure` procedures, so the dummy `f` must be known-pure, which it inherits from a `pure` abstract interface. `pure` promises no side effects and no hidden state, which lets the optimizer reorder, vectorize, and even parallelize the evaluations of `f` across the panels (the payoff previewed in [Chapter 27](../../part-07-performance/chapter-27-why-fortran-is-fast/index.md)).

2. (Ch. 6) Why is passing the integrand as a procedure argument better than hard-coding one function inside simpson, and what modern feature makes it possible to pass an internal function like sq?
   

   Answer

   A procedure argument makes `simpson` reusable for *any* integrand without editing it — the same separation-of-concerns that `intent` and modules serve. Fortran 2008+ lets you pass an *internal* procedure (one defined in a program's or procedure's `contains`) as an actual argument, so `sq` and `quartic` can live beside `main` and still be handed to the integrator.

3. (Ch. 20) In one sentence, why does shrinking $h$ below about $10^{-8}$ make a forward-difference derivative less accurate in double precision?
   

   Answer

   Because $f(x+h)$ and $f(x)$ then agree to nearly all 16 digits, so their subtraction is catastrophic cancellation — round-off of size $\varepsilon|f|$ survives and is amplified by the division by the tiny $h$, and this $\varepsilon/h$ round-off error overtakes the shrinking $\tfrac{h}{2}|f''|$ truncation error.

4. (Ch. 20) The convergence table's central-difference error ratio is a clean $4.00$ at moderate $h$ but would drift and then collapse if you continued to $h = 10^{-9}$. Which recurring theme, and which two competing errors, does that illustrate?
   

   Answer

   It illustrates that a `real(dp)` is a grid point, not a real number: truncation error ($\propto h^2$, shrinking) and round-off error from cancellation ($\propto \varepsilon/h$, growing) compete, and once round-off dominates the ratio stops being $4$ and the "error" grows — the U-shaped total error of §22.4, straight from Chapter 20's arithmetic.

What's Next

You can now approximate the two operations of calculus from nothing but function samples, and — more importantly — bound and measure the error when you do. That is precisely the toolkit the next two chapters demand. Chapter 23 marches ordinary differential equations forward in time with Euler's method and Runge-Kutta, and every one of those steppers is a finite-difference formula in disguise, carrying an order of accuracy you will verify with the halving test you just learned. Then Chapter 24 assembles the three-point second difference of this chapter into the five-point Laplacian stencil and turns your plate into a living simulation — the real PDE core of the whole project. The derivative you learned to compute without its formula is about to become the beating heart of the heat equation. Let's march some equations through time.