libs/math/doc/quadrature/trapezoidal.qbk

   1 [/
   2 Copyright (c) 2017 Nick Thompson
   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:trapezoidal Trapezoidal Quadrature]
   9
  10 [heading Synopsis]
  11
  12     #include <boost/math/quadrature/trapezoidal.hpp>
  13     namespace boost{ namespace math{ namespace quadrature {
  14
  15     template<class F, class Real>
  16     auto trapezoidal(F f, Real a, Real b,
  17                      Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
  18                      size_t max_refinements = 12,
  19                      Real* error_estimate = nullptr,
  20                      Real* L1 = nullptr);
  21
  22     template<class F, class Real, class ``__Policy``>
  23     auto trapezoidal(F f, Real a, Real b, Real tol, size_t max_refinements,
  24                      Real* error_estimate, Real* L1, const ``__Policy``& pol);
  25
  26     }}} // namespaces
  27
  28 [heading Description]
  29
  30 The functional `trapezoidal` calculates the integral of a function /f/ using the surprisingly simple trapezoidal rule.
  31 If we assume only that the integrand is twice continuously differentiable,
  32 we can prove that the error of the composite trapezoidal rule
  33 is [bigo](h[super 2]).
  34 Hence halving the interval only cuts the error by about a fourth,
  35 which in turn implies that we must evaluate the function many times before an acceptable accuracy can be achieved.
  36
  37 However, the trapezoidal rule has an astonishing property:
  38 If the integrand is periodic, and we integrate it over a period,
  39 then the trapezoidal rule converges faster than any power of the step size /h/.
  40 This can be seen by examination of the [@https://en.wikipedia.org/wiki/Euler-Maclaurin_formula Euler-Maclaurin summation formula],
  41 which relates a definite integral to its trapezoidal sum and error terms proportional to the derivatives of the function at the endpoints and the Bernoulli numbers.
  42 If the derivatives at the endpoints are the same or vanish, then the error very nearly vanishes.
  43 Hence the trapezoidal rule is essentially optimal for periodic integrands.
  44
  45 Other classes of integrands which are integrated efficiently by this method are the C[sub 0][super \u221E](\u221D) [@https://en.wikipedia.org/wiki/Bump_function bump functions] and bell-shaped integrals over the infinite interval.
  46 For details, see [@http://epubs.siam.org/doi/pdf/10.1137/130932132 Trefethen's] SIAM review.
  47
  48
  49 In its simplest form, an integration can be performed by the following code
  50
  51     using boost::math::quadrature::trapezoidal;
  52     auto f = [](double x) { return 1/(5 - 4*cos(x)); };
  53     double I = trapezoidal(f, 0.0, boost::math::constants::two_pi<double>());
  54
  55 The integrand must accept a real number argument, but can return a complex number.
  56 This is useful for contour integrals (which are manifestly periodic) and high-order numerical differentiation of analytic functions.
  57 An example using the integral definition of the complex Bessel function is shown here:
  58
  59     auto bessel_integrand = [&n, &z](double theta)->std::complex<double>
  60     {
  61         std::complex<double> z{2, 3};
  62         using std::cos;
  63         using std::sin;
  64         return cos(z*sin(theta) - 2*theta)/pi<double>();
  65     };
  66
  67     using boost::math::quadrature::trapezoidal;
  68     std::complex<double> Jnz = trapezoidal(bessel_integrand, 0.0, pi<Real>());
  69
  70 Other special functions which are efficiently evaluated in the complex plane by trapezoidal quadrature are modified Bessel functions and the complementary error function. Another application of complex-valued trapezoidal quadrature is computation of high-order numerical derivatives; see Lyness and Moler for details.
  71
  72
  73 Since the routine is adaptive, step sizes are halved continuously until a tolerance is reached.
  74 In order to control this tolerance, simply call the routine with an additional argument
  75
  76     double I = trapezoidal(f, 0.0, two_pi<double>(), 1e-6);
  77
  78 The routine stops when successive estimates of the integral `I1` and `I0` differ by less than the tolerance multiplied by the estimated L[sub 1] norm of the function.
  79 A good choice for the tolerance is [radic][epsilon], which is the default.
  80 If the integrand is periodic, then the number of correct digits should double on each interval halving.
  81 Hence, once the integration routine has estimated that the error is [radic][epsilon], then the actual error should be ~[epsilon].
  82 If the integrand is *not* periodic, then reducing the error to [radic][epsilon] takes much longer, but is nonetheless possible without becoming a major performance bug.
  83
  84 A question arises as to what to do when successive estimates never pass below the tolerance threshold.
  85 The stepsize would be halved repeatedly, generating an exponential explosion in function evaluations.
  86 As such, you may pass an optional argument `max_refinements` which controls how many times the interval may be halved before giving up.
  87 By default, this maximum number of refinement steps is 12,
  88 which should never be hit in double precision unless the function is not periodic.
  89 However, for higher-precision types,
  90 it may be of interest to allow the algorithm to compute more refinements:
  91
  92     size_t max_refinements = 15;
  93     long double I = trapezoidal(f, 0.0L, two_pi<long double>(), 1e-9L, max_refinements);
  94
  95 Note that the maximum allowed compute time grows exponentially with `max_refinements`.
  96 The routine will not throw an exception if the maximum refinements is achieved without the requested tolerance being attained.
  97 This is because the value calculated is more often than not still usable.
  98 However, for applications with high-reliability requirements,
  99 the error estimate should be queried.
 100 This is achieved by passing additional pointers into the routine:
 101
 102     double error_estimate;
 103     double L1;
 104     double I = trapezoidal(f, 0.0, two_pi<double>(), tolerance, max_refinements, &error_estimate, &L1);
 105     if (error_estimate > tolerance*L1)
 106     {
 107          double I = some_other_quadrature_method(f, 0, two_pi<double>());
 108     }
 109
 110 The final argument is the L[sub 1] norm of the integral.
 111 This is computed along with the integral, and is an essential component of the algorithm.
 112 First, the L[sub 1] norm establishes a scale against which the error can be measured.
 113 Second, the L[sub 1] norm can be used to evaluate the stability of the computation.
 114 This can be formulated in a rigorous manner by defining the *condition number of summation*.
 115 The condition number of summation is defined by
 116
 117 [expression ['[kappa](S[sub n]) := [Sigma][sub i][super n] |x[sub i]|/|[Sigma][sub i][super n] x[sub i]|]]
 118
 119 If this number of ~10[super k],
 120 then /k/ additional digits are expected to be lost in addition to digits lost due to floating point rounding error.
 121 As all numerical quadrature methods reduce to summation,
 122 their stability is controlled by the ratio \u222B |f| dt/|\u222B f dt |,
 123 which is easily seen to be equivalent to condition number of summation when evaluated numerically.
 124 Hence both the error estimate and the condition number of summation should be analyzed in applications requiring very high precision and reliability.
 125
 126 As an example, we consider evaluation of Bessel functions by trapezoidal quadrature.
 127 The Bessel function of the first kind is defined via
 128
 129 [expression ['J[sub n](x) = 1/2\u03A0 \u222B[sub -\u03A0][super \u03A0] cos(n t - x sin(t)) dt]]
 130
 131 The integrand is periodic, so the Euler-Maclaurin summation formula guarantees exponential convergence via the trapezoidal quadrature.
 132 Without careful consideration, it seems this would be a very attractive method to compute Bessel functions.
 133 However, we see that for large /n/ the integrand oscillates rapidly,
 134 taking on positive and negative values,
 135 and hence the trapezoidal sums become ill-conditioned.
 136 In double precision, /x = 17/ and /n = 25/ gives a sum which is so poorly conditioned that zero correct digits are obtained.
 137
 138 [optional_policy]
 139
 140 References:
 141
 142 Trefethen, Lloyd N., Weideman, J.A.C., ['The Exponentially Convergent Trapezoidal Rule], SIAM Review, Vol. 56, No. 3, 2014.
 143
 144 Stoer, Josef, and Roland Bulirsch. ['Introduction to numerical analysis. Vol. 12.], Springer Science & Business Media, 2013.
 145
 146 Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Society for industrial and applied mathematics, 2002.
 147
 148 Lyness, James N., and Cleve B. Moler. ['Numerical differentiation of analytic functions.] SIAM Journal on Numerical Analysis 4.2 (1967): 202-210.
 149
 150 Gil, Amparo, Javier Segura, and Nico M. Temme. ['Computing special functions by using quadrature rules.] Numerical Algorithms 33.1-4 (2003): 265-275.
 151 [endsect] [/section:trapezoidal Trapezoidal Quadrature]
 152