townforge/math
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
cc
math/doc/html/math_toolkit/gauss_kronrod.html
jzmaddock 72e469da0a Update version number in docs. 
[CI SKIP]
2019-10-31 17:55:35 +00:00

289 lines
26 KiB
HTML
Raw Permalink Blame History

<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>
Reference in New Issue View Git Blame Copy Permalink