The Trapezoidal Rule, Worked

Estimate an integral by summing trapezoids under the curve, and understand its second-order error and one surprising strength.

The rule

Approximate the area under f between sample points by trapezoids. With n equal intervals of width h, the integral is approximately h[(f0 + fn)/2 + f1 + f2 + ... + f(n-1)]: interior points full weight, endpoints half.

import numpy as np
def trap(f,a,b,n):
    x=np.linspace(a,b,n+1); h=(b-a)/n; y=f(x)
    return h*(y[0]/2+y[-1]/2+y[1:-1].sum())
for n in [4,8,16,32]:
    print(n, round(abs(trap(np.sin,0,np.pi,n)-2),6))  # error ~ 1/n^2

Error

The trapezoidal error scales as h^2 and is proportional to the second derivative of f. Concave-down functions are underestimated, concave-up overestimated - a useful intuition for bounding the error.

The periodic surprise

For a smooth periodic function integrated over a full period, the trapezoidal rule is spectacularly accurate - error falls faster than any power of h. This is why it is the natural quadrature for Fourier coefficients and for integrals around closed contours.

Richardson extrapolation

Combining trapezoidal estimates at h and h/2 to cancel the leading h^2 error yields Simpson's rule; iterating this idea produces Romberg integration, which climbs to high order automatically. The humble trapezoid is thus the foundation of a whole family of accurate quadrature methods.