libs/math/doc/html/math_toolkit/trapezoidal.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   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.11.0">
   8 <link rel="up" href="../quadrature.html" title="Chapter&#160;13.&#160;Quadrature and Differentiation">
   9 <link rel="prev" href="../quadrature.html" title="Chapter&#160;13.&#160;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">&lt;</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">&gt;</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">&lt;</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">&gt;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;::</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">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;</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 &#119926;(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>&#8734;</sup>(&#8733;) <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">&lt;</span><span class="keyword">double</span><span class="special">&gt;());</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">[&amp;</span><span class="identifier">n</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">z</span><span class="special">](</span><span class="keyword">double</span> <span class="identifier">theta</span><span class="special">)-&gt;</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span>
  94 <span class="special">{</span>
  95     <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;();</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;());</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(),</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 &#8730;&#949;, 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 &#8730;&#949;, then the actual
 125       error should be ~&#949;. If the integrand is <span class="bold"><strong>not</strong></span>
 126       periodic, then reducing the error to &#8730;&#949; 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">&lt;</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">&gt;(),</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(),</span> <span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">max_refinements</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error_estimate</span><span class="special">,</span> <span class="special">&amp;</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">&gt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;());</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>&#954;(S<sub>n</sub>) := &#931;<sub>i</sub><sup>n</sup> |x<sub>i</sub>|/|&#931;<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 &#8747; |f| dt/|&#8747; 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&#928; &#8747;<sub>-&#928;</sub><sup>&#928;</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&#160;20.&#160;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&#160;20.&#160;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 &amp; 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 &#169; 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&#229;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>