Trapezoidal rule for numerical integration

The trapezoidal rule is the simplest method for approximating definite integrals. Based on approximating f with line segments between consecutive points, forming trapezoids whose total area estimates the integral. Known since antiquity (Newton-Cotes 1, already in Heron and Archimedes).

Basic formula

For interval [a,b] with two points a and b:

integral_{a}^{b} f(x) dx ~ (b-a)/2 * [f(a) + f(b)]

Exact for linear functions. Error: -(b-a)^3/12 * f''(xi) for some xi in (a,b).

Composite trapezoidal

Divides [a,b] in n sub-intervals of width h = (b-a)/n and applies trapezoidal to each:

integral ~ h [f(x_0)/2 + f(x_1) + f(x_2) + ... + f(x_{n-1}) + f(x_n)/2]

where x_i = a + i h. Global error: -(b-a) h^2/12 * max|f''|. Order 2 convergence (error /4 halving h).

Linear polynomial exactness

Trapezoidal is EXACT for linear f (ax + b). For quadratic f = x^2 on [0,1] (exact = 1/3): trapezoidal n=10 gives 0.335 (error 1.5e-3); n=100 gives 0.33335 (error 1.7e-5). For polynomials degree >= 2 error proportional to max|f''|.

Richardson extrapolation

Idea: if T_h is trapezoidal approximation with step h, then true value I is given by T_h - C h^2 - C2 h^4 - ... Combining T_h and T_{h/2} eliminates C h^2:

I ~ (4 T_{h/2} - T_h) / 3 (= refined Simpson 1/3)

Iterating eliminates successive orders. Result: Romberg table.

Romberg table

Romberg algorithm builds matrix R[k][j]:

• R[k][0] = trapezoidal with 2^k intervals.
• R[k][j] = R[k][j-1] + (R[k][j-1] - R[k-1][j-1]) / (4^j - 1).

R[k][k] (diagonal) has order 2(k+1). For k=4 order 10! Exponential convergence for regular (analytic) functions.

Numerical examples

Integral of sin(x) on [0, pi] (exact = 2):

• Trapezoidal n=10: 1.9835 (error 1.6e-2)
• Trapezoidal n=20: 1.9959 (error 4.1e-3, /4 of above)
• Trapezoidal n=40: 1.9990 (error 1.0e-3, /4)
• Romberg 4 levels: 1.99999992 (error 8e-8)
• Romberg 6 levels: machine precision

Trapezoidal advantages

• Very simple implementation, closed formula.
• Cost: 1 f evaluation per point (n+1 total).
• Exact for linear (useful for adaptive grid).
• Numerically stable (formula with positive coefficients).
• Basis for more sophisticated methods (Romberg, Simpson).

Disadvantages

• Slow convergence: order 2 (error /4 halving h).
• Inefficient for extreme accuracy: needs huge n.
• Overestimates for convex functions (positive curvature).
• Underestimates for concave functions (negative curvature).
• For highly oscillating integrals: needs adaptive sampling.

Trapezoidal vs Simpson

Trapezoidal (order 2) and Simpson (order 4) have similar costs but very different accuracies. For same n:

• Trapezoidal: error O(h^2), approximating polynomial degree 1.
• Simpson 1/3: error O(h^4), approximating polynomial degree 2.
• For smooth functions, Simpson typically 100-1000x more accurate.
• For discontinuous functions or few points, trapezoidal may be preferable (less overfitting).

Adaptive quadrature

Strategy: divide [a,b] in sub-intervals and apply trapezoidal (or Simpson) until estimated error < tolerance. Error estimate via difference T_h and T_{h/2}: error_h ~ T_h - T_{h/2}. Result: small steps where f varies rapidly, large steps where flat.

History

Known to ancients (Babylonians, Archimedes). Newton-Cotes 1 (closed formula with 2 points) is simplest case. Roger Cotes 1722 systematizes. Werner Romberg 1955 introduces systematic Richardson extrapolation. Today basis of professional libraries (scipy.integrate.romb, MATLAB trapz, numpy.trapz).

The trapezoidal error and the Euler-Maclaurin formula

The error of the trapezoidal rule has a precise form: for a single interval it is -(b-a)^3/12 f''(xi) for some interior point xi, and in the composite version it becomes -(b-a) h^2/12 f''(xi), hence order two. The sign of the second derivative determines the direction of the error: for convex functions (f'' > 0) the trapezoids systematically overestimate the integral, because the chord lies above the curve; for concave functions they underestimate it.

The Euler-Maclaurin formula gives the full expansion of the error in even powers of h, with coefficients given by the Bernoulli numbers and the values of the derivatives at the interval endpoints. It is this structure, with only even terms, that makes Richardson extrapolation so effective: eliminating one term at a time yields the Romberg table, which rapidly derives from the trapezoidal rule approximations of ever higher order.

The surprise: spectral convergence on periodic functions

The trapezoidal rule has a surprising property: applied to a smooth periodic function, integrated over a full period, it converges exponentially instead of as h^2. The reason emerges from the Euler-Maclaurin formula: for a periodic integrand the values of the derivatives at the two endpoints coincide and all the correction terms cancel in pairs, leaving an error that decreases faster than any power of h.

This spectral convergence explains why the trapezoidal rule, apparently the crudest method, is in fact the optimal tool for Fourier analysis and the discrete transform (DFT), where signals are periodic by construction. The same idea makes the trapezoidal rule the method of choice for contour integrals in the complex plane, where the closed path is naturally periodic. An elegant example of how an elementary method can excel in the right context.

Updated July 2026