Trapezoidal rule

right|thumb|The function f(x) (in blue) is approximated by a linear function (in red).

In calculus, the trapezoidal rule (informally trapezoid rule; or in British English trapezium rule) is a technique for numerical integration, i.e. approximating the definite integral:

<!--

Note: the "animation" below is not animated! **

-->

thumb|right|An animation that shows what the trapezoidal rule is and how the error in approximation decreases as the step size decreases

<!--

Note: the "animation" above is not animated! **

-->

The trapezoidal rule works by approximating the region under the graph of the function

<math>f(x)</math> as a trapezoid and calculating its area. This is easily calculated by noting that the area of the region is made up of a rectangle with width <math>(b - a)</math> and height <math>f(a)</math>, and a triangle of width <math>(b-a)</math> and height <math>f(b) - f(a)</math>.

Therefore,

\begin{align}

\int_a^b f(x) \, dx

&\approx \underbrace{(b - a) \cdot f(a)}_{\text{Area of rectangle + \underbrace{\tfrac{1}{2} (b - a) \cdot [f(b) - f(a)]}_{\text{Area of triangle \\

&= (b - a) \cdot \left(f(a) + \tfrac{1}{2} f(b) - \tfrac{1}{2} f(a)\right) \\

&= (b - a) \cdot \left(\tfrac{1}{2} f(a) + \tfrac{1}{2} f(b)\right) \\

&= \frac{1}{2} (b - a) [f(a) + f(b)].

\end{align}

</math>

The rule can also be derived by replacing the integrand with the equation of the line joining points <math>\big(a, f(a)\big)</math> and <math>\big(b, f(b)\big)</math>, which using the two point form of the equation of a line, is

y = (x - a)\,\frac{f(b) - f(a)}{b - a} + f(a).

</math>

Therefore,

\begin{align}

\int_a^b f(x) \, dx

&\approx \int_a^b (x - a)\,\frac{f(b) - f(a)}{b - a} + f(a) \, dx \\

&= \left[\frac{1}{2} (x - a)^2 \,\frac{f(b) - f(a)}{b - a} + xf(a) \right]_{x=a}^{x=b} \\

&= \left[\frac{1}{2} (b - a)^2 \,\frac{f(b) - f(a)}{b - a} + bf(a) \right] - af(a) \\

&= \frac{1}{2} (b - a) [f(b) - f(a)] + (b - a) f(a) \\

&= \frac{1}{2} (b - a) [f(a) + f(b)],

\end{align}

</math>

as before.

[[File:Integration num trapezes notation.svg|thumb|Illustration of "chained trapezoidal rule" used on an irregularly spaced partition of <math>[a, b]</math>]]

The integral can be even better approximated by partitioning the integration interval, applying the trapezoidal rule to each subinterval and summing the results. In practice, this "chained" (or "composite") trapezoidal rule is usually what is meant by "integrating with the trapezoidal rule". Let <math>\{x_k\}</math> be a partition of <math>[a, b]</math> such that <math>a = x_0 < x_1 < \cdots < x_{N-1} < x_N = b,</math> and <math>\Delta x_k</math> be the length of the <math>k</math>-th subinterval (that is, <math>\Delta x_k = x_k - x_{k-1}</math>), then

\int_a^b f(x) \, dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k.

</math>

The trapezoidal rule may be viewed as the result obtained by averaging the left and right Riemann sums and is sometimes defined this way.

The approximation becomes more accurate as the resolution of the partition increases (that is, for larger <math>N</math>, all <math>\Delta x_k</math> decrease).

When the partition has a regular spacing, as is often the case, that is, when all the <math>\Delta x_k</math> have the same value <math>\Delta x,</math> the formula can be simplified for calculation efficiency by factoring <math>\Delta x</math> out:

\int_a^b f(x) \, dx \approx \Delta x \left(\frac{f(x_0) + f(x_N)} 2 + \sum_{k=1}^{N-1} f(x_k) \right).

</math>

As discussed below, it is also possible to place error bounds on the accuracy of the value of a definite integral estimated using a trapezoidal rule.

History

Evidence discovered in Babylonian cuniform tablets circa 200 BCE predicting the motion of Jupiter along the ecliptic suggests that the trapezoid rule was in use for numerical approximation long before calculus itself was invented.

Numerical implementation

Non-uniform grid

When the grid spacing is non-uniform, one can use the formula

\int_a^b f(x) \,dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k,

</math>

where <math>\Delta x_k = x_k - x_{k-1},</math> or more a computationally efficient formula

\int_a^b f(x) \,dx \approx \frac{1}{2} \biggl(f(x_0) \Delta_{+1} x_0 + f(x_N) \Delta_{-1} x_N + \sum_{k=1}^{N-1} f(x_k) \Delta_{\pm1} x_k\biggr),

</math>

where <math>\Delta_{+1} x_0 = x_1 - x_0,</math> <math>\Delta_{-1} x_N = x_N - x_{N-1},</math> <math>\Delta_{\pm1} x_k = x_{k+1} - x_{k-1}</math> are the corresponding forward, backward, and central differences.

Uniform grid

For a domain partitioned by <math>N</math> equally spaced points, considerable simplification may occur.

Let <math>\Delta x = \frac{b - a}{N}</math> and <math>x_k = a + k \Delta x</math> for <math>k = 0, 1, \ldots, N.</math>

The approximation to the integral becomes

<math display="block">\begin{align}

\int_a^b f(x) \,dx

&\approx \frac{\Delta x}{2} \sum_{k=1}^N [f(x_{k-1}) + f(x_{k})] \\

&= \Delta x \biggl( \tfrac12f(x_0) + \tfrac12f(x_N) + \sum_{k=1}^{N-1} f(x_k) \biggr).

\end{align}</math>

Sometimes this expression is written as

\Delta x \!\! \mathop{\ \sum\vphantom\big)'}_{k=0}^{N} f(x_k),

</math>

where the symbol indicates that the first and last terms are halved.

Error analysis

right|thumb|An animation showing how the trapezoidal rule approximation improves with more strips for an interval with <math>a=2</math> and <math>b=8</math>. As the number of intervals <math>N</math> increases, so too does the accuracy of the result.

The error of the composite trapezoidal rule is the difference between the value of the integral and the numerical result:

<math display="block"> \text{E} = \int_a^b f(x)\,dx - \frac{b-a}{N} \left[ {f(a) + f(b) \over 2} + \sum_{k=1}^{N-1} f \left( a+k \frac{b-a}{N} \right) \right]</math>

There exists a number ξ between a and b, such that

It follows that if the integrand is concave up (and thus has a positive second derivative), then the error is negative and the trapezoidal rule overestimates the true value. This can also be seen from the geometric picture: the trapezoids include all of the area under the curve and extend over it. Similarly, a concave-down function yields an underestimate because area is unaccounted for under the curve, but none is counted above. If the interval of the integral being approximated includes an inflection point, the sign of the error is harder to identify.

An asymptotic error estimate for N → ∞ is given by

Further terms in this error estimate are given by the Euler–Maclaurin summation formula.

Several techniques can be used to analyze the error, including:

Fourier series
Residue calculus
Euler–Maclaurin summation formula
Polynomial interpolation

It is argued that the speed of convergence of the trapezoidal rule reflects and can be used as a definition of classes of smoothness of the functions.

Proof

First suppose that <math>h=\frac{b-a}{N}</math> and <math>a_k=a+(k-1)h</math>. Let <math display="block"> g_k(t) = \frac{1}{2} t[f(a_k)+f(a_k+t)] - \int_{a_k}^{a_k+t} f(x) \, dx</math> be the function such that <math> |g_k(h)| </math> is the error of the trapezoidal rule on one of the intervals, <math> [a_k, a_k+h] </math>. Then

and

Now suppose that <math> \left| f(x) \right| \leq \left| f(\xi) \right|, </math> which holds if <math> f </math> is sufficiently smooth. It then follows that

<math display="block"> \left| f(a_k+t) \right| \leq f(\xi)</math>

which is equivalent to

Since <math> g_k'(0)=0</math> and <math> g_k(0)=0</math>,

Using these results, we find

and

Letting <math> t = h </math> we find

Summing all of the local error terms we find

<math display="block"> \sum_{k=1}^{N} g_k(h) = \frac{b-a}{N} \left[ {f(a) + f(b) \over 2} + \sum_{k=1}^{N-1} f \left( a+k \frac{b-a}{N} \right) \right] - \int_a^b f(x)dx.</math>

But we also have

and

so that

<math display="block"> -\frac{f(\xi)h^3N}{12} \leq \frac{b-a}{N} \left[ {f(a) + f(b) \over 2} + \sum_{k=1}^{N-1} f \left( a+k \frac{b-a}{N} \right) \right]-\int_a^bf(x)dx \leq \frac{f(\xi)h^3N}{12}.</math>

Therefore the total error is bounded by

<math display="block"> \text{error} = \int_a^b f(x)\,dx - \frac{b-a}{N} \left[ {f(a) + f(b) \over 2} + \sum_{k=1}^{N-1} f \left( a+k \frac{b-a}{N} \right) \right] = \frac{f(\xi)h^3N}{12}=\frac{f(\xi)(b-a)^3}{12N^2}.</math>

Periodic and peak functions

The trapezoidal rule converges rapidly for periodic functions. This is an easy consequence of the Euler–Maclaurin summation formula, which says that

if <math>f</math> is <math>p</math> times continuously differentiable with period <math>T,</math> then

\sum_{k=0}^{N-1} f(kh)h =

\int_0^T f(x)\,dx +

\sum_{k=1}^{\lfloor p/2\rfloor} \frac{B_{2k{(2k)!} \left(f^{(2k - 1)}(T) - f^{(2k - 1)}(0)\right) - (-1)^p h^p \int_0^T\tilde{B}_{p}(x/T)f^{(p)}(x) \, dx,

</math>

where <math>h := T/N,</math> and <math>\tilde{B}_p</math> is the periodic extension of the <math>p</math>-th Bernoulli polynomial. Due to the periodicity, the derivatives at the endpoint cancel, and we see that the error is <math>O(h^p)</math>.

A similar effect is available for peak-like functions, such as Gaussian, Exponentially modified Gaussian and other functions with derivatives at integration limits that can be neglected. The evaluation of the full integral of a Gaussian function by trapezoidal rule with 1% accuracy can be made using just 4 points. Simpson's rule requires 1.8 times more points to achieve the same accuracy.