devn00b
/
EQ2EMu


			
				
					
						
						
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246
							<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Trapezoidal Quadrature</title>
<link rel="stylesheet" href="../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../index.html" title="Math Toolkit 2.11.0">
<link rel="up" href="../quadrature.html" title="Chapter&#160;13.&#160;Quadrature and Differentiation">
<link rel="prev" href="../quadrature.html" title="Chapter&#160;13.&#160;Quadrature and Differentiation">
<link rel="next" href="gauss.html" title="Gauss-Legendre quadrature">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
<td align="center"><a href="../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<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>
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="math_toolkit.trapezoidal"></a><a class="link" href="trapezoidal.html" title="Trapezoidal Quadrature">Trapezoidal Quadrature</a>
</h2></div></div></div>
<h4>
<a name="math_toolkit.trapezoidal.h0"></a>
      <span class="phrase"><a name="math_toolkit.trapezoidal.synopsis"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.synopsis">Synopsis</a>
    </h4>
<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>
<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>

<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>
<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">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>
                 <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>
                 <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>
                 <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>

<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>
<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>
                 <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>

<span class="special">}}}</span> <span class="comment">// namespaces</span>
</pre>
<h4>
<a name="math_toolkit.trapezoidal.h1"></a>
      <span class="phrase"><a name="math_toolkit.trapezoidal.description"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.description">Description</a>
    </h4>
<p>
      The functional <code class="computeroutput"><span class="identifier">trapezoidal</span></code>
      calculates the integral of a function <span class="emphasis"><em>f</em></span> using the surprisingly
      simple trapezoidal rule. If we assume only that the integrand is twice continuously
      differentiable, we can prove that the error of the composite trapezoidal rule
      is &#119926;(h<sup>2</sup>). Hence halving the interval only cuts the error by about a fourth,
      which in turn implies that we must evaluate the function many times before
      an acceptable accuracy can be achieved.
    </p>
<p>
      However, the trapezoidal rule has an astonishing property: If the integrand
      is periodic, and we integrate it over a period, then the trapezoidal rule converges
      faster than any power of the step size <span class="emphasis"><em>h</em></span>. This can be
      seen by examination of the <a href="https://en.wikipedia.org/wiki/Euler-Maclaurin_formula" target="_top">Euler-Maclaurin
      summation formula</a>, which relates a definite integral to its trapezoidal
      sum and error terms proportional to the derivatives of the function at the
      endpoints and the Bernoulli numbers. If the derivatives at the endpoints are
      the same or vanish, then the error very nearly vanishes. Hence the trapezoidal
      rule is essentially optimal for periodic integrands.
    </p>
<p>
      Other classes of integrands which are integrated efficiently by this method
      are the C<sub>0</sub><sup>&#8734;</sup>(&#8733;) <a href="https://en.wikipedia.org/wiki/Bump_function" target="_top">bump
      functions</a> and bell-shaped integrals over the infinite interval. For
      details, see <a href="http://epubs.siam.org/doi/pdf/10.1137/130932132" target="_top">Trefethen's</a>
      SIAM review.
    </p>
<p>
      In its simplest form, an integration can be performed by the following code
    </p>
<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>
<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>
<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>
</pre>
<p>
      The integrand must accept a real number argument, but can return a complex
      number. This is useful for contour integrals (which are manifestly periodic)
      and high-order numerical differentiation of analytic functions. An example
      using the integral definition of the complex Bessel function is shown here:
    </p>
<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>
<span class="special">{</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> <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>
    <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span><span class="special">;</span>
    <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">;</span>
    <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>
<span class="special">};</span>

<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>
<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>
</pre>
<p>
      Other special functions which are efficiently evaluated in the complex plane
      by trapezoidal quadrature are modified Bessel functions and the complementary
      error function. Another application of complex-valued trapezoidal quadrature
      is computation of high-order numerical derivatives; see Lyness and Moler for
      details.
    </p>
<p>
      Since the routine is adaptive, step sizes are halved continuously until a tolerance
      is reached. In order to control this tolerance, simply call the routine with
      an additional argument
    </p>
<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>
</pre>
<p>
      The routine stops when successive estimates of the integral <code class="computeroutput"><span class="identifier">I1</span></code>
      and <code class="computeroutput"><span class="identifier">I0</span></code> differ by less than
      the tolerance multiplied by the estimated L<sub>1</sub> norm of the function. A good choice
      for the tolerance is &#8730;&#949;, which is the default. If the integrand is periodic,
      then the number of correct digits should double on each interval halving. Hence,
      once the integration routine has estimated that the error is &#8730;&#949;, then the actual
      error should be ~&#949;. If the integrand is <span class="bold"><strong>not</strong></span>
      periodic, then reducing the error to &#8730;&#949; takes much longer, but is nonetheless
      possible without becoming a major performance bug.
    </p>
<p>
      A question arises as to what to do when successive estimates never pass below
      the tolerance threshold. The stepsize would be halved repeatedly, generating
      an exponential explosion in function evaluations. As such, you may pass an
      optional argument <code class="computeroutput"><span class="identifier">max_refinements</span></code>
      which controls how many times the interval may be halved before giving up.
      By default, this maximum number of refinement steps is 12, which should never
      be hit in double precision unless the function is not periodic. However, for
      higher-precision types, it may be of interest to allow the algorithm to compute
      more refinements:
    </p>
<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>
<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>
</pre>
<p>
      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
      an exception if the maximum refinements is achieved without the requested tolerance
      being attained. This is because the value calculated is more often than not
      still usable. However, for applications with high-reliability requirements,
      the error estimate should be queried. This is achieved by passing additional
      pointers into the routine:
    </p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error_estimate</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</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.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>
<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>
<span class="special">{</span>
     <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>
<span class="special">}</span>
</pre>
<p>
      The final argument is the L<sub>1</sub> norm of the integral. This is computed along with
      the integral, and is an essential component of the algorithm. First, the L<sub>1</sub> norm
      establishes a scale against which the error can be measured. Second, the L<sub>1</sub> norm
      can be used to evaluate the stability of the computation. This can be formulated
      in a rigorous manner by defining the <span class="bold"><strong>condition number
      of summation</strong></span>. The condition number of summation is defined by
    </p>
<div class="blockquote"><blockquote class="blockquote"><p>
        <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>
      </p></blockquote></div>
<p>
      If this number of ~10<sup>k</sup>, then <span class="emphasis"><em>k</em></span> additional digits are expected
      to be lost in addition to digits lost due to floating point rounding error.
      As all numerical quadrature methods reduce to summation, their stability is
      controlled by the ratio &#8747; |f| dt/|&#8747; f dt |, which is easily seen
      to be equivalent to condition number of summation when evaluated numerically.
      Hence both the error estimate and the condition number of summation should
      be analyzed in applications requiring very high precision and reliability.
    </p>
<p>
      As an example, we consider evaluation of Bessel functions by trapezoidal quadrature.
      The Bessel function of the first kind is defined via
    </p>
<div class="blockquote"><blockquote class="blockquote"><p>
        <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
        t - x sin(t)) dt</em></span></span>
      </p></blockquote></div>
<p>
      The integrand is periodic, so the Euler-Maclaurin summation formula guarantees
      exponential convergence via the trapezoidal quadrature. Without careful consideration,
      it seems this would be a very attractive method to compute Bessel functions.
      However, we see that for large <span class="emphasis"><em>n</em></span> the integrand oscillates
      rapidly, taking on positive and negative values, and hence the trapezoidal
      sums become ill-conditioned. In double precision, <span class="emphasis"><em>x = 17</em></span>
      and <span class="emphasis"><em>n = 25</em></span> gives a sum which is so poorly conditioned
      that zero correct digits are obtained.
    </p>
<p>
      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
      be used to control the behaviour of the function: how it handles errors, what
      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
      for more details</a>.
    </p>
<p>
      References:
    </p>
<p>
      Trefethen, Lloyd N., Weideman, J.A.C., <span class="emphasis"><em>The Exponentially Convergent
      Trapezoidal Rule</em></span>, SIAM Review, Vol. 56, No. 3, 2014.
    </p>
<p>
      Stoer, Josef, and Roland Bulirsch. <span class="emphasis"><em>Introduction to numerical analysis.
      Vol. 12.</em></span>, Springer Science &amp; Business Media, 2013.
    </p>
<p>
      Higham, Nicholas J. <span class="emphasis"><em>Accuracy and stability of numerical algorithms.</em></span>
      Society for industrial and applied mathematics, 2002.
    </p>
<p>
      Lyness, James N., and Cleve B. Moler. <span class="emphasis"><em>Numerical differentiation of
      analytic functions.</em></span> SIAM Journal on Numerical Analysis 4.2 (1967):
      202-210.
    </p>
<p>
      Gil, Amparo, Javier Segura, and Nico M. Temme. <span class="emphasis"><em>Computing special
      functions by using quadrature rules.</em></span> Numerical Algorithms 33.1-4
      (2003): 265-275.
    </p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
      Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
      Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
      R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
      Daryle Walker and Xiaogang Zhang<p>
        Distributed under the Boost Software License, Version 1.0. (See accompanying
        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>)
      </p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<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>
</div>
</body>
</html>