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>
test_complex_bessel()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>
test_I0_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>
test_erfc()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>
test_constant()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 Q = trapezoidal<decltype(f), Real>(f, (Real) 10.0, (Real) 0.0);
134 BOOST_CHECK_CLOSE(Q, -5.0, 100*std::numeric_limits<Real>::epsilon());
135
136 Q = trapezoidal<decltype(f), Real>(f, (Real) 10.0, (Real) 10.0);
137 BOOST_CHECK_CLOSE(Q, Real(0), 100*std::numeric_limits<Real>::epsilon());
138 }
139
140
141 template<class Real>
test_rational_periodic()142 void test_rational_periodic()
143 {
144 using boost::math::constants::two_pi;
145 using boost::math::constants::third;
146 std::cout << "Testing that rational periodic functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
147
148 auto f = [](Real x)->Real { return 1/(5 - 4*cos(x)); };
149
150 Real tol = 100*boost::math::tools::epsilon<Real>();
151 Real Q = trapezoidal(f, (Real) 0.0, two_pi<Real>(), tol);
152
153 BOOST_CHECK_CLOSE_FRACTION(Q, two_pi<Real>()*third<Real>(), 10*tol);
154 }
155
156 template<class Real>
test_bump_function()157 void test_bump_function()
158 {
159 std::cout << "Testing that bump functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
160 auto f = [](Real x)->Real {
161 if( x>= 1 || x <= -1)
162 {
163 return (Real) 0;
164 }
165 return (Real) exp(-(Real) 1/(1-x*x));
166 };
167 Real tol = boost::math::tools::epsilon<Real>();
168 Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
169 // 2*NIntegrate[Exp[-(1/(1 - x^2))], {x, 0, 1}, WorkingPrecision -> 210]
170 Real Q_exp = boost::lexical_cast<Real>("0.44399381616807943782304892117055266376120178904569749730748455394704");
171 BOOST_CHECK_CLOSE_FRACTION(Q, Q_exp, 30*tol);
172 }
173
174 template<class Real>
test_zero_function()175 void test_zero_function()
176 {
177 std::cout << "Testing that zero functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
178 auto f = [](Real)->Real { return (Real) 0;};
179 Real tol = 100* boost::math::tools::epsilon<Real>();
180 Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
181 BOOST_CHECK_SMALL(Q, 100*tol);
182 }
183
184 template<class Real>
test_sinsq()185 void test_sinsq()
186 {
187 std::cout << "Testing that sin(x)^2 is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
188 auto f = [](Real x)->Real { return sin(10*x)*sin(10*x); };
189 Real tol = 100* boost::math::tools::epsilon<Real>();
190 Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi<Real>(), tol);
191 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>(), tol);
192
193 }
194
195 template<class Real>
test_slowly_converging()196 void test_slowly_converging()
197 {
198 using std::sqrt;
199 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";
200 // This function is not periodic, so it should not be fast to converge:
201 auto f = [](Real x)->Real { using std::sqrt; return sqrt(1 - x*x); };
202
203 Real tol = sqrt(sqrt(boost::math::tools::epsilon<Real>()));
204 Real error_estimate;
205 Real Q = trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate);
206 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>()/2, 10*tol);
207 }
208
209 template<class Real>
test_rational_sin()210 void test_rational_sin()
211 {
212 using std::pow;
213 using std::sin;
214 using boost::math::constants::two_pi;
215 using boost::math::constants::half;
216 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";
217 Real a = 5;
218 auto f = [=](Real x)->Real { using std::sin; Real t = a + sin(x); return 1.0f / (t*t); };
219
220 Real expected = two_pi<Real>()*a/pow(a*a - 1, 3*half<Real>());
221 Real tol = 100* boost::math::tools::epsilon<Real>();
222 Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi<Real>(), tol);
223 BOOST_CHECK_CLOSE_FRACTION(Q, expected, tol);
224 }
225
BOOST_AUTO_TEST_CASE(trapezoidal_quadrature)226 BOOST_AUTO_TEST_CASE(trapezoidal_quadrature)
227 {
228 test_constant<float>();
229 test_constant<double>();
230 test_constant<long double>();
231 test_constant<boost::math::concepts::real_concept>();
232 test_constant<cpp_bin_float_50>();
233 test_constant<cpp_bin_float_100>();
234
235 test_rational_periodic<float>();
236 test_rational_periodic<double>();
237 test_rational_periodic<long double>();
238 test_rational_periodic<boost::math::concepts::real_concept>();
239 test_rational_periodic<cpp_bin_float_50>();
240 test_rational_periodic<cpp_bin_float_100>();
241
242 test_bump_function<float>();
243 test_bump_function<double>();
244 test_bump_function<long double>();
245 test_rational_periodic<boost::math::concepts::real_concept>();
246 test_rational_periodic<cpp_bin_float_50>();
247
248 test_zero_function<float>();
249 test_zero_function<double>();
250 test_zero_function<long double>();
251 test_zero_function<boost::math::concepts::real_concept>();
252 test_zero_function<cpp_bin_float_50>();
253 test_zero_function<cpp_bin_float_100>();
254
255 test_sinsq<float>();
256 test_sinsq<double>();
257 test_sinsq<long double>();
258 test_sinsq<boost::math::concepts::real_concept>();
259 test_sinsq<cpp_bin_float_50>();
260 test_sinsq<cpp_bin_float_100>();
261
262 test_slowly_converging<float>();
263 test_slowly_converging<double>();
264 test_slowly_converging<long double>();
265 test_slowly_converging<boost::math::concepts::real_concept>();
266
267 test_rational_sin<float>();
268 test_rational_sin<double>();
269 test_rational_sin<long double>();
270 //test_rational_sin<boost::math::concepts::real_concept>();
271 test_rational_sin<cpp_bin_float_50>();
272
273 test_complex_bessel<std::complex<float>>();
274 test_complex_bessel<std::complex<double>>();
275 test_complex_bessel<std::complex<long double>>();
276 //test_complex_bessel<boost::multiprecision::mpc_complex_100>();
277 test_I0_complex<std::complex<float>>();
278 test_I0_complex<std::complex<double>>();
279 test_I0_complex<std::complex<long double>>();
280 //test_I0_complex<boost::multiprecision::mpc_complex_100>();
281 test_erfc<std::complex<float>>();
282 test_erfc<std::complex<double>>();
283 test_erfc<std::complex<long double>>();
284 //test_erfc<boost::multiprecision::number<boost::multiprecision::mpc_complex_backend<20>>>();
285 //test_erfc<boost::multiprecision::number<boost::multiprecision::mpc_complex_backend<30>>>();
286 //test_erfc<boost::multiprecision::mpc_complex_50>();
287 //test_erfc<boost::multiprecision::mpc_complex_100>();
288
289 #ifdef BOOST_HAS_FLOAT128
290 test_complex_bessel<boost::multiprecision::complex128>();
291 test_I0_complex<boost::multiprecision::complex128>();
292 test_erfc<boost::multiprecision::complex128>();
293 #endif
294
295 }
296