1[/ 2Copyright (c) 2017 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:gauss_kronrod Gauss-Kronrod Quadrature] 9 10[heading Overview] 11 12Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides an a posteriori error estimate for the integral. 13 14The idea behind Gaussian quadrature is to choose /n/ nodes and weights in such a way that polynomials of order /2n-1/ are integrated exactly. 15However, integration of polynomials is trivial, so it is rarely done via numerical methods. 16Instead, transcendental and numerically defined functions are integrated via Gaussian quadrature, and the defining problem becomes how to estimate the remainder. 17Gaussian quadrature alone (without some form of interval splitting) cannot answer this question. 18 19It is possible to compute a Gaussian quadrature of order /n/ and another of order (say) /2n+1/, and use the difference as an error estimate. 20However, this is not optimal, as the zeros of the Legendre polynomials (nodes of the Gaussian quadrature) are never the same for different orders, so /3n+1/ function evaluations must be performed. 21Kronrod 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. 22This allows an a posteriori error estimate to be provided while still preserving exponential convergence. 23Kronrod discovered that by adding /n+1/ nodes (computed from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature of order /n/, he could integrate polynomials of order /3n+1/. 24 25The integration routines provided here will perform either adaptive or non-adaptive quadrature, they should be chosen for the integration of smooth functions 26with no end point singularities. For difficult functions, or those with end point singularities, please refer to the [link math_toolkit.double_exponential double-exponential integration schemes]. 27 28 #include <boost/math/quadrature/gauss_kronrod.hpp> 29 30 template <class Real, unsigned N, class ``__Policy`` = boost::math::policies::policy<> > 31 class gauss_kronrod 32 { 33 static const RandomAccessContainer& abscissa(); 34 static const RandomAccessContainer& weights(); 35 36 template <class F> 37 static auto integrate(F f, 38 Real a, Real b, 39 unsigned max_depth = 15, 40 Real tol = tools::root_epsilon<Real>(), 41 Real* error = nullptr, 42 Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>())); 43 }; 44 45[heading Description] 46 47 static const RandomAccessContainer& abscissa(); 48 static const RandomAccessContainer& weights(); 49 50These functions provide direct access to the abscissa and weights used to perform the quadrature: the return type depends on the 51/Points/ template parameter, but is always a RandomAccessContainer type. Note that only positive (or zero) abscissa and weights 52are stored, and that they contain both the Gauss and Kronrod points. 53 54 template <class F> 55 static auto integrate(F f, 56 Real a, Real b, 57 unsigned max_depth = 15, 58 Real tol = tools::root_epsilon<Real>(), 59 Real* error = nullptr, 60 Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>())); 61 62Performs adaptive Gauss-Kronrod quadrature on function /f/ over the range (a,b). 63 64['max_depth] sets the maximum number of interval splittings permitted before stopping. Set this to zero for non-adaptive quadrature. Note that the 65algorithm descends the tree depth first, so only "difficult" areas of the integral result in interval splitting. 66 67['tol] Sets the maximum relative error in the result, this should not be set too close to machine epsilon or the function 68will simply thrash and probably not return accurate results. On the other hand the default may be overly-pessimistic. 69 70['error] When non-null, `*error` is set to the difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation. 71 72['pL1] When non-null, `*pL1` 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 73likely to be ill-conditioned. 74 75[heading Choosing the number of points] 76 77The number of points specified in the ['Points] template parameter must be an odd number: giving a (N-1)/2 Gauss quadrature as the comparison for error estimation. 78 79Internally class `gauss_kronrod` has pre-computed tables of abscissa and weights for 15, 31, 41, 51 and 61 Gauss-Kronrod points at up to 100-decimal 80digit precision. That means that using for example, `gauss_kronrod<double, 31>::integrate` incurs absolutely zero set-up overhead from 81computing the abscissa/weight pairs. When using multiprecision types with less than 100 digits of precision, then there is a small 82initial one time cost, while the abscissa/weight pairs are constructed from strings. 83 84However, for types with higher precision, or numbers of points other than those given above, the abscissa/weight pairs are computed 85when first needed and then cached for future use, which does incur a noticeable overhead. If this is likely to be an issue, then: 86 87* Defining BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND will result in a compile-time error, whenever a combination of number type 88and number of points is used which does not have pre-computed values. 89* There is a program [@../../tools/gauss_kronrod_constants.cpp gauss_kronrod_constants.cpp] which was used to provide the 90pre-computed values already in gauss_kronrod.hpp. The program can be trivially modified to generate code and constants for other precisions 91and numbers of points. 92 93[heading Complex Quadrature] 94 95The Gauss-Kronrod quadrature support integrands defined on the real line and returning complex values. 96In this case, the template argument is the real type, and the complex type is deduced via the return type of the function. 97 98[heading Examples] 99 100[import ../../example/gauss_example.cpp] 101 102[gauss_kronrod_example] 103 104[heading Caveats] 105 106For essentially all analytic integrands bounded on the domain, the error estimates provided by the routine are woefully pessimistic. 107In fact, in this case the error is very nearly equal to the error of the Gaussian quadrature formula of order `(N-1)/2`, 108whereas the Kronrod extension converges exponentially beyond the Gaussian estimate. 109Very sophisticated method exist to estimate the error, but all require the integrand to lie in a particular function space. 110A more sophisticated a posteriori error estimate for an element of a particular function space is left to the user. 111 112These routines are deliberately kept relatively simple: when they work, they work very well and very rapidly. However, no effort 113has been made to make these routines work well with end-point singularities or other "difficult" integrals. In such cases please 114use one of the [link math_toolkit.double_exponential double-exponential integration schemes] which are generally much more robust. 115 116[heading References] 117 118* Kronrod, Aleksandr Semenovish (1965), ['Nodes and weights of quadrature formulas. Sixteen-place tables], New York: Consultants Bureau 119* Dirk P. Laurie, ['Calculation of Gauss-Kronrod Quadrature Rules], Mathematics of Computation, Volume 66, Number 219, 1997 120* Gonnet, Pedro, ['A Review of Error Estimation in Adaptive Quadrature], https://arxiv.org/pdf/1003.4629.pdf 121 122[endsect] [/section:gauss_kronrod Gauss-Kronrod Quadrature] 123 124