Chapter 22 Exercises: Numerical Integration and Differentiation

These problems build the two reflexes this chapter exists to instill: name a method's order before you trust it, and measure that order on your own machine. Work them with a compiler open — predict every output first, then compile and check. When your number matches the hand-computed one, you have understood the method; when it does not, you have found something worth understanding.

Difficulty tiers. ⭐ warm-up (minutes), ⭐⭐ standard (the core skill), ⭐⭐⭐ challenge (synthesis or open-ended). Solutions to the †-marked and all odd-numbered problems are in appendices/answers-to-selected.md; the computational ones are also provided as compilable, hand-checked code in code/exercise-solutions.f90. Never run code to find an answer — reason it out, then verify.

Part A — Type, Compile, and Run

Predict the printed output before compiling.

22.1 † ⭐ A forward-difference program approximates $f'(3)$ for $f(x) = x^2$ (exact answer $6$) with $h = 0.5$:

fwd = (sq(3.0_dp + 0.5_dp) - sq(3.0_dp)) / 0.5_dp
print '(a, f8.4)', 'forward estimate = ', fwd

What does it print, and what is the error? From the Taylor expansion, what is the exact truncation error of the forward difference for a function with constant $f''$, and does the printed error match it?

22.2 ⭐ Predict the output of the trapezoidal rule applied to $\int_0^2 x\,dx$ (exact answer $2$) with $n = 2$ panels. Why is the trapezoidal rule exact for this integrand, and for what class of functions is that always true?

22.3 † ⭐⭐ Simpson's rule is applied to $\int_0^1 x^3\,dx$ (exact answer $\tfrac14$) with $n = 2$ panels. Predict the printed value to eight decimals, and explain in one sentence why Simpson returns the exact answer despite $x^3$ being a cubic, not a quadratic.

22.4 ⭐⭐ The program in code/exercise-solutions.f90 (routine ex04_convergence) prints a convergence table for the central difference of $f(x) = x^4$ at $x = 1$ (exact $f'(1) = 4$) on grids $h = 0.1, 0.05, 0.025$. Predict the three error values and the two ratios, and state the order of accuracy you would read off.

Part B — By Hand

Paper and pencil; no compiler needed.

22.5 † ⭐ Starting from the Taylor expansion of $f(x-h)$ about $x$, derive the truncation error of the backward difference $\big(f(x) - f(x-h)\big)/h$. Show its leading term is $-\tfrac{h}{2}f''(x)$, hence $O(h)$.

22.6 ⭐⭐ Add the Taylor expansions of $f(x+h)$ and $f(x-h)$ and rearrange to derive the three-point second difference $\big(f(x+h) - 2f(x) + f(x-h)\big)/h^2 \approx f''(x)$. Show its leading error term is $\tfrac{h^2}{12}f^{(4)}(x)$, and state the order.

22.7 † ⭐⭐ Apply Simpson's rule with $n = 2$ to $\int_0^1 e^x\,dx$ by hand, using the rounded values $e^0 = 1$, $e^{0.5} = 1.64872$, $e^1 = 2.71828$. Compare with the exact value $e - 1 = 1.71828$. Is the error consistent with an $O(h^4)$ method?

22.8 ⭐⭐⭐ Apply two-point Gauss-Legendre to the same $\int_0^1 e^x\,dx$, using the mapped nodes $x \approx 0.21132$ and $0.78868$ with $e^{0.21132} = 1.23533$ and $e^{0.78868} = 2.20050$. Compare your result with Simpson's from 22.7 and with the exact value. Two function evaluations versus Simpson's three — which is more accurate here, and why is that the expected outcome?

Part C — Find the Bug

Each snippet compiles but is wrong (or is right code used wrongly). Diagnose and fix it.

22.9 † ⭐⭐ This Simpson routine gives the wrong answer for every non-trivial integrand. The weights are misapplied — find the error.

s = f(a) + f(b)
do i = 1, n - 1
  if (mod(i, 2) == 0) then
    s = s + 4.0_dp * f(a + i*h)     ! even node
  else
    s = s + 2.0_dp * f(a + i*h)     ! odd node
  end if
end do
s = s * h / 3.0_dp

22.10 ⭐⭐ A colleague reports that their central-difference derivative "gets worse the smaller I make h" and shows you a loop that drives $h$ down to 1.0e-14_dp. The code is otherwise correct. What is happening, at what $h$ should they have stopped, and which chapter explains it?

22.11 † ⭐⭐ This trapezoidal routine is systematically too large. Find the one missing factor.

s = f(a) + f(b)
do i = 1, n - 1
  s = s + f(a + i*h)
end do
s = s * h

Part D — Port It

Translate the given snippet to modern Fortran, and reason about relative speed.

22.12 ⭐⭐ Port this NumPy central-difference one-liner to a Fortran whole-array section expression (no explicit loop). Assume u(1:n) holds samples on a uniform grid of spacing h; produce the interior derivatives d(2:n-1).

d[1:-1] = (u[2:] - u[:-2]) / (2*h)     # central difference, interior points

22.13 † ⭐⭐ Port this pure-Python Simpson loop to a Fortran function, and estimate the speedup for $n = 10^{7}$ panels of a cheap integrand — is the win closer to $2\times$, $10\times$, or $50\times$, and what determines it?

s = f(a) + f(b)
for i in range(1, n):
    s += (4 if i % 2 else 2) * f(a + i*h)
s *= h/3

22.14 ⭐⭐⭐ MATLAB's trapz(x, y) integrates sampled data y at abscissae x with the trapezoidal rule, allowing non-uniform spacing. Write a Fortran function trapz(x, y) taking two assumed-shape arrays and returning $\sum_i \tfrac{1}{2}(x_{i+1}-x_i)(y_i + y_{i+1})$. Why does non-uniform spacing rule out the tidy "half the endpoints" formula?

Part E — Back of the Envelope

Estimates, not exact answers. Show your reasoning.

22.15 † ⭐⭐ You refine a quadrature rule by halving $h$ three times and measure error ratios of $3.9, 3.98, 4.0$. What is the method's order of accuracy, and which of the chapter's rules is it? If instead the ratios were $15.8, 16.0, 16.1$, which rule?

22.16 ⭐⭐ Estimate the optimal step $h^{*}$ and the best achievable error for a central-difference derivative in double precision ($\varepsilon \approx 2.2\times10^{-16}$), given that truncation error is $\sim h^2$ and round-off error is $\sim \varepsilon/h$. Compare with the forward difference's $h^{*}\sim \sqrt\varepsilon$.

22.17 † ⭐⭐ To integrate a smooth function to a relative accuracy of $10^{-6}$, roughly how many panels does the trapezoidal rule ($O(h^2)$) need compared with Simpson's rule ($O(h^4)$)? Take the error constants as comparable and reason from the orders alone.

22.18 ⭐⭐⭐ A product Gauss rule uses $5$ nodes per dimension. How many function evaluations does it need in $10$ dimensions? A Monte Carlo estimate reaches error $\sim 0.1/\sqrt{N}$. Roughly how many samples $N$ would Monte Carlo need to match a $1\%$ accuracy, and how does that compare with the product rule's cost? What does this say about when to abandon grids?

Part F — Design It

Extend the chapter's code or the heat solver. Compile everything.

22.19 † ⭐⭐ Extend the quadrature module with a three-point Gauss-Legendre rule gauss3(f, a, b), using nodes $0, \pm\sqrt{3/5}$ and weights $\tfrac89, \tfrac59, \tfrac59$ on $[-1,1]$ (mapped to $[a,b]$). Test it on $\int_{-1}^{1} x^4\,dx = \tfrac25$, which it must return exactly (degree $4 \le 2n-1 = 5$). The solution is in code/exercise-solutions.f90 (ex19_gauss3).

22.20 ⭐⭐ Write a richardson(coarse, fine, p) function that combines two estimates of orders $p$ and returns the extrapolated value $\big(2^p\,\text{fine} - \text{coarse}\big)/(2^p - 1)$. Feed it the Simpson values $S_2 = 0.20833333$ and $S_4 = 0.20052083$ for $\int_0^1 x^4$ with $p = 4$; it should return the exact $0.2$. (ex20_richardson in the code.)

22.21 † ⭐⭐⭐ Design a reusable subroutine measure_order(op, x0, exact, h0, levels) that applies a difference operator op (a procedure argument) on a halving sequence of steps and prints the measured order at each level, generalizing the Project Checkpoint. What does it print for a correct $O(h^2)$ stencil, and what would a coding bug that made it accidentally $O(h)$ look like in the output?

22.22 ⭐⭐⭐ Add a diagnostic to the heat solver: a two-dimensional trapezoidal rule integrate_field(u, dx, dy) that returns $\int\!\!\int u\,dx\,dy$ over the plate (a double sum with half-weighted edges and quarter-weighted corners). Verify it on $u(x,y) = xy$ over the unit square, where the exact integral is $\tfrac14$. (ex22_field_integral in the code.)

Part G — Interleaved (earlier chapters)

22.23 † ⭐ (Ch. 6) The trapezoid and simpson functions of §22.2 are declared pure and take the integrand as a procedure(scalar_fn) argument. Why must the abstract interface scalar_fn itself be declared pure function? What compiler error would you get if the actual integrand you passed were not pure?

22.24 ⭐⭐ (Ch. 20) Numerical differentiation has a round-off floor (a smallest useful $h$), but numerical integration does not — you can keep adding panels and keep improving. Explain the asymmetry in terms of the operations each performs on nearly equal floating-point values.

22.25 † ⭐⭐ (Ch. 21) The three-point second difference, applied at every interior node of a 1-D grid, assembles into a tridiagonal matrix with $-2$ on the diagonal and $1$ on the off-diagonals (times $1/h^2$). Which LAPACK routine from Chapter 21 would solve a system with this matrix in $O(n)$ time, and why is the dense dgesv the wrong tool for it?

22.26 ⭐⭐ (Ch. 5) Write a single whole-array expression that computes the central-difference first derivative of the interior of real(dp) :: u(:) on a uniform grid of spacing h, using array sections. Why does this hand the compiler more freedom to vectorize than an explicit do loop would?

Check yourself: the two skills you should own after these problems are (1) reading a method's order off its Taylor expansion and confirming it with a halving experiment, and (2) knowing that a finite-difference derivative has a best step near $\sqrt\varepsilon$ — smaller is worse. If either felt shaky, re-read §22.4 and redo 22.4, 22.15, and 22.16.