Trapezoidal Rule

1. Basic Trapezoidal Rule

We approximate the definite integral of a function $f(x)$ over an interval $[a,b]$ by the area of the trapezoid bounded by the straight line joining the two endpoints.

This amounts to interpolating $f$ on $[a,b]$ by a degree-1 polynomial (a straight line) and taking the area beneath it. It has first-order accuracy and is exact when $f$ is linear.

In the animation below we divide $[0,1]$ into $n$ equal parts, approximate each subinterval by a trapezoid, and write the resulting sum as the approximation $T_n$ (here $n$ is the number of trapezoids the interval is split into). The larger $n$ is (the finer the subdivision), the closer $T_n$ gets to the true value of the integral. The exact formula for $T_n$ is given in Section 2 (the composite trapezoidal rule).

2. Composite Trapezoidal Rule

Divide $[a,b]$ into $n$ equal parts and sum the trapezoidal rule applied to each subinterval. With step size $h = (b-a)/n$ and nodes $x_k = a + kh$,

The two endpoints are summed with weight $1/2$ and the interior nodes with weight $1$.

3. Error Estimate

Halving the step size $h$ reduces the error by roughly $1/4$. This property underlies Richardson extrapolation and Romberg integration (see also Section 4).

4. Endpoint Correction

By the Euler-Maclaurin formula, the error of the composite trapezoidal rule can be expanded as a power series in $h$:

Here the $B_{2k}$ are the Bernoulli numbers ($B_2 = \tfrac{1}{6},\ B_4 = -\tfrac{1}{30},\ \ldots$). If the endpoint derivative values $f^{(2k-1)}(a),\,f^{(2k-1)}(b)$ are known, one can subtract the correction terms of this expansion from $T_n$ to raise the accuracy. This is called endpoint correction (gradient correction).

First-order endpoint correction

Use only the $k=1$ term. Since $\dfrac{B_2}{2!} = \dfrac{1}{12}$, the leading error term is $\dfrac{h^2}{12}\bigl[f'(b) - f'(a)\bigr]$. Subtracting it from $T_n$,

the $h^2$ term cancels and the error of $T_n^{\text{corr}}$ improves to $O(h^4)$. The advantage is that, knowing only the first derivatives $f'(a),\,f'(b)$ at the endpoints, the accuracy gains two orders without any additional function evaluations.

Second-order endpoint correction

Correcting further, up to the $k=2$ term (which uses $f^{(3)}$), gives $O(h^6)$. Since $\dfrac{B_4}{4!} = \dfrac{-1/30}{24} = -\dfrac{1}{720}$, the sign of the correction term becomes positive, and

Both the $h^2$ and $h^4$ terms cancel, so the error of $T_n^{\text{corr2}}$ is $O(h^6)$. In general, the larger $k$ is the higher the endpoint derivatives $f^{(2k-1)}(a),\,f^{(2k-1)}(b)$ required, but each correction raises the accuracy by two more orders.

Worked example: the effect of the correction

Take $\displaystyle\int_0^1 e^x\,dx$. For $f(x)=e^x$ we have $f'(x)=f^{(3)}(x)=e^x$, so $f'(1)-f'(0) = f^{(3)}(1)-f^{(3)}(0) = e - 1$. We correct the value $T_4 = 1.727222$ from Figure 1 ($n=4$, error $8.94\times10^{-3}$), taking $h=\tfrac14$.

First-order correction ($O(h^4)$):

The error from the true value $e-1 = 1.718282$ is $9.3\times10^{-6}$, about $1/1000$ of the uncorrected error.

Second-order correction ($O(h^6)$): adding the $f^{(3)}$ term, with $\dfrac{h^4}{720}(e-1) = 0.0000093$,

The error drops to $1.4\times10^{-8}$, a further reduction of about $1/700$ over the first-order correction (roughly $1/600{,}000$ of the uncorrected error).

Adding one endpoint derivative at a time improves the accuracy by orders of magnitude without increasing $n$.

5. Super-convergence for Periodic Functions

If $f$ is a smooth periodic function of period $T = b - a$, then $f^{(k)}(a) = f^{(k)}(b)$ holds for every $k$, so all the correction terms of the Euler-Maclaurin expansion vanish. As a result, the trapezoidal rule converges exponentially.

This is extremely important in practice: for integrating periodic functions the trapezoidal rule is one of the best methods available.

Example 1: trigonometric polynomials are integrated exactly

$\sin^2 x = \dfrac{1 - \cos 2x}{2}$ is a trigonometric polynomial of frequency up to 2. Its true value over one period $[0,2\pi]$ is

The composite trapezoidal rule ($h = 2\pi/n$) gives the following values.

For $n \ge 3$ it equals the true value $\pi$ exactly. This is because, over one full period, the trapezoidal rule integrates a trigonometric polynomial with zero error once the subdivision count $n$ exceeds the highest frequency of the integrand (here the frequency is 2, so $n \ge 3$). For $n = 1, 2$ the frequency-2 component is aliased onto the mean, producing an error. Likewise the simplest case $\displaystyle\int_0^{2\pi}\sin x\,dx = 0$ is obtained exactly as $0$ for every $n$.

This exactness, however, is a privilege of integrating over one complete period. Computing the half-period integral $\displaystyle\int_0^{\pi}\sin^2 x\,dx = \dfrac{\pi}{2}$ with the trapezoidal rule does not benefit from the periodicity and reverts to the usual $O(h^2)$ error.

Example 2: exponential convergence for a non-polynomial integrand

If the integrand is not a trigonometric polynomial the result is never exact, but for a smooth periodic function it still converges to the true value exponentially. The Fourier series of $e^{\cos x}$ is infinite, so it is never exact for finite $n$; the true value can be written using the modified Bessel function of the first kind $I_0$ as

The composite trapezoidal rule gives the following values $T_n$ and errors.

Simply going from $n = 4$ to $n = 8$ cuts the error from $3.4\times10^{-2}$ to $1.3\times10^{-6}$, a factor of about $1/20{,}000$. This far exceeds $O(h^2)$ (where doubling $n$ reduces the error only by $1/4$): the error decreases geometrically in $n$, which is super-convergence. The smoother and more periodic the integrand, the more the Euler-Maclaurin correction terms vanish one after another, giving this fast convergence.

6. Singular Endpoints

The composite trapezoidal rule $T_n = h\left[\dfrac{f(a)}{2} + \displaystyle\sum_{k=1}^{n-1} f(x_k) + \dfrac{f(b)}{2}\right]$ contains the endpoint values $f(a),\,f(b)$ directly. So when $f$ is undefined at an endpoint (its value does not exist, or it diverges), the rule cannot be applied as is. This is a question about improper integrals; we organize it into three cases.

Case 1: the limit at the endpoint is finite (a removable singularity)

For example, $f(x)=\dfrac{\sin x}{x}$ is undefined at $x=0$, yet $\displaystyle\lim_{x\to 0}\frac{\sin x}{x}=1$. For such a removable singularity, redefining $f(a)$ by the limit value ($f(0):=1$) makes the function continuous and the trapezoidal rule applies directly. Once $f$ is again twice continuously differentiable up to the endpoints, the $O(h^2)$ error estimate holds as well.

Case 2: the function diverges at the endpoint but the integral converges

Example: $\displaystyle\int_0^1 \frac{1}{\sqrt{x}}\,dx = 2$ converges, but $f(0)=+\infty$, so $f(a)$ cannot be written down. This is the situation the trapezoidal rule handles worst; the following remedies are available.

• Use an open quadrature rule: switch to a rule that does not evaluate the endpoints, such as the midpoint rule, so the singular point is never touched.
• Remove the singularity by a change of variable: a suitable substitution makes the integrand smooth. For endpoint singularities the double-exponential (DE, tanh-sinh) formula of Takahashi and Mori is especially powerful and is widely used in practice.
• Separate the singular part analytically: factor out the singular factor, as in $f(x)=\dfrac{g(x)}{\sqrt{x}}$, integrate it by hand, and apply numerical integration only to the remaining smooth part.
• Use weighted quadrature: absorb the endpoint singularity into a weight function, as in Gauss-Jacobi quadrature.

Case 3: the integral itself diverges

Example: $\displaystyle\int_0^1 \frac{1}{x}\,dx = \infty$. Since the value of the integral does not exist, no numerical method can produce a correct finite value (the value keeps growing as the mesh is refined). Here the rule is used to detect the divergence numerically.

Summary of the three cases

Cases 1 to 3 above can be summarized as follows.

The key point is that the $O(h^2)$ error of the trapezoidal rule assumes $f \in C^2[a,b]$ (twice continuously differentiable up to the endpoints). If $f$ or its derivatives blow up at an endpoint, then even when a value can be filled in and the rule computed, convergence becomes extremely slow and the accuracy guarantee is lost. The basic strategy for endpoint singularities is therefore to make the integrand smooth first, by a change of variable or by separating the singular part, and only then apply the trapezoidal rule.

7. FAQ

8. References

• Wikipedia, "Trapezoidal rule"
• L. N. Trefethen & J. A. C. Weideman, "The Exponentially Convergent Trapezoidal Rule," SIAM Review, 56(3), 2014.
• R. L. Burden & J. D. Faires, Numerical Analysis, 10th ed., Cengage, 2016.