libs/math/test/adaptive_gauss_kronrod_quadrature_test.cpp

   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 #define BOOST_TEST_MODULE adaptive_gauss_kronrod_quadrature_test
   8
   9 #include <boost/config.hpp>
  10 #include <boost/detail/workaround.hpp>
  11
  12 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
  13
  14 #include <boost/math/concepts/real_concept.hpp>
  15 #include <boost/test/included/unit_test.hpp>
  16 #include <boost/test/tools/floating_point_comparison.hpp>
  17 #include <boost/math/quadrature/gauss_kronrod.hpp>
  18 #include <boost/math/special_functions/sinc.hpp>
  19 #include <boost/multiprecision/cpp_bin_float.hpp>
  20
  21 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
  22 #  define TEST1
  23 #  define TEST1A
  24 #  define TEST2
  25 #  define TEST3
  26 #endif
  27
  28 #ifdef _MSC_VER
  29 #pragma warning(disable:4127)  // Conditional expression is constant
  30 #endif
  31
  32 using std::expm1;
  33 using std::atan;
  34 using std::tan;
  35 using std::log;
  36 using std::log1p;
  37 using std::asinh;
  38 using std::atanh;
  39 using std::sqrt;
  40 using std::isnormal;
  41 using std::abs;
  42 using std::sinh;
  43 using std::tanh;
  44 using std::cosh;
  45 using std::pow;
  46 using std::exp;
  47 using std::sin;
  48 using std::cos;
  49 using std::string;
  50 using boost::math::quadrature::gauss_kronrod;
  51 using boost::math::constants::pi;
  52 using boost::math::constants::half_pi;
  53 using boost::math::constants::two_div_pi;
  54 using boost::math::constants::two_pi;
  55 using boost::math::constants::half;
  56 using boost::math::constants::third;
  57 using boost::math::constants::half;
  58 using boost::math::constants::third;
  59 using boost::math::constants::catalan;
  60 using boost::math::constants::ln_two;
  61 using boost::math::constants::root_two;
  62 using boost::math::constants::root_two_pi;
  63 using boost::math::constants::root_pi;
  64 using boost::multiprecision::cpp_bin_float_quad;
  65
  66 template <class Real>
  67 Real get_termination_condition()
  68 {
  69    return boost::math::tools::epsilon<Real>() * 1000;
  70 }
  71
  72
  73 template<class Real, unsigned Points>
  74 void test_linear()
  75 {
  76     std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  77     Real tol = boost::math::tools::epsilon<Real>() * 10;
  78     Real error;
  79     auto f = [](const Real& x)
  80     {
  81        return 5*x + 7;
  82     };
  83     Real L1;
  84     Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 15, get_termination_condition<Real>(), &error, &L1);
  85     BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
  86     BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
  87     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
  88     BOOST_CHECK_GE(fabs(error), fabs(Q - 9.5));
  89 }
  90
  91 template<class Real, unsigned Points>
  92 void test_quadratic()
  93 {
  94     std::cout << "Testing quadratic functions are integrated properly by gauss-kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  95     Real tol = boost::math::tools::epsilon<Real>() * 10;
  96     Real error;
  97
  98     auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
  99     Real L1;
 100     Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
 101     BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
 102     BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
 103     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 104     BOOST_CHECK_GE(fabs(error), fabs(Q - ((Real)17 + half<Real>()*third<Real>())));
 105 }
 106
 107 // Examples taken from
 108 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
 109 template<class Real, unsigned Points>
 110 void test_ca()
 111 {
 112     std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 113     Real tol = boost::math::tools::epsilon<Real>() * 10;
 114     Real L1;
 115     Real error;
 116
 117     auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
 118     Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
 119     Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
 120     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 121     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 122     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 123     BOOST_CHECK_GE(fabs(error), fabs(Q - Q_expected));
 124
 125     auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
 126     Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 15, get_termination_condition<Real>(), &error, &L1);
 127     Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
 128     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 129     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 130     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 131     BOOST_CHECK_GE(fabs(error), fabs(Q - Q_expected));
 132
 133     auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
 134     Q = gauss_kronrod<Real, Points>::integrate(f5, 0, 1, 25);
 135     Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
 136     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 100 * tol);
 137 }
 138
 139 template<class Real, unsigned Points>
 140 void test_three_quadrature_schemes_examples()
 141 {
 142     std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 143     Real tol = boost::math::tools::epsilon<Real>() * 10;
 144     Real Q;
 145     Real Q_expected;
 146
 147     // Example 1:
 148     auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
 149     Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1);
 150     Q_expected = half<Real>()*half<Real>();
 151     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 152
 153
 154     // Example 2:
 155     auto f2 = [](const Real& t) { return t*t*atan(t); };
 156     Q = gauss_kronrod<Real, Points>::integrate(f2, 0, 1);
 157     Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
 158     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
 159
 160     // Example 3:
 161     auto f3 = [](const Real& t) { return exp(t)*cos(t); };
 162     Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>());
 163     Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
 164     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 165
 166     // Example 4:
 167     auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
 168     Q = gauss_kronrod<Real, Points>::integrate(f4, 0, 1);
 169     Q_expected = 5*pi<Real>()*pi<Real>()/96;
 170     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 171 }
 172
 173
 174 template<class Real, unsigned Points>
 175 void test_integration_over_real_line()
 176 {
 177     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";
 178     Real tol = boost::math::tools::epsilon<Real>() * 10;
 179     Real Q;
 180     Real Q_expected;
 181     Real L1;
 182     Real error;
 183
 184     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 185     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);
 186     Q_expected = pi<Real>();
 187     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 188     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 189     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 190
 191     auto f4 = [](const Real& t) { return 1/cosh(t);};
 192     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);
 193     Q_expected = pi<Real>();
 194     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 195     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 196     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 197
 198 }
 199
 200 template<class Real, unsigned Points>
 201 void test_right_limit_infinite()
 202 {
 203     std::cout << "Testing right limit infinite for Gauss-Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 204     Real tol = boost::math::tools::epsilon<Real>() * 10;
 205     Real Q;
 206     Real Q_expected;
 207     Real L1;
 208     Real error;
 209
 210     // Example 11:
 211     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 212     Q = gauss_kronrod<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
 213     Q_expected = half_pi<Real>();
 214     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 215     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 216
 217     auto f4 = [](const Real& t) { return 1/(1+t*t); };
 218     Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
 219     Q_expected = pi<Real>()/4;
 220     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 221     BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
 222 }
 223
 224 template<class Real, unsigned Points>
 225 void test_left_limit_infinite()
 226 {
 227     std::cout << "Testing left limit infinite for Gauss-Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 228     Real tol = boost::math::tools::epsilon<Real>() * 10;
 229     Real Q;
 230     Real Q_expected;
 231
 232     // Example 11:
 233     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 234     Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), 0);
 235     Q_expected = half_pi<Real>();
 236     BOOST_CHECK_CLOSE(Q, Q_expected, 300*tol);
 237 }
 238
 239 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
 240 {
 241 #ifdef TEST1
 242     std::cout << "Testing with 15 point Gauss-Kronrod rule:\n";
 243     test_linear<double, 15>();
 244     test_quadratic<double, 15>();
 245     test_ca<double, 15>();
 246     test_three_quadrature_schemes_examples<double, 15>();
 247     test_integration_over_real_line<double, 15>();
 248     test_right_limit_infinite<double, 15>();
 249     test_left_limit_infinite<double, 15>();
 250
 251     //  test one case where we do not have pre-computed constants:
 252     std::cout << "Testing with 17 point Gauss-Kronrod rule:\n";
 253     test_linear<double, 17>();
 254     test_quadratic<double, 17>();
 255     test_ca<double, 17>();
 256     test_three_quadrature_schemes_examples<double, 17>();
 257     test_integration_over_real_line<double, 17>();
 258     test_right_limit_infinite<double, 17>();
 259     test_left_limit_infinite<double, 17>();
 260 #endif
 261 #ifdef TEST1A
 262     std::cout << "Testing with 21 point Gauss-Kronrod rule:\n";
 263     test_linear<cpp_bin_float_quad, 21>();
 264     test_quadratic<cpp_bin_float_quad, 21>();
 265     test_ca<cpp_bin_float_quad, 21>();
 266     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
 267     test_integration_over_real_line<cpp_bin_float_quad, 21>();
 268     test_right_limit_infinite<cpp_bin_float_quad, 21>();
 269     test_left_limit_infinite<cpp_bin_float_quad, 21>();
 270
 271     std::cout << "Testing with 31 point Gauss-Kronrod rule:\n";
 272     test_linear<cpp_bin_float_quad, 31>();
 273     test_quadratic<cpp_bin_float_quad, 31>();
 274     test_ca<cpp_bin_float_quad, 31>();
 275     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
 276     test_integration_over_real_line<cpp_bin_float_quad, 31>();
 277     test_right_limit_infinite<cpp_bin_float_quad, 31>();
 278     test_left_limit_infinite<cpp_bin_float_quad, 31>();
 279 #endif
 280 #ifdef TEST2
 281     std::cout << "Testing with 41 point Gauss-Kronrod rule:\n";
 282     test_linear<cpp_bin_float_quad, 41>();
 283     test_quadratic<cpp_bin_float_quad, 41>();
 284     test_ca<cpp_bin_float_quad, 41>();
 285     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
 286     test_integration_over_real_line<cpp_bin_float_quad, 41>();
 287     test_right_limit_infinite<cpp_bin_float_quad, 41>();
 288     test_left_limit_infinite<cpp_bin_float_quad, 41>();
 289
 290     std::cout << "Testing with 51 point Gauss-Kronrod rule:\n";
 291     test_linear<cpp_bin_float_quad, 51>();
 292     test_quadratic<cpp_bin_float_quad, 51>();
 293     test_ca<cpp_bin_float_quad, 51>();
 294     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 51>();
 295     test_integration_over_real_line<cpp_bin_float_quad, 51>();
 296     test_right_limit_infinite<cpp_bin_float_quad, 51>();
 297     test_left_limit_infinite<cpp_bin_float_quad, 51>();
 298 #endif
 299 #ifdef TEST3
 300     std::cout << "Testing with 61 point Gauss-Kronrod rule:\n";
 301     test_linear<cpp_bin_float_quad, 61>();
 302     test_quadratic<cpp_bin_float_quad, 61>();
 303     test_ca<cpp_bin_float_quad, 61>();
 304     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 61>();
 305     test_integration_over_real_line<cpp_bin_float_quad, 61>();
 306     test_right_limit_infinite<cpp_bin_float_quad, 61>();
 307     test_left_limit_infinite<cpp_bin_float_quad, 61>();
 308 #endif
 309 }
 310
 311 #else
 312
 313 int main() { return 0; }
 314
 315 #endif