libs/math/test/test_trapezoidal.cpp

   1 /*
   2  * Copyright Nick Thompson, 2017
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7 #define BOOST_TEST_MODULE trapezoidal_quadrature
   8
   9 #include <complex>
  10 #include <boost/config.hpp>
  11 //#include <boost/multiprecision/mpc.hpp>
  12 #include <boost/test/included/unit_test.hpp>
  13 #include <boost/test/tools/floating_point_comparison.hpp>
  14 #include <boost/math/concepts/real_concept.hpp>
  15 #include <boost/math/special_functions/bessel.hpp>
  16 #include <boost/math/quadrature/trapezoidal.hpp>
  17 #include <boost/multiprecision/cpp_bin_float.hpp>
  18 #include <boost/multiprecision/cpp_dec_float.hpp>
  19 #ifdef BOOST_HAS_FLOAT128
  20 #include <boost/multiprecision/complex128.hpp>
  21 #endif
  22 using boost::multiprecision::cpp_bin_float_50;
  23 using boost::multiprecision::cpp_bin_float_100;
  24 using boost::math::quadrature::trapezoidal;
  25
  26 // These tests come from:
  27 // https://doi.org/10.1023/A:1025524324969
  28 // "Computing special functions by using quadrature rules",  Gil, Segura, and Temme.
  29 template<class Complex>
  30 void test_complex_bessel()
  31 {
  32     std::cout << "Testing that complex-valued integrands are integrated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Complex>().pretty_name()  << "\n";
  33     typedef typename Complex::value_type Real;
  34     Complex z{2, 3};
  35     int n = 2;
  36     using boost::math::constants::pi;
  37     auto bessel_integrand = [&n, &z](Real theta)->Complex
  38     {
  39         using std::cos;
  40         using std::sin;
  41         Real t1 = sin(theta);
  42         Real t2 = - n*theta;
  43         Complex arg = z*t1 + t2;
  44         return cos(arg)/pi<Real>();
  45     };
  46
  47     using boost::math::quadrature::trapezoidal;
  48
  49     Real a = 0;
  50     Real b = pi<Real>();
  51     Complex Jnz = trapezoidal<decltype(bessel_integrand), Real>(bessel_integrand, a, b);
  52     // N[BesselJ[2, 2 + 3 I], 143]
  53     // 1.257674591970511077630764085052638490387449039392695959943027966195657681586539389134094087028482099931927725892... +
  54     // 2.318771368505683055818032722011594415038779144567369903204833213112457210243098545874099591376455981793627257060... i
  55     Real Jnzx = boost::lexical_cast<Real>("1.257674591970511077630764085052638490387449039392695959943027966195657681586539389134094087028482099931927725892");
  56     Real Jnzy = boost::lexical_cast<Real>("2.318771368505683055818032722011594415038779144567369903204833213112457210243098545874099591376455981793627257060");
  57     Real tol = 10*std::numeric_limits<Real>::epsilon();
  58     BOOST_CHECK_CLOSE_FRACTION(Jnz.real(), Jnzx, tol);
  59     BOOST_CHECK_CLOSE_FRACTION(Jnz.imag(), Jnzy, tol);
  60 }
  61
  62 template<class Complex>
  63 void test_I0_complex()
  64 {
  65     std::cout << "Testing that complex-argument I0 is calculated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Complex>().pretty_name()  << "\n";
  66     typedef typename Complex::value_type Real;
  67     Complex z{2, 3};
  68     using boost::math::constants::pi;
  69     auto I0 = [&z](Real theta)->Complex
  70     {
  71         using std::cos;
  72         using std::exp;
  73         return exp(z*cos(theta))/pi<Real>();
  74     };
  75
  76     using boost::math::quadrature::trapezoidal;
  77
  78     Real a = 0;
  79     Real b = pi<Real>();
  80     Complex I0z = trapezoidal<decltype(I0), Real>(I0, a, b);
  81     // N[BesselI[0, 2 + 3 I], 143]
  82     // -1.24923487960742219637619681391438589436703710701063561548156438052154090067526565701278826317992172207565649925713468090525951417141982808439560899101
  83     // 0.947983792057734776114060623981442199525094227418764823692296622398838765371662384207319492908490909109393495109183270208372778907692930675595924819922 i
  84     Real I0zx = boost::lexical_cast<Real>("-1.24923487960742219637619681391438589436703710701063561548156438052154090067526565701278826317992172207565649925713468090525951417141982808439560899101");
  85     Real I0zy = boost::lexical_cast<Real>("0.947983792057734776114060623981442199525094227418764823692296622398838765371662384207319492908490909109393495109183270208372778907692930675595924819922");
  86     Real tol = 10*std::numeric_limits<Real>::epsilon();
  87     BOOST_CHECK_CLOSE_FRACTION(I0z.real(), I0zx, tol);
  88     BOOST_CHECK_CLOSE_FRACTION(I0z.imag(), I0zy, tol);
  89 }
  90
  91
  92 template<class Complex>
  93 void test_erfc()
  94 {
  95     std::cout << "Testing that complex-argument erfc is calculated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Complex>().pretty_name()  << "\n";
  96     typedef typename Complex::value_type Real;
  97     Complex z{2, -1};
  98     using boost::math::constants::pi;
  99     using boost::math::constants::two_pi;
 100     auto erfc = [&z](Real theta)->Complex
 101     {
 102         using std::exp;
 103         using std::tan;
 104         Real t = tan(theta/2);
 105         Complex arg = -z*z*(1+t*t);
 106         return exp(arg)/two_pi<Real>();
 107     };
 108
 109     using boost::math::quadrature::trapezoidal;
 110
 111     Real a = -pi<Real>();
 112     Real b = pi<Real>();
 113     Complex erfcz = trapezoidal<decltype(erfc), Real>(erfc, a, b, boost::math::tools::root_epsilon<Real>(), 17);
 114     // N[Erfc[2-i], 150]
 115     //-0.00360634272565175091291182820541914235532928536595056623793472801084629874817202107825472707423984408881473019087931573313969503235634965264302640170177
 116     // - 0.0112590060288150250764009156316482248536651598819882163212627394923365188251633729432967232423246312345152595958230197778555210858871376231770868078020 i
 117     Real erfczx = boost::lexical_cast<Real>("-0.00360634272565175091291182820541914235532928536595056623793472801084629874817202107825472707423984408881473019087931573313969503235634965264302640170177");
 118     Real erfczy = boost::lexical_cast<Real>("-0.0112590060288150250764009156316482248536651598819882163212627394923365188251633729432967232423246312345152595958230197778555210858871376231770868078020");
 119     Real tol = 5000*std::numeric_limits<Real>::epsilon();
 120     BOOST_CHECK_CLOSE_FRACTION(erfcz.real(), erfczx, tol);
 121     BOOST_CHECK_CLOSE_FRACTION(erfcz.imag(), erfczy, tol);
 122 }
 123
 124
 125 template<class Real>
 126 void test_constant()
 127 {
 128     std::cout << "Testing constants are integrated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 129
 130     auto f = [](Real)->Real { return boost::math::constants::half<Real>(); };
 131     Real Q = trapezoidal<decltype(f), Real>(f, (Real) 0.0, (Real) 10.0);
 132     BOOST_CHECK_CLOSE(Q, 5.0, 100*std::numeric_limits<Real>::epsilon());
 133 }
 134
 135
 136 template<class Real>
 137 void test_rational_periodic()
 138 {
 139     using boost::math::constants::two_pi;
 140     using boost::math::constants::third;
 141     std::cout << "Testing that rational periodic functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 142
 143     auto f = [](Real x)->Real { return 1/(5 - 4*cos(x)); };
 144
 145     Real tol = 100*boost::math::tools::epsilon<Real>();
 146     Real Q = trapezoidal(f, (Real) 0.0, two_pi<Real>(), tol);
 147
 148     BOOST_CHECK_CLOSE_FRACTION(Q, two_pi<Real>()*third<Real>(), 10*tol);
 149 }
 150
 151 template<class Real>
 152 void test_bump_function()
 153 {
 154     std::cout << "Testing that bump functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 155     auto f = [](Real x)->Real {
 156         if( x>= 1 || x <= -1)
 157         {
 158             return (Real) 0;
 159         }
 160         return (Real) exp(-(Real) 1/(1-x*x));
 161     };
 162     Real tol = boost::math::tools::epsilon<Real>();
 163     Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
 164     // 2*NIntegrate[Exp[-(1/(1 - x^2))], {x, 0, 1}, WorkingPrecision -> 210]
 165     Real Q_exp = boost::lexical_cast<Real>("0.44399381616807943782304892117055266376120178904569749730748455394704");
 166     BOOST_CHECK_CLOSE_FRACTION(Q, Q_exp, 30*tol);
 167 }
 168
 169 template<class Real>
 170 void test_zero_function()
 171 {
 172     std::cout << "Testing that zero functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 173     auto f = [](Real)->Real { return (Real) 0;};
 174     Real tol = 100* boost::math::tools::epsilon<Real>();
 175     Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
 176     BOOST_CHECK_SMALL(Q, 100*tol);
 177 }
 178
 179 template<class Real>
 180 void test_sinsq()
 181 {
 182     std::cout << "Testing that sin(x)^2 is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 183     auto f = [](Real x)->Real { return sin(10*x)*sin(10*x); };
 184     Real tol = 100* boost::math::tools::epsilon<Real>();
 185     Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi<Real>(), tol);
 186     BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>(), tol);
 187
 188 }
 189
 190 template<class Real>
 191 void test_slowly_converging()
 192 {
 193     using std::sqrt;
 194     std::cout << "Testing that non-periodic functions are correctly integrated by the trapezoidal rule, even if slowly, on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 195     // This function is not periodic, so it should not be fast to converge:
 196     auto f = [](Real x)->Real { using std::sqrt;  return sqrt(1 - x*x); };
 197
 198     Real tol = sqrt(sqrt(boost::math::tools::epsilon<Real>()));
 199     Real error_estimate;
 200     Real Q = trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate);
 201     BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>()/2, 10*tol);
 202 }
 203
 204 template<class Real>
 205 void test_rational_sin()
 206 {
 207     using std::pow;
 208     using std::sin;
 209     using boost::math::constants::two_pi;
 210     using boost::math::constants::half;
 211     std::cout << "Testing that a rational sin function is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 212     Real a = 5;
 213     auto f = [=](Real x)->Real { using std::sin;  Real t = a + sin(x); return 1.0f / (t*t); };
 214
 215     Real expected = two_pi<Real>()*a/pow(a*a - 1, 3*half<Real>());
 216     Real tol = 100* boost::math::tools::epsilon<Real>();
 217     Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi<Real>(), tol);
 218     BOOST_CHECK_CLOSE_FRACTION(Q, expected, tol);
 219 }
 220
 221 BOOST_AUTO_TEST_CASE(trapezoidal_quadrature)
 222 {
 223     test_constant<float>();
 224     test_constant<double>();
 225     test_constant<long double>();
 226     test_constant<boost::math::concepts::real_concept>();
 227     test_constant<cpp_bin_float_50>();
 228     test_constant<cpp_bin_float_100>();
 229
 230     test_rational_periodic<float>();
 231     test_rational_periodic<double>();
 232     test_rational_periodic<long double>();
 233     test_rational_periodic<boost::math::concepts::real_concept>();
 234     test_rational_periodic<cpp_bin_float_50>();
 235     test_rational_periodic<cpp_bin_float_100>();
 236
 237     test_bump_function<float>();
 238     test_bump_function<double>();
 239     test_bump_function<long double>();
 240     test_rational_periodic<boost::math::concepts::real_concept>();
 241     test_rational_periodic<cpp_bin_float_50>();
 242
 243     test_zero_function<float>();
 244     test_zero_function<double>();
 245     test_zero_function<long double>();
 246     test_zero_function<boost::math::concepts::real_concept>();
 247     test_zero_function<cpp_bin_float_50>();
 248     test_zero_function<cpp_bin_float_100>();
 249
 250     test_sinsq<float>();
 251     test_sinsq<double>();
 252     test_sinsq<long double>();
 253     test_sinsq<boost::math::concepts::real_concept>();
 254     test_sinsq<cpp_bin_float_50>();
 255     test_sinsq<cpp_bin_float_100>();
 256
 257     test_slowly_converging<float>();
 258     test_slowly_converging<double>();
 259     test_slowly_converging<long double>();
 260     test_slowly_converging<boost::math::concepts::real_concept>();
 261
 262     test_rational_sin<float>();
 263     test_rational_sin<double>();
 264     test_rational_sin<long double>();
 265     //test_rational_sin<boost::math::concepts::real_concept>();
 266     test_rational_sin<cpp_bin_float_50>();
 267
 268     test_complex_bessel<std::complex<float>>();
 269     test_complex_bessel<std::complex<double>>();
 270     test_complex_bessel<std::complex<long double>>();
 271     //test_complex_bessel<boost::multiprecision::mpc_complex_100>();
 272     test_I0_complex<std::complex<float>>();
 273     test_I0_complex<std::complex<double>>();
 274     test_I0_complex<std::complex<long double>>();
 275     //test_I0_complex<boost::multiprecision::mpc_complex_100>();
 276     test_erfc<std::complex<float>>();
 277     test_erfc<std::complex<double>>();
 278     test_erfc<std::complex<long double>>();
 279     //test_erfc<boost::multiprecision::number<boost::multiprecision::mpc_complex_backend<20>>>();
 280     //test_erfc<boost::multiprecision::number<boost::multiprecision::mpc_complex_backend<30>>>();
 281     //test_erfc<boost::multiprecision::mpc_complex_50>();
 282     //test_erfc<boost::multiprecision::mpc_complex_100>();
 283
 284 #ifdef BOOST_HAS_FLOAT128
 285     test_complex_bessel<boost::multiprecision::complex128>();
 286     test_I0_complex<boost::multiprecision::complex128>();
 287     test_erfc<boost::multiprecision::complex128>();
 288 #endif
 289
 290 }