devn00b
/
EQ2EMu


			
				
					
						
						
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288
							<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Gauss-Kronrod 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="gauss.html" title="Gauss-Legendre quadrature">
<link rel="next" href="double_exponential.html" title="Double-exponential 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="gauss.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="double_exponential.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.gauss_kronrod"></a><a class="link" href="gauss_kronrod.html" title="Gauss-Kronrod Quadrature">Gauss-Kronrod Quadrature</a>
</h2></div></div></div>
<h4>
<a name="math_toolkit.gauss_kronrod.h0"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.overview"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.overview">Overview</a>
    </h4>
<p>
      Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides
      an a posteriori error estimate for the integral.
    </p>
<p>
      The idea behind Gaussian quadrature is to choose <span class="emphasis"><em>n</em></span> nodes
      and weights in such a way that polynomials of order <span class="emphasis"><em>2n-1</em></span>
      are integrated exactly. However, integration of polynomials is trivial, so
      it is rarely done via numerical methods. Instead, transcendental and numerically
      defined functions are integrated via Gaussian quadrature, and the defining
      problem becomes how to estimate the remainder. Gaussian quadrature alone (without
      some form of interval splitting) cannot answer this question.
    </p>
<p>
      It is possible to compute a Gaussian quadrature of order <span class="emphasis"><em>n</em></span>
      and another of order (say) <span class="emphasis"><em>2n+1</em></span>, and use the difference
      as an error estimate. However, this is not optimal, as the zeros of the Legendre
      polynomials (nodes of the Gaussian quadrature) are never the same for different
      orders, so <span class="emphasis"><em>3n+1</em></span> function evaluations must be performed.
      Kronrod considered the problem of how to interleave nodes into a Gaussian quadrature
      in such a way that all previous function evaluations can be reused, while increasing
      the order of polynomials that can be integrated exactly. This allows an a posteriori
      error estimate to be provided while still preserving exponential convergence.
      Kronrod discovered that by adding <span class="emphasis"><em>n+1</em></span> nodes (computed
      from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature
      of order <span class="emphasis"><em>n</em></span>, he could integrate polynomials of order <span class="emphasis"><em>3n+1</em></span>.
    </p>
<p>
      The integration routines provided here will perform either adaptive or non-adaptive
      quadrature, they should be chosen for the integration of smooth functions with
      no end point singularities. For difficult functions, or those with end point
      singularities, please refer to the <a class="link" href="double_exponential.html" title="Double-exponential quadrature">double-exponential
      integration schemes</a>.
    </p>
<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">gauss_kronrod</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span>

<span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">unsigned</span> <span class="identifier">N</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">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">&lt;&gt;</span> <span class="special">&gt;</span>
<span class="keyword">class</span> <span class="identifier">gauss_kronrod</span>
<span class="special">{</span>
   <span class="keyword">static</span> <span class="keyword">const</span> <span class="identifier">RandomAccessContainer</span><span class="special">&amp;</span> <span class="identifier">abscissa</span><span class="special">();</span>
   <span class="keyword">static</span> <span class="keyword">const</span> <span class="identifier">RandomAccessContainer</span><span class="special">&amp;</span> <span class="identifier">weights</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">&gt;</span>
   <span class="keyword">static</span> <span class="keyword">auto</span> <span class="identifier">integrate</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="keyword">unsigned</span> <span class="identifier">max_depth</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span>
                         <span class="identifier">Real</span> <span class="identifier">tol</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
                         <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</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">pL1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()));</span>
<span class="special">};</span>
</pre>
<h4>
<a name="math_toolkit.gauss_kronrod.h1"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.description"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.description">Description</a>
    </h4>
<pre class="programlisting"><span class="keyword">static</span> <span class="keyword">const</span> <span class="identifier">RandomAccessContainer</span><span class="special">&amp;</span> <span class="identifier">abscissa</span><span class="special">();</span>
<span class="keyword">static</span> <span class="keyword">const</span> <span class="identifier">RandomAccessContainer</span><span class="special">&amp;</span> <span class="identifier">weights</span><span class="special">();</span>
</pre>
<p>
      These functions provide direct access to the abscissa and weights used to perform
      the quadrature: the return type depends on the <span class="emphasis"><em>Points</em></span>
      template parameter, but is always a RandomAccessContainer type. Note that only
      positive (or zero) abscissa and weights are stored, and that they contain both
      the Gauss and Kronrod points.
    </p>
<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
<span class="keyword">static</span> <span class="keyword">auto</span> <span class="identifier">integrate</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="keyword">unsigned</span> <span class="identifier">max_depth</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span>
                            <span class="identifier">Real</span> <span class="identifier">tol</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
                            <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</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">pL1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()));</span>
</pre>
<p>
      Performs adaptive Gauss-Kronrod quadrature on function <span class="emphasis"><em>f</em></span>
      over the range (a,b).
    </p>
<p>
      <span class="emphasis"><em>max_depth</em></span> sets the maximum number of interval splittings
      permitted before stopping. Set this to zero for non-adaptive quadrature. Note
      that the algorithm descends the tree depth first, so only "difficult"
      areas of the integral result in interval splitting.
    </p>
<p>
      <span class="emphasis"><em>tol</em></span> Sets the maximum relative error in the result, this
      should not be set too close to machine epsilon or the function will simply
      thrash and probably not return accurate results. On the other hand the default
      may be overly-pressimistic.
    </p>
<p>
      <span class="emphasis"><em>error</em></span> When non-null, <code class="computeroutput"><span class="special">*</span><span class="identifier">error</span></code> is set to the difference between the
      (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation.
    </p>
<p>
      <span class="emphasis"><em>pL1</em></span> When non-null, <code class="computeroutput"><span class="special">*</span><span class="identifier">pL1</span></code> is set to the L1 norm of the result,
      if there is a significant difference between this and the returned value, then
      the result is likely to be ill-conditioned.
    </p>
<h4>
<a name="math_toolkit.gauss_kronrod.h2"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.choosing_the_number_of_points"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.choosing_the_number_of_points">Choosing
      the number of points</a>
    </h4>
<p>
      The number of points specified in the <span class="emphasis"><em>Points</em></span> template
      parameter must be an odd number: giving a (N-1)/2 Gauss quadrature as the comparison
      for error estimation.
    </p>
<p>
      Internally class <code class="computeroutput"><span class="identifier">gauss_kronrod</span></code>
      has pre-computed tables of abscissa and weights for 15, 31, 41, 51 and 61 Gauss-Kronrod
      points at up to 100-decimal digit precision. That means that using for example,
      <code class="computeroutput"><span class="identifier">gauss_kronrod</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="number">31</span><span class="special">&gt;::</span><span class="identifier">integrate</span></code>
      incurs absolutely zero set-up overhead from computing the abscissa/weight pairs.
      When using multiprecision types with less than 100 digits of precision, then
      there is a small initial one time cost, while the abscissa/weight pairs are
      constructed from strings.
    </p>
<p>
      However, for types with higher precision, or numbers of points other than those
      given above, the abscissa/weight pairs are computed when first needed and then
      cached for future use, which does incur a noticeable overhead. If this is likely
      to be an issue, then:
    </p>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
<li class="listitem">
          Defining BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND will result in a compile-time
          error, whenever a combination of number type and number of points is used
          which does not have pre-computed values.
        </li>
<li class="listitem">
          There is a program <a href="../../../tools/gauss_kronrod_constants.cpp" target="_top">gauss_kronrod_constants.cpp</a>
          which was used to provide the pre-computed values already in gauss_kronrod.hpp.
          The program can be trivially modified to generate code and constants for
          other precisions and numbers of points.
        </li>
</ul></div>
<h4>
<a name="math_toolkit.gauss_kronrod.h3"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.complex_quadrature"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.complex_quadrature">Complex
      Quadrature</a>
    </h4>
<p>
      The Gauss-Kronrod quadrature support integrands defined on the real line and
      returning complex values. In this case, the template argument is the real type,
      and the complex type is deduced via the return type of the function.
    </p>
<h4>
<a name="math_toolkit.gauss_kronrod.h4"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.examples"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.examples">Examples</a>
    </h4>
<p>
      We'll begin by integrating exp(-t<sup>2</sup>/2) over (0,+&#8734;) using a 7 term Gauss rule
      and 15 term Kronrod rule, and begin by defining the function to integrate as
      a C++ lambda expression:
    </p>
<pre class="programlisting"><span class="keyword">using</span> <span class="keyword">namespace</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="keyword">auto</span> <span class="identifier">f1</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">t</span><span class="special">*</span><span class="identifier">t</span> <span class="special">/</span> <span class="number">2</span><span class="special">);</span> <span class="special">};</span>
</pre>
<p>
      We'll start off with a one shot (ie non-adaptive) integration, and keep track
      of the estimated error:
    </p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">gauss_kronrod</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">&gt;::</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f1</span><span class="special">,</span> <span class="number">0</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="keyword">double</span><span class="special">&gt;::</span><span class="identifier">infinity</span><span class="special">(),</span> <span class="number">0</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">);</span>
</pre>
<p>
      This yields Q = 1.25348207361, which has an absolute error of 1e-4 compared
      to the estimated error of 5e-3: this is fairly typical, with the difference
      between Gauss and Gauss-Kronrod schemes being much higher than the actual error.
      Before moving on to adaptive quadrature, lets try again with more points, in
      fact with the largest Gauss-Kronrod scheme we have cached (30/61):
    </p>
<pre class="programlisting"><span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">gauss_kronrod</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="number">61</span><span class="special">&gt;::</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f1</span><span class="special">,</span> <span class="number">0</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="keyword">double</span><span class="special">&gt;::</span><span class="identifier">infinity</span><span class="special">(),</span> <span class="number">0</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">);</span>
</pre>
<p>
      This yields an absolute error of 3e-15 against an estimate of 1e-8, which is
      about as good as we're going to get at double precision
    </p>
<p>
      However, instead of continuing with ever more points, lets switch to adaptive
      integration, and set the desired relative error to 1e-14 against a maximum
      depth of 5:
    </p>
<pre class="programlisting"><span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">gauss_kronrod</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">&gt;::</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f1</span><span class="special">,</span> <span class="number">0</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="keyword">double</span><span class="special">&gt;::</span><span class="identifier">infinity</span><span class="special">(),</span> <span class="number">5</span><span class="special">,</span> <span class="number">1e-14</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">);</span>
</pre>
<p>
      This yields an actual error of zero, against an estimate of 4e-15. In fact
      in this case the requested tolerance was almost certainly set too low: as we've
      seen above, for smooth functions, the precision achieved is often double that
      of the estimate, so if we integrate with a tolerance of 1e-9:
    </p>
<pre class="programlisting"><span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">gauss_kronrod</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">&gt;::</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f1</span><span class="special">,</span> <span class="number">0</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="keyword">double</span><span class="special">&gt;::</span><span class="identifier">infinity</span><span class="special">(),</span> <span class="number">5</span><span class="special">,</span> <span class="number">1e-9</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">);</span>
</pre>
<p>
      We still achieve 1e-15 precision, with an error estimate of 1e-10.
    </p>
<h4>
<a name="math_toolkit.gauss_kronrod.h5"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.caveats"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.caveats">Caveats</a>
    </h4>
<p>
      For essentially all analytic integrands bounded on the domain, the error estimates
      provided by the routine are woefully pessimistic. In fact, in this case the
      error is very nearly equal to the error of the Gaussian quadrature formula
      of order <code class="computeroutput"><span class="special">(</span><span class="identifier">N</span><span class="special">-</span><span class="number">1</span><span class="special">)/</span><span class="number">2</span></code>, whereas the Kronrod extension converges exponentially
      beyond the Gaussian estimate. Very sophisticated method exist to estimate the
      error, but all require the integrand to lie in a particular function space.
      A more sophisticated a posteriori error estimate for an element of a particular
      function space is left to the user.
    </p>
<p>
      These routines are deliberately kept relatively simple: when they work, they
      work very well and very rapidly. However, no effort has been made to make these
      routines work well with end-point singularities or other "difficult"
      integrals. In such cases please use one of the <a class="link" href="double_exponential.html" title="Double-exponential quadrature">double-exponential
      integration schemes</a> which are generally much more robust.
    </p>
<h4>
<a name="math_toolkit.gauss_kronrod.h6"></a>
      <span class="phrase"><a name="math_toolkit.gauss_kronrod.references"></a></span><a class="link" href="gauss_kronrod.html#math_toolkit.gauss_kronrod.references">References</a>
    </h4>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
<li class="listitem">
          Kronrod, Aleksandr Semenovish (1965), <span class="emphasis"><em>Nodes and weights of quadrature
          formulas. Sixteen-place tables</em></span>, New York: Consultants Bureau
        </li>
<li class="listitem">
          Dirk P. Laurie, <span class="emphasis"><em>Calculation of Gauss-Kronrod Quadrature Rules</em></span>,
          Mathematics of Computation, Volume 66, Number 219, 1997
        </li>
<li class="listitem">
          Gonnet, Pedro, <span class="emphasis"><em>A Review of Error Estimation in Adaptive Quadrature</em></span>,
          https://arxiv.org/pdf/1003.4629.pdf
        </li>
</ul></div>
</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="gauss.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="double_exponential.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>