math/quadrature/tanh_sinh.hpp

0001 // Copyright Nick Thompson, 2017
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 /*
0008  * This class performs tanh-sinh quadrature on the real line.
0009  * Tanh-sinh quadrature is exponentially convergent for integrands in Hardy spaces,
0010  * (see https://en.wikipedia.org/wiki/Hardy_space for a formal definition), and is optimal for a random function from that class.
0011  *
0012  * The tanh-sinh quadrature is one of a class of so called "double exponential quadratures"-there is a large family of them,
0013  * but this one seems to be the most commonly used.
0014  *
0015  * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk,
0016  * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not
0017  * require the function to be holomorphic, only differentiable up to some order.
0018  *
0019  * In addition, if you are integrating a periodic function over a period, the trapezoidal rule is better.
0020  *
0021  * References:
0022  *
0023  * 1) Mori, Masatake. "Quadrature formulas obtained by variable transformation and the DE-rule." Journal of Computational and Applied Mathematics 12 (1985): 119-130.
0024  * 2) Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. "A comparison of three high-precision quadrature schemes." Experimental Mathematics 14.3 (2005): 317-329.
0025  * 3) Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
0026  *
0027  */
0028
0029 #ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP
0030 #define BOOST_MATH_QUADRATURE_TANH_SINH_HPP
0031
0032 #include <cmath>
0033 #include <limits>
0034 #include <memory>
0035 #include <boost/math/quadrature/detail/tanh_sinh_detail.hpp>
0036
0037 namespace boost{ namespace math{ namespace quadrature {
0038
0039 template<class Real, class Policy = policies::policy<> >
0040 class tanh_sinh
0041 {
0042 public:
0043     tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
0044     : m_imp(std::make_shared<detail::tanh_sinh_detail<Real, Policy>>(max_refinements, min_complement)) {}
0045
0046     template<class F>
0047     auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
0048     template<class F>
0049     auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));
0050
0051     template<class F>
0052     auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
0053     template<class F>
0054     auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));
0055
0056 private:
0057     std::shared_ptr<detail::tanh_sinh_detail<Real, Policy>> m_imp;
0058 };
0059
0060 template<class Real, class Policy>
0061 template<class F>
0062 auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
0063 {
0064     BOOST_MATH_STD_USING
0065     using boost::math::constants::half;
0066     using boost::math::quadrature::detail::tanh_sinh_detail;
0067
0068     static const char* function = "tanh_sinh<%1%>::integrate";
0069
0070     typedef decltype(std::declval<F>()(std::declval<Real>())) result_type;
0071     static_assert(!std::is_integral<result_type>::value,
0072                   "The return type cannot be integral, it must be either a real or complex floating point type.");
0073     if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
0074     {
0075
0076        // Infinite limits:
0077        if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
0078        {
0079           auto u = [&](const Real& t, const Real& tc)->result_type
0080           {
0081              Real t_sq = t*t;
0082              Real inv;
0083              if (t > 0.5f)
0084                 inv = 1 / ((2 - tc) * tc);
0085              else if(t < -0.5)
0086                 inv = 1 / ((2 + tc) * -tc);
0087              else
0088                 inv = 1 / (1 - t_sq);
0089              return f(t*inv)*(1 + t_sq)*inv*inv;
0090           };
0091           Real limit = sqrt(tools::min_value<Real>()) * 4;
0092           return m_imp->integrate(u, error, L1, function, limit, limit, tolerance, levels);
0093        }
0094
0095        // Right limit is infinite:
0096        if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
0097        {
0098           auto u = [&](const Real& t, const Real& tc)->result_type
0099           {
0100              Real z, arg;
0101              if (t > -0.5f)
0102                 z = 1 / (t + 1);
0103              else
0104                 z = -1 / tc;
0105              if (t < 0.5)
0106                 arg = 2 * z + a - 1;
0107              else
0108                 arg = a + tc / (2 - tc);
0109              return f(arg)*z*z;
0110           };
0111           Real left_limit = sqrt(tools::min_value<Real>()) * 4;
0112           result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
0113           if (L1)
0114           {
0115              *L1 *= 2;
0116           }
0117           if (error)
0118           {
0119              *error *= 2;
0120           }
0121
0122           return Q;
0123        }
0124
0125        if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
0126        {
0127           auto v = [&](const Real& t, const Real& tc)->result_type
0128           {
0129              Real z;
0130              if (t > -0.5)
0131                 z = 1 / (t + 1);
0132              else
0133                 z = -1 / tc;
0134              Real arg;
0135              if (t < 0.5)
0136                 arg = 2 * z - 1;
0137              else
0138                 arg = tc / (2 - tc);
0139              return f(b - arg) * z * z;
0140           };
0141
0142           Real left_limit = sqrt(tools::min_value<Real>()) * 4;
0143           result_type Q = Real(2) * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
0144           if (L1)
0145           {
0146              *L1 *= 2;
0147           }
0148           if (error)
0149           {
0150              *error *= 2;
0151           }
0152           return Q;
0153        }
0154
0155        if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
0156        {
0157           if (a == b)
0158           {
0159              return result_type(0);
0160           }
0161           if (b < a)
0162           {
0163              return -this->integrate(f, b, a, tolerance, error, L1, levels);
0164           }
0165           Real avg = (a + b)*half<Real>();
0166           Real diff = (b - a)*half<Real>();
0167           Real avg_over_diff_m1 = a / diff;
0168           Real avg_over_diff_p1 = b / diff;
0169           bool have_small_left = fabs(a) < 0.5f;
0170           bool have_small_right = fabs(b) < 0.5f;
0171           Real left_min_complement = float_next(avg_over_diff_m1) - avg_over_diff_m1;
0172           Real min_complement_limit = (std::max)(tools::min_value<Real>(), float_next(Real(tools::min_value<Real>() / diff)));
0173           if (left_min_complement < min_complement_limit)
0174              left_min_complement = min_complement_limit;
0175           Real right_min_complement = avg_over_diff_p1 - float_prior(avg_over_diff_p1);
0176           if (right_min_complement < min_complement_limit)
0177              right_min_complement = min_complement_limit;
0178           //
0179           // These asserts will fail only if rounding errors on
0180           // type Real have accumulated so much error that it's
0181           // broken our internal logic.  Should that prove to be
0182           // a persistent issue, we might need to add a bit of fudge
0183           // factor to move left_min_complement and right_min_complement
0184           // further from the end points of the range.
0185           //
0186           BOOST_MATH_ASSERT((left_min_complement * diff + a) > a);
0187           BOOST_MATH_ASSERT((b - right_min_complement * diff) < b);
0188           auto u = [&](Real z, Real zc)->result_type
0189           {
0190              Real position;
0191              if (z < -0.5)
0192              {
0193                 if(have_small_left)
0194                   return f(diff * (avg_over_diff_m1 - zc));
0195                 position = a - diff * zc;
0196              }
0197              else if (z > 0.5)
0198              {
0199                 if(have_small_right)
0200                   return f(diff * (avg_over_diff_p1 - zc));
0201                 position = b - diff * zc;
0202              }
0203              else
0204                 position = avg + diff*z;
0205              BOOST_MATH_ASSERT(position != a);
0206              BOOST_MATH_ASSERT(position != b);
0207              return f(position);
0208           };
0209           result_type Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
0210
0211           if (L1)
0212           {
0213              *L1 *= diff;
0214           }
0215           if (error)
0216           {
0217              *error *= diff;
0218           }
0219           return Q;
0220        }
0221     }
0222     return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
0223 }
0224
0225 template<class Real, class Policy>
0226 template<class F>
0227 auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()))
0228 {
0229    BOOST_MATH_STD_USING
0230       using boost::math::constants::half;
0231    using boost::math::quadrature::detail::tanh_sinh_detail;
0232
0233    static const char* function = "tanh_sinh<%1%>::integrate";
0234
0235    if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
0236    {
0237       if (b <= a)
0238       {
0239          return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
0240       }
0241       auto u = [&](Real z, Real zc)->Real
0242       {
0243          if (z < 0)
0244             return f((a - b) * zc / 2 + a, (b - a) * zc / 2);
0245          else
0246             return f((a - b) * zc / 2 + b, (b - a) * zc / 2);
0247       };
0248       Real diff = (b - a)*half<Real>();
0249       Real left_min_complement = tools::min_value<Real>() * 4;
0250       Real right_min_complement = tools::min_value<Real>() * 4;
0251       Real Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
0252
0253       if (L1)
0254       {
0255          *L1 *= diff;
0256       }
0257       if (error)
0258       {
0259          *error *= diff;
0260       }
0261       return Q;
0262    }
0263    return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
0264 }
0265
0266 template<class Real, class Policy>
0267 template<class F>
0268 auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
0269 {
0270    using boost::math::quadrature::detail::tanh_sinh_detail;
0271    static const char* function = "tanh_sinh<%1%>::integrate";
0272    Real min_complement = tools::epsilon<Real>();
0273    return m_imp->integrate([&](const Real& arg, const Real&) { return f(arg); }, error, L1, function, min_complement, min_complement, tolerance, levels);
0274 }
0275
0276 template<class Real, class Policy>
0277 template<class F>
0278 auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()))
0279 {
0280    using boost::math::quadrature::detail::tanh_sinh_detail;
0281    static const char* function = "tanh_sinh<%1%>::integrate";
0282    Real min_complement = tools::min_value<Real>() * 4;
0283    return m_imp->integrate(f, error, L1, function, min_complement, min_complement, tolerance, levels);
0284 }
0285
0286 }
0287 }
0288 }
0289 #endif