Case Study 2: A Quadrature Toolkit You Can Trust

"Testing shows the presence, not the absence of bugs." — Edsger W. Dijkstra

Executive Summary

You need to integrate — repeatedly, in production — both analytic functions (a beam-intensity model) and sampled data (a measured profile), and you need to trust the results. Rather than sprinkle ad-hoc sums through the code, you will build one small, reusable quadrature module and, crucially, validate it against integrals whose answers you already know before pointing it at anything real. The deliverable is a module exposing the trapezoidal rule, Simpson's rule, two- and three-point Gauss-Legendre, and a convergence-order tester — each verified to converge at its advertised order — plus a demonstration that uses the right tool for each job: Gauss for the smooth analytic integrand (maximum accuracy per evaluation) and the trapezoidal rule for the sampled data (where you cannot choose the nodes).

Skills applied - Designing quadrature routines as pure functions with procedure arguments (§22.2, §22.3, Chapter 6). - Validating a numerical method by measuring its order of accuracy (§22.4). - Choosing a rule by accuracy per function evaluation (§22.3) and by whether nodes are yours to place. - Integrating sampled data, where only the trapezoidal rule (or data-Simpson) applies (§22.2, §22.5).

Background

The instrument produces a one-dimensional intensity profile $I(x)$. Two numbers matter downstream: the total power $P = \int I\,dx$ (a normalization) and the centroid $\bar{x} = \int x\,I\,dx / \int I\,dx$ (where the beam is centered). Sometimes $I$ is an analytic model you can evaluate anywhere; sometimes it is a table of detector samples. A single validated toolkit should serve both, and — because these numbers feed a calibration — it must be demonstrably correct, not merely plausible.

The engineering principle, borrowed from Case Study 1's hard-won lesson, is: never trust a numerical routine you have not watched converge at the right order on a problem whose answer you know.

Phase 1 — Design the Module

The interface mirrors §22.2–§22.3: each rule is a pure function taking the integrand as a procedure(scalar_fn) argument. We add gauss3 (from Exercise 22.19) and a measure_order routine that runs the §22.4 halving experiment automatically.

module quadrature
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: dp, trapezoid, simpson, gauss2, gauss3, trapz_data

  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
  ! trapezoid, simpson, gauss2, gauss3 as in sections 22.2-22.3 and Exercise 22.19 ...

  pure function trapz_data(x, y) result(area)     ! integrate SAMPLED data (nodes not ours to choose)
    real(dp), intent(in) :: x(:), y(:)            ! possibly non-uniform abscissae
    real(dp)             :: area
    integer              :: i
    area = 0.0_dp
    do i = 1, size(x) - 1
      area = area + 0.5_dp * (x(i+1) - x(i)) * (y(i) + y(i+1))
    end do
  end function trapz_data
end module quadrature

The key design decision is the split between the function integrators (which may sample wherever they like, so Gauss is available) and trapz_data (which must accept the abscissae it is given). You cannot run Gauss-Legendre on a fixed table of samples — the nodes are not yours to move — so sampled data gets the trapezoidal rule, and that is not a compromise but the correct tool.

Phase 2 — Validate the Function Integrators

Before any real use, point every rule at $\int_0^1 x^4\,dx = \tfrac15 = 0.2$, whose value and whose convergence we can check to the digit.

integral_0^1 x^4 dx  (exact 0.20000000)
  trapezoid n=2 :  0.28125000    error  8.13e-02
  trapezoid n=4 :  0.22070313    error  2.07e-02     ratio 3.93  -> ~2nd order
  simpson   n=2 :  0.20833333    error  8.33e-03
  simpson   n=4 :  0.20052083    error  5.21e-04     ratio 16.0  -> 4th order
  gauss2        :  0.19444444    error  5.56e-03     (2 evals)
  gauss3        :  0.20000000    error  0            (3 evals, exact: degree 4 <= 5)

Read this table as an acceptance test. Trapezoid's error ratio $\approx 4$ confirms $O(h^2)$; Simpson's ratio of exactly $16$ confirms $O(h^4)$; and gauss3 returns the exact answer because $x^4$ has degree $4 \le 2n-1 = 5$. Every number is hand-checkable: Simpson $n=4$ is $\tfrac{0.25}{3}[0 + 4(0.00390625) + 2(0.0625) + 4(0.31640625) + 1] = \tfrac{2.40625}{12} = 0.20052083$, and `gauss2` is $\tfrac12(x_1^4 + x_2^4) = \tfrac12(0.38888889) = 0.19444444$ (from Exercise 22.8's algebra). A rule that failed to show its order here would be rejected before it ever touched real data.

The measured orders are the contract. If a later optimization (say, restructuring Simpson's loop for vectorization) accidentally broke a weight, the ratio would drift from $16$ and this table would fail. Keep it as a regression test.

Phase 3 — Optimize: Accuracy per Evaluation

The validation table already reveals the performance story. gauss2 used two function evaluations and beat simpson n=2's three ($0.0056$ versus $0.0083$ error); gauss3 used three and was exact. For a smooth analytic integrand — where each evaluation might be expensive — Gauss-Legendre extracts more accuracy per call than any equally-spaced rule, which is the §22.3 lesson made quantitative. Performance is not accidental: when function evaluations dominate the cost, the rule that needs fewer of them wins, and you choose it deliberately.

The corollary matters just as much: for the sampled data below, none of this applies. The detector fixed the abscissae, so Gauss is off the table and the trapezoidal rule is both the correct and the only honest choice. Optimization is not "always use the fanciest rule" — it is "use the rule that fits the constraints," and the binding constraint here is who chose the nodes.

Phase 4 — Apply to the Real Profile

Now the payoff. The detector returns five samples of $I(x)$ on a uniform grid $x = 0,1,2,3,4$, with intensities $I = [\,0,\ 1,\ 4,\ 1,\ 0\,]$ — a symmetric peak. We compute the total power $P = \int I\,dx$ and the centroid $\bar{x}$ with trapz_data, sampling $x\,I(x) = [\,0,\ 1,\ 8,\ 3,\ 0\,]$ for the moment:

total power   P     = 6.00000000     (trapezoid on I)
first moment  int xI = 12.0000000    (trapezoid on x*I)
centroid      xbar  = 2.00000000     (= 12 / 6)

Every number is exact by hand: with unit spacing, $P = \tfrac12(0) + 1 + 4 + 1 + \tfrac12(0) = 6$; the first moment is $\tfrac12(0) + 1 + 8 + 3 + \tfrac12(0) = 12$; and $\bar{x} = 12/6 = 2$, dead center, as the symmetric profile demands. The centroid landing exactly at the middle sample is a sanity check the physics hands you for free — a symmetric profile must be centered — and it confirms both the arithmetic and the half-weighted-endpoint convention are right.

Phase 5 — Harden and Hand Off

Three finishing moves make the toolkit production-grade. (1) Wrap the validation table (Phase 2) as an automated test that asserts each order ratio is within tolerance of $4$ or $16$ — the Chapter 37 discipline in miniature. (2) Document, in each routine's header, its order and its exactness degree, and note that simpson requires an even $n$ (guard it with an error stop on odd $n$, the Chapter 13 way). (3) State the selection rule at the top of the module: analytic integrand → Gauss; sampled data → trapz_data; need an error estimate → run two resolutions and compare (§22.3). A colleague who reads that header knows which tool to reach for without re-deriving this case study.

Discussion Questions

  1. Why can gauss3 integrate $x^4$ exactly but not $x^6$? State the degree rule and compute what gauss3 returns for $\int_{-1}^1 x^6\,dx$ (true value $\tfrac27$).
  2. The sampled profile used trapz_data, not Simpson. Under what condition on the number and spacing of samples could you legitimately apply Simpson's rule to the data, and what would you gain?
  3. Phase 2 treats the order ratios as an acceptance test. Design one more test case that would catch a bug where someone swapped the $4$ and $2$ Simpson weights (odd vs even). What integrand exposes it fastest?
  4. The centroid came out exactly at the middle sample. Why is that a weak test on its own, and what asymmetric profile would you add to test the centroid computation properly?

Your Turn: Extensions

  • Option A (build). Add simpson_data(x, y) for an odd number of uniformly spaced samples, and compare its accuracy with trapz_data on a finely sampled smooth profile.
  • Option B (optimize). Add an adaptive integrator that recursively bisects any subinterval whose Simpson error estimate $(S_{2}-S_{1})/15$ exceeds a tolerance, and count how many evaluations it saves versus a uniform Simpson rule on a profile with one sharp spike.
  • Option C (extend). Generalize trapz_data to two dimensions for a sampled image $I(x,y)$ — a double sum with half-weighted edges and quarter-weighted corners (the §22.5 product rule) — and compute a 2-D centroid $(\bar{x}, \bar{y})$.

Key Takeaways

  • Build quadrature once, as a reusable pure-function module, and validate it against known integrals by measuring the order-of-accuracy ratio before trusting it on real data.
  • Choose the rule by the constraints: Gauss for smooth analytic integrands (most accuracy per evaluation), the trapezoidal rule for sampled data (you cannot move the nodes).
  • Convergence-ratio checks ($4$ for $O(h^2)$, $16$ for $O(h^4)$) are cheap, decisive regression tests that catch weight and loop-bound bugs.
  • Let the physics sanity-check you: a symmetric profile's centroid must sit at the center, and a result that violates a known symmetry is a bug you were handed for free.