1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Gauss-Kronrod Quadrature</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation"> 9<link rel="prev" href="gauss.html" title="Gauss-Legendre quadrature"> 10<link rel="next" href="double_exponential.html" title="Double-exponential quadrature"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="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> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.gauss_kronrod"></a><a class="link" href="gauss_kronrod.html" title="Gauss-Kronrod Quadrature">Gauss-Kronrod Quadrature</a> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.gauss_kronrod.h0"></a> 31 <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> 32 </h4> 33<p> 34 Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides 35 an a posteriori error estimate for the integral. 36 </p> 37<p> 38 The idea behind Gaussian quadrature is to choose <span class="emphasis"><em>n</em></span> nodes 39 and weights in such a way that polynomials of order <span class="emphasis"><em>2n-1</em></span> 40 are integrated exactly. However, integration of polynomials is trivial, so 41 it is rarely done via numerical methods. Instead, transcendental and numerically 42 defined functions are integrated via Gaussian quadrature, and the defining 43 problem becomes how to estimate the remainder. Gaussian quadrature alone (without 44 some form of interval splitting) cannot answer this question. 45 </p> 46<p> 47 It is possible to compute a Gaussian quadrature of order <span class="emphasis"><em>n</em></span> 48 and another of order (say) <span class="emphasis"><em>2n+1</em></span>, and use the difference 49 as an error estimate. However, this is not optimal, as the zeros of the Legendre 50 polynomials (nodes of the Gaussian quadrature) are never the same for different 51 orders, so <span class="emphasis"><em>3n+1</em></span> function evaluations must be performed. 52 Kronrod considered the problem of how to interleave nodes into a Gaussian quadrature 53 in such a way that all previous function evaluations can be reused, while increasing 54 the order of polynomials that can be integrated exactly. This allows an a posteriori 55 error estimate to be provided while still preserving exponential convergence. 56 Kronrod discovered that by adding <span class="emphasis"><em>n+1</em></span> nodes (computed 57 from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature 58 of order <span class="emphasis"><em>n</em></span>, he could integrate polynomials of order <span class="emphasis"><em>3n+1</em></span>. 59 </p> 60<p> 61 The integration routines provided here will perform either adaptive or non-adaptive 62 quadrature, they should be chosen for the integration of smooth functions with 63 no end point singularities. For difficult functions, or those with end point 64 singularities, please refer to the <a class="link" href="double_exponential.html" title="Double-exponential quadrature">double-exponential 65 integration schemes</a>. 66 </p> 67<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> 68 69<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 21. 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> 70<span class="keyword">class</span> <span class="identifier">gauss_kronrod</span> 71<span class="special">{</span> 72 <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> 73 <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> 74 75 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></span> 76 <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> 77 <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> 78 <span class="keyword">unsigned</span> <span class="identifier">max_depth</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span> 79 <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> 80 <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> 81 <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> 82<span class="special">};</span> 83</pre> 84<h4> 85<a name="math_toolkit.gauss_kronrod.h1"></a> 86 <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> 87 </h4> 88<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> 89<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> 90</pre> 91<p> 92 These functions provide direct access to the abscissa and weights used to perform 93 the quadrature: the return type depends on the <span class="emphasis"><em>Points</em></span> 94 template parameter, but is always a RandomAccessContainer type. Note that only 95 positive (or zero) abscissa and weights are stored, and that they contain both 96 the Gauss and Kronrod points. 97 </p> 98<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> 99<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> 100 <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> 101 <span class="keyword">unsigned</span> <span class="identifier">max_depth</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span> 102 <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> 103 <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> 104 <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> 105</pre> 106<p> 107 Performs adaptive Gauss-Kronrod quadrature on function <span class="emphasis"><em>f</em></span> 108 over the range (a,b). 109 </p> 110<p> 111 <span class="emphasis"><em>max_depth</em></span> sets the maximum number of interval splittings 112 permitted before stopping. Set this to zero for non-adaptive quadrature. Note 113 that the algorithm descends the tree depth first, so only "difficult" 114 areas of the integral result in interval splitting. 115 </p> 116<p> 117 <span class="emphasis"><em>tol</em></span> Sets the maximum relative error in the result, this 118 should not be set too close to machine epsilon or the function will simply 119 thrash and probably not return accurate results. On the other hand the default 120 may be overly-pessimistic. 121 </p> 122<p> 123 <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 124 (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation. 125 </p> 126<p> 127 <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, 128 if there is a significant difference between this and the returned value, then 129 the result is likely to be ill-conditioned. 130 </p> 131<h4> 132<a name="math_toolkit.gauss_kronrod.h2"></a> 133 <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 134 the number of points</a> 135 </h4> 136<p> 137 The number of points specified in the <span class="emphasis"><em>Points</em></span> template 138 parameter must be an odd number: giving a (N-1)/2 Gauss quadrature as the comparison 139 for error estimation. 140 </p> 141<p> 142 Internally class <code class="computeroutput"><span class="identifier">gauss_kronrod</span></code> 143 has pre-computed tables of abscissa and weights for 15, 31, 41, 51 and 61 Gauss-Kronrod 144 points at up to 100-decimal digit precision. That means that using for example, 145 <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> 146 incurs absolutely zero set-up overhead from computing the abscissa/weight pairs. 147 When using multiprecision types with less than 100 digits of precision, then 148 there is a small initial one time cost, while the abscissa/weight pairs are 149 constructed from strings. 150 </p> 151<p> 152 However, for types with higher precision, or numbers of points other than those 153 given above, the abscissa/weight pairs are computed when first needed and then 154 cached for future use, which does incur a noticeable overhead. If this is likely 155 to be an issue, then: 156 </p> 157<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 158<li class="listitem"> 159 Defining BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND will result in a compile-time 160 error, whenever a combination of number type and number of points is used 161 which does not have pre-computed values. 162 </li> 163<li class="listitem"> 164 There is a program <a href="../../../tools/gauss_kronrod_constants.cpp" target="_top">gauss_kronrod_constants.cpp</a> 165 which was used to provide the pre-computed values already in gauss_kronrod.hpp. 166 The program can be trivially modified to generate code and constants for 167 other precisions and numbers of points. 168 </li> 169</ul></div> 170<h4> 171<a name="math_toolkit.gauss_kronrod.h3"></a> 172 <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 173 Quadrature</a> 174 </h4> 175<p> 176 The Gauss-Kronrod quadrature support integrands defined on the real line and 177 returning complex values. In this case, the template argument is the real type, 178 and the complex type is deduced via the return type of the function. 179 </p> 180<h4> 181<a name="math_toolkit.gauss_kronrod.h4"></a> 182 <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> 183 </h4> 184<p> 185 We'll begin by integrating exp(-t<sup>2</sup>/2) over (0,+∞) using a 7 term Gauss rule 186 and 15 term Kronrod rule, and begin by defining the function to integrate as 187 a C++ lambda expression: 188 </p> 189<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> 190 191<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> 192</pre> 193<p> 194 We'll start off with a one shot (ie non-adaptive) integration, and keep track 195 of the estimated error: 196 </p> 197<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error</span><span class="special">;</span> 198<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> 199</pre> 200<p> 201 This yields Q = 1.25348207361, which has an absolute error of 1e-4 compared 202 to the estimated error of 5e-3: this is fairly typical, with the difference 203 between Gauss and Gauss-Kronrod schemes being much higher than the actual error. 204 Before moving on to adaptive quadrature, lets try again with more points, in 205 fact with the largest Gauss-Kronrod scheme we have cached (30/61): 206 </p> 207<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> 208</pre> 209<p> 210 This yields an absolute error of 3e-15 against an estimate of 1e-8, which is 211 about as good as we're going to get at double precision 212 </p> 213<p> 214 However, instead of continuing with ever more points, lets switch to adaptive 215 integration, and set the desired relative error to 1e-14 against a maximum 216 depth of 5: 217 </p> 218<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> 219</pre> 220<p> 221 This yields an actual error of zero, against an estimate of 4e-15. In fact 222 in this case the requested tolerance was almost certainly set too low: as we've 223 seen above, for smooth functions, the precision achieved is often double that 224 of the estimate, so if we integrate with a tolerance of 1e-9: 225 </p> 226<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> 227</pre> 228<p> 229 We still achieve 1e-15 precision, with an error estimate of 1e-10. 230 </p> 231<h4> 232<a name="math_toolkit.gauss_kronrod.h5"></a> 233 <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> 234 </h4> 235<p> 236 For essentially all analytic integrands bounded on the domain, the error estimates 237 provided by the routine are woefully pessimistic. In fact, in this case the 238 error is very nearly equal to the error of the Gaussian quadrature formula 239 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 240 beyond the Gaussian estimate. Very sophisticated method exist to estimate the 241 error, but all require the integrand to lie in a particular function space. 242 A more sophisticated a posteriori error estimate for an element of a particular 243 function space is left to the user. 244 </p> 245<p> 246 These routines are deliberately kept relatively simple: when they work, they 247 work very well and very rapidly. However, no effort has been made to make these 248 routines work well with end-point singularities or other "difficult" 249 integrals. In such cases please use one of the <a class="link" href="double_exponential.html" title="Double-exponential quadrature">double-exponential 250 integration schemes</a> which are generally much more robust. 251 </p> 252<h4> 253<a name="math_toolkit.gauss_kronrod.h6"></a> 254 <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> 255 </h4> 256<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 257<li class="listitem"> 258 Kronrod, Aleksandr Semenovish (1965), <span class="emphasis"><em>Nodes and weights of quadrature 259 formulas. Sixteen-place tables</em></span>, New York: Consultants Bureau 260 </li> 261<li class="listitem"> 262 Dirk P. Laurie, <span class="emphasis"><em>Calculation of Gauss-Kronrod Quadrature Rules</em></span>, 263 Mathematics of Computation, Volume 66, Number 219, 1997 264 </li> 265<li class="listitem"> 266 Gonnet, Pedro, <span class="emphasis"><em>A Review of Error Estimation in Adaptive Quadrature</em></span>, 267 https://arxiv.org/pdf/1003.4629.pdf 268 </li> 269</ul></div> 270</div> 271<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 272<td align="left"></td> 273<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 274 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 275 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 276 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 277 Daryle Walker and Xiaogang Zhang<p> 278 Distributed under the Boost Software License, Version 1.0. (See accompanying 279 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>) 280 </p> 281</div></td> 282</tr></table> 283<hr> 284<div class="spirit-nav"> 285<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> 286</div> 287</body> 288</html> 289
