1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Trapezoidal Quadrature</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation"> 9<link rel="prev" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation"> 10<link rel="next" href="gauss.html" title="Gauss-Legendre quadrature"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="../quadrature.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="gauss.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.trapezoidal"></a><a class="link" href="trapezoidal.html" title="Trapezoidal Quadrature">Trapezoidal Quadrature</a> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.trapezoidal.h0"></a> 31 <span class="phrase"><a name="math_toolkit.trapezoidal.synopsis"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.synopsis">Synopsis</a> 32 </h4> 33<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">quadrature</span><span class="special">/</span><span class="identifier">trapezoidal</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34<span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">quadrature</span> <span class="special">{</span> 35 36<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span> 37<span class="keyword">auto</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> 38 <span class="identifier">Real</span> <span class="identifier">tol</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()),</span> 39 <span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">12</span><span class="special">,</span> 40 <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error_estimate</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span> 41 <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">);</span> 42 43<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 44<span class="keyword">auto</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">tol</span><span class="special">,</span> <span class="identifier">size_t</span> <span class="identifier">max_refinements</span><span class="special">,</span> 45 <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error_estimate</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&</span> <span class="identifier">pol</span><span class="special">);</span> 46 47<span class="special">}}}</span> <span class="comment">// namespaces</span> 48</pre> 49<h4> 50<a name="math_toolkit.trapezoidal.h1"></a> 51 <span class="phrase"><a name="math_toolkit.trapezoidal.description"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.description">Description</a> 52 </h4> 53<p> 54 The functional <code class="computeroutput"><span class="identifier">trapezoidal</span></code> 55 calculates the integral of a function <span class="emphasis"><em>f</em></span> using the surprisingly 56 simple trapezoidal rule. If we assume only that the integrand is twice continuously 57 differentiable, we can prove that the error of the composite trapezoidal rule 58 is ��(h<sup>2</sup>). Hence halving the interval only cuts the error by about a fourth, 59 which in turn implies that we must evaluate the function many times before 60 an acceptable accuracy can be achieved. 61 </p> 62<p> 63 However, the trapezoidal rule has an astonishing property: If the integrand 64 is periodic, and we integrate it over a period, then the trapezoidal rule converges 65 faster than any power of the step size <span class="emphasis"><em>h</em></span>. This can be 66 seen by examination of the <a href="https://en.wikipedia.org/wiki/Euler-Maclaurin_formula" target="_top">Euler-Maclaurin 67 summation formula</a>, which relates a definite integral to its trapezoidal 68 sum and error terms proportional to the derivatives of the function at the 69 endpoints and the Bernoulli numbers. If the derivatives at the endpoints are 70 the same or vanish, then the error very nearly vanishes. Hence the trapezoidal 71 rule is essentially optimal for periodic integrands. 72 </p> 73<p> 74 Other classes of integrands which are integrated efficiently by this method 75 are the C<sub>0</sub><sup>∞</sup>(∝) <a href="https://en.wikipedia.org/wiki/Bump_function" target="_top">bump 76 functions</a> and bell-shaped integrals over the infinite interval. For 77 details, see <a href="http://epubs.siam.org/doi/pdf/10.1137/130932132" target="_top">Trefethen's</a> 78 SIAM review. 79 </p> 80<p> 81 In its simplest form, an integration can be performed by the following code 82 </p> 83<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">quadrature</span><span class="special">::</span><span class="identifier">trapezoidal</span><span class="special">;</span> 84<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">5</span> <span class="special">-</span> <span class="number">4</span><span class="special">*</span><span class="identifier">cos</span><span class="special">(</span><span class="identifier">x</span><span class="special">));</span> <span class="special">};</span> 85<span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>());</span> 86</pre> 87<p> 88 The integrand must accept a real number argument, but can return a complex 89 number. This is useful for contour integrals (which are manifestly periodic) 90 and high-order numerical differentiation of analytic functions. An example 91 using the integral definition of the complex Bessel function is shown here: 92 </p> 93<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">bessel_integrand</span> <span class="special">=</span> <span class="special">[&</span><span class="identifier">n</span><span class="special">,</span> <span class="special">&</span><span class="identifier">z</span><span class="special">](</span><span class="keyword">double</span> <span class="identifier">theta</span><span class="special">)-></span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> 94<span class="special">{</span> 95 <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">z</span><span class="special">{</span><span class="number">2</span><span class="special">,</span> <span class="number">3</span><span class="special">};</span> 96 <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span><span class="special">;</span> 97 <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">;</span> 98 <span class="keyword">return</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">z</span><span class="special">*</span><span class="identifier">sin</span><span class="special">(</span><span class="identifier">theta</span><span class="special">)</span> <span class="special">-</span> <span class="number">2</span><span class="special">*</span><span class="identifier">theta</span><span class="special">)/</span><span class="identifier">pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span> 99<span class="special">};</span> 100 101<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">quadrature</span><span class="special">::</span><span class="identifier">trapezoidal</span><span class="special">;</span> 102<span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">Jnz</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">bessel_integrand</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">pi</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>());</span> 103</pre> 104<p> 105 Other special functions which are efficiently evaluated in the complex plane 106 by trapezoidal quadrature are modified Bessel functions and the complementary 107 error function. Another application of complex-valued trapezoidal quadrature 108 is computation of high-order numerical derivatives; see Lyness and Moler for 109 details. 110 </p> 111<p> 112 Since the routine is adaptive, step sizes are halved continuously until a tolerance 113 is reached. In order to control this tolerance, simply call the routine with 114 an additional argument 115 </p> 116<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(),</span> <span class="number">1e-6</span><span class="special">);</span> 117</pre> 118<p> 119 The routine stops when successive estimates of the integral <code class="computeroutput"><span class="identifier">I1</span></code> 120 and <code class="computeroutput"><span class="identifier">I0</span></code> differ by less than 121 the tolerance multiplied by the estimated L<sub>1</sub> norm of the function. A good choice 122 for the tolerance is √ε, which is the default. If the integrand is periodic, 123 then the number of correct digits should double on each interval halving. Hence, 124 once the integration routine has estimated that the error is √ε, then the actual 125 error should be ~ε. If the integrand is <span class="bold"><strong>not</strong></span> 126 periodic, then reducing the error to √ε takes much longer, but is nonetheless 127 possible without becoming a major performance bug. 128 </p> 129<p> 130 A question arises as to what to do when successive estimates never pass below 131 the tolerance threshold. The stepsize would be halved repeatedly, generating 132 an exponential explosion in function evaluations. As such, you may pass an 133 optional argument <code class="computeroutput"><span class="identifier">max_refinements</span></code> 134 which controls how many times the interval may be halved before giving up. 135 By default, this maximum number of refinement steps is 12, which should never 136 be hit in double precision unless the function is not periodic. However, for 137 higher-precision types, it may be of interest to allow the algorithm to compute 138 more refinements: 139 </p> 140<pre class="programlisting"><span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">15</span><span class="special">;</span> 141<span class="keyword">long</span> <span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0L</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">>(),</span> <span class="number">1e-9L</span><span class="special">,</span> <span class="identifier">max_refinements</span><span class="special">);</span> 142</pre> 143<p> 144 Note that the maximum allowed compute time grows exponentially with <code class="computeroutput"><span class="identifier">max_refinements</span></code>. The routine will not throw 145 an exception if the maximum refinements is achieved without the requested tolerance 146 being attained. This is because the value calculated is more often than not 147 still usable. However, for applications with high-reliability requirements, 148 the error estimate should be queried. This is achieved by passing additional 149 pointers into the routine: 150 </p> 151<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error_estimate</span><span class="special">;</span> 152<span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</span> 153<span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(),</span> <span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">max_refinements</span><span class="special">,</span> <span class="special">&</span><span class="identifier">error_estimate</span><span class="special">,</span> <span class="special">&</span><span class="identifier">L1</span><span class="special">);</span> 154<span class="keyword">if</span> <span class="special">(</span><span class="identifier">error_estimate</span> <span class="special">></span> <span class="identifier">tolerance</span><span class="special">*</span><span class="identifier">L1</span><span class="special">)</span> 155<span class="special">{</span> 156 <span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">some_other_quadrature_method</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>());</span> 157<span class="special">}</span> 158</pre> 159<p> 160 The final argument is the L<sub>1</sub> norm of the integral. This is computed along with 161 the integral, and is an essential component of the algorithm. First, the L<sub>1</sub> norm 162 establishes a scale against which the error can be measured. Second, the L<sub>1</sub> norm 163 can be used to evaluate the stability of the computation. This can be formulated 164 in a rigorous manner by defining the <span class="bold"><strong>condition number 165 of summation</strong></span>. The condition number of summation is defined by 166 </p> 167<div class="blockquote"><blockquote class="blockquote"><p> 168 <span class="serif_italic"><span class="emphasis"><em>κ(S<sub>n</sub>) := Σ<sub>i</sub><sup>n</sup> |x<sub>i</sub>|/|Σ<sub>i</sub><sup>n</sup> x<sub>i</sub>|</em></span></span> 169 </p></blockquote></div> 170<p> 171 If this number of ~10<sup>k</sup>, then <span class="emphasis"><em>k</em></span> additional digits are expected 172 to be lost in addition to digits lost due to floating point rounding error. 173 As all numerical quadrature methods reduce to summation, their stability is 174 controlled by the ratio ∫ |f| dt/|∫ f dt |, which is easily seen 175 to be equivalent to condition number of summation when evaluated numerically. 176 Hence both the error estimate and the condition number of summation should 177 be analyzed in applications requiring very high precision and reliability. 178 </p> 179<p> 180 As an example, we consider evaluation of Bessel functions by trapezoidal quadrature. 181 The Bessel function of the first kind is defined via 182 </p> 183<div class="blockquote"><blockquote class="blockquote"><p> 184 <span class="serif_italic"><span class="emphasis"><em>J<sub>n</sub>(x) = 1/2Π ∫<sub>-Π</sub><sup>Π</sup> cos(n 185 t - x sin(t)) dt</em></span></span> 186 </p></blockquote></div> 187<p> 188 The integrand is periodic, so the Euler-Maclaurin summation formula guarantees 189 exponential convergence via the trapezoidal quadrature. Without careful consideration, 190 it seems this would be a very attractive method to compute Bessel functions. 191 However, we see that for large <span class="emphasis"><em>n</em></span> the integrand oscillates 192 rapidly, taking on positive and negative values, and hence the trapezoidal 193 sums become ill-conditioned. In double precision, <span class="emphasis"><em>x = 17</em></span> 194 and <span class="emphasis"><em>n = 25</em></span> gives a sum which is so poorly conditioned 195 that zero correct digits are obtained. 196 </p> 197<p> 198 The final <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can 199 be used to control the behaviour of the function: how it handles errors, what 200 level of precision to use etc. Refer to the <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">policy documentation 201 for more details</a>. 202 </p> 203<p> 204 References: 205 </p> 206<p> 207 Trefethen, Lloyd N., Weideman, J.A.C., <span class="emphasis"><em>The Exponentially Convergent 208 Trapezoidal Rule</em></span>, SIAM Review, Vol. 56, No. 3, 2014. 209 </p> 210<p> 211 Stoer, Josef, and Roland Bulirsch. <span class="emphasis"><em>Introduction to numerical analysis. 212 Vol. 12.</em></span>, Springer Science & Business Media, 2013. 213 </p> 214<p> 215 Higham, Nicholas J. <span class="emphasis"><em>Accuracy and stability of numerical algorithms.</em></span> 216 Society for industrial and applied mathematics, 2002. 217 </p> 218<p> 219 Lyness, James N., and Cleve B. Moler. <span class="emphasis"><em>Numerical differentiation of 220 analytic functions.</em></span> SIAM Journal on Numerical Analysis 4.2 (1967): 221 202-210. 222 </p> 223<p> 224 Gil, Amparo, Javier Segura, and Nico M. Temme. <span class="emphasis"><em>Computing special 225 functions by using quadrature rules.</em></span> Numerical Algorithms 33.1-4 226 (2003): 265-275. 227 </p> 228</div> 229<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 230<td align="left"></td> 231<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 232 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 233 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 234 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 235 Daryle Walker and Xiaogang Zhang<p> 236 Distributed under the Boost Software License, Version 1.0. (See accompanying 237 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 238 </p> 239</div></td> 240</tr></table> 241<hr> 242<div class="spirit-nav"> 243<a accesskey="p" href="../quadrature.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="gauss.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 244</div> 245</body> 246</html> 247
