1 /*
2 * Copyright Nick Thompson, 2017
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 * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period,
8 * or to integrate a function whose derivative vanishes at the endpoints.
9 *
10 * If your function does not satisfy these conditions, and instead is simply continuous and bounded
11 * over the whole interval, then this routine will still converge, albeit slowly. However, there
12 * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.
13 */
14
15 #ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
16 #define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
17
18 #include <cmath>
19 #include <limits>
20 #include <stdexcept>
21 #include <boost/math/constants/constants.hpp>
22 #include <boost/math/special_functions/fpclassify.hpp>
23 #include <boost/math/policies/error_handling.hpp>
24 #include <boost/math/tools/cxx03_warn.hpp>
25
26 namespace boost{ namespace math{ namespace quadrature {
27
28 template<class F, class Real, class Policy>
trapezoidal(F f,Real a,Real b,Real tol,std::size_t max_refinements,Real * error_estimate,Real * L1,const Policy & pol)29 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>()))
30 {
31 static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)";
32 using std::abs;
33 using boost::math::constants::half;
34 // In many math texts, K represents the field of real or complex numbers.
35 // Too bad we can't put blackboard bold into C++ source!
36 typedef decltype(f(a)) K;
37 if (!(boost::math::isfinite)(a))
38 {
39 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));
40 }
41 if (!(boost::math::isfinite)(b))
42 {
43 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));
44 }
45
46 if (a == b)
47 {
48 return static_cast<K>(0);
49 }
50 if(a > b)
51 {
52 return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol);
53 }
54
55
56 K ya = f(a);
57 K yb = f(b);
58 Real h = (b - a)*half<Real>();
59 K I0 = (ya + yb)*h;
60 Real IL0 = (abs(ya) + abs(yb))*h;
61
62 K yh = f(a + h);
63 K I1;
64 I1 = I0*half<Real>() + yh*h;
65 Real IL1 = IL0*half<Real>() + abs(yh)*h;
66
67 // The recursion is:
68 // 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)
69 std::size_t k = 2;
70 // We want to go through at least 5 levels so we have sampled the function at least 20 times.
71 // Otherwise, we could terminate prematurely and miss essential features.
72 // This is of course possible anyway, but 20 samples seems to be a reasonable compromise.
73 Real error = abs(I0 - I1);
74 // I take k < 5, rather than k < 4, or some other smaller minimum number,
75 // because I hit a truly exceptional bug where the k = 2 and k =3 refinement were bitwise equal,
76 // but the quadrature had not yet converged.
77 while (k < 5 || (k < max_refinements && error > tol*IL1) )
78 {
79 I0 = I1;
80 IL0 = IL1;
81
82 I1 = I0*half<Real>();
83 IL1 = IL0*half<Real>();
84 std::size_t p = static_cast<std::size_t>(1u) << k;
85 h *= half<Real>();
86 K sum = 0;
87 Real absum = 0;
88
89 for(std::size_t j = 1; j < p; j += 2)
90 {
91 K y = f(a + j*h);
92 sum += y;
93 absum += abs(y);
94 }
95
96 I1 += sum*h;
97 IL1 += absum*h;
98 ++k;
99 error = abs(I0 - I1);
100 }
101
102 if (error_estimate)
103 {
104 *error_estimate = error;
105 }
106
107 if (L1)
108 {
109 *L1 = IL1;
110 }
111
112 return static_cast<K>(I1);
113 }
114 #if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
115 // Template argument deduction failure otherwise:
116 template<class F, class Real>
trapezoidal(F f,Real a,Real b,Real tol=0,std::size_t max_refinements=12,Real * error_estimate=0,Real * L1=0)117 auto trapezoidal(F f, Real a, Real b, Real tol = 0, std::size_t max_refinements = 12, Real* error_estimate = 0, Real* L1 = 0)->decltype(std::declval<F>()(std::declval<Real>()))
118 #elif !defined(BOOST_NO_CXX11_NULLPTR)
119 template<class F, class Real>
120 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>()))
121 #else
122 template<class F, class Real>
123 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 = 0, Real* L1 = 0)->decltype(std::declval<F>()(std::declval<Real>()))
124 #endif
125 {
126 #if BOOST_WORKAROUND(BOOST_MSVC, <= 1600)
127 if (tol == 0)
128 tol = boost::math::tools::root_epsilon<Real>();
129 #endif
130 return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>());
131 }
132
133 }}}
134 #endif
135