From 136e7411f59745eecbb041ce6f7c4ccb83f2e182 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 5 Mar 2017 19:05:41 -0600 Subject: [PATCH] Adaptive Trapezoidal Quadrature This routine estimates the definite integral of a function f. Assuming that f is periodic, it can be shown that this routine converges exponentially fast. In fact, the test cases given exhibit exponential convergence with decreasing stepsize. A potential improvement is using the Bulirsch sequence rather than the Romberg sequence to schedule the refinements. However, the convergence is so rapid for functions of the class specified above that there seems to be no need at present. This code is cppcheck clean, and runs successfully under AddressSanitizer and UndefinedBehaviorSanitizer. --- doc/quadrature/adaptive_trapezoidal.qbk | 61 +++++++ example/adaptive_trapezoidal_example.cpp | 25 +++ .../math/quadrature/adaptive_trapezoidal.hpp | 90 ++++++++++ test/test_adaptive_trapezoidal.cpp | 154 ++++++++++++++++++ 4 files changed, 330 insertions(+) create mode 100644 doc/quadrature/adaptive_trapezoidal.qbk create mode 100644 example/adaptive_trapezoidal_example.cpp create mode 100644 include/boost/math/quadrature/adaptive_trapezoidal.hpp create mode 100644 test/test_adaptive_trapezoidal.cpp diff --git a/doc/quadrature/adaptive_trapezoidal.qbk b/doc/quadrature/adaptive_trapezoidal.qbk new file mode 100644 index 000000000..548de93f2 --- /dev/null +++ b/doc/quadrature/adaptive_trapezoidal.qbk @@ -0,0 +1,61 @@ +[/ +Copyright (c) 2017 Nick Thompson +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) +] + +[section:adaptive_trapezoidal Adaptive Trapezoidal Quadrature] + +[heading Synopsis] + +`` + #include +`` + +[heading Description] + +The function __adaptive_trapezoidal calculates the integral of a function /f/ using the surprisingly simple trapezoidal rule. +If we assume only that the integrand is twice continuously differentiable, we can prove that the error of the composite trapezoidal rule +is [bigo](h^2). +Hence halving the interval only cuts the error by about a fourth, which in turn implies that we must evaluate the function many times before an acceptable accuracy can be achieved. + +However, the trapezoidal rule has an astonishing property: If the integrand is periodic, and we integrate it over a period, then the trapezoidal rule converges exponentially fast. +This can be seen by examination of the Euler-Maclaurin formula, which relates a definite integral to its trapezoidal sum and error terms proportional to the derivatives of the function at the endpoints. +If the derivatives at the endpoints are the same or vanish, then the error very nearly vanishes. +Hence the trapezoidal rule is essentially optimal for periodic integrands. + +In it's simplest form, an integration can be performed by the following code + + __auto f = [](double x) { return 1/(5 - 4*cos(x)); }; + __double I = boost::math::adaptive_trapezoidal(f, 0, two_pi()); + +Since the routine is adaptive, step sizes are halved continuously until a tolerance is reached. +In order to control this tolerance, simply call the routine with an additional argument + + __double I = adaptive_trapezoidal(f, 0, two_pi(), 1e-6); + +The routine stops when successive estimates of the integral /I1/ and /I0/ differ by less than the tolerance multiplied by the estimated L1 norm of the function. +A question arises as to what to do when successive estimates never pass below this threshold. +The stepsize would be halved until it eventually would be flushed to zero, leading to an infinite loop. +As such, you may pass an optional argument __max_refinements which controls how many times the interval may be halved before giving up. +By default, this maximum number of refinement steps is 15, leading to ~30,000 function evaluations. +For certain problems, this is clearly excessive, so specifying a smaller number is reasonable + + __size_t max_refinements = 10; + __double I = adaptive_trapezoidal(f, 0, two_pi(), 1e-6, max_refinements); + +Finally, you might wonder what the final error estimate actually was: This is achieved by passing a pointer into the last argument + + __double error_estimate; + __double I = adaptive_trapezoidal(f, 0, two_pi(), 1e-6, 10, &error_estimate); + +If need be, you can query this value, and if it is unacceptably large, use another method. + + +References: + +Stoer, Josef, and Roland Bulirsch. Introduction to numerical analysis. Vol. 12. Springer Science & Business Media, 2013. + + +[endsect] diff --git a/example/adaptive_trapezoidal_example.cpp b/example/adaptive_trapezoidal_example.cpp new file mode 100644 index 000000000..af833dfd6 --- /dev/null +++ b/example/adaptive_trapezoidal_example.cpp @@ -0,0 +1,25 @@ +/* + * Copyright Nick Thompson, 2017 + * + * This example shows to to numerically integrate a periodic function using the adaptive_trapezoidal routine provided by boost. + */ + +#include +#include +#include +#include + +int main() +{ + using boost::math::constants::two_pi; + using boost::math::constants::third; + // This function has an analytic form for its integral over a period: 2pi/3. + auto f = [](double x) { return 1/(5 - 4*cos(x)); }; + + double Q = boost::math::adaptive_trapezoidal(f, (double) 0, two_pi()); + + std::cout << std::setprecision(std::numeric_limits::digits10); + std::cout << "The adaptive trapezoidal rule gives the integral of our function as " << Q << "\n"; + std::cout << "The exact result is " << two_pi()*third() << "\n"; + +} diff --git a/include/boost/math/quadrature/adaptive_trapezoidal.hpp b/include/boost/math/quadrature/adaptive_trapezoidal.hpp new file mode 100644 index 000000000..42bb41ffe --- /dev/null +++ b/include/boost/math/quadrature/adaptive_trapezoidal.hpp @@ -0,0 +1,90 @@ +/* + * Copyright Nick Thompson, 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) + * + * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period, + * or to integrate a function whose derivative vanishes at the endpoints. + * + * If your function does not satisfy these conditions, and instead is simply continuous and bounded + * over the whole interval, then this routine will still converge, albeit slowly. However, there + * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature. + */ + +#ifndef BOOST_MATH_QUADRATURE_ADAPTIVE_TRAPEZOIDAL_HPP +#define BOOST_MATH_QUADRATURE_ADAPTIVE_TRAPEZOIDAL_HPP + +#include +#include +#include + +namespace boost{ namespace math{ + +template +Real adaptive_trapezoidal(F f, Real a, Real b, Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 15, Real* error_estimate = nullptr) +{ + using std::abs; + using std::isfinite; + using boost::math::constants::half; + if(a >= b) + { + throw std::domain_error("a < b for integration over the region [a, b] is required.\n"); + } + if (!isfinite(a)) + { + throw std::domain_error("Left endpoint of integration must be finite for adaptive trapezoidal integration.\n"); + } + if (!isfinite(b)) + { + throw std::domain_error("Right endpoint of integration must be finite for adaptive trapedzoidal integration.\n"); + } + + Real ya = f(a); + Real yb = f(b); + Real h = (b - a)*half(); + Real I0 = (ya + yb)*h; + Real IL0 = (abs(ya) + abs(yb))*h; + + Real yh = f(a + h); + Real I1 = half()*I0 + yh*h; + Real IL1 = half()*IL0 + abs(yh)*h; + + // The recursion is: + // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k) + size_t k = 2; + // We want to go through at least 5 levels so we have sampled the function at least 33 times. + // Otherwise, we could terminate prematurely and miss essential features. + // This is of course possible anyway, but 33 samples seems to be a reasonable compromise. + Real error = abs(I0 - I1); + while (k < 5 || (k < max_refinements && error > tol*IL1) ) + { + I0 = I1; + IL0 = IL1; + + I1 = half()*I0; + IL1 = half()*IL0; + size_t p = 1 << k; + h *= half(); + Real sum = 0; + Real absum = 0; + for(size_t j = 1; j < p; j += 2) + { + Real y = f(a + j*h); + sum += y; + absum += abs(y); + } + I1 += sum*h; + IL1 += absum*h; + ++k; + error = abs(I0 - I1); + if(error_estimate) + { + *error_estimate = error; + } + } + return I1; +} + +}} +#endif diff --git a/test/test_adaptive_trapezoidal.cpp b/test/test_adaptive_trapezoidal.cpp new file mode 100644 index 000000000..93fd6cb6a --- /dev/null +++ b/test/test_adaptive_trapezoidal.cpp @@ -0,0 +1,154 @@ +#define BOOST_TEST_MODULE adaptive_trapezoidal + +#include +#include +#include +#include +#include +#include +#ifdef __GNUC__ +#ifndef __clang__ +#include +#endif +#endif + +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; + +template +void test_constant() +{ + std::cout << "Testing constants are integrated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + auto f = [](Real x) { return boost::math::constants::half(); }; + + Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0.0, (Real) 10.0); + + BOOST_CHECK_CLOSE(Q, 5.0, 100*std::numeric_limits::epsilon()); +} + + +template +void test_rational_periodic() +{ + using boost::math::constants::two_pi; + using boost::math::constants::third; + std::cout << "Testing that rational periodic functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + auto f = [](Real x) { return 1/(5 - 4*cos(x)); }; + + Real tol = 100*std::numeric_limits::epsilon(); + Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0.0, two_pi(), tol); + + BOOST_CHECK_CLOSE(Q, two_pi()*third(), 200*tol); +} + +template +void test_bump_function() +{ + std::cout << "Testing that bump functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id().pretty_name() << "\n"; + auto f = [](Real x) { + if( x>= 1 || x <= -1) + { + return (Real) 0; + } + return (Real) exp(-(Real) 1/(1-x*x)); + }; + Real tol = 100*std::numeric_limits::epsilon(); + Real Q = boost::math::adaptive_trapezoidal(f, (Real) -1, (Real) 1, tol); + // This is all the digits Mathematica gave me! + BOOST_CHECK_CLOSE(Q, 0.443994, 1e-4); +} + +template +void test_zero_function() +{ + std::cout << "Testing that zero functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id().pretty_name() << "\n"; + auto f = [](Real x) { return (Real) 0;}; + Real tol = 100*std::numeric_limits::epsilon(); + Real Q = boost::math::adaptive_trapezoidal(f, (Real) -1, (Real) 1, tol); + BOOST_CHECK_SMALL(Q, 100*tol); +} + +template +void test_sinsq() +{ + std::cout << "Testing that sin(x)^2 is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id().pretty_name() << "\n"; + auto f = [](Real x) { return sin(10*x)*sin(10*x); }; + Real tol = 100*std::numeric_limits::epsilon(); + Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi(), tol); + BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi(), 100*tol); + +} + +template +void test_slowly_converging() +{ + std::cout << "Testing that non-periodic functions are correctly integrated by the trapezoidal rule, even if slowly, on type " << boost::typeindex::type_id().pretty_name() << "\n"; + // This function is not periodic, so it should not be fast to converge: + auto f = [](Real x) { return sqrt(1 - x*x); }; + + Real tol = sqrt(sqrt(std::numeric_limits::epsilon())); + Real error_estimate; + Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate); + BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi()/2, 1000*tol); +} + +template +void test_rational_sin() +{ + using std::pow; + using boost::math::constants::two_pi; + std::cout << "Testing that a rational sin function is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real a = 5; + auto f = [=](Real x) { Real t = a + sin(x); return 1.0/(t*t); }; + + Real expected = two_pi()*a/pow(a*a - 1, 1.5); + Real tol = 100*std::numeric_limits::epsilon(); + Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi(), tol); + BOOST_CHECK_CLOSE(Q, expected, 100*tol); +} + +BOOST_AUTO_TEST_CASE(adaptive_trapezoidal) +{ + test_constant(); + test_constant(); + test_constant(); + test_constant(); + test_constant(); + + test_rational_periodic(); + test_rational_periodic(); + test_rational_periodic(); + test_rational_periodic(); + test_rational_periodic(); + + test_bump_function(); + + test_zero_function(); + test_zero_function(); + test_zero_function(); + test_zero_function(); + test_zero_function(); + + test_sinsq(); + test_sinsq(); + test_sinsq(); + test_sinsq(); + test_sinsq(); + + test_slowly_converging(); + test_slowly_converging(); + test_slowly_converging(); + + test_rational_sin(); + test_rational_sin(); + test_rational_sin(); + +#ifdef __GNUC__ +#ifndef __clang__ + test_constant(); + test_rational_periodic(); +#endif +#endif +}