Trapezoidal Rule

The straight-line idea

You want a number for ∫_a^b f(x) dx but f is messy enough that no antiderivative is in reach. Approximate f by something easier — the simplest reasonable choice is the straight line L(x) that agrees with f at the two endpoints. That line is

L(x) = f(a) + (f(b) − f(a)) · (x − a) / (b − a)

The region between L and the x-axis on [a, b] is a trapezoid with parallel sides f(a) and f(b) and width h = b − a. Its area is

I ≈ (h / 2) · (f(a) + f(b))

That is the trapezoidal rule. One function evaluation at each endpoint, one multiplication, one addition. The shape of f between the endpoints is replaced wholesale by a line — the cheapest possible interpolation. Whatever the line misses is the truncation error, and bounding that error is the entire technical content of the rule.

The composite rule

Stretching a single trapezoid across a long interval is wasteful when f curves. The fix is to subdivide [a, b] into n equal strips of width h = (b − a)/n with nodes x_k = a + k·h for k = 0, 1, …, n. Apply the rule to each strip and add the results. Adjacent strips share an endpoint, so the interior nodes get weight 1 and the two boundary nodes get weight 1/2:

I ≈ h · [ f(x_0)/2 + f(x_1) + f(x_2) + … + f(x_{n-1}) + f(x_n)/2 ]

That is the composite trapezoidal rule. The cost is n + 1 evaluations of f. In NumPy: np.trapz(f(x), x) for any array x. In SciPy: scipy.integrate.trapezoid. Every introductory engineering course shows the formula, every CFD code has a function that computes it.

Why the error is O(h²)

On a single strip [x_k, x_k + h], interpolation theory gives

f(x) − L(x) = ½ · (x − x_k) · (x − x_k − h) · f''(ξ)        for some ξ in the strip

Integrating from x_k to x_k + h, the polynomial factor integrates to −h³/6, and the f''(ξ) factor is bounded by max|f''| on the strip. So

strip error = −(h³ / 12) · f''(ξ_k)

Summing across n strips and absorbing the average of f''(ξ_k) into a single representative f''(ξ) by the mean value theorem,

global error E_n = −(b − a) · h² · f''(ξ) / 12        for some ξ in [a, b]

That is the central error bound for the composite trapezoidal rule. It is second-order: the error scales as h², so halving the step size shrinks the error by a factor of four, and going from 10 to 100 strips drops the error by a factor of one hundred. The sign of the error follows the sign of f'' — for a convex integrand (f'' > 0) the trapezoidal rule overestimates the area, because the secant always lies above a convex curve.

Worked example: ∫₀¹ exp(−x²) dx

This is the integrand whose values define erf(1) — no elementary antiderivative. The true value is √π/2 · erf(1) ≈ 0.74682413281…. Try the trapezoidal rule with various n:

n       h         trapezoidal value    error
─       ──        ──────────────────   ────────────
1       1.000     0.683939720586       6.28 × 10⁻²
2       0.500     0.731370251828       1.55 × 10⁻²
4       0.250     0.742984097800       3.84 × 10⁻³
8       0.125     0.745865614845       9.59 × 10⁻⁴
16      0.0625    0.746584596788       2.40 × 10⁻⁴
32      0.03125   0.746764150605       5.98 × 10⁻⁵
64      0.015625  0.746809138114       1.50 × 10⁻⁵
128     0.007812  0.746820376113       3.76 × 10⁻⁶
1024    0.000977  0.746824073926       5.87 × 10⁻⁸

The pattern is exact: each doubling of n shrinks the error by close to a factor of four. From n = 32 to n = 64 the error drops from 5.98 × 10⁻⁵ to 1.50 × 10⁻⁵ — ratio 3.99. From n = 64 to n = 128, from 1.50 × 10⁻⁵ to 3.76 × 10⁻⁶ — ratio 3.99. That is the O(h²) signature. Reaching 10⁻⁸ accuracy demands roughly 1024 strips; Simpson's rule, by contrast, hits the same accuracy with just 32 strips because its error scales as h⁴.

For this integrand, the second derivative is f''(x) = (4x² − 2) e^{−x²}, with max value about 2 at x = 0. Plugging into the bound: |E_n| ≤ (1 − 0) · (1/n)² · 2 / 12 = 1/(6n²). At n = 32, the bound predicts |E| ≤ 0.000163, and the observed error 5.98 × 10⁻⁵ sits well below — typical, since f''(ξ) at some interior point is smaller than the maximum.

Trapezoidal vs Simpson vs Gauss

	Trapezoidal	Simpson's 1/3	Gauss-Legendre n-pt	Romberg (k levels)	Tanh-Sinh	Midpoint
Order of accuracy	2	4	2n (exact for poly degree ≤ 2n−1)	2k (Richardson-improved)	spectral / exponential	2
Error formula	−(b−a)h²·f''(ξ)/12	−(b−a)h⁴·f⁽⁴⁾(ξ)/180	varies; ∝ f⁽²ⁿ⁾	cancels low powers of h	doubly exponential	(b−a)h²·f''(ξ)/24
Function evals (n=8)	9	9	8	9 + extrapolation	~20 typical	8
Endpoint values needed	yes	yes	no (interior nodes)	yes	no	no
Periodic smooth f	exponentially convergent	fourth-order only	fourth-order on a period	same as trapezoidal	fast	exponentially convergent
Handles endpoint singularity	poorly	poorly	poorly	poorly	extremely well	moderately
Best when	simple, low accuracy	smooth f, modest accuracy	polynomial-like f	smooth, many digits needed	singular endpoints	can't evaluate at endpoints

The hierarchy: trapezoidal is the cheapest possible non-trivial quadrature. Simpson buys two more orders of accuracy at no extra evaluations (just a different weighting). Gauss-Legendre takes the same n evaluations and chooses both weights and node positions to maximise accuracy — for polynomial integrands of degree up to 2n − 1, Gauss is exact. Romberg uses trapezoidal at doubling refinements and Richardson-extrapolates the entries to cancel successive powers of h. SciPy's quad defaults to adaptive Gauss-Kronrod (a Gauss rule with extra nodes for error estimation), which is the modern production choice for general one-dimensional integration.

The miracle case: smooth periodic functions

For a smooth function periodic with period (b − a), the trapezoidal rule converges exponentially fast. The Euler-Maclaurin formula expands the error as

E_n = -Σ_{k=1}^∞ (B_{2k} / (2k)!) · h^{2k} · [f^{(2k-1)}(b) - f^{(2k-1)}(a)] + R

Each boundary term involves an odd derivative at the endpoints. If f is periodic, all those derivatives agree at a and b, so every boundary term cancels exactly. What remains is a remainder R that decays faster than any power of h — spectrally fast, like e^{−c/h}. Concretely: integrating cos²(x) over [0, 2π] by the trapezoidal rule with n = 4 nodes already gives machine precision. This is why discrete Fourier transforms work: trapezoidal sampling of e^{−2πikt} f(t) on [0, 1] is essentially exact for any smooth periodic f, which is the entire reason FFTs make signal processing tractable.

When endpoints cause trouble

Integrands with singularities or infinite derivatives at endpoints behave badly under trapezoidal. The bound E = −(b − a)h²f''(ξ)/12 assumes f'' is bounded — when f'' → ∞ near a or b, the constant explodes and convergence stalls. Classic example: ∫₀¹ √x dx. The true value is 2/3, and f'' = −¼ x^{−3/2} blows up at zero. The trapezoidal error decays only as O(h^{3/2}) here instead of O(h²).

Three fixes. Variable transformation: substitute x = t² to map the singularity onto a smooth integrand. Open rules: the midpoint rule and Gauss quadrature don't evaluate f at endpoints, dodging endpoint blowups directly. Tanh-sinh quadrature: a transformation x = ½ + ½·tanh(½π sinh t) packs nodes ultra-densely near endpoints, achieving doubly-exponential convergence for typical singularities. Mathematica and the mpmath library default to tanh-sinh for high-precision quadrature with endpoint singularities.

Romberg integration: trapezoidal squared

Apply the composite trapezoidal rule with n = 1, 2, 4, 8, 16, … strips, getting estimates T_0, T_1, T_2, …. Each successive T_k reuses all the function values from T_{k-1} plus n new points. By the Euler-Maclaurin formula, T_k has known error structure:

T_k = I + c_2·h_k² + c_4·h_k⁴ + c_6·h_k⁶ + ...

Richardson extrapolation eliminates the leading O(h²) term by combining T_{k-1} and T_k:

T_{k,1} = (4 T_k − T_{k-1}) / 3

This combination is exactly Simpson's rule on the same nodes. Iterating once more cancels the O(h⁴) term, giving order O(h⁶), and so on — each Romberg level cancels the next power of h. With k levels and 2^k + 1 trapezoidal evaluations you reach order O(h^{2k}). For smooth integrands, a 6-level Romberg table reaches machine precision in 65 function evaluations — what would take millions of strips with raw trapezoidal.

Where the trapezoidal rule shows up

• Finite-element matrix assembly. Computing element stiffness and mass matrices in FEM codes (ABAQUS, ANSYS, Code_Aster) requires integration of basis function products over each element. Linear-element matrices are computed analytically, but higher-order elements use composite trapezoidal or Gauss quadrature inside each element.
• Signal processing and DFT. Sampling a periodic signal at n equally spaced points and applying the trapezoidal rule is exactly the discrete Fourier transform up to a normalisation constant. The exponential convergence on periodic smooth functions is what makes spectral methods sharp.
• Image processing — cumulative distribution functions. Histogram equalisation, gamma correction lookup tables, and many tone-mapping algorithms compute CDFs by trapezoidal integration of pixel histograms. ImageMagick, OpenCV, and FFmpeg all do this.
• Financial pricing. Black-Scholes options pricing reduces to integrals of cumulative normal distributions. The trapezoidal rule applied to log-strike grids underlies the Carr-Madan FFT pricing method for European options.
• Pharmacokinetics — AUC. The area under the plasma-concentration-versus-time curve (AUC) is the standard pharmacokinetic exposure measure and is computed by the linear trapezoidal rule in every PK analysis package (WinNonlin/Phoenix, NONMEM, Pumas). FDA bioequivalence guidance specifically prescribes the trapezoidal AUC.
• GPS and motion sensors. Inertial navigation integrates accelerometer readings to velocity and position. Strap-down INS systems use composite trapezoidal or low-order Runge-Kutta on accelerometer time series at hundreds of hertz.

Variants and extensions

• Midpoint rule. Evaluate f only at strip centres x_k + h/2 with weight h: I ≈ h Σ f(x_k + h/2). Same O(h²) accuracy but error is +h²·f''(ξ)/24 — opposite sign and half the magnitude of trapezoidal. Robust to endpoint singularities since it never evaluates at the endpoints.
• Simpson's 1/3 rule. Fit a quadratic through three equally spaced points and integrate it: I ≈ (h/3)(f(x_0) + 4f(x_1) + f(x_2)). Fourth-order: error scales as h⁴. The next member of the Newton-Cotes family after trapezoidal.
• Simpson's 3/8 rule. Fit a cubic through four equally spaced points. Same fourth order as Simpson's 1/3 but uses three subintervals per fit; useful when the number of strips is divisible by 3.
• Bode's rule. Fit a quartic through five points. Sixth-order. Higher Newton-Cotes rules quickly suffer from negative weights and Runge phenomenon — degree 8 is the practical ceiling.
• Gauss-Legendre quadrature. Choose nodes and weights to maximise polynomial accuracy: an n-point rule is exact for polynomials of degree up to 2n − 1. Dominant in finite-element matrix assembly.
• Romberg integration. Trapezoidal at doubling refinements + Richardson extrapolation. Achieves order 2k after k extrapolation levels. Standard in textbooks; less common in production code than Gauss-Kronrod.
• Adaptive quadrature. Apply the trapezoidal (or Simpson, or Gauss) rule on the whole interval, then on each half, compare results, and recursively subdivide subintervals where the error estimate exceeds tolerance. SciPy's scipy.integrate.quad wraps QUADPACK, which uses adaptive Gauss-Kronrod.

Common pitfalls

• Forgetting the half-weight on endpoints. The composite formula gives endpoints weight 1/2, not 1. Treating them as 1 doubles their contribution and produces a consistent bias of order h.
• Refining when O(h²) won't reach your tolerance. Demanding 10⁻¹⁰ accuracy from raw trapezoidal on a non-periodic integrand means roughly 10⁵ strips. Switch to Simpson, Romberg, or Gauss before stamping out trapezoidal strips by the million.
• Ignoring endpoint singularities. A function with √x or log(x) at the lower endpoint slows convergence below O(h²). Use a transformation, an open rule, or tanh-sinh quadrature.
• Round-off saturation. At very small h, round-off in the sum f(x_0) + f(x_1) + … overwhelms truncation error. Use Kahan summation or compensated summation for n > 10⁶.
• Using trapezoidal in 2D and 3D as a default. Tensor-product trapezoidal in 2D needs O(n²) evaluations for O(h²) accuracy. In d dimensions, n^d evaluations becomes prohibitive — Monte Carlo, sparse grids, or quasi-Monte Carlo dominate for d ≥ 4.

Why the simplest rule still earns its keep

The trapezoidal rule is rarely the fastest quadrature on any specific problem. Simpson is more accurate per evaluation on smooth f, Gauss-Legendre dominates polynomial-like integrands, tanh-sinh sweeps up endpoint singularities. Yet every textbook starts here, every numerical library exposes a trapz, and every production code that integrates anything has a trapezoidal-sum loop somewhere — because the rule has three properties no other quadrature combines as cleanly: it is unconditionally stable, it requires nothing about f beyond function values at the chosen points, and the error is given in closed form by a single bound that any engineer can write down from memory. For uncertain integrands of modest accuracy needs, that is the right starting point. Everything fancier was built by improving on it.