126 lines · plain
1/*2 * Copyright Nick Thompson, 20173 * Use, modification and distribution are subject to the4 * Boost Software License, Version 1.0. (See accompanying file5 * 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 bounded11 * over the whole interval, then this routine will still converge, albeit slowly. However, there12 * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.13 */14 15#ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP16#define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP17 18#include <cmath>19#include <limits>20#include <utility>21#include <stdexcept>22#include <boost/math/constants/constants.hpp>23#include <boost/math/special_functions/fpclassify.hpp>24#include <boost/math/policies/error_handling.hpp>25#include <boost/math/tools/cxx03_warn.hpp>26 27namespace boost{ namespace math{ namespace quadrature {28 29template<class F, class Real, class Policy>30auto 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>()))31{32 static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)";33 using std::abs;34 using boost::math::constants::half;35 // In many math texts, K represents the field of real or complex numbers.36 // Too bad we can't put blackboard bold into C++ source!37 typedef decltype(f(a)) K;38 static_assert(!std::is_integral<K>::value,39 "The return type cannot be integral, it must be either a real or complex floating point type.");40 if (!(boost::math::isfinite)(a))41 {42 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));43 }44 if (!(boost::math::isfinite)(b))45 {46 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));47 }48 49 if (a == b)50 {51 return static_cast<K>(0);52 }53 if(a > b)54 {55 return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol);56 }57 58 59 K ya = f(a);60 K yb = f(b);61 Real h = (b - a)*half<Real>();62 K I0 = (ya + yb)*h;63 Real IL0 = (abs(ya) + abs(yb))*h;64 65 K yh = f(a + h);66 K I1;67 I1 = I0*half<Real>() + yh*h;68 Real IL1 = IL0*half<Real>() + abs(yh)*h;69 70 // The recursion is:71 // 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)72 std::size_t k = 2;73 // We want to go through at least 5 levels so we have sampled the function at least 20 times.74 // Otherwise, we could terminate prematurely and miss essential features.75 // This is of course possible anyway, but 20 samples seems to be a reasonable compromise.76 Real error = abs(I0 - I1);77 // I take k < 5, rather than k < 4, or some other smaller minimum number,78 // because I hit a truly exceptional bug where the k = 2 and k =3 refinement were bitwise equal,79 // but the quadrature had not yet converged.80 while (k < 5 || (k < max_refinements && error > tol*IL1) )81 {82 I0 = I1;83 IL0 = IL1;84 85 I1 = I0*half<Real>();86 IL1 = IL0*half<Real>();87 std::size_t p = static_cast<std::size_t>(1u) << k;88 h *= half<Real>();89 K sum = 0;90 Real absum = 0;91 92 for(std::size_t j = 1; j < p; j += 2)93 {94 K y = f(a + j*h);95 sum += y;96 absum += abs(y);97 }98 99 I1 += sum*h;100 IL1 += absum*h;101 ++k;102 error = abs(I0 - I1);103 }104 105 if (error_estimate)106 {107 *error_estimate = error;108 }109 110 if (L1)111 {112 *L1 = IL1;113 }114 115 return static_cast<K>(I1);116}117 118template<class F, class Real>119auto 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>()))120{121 return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>());122}123 124}}}125#endif126