From 11c45c0ec0d29bf722c11c6237b6c40b9e353d67 Mon Sep 17 00:00:00 2001 From: Nick Date: Sat, 21 Mar 2020 09:22:41 -0400 Subject: [PATCH] Accept b<= a in quadrature. --- include/boost/math/quadrature/gauss.hpp | 8 ++++++-- include/boost/math/quadrature/gauss_kronrod.hpp | 8 ++++++-- include/boost/math/quadrature/tanh_sinh.hpp | 8 ++++++-- include/boost/math/quadrature/trapezoidal.hpp | 13 +++++++++---- test/gauss_kronrod_quadrature_test.cpp | 7 +++++++ test/gauss_quadrature_test.cpp | 4 ++++ test/tanh_sinh_quadrature_test.cpp | 5 +++++ test/test_trapezoidal.cpp | 5 +++++ 8 files changed, 48 insertions(+), 10 deletions(-) diff --git a/include/boost/math/quadrature/gauss.hpp b/include/boost/math/quadrature/gauss.hpp index 2f67666f3..b8caee11a 100644 --- a/include/boost/math/quadrature/gauss.hpp +++ b/include/boost/math/quadrature/gauss.hpp @@ -1258,9 +1258,13 @@ public: if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) { - if (b <= a) + if (a == b) { - return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy()); + return K(0); + } + if (b < a) + { + return -integrate(f, b, a, pL1); } Real avg = (a + b)*constants::half(); Real scale = (b - a)*constants::half(); diff --git a/include/boost/math/quadrature/gauss_kronrod.hpp b/include/boost/math/quadrature/gauss_kronrod.hpp index fe6ba92c9..3deeb8bef 100644 --- a/include/boost/math/quadrature/gauss_kronrod.hpp +++ b/include/boost/math/quadrature/gauss_kronrod.hpp @@ -1927,11 +1927,15 @@ public: if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) { - if (b <= a) + if (a==b) { - return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy()); + return K(0); } recursive_info info = { f, tol }; + if (b < a) + { + return -recursive_adaptive_integrate(&info, b, a, max_depth, Real(0), error, pL1); + } return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1); } } diff --git a/include/boost/math/quadrature/tanh_sinh.hpp b/include/boost/math/quadrature/tanh_sinh.hpp index 25c65111d..056ecfa88 100644 --- a/include/boost/math/quadrature/tanh_sinh.hpp +++ b/include/boost/math/quadrature/tanh_sinh.hpp @@ -145,9 +145,13 @@ auto tanh_sinh::integrate(const F f, Real a, Real b, Real toleranc if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) { - if (b <= a) + if (a == b) { - return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy()); + return result_type(0); + } + if (b < a) + { + return -this->integrate(f, b, a, tolerance, error, L1, levels); } Real avg = (a + b)*half(); Real diff = (b - a)*half(); diff --git a/include/boost/math/quadrature/trapezoidal.hpp b/include/boost/math/quadrature/trapezoidal.hpp index 4899e356c..935b01a6b 100644 --- a/include/boost/math/quadrature/trapezoidal.hpp +++ b/include/boost/math/quadrature/trapezoidal.hpp @@ -33,10 +33,6 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea // In many math texts, K represents the field of real or complex numbers. // Too bad we can't put blackboard bold into C++ source! typedef decltype(f(a)) K; - if(a >= b) - { - return static_cast(boost::math::policies::raise_domain_error(function, "a < b for integration over the region [a, b] is required, but got a = %1%.\n", a, pol)); - } if (!(boost::math::isfinite)(a)) { return static_cast(boost::math::policies::raise_domain_error(function, "Left endpoint of integration must be finite for adaptive trapezoidal integration but got a = %1%.\n", a, pol)); @@ -46,6 +42,15 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea return static_cast(boost::math::policies::raise_domain_error(function, "Right endpoint of integration must be finite for adaptive trapezoidal integration but got b = %1%.\n", b, pol)); } + if (a == b) + { + return static_cast(0); + } + if(a > b) + { + return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol); + } + K ya = f(a); K yb = f(b); diff --git a/test/gauss_kronrod_quadrature_test.cpp b/test/gauss_kronrod_quadrature_test.cpp index 3518ffcaa..c80ffe854 100644 --- a/test/gauss_kronrod_quadrature_test.cpp +++ b/test/gauss_kronrod_quadrature_test.cpp @@ -264,6 +264,13 @@ void test_linear() Real Q = gauss_kronrod::integrate(f, (Real) 0, (Real) 1, 0, 0, &error, &L1); BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol); BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); + + Q = gauss_kronrod::integrate(f, (Real) 1, (Real) 0, 0, 0, &error, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); + + Q = gauss_kronrod::integrate(f, (Real) 0, (Real) 0, 0, 0, &error, &L1); + BOOST_CHECK_CLOSE(Q, Real(0), tol); } template diff --git a/test/gauss_quadrature_test.cpp b/test/gauss_quadrature_test.cpp index 61145db16..b8b7f3b72 100644 --- a/test/gauss_quadrature_test.cpp +++ b/test/gauss_quadrature_test.cpp @@ -259,6 +259,10 @@ void test_linear() Real Q = gauss::integrate(f, (Real) 0, (Real) 1, &L1); BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol); BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); + Q = gauss::integrate(f, (Real) 0, (Real) 0, &L1); + BOOST_CHECK_CLOSE(Q, 0, tol); + Q = gauss::integrate(f, (Real) 1, (Real) 0, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol); } template diff --git a/test/tanh_sinh_quadrature_test.cpp b/test/tanh_sinh_quadrature_test.cpp index e1473a872..6868ee933 100644 --- a/test/tanh_sinh_quadrature_test.cpp +++ b/test/tanh_sinh_quadrature_test.cpp @@ -218,6 +218,11 @@ void test_linear() Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance(), &error, &L1); BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol); BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); + Q = integrator.integrate(f, (Real) 1, (Real) 0, get_convergence_tolerance(), &error, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); + Q = integrator.integrate(f, (Real) 1, (Real) 1, get_convergence_tolerance(), &error, &L1); + BOOST_CHECK_EQUAL(Q, Real(0)); } diff --git a/test/test_trapezoidal.cpp b/test/test_trapezoidal.cpp index 6f78c8b83..c99d9208e 100644 --- a/test/test_trapezoidal.cpp +++ b/test/test_trapezoidal.cpp @@ -130,6 +130,11 @@ void test_constant() auto f = [](Real)->Real { return boost::math::constants::half(); }; Real Q = trapezoidal(f, (Real) 0.0, (Real) 10.0); BOOST_CHECK_CLOSE(Q, 5.0, 100*std::numeric_limits::epsilon()); + Q = trapezoidal(f, (Real) 10.0, (Real) 0.0); + BOOST_CHECK_CLOSE(Q, -5.0, 100*std::numeric_limits::epsilon()); + + Q = trapezoidal(f, (Real) 10.0, (Real) 10.0); + BOOST_CHECK_CLOSE(Q, Real(0), 100*std::numeric_limits::epsilon()); }