12. Numerical integration

Some integrals have no elementary antiderivative — \int e^{-x^2}dx from §4.10's error function is one. Others have an antiderivative you simply don't want to find by hand. Either way, you need the number, not the formula, and §4.1's Riemann sums already showed the crude approach: chop the interval, sum up rectangles. This lesson makes that idea precise and then makes it fast.

The trapezoid rule

Riemann sums approximate each strip with a flat-topped rectangle. A trapezoid does better — it approximates each strip with a straight line between the two endpoint values, which tracks a curving function more closely.

\int_a^bf(x)\,dx\approx\frac{h}{2}\Big[f(x_0)+2f(x_1)+2f(x_2)+\cdots+2f(x_{n-1})+f(x_n)\Big]

where h=(b-a)/n and x_i=a+ih. Every interior point counts twice — once as the right end of one trapezoid, once as the left end of the next — and the two endpoints count once each.

Error. For a function with bounded second derivative M on [a,b],

|E_{\text{trap}}|\le\frac{(b-a)^3}{12n^2}M

Halving h (doubling n) cuts the error by a factor of about 4 — the error scales as h^2. That is a direct consequence of the second derivative appearing in the bound: a trapezoid is exact for straight lines, so the error comes entirely from curvature, and h^2 is how fast a Taylor remainder with a second-derivative term shrinks (§8.1).

The midpoint rule

Instead of connecting endpoints, evaluate at the middle of each strip and use that as the rectangle's height:

\int_a^bf(x)\,dx\approx h\sum_{i=0}^{n-1}f\!\left(a+\left(i+\tfrac12\right)h\right)

This is the rule you've been using in every "Doing it in Python" numerical check since §4.2 — it's simple, and it turns out to already carry the same h^2 error order as the trapezoid rule, with a bound half the size:

|E_{\text{mid}}|\le\frac{(b-a)^3}{24n^2}M

Same order, smaller constant, opposite sign of error. A convex function sits above its tangent lines, so the midpoint rectangle underestimates while the trapezoid — cutting the corner off that convexity — overestimates. That opposite-sign relationship is what the next rule exploits.

Simpson's rule

If the trapezoid and midpoint errors have opposite signs, a weighted average of the two should cancel most of the error. Simpson's rule is exactly that average (in a 1:2 ratio), and it's equivalent to fitting a parabola through each group of three points instead of a line through two:

\int_a^bf(x)\,dx\approx\frac h3\Big[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\cdots+4f(x_{n-1})+f(x_n)\Big]

(requires n even) — weights alternate 4,2,4,2,\ldots,4 between the fixed endpoint weights of 1.

Error. The cancellation is not partial — it's a full order of h better:

|E_{\text{Simp}}|\le\frac{(b-a)^5}{180n^4}M_4

where M_4 bounds the fourth derivative. Halving h now cuts the error by roughly a factor of 16. A parabola is exact for cubics too (the odd-degree term integrates to zero by symmetry within each panel), which is why the order jumps by two rather than one.

Why the order jump matters

Trapezoid and midpoint are O(h^2): to get one more correct digit you need roughly \sqrt{10}\approx3.16\times as many points. Simpson is O(h^4): one more digit needs only 10^{1/4}\approx1.78\times as many points. For a smooth function, Simpson's rule with n=32 routinely beats the trapezoid rule with n=10{,}000. This is the same lesson as §8's Taylor series and §3.2's Newton's method: matching more derivatives at each step buys a higher order of convergence, and higher order compounds fast.

Doing it in Python

All three rules on \int_0^1e^x\,dx=e-1, watching the error shrink as n grows:

from math import exp, e

def trapezoid(f, a, b, n):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    total += sum(f(a + i * h) for i in range(1, n))
    return total * h

def midpoint(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + (i + 0.5) * h) for i in range(n))

def simpson(f, a, b, n):
    h = (b - a) / n
    total = f(a) + f(b)
    total += 4 * sum(f(a + i * h) for i in range(1, n, 2))
    total += 2 * sum(f(a + i * h) for i in range(2, n, 2))
    return total * h / 3

exact = e - 1
print(f"{'n':>5} {'trapezoid err':>16} {'midpoint err':>16} {'simpson err':>16}")
for n in (2, 4, 8, 16, 32):
    et = abs(trapezoid(exp, 0, 1, n) - exact)
    em = abs(midpoint(exp, 0, 1, n) - exact)
    es = abs(simpson(exp, 0, 1, n) - exact)
    print(f"{n:>5} {et:>16.3e} {em:>16.3e} {es:>16.3e}")

Doubling n each row: trapezoid and midpoint errors drop by roughly 4\times, Simpson's by roughly 16\times — exactly the O(h^2) versus O(h^4) the error bounds predicted.

An integral with no elementary antiderivative, computed to high precision:

from math import exp

def simpson(f, a, b, n):
    h = (b - a) / n
    total = f(a) + f(b)
    total += 4 * sum(f(a + i * h) for i in range(1, n, 2))
    total += 2 * sum(f(a + i * h) for i in range(2, n, 2))
    return total * h / 3

f = lambda x: exp(-x**2)
approx = simpson(f, 0, 1, 100)
print(f"int_0^1 e^(-x^2) dx = {approx:.10f}")
print("no elementary antiderivative -- Simpson gets 10 digits from n=100")

scipy.integrate.quad, the production-grade version of this idea (it uses an adaptive scheme related to Simpson's rule, refining wherever the error estimate is largest):

from scipy import integrate
from math import exp

result, error_estimate = integrate.quad(lambda x: exp(-x**2), 0, 1)
print(f"quad result         : {result:.10f}")
print(f"quad's own estimate : {error_estimate:.2e}")
print("this is what runs under the hood whenever a real system")
print("needs a numerical integral")

Worked example

Estimate \displaystyle\int_0^{\pi/2}\sqrt{1-0.5\sin^2x}\,dx — an elliptic integral (§9's ellipse arc-length family) with no elementary antiderivative — using Simpson's rule with n=4, and bound the error using n=8 as a reference.

h=\frac{\pi/2}{4}=\frac{\pi}{8}. Nodes: x_0=0,\,x_1=\frac\pi8,\,x_2=\frac\pi4,\,x_3=\frac{3\pi}8,\,x_4=\frac\pi2.

f(x)=\sqrt{1-0.5\sin^2x}

f(x_0)=1,\quad f(x_1)\approx0.962692,\quad f(x_2)\approx0.866025,\quad f(x_3)\approx0.757115,\quad f(x_4)\approx0.707107

\int\approx\frac h3\Big[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+f(x_4)\Big]

=\frac{\pi/8}{3}\Big[1+4(0.962692)+2(0.866025)+4(0.757115)+0.707107\Big]

=0.1308997\times10.318385\approx\boxed{1.35067}

Sanity check. The integrand is between \sqrt{0.5}\approx0.7071 and 1 across an interval of length \pi/2\approx1.5708, so the true value must lie between 1.1107 and 1.5708. 1.35067 sits comfortably inside that range. ✓

Now the promised error estimate, from n=8. Rerun with h=\frac{\pi}{16}:

S_8\approx1.3506439

The two estimates differ by |S_8-S_4|\approx3.0\times10^{-5}. Because Simpson is O(h^4), halving h should shrink the error about 16-fold, so almost all of that gap is S_4's error, and S_8's own error is roughly \frac1{15} of it — the standard Richardson argument:

|S_4-\text{true}|\approx|S_8-S_4|=3.0\times10^{-5},\qquad |S_8-\text{true}|\lesssim\frac{|S_8-S_4|}{15}=2.0\times10^{-6}

Checking against the true value E(0.5)=1.3506439: S_4 is off by 3.0\times10^{-5} — exactly what the gap predicted — and S_8 is off by 8.7\times10^{-9}, comfortably inside its 2\times10^{-6} estimate. The S_8 figure is very conservative here, because this integrand has vanishing derivatives at both endpoints and so converges faster than the generic O(h^4); an error estimate erring on the safe side is the right failure direction.

The point is the method, not the numbers. Two runs at different n give you an error bar without ever knowing the answer — which is the only kind of error estimate available when there is no antiderivative to compare against.

Your turn

1. Using the trapezoid rule with n=2, estimate \displaystyle\int_0^2x^2\,dx (exact value 8/3).

2. Using the midpoint rule with n=2 on the same integral, estimate it and compare the sign of the error to part 1.

3. True or false: doubling n in Simpson's rule is guaranteed to reduce the error by a factor of exactly 16.

Solutions

1. h=1, nodes x_0=0,x_1=1,x_2=2, f(x)=x^2:

\frac12\Big[f(0)+2f(1)+f(2)\Big]=\frac12\big[0+2+4\big]=3

Error: 3-\frac83=\frac13\approx0.333, an overestimate — x^2 is convex, so the trapezoid's straight-line tops sit above the curve.

2. h=1, midpoints 0.5 and 1.5:

1\cdot\big[f(0.5)+f(1.5)\big]=0.25+2.25=2.5

Error: 2.5-\frac83\approx-0.167, an underestimate, opposite sign from part 1 and about half the magnitude — exactly the relationship the error bounds predicted, and why Simpson's weighted blend of the two cancels so much error.

3. False. The bound guarantees the error shrinks by a factor of up to 16 for smooth functions with a well-behaved fourth derivative — it's an asymptotic rate, not an exact multiplier for any particular n. Real error ratios cluster near 16 for smooth f but can differ, especially at small n or when f^{(4)} varies a lot across the interval.

Check yourself in code

Compute \int_0^1e^x\,dx with the trapezoid, midpoint, and Simpson rules for n=2,4,8,16,32, printing the absolute error against the exact value e-1 for each.

Print exactly this:

    n    trapezoid err     midpoint err      simpson err
    2        3.565e-02        1.777e-02        5.793e-04
    4        8.940e-03        4.467e-03        3.701e-05
    8        2.237e-03        1.118e-03        2.326e-06
   16        5.593e-04        2.796e-04        1.456e-07
   32        1.398e-04        6.992e-05        9.103e-09

from math import exp, e

def trapezoid(f, a, b, n):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    total += sum(f(a + i * h) for i in range(1, n))
    return total * h

def midpoint(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + (i + 0.5) * h) for i in range(n))

def simpson(f, a, b, n):
    h = (b - a) / n
    total = f(a) + f(b)
    total += 4 * sum(f(a + i * h) for i in range(1, n, 2))
    total += 2 * sum(f(a + i * h) for i in range(2, n, 2))
    return total * h / 3

exact = e - 1
print(f"{'n':>5} {'trapezoid err':>16} {'midpoint err':>16} {'simpson err':>16}")
# for each n, print n, trapezoid error, midpoint error, simpson error

from math import exp, e

def trapezoid(f, a, b, n):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    total += sum(f(a + i * h) for i in range(1, n))
    return total * h

def midpoint(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + (i + 0.5) * h) for i in range(n))

def simpson(f, a, b, n):
    h = (b - a) / n
    total = f(a) + f(b)
    total += 4 * sum(f(a + i * h) for i in range(1, n, 2))
    total += 2 * sum(f(a + i * h) for i in range(2, n, 2))
    return total * h / 3

exact = e - 1
print(f"{'n':>5} {'trapezoid err':>16} {'midpoint err':>16} {'simpson err':>16}")
for n in (2, 4, 8, 16, 32):
    et = abs(trapezoid(exp, 0, 1, n) - exact)
    em = abs(midpoint(exp, 0, 1, n) - exact)
    es = abs(simpson(exp, 0, 1, n) - exact)
    print(f"{n:>5} {et:>16.3e} {em:>16.3e} {es:>16.3e}")

---

The trapezoid rule fits a line through each pair of points, the midpoint rule evaluates at each strip's center, and both carry O(h^2) error — halving h quarters the error. Simpson's rule averages the two, equivalent to fitting a parabola through each triple, and its errors cancel to one full order better: O(h^4), where halving h shrinks the error sixteenfold. That order jump is why Simpson's rule with a few dozen points can outperform crude sums with thousands, and it's the same "match more derivatives, buy a higher order" trade you saw in §3.2's Newton's method and §8's Taylor series.

Next: with three techniques and three numerical rules on the table, how do you pick?