/* * Copyright John Maddock, 2017 * Use, modification and distribution are subject to the * Boost Software License, Version 1.0. (See accompanying file * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) * * This example Illustrates numerical integration via Gauss and Gauss-Kronrod quadrature. */ #include #include #include #include #include #include #include #include void gauss_examples() { //[gauss_example /*` We'll begin by integrating t[super 2] atan(t) over (0,1) using a 7 term Gauss-Legendre rule, and begin by defining the function to integrate as a C++ lambda expression: */ using namespace boost::math::quadrature; auto f = [](const double& t) { return t * t * std::atan(t); }; /*` Integration is simply a matter of calling the `gauss::integrate` method: */ double Q = gauss::integrate(f, 0, 1); /*` Which yields a value 0.2106572512 accurate to 1e-10. For more accurate evaluations, we'll move to a multiprecision type and use a 20-point integration scheme: */ using boost::multiprecision::cpp_bin_float_quad; auto f2 = [](const cpp_bin_float_quad& t) { return t * t * atan(t); }; cpp_bin_float_quad Q2 = gauss::integrate(f2, 0, 1); /*` Which yields 0.2106572512258069881080923020669, which is accurate to 5e-28. */ //] std::cout << std::setprecision(18) << Q << std::endl; std::cout << boost::math::relative_difference(Q, (boost::math::constants::pi() - 2 + 2 * boost::math::constants::ln_two()) / 12) << std::endl; std::cout << std::setprecision(34) << Q2 << std::endl; std::cout << boost::math::relative_difference(Q2, (boost::math::constants::pi() - 2 + 2 * boost::math::constants::ln_two()) / 12) << std::endl; } void gauss_kronrod_examples() { //[gauss_kronrod_example /*` We'll begin by integrating exp(-t[super 2]/2) over (0,+[infin]) using a 7 term Gauss rule and 15 term Kronrod rule, and begin by defining the function to integrate as a C++ lambda expression: */ using namespace boost::math::quadrature; auto f1 = [](double t) { return std::exp(-t*t / 2); }; //<- double Q_expected = sqrt(boost::math::constants::half_pi()); //-> /*` We'll start off with a one shot (ie non-adaptive) integration, and keep track of the estimated error: */ double error; double Q = gauss_kronrod::integrate(f1, 0, std::numeric_limits::infinity(), 0, 0, &error); /*` This yields Q = 1.25348207361, which has an absolute error of 1e-4 compared to the estimated error of 5e-3: this is fairly typical, with the difference between Gauss and Gauss-Kronrod schemes being much higher than the actual error. Before moving on to adaptive quadrature, lets try again with more points, in fact with the largest Gauss-Kronrod scheme we have cached (30/61): */ //<- std::cout << std::setprecision(16) << Q << std::endl; std::cout << boost::math::relative_difference(Q, Q_expected) << std::endl; std::cout << fabs(Q - Q_expected) << std::endl; std::cout << error << std::endl; //-> Q = gauss_kronrod::integrate(f1, 0, std::numeric_limits::infinity(), 0, 0, &error); //<- std::cout << std::setprecision(16) << Q << std::endl; std::cout << boost::math::relative_difference(Q, Q_expected) << std::endl; std::cout << fabs(Q - Q_expected) << std::endl; std::cout << error << std::endl; //-> /*` This yields an absolute error of 3e-15 against an estimate of 1e-8, which is about as good as we're going to get at double precision However, instead of continuing with ever more points, lets switch to adaptive integration, and set the desired relative error to 1e-14 against a maximum depth of 5: */ Q = gauss_kronrod::integrate(f1, 0, std::numeric_limits::infinity(), 5, 1e-14, &error); //<- std::cout << std::setprecision(16) << Q << std::endl; std::cout << boost::math::relative_difference(Q, Q_expected) << std::endl; std::cout << fabs(Q - Q_expected) << std::endl; std::cout << error << std::endl; //-> /*` This yields an actual error of zero, against an estimate of 4e-15. In fact in this case the requested tolerance was almost certainly set too low: as we've seen above, for smooth functions, the precision achieved is often double that of the estimate, so if we integrate with a tolerance of 1e-9: */ Q = gauss_kronrod::integrate(f1, 0, std::numeric_limits::infinity(), 5, 1e-9, &error); //<- std::cout << std::setprecision(16) << Q << std::endl; std::cout << boost::math::relative_difference(Q, Q_expected) << std::endl; std::cout << fabs(Q - Q_expected) << std::endl; std::cout << error << std::endl; //-> /*` We still achieve 1e-15 precision, with an error estimate of 1e-10. */ //] } int main() { gauss_examples(); gauss_kronrod_examples(); return 0; }