libs/math/doc/quadrature/gauss.qbk

   1 [/
   2 Copyright (c) 2017 John Maddock
   3 Use, modification and distribution are subject to the
   4 Boost Software License, Version 1.0. (See accompanying file
   5 LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 ]
   7
   8 [section:gauss Gauss-Legendre quadrature]
   9
  10 [heading Synopsis]
  11
  12 `#include <boost/math/quadrature/gauss.hpp>`
  13
  14    namespace boost{ namespace math{ namespace quadrature{
  15
  16    template <class Real, unsigned Points, class ``__Policy`` = boost::math::policies::policy<> >
  17    struct gauss
  18    {
  19       static const RandomAccessContainer& abscissa();
  20       static const RandomAccessContainer& weights();
  21
  22       template <class F>
  23       static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
  24
  25       template <class F>
  26       static auto integrate(F f, Real a, Real b, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
  27    };
  28
  29    }}} // namespaces
  30
  31 [heading description]
  32
  33 The `gauss` class template performs "one shot" non-adaptive Gauss-Legendre integration on some arbitrary function /f/ using the
  34 number of evaluation points as specified by /Points/.
  35
  36 This is intentionally a very simple quadrature routine, it obtains no estimate of the error, and is not adaptive, but is very efficient
  37 in simple cases that involve integrating smooth "bell like" functions and functions with rapidly convergent power series.
  38
  39    static const RandomAccessContainer& abscissa();
  40    static const RandomAccessContainer& weights();
  41
  42 These functions provide direct access to the abscissa and weights used to perform the quadrature: the return type depends on the
  43 /Points/ template parameter, but is always a RandomAccessContainer type.  Note that only positive (or zero) abscissa and weights
  44 are stored.
  45
  46     template <class F>
  47     static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
  48
  49 Integrates /f/ over (-1,1), and optionally sets `*pL1` to the L1 norm of the returned value: if this is substantially larger
  50 than the return value, then the sum was ill-conditioned.  Note however, that no error estimate is available.
  51
  52     template <class F>
  53     static auto integrate(F f, Real a, Real b, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
  54
  55 Integrates /f/ over (a,b), and optionally sets `*pL1` to the L1 norm of the returned value: if this is substantially larger
  56 than the return value, then the sum was ill-conditioned.  Note however, that no error estimate is available.  This function
  57 supports both finite and infinite /a/ and /b/, as long as `a < b`.
  58
  59 The Gaussian quadrature routine support both real and complex-valued quadrature.
  60 For example, the Lambert-W function admits the integral representation
  61
  62 [expression ['W(z) = 1/2\u03A0 \u222B[sub -\u03A0][super \u03A0]  ((1- v cot(v) )^2 + v^2)/(z + v csc(v) exp(-v cot(v))) dv]]
  63
  64 so it can be effectively computed via Gaussian quadrature using the following code:
  65
  66     Complex z{2, 3};
  67     auto lw = [&z](Real v)->Complex {
  68       using std::cos;
  69       using std::sin;
  70       using std::exp;
  71       Real sinv = sin(v);
  72       Real cosv = cos(v);
  73
  74       Real cotv = cosv/sinv;
  75       Real cscv = 1/sinv;
  76       Real t = (1-v*cotv)*(1-v*cotv) + v*v;
  77       Real x = v*cscv*exp(-v*cotv);
  78       Complex den = z + x;
  79       Complex num = t*(z/pi<Real>());
  80       Complex res = num/den;
  81       return res;
  82     };
  83
  84     boost::math::quadrature::gauss<Real, 30> integrator;
  85     Complex W = integrator.integrate(lw, (Real) 0, pi<Real>());
  86
  87
  88 [heading Choosing the number of points]
  89
  90 Internally class `gauss` has pre-computed tables of abscissa and weights for 7, 15, 20, 25 and 30 points at up to 100-decimal
  91 digit precision.  That means that using for example, `gauss<double, 30>::integrate` incurs absolutely zero set-up overhead from
  92 computing the abscissa/weight pairs.  When using multiprecision types with less than 100 digits of precision, then there is a small
  93 initial one time cost, while the abscissa/weight pairs are constructed from strings.
  94
  95 However, for types with higher precision, or numbers of points other than those given above, the abscissa/weight pairs are computed
  96 when first needed and then cached for future use, which does incur a noticeable overhead.  If this is likely to be an issue, then
  97
  98 * Defining BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND will result in a compile-time error, whenever a combination of number type
  99 and number of points is used which does not have pre-computed values.
 100 * There is a program [@../../tools/gauss_kronrod_constants.cpp gauss_kronrod_constants.cpp] which was used to provide the
 101 pre-computed values already in gauss.hpp.  The program can be trivially modified to generate code and constants for other precisions
 102 and numbers of points.
 103
 104 [heading Examples]
 105
 106 [import ../../example/gauss_example.cpp]
 107
 108 [gauss_example]
 109
 110 [endsect] [/section:gauss Gauss-Legendre quadrature]
 111