boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-26 06:42:12 +00:00
Code Issues Projects Releases Wiki Activity

Quadrature: add gauss and gauss-kronrod quadrature.

Browse Source
This commit is contained in:
jzmaddock
2017-08-31 19:42:26 +01:00
parent e7e915816f
commit 7d2002db80
5 changed files with 1472 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

183
include/boost/math/quadrature/gauss.hpp Normal file
View File

@@ -0,0 +1,183 @@
// Copyright John Maddock 2015.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_QUADRATURE_GAUSS_HPP
#define BOOST_MATH_QUADRATURE_GAUSS_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include <vector>
#include <boost/math/special_functions/legendre.hpp>
namespace boost { namespace math{ namespace quadrature{ namespace detail{
template <class Real, unsigned N>
class gauss_detail
{
static std::vector<Real> calculate_weights()
{
std::vector<Real> result(abscissa().size(), 0);
for (unsigned i = 0; i < abscissa().size(); ++i)
{
Real x = abscissa()[i];
Real p = boost::math::legendre_p_prime(N, x);
result[i] = 2 / ((1 - x * x) * p * p);
}
return result;
}
public:
static const std::vector<Real>& abscissa()
{
static std::vector<Real> data = boost::math::legendre_p_zeros<Real>(N);
return data;
}
static const std::vector<Real>& weights()
{
static std::vector<Real> data = calculate_weights();
return data;
}
};
template <>
class gauss_detail<double, 7>
{
public:
static constexpr std::array<double, 4> const & abscissa()
{
static constexpr std::array<double, 4> data = {
0.000000000000000000000000000000000e+00,
4.058451513773971669066064120769615e-01,
7.415311855993944398638647732807884e-01,
9.491079123427585245261896840478513e-01,
};
return data;
}
static constexpr std::array<double, 4> const & weights()
{
static constexpr std::array<double, 4> data = {
4.179591836734693877551020408163265e-01,
3.818300505051189449503697754889751e-01,
2.797053914892766679014677714237796e-01,
1.294849661688696932706114326790820e-01,
};
return data;
}
};
}
template <class Real, unsigned N, class Policy = boost::math::policies::policy<> >
class gauss : public detail::gauss_detail<Real, N>
{
public:
typedef Real value_type;
template <class F>
static value_type integrate(F f, Real* pL1 = nullptr)
{
using std::fabs;
unsigned non_zero_start = 1;
value_type result = 0;
if (N & 1)
result = f(value_type(0)) * weights()[0];
else
non_zero_start = 0;
value_type L1 = fabs(result);
for (unsigned i = non_zero_start; i < abscissa().size(); ++i)
{
value_type fp = f(abscissa()[i]);
value_type fm = f(-abscissa()[i]);
result += (fp + fm) * weights()[i];
L1 += (fabs(fp) + fabs(fm)) * weights()[i];
}
if (pL1)
*pL1 = L1;
return result;
}
template <class F>
static value_type integrate(F f, Real a, Real b, Real* pL1 = nullptr)
{
static const char* function = "boost::math::quadrature::gauss<%1%>::integrate(f, %1%, %1%)";
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
// Infinite limits:
if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t)->Real
{
Real t_sq = t*t;
Real inv = 1 / (1 - t_sq);
return f(t*inv)*(1 + t_sq)*inv*inv;
};
return integrate(u, pL1);
}
// Right limit is infinite:
if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t)->Real
{
Real z = 1 / (t + 1);
Real arg = 2 * z + a - 1;
return f(arg)*z*z;
};
Real Q = 2 * integrate(u, pL1);
if (pL1)
{
*pL1 *= 2;
}
return Q;
}
if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
{
auto v = [&](const Real& t)->Real
{
Real z = 1 / (t + 1);
Real arg = 2 * z - 1;
return f(b - arg) * z * z;
};
Real Q = 2 * integrate(v, pL1);
if (pL1)
{
*pL1 *= 2;
}
return Q;
}
if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
{
if (b <= a)
{
return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
}
Real avg = (a + b)*half<Real>();
Real scale = (b - a)*half<Real>();
auto u = [&](Real z)->Real
{
return f(avg + scale*z);
};
Real Q = scale*integrate(u, pL1);
if (pL1)
{
*pL1 *= scale;
}
return Q;
}
}
return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
}
};
} // namespace quadrature
} // namespace math
} // namespace boost
#endif // BOOST_MATH_QUADRATURE_GAUSS_HPP

281
include/boost/math/quadrature/gauss_kronrod.hpp Normal file
View File

@@ -0,0 +1,281 @@
// Copyright John Maddock 2017.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP
#define BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP
#ifdef _MSC_VER
#pragma once
#pragma warning(push)
#pragma warning(disable: 4127)
#endif
#include <vector>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/legendre_stieltjes.hpp>
#include <boost/math/quadrature/gauss.hpp>
namespace boost { namespace math{ namespace quadrature{ namespace detail{
template <class Real, unsigned N>
class gauss_kronrod_detail
{
static legendre_stieltjes<Real> const& get_legendre_stieltjes()
{
static const legendre_stieltjes<Real> data((N - 1) / 2 + 1);
return data;
}
static std::vector<Real> calculate_abscissa()
{
static std::vector<Real> result = boost::math::legendre_p_zeros<Real>((N - 1) / 2);
const legendre_stieltjes<Real> E = get_legendre_stieltjes();
std::vector<Real> ls_zeros = E.zeros();
result.insert(result.end(), ls_zeros.begin(), ls_zeros.end());
std::sort(result.begin(), result.end());
return result;
}
static std::vector<Real> calculate_weights()
{
std::vector<Real> result(abscissa().size(), 0);
unsigned gauss_order = (N - 1) / 2;
unsigned gauss_start = gauss_order & 1 ? 0 : 1;
const legendre_stieltjes<Real>& E = get_legendre_stieltjes();
for (unsigned i = gauss_start; i < abscissa().size(); i += 2)
{
Real x = abscissa()[i];
Real p = boost::math::legendre_p_prime(gauss_order, x);
Real gauss_weight = 2 / ((1 - x * x) * p * p);
result[i] = gauss_weight + static_cast<Real>(2) / (static_cast<Real>(gauss_order + 1) * legendre_p_prime(gauss_order, x) * E(x));
}
for (unsigned i = gauss_start ? 0 : 1; i < abscissa().size(); i += 2)
{
Real x = abscissa()[i];
result[i] = static_cast<Real>(2) / (static_cast<Real>(gauss_order + 1) * legendre_p(gauss_order, x) * E.prime(x));
}
return result;
}
public:
static const std::vector<Real>& abscissa()
{
static std::vector<Real> data = calculate_abscissa();
return data;
}
static const std::vector<Real>& weights()
{
static std::vector<Real> data = calculate_weights();
return data;
}
};
template <>
class gauss_kronrod_detail<double, 15>
{
public:
static constexpr std::array<double, 8> const & abscissa()
{
static constexpr std::array<double, 8> data = {
0.000000000000000000000000000000000e+00,
2.077849550078984676006894037732449e-01,
4.058451513773971669066064120769615e-01,
5.860872354676911302941448382587296e-01,
7.415311855993944398638647732807884e-01,
8.648644233597690727897127886409262e-01,
9.491079123427585245261896840478513e-01,
9.914553711208126392068546975263285e-01,
};
return data;
}
static constexpr std::array<double, 8> const & weights()
{
static constexpr std::array<double, 8> data = {
2.094821410847278280129991748917143e-01,
2.044329400752988924141619992346491e-01,
1.903505780647854099132564024210137e-01,
1.690047266392679028265834265985503e-01,
1.406532597155259187451895905102379e-01,
1.047900103222501838398763225415180e-01,
6.309209262997855329070066318920429e-02,
2.293532201052922496373200805896959e-02,
};
return data;
}
};
}
template <class Real, unsigned N, class Policy = boost::math::policies::policy<> >
class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N>
{
public:
typedef Real value_type;
private:
template <class F>
static value_type integrate_non_adaptive_m1_1(F f, Real* error = nullptr, Real* pL1 = nullptr)
{
using std::fabs;
unsigned gauss_start = 2;
unsigned kronrod_start = 1;
unsigned gauss_order = (N - 1) / 2;
value_type kronrod_result = 0;
value_type gauss_result = 0;
value_type fp, fm;
if (gauss_order & 1)
{
fp = f(value_type(0));
kronrod_result = fp * weights()[0];
gauss_result += fp * gauss<Real, (N - 1) / 2>::weights()[0];
}
else
{
fp = f(value_type(0));
kronrod_result = fp * weights()[0];
gauss_start = 1;
kronrod_start = 2;
}
value_type L1 = fabs(kronrod_result);
for (unsigned i = gauss_start; i < abscissa().size(); i += 2)
{
fp = f(abscissa()[i]);
fm = f(-abscissa()[i]);
kronrod_result += (fp + fm) * weights()[i];
L1 += (fabs(fp) + fabs(fm)) * weights()[i];
gauss_result += (fp + fm) * gauss<Real, (N - 1) / 2>::weights()[i / 2];
}
for (unsigned i = kronrod_start; i < abscissa().size(); i += 2)
{
fp = f(abscissa()[i]);
fm = f(-abscissa()[i]);
kronrod_result += (fp + fm) * weights()[i];
L1 += (fabs(fp) + fabs(fm)) * weights()[i];
}
if (pL1)
*pL1 = L1;
if (error)
*error = fabs(kronrod_result - gauss_result);
return kronrod_result;
}
template <class F>
struct recursive_info
{
F f;
Real tol;
};
template <class F>
static value_type recursive_adaptive_integrate(const recursive_info<F>* info, Real a, Real b, unsigned max_levels, Real abs_tol, Real* error, Real* L1)
{
Real error_local;
Real mean = (b + a) / 2;
Real scale = (b - a) / 2;
auto ff = [&](const Real& x)
{
return info->f(scale * x + mean);
};
Real estimate = scale * integrate_non_adaptive_m1_1(ff, &error_local, L1);
Real abs_tol1 = fabs(estimate * info->tol);
if (abs_tol == 0)
abs_tol = abs_tol1;
if (max_levels && (abs_tol1 < error_local) && (abs_tol < error_local))
{
Real mid = (a + b) / 2;
Real L1_local;
estimate = recursive_adaptive_integrate(info, a, mid, max_levels - 1, abs_tol / 2, error, L1);
estimate += recursive_adaptive_integrate(info, mid, b, max_levels - 1, abs_tol / 2, &error_local, &L1_local);
if (error)
*error += error_local;
if (L1)
*L1 += L1_local;
return estimate;
}
if(L1)
*L1 *= scale;
if (error)
*error = error_local;
return estimate;
}
public:
template <class F>
static value_type integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr)
{
static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)";
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
// Infinite limits:
if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t)->Real
{
Real t_sq = t*t;
Real inv = 1 / (1 - t_sq);
return f(t*inv)*(1 + t_sq)*inv*inv;
};
recursive_info<decltype(u)> info = { u, tol };
return recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
}
// Right limit is infinite:
if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t)->Real
{
Real z = 1 / (t + 1);
Real arg = 2 * z + a - 1;
return f(arg)*z*z;
};
recursive_info<decltype(u)> info = { u, tol };
Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
if (pL1)
{
*pL1 *= 2;
}
return Q;
}
if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
{
auto v = [&](const Real& t)->Real
{
Real z = 1 / (t + 1);
Real arg = 2 * z - 1;
return f(b - arg) * z * z;
};
recursive_info<decltype(v)> info = { v, tol };
Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
if (pL1)
{
*pL1 *= 2;
}
return Q;
}
if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
{
if (b <= a)
{
return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
}
recursive_info<F> info = { f, tol };
return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1);
}
}
return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
}
};
} // namespace quadrature
} // namespace math
} // namespace boost
#ifdef _MSC_VER
#pragma warning(pop)
#endif
#endif // BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP

259
test/adaptive_gauss_kronrod_quadrature_test.cpp Normal file
View File

@@ -0,0 +1,259 @@
// Copyright Nick Thompson, 2017
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#define BOOST_TEST_MODULE tanh_sinh_quadrature_test
#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>
#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#ifdef _MSC_VER
#pragma warning(disable:4127) // Conditional expression is constant
#endif
using std::expm1;
using std::atan;
using std::tan;
using std::log;
using std::log1p;
using std::asinh;
using std::atanh;
using std::sqrt;
using std::isnormal;
using std::abs;
using std::sinh;
using std::tanh;
using std::cosh;
using std::pow;
using std::exp;
using std::sin;
using std::cos;
using std::string;
using boost::math::quadrature::gauss_kronrod;
using boost::math::constants::pi;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::math::constants::two_pi;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::catalan;
using boost::math::constants::ln_two;
using boost::math::constants::root_two;
using boost::math::constants::root_two_pi;
using boost::math::constants::root_pi;
using boost::multiprecision::cpp_bin_float_quad;
template <class Real>
Real get_termination_condition()
{
return boost::math::tools::epsilon<Real>() * 1000;
}
template<class Real, unsigned Points>
void test_linear()
{
std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real error;
auto f = [](const Real& x)
{
return 5*x + 7;
};
Real L1;
Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 15, get_termination_condition<Real>(), &error, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
}
template<class Real, unsigned Points>
void test_quadratic()
{
std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real error;
auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
Real L1;
Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
}
// Examples taken from
//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
template<class Real, unsigned Points>
void test_ca()
{
std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real L1;
Real error;
auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 15, get_termination_condition<Real>(), &error, &L1);
Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
Q = gauss_kronrod<Real, Points>::integrate(f5, 0, 1, 25);
Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 100 * tol);
}
template<class Real, unsigned Points>
void test_three_quadrature_schemes_examples()
{
std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real Q;
Real Q_expected;
// Example 1:
auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1);
Q_expected = half<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 2:
auto f2 = [](const Real& t) { return t*t*atan(t); };
Q = gauss_kronrod<Real, Points>::integrate(f2, 0, 1);
Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
// Example 3:
auto f3 = [](const Real& t) { return exp(t)*cos(t); };
Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>());
Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 4:
auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
Q = gauss_kronrod<Real, Points>::integrate(f4, 0, 1);
Q_expected = 5*pi<Real>()*pi<Real>()/96;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_integration_over_real_line()
{
std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real Q;
Real Q_expected;
Real L1;
Real error;
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f4 = [](const Real& t) { return 1/cosh(t);};
Q = gauss_kronrod<Real, Points>::integrate(f4, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_right_limit_infinite()
{
std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real Q;
Real Q_expected;
Real L1;
Real error;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
auto f4 = [](const Real& t) { return 1/(1+t*t); };
Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
Q_expected = pi<Real>()/4;
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
template<class Real, unsigned Points>
void test_left_limit_infinite()
{
std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real Q;
Real Q_expected;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), 0);
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
{
test_linear<double, 15>();
test_quadratic<double, 15>();
test_ca<double, 15>();
test_three_quadrature_schemes_examples<double, 15>();
test_integration_over_real_line<double, 15>();
test_right_limit_infinite<double, 15>();
test_left_limit_infinite<double, 15>();
test_linear<cpp_bin_float_quad, 21>();
test_quadratic<cpp_bin_float_quad, 21>();
test_ca<cpp_bin_float_quad, 21>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
test_integration_over_real_line<cpp_bin_float_quad, 21>();
test_right_limit_infinite<cpp_bin_float_quad, 21>();
test_left_limit_infinite<cpp_bin_float_quad, 21>();
test_linear<cpp_bin_float_quad, 31>();
test_quadratic<cpp_bin_float_quad, 31>();
test_ca<cpp_bin_float_quad, 31>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
test_integration_over_real_line<cpp_bin_float_quad, 31>();
test_right_limit_infinite<cpp_bin_float_quad, 31>();
test_left_limit_infinite<cpp_bin_float_quad, 31>();
test_linear<cpp_bin_float_quad, 41>();
test_quadratic<cpp_bin_float_quad, 41>();
test_ca<cpp_bin_float_quad, 41>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
test_integration_over_real_line<cpp_bin_float_quad, 41>();
test_right_limit_infinite<cpp_bin_float_quad, 41>();
test_left_limit_infinite<cpp_bin_float_quad, 41>();
}
#else
int main() { return 0; }
#endif

377
test/gauss_kronrod_quadrature_test.cpp Normal file
View File

@@ -0,0 +1,377 @@
// Copyright Nick Thompson, 2017
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#define BOOST_TEST_MODULE tanh_sinh_quadrature_test
#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>
#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#ifdef _MSC_VER
#pragma warning(disable:4127) // Conditional expression is constant
#endif
using std::expm1;
using std::atan;
using std::tan;
using std::log;
using std::log1p;
using std::asinh;
using std::atanh;
using std::sqrt;
using std::isnormal;
using std::abs;
using std::sinh;
using std::tanh;
using std::cosh;
using std::pow;
using std::exp;
using std::sin;
using std::cos;
using std::string;
using boost::math::quadrature::gauss_kronrod;
using boost::math::constants::pi;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::math::constants::two_pi;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::catalan;
using boost::math::constants::ln_two;
using boost::math::constants::root_two;
using boost::math::constants::root_two_pi;
using boost::math::constants::root_pi;
using boost::multiprecision::cpp_bin_float_quad;
//
// Error rates depend only on the number of points in the approximation, not the type being tested,
// define all our expected errors here:
//
enum
{
test_ca_error_id,
test_ca_error_id_2,
test_three_quad_error_id,
test_three_quad_error_id_2,
test_integration_over_real_line_error_id,
test_right_limit_infinite_error_id,
test_left_limit_infinite_error_id
};
template <unsigned Points>
double expected_error(unsigned)
{
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<15>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-7;
case test_ca_error_id_2:
return 2e-5;
case test_three_quad_error_id:
return 1e-8;
case test_three_quad_error_id_2:
return 3.5e-3;
case test_integration_over_real_line_error_id:
return 6e-3;
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 1e-5;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<21>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-12;
case test_ca_error_id_2:
return 3e-6;
case test_three_quad_error_id:
return 2e-13;
case test_three_quad_error_id_2:
return 2e-3;
case test_integration_over_real_line_error_id:
return 6e-3; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 5e-8;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<31>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 6e-20;
case test_ca_error_id_2:
return 3e-7;
case test_three_quad_error_id:
return 1e-19;
case test_three_quad_error_id_2:
return 6e-4;
case test_integration_over_real_line_error_id:
return 6e-3; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 5e-11;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<41>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-26;
case test_ca_error_id_2:
return 1e-7;
case test_three_quad_error_id:
return 3e-27;
case test_three_quad_error_id_2:
return 3e-4;
case test_integration_over_real_line_error_id:
return 5e-5; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 1e-15;
}
return 0; // placeholder, all tests will fail
}
template<class Real, unsigned Points>
void test_linear()
{
std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real error;
auto f = [](const Real& x)
{
return 5*x + 7;
};
Real L1;
Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 0, 0, &error, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
}
template<class Real, unsigned Points>
void test_quadratic()
{
std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
Real error;
auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
Real L1;
Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 0, 0, &error, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
}
// Examples taken from
//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
template<class Real, unsigned Points>
void test_ca()
{
std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_ca_error_id);
Real L1;
Real error;
auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 0, 0, &error, &L1);
Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 0, 0, &error, &L1);
Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
tol = expected_error<Points>(test_ca_error_id_2);
auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
Q = gauss_kronrod<Real, Points>::integrate(f5, 0, 1, 0);
Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_three_quadrature_schemes_examples()
{
std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_three_quad_error_id);
Real Q;
Real Q_expected;
// Example 1:
auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1, 0);
Q_expected = half<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 2:
auto f2 = [](const Real& t) { return t*t*atan(t); };
Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 0);
Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
// Example 3:
auto f3 = [](const Real& t) { return exp(t)*cos(t); };
Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>(), 0);
Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 4:
auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
Q = gauss_kronrod<Real, Points>::integrate(f4, 0 , 1, 0);
Q_expected = 5*pi<Real>()*pi<Real>()/96;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
tol = expected_error<Points>(test_three_quad_error_id_2);
// Example 5:
auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
Q = gauss_kronrod<Real, Points>::integrate(f5, 0 , 1, 0);
Q_expected = -4/ (Real) 9;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 6:
auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
Q = gauss_kronrod<Real, Points>::integrate(f6, 0 , 1, 0);
Q_expected = pi<Real>()/4;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_integration_over_real_line()
{
std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
Real Q;
Real Q_expected;
Real L1;
Real error;
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f4 = [](const Real& t) { return 1/cosh(t);};
Q = gauss_kronrod<Real, Points>::integrate(f4, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_right_limit_infinite()
{
std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
Real Q;
Real Q_expected;
Real L1;
Real error;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
auto f4 = [](const Real& t) { return 1/(1+t*t); };
Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
Q_expected = pi<Real>()/4;
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
template<class Real, unsigned Points>
void test_left_limit_infinite()
{
std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
Real Q;
Real Q_expected;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0), 0);
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
{
test_linear<double, 15>();
test_quadratic<double, 15>();
test_ca<double, 15>();
test_three_quadrature_schemes_examples<double, 15>();
test_integration_over_real_line<double, 15>();
test_right_limit_infinite<double, 15>();
test_left_limit_infinite<double, 15>();
test_linear<cpp_bin_float_quad, 21>();
test_quadratic<cpp_bin_float_quad, 21>();
test_ca<cpp_bin_float_quad, 21>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
test_integration_over_real_line<cpp_bin_float_quad, 21>();
test_right_limit_infinite<cpp_bin_float_quad, 21>();
test_left_limit_infinite<cpp_bin_float_quad, 21>();
test_linear<cpp_bin_float_quad, 31>();
test_quadratic<cpp_bin_float_quad, 31>();
test_ca<cpp_bin_float_quad, 31>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
test_integration_over_real_line<cpp_bin_float_quad, 31>();
test_right_limit_infinite<cpp_bin_float_quad, 31>();
test_left_limit_infinite<cpp_bin_float_quad, 31>();
test_linear<cpp_bin_float_quad, 41>();
test_quadratic<cpp_bin_float_quad, 41>();
test_ca<cpp_bin_float_quad, 41>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
test_integration_over_real_line<cpp_bin_float_quad, 41>();
test_right_limit_infinite<cpp_bin_float_quad, 41>();
test_left_limit_infinite<cpp_bin_float_quad, 41>();
}
#else
int main() { return 0; }
#endif

372
test/gauss_quadrature_test.cpp Normal file
View File

@@ -0,0 +1,372 @@
// Copyright Nick Thompson, 2017
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#define BOOST_TEST_MODULE tanh_sinh_quadrature_test
#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>
#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/quadrature/gauss.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#ifdef _MSC_VER
#pragma warning(disable:4127) // Conditional expression is constant
#endif
using std::expm1;
using std::atan;
using std::tan;
using std::log;
using std::log1p;
using std::asinh;
using std::atanh;
using std::sqrt;
using std::isnormal;
using std::abs;
using std::sinh;
using std::tanh;
using std::cosh;
using std::pow;
using std::exp;
using std::sin;
using std::cos;
using std::string;
using boost::math::quadrature::gauss;
using boost::math::constants::pi;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::math::constants::two_pi;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::catalan;
using boost::math::constants::ln_two;
using boost::math::constants::root_two;
using boost::math::constants::root_two_pi;
using boost::math::constants::root_pi;
using boost::multiprecision::cpp_bin_float_quad;
//
// Error rates depend only on the number of points in the approximation, not the type being tested,
// define all our expected errors here:
//
enum
{
test_ca_error_id,
test_ca_error_id_2,
test_three_quad_error_id,
test_three_quad_error_id_2,
test_integration_over_real_line_error_id,
test_right_limit_infinite_error_id,
test_left_limit_infinite_error_id
};
template <unsigned Points>
double expected_error(unsigned)
{
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<7>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-7;
case test_ca_error_id_2:
return 2e-5;
case test_three_quad_error_id:
return 1e-8;
case test_three_quad_error_id_2:
return 3.5e-3;
case test_integration_over_real_line_error_id:
return 6e-3;
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 1e-5;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<10>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-12;
case test_ca_error_id_2:
return 3e-6;
case test_three_quad_error_id:
return 2e-13;
case test_three_quad_error_id_2:
return 2e-3;
case test_integration_over_real_line_error_id:
return 6e-3; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 5e-8;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<15>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 6e-20;
case test_ca_error_id_2:
return 3e-7;
case test_three_quad_error_id:
return 1e-19;
case test_three_quad_error_id_2:
return 6e-4;
case test_integration_over_real_line_error_id:
return 6e-3; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 5e-11;
}
return 0; // placeholder, all tests will fail
}
template <>
double expected_error<20>(unsigned id)
{
switch (id)
{
case test_ca_error_id:
return 1e-26;
case test_ca_error_id_2:
return 1e-7;
case test_three_quad_error_id:
return 3e-27;
case test_three_quad_error_id_2:
return 3e-4;
case test_integration_over_real_line_error_id:
return 5e-5; // doesn't get any better with more points!
case test_right_limit_infinite_error_id:
case test_left_limit_infinite_error_id:
return 1e-15;
}
return 0; // placeholder, all tests will fail
}
template<class Real, unsigned Points>
void test_linear()
{
std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
auto f = [](const Real& x)
{
return 5*x + 7;
};
Real L1;
Real Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 1, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
}
template<class Real, unsigned Points>
void test_quadratic()
{
std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = boost::math::tools::epsilon<Real>() * 10;
auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
Real L1;
Real Q = gauss<Real, Points>::integrate(f, 0, 1, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
}
// Examples taken from
//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
template<class Real, unsigned Points>
void test_ca()
{
std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_ca_error_id);
Real L1;
auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
Q = gauss<Real, Points>::integrate(f2, 0 , 1, &L1);
Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
tol = expected_error<Points>(test_ca_error_id_2);
auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
Q = gauss<Real, Points>::integrate(f5, 0 , 1);
Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_three_quadrature_schemes_examples()
{
std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_three_quad_error_id);
Real Q;
Real Q_expected;
// Example 1:
auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
Q = gauss<Real, Points>::integrate(f1, 0 , 1);
Q_expected = half<Real>()*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 2:
auto f2 = [](const Real& t) { return t*t*atan(t); };
Q = gauss<Real, Points>::integrate(f2, 0 , 1);
Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
// Example 3:
auto f3 = [](const Real& t) { return exp(t)*cos(t); };
Q = gauss<Real, Points>::integrate(f3, 0, half_pi<Real>());
Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 4:
auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
Q = gauss<Real, Points>::integrate(f4, 0 , 1);
Q_expected = 5*pi<Real>()*pi<Real>()/96;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
tol = expected_error<Points>(test_three_quad_error_id_2);
// Example 5:
auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
Q = gauss<Real, Points>::integrate(f5, 0 , 1);
Q_expected = -4/ (Real) 9;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
// Example 6:
auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
Q = gauss<Real, Points>::integrate(f6, 0 , 1);
Q_expected = pi<Real>()/4;
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_integration_over_real_line()
{
std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
Real Q;
Real Q_expected;
Real L1;
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
auto f4 = [](const Real& t) { return 1/cosh(t);};
Q = gauss<Real, Points>::integrate(f4, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
Q_expected = pi<Real>();
BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}
template<class Real, unsigned Points>
void test_right_limit_infinite()
{
std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
Real Q;
Real Q_expected;
Real L1;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), &L1);
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
auto f4 = [](const Real& t) { return 1/(1+t*t); };
Q = gauss<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), &L1);
Q_expected = pi<Real>()/4;
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
template<class Real, unsigned Points>
void test_left_limit_infinite()
{
std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
Real Q;
Real Q_expected;
// Example 11:
auto f1 = [](const Real& t) { return 1/(1+t*t);};
Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0));
Q_expected = half_pi<Real>();
BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}
BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
{
test_linear<double, 7>();
test_quadratic<double, 7>();
test_ca<double, 7>();
test_three_quadrature_schemes_examples<double, 7>();
test_integration_over_real_line<double, 7>();
test_right_limit_infinite<double, 7>();
test_left_limit_infinite<double, 7>();
test_linear<cpp_bin_float_quad, 10>();
test_quadratic<cpp_bin_float_quad, 10>();
test_ca<cpp_bin_float_quad, 10>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 10>();
test_integration_over_real_line<cpp_bin_float_quad, 10>();
test_right_limit_infinite<cpp_bin_float_quad, 10>();
test_left_limit_infinite<cpp_bin_float_quad, 10>();
test_linear<cpp_bin_float_quad, 15>();
test_quadratic<cpp_bin_float_quad, 15>();
test_ca<cpp_bin_float_quad, 15>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 15>();
test_integration_over_real_line<cpp_bin_float_quad, 15>();
test_right_limit_infinite<cpp_bin_float_quad, 15>();
test_left_limit_infinite<cpp_bin_float_quad, 15>();
test_linear<cpp_bin_float_quad, 20>();
test_quadratic<cpp_bin_float_quad, 20>();
test_ca<cpp_bin_float_quad, 20>();
test_three_quadrature_schemes_examples<cpp_bin_float_quad, 20>();
test_integration_over_real_line<cpp_bin_float_quad, 20>();
test_right_limit_infinite<cpp_bin_float_quad, 20>();
test_left_limit_infinite<cpp_bin_float_quad, 20>();
}
#else
int main() { return 0; }
#endif