boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-24 06:02:08 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #326 from boostorg/quadrature_bounds

Browse Source 
Accept b<= a in quadrature.
This commit is contained in:
Nick
2020-03-22 08:11:25 -04:00
committed by GitHub
parent abbe41b3b3 11c45c0ec0
commit 5800a6493f
8 changed files with 48 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
include/boost/math/quadrature/gauss.hpp
View File

@@ -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>();
Real scale = (b - a)*constants::half<Real>();

8
include/boost/math/quadrature/gauss_kronrod.hpp
View File

@@ -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<F> 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);
}
}

8
include/boost/math/quadrature/tanh_sinh.hpp
View File

@@ -145,9 +145,13 @@ auto tanh_sinh<Real, Policy>::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>();
Real diff = (b - a)*half<Real>();

13
include/boost/math/quadrature/trapezoidal.hpp
View File

@@ -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<K>(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<K>(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<K>(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<K>(0);
}
if(a > b)
{
return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol);
}
K ya = f(a);
K yb = f(b);

7
test/gauss_kronrod_quadrature_test.cpp
View File

@@ -264,6 +264,13 @@ void test_linear()
Real Q = gauss_kronrod<Real, Points>::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<Real, Points>::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<Real, Points>::integrate(f, (Real) 0, (Real) 0, 0, 0, &error, &L1);
BOOST_CHECK_CLOSE(Q, Real(0), tol);
}
template<class Real, unsigned Points>

4
test/gauss_quadrature_test.cpp
View File

@@ -259,6 +259,10 @@ void test_linear()
Real Q = gauss<Real, Points>::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<Real, Points>::integrate(f, (Real) 0, (Real) 0, &L1);
BOOST_CHECK_CLOSE(Q, 0, tol);
Q = gauss<Real, Points>::integrate(f, (Real) 1, (Real) 0, &L1);
BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
}
template<class Real, unsigned Points>

5
test/tanh_sinh_quadrature_test.cpp
View File

@@ -218,6 +218,11 @@ void test_linear()
Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &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<Real>(), &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<Real>(), &error, &L1);
BOOST_CHECK_EQUAL(Q, Real(0));
}

5
test/test_trapezoidal.cpp
View File

@@ -130,6 +130,11 @@ void test_constant()
auto f = [](Real)->Real { return boost::math::constants::half<Real>(); };
Real Q = trapezoidal<decltype(f), Real>(f, (Real) 0.0, (Real) 10.0);
BOOST_CHECK_CLOSE(Q, 5.0, 100*std::numeric_limits<Real>::epsilon());
Q = trapezoidal<decltype(f), Real>(f, (Real) 10.0, (Real) 0.0);
BOOST_CHECK_CLOSE(Q, -5.0, 100*std::numeric_limits<Real>::epsilon());
Q = trapezoidal<decltype(f), Real>(f, (Real) 10.0, (Real) 10.0);
BOOST_CHECK_CLOSE(Q, Real(0), 100*std::numeric_limits<Real>::epsilon());
}