123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123 |
- [/
- Copyright (c) 2017 Nick Thompson
- Use, modification and distribution are subject to the
- Boost Software License, Version 1.0. (See accompanying file
- LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- ]
- [section:gauss_kronrod Gauss-Kronrod Quadrature]
- [heading Overview]
- Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides an a posteriori error estimate for the integral.
- The idea behind Gaussian quadrature is to choose /n/ nodes and weights in such a way that polynomials of order /2n-1/ 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.
- It 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.
- 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 /3n+1/ 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 /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/.
- 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 [link math_toolkit.double_exponential double-exponential integration schemes].
- #include <boost/math/quadrature/gauss_kronrod.hpp>
- template <class Real, unsigned N, class ``__Policy`` = boost::math::policies::policy<> >
- class gauss_kronrod
- {
- static const RandomAccessContainer& abscissa();
- static const RandomAccessContainer& weights();
- template <class F>
- static auto integrate(F f,
- Real a, Real b,
- unsigned max_depth = 15,
- Real tol = tools::root_epsilon<Real>(),
- Real* error = nullptr,
- Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
- };
- [heading Description]
- static const RandomAccessContainer& abscissa();
- static const RandomAccessContainer& weights();
- These functions provide direct access to the abscissa and weights used to perform the quadrature: the return type depends on the
- /Points/ 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.
- template <class F>
- static auto integrate(F f,
- Real a, Real b,
- unsigned max_depth = 15,
- Real tol = tools::root_epsilon<Real>(),
- Real* error = nullptr,
- Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
- Performs adaptive Gauss-Kronrod quadrature on function /f/ over the range (a,b).
- ['max_depth] 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.
- ['tol] 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.
- ['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.
- ['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
- likely to be ill-conditioned.
- [heading Choosing the number of points]
- The 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.
- Internally 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
- digit precision. That means that using for example, `gauss_kronrod<double, 31>::integrate` 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.
- 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:
- * 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.
- * There is a program [@../../tools/gauss_kronrod_constants.cpp gauss_kronrod_constants.cpp] 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.
- [heading Complex Quadrature]
- 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.
- [heading Examples]
- [import ../../example/gauss_example.cpp]
- [gauss_kronrod_example]
- [heading Caveats]
- 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 `(N-1)/2`,
- 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.
- 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 [link math_toolkit.double_exponential double-exponential integration schemes] which are generally much more robust.
- [heading References]
- * Kronrod, Aleksandr Semenovish (1965), ['Nodes and weights of quadrature formulas. Sixteen-place tables], New York: Consultants Bureau
- * Dirk P. Laurie, ['Calculation of Gauss-Kronrod Quadrature Rules], Mathematics of Computation, Volume 66, Number 219, 1997
- * Gonnet, Pedro, ['A Review of Error Estimation in Adaptive Quadrature], https://arxiv.org/pdf/1003.4629.pdf
- [endsect] [/section:gauss_kronrod Gauss-Kronrod Quadrature]
|