EQ2EMu/EQ2/source/depends/boost_1_72_0/libs/math/doc/quadrature/gauss_kronrod.qbk

[/
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]