math/quadrature/trapezoidal.hpp

0001 /*
0002  * Copyright Nick Thompson, 2017
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  *
0007  * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period,
0008  * or to integrate a function whose derivative vanishes at the endpoints.
0009  *
0010  * If your function does not satisfy these conditions, and instead is simply continuous and bounded
0011  * over the whole interval, then this routine will still converge, albeit slowly. However, there
0012  * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.
0013  */
0014
0015 #ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
0016 #define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
0017
0018 #include <cmath>
0019 #include <limits>
0020 #include <utility>
0021 #include <stdexcept>
0022 #include <boost/math/constants/constants.hpp>
0023 #include <boost/math/special_functions/fpclassify.hpp>
0024 #include <boost/math/policies/error_handling.hpp>
0025 #include <boost/math/tools/cxx03_warn.hpp>
0026
0027 namespace boost{ namespace math{ namespace quadrature {
0028
0029 template<class F, class Real, class Policy>
0030 auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Real* error_estimate, Real* L1, const Policy& pol)->decltype(std::declval<F>()(std::declval<Real>()))
0031 {
0032     static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)";
0033     using std::abs;
0034     using boost::math::constants::half;
0035     // In many math texts, K represents the field of real or complex numbers.
0036     // Too bad we can't put blackboard bold into C++ source!
0037     typedef decltype(f(a)) K;
0038     static_assert(!std::is_integral<K>::value,
0039                   "The return type cannot be integral, it must be either a real or complex floating point type.");
0040     if (!(boost::math::isfinite)(a))
0041     {
0042        return static_cast<K>(boost::math::policies::raise_domain_error(function, "Left endpoint of integration must be finite for adaptive trapezoidal integration but got a = %1%.\n", a, pol));
0043     }
0044     if (!(boost::math::isfinite)(b))
0045     {
0046        return static_cast<K>(boost::math::policies::raise_domain_error(function, "Right endpoint of integration must be finite for adaptive trapezoidal integration but got b = %1%.\n", b, pol));
0047     }
0048
0049     if (a == b)
0050     {
0051         return static_cast<K>(0);
0052     }
0053     if(a > b)
0054     {
0055         return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol);
0056     }
0057
0058
0059     K ya = f(a);
0060     K yb = f(b);
0061     Real h = (b - a)*half<Real>();
0062     K I0 = (ya + yb)*h;
0063     Real IL0 = (abs(ya) + abs(yb))*h;
0064
0065     K yh = f(a + h);
0066     K I1;
0067     I1 = I0*half<Real>() + yh*h;
0068     Real IL1 = IL0*half<Real>() + abs(yh)*h;
0069
0070     // The recursion is:
0071     // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k)
0072     std::size_t k = 2;
0073     // We want to go through at least 5 levels so we have sampled the function at least 20 times.
0074     // Otherwise, we could terminate prematurely and miss essential features.
0075     // This is of course possible anyway, but 20 samples seems to be a reasonable compromise.
0076     Real error = abs(I0 - I1);
0077     // I take k < 5, rather than k < 4, or some other smaller minimum number,
0078     // because I hit a truly exceptional bug where the k = 2 and k =3 refinement were bitwise equal,
0079     // but the quadrature had not yet converged.
0080     while (k < 5 || (k < max_refinements && error > tol*IL1) )
0081     {
0082         I0 = I1;
0083         IL0 = IL1;
0084
0085         I1 = I0*half<Real>();
0086         IL1 = IL0*half<Real>();
0087         std::size_t p = static_cast<std::size_t>(1u) << k;
0088         h *= half<Real>();
0089         K sum = 0;
0090         Real absum = 0;
0091
0092         for(std::size_t j = 1; j < p; j += 2)
0093         {
0094             K y = f(a + j*h);
0095             sum += y;
0096             absum += abs(y);
0097         }
0098
0099         I1 += sum*h;
0100         IL1 += absum*h;
0101         ++k;
0102         error = abs(I0 - I1);
0103     }
0104
0105     if (error_estimate)
0106     {
0107         *error_estimate = error;
0108     }
0109
0110     if (L1)
0111     {
0112         *L1 = IL1;
0113     }
0114
0115     return static_cast<K>(I1);
0116 }
0117
0118 template<class F, class Real>
0119 auto trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), std::size_t max_refinements = 12, Real* error_estimate = nullptr, Real* L1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
0120 {
0121    return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>());
0122 }
0123
0124 }}}
0125 #endif