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