1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<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 4.2.1">
<link rel="up" href="../quadrature.html" title="Chapter 14. Quadrature and Differentiation">
<link rel="prev" href="gauss.html" title="Gauss-Legendre quadrature">
<link rel="next" href="double_exponential.html" title="Double-exponential quadrature">
<meta name="viewport" content="width=device-width, initial-scale=1">
</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"><</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">></span>
<span class="keyword">template</span> <span class="special"><</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 22. 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"><></span> <span class="special">></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">&</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">&</span> <span class="identifier">weights</span><span class="special">();</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></span>
<span class="keyword">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"><</span><span class="identifier">Real</span><span class="special">>(),</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">)-></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"><</span><span class="identifier">F</span><span class="special">>()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>()));</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">&</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">&</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"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></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"><</span><span class="identifier">Real</span><span class="special">>(),</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">)-></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"><</span><span class="identifier">F</span><span class="special">>()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>()));</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-pessimistic.
</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"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">31</span><span class="special">>::</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,+∞) 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"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">>::</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"><</span><span class="keyword">double</span><span class="special">>::</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">&</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"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">61</span><span class="special">>::</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"><</span><span class="keyword">double</span><span class="special">>::</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">&</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"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">>::</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"><</span><span class="keyword">double</span><span class="special">>::</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">&</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"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">15</span><span class="special">>::</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"><</span><span class="keyword">double</span><span class="special">>::</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">&</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>
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Rå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>
<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>
|
