libs/math/test/gauss_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 tanh_sinh_quadrature_test
   8
   9 #include <complex>
  10 //#include <boost/multiprecision/mpc.hpp>
  11 #include <boost/config.hpp>
  12 #include <boost/detail/workaround.hpp>
  13
  14 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
  15
  16 #include <boost/math/concepts/real_concept.hpp>
  17 #include <boost/test/included/unit_test.hpp>
  18 #include <boost/test/tools/floating_point_comparison.hpp>
  19 #include <boost/math/quadrature/gauss.hpp>
  20 #include <boost/math/special_functions/sinc.hpp>
  21 #include <boost/multiprecision/cpp_bin_float.hpp>
  22 #include <boost/multiprecision/cpp_complex.hpp>
  23
  24 #ifdef BOOST_HAS_FLOAT128
  25 #include <boost/multiprecision/complex128.hpp>
  26 #endif
  27
  28 #ifdef _MSC_VER
  29 #pragma warning(disable:4127)  // Conditional expression is constant
  30 #endif
  31
  32 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3)
  33 #  define TEST1
  34 #  define TEST2
  35 #  define TEST3
  36 #endif
  37
  38 using std::expm1;
  39 using std::atan;
  40 using std::tan;
  41 using std::log;
  42 using std::log1p;
  43 using std::asinh;
  44 using std::atanh;
  45 using std::sqrt;
  46 using std::isnormal;
  47 using std::abs;
  48 using std::sinh;
  49 using std::tanh;
  50 using std::cosh;
  51 using std::pow;
  52 using std::exp;
  53 using std::sin;
  54 using std::cos;
  55 using std::string;
  56 using boost::math::quadrature::gauss;
  57 using boost::math::constants::pi;
  58 using boost::math::constants::half_pi;
  59 using boost::math::constants::two_div_pi;
  60 using boost::math::constants::two_pi;
  61 using boost::math::constants::half;
  62 using boost::math::constants::third;
  63 using boost::math::constants::half;
  64 using boost::math::constants::third;
  65 using boost::math::constants::catalan;
  66 using boost::math::constants::ln_two;
  67 using boost::math::constants::root_two;
  68 using boost::math::constants::root_two_pi;
  69 using boost::math::constants::root_pi;
  70 using boost::multiprecision::cpp_bin_float_quad;
  71
  72 //
  73 // Error rates depend only on the number of points in the approximation, not the type being tested,
  74 // define all our expected errors here:
  75 //
  76
  77 enum
  78 {
  79    test_ca_error_id,
  80    test_ca_error_id_2,
  81    test_three_quad_error_id,
  82    test_three_quad_error_id_2,
  83    test_integration_over_real_line_error_id,
  84    test_right_limit_infinite_error_id,
  85    test_left_limit_infinite_error_id
  86 };
  87
  88 template <unsigned Points>
  89 double expected_error(unsigned)
  90 {
  91    return 0; // placeholder, all tests will fail
  92 }
  93
  94 template <>
  95 double expected_error<7>(unsigned id)
  96 {
  97    switch (id)
  98    {
  99    case test_ca_error_id:
 100       return 1e-7;
 101    case test_ca_error_id_2:
 102       return 2e-5;
 103    case test_three_quad_error_id:
 104       return 1e-8;
 105    case test_three_quad_error_id_2:
 106       return 3.5e-3;
 107    case test_integration_over_real_line_error_id:
 108       return 6e-3;
 109    case test_right_limit_infinite_error_id:
 110    case test_left_limit_infinite_error_id:
 111       return 1e-5;
 112    }
 113    return 0;  // placeholder, all tests will fail
 114 }
 115
 116 template <>
 117 double expected_error<9>(unsigned id)
 118 {
 119    switch (id)
 120    {
 121    case test_ca_error_id:
 122       return 1e-7;
 123    case test_ca_error_id_2:
 124       return 2e-5;
 125    case test_three_quad_error_id:
 126       return 1e-8;
 127    case test_three_quad_error_id_2:
 128       return 3.5e-3;
 129    case test_integration_over_real_line_error_id:
 130       return 6e-3;
 131    case test_right_limit_infinite_error_id:
 132    case test_left_limit_infinite_error_id:
 133       return 1e-5;
 134    }
 135    return 0;  // placeholder, all tests will fail
 136 }
 137
 138 template <>
 139 double expected_error<10>(unsigned id)
 140 {
 141    switch (id)
 142    {
 143    case test_ca_error_id:
 144       return 1e-12;
 145    case test_ca_error_id_2:
 146       return 3e-6;
 147    case test_three_quad_error_id:
 148       return 2e-13;
 149    case test_three_quad_error_id_2:
 150       return 2e-3;
 151    case test_integration_over_real_line_error_id:
 152       return 6e-3;  // doesn't get any better with more points!
 153    case test_right_limit_infinite_error_id:
 154    case test_left_limit_infinite_error_id:
 155       return 5e-8;
 156    }
 157    return 0;  // placeholder, all tests will fail
 158 }
 159
 160 template <>
 161 double expected_error<15>(unsigned id)
 162 {
 163    switch (id)
 164    {
 165    case test_ca_error_id:
 166       return 6e-20;
 167    case test_ca_error_id_2:
 168       return 3e-7;
 169    case test_three_quad_error_id:
 170       return 1e-19;
 171    case test_three_quad_error_id_2:
 172       return 6e-4;
 173    case test_integration_over_real_line_error_id:
 174       return 6e-3;  // doesn't get any better with more points!
 175    case test_right_limit_infinite_error_id:
 176    case test_left_limit_infinite_error_id:
 177       return 5e-11;
 178    }
 179    return 0;  // placeholder, all tests will fail
 180 }
 181
 182 template <>
 183 double expected_error<20>(unsigned id)
 184 {
 185    switch (id)
 186    {
 187    case test_ca_error_id:
 188       return 1e-26;
 189    case test_ca_error_id_2:
 190       return 1e-7;
 191    case test_three_quad_error_id:
 192       return 3e-27;
 193    case test_three_quad_error_id_2:
 194       return 3e-4;
 195    case test_integration_over_real_line_error_id:
 196       return 5e-5;  // doesn't get any better with more points!
 197    case test_right_limit_infinite_error_id:
 198    case test_left_limit_infinite_error_id:
 199       return 1e-15;
 200    }
 201    return 0;  // placeholder, all tests will fail
 202 }
 203
 204 template <>
 205 double expected_error<25>(unsigned id)
 206 {
 207    switch (id)
 208    {
 209    case test_ca_error_id:
 210       return 5e-33;
 211    case test_ca_error_id_2:
 212       return 1e-8;
 213    case test_three_quad_error_id:
 214       return 1e-32;
 215    case test_three_quad_error_id_2:
 216       return 3e-4;
 217    case test_integration_over_real_line_error_id:
 218       return 1e-14;
 219    case test_right_limit_infinite_error_id:
 220    case test_left_limit_infinite_error_id:
 221       return 3e-19;
 222    }
 223    return 0;  // placeholder, all tests will fail
 224 }
 225
 226 template <>
 227 double expected_error<30>(unsigned id)
 228 {
 229    switch (id)
 230    {
 231    case test_ca_error_id:
 232       return 2e-34;
 233    case test_ca_error_id_2:
 234       return 5e-9;
 235    case test_three_quad_error_id:
 236       return 4e-34;
 237    case test_three_quad_error_id_2:
 238       return 1e-4;
 239    case test_integration_over_real_line_error_id:
 240       return 1e-16;
 241    case test_right_limit_infinite_error_id:
 242    case test_left_limit_infinite_error_id:
 243       return 3e-23;
 244    }
 245    return 0;  // placeholder, all tests will fail
 246 }
 247
 248
 249 template<class Real, unsigned Points>
 250 void test_linear()
 251 {
 252     std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 253     Real tol = boost::math::tools::epsilon<Real>() * 10;
 254     auto f = [](const Real& x)
 255     {
 256        return 5*x + 7;
 257     };
 258     Real L1;
 259     Real Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 1, &L1);
 260     BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
 261     BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
 262 }
 263
 264 template<class Real, unsigned Points>
 265 void test_quadratic()
 266 {
 267     std::cout << "Testing quadratic functions are integrated properly by Gaussian quadrature on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 268     Real tol = boost::math::tools::epsilon<Real>() * 10;
 269
 270     auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
 271     Real L1;
 272     Real Q = gauss<Real, Points>::integrate(f, 0, 1, &L1);
 273     BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
 274     BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
 275 }
 276
 277 // Examples taken from
 278 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
 279 template<class Real, unsigned Points>
 280 void test_ca()
 281 {
 282     std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 283     Real tol = expected_error<Points>(test_ca_error_id);
 284     Real L1;
 285
 286     auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
 287     Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
 288     Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
 289     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 290     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 291
 292     auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
 293     Q = gauss<Real, Points>::integrate(f2, 0 , 1, &L1);
 294     Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
 295     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 296     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 297
 298     tol = expected_error<Points>(test_ca_error_id_2);
 299     auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
 300     Q = gauss<Real, Points>::integrate(f5, 0 , 1);
 301     Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
 302     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 303 }
 304
 305 template<class Real, unsigned Points>
 306 void test_three_quadrature_schemes_examples()
 307 {
 308     std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 309     Real tol = expected_error<Points>(test_three_quad_error_id);
 310     Real Q;
 311     Real Q_expected;
 312
 313     // Example 1:
 314     auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
 315     Q = gauss<Real, Points>::integrate(f1, 0 , 1);
 316     Q_expected = half<Real>()*half<Real>();
 317     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 318
 319
 320     // Example 2:
 321     auto f2 = [](const Real& t) { return t*t*atan(t); };
 322     Q = gauss<Real, Points>::integrate(f2, 0 , 1);
 323     Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
 324     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
 325
 326     // Example 3:
 327     auto f3 = [](const Real& t) { return exp(t)*cos(t); };
 328     Q = gauss<Real, Points>::integrate(f3, 0, half_pi<Real>());
 329     Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
 330     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 331
 332     // Example 4:
 333     auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
 334     Q = gauss<Real, Points>::integrate(f4, 0 , 1);
 335     Q_expected = 5*pi<Real>()*pi<Real>()/96;
 336     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 337
 338     tol = expected_error<Points>(test_three_quad_error_id_2);
 339     // Example 5:
 340     auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
 341     Q = gauss<Real, Points>::integrate(f5, 0 , 1);
 342     Q_expected = -4/ (Real) 9;
 343     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 344
 345     // Example 6:
 346     auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
 347     Q = gauss<Real, Points>::integrate(f6, 0 , 1);
 348     Q_expected = pi<Real>()/4;
 349     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 350 }
 351
 352
 353 template<class Real, unsigned Points>
 354 void test_integration_over_real_line()
 355 {
 356     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";
 357     Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
 358     Real Q;
 359     Real Q_expected;
 360     Real L1;
 361
 362     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 363     Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
 364     Q_expected = pi<Real>();
 365     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 366     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 367 }
 368
 369 template<class Real, unsigned Points>
 370 void test_right_limit_infinite()
 371 {
 372     std::cout << "Testing right limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 373     Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
 374     Real Q;
 375     Real Q_expected;
 376     Real L1;
 377
 378     // Example 11:
 379     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 380     Q = gauss<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), &L1);
 381     Q_expected = half_pi<Real>();
 382     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 383
 384     auto f4 = [](const Real& t) { return 1/(1+t*t); };
 385     Q = gauss<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), &L1);
 386     Q_expected = pi<Real>()/4;
 387     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 388 }
 389
 390 template<class Real, unsigned Points>
 391 void test_left_limit_infinite()
 392 {
 393     std::cout << "Testing left limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 394     Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
 395     Real Q;
 396     Real Q_expected;
 397
 398     // Example 11:
 399     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 400     Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0));
 401     Q_expected = half_pi<Real>();
 402     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 403 }
 404
 405 template<class Complex>
 406 void test_complex_lambert_w()
 407 {
 408     std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
 409     typedef typename Complex::value_type Real;
 410     Real tol = 10e-9;
 411     using boost::math::constants::pi;
 412     Complex z{2, 3};
 413     auto lw = [&z](Real v)->Complex {
 414       using std::cos;
 415       using std::sin;
 416       using std::exp;
 417       Real sinv = sin(v);
 418       Real cosv = cos(v);
 419
 420       Real cotv = cosv/sinv;
 421       Real cscv = 1/sinv;
 422       Real t = (1-v*cotv)*(1-v*cotv) + v*v;
 423       Real x = v*cscv*exp(-v*cotv);
 424       Complex den = z + x;
 425       Complex num = t*(z/pi<Real>());
 426       Complex res = num/den;
 427       return res;
 428     };
 429
 430     //N[ProductLog[2+3*I], 150]
 431     Complex Q = gauss<Real, 30>::integrate(lw, (Real) 0, pi<Real>());
 432     BOOST_CHECK_CLOSE_FRACTION(Q.real(), boost::lexical_cast<Real>("1.09007653448579084630177782678166964987102108635357778056449870727913321296238687023915522935120701763447787503167111962008709116746523970476893277703"), tol);
 433     BOOST_CHECK_CLOSE_FRACTION(Q.imag(), boost::lexical_cast<Real>("0.530139720774838801426860213574121741928705631382703178297940568794784362495390544411799468140433404536019992695815009036975117285537382995180319280835"), tol);
 434 }
 435
 436 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
 437 {
 438
 439 #ifdef TEST1
 440     test_linear<double, 7>();
 441     test_quadratic<double, 7>();
 442     test_ca<double, 7>();
 443     test_three_quadrature_schemes_examples<double, 7>();
 444     test_integration_over_real_line<double, 7>();
 445     test_right_limit_infinite<double, 7>();
 446     test_left_limit_infinite<double, 7>();
 447
 448     test_linear<double, 9>();
 449     test_quadratic<double, 9>();
 450     test_ca<double, 9>();
 451     test_three_quadrature_schemes_examples<double, 9>();
 452     test_integration_over_real_line<double, 9>();
 453     test_right_limit_infinite<double, 9>();
 454     test_left_limit_infinite<double, 9>();
 455
 456     test_linear<cpp_bin_float_quad, 10>();
 457     test_quadratic<cpp_bin_float_quad, 10>();
 458     test_ca<cpp_bin_float_quad, 10>();
 459     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 10>();
 460     test_integration_over_real_line<cpp_bin_float_quad, 10>();
 461     test_right_limit_infinite<cpp_bin_float_quad, 10>();
 462     test_left_limit_infinite<cpp_bin_float_quad, 10>();
 463 #endif
 464 #ifdef TEST2
 465     test_linear<cpp_bin_float_quad, 15>();
 466     test_quadratic<cpp_bin_float_quad, 15>();
 467     test_ca<cpp_bin_float_quad, 15>();
 468     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 15>();
 469     test_integration_over_real_line<cpp_bin_float_quad, 15>();
 470     test_right_limit_infinite<cpp_bin_float_quad, 15>();
 471     test_left_limit_infinite<cpp_bin_float_quad, 15>();
 472
 473     test_linear<cpp_bin_float_quad, 20>();
 474     test_quadratic<cpp_bin_float_quad, 20>();
 475     test_ca<cpp_bin_float_quad, 20>();
 476     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 20>();
 477     test_integration_over_real_line<cpp_bin_float_quad, 20>();
 478     test_right_limit_infinite<cpp_bin_float_quad, 20>();
 479     test_left_limit_infinite<cpp_bin_float_quad, 20>();
 480
 481     test_linear<cpp_bin_float_quad, 25>();
 482     test_quadratic<cpp_bin_float_quad, 25>();
 483     test_ca<cpp_bin_float_quad, 25>();
 484     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 25>();
 485     test_integration_over_real_line<cpp_bin_float_quad, 25>();
 486     test_right_limit_infinite<cpp_bin_float_quad, 25>();
 487     test_left_limit_infinite<cpp_bin_float_quad, 25>();
 488
 489     test_linear<cpp_bin_float_quad, 30>();
 490     test_quadratic<cpp_bin_float_quad, 30>();
 491     test_ca<cpp_bin_float_quad, 30>();
 492     test_three_quadrature_schemes_examples<cpp_bin_float_quad, 30>();
 493     test_integration_over_real_line<cpp_bin_float_quad, 30>();
 494     test_right_limit_infinite<cpp_bin_float_quad, 30>();
 495     test_left_limit_infinite<cpp_bin_float_quad, 30>();
 496
 497
 498 #endif
 499 #ifdef TEST3
 500     test_left_limit_infinite<cpp_bin_float_quad, 30>();
 501     test_complex_lambert_w<std::complex<double>>();
 502     test_complex_lambert_w<std::complex<long double>>();
 503 #ifdef BOOST_HAS_FLOAT128
 504     test_left_limit_infinite<boost::multiprecision::float128, 30>();
 505     test_complex_lambert_w<boost::multiprecision::complex128>();
 506 #endif
 507     test_complex_lambert_w<boost::multiprecision::cpp_complex_quad>();
 508 #endif
 509 }
 510
 511 #else
 512
 513 int main() { return 0; }
 514
 515 #endif