Jump to content

Trapezoidal rule: Difference between revisions

From Wikipedia, the free encyclopedia
Content deleted Content added
Citation bot (talk | contribs)
Added bibcode. | Use this bot. Report bugs. | Suggested by Dominic3203 | Category:Numerical integration (quadrature) | #UCB_Category 14/27
 
(14 intermediate revisions by 11 users not shown)
Line 3: Line 3:
[[File:Trapezoidal rule illustration.svg|right|thumb|The function ''f''(''x'') (in blue) is approximated by a linear function (in red).]]
[[File:Trapezoidal rule illustration.svg|right|thumb|The function ''f''(''x'') (in blue) is approximated by a linear function (in red).]]


In [[calculus]], the '''trapezoidal rule''' (also known as the '''trapezoid rule''' or '''trapezium rule'''; see [[Trapezoid]] for more information on terminology) is a technique for approximating the [[integral|definite integral]].
In [[calculus]], the '''trapezoidal rule''' (also known as the '''trapezoid rule''' or '''trapezium rule'''){{efn|See [[Trapezoid]] for more information on terminology.}} is a technique for [[numerical integration]], i.e., approximating the [[integral|definite integral]]:
<math display="block">\int_a^b f(x) \, dx.</math>

:<math>\int_a^b f(x) \, dx.</math>


The trapezoidal rule works by approximating the region under the graph of the function
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. It follows that
<math>f(x)</math> as a [[trapezoid]] and calculating its area. It follows that
<math display="block">\int_{a}^{b} f(x) \, dx \approx (b-a) \cdot \tfrac{1}{2}(f(a)+f(b)).</math>

:<math>\int_{a}^{b} f(x) \, dx \approx (b-a) \cdot \tfrac{1}{2}(f(a)+f(b)).</math>


[[File:WikiTrap.gif|thumb|right|An animation that shows what the trapezoidal rule is and how the error in approximation decreases as the step size decreases]]
[[File:WikiTrap.gif|thumb|right|An animation that shows what the trapezoidal rule is and how the error in approximation decreases as the step size decreases]]
The trapezoidal rule may be viewed as the result obtained by averaging the [[Riemann sum#Left Riemann sum|left]] and [[Riemann sum#Right Riemann sum|right]] [[Riemann sum]]s, and is sometimes defined this way. The integral can be even better approximated by [[Partition of an interval|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
The trapezoidal rule may be viewed as the result obtained by averaging the [[Riemann sum#Left rule|left]] and [[Riemann sum#Right rule|right]] [[Riemann sum]]s, and is sometimes defined this way. The integral can be even better approximated by [[Partition of an interval|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
:<math>\int_a^b f(x) \, dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k.</math>
<math display="block">\int_a^b f(x) \, dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k.</math>
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:.
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:.
:<math>\int_a^b f(x) \, dx \approx \frac{\Delta x}{2} \left(f(x_0) + 2f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \cdots + 2f(x_{N-1}) + f(x_N)\right).</math>
<math display="block">\int_a^b f(x) \, dx \approx \frac{\Delta x}{2} \left(f(x_0) + 2f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \cdots + 2f(x_{N-1}) + f(x_N)\right).</math>


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).
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).
Line 22: Line 20:
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.
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.
[[File:Integration num trapezes notation.svg|thumb|Illustration of "chained trapezoidal rule" used on an irregularly-spaced partition of <math>[a,b]</math>.]]
[[File:Integration num trapezes notation.svg|thumb|Illustration of "chained trapezoidal rule" used on an irregularly-spaced partition of <math>[a,b]</math>.]]



== History ==
== History ==
A 2016 [[Science (journal)|Science]] paper reports that the trapezoid rule was in use in [[Babylon]] before 50 BCE for integrating the velocity of [[Jupiter]] along the [[ecliptic]].<ref>{{cite journal |last=Ossendrijver |first=Mathieu |date=Jan 29, 2016 |title=Ancient Babylonian astronomers calculated Jupiter's position from the area under a time-velocity graph |journal=Science |doi=10.1126/science.aad8085 |pmid=26823423 |volume=351 |issue=6272 |pages=482–484 |s2cid=206644971 |url=https://www.science.org/doi/full/10.1126/science.aad8085}}</ref>
A 2016 ''[[Science (journal)|Science]]'' paper reports that the trapezoid rule was in use in [[Babylon]] before 50 BCE for integrating the velocity of [[Jupiter]] along the [[ecliptic]].<ref>{{cite journal |last=Ossendrijver |first=Mathieu |date=Jan 29, 2016 |title=Ancient Babylonian astronomers calculated Jupiter's position from the area under a time-velocity graph |journal=Science |doi=10.1126/science.aad8085 |pmid=26823423 |volume=351 |issue=6272 |pages=482–484 |bibcode=2016Sci...351..482O |s2cid=206644971 |url=https://www.science.org/doi/full/10.1126/science.aad8085}}</ref>


== Numerical implementation ==
== Numerical implementation ==
Line 31: Line 28:
=== Non-uniform grid ===
=== Non-uniform grid ===
When the grid spacing is non-uniform, one can use the formula
When the grid spacing is non-uniform, one can use the formula
<math display="block"> \int_{a}^{b} f(x)\, dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k ,</math>

:<math> \int_{a}^{b} f(x)\, dx \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k ,</math>
wherein <math>\Delta x_k = x_{k} - x_{k-1} .</math>
wherein <math>\Delta x_k = x_{k}-x_{k-1} .</math>



=== Uniform grid ===
=== Uniform grid ===
For a domain discretized into <math>N</math> equally spaced panels, considerable simplification may occur. Let
For a domain discretized into <math>N</math> equally spaced panels, considerable simplification may occur. Let
:<math>\Delta x_k = \Delta x = \frac{b-a}{N}</math>
<math display="block">\Delta x_k = \Delta x = \frac{b-a}{N}</math>
the approximation to the integral becomes
the approximation to the integral becomes
:<math>\begin{align}
<math display="block">\begin{align}
\int_{a}^{b} f(x)\, dx
\int_{a}^{b} f(x)\, dx
&\approx \frac{\Delta x}{2} \sum_{k=1}^{N} \left( f(x_{k-1}) + f(x_{k}) \right) \\[6pt]
&\approx \frac{\Delta x}{2} \sum_{k=1}^{N} \left( f(x_{k-1}) + f(x_{k}) \right) \\[1ex]
&= \frac{\Delta x}{2} ( f(x_0) + 2f(x_1) + 2f(x_2) + 2f(x_3) + \dotsb + 2f(x_{N-1}) + f(x_N) ) \\[6pt]
&= \frac{\Delta x}{2} \Biggl( f(x_0) + 2f(x_1) + 2f(x_2) + 2f(x_3) + \dotsb + 2f(x_{N-1}) + f(x_N) \Biggr) \\[1ex]
&= \Delta x \left( \sum_{k=1}^{N-1} f(x_k) + \frac{f(x_N) + f(x_0) }{2} \right).
&= \Delta x \left( \frac{f(x_N) + f(x_0) }{2} + \sum_{k=1}^{N-1} f(x_k) \right) .
\end{align}</math>
\end{align}</math>


Line 50: Line 45:
[[File:Trapezium2.gif|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.]]
[[File:Trapezium2.gif|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:
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>

:<math> \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<ref>{{harvtxt|Atkinson|1989|loc=equation (5.1.7)}}</ref>
There exists a number ''ξ'' between ''a'' and ''b'', such that<ref>{{harvtxt|Atkinson|1989|loc=equation (5.1.7)}}</ref>
<math display="block"> \text{E} = -\frac{(b-a)^3}{12N^2} f''(\xi)</math>


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.
:<math> \text{E} = -\frac{(b-a)^3}{12N^2} f''(\xi)</math>

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 error is harder to identify.


An asymptotic error estimate for ''N'' → ∞ is given by
An asymptotic error estimate for ''N'' → ∞ is given by
:<math> \text{E} = -\frac{(b-a)^2}{12N^2} \big[ f'(b)-f'(a) \big] + O(N^{-3}). </math>
<math display="block"> \text{E} = -\frac{(b-a)^2}{12N^2} \big[ f'(b)-f'(a) \big] + O(N^{-3}). </math>
Further terms in this error estimate are given by the Euler–Maclaurin summation formula.
Further terms in this error estimate are given by the Euler–Maclaurin summation formula.


Line 67: Line 60:
#[[Residue calculus]]
#[[Residue calculus]]
#[[Euler–Maclaurin summation formula]]<ref>{{harvtxt|Atkinson|1989|loc=equation (5.1.9)}}</ref><ref>{{harvtxt|Atkinson|1989|loc=p. 285}}</ref>
#[[Euler–Maclaurin summation formula]]<ref>{{harvtxt|Atkinson|1989|loc=equation (5.1.9)}}</ref><ref>{{harvtxt|Atkinson|1989|loc=p. 285}}</ref>
#[[Polynomial interpolation]]<ref>{{harvtxt|Burden|Faires|2011|loc=p. 194}} </ref>
#[[Polynomial interpolation]]<ref>{{harvtxt|Burden|Faires|2011|p=194}} </ref>



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.<ref name="rs90" />
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.<ref name="rs90" />


=== Proof ===
=== Proof ===
First suppose that <math>h=\frac{b-a}{N}</math> and <math>a_k=a+(k-1)h</math>. Let <math> 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
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

<math display="block"> {dg_k \over dt}={1 \over 2}[f(a_k)+f(a_k+t)]+{1\over2}t\cdot f'(a_k+t)-f(a_k+t),</math>
<math display="block"> {dg_k \over dt}={1 \over 2}[f(a_k)+f(a_k+t)]+{1\over2}t\cdot f'(a_k+t)-f(a_k+t),</math>

and
and

<math display="block"> {d^2g_k \over dt^2}={1\over 2}t\cdot f''(a_k+t).</math>
<math display="block"> {d^2g_k \over dt^2}={1\over 2}t\cdot f''(a_k+t).</math>


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

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

which is equivalent to
which is equivalent to
<math> -f''(\xi) \leq f''(a_k+t) \leq f''(\xi)</math>, or <math> -\frac{f''(\xi)t}{2} \leq g_k''(t) \leq \frac{f''(\xi)t}{2}.</math>
<math> -f''(\xi) \leq f''(a_k+t) \leq f''(\xi)</math>, or <math> -\frac{f''(\xi)t}{2} \leq g_k''(t) \leq \frac{f''(\xi)t}{2}.</math>


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

<math display="block"> \int_0^t g_k''(x) dx = g_k'(t)</math> and <math display="block"> \int_0^t g_k'(x) dx = g_k(t).</math>
<math display="block"> \int_0^t g_k''(x) dx = g_k'(t)</math> and <math display="block"> \int_0^t g_k'(x) dx = g_k(t).</math>


Using these results, we find
Using these results, we find

<math display="block"> -\frac{f''(\xi)t^2}{4} \leq g_k'(t) \leq \frac{f''(\xi)t^2}{4}</math>
<math display="block"> -\frac{f''(\xi)t^2}{4} \leq g_k'(t) \leq \frac{f''(\xi)t^2}{4}</math>

and
and

<math display="block"> -\frac{f''(\xi)t^3}{12} \leq g_k(t) \leq \frac{f''(\xi)t^3}{12}</math>
<math display="block"> -\frac{f''(\xi)t^3}{12} \leq g_k(t) \leq \frac{f''(\xi)t^3}{12}</math>


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

<math display="block"> -\frac{f''(\xi)h^3}{12} \leq g_k(h) \leq \frac{f''(\xi)h^3}{12}.</math>
<math display="block"> -\frac{f''(\xi)h^3}{12} \leq g_k(h) \leq \frac{f''(\xi)h^3}{12}.</math>


Summing all of the local error terms 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>
<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
But we also have

<math display="block"> - \sum_{k=1}^N \frac{f''(\xi)h^3}{12} \leq \sum_{k=1}^N g_k(h) \leq \sum_{k=1}^N \frac{f''(\xi)h^3}{12}</math>
<math display="block"> - \sum_{k=1}^N \frac{f''(\xi)h^3}{12} \leq \sum_{k=1}^N g_k(h) \leq \sum_{k=1}^N \frac{f''(\xi)h^3}{12}</math>

and
and

<math display="block"> \sum_{k=1}^N \frac{f''(\xi)h^3}{12}=\frac{f''(\xi)h^3N}{12},</math>
<math display="block"> \sum_{k=1}^N \frac{f''(\xi)h^3}{12}=\frac{f''(\xi)h^3N}{12},</math>


Line 122: Line 100:
Therefore the total error is bounded by
Therefore the total error is bounded by


<math display="block"> \text{error}=\int_a^bf(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>
<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 ===
=== 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
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>
if <math>f</math> is <math>p</math> times continuously differentiable with period <math>T</math>
:<math>\sum_{k=0}^{N-1} f(kh)h =
<math display="block">\sum_{k=0}^{N-1} f(kh)h =
\int_{0}^{T} f(x)\,dx +
\int_0^T f(x)\,dx +
\sum_{k=1}^{\lfloor p/2\rfloor} \frac{B_{2k}}{(2k)!} (f^{(2k - 1)}(T) - f^{(2k - 1)}(0)) - (-1)^{p}h^{p}\int_{0}^{T}\tilde{B}_{p}(x/T)f^{(p)}(x) \, \mathrm{d}x
\sum_{k=1}^{\lfloor p/2\rfloor} \frac{B_{2k}}{(2k)!} (f^{(2k - 1)}(T) - f^{(2k - 1)}(0)) - (-1)^p h^p \int_0^T\tilde{B}_{p}(x/T)f^{(p)}(x) \, dx
</math>
</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.<ref>{{cite book|title=Numerical Analysis, volume 181 of Graduate Texts in Mathematics|first=Rainer|last=Kress |year=1998 |publisher=Springer-Verlag}}</ref> 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 function|Gaussian]], [[Exponentially modified Gaussian distribution|Exponentially modified Gaussian]] and other functions with derivatives at integration limits that can be neglected.<ref>{{Cite journal |last=Goodwin|first=E. T. |date=1949 |title=The evaluation of integrals of the form | journal=Mathematical Proceedings of the Cambridge Philosophical Society |language=en |volume=45 |issue=2 |pages=241–245 |doi=10.1017/S0305004100024786 |bibcode=1949PCPS...45..241G |issn=1469-8064}}</ref> The evaluation of the full integral of a Gaussian function by trapezoidal rule with 1% accuracy can be made using just 4 points.<ref name=":0">{{Cite journal| last1=Kalambet|first1=Yuri |last2=Kozmin|first2=Yuri |last3=Samokhin|first3=Andrey |date=2018 |title=Comparison of integration rules in the case of very narrow chromatographic peaks |journal=Chemometrics and Intelligent Laboratory Systems|volume=179 |pages=22–30 |doi=10.1016/j.chemolab.2018.06.001|issn=0169-7439}}</ref> [[Simpson's rule]] requires 1.8 times more points to achieve the same accuracy.<ref name=":0" /><ref name="w02" />
where <math>h:=T/N</math> and <math>\tilde{B}_{p}</math> is the periodic extension of the <math>p</math>th Bernoulli polynomial.<ref>{{cite book|title=Numerical Analysis, volume 181 of Graduate Texts in Mathematics|first=Rainer|last=Kress|year=1998|publisher=Springer-Verlag}}</ref> Due to the periodicity, the derivatives at the endpoint cancel and we see that the error is <math>O(h^p)</math>.


Although some effort has been made to extend the Euler-Maclaurin summation formula to higher dimensions,<ref>{{cite web |url=https://math.stackexchange.com/q/30384 |website=math.stackexchange.com |title=Euler-Maclaurin Summation Formula for Multiple Sums}}</ref> the most straightforward proof of the rapid convergence of the trapezoidal rule in higher dimensions is to reduce the problem to that of convergence of Fourier series. This line of reasoning shows that if <math>f</math> is periodic on a <math>n</math>-dimensional space with <math>p</math> continuous derivatives, the speed of convergence is <math>O(h^{p/d})</math>. For very large dimension, the shows that Monte-Carlo integration is most likely a better choice, but for 2 and 3 dimensions, equispaced sampling is efficient. This is exploited in computational solid state physics where equispaced sampling over primitive cells in the reciprocal lattice is known as ''Monkhorst-Pack integration''.<ref>{{cite web |url=https://bandgap.io/blog/numerical_integration_bz/| first=Nick|last=Thompson |website=bandgap.io | access-date=19 December 2017 |title=Numerical Integration over Brillouin Zones}}</ref>
A similar effect is available for peak-like functions, such as [[Gaussian function|Gaussian]], [[Exponentially modified Gaussian distribution|Exponentially modified Gaussian]] and other functions with derivatives at integration limits that can be neglected.<ref>{{Cite journal|last=Goodwin|first=E. T.|date=1949|title=The evaluation of integrals of the form|journal=Mathematical Proceedings of the Cambridge Philosophical Society|language=en|volume=45|issue=2|pages=241–245|doi=10.1017/S0305004100024786|issn=1469-8064}}</ref> The evaluation of the full integral of a Gaussian function by trapezoidal rule with 1% accuracy can be made using just 4 points.<ref name=":0">{{Cite journal|last1=Kalambet|first1=Yuri|last2=Kozmin|first2=Yuri|last3=Samokhin|first3=Andrey|date=2018|title=Comparison of integration rules in the case of very narrow chromatographic peaks|journal=Chemometrics and Intelligent Laboratory Systems|volume=179|pages=22–30|doi=10.1016/j.chemolab.2018.06.001|issn=0169-7439}}</ref> [[Simpson's rule]] requires 1.8 times more points to achieve the same accuracy.<ref name=":0" /><ref name="w02" />

Although some effort has been made to extend the Euler-Maclaurin summation formula to higher dimensions,<ref>{{cite web|url=https://math.stackexchange.com/q/30384|website=math.stackexchange.com|title=Euler-Maclaurin Summation Formula for Multiple Sums}}</ref> the most straightforward proof of the rapid convergence of the trapezoidal rule in higher dimensions is to reduce the problem to that of convergence of Fourier series. This line of reasoning shows that if <math>f</math> is periodic on a <math>n</math>-dimensional space with <math>p</math> continuous derivatives, the speed of convergence is <math>O(h^{p/d})</math>. For very large dimension, the shows that Monte-Carlo integration is most likely a better choice, but for 2 and 3 dimensions, equispaced sampling is efficient. This is exploited in computational solid state physics where equispaced sampling over primitive cells in the reciprocal lattice is known as ''Monkhorst-Pack integration''.<ref>{{cite web|url=https://bandgap.io/blog/numerical_integration_bz/|first=Nick|last=Thompson|website=bandgap.io|access-date=19 December 2017|title=Numerical Integration over Brillouin Zones}}</ref>


=== "Rough" functions ===
=== "Rough" functions ===
Line 150: Line 127:
==Example==
==Example==
The following integral is given:
The following integral is given:
<math display="block"> \int_{0.1}^{1.3}{5xe^{- 2x}{dx}} </math>
<math display="block"> \int_{0.1}^{1.3}{5xe^{- 2x}{dx}} </math>
{{ordered list| list-style-type = lower-alpha
| Use the composite trapezoidal rule to estimate the value of this integral. Use three segments.
| Find the true error <math display="inline"> E_{t} </math> for part (a).
a) Use the composite trapezoidal rule to estimate the value of this integral. Use three segments.
b) Find the true error <math display="inline"> E_{t} </math> for part (a).
c) Find the absolute relative true error <math display="inline"> \left| \varepsilon_{t} \right| </math> for part (a).
| Find the absolute relative true error <math display="inline"> \left| \varepsilon_{t} \right| </math> for part (a).
}}
'''Solution'''
'''Solution'''
{{ordered list| list-style-type = lower-alpha
a)  The solution using the composite trapezoidal rule with 3 segments is applied as follows.
| The solution using the composite trapezoidal rule with 3 segments is applied as follows.
<math display="block"> \int_{a}^{b}{f(x){dx}} \approx \frac{b - a}{2n}\left\lbrack f(a) + 2\sum_{i = 1}^{n - 1}{f(a + {ih})} + f(b) \right\rbrack </math>
<math display="block"> \int_{a}^{b}{f(x){dx}} \approx \frac{b - a}{2n}\left\lbrack f(a) + 2\sum_{i = 1}^{n - 1}{f(a + {ih})} + f(b) \right\rbrack </math>
<math display="block"> n = 3 </math>
<math display="block">\begin{align}
n &= 3 \\
a &= 0.1 \\
<math display="block"> a = 0.1 </math>
b &= 1.3 \\
h &= \frac{b - a}{n}
<math display="block"> b = 1.3 </math>
= \frac{1.3 - 0.1}{3}
= 0.4
<math display="block"> \begin{align} h &= \frac{b - a}{n}\\ &= \frac{1.3 - 0.1}{3}\\ &= 0.4 \end{align} </math>
\end{align} </math>
Using the composite trapezoidal rule formula
Using the composite trapezoidal rule formula
<math display="block"> \begin{align} \int_a^b {f(x){dx}} \approx \frac{b - a}{2n} \left\lbrack f(a) + 2\left\{ \sum_{i = 1}^{n - 1}{f(a + {ih})} \right\} + f(b) \right\rbrack\;\;\;\;\;\;\;\;\;\;\;\; (3) \end{align} </math>

<math display="block"> \begin{align} \int_{a}^{b}{f(x){dx}} \approx \frac{b - a}{2n}\left\lbrack f(a) + 2\left\{ \sum_{i = 1}^{n - 1}{f(a + {ih})} \right\} + f(b) \right\rbrack\;\;\;\;\;\;\;\;\;\;\;\; (3) \end{align} </math>
<math display="block"> \begin{align} I &\approx \frac{1.3 - 0.1}{6}\left\lbrack f(0.1) + 2\sum_{i = 1}^{3 - 1}{f(0.1 + 0.4i)} + f(1.3) \right\rbrack\\ I &\approx \frac{1.3 - 0.1}{6}\left\lbrack f(0.1) + 2\sum_{i = 1}^{2}{f(0.1 + 0.4i)} + f(1.3) \right\rbrack\\ &= 0.2\lbrack f(0.1) + 2f(0.5) + 2f(0.9) + f(1.3)\rbrack\\ &= 0.2[5 \times 0.1 \times e^{- 2(0.1)}+2(5 \times 0.5 \times e^{- 2(0.5)})+2(5 \times 0.9 \times e^{- 2(0.9)}) + 5 \times 1.3 \times e^{- 2(1.3)}\rbrack\\ &= 0.84385 \end{align} </math>
b)  The exact value of the above integral can be found by integration by parts and is
<math display="block"> \int_{0.1}^{1.3}{5xe^{- 2x}{dx} = 0.89387} </math>
<math display="block"> \begin{align}
So the true error is
I &\approx \frac{1.3 - 0.1}{6}\left\lbrack f(0.1) + 2\sum_{i = 1}^{3 - 1}{f(0.1 + 0.4i)} + f(1.3) \right\rbrack\\
I &\approx \frac{1.3 - 0.1}{6}\left\lbrack f(0.1) + 2\sum_{i = 1}^{2}{f(0.1 + 0.4i)} + f(1.3) \right\rbrack\\
&= 0.2\lbrack f(0.1) + 2f(0.5) + 2f(0.9) + f(1.3)\rbrack\\
&= 0.2[5 \times 0.1 \times e^{- 2(0.1)}+2(5 \times 0.5 \times e^{- 2(0.5)})+2(5 \times 0.9 \times e^{- 2(0.9)}) + 5 \times 1.3 \times e^{- 2(1.3)}\rbrack\\
&= 0.84385
\end{align} </math>
| The exact value of the above integral can be found by integration by parts and is
<math display="block"> \begin{align} E_{t} &= \text{True Value} - \text{Approximate Value}\\ &= 0.89387 - 0.84385\\ &= 0.05002 \end{align} </math>
<math display="block"> \int_{0.1}^{1.3} 5xe^{- 2x}{dx} = 0.89387 </math>
c)  The absolute relative true error is
So the true error is
<math display="block"> \begin{align}
E_{t} &= \text{True Value} - \text{Approximate Value}\\
&= 0.89387 - 0.84385\\
&= 0.05002
\end{align} </math>
| The absolute relative true error is
<math display="block"> \displaystyle \begin{align}\left| \varepsilon_{t} \right| &= \left| \frac{\text{True Error}}{\text{True Value}} \right| \times 100\\ &= \left| \frac{0.05002}{0.89387} \right| \times 100\\ &= 5.5959\% \end{align} </math>
<math display="block"> \displaystyle \begin{align}\left| \varepsilon_{t} \right| &= \left| \frac{\text{True Error}}{\text{True Value}} \right| \times 100\%\\
&= \left| \frac{0.05002}{0.89387} \right| \times 100\%\\
&= 5.5959\%
\end{align} </math>
}}


==See also==
==See also==
Line 197: Line 183:
* [[Romberg's method]]
* [[Romberg's method]]
* [[Simpson's rule]]
* [[Simpson's rule]]
*[[Volterra integral equation#Numerical Solution using Trapezoidal Rule]]
* {{slink|Volterra integral equation#Numerical Solution using Trapezoidal Rule}}


==Notes==
==Notes==
{{notelist}}
<references/>
<references/>


Line 206: Line 193:
*{{citation |last=Atkinson |first=Kendall E. |year=1989 |title=An Introduction to Numerical Analysis |edition=2nd |publisher=[[John Wiley & Sons]] |location=New York |isbn=978-0-471-50023-0}}
*{{citation |last=Atkinson |first=Kendall E. |year=1989 |title=An Introduction to Numerical Analysis |edition=2nd |publisher=[[John Wiley & Sons]] |location=New York |isbn=978-0-471-50023-0}}
*{{citation |last1=Rahman |first1=Qazi I. |last2=Schmeisser |first2=Gerhard |date=December 1990 |title=Characterization of the speed of convergence of the trapezoidal rule |journal=Numerische Mathematik |issn=0945-3245 |doi=10.1007/BF01386402 |volume=57 |issue=1 |pages=123–138|s2cid=122245944 }}
*{{citation |last1=Rahman |first1=Qazi I. |last2=Schmeisser |first2=Gerhard |date=December 1990 |title=Characterization of the speed of convergence of the trapezoidal rule |journal=Numerische Mathematik |issn=0945-3245 |doi=10.1007/BF01386402 |volume=57 |issue=1 |pages=123–138|s2cid=122245944 }}
*{{citation |last1=Burden |first1=Richard L. |last2=Faires |first2=J. Douglas |year=2000 |title=Numerical Analysis |edition=7th |publisher=Brooks/Cole |isbn=978-0-534-38216-2 |url-access=registration |url=https://archive.org/details/numericalanalysi0000burd }}
*{{citation |last1=Burden |first1=Richard L. |last2=Faires |first2=J. Douglas |year=2011|title=Numerical Analysis |edition=9th |publisher=Brooks/Cole }}
*{{citation |last=Weideman |first=J. A. C. |date=January 2002 |title=Numerical Integration of Periodic Functions: A Few Examples |journal=[[The American Mathematical Monthly]] |doi=10.2307/2695765 |jstor=2695765 |volume=109 |issue=1 |pages=21–36}}
*{{citation |last=Weideman |first=J. A. C. |date=January 2002 |title=Numerical Integration of Periodic Functions: A Few Examples |journal=[[The American Mathematical Monthly]] |doi=10.2307/2695765 |jstor=2695765 |volume=109 |issue=1 |pages=21–36}}
*{{citation |last1=Cruz-Uribe |first1=D. |last2=Neugebauer |first2=C. J. |year=2002 |title=Sharp Error Bounds for the Trapezoidal Rule and Simpson's Rule |journal=Journal of Inequalities in Pure and Applied Mathematics |volume=3 |issue=4 |url=http://www.emis.de/journals/JIPAM/images/031_02_JIPAM/031_02.pdf }}
*{{citation |last1=Cruz-Uribe |first1=D. |last2=Neugebauer |first2=C. J. |year=2002 |title=Sharp Error Bounds for the Trapezoidal Rule and Simpson's Rule |journal=Journal of Inequalities in Pure and Applied Mathematics |volume=3 |issue=4 |url=http://www.emis.de/journals/JIPAM/images/031_02_JIPAM/031_02.pdf }}
Line 217: Line 204:
*[http://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/trapezoidal.html An implementation of trapezoidal quadrature provided by Boost.Math]
*[http://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/trapezoidal.html An implementation of trapezoidal quadrature provided by Boost.Math]


{{Numerical integration}}
{{Calculus topics}}
{{Calculus topics}}



Latest revision as of 05:02, 16 December 2024

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

In calculus, the trapezoidal rule (also known as the trapezoid rule or trapezium rule)[a] is a technique for numerical integration, i.e., approximating the definite integral:

The trapezoidal rule works by approximating the region under the graph of the function as a trapezoid and calculating its area. It follows that

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

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 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 be a partition of such that and be the length of the -th subinterval (that is, ), then When the partition has a regular spacing, as is often the case, that is, when all the have the same value the formula can be simplified for calculation efficiency by factoring out:.

The approximation becomes more accurate as the resolution of the partition increases (that is, for larger , all decrease).

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.

Illustration of "chained trapezoidal rule" used on an irregularly-spaced partition of .

History

[edit]

A 2016 Science paper reports that the trapezoid rule was in use in Babylon before 50 BCE for integrating the velocity of Jupiter along the ecliptic.[1]

Numerical implementation

[edit]

Non-uniform grid

[edit]

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

Uniform grid

[edit]

For a domain discretized into equally spaced panels, considerable simplification may occur. Let the approximation to the integral becomes

Error analysis

[edit]
An animation showing how the trapezoidal rule approximation improves with more strips for an interval with and . As the number of intervals 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:

There exists a number ξ between a and b, such that[2]

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:[3]

  1. Fourier series
  2. Residue calculus
  3. Euler–Maclaurin summation formula[4][5]
  4. Polynomial interpolation[6]

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.[7]

Proof

[edit]

First suppose that and . Let be the function such that is the error of the trapezoidal rule on one of the intervals, . Then and

Now suppose that which holds if is sufficiently smooth. It then follows that which is equivalent to , or

Since and , and

Using these results, we find and

Letting we find

Summing all of the local error terms we find

But we also have and

so that

Therefore the total error is bounded by

Periodic and peak functions

[edit]

The trapezoidal rule converges rapidly for periodic functions. This is an easy consequence of the Euler-Maclaurin summation formula, which says that if is times continuously differentiable with period where and is the periodic extension of the th Bernoulli polynomial.[8] Due to the periodicity, the derivatives at the endpoint cancel and we see that the error is .

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.[9] The evaluation of the full integral of a Gaussian function by trapezoidal rule with 1% accuracy can be made using just 4 points.[10] Simpson's rule requires 1.8 times more points to achieve the same accuracy.[10][11]

Although some effort has been made to extend the Euler-Maclaurin summation formula to higher dimensions,[12] the most straightforward proof of the rapid convergence of the trapezoidal rule in higher dimensions is to reduce the problem to that of convergence of Fourier series. This line of reasoning shows that if is periodic on a -dimensional space with continuous derivatives, the speed of convergence is . For very large dimension, the shows that Monte-Carlo integration is most likely a better choice, but for 2 and 3 dimensions, equispaced sampling is efficient. This is exploited in computational solid state physics where equispaced sampling over primitive cells in the reciprocal lattice is known as Monkhorst-Pack integration.[13]

"Rough" functions

[edit]

For functions that are not in C2, the error bound given above is not applicable. Still, error bounds for such rough functions can be derived, which typically show a slower convergence with the number of function evaluations than the behaviour given above. Interestingly, in this case the trapezoidal rule often has sharper bounds than Simpson's rule for the same number of function evaluations.[14]

Applicability and alternatives

[edit]

The trapezoidal rule is one of a family of formulas for numerical integration called Newton–Cotes formulas, of which the midpoint rule is similar to the trapezoid rule. Simpson's rule is another member of the same family, and in general has faster convergence than the trapezoidal rule for functions which are twice continuously differentiable, though not in all specific cases. However, for various classes of rougher functions (ones with weaker smoothness conditions), the trapezoidal rule has faster convergence in general than Simpson's rule.[14]

Moreover, the trapezoidal rule tends to become extremely accurate when periodic functions are integrated over their periods, which can be analyzed in various ways.[7][11] A similar effect is available for peak functions.[10][11]

For non-periodic functions, however, methods with unequally spaced points such as Gaussian quadrature and Clenshaw–Curtis quadrature are generally far more accurate; Clenshaw–Curtis quadrature can be viewed as a change of variables to express arbitrary integrals in terms of periodic integrals, at which point the trapezoidal rule can be applied accurately.

Example

[edit]

The following integral is given:

  1.  Use the composite trapezoidal rule to estimate the value of this integral. Use three segments.
  2.  Find the true error for part (a).
  3.  Find the absolute relative true error for part (a).

Solution

  1. The solution using the composite trapezoidal rule with 3 segments is applied as follows.

    Using the composite trapezoidal rule formula

  2. The exact value of the above integral can be found by integration by parts and is So the true error is
  3. The absolute relative true error is

See also

[edit]

Notes

[edit]
  1. ^ See Trapezoid for more information on terminology.
  1. ^ Ossendrijver, Mathieu (Jan 29, 2016). "Ancient Babylonian astronomers calculated Jupiter's position from the area under a time-velocity graph". Science. 351 (6272): 482–484. Bibcode:2016Sci...351..482O. doi:10.1126/science.aad8085. PMID 26823423. S2CID 206644971.
  2. ^ Atkinson (1989, equation (5.1.7))
  3. ^ (Weideman 2002, p. 23, section 2)
  4. ^ Atkinson (1989, equation (5.1.9))
  5. ^ Atkinson (1989, p. 285)
  6. ^ Burden & Faires (2011, p. 194)
  7. ^ Jump up to: a b (Rahman & Schmeisser 1990)
  8. ^ Kress, Rainer (1998). Numerical Analysis, volume 181 of Graduate Texts in Mathematics. Springer-Verlag.
  9. ^ Goodwin, E. T. (1949). "The evaluation of integrals of the form". Mathematical Proceedings of the Cambridge Philosophical Society. 45 (2): 241–245. Bibcode:1949PCPS...45..241G. doi:10.1017/S0305004100024786. ISSN 1469-8064.
  10. ^ Jump up to: a b c Kalambet, Yuri; Kozmin, Yuri; Samokhin, Andrey (2018). "Comparison of integration rules in the case of very narrow chromatographic peaks". Chemometrics and Intelligent Laboratory Systems. 179: 22–30. doi:10.1016/j.chemolab.2018.06.001. ISSN 0169-7439.
  11. ^ Jump up to: a b c (Weideman 2002)
  12. ^ "Euler-Maclaurin Summation Formula for Multiple Sums". math.stackexchange.com.
  13. ^ Thompson, Nick. "Numerical Integration over Brillouin Zones". bandgap.io. Retrieved 19 December 2017.
  14. ^ Jump up to: a b (Cruz-Uribe & Neugebauer 2002)

References

[edit]
[edit]