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)
7 #define BOOST_TEST_MODULE adaptive_gauss_kronrod_quadrature_test
9 #include <boost/config.hpp>
10 #include <boost/detail/workaround.hpp>
12 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
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>
21 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
29 #pragma warning(disable:4127) // Conditional expression is constant
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
;
67 Real
get_termination_condition()
69 return boost::math::tools::epsilon
<Real
>() * 1000;
73 template<class Real
, unsigned Points
>
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;
79 auto f
= [](const Real
& x
)
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));
91 template<class Real
, unsigned Points
>
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;
98 auto f
= [](const Real
& x
) { return 5*x
*x
+ 7*x
+ 12; };
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
>())));
107 // Examples taken from
108 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
109 template<class Real
, unsigned Points
>
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;
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
));
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
));
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
);
139 template<class Real
, unsigned Points
>
140 void test_three_quadrature_schemes_examples()
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;
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
);
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
);
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
);
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
);
174 template<class Real
, unsigned Points
>
175 void test_integration_over_real_line()
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;
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
>());
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
>());
200 template<class Real
, unsigned Points
>
201 void test_right_limit_infinite()
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;
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
>());
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
>());
224 template<class Real
, unsigned Points
>
225 void test_left_limit_infinite()
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;
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
);
239 BOOST_AUTO_TEST_CASE(gauss_quadrature_test
)
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>();
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>();
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>();
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>();
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>();
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>();
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>();
313 int main() { return 0; }