]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Copyright Nick Thompson, 2017 |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. | |
4 | // (See accompanying file LICENSE_1_0.txt | |
5 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | /* | |
8 | * This class performs exp-sinh quadrature on half infinite intervals. | |
9 | * | |
10 | * References: | |
11 | * | |
12 | * 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655. | |
13 | */ | |
14 | ||
15 | #ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP | |
16 | #define BOOST_MATH_QUADRATURE_EXP_SINH_HPP | |
17 | ||
18 | #include <cmath> | |
19 | #include <limits> | |
20 | #include <memory> | |
1e59de90 | 21 | #include <string> |
b32b8144 FG |
22 | #include <boost/math/quadrature/detail/exp_sinh_detail.hpp> |
23 | ||
24 | namespace boost{ namespace math{ namespace quadrature { | |
25 | ||
26 | template<class Real, class Policy = policies::policy<> > | |
27 | class exp_sinh | |
28 | { | |
29 | public: | |
30 | exp_sinh(size_t max_refinements = 9) | |
31 | : m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {} | |
32 | ||
33 | template<class F> | |
92f5a8d4 | 34 | auto integrate(const F& f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const; |
b32b8144 | 35 | template<class F> |
92f5a8d4 | 36 | auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const; |
b32b8144 FG |
37 | |
38 | private: | |
39 | std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp; | |
40 | }; | |
41 | ||
42 | template<class Real, class Policy> | |
43 | template<class F> | |
92f5a8d4 | 44 | auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) const |
b32b8144 | 45 | { |
92f5a8d4 | 46 | typedef decltype(f(a)) K; |
1e59de90 TL |
47 | static_assert(!std::is_integral<K>::value, |
48 | "The return type cannot be integral, it must be either a real or complex floating point type."); | |
b32b8144 FG |
49 | using std::abs; |
50 | using boost::math::constants::half; | |
51 | using boost::math::quadrature::detail::exp_sinh_detail; | |
52 | ||
53 | static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate"; | |
54 | ||
55 | // Neither limit may be a NaN: | |
56 | if((boost::math::isnan)(a) || (boost::math::isnan)(b)) | |
92f5a8d4 TL |
57 | { |
58 | return static_cast<K>(policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy())); | |
59 | } | |
b32b8144 FG |
60 | // Right limit is infinite: |
61 | if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>())) | |
62 | { | |
63 | // If a = 0, don't use an additional level of indirection: | |
1e59de90 | 64 | if (a == static_cast<Real>(0)) |
b32b8144 FG |
65 | { |
66 | return m_imp->integrate(f, error, L1, function, tolerance, levels); | |
67 | } | |
f67539c2 | 68 | const auto u = [&](Real t)->K { return f(t + a); }; |
b32b8144 FG |
69 | return m_imp->integrate(u, error, L1, function, tolerance, levels); |
70 | } | |
71 | ||
72 | if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>()) | |
73 | { | |
f67539c2 | 74 | const auto u = [&](Real t)->K { return f(b-t);}; |
b32b8144 FG |
75 | return m_imp->integrate(u, error, L1, function, tolerance, levels); |
76 | } | |
77 | ||
78 | // Infinite limits: | |
79 | if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>())) | |
80 | { | |
92f5a8d4 | 81 | return static_cast<K>(policies::raise_domain_error(function, "Use sinh_sinh quadrature for integration over the whole real line; exp_sinh is for half infinite integrals.", a, Policy())); |
b32b8144 FG |
82 | } |
83 | // If we get to here then both ends must necessarily be finite: | |
92f5a8d4 | 84 | return static_cast<K>(policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy())); |
b32b8144 FG |
85 | } |
86 | ||
87 | template<class Real, class Policy> | |
88 | template<class F> | |
92f5a8d4 | 89 | auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) const |
b32b8144 | 90 | { |
b32b8144 | 91 | static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate"; |
1e59de90 TL |
92 | using std::abs; |
93 | if (abs(tolerance) > 1) { | |
94 | std::string msg = std::string(__FILE__) + ":" + std::to_string(__LINE__) + ":" + std::string(function) + ": The tolerance provided is unusually large; did you confuse it with a domain bound?"; | |
95 | throw std::domain_error(msg); | |
96 | } | |
b32b8144 FG |
97 | return m_imp->integrate(f, error, L1, function, tolerance, levels); |
98 | } | |
99 | ||
100 | ||
101 | }}} | |
102 | #endif |