boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Add namespace boost::math::quadrature. Remove throw when condition number of summation exceeds precision of type; how to properly mollify the condition number is not clear and should be done consistently rather than ad-hoc.

Browse Source
This commit is contained in:
Nick Thompson
2017-05-11 21:03:00 -06:00
parent 48a0cf714b
commit 4801e2d8bf
4 changed files with 60 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
doc/quadrature/trapezoidal.qbk
View File

@@ -12,7 +12,7 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
``
#include <boost/math/quadrature/trapezoidal.hpp>
namespace boost{ namespace math{
namespace boost{ namespace math{ namespace quadrature {
template<class F, class Real>
Real trapezoidal(F f, Real a, Real b,
@@ -20,18 +20,21 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
size_t max_refinements = 10,
Real* error_estimate = nullptr,
Real* L1 = nullptr);
}} // namespaces
}}} // namespaces
``
[heading Description]
The functional `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
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[super 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.
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 faster than any power of the step size /h/.
If the integrand is periodic, and we integrate it over a period,
then the trapezoidal rule converges faster than any power of the step size /h/.
This can be seen by examination of the [@https://en.wikipedia.org/wiki/Euler-Maclaurin_formula Euler-Maclaurin summation 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.
@@ -43,8 +46,9 @@ For details, see [@http://epubs.siam.org/doi/pdf/10.1137/130932132 Trefethen's]
In its simplest form, an integration can be performed by the following code
using boost::math::quadrature::trapezoidal;
auto f = [](double x) { return 1/(5 - 4*cos(x)); };
double I = boost::math::trapezoidal(f, 0, boost::math::constants::two_pi<double>());
double I = trapezoidal(f, 0, boost::math::constants::two_pi<double>());
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
@@ -60,16 +64,19 @@ If the integrand is *not* periodic, then reducing the error to [radic][epsilon]
A question arises as to what to do when successive estimates never pass below the tolerance 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 10, which should never be hit in double precision unless the function is not periodic.
However, for higher-precision types, it may be of interest to allow the algorithm to compute more refinements:
By default, this maximum number of refinement steps is 10,
which should never be hit in double precision unless the function is not periodic.
However, for higher-precision types,
it may be of interest to allow the algorithm to compute more refinements:
size_t max_refinements = 15;
long double I = trapezoidal(f, 0, two_pi<long double>(), 1e-6L, max_refinements);
long double I = trapezoidal(f, 0, two_pi<long double>(), 1e-9L, max_refinements);
Note that the maximum allowed compute time grows exponentially with `max_refinements`.
The routine will not throw an exception if the maximum refinements is achieved without the requested tolerance being attained.
This is because the value calculated is more often than not still usable.
However, for applications with high-reliability requirements, the error estimate should be queried.
However, for applications with high-reliability requirements,
the error estimate should be queried.
This is achieved by passing additional pointers into the routine:
double error_estimate;
@@ -89,10 +96,12 @@ The condition number of summation is defined by
[kappa](S[sub n]) := [Sigma][sub i][super n] |x[sub i]|/|[Sigma][sub i][super n] x[sub i]|
If this number of ~10[super k], then /k/ additional digits are expected to be lost in addition to digits lost due to floating point rounding error.
As all numerical quadrature methods reduce to summation, their stability is controlled by the ratio \u222B |f| dt/|\u222B f dt |,
If this number of ~10[super k],
then /k/ additional digits are expected to be lost in addition to digits lost due to floating point rounding error.
As all numerical quadrature methods reduce to summation,
their stability is controlled by the ratio \u222B |f| dt/|\u222B f dt |,
which is easily seen to be equivalent to condition number of summation when evaluated numerically.
Hence both the error estimate and the condition number of summation should be analyzed in applications requiring very high precision.
Hence both the error estimate and the condition number of summation should be analyzed in applications requiring very high precision and reliability.
As an example, we consider evaluation of Bessel functions by trapezoidal quadrature.
The Bessel function of the first kind is defined via
@@ -101,10 +110,10 @@ J[sub n](x) = 1/2\u03A0 \u222B[sub -\u03A0][super \u03A0] cos(n t - x sin(t)) dt
The integrand is periodic, so the Euler-Maclaurin summation formula guarantees exponential convergence via the trapezoidal quadrature.
Without careful consideration, it seems this would be a very attractive method to compute Bessel functions.
However, we see that for large /n/ the integrand oscillates rapidly, taking on positive and negative values, and hence the trapezoidal sums become ill-conditioned.
However, we see that for large /n/ the integrand oscillates rapidly,
taking on positive and negative values,
and hence the trapezoidal sums become ill-conditioned.
In double precision, /x = 17/ and /n = 25/ gives a sum which is so poorly conditioned that zero correct digits are obtained.
In this case, the routine throw a `boost::math::evaluation_error` which explains the poor conditioning.
However, the error is only thrown when *all* digits are expected to be incorrect, so the user is still responsible for error analysis in many circumstances.
References:
@@ -116,9 +125,3 @@ Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Society f
[endsect]
[/
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)
]

8
example/adaptive_trapezoidal_example.cpp → example/trapezoidal_example.cpp
View File

@@ -1,5 +1,8 @@
/*
* 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)
*
* This example shows to to numerically integrate a periodic function using the adaptive_trapezoidal routine provided by boost.
*/
@@ -7,16 +10,17 @@
#include <iostream>
#include <cmath>
#include <limits>
#include <boost/math/quadrature/adaptive_trapezoidal.hpp>
#include <boost/math/quadrature/trapezoidal.hpp>
int main()
{
using boost::math::constants::two_pi;
using boost::math::constants::third;
using boost::math::quadrature::trapezoidal;
// 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<double>());
double Q = trapezoidal(f, (double) 0, two_pi<double>());
std::cout << std::setprecision(std::numeric_limits<double>::digits10);
std::cout << "The adaptive trapezoidal rule gives the integral of our function as " << Q << "\n";

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

@@ -12,16 +12,15 @@
* 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
#ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
#define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
#include <cmath>
#include <limits>
#include <stdexcept>
#include <boost/format.hpp>
#include <boost/math/constants/constants.hpp>
namespace boost{ namespace math{
namespace boost{ namespace math{ namespace quadrature {
template<class F, class Real>
Real trapezoidal(F f, Real a, Real b, Real tol = sqrt(std::numeric_limits<Real>::epsilon()), size_t max_refinements = 10, Real* error_estimate = nullptr, Real* L1 = nullptr)
@@ -94,18 +93,8 @@ Real trapezoidal(F f, Real a, Real b, Real tol = sqrt(std::numeric_limits<Real>:
*L1 = IL1;
}
// The condition abs(I) > eps is a poor man's way of mollifying the condition number.
// Recall that no number is close to zero, so it makes no sense to throw when the I1 is very close to zero.
if (abs(I1) > std::numeric_limits<Real>::epsilon() && IL1/abs(I1) > 1.0/std::numeric_limits<Real>::epsilon())
{
Real cond = IL1/abs(I1);
Real inv_prec = 1.0/std::numeric_limits<Real>::epsilon();
boost::basic_format<char> err_msg = boost::format("\nThe condition number of the quadrature sum is %1%, which exceeds the inverse of the precision of the type %2%.\nNo correct digits can be expected for the integral.\n") % cond % inv_prec;
throw boost::math::evaluation_error(err_msg.str());
}
return I1;
}
}}
}}}
#endif

42
test/test_trapezoidal.cpp
View File

@@ -4,7 +4,7 @@
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#define BOOST_TEST_MODULE trapezoidal
#define BOOST_TEST_MODULE trapezoidal_quadrature
#include <boost/type_index.hpp>
#include <boost/test/included/unit_test.hpp>
@@ -16,6 +16,7 @@
using boost::multiprecision::cpp_bin_float_50;
using boost::multiprecision::cpp_bin_float_100;
using boost::math::quadrature::trapezoidal;
template<class Real>
void test_constant()
@@ -23,7 +24,7 @@ void test_constant()
std::cout << "Testing constants are integrated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
auto f = [](Real x) { return boost::math::constants::half<Real>(); };
Real Q = boost::math::trapezoidal<decltype(f), Real>(f, (Real) 0.0, (Real) 10.0);
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());
}
@@ -38,9 +39,9 @@ void test_rational_periodic()
auto f = [](Real x) { return 1/(5 - 4*cos(x)); };
Real tol = 100*std::numeric_limits<Real>::epsilon();
Real Q = boost::math::trapezoidal(f, (Real) 0.0, two_pi<Real>(), tol);
Real Q = trapezoidal(f, (Real) 0.0, two_pi<Real>(), tol);
BOOST_CHECK_CLOSE(Q, two_pi<Real>()*third<Real>(), 200*tol);
BOOST_CHECK_CLOSE_FRACTION(Q, two_pi<Real>()*third<Real>(), 10*tol);
}
template<class Real>
@@ -54,10 +55,11 @@ void test_bump_function()
}
return (Real) exp(-(Real) 1/(1-x*x));
};
Real tol = 100*std::numeric_limits<Real>::epsilon();
Real Q = boost::math::trapezoidal(f, (Real) -1, (Real) 1, tol);
// This is all the digits Mathematica gave me!
BOOST_CHECK_CLOSE(Q, 0.443994, 1e-4);
Real tol = std::numeric_limits<Real>::epsilon();
Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
// 2*NIntegrate[Exp[-(1/(1 - x^2))], {x, 0, 1}, WorkingPrecision -> 210]
Real Q_exp = boost::lexical_cast<Real>("0.44399381616807943782304892117055266376120178904569749730748455394704");
BOOST_CHECK_CLOSE_FRACTION(Q, Q_exp, 15*tol);
}
template<class Real>
@@ -66,7 +68,7 @@ void test_zero_function()
std::cout << "Testing that zero functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
auto f = [](Real x) { return (Real) 0;};
Real tol = 100*std::numeric_limits<Real>::epsilon();
Real Q = boost::math::trapezoidal(f, (Real) -1, (Real) 1, tol);
Real Q = trapezoidal(f, (Real) -1, (Real) 1, tol);
BOOST_CHECK_SMALL(Q, 100*tol);
}
@@ -76,8 +78,8 @@ void test_sinsq()
std::cout << "Testing that sin(x)^2 is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
auto f = [](Real x) { return sin(10*x)*sin(10*x); };
Real tol = 100*std::numeric_limits<Real>::epsilon();
Real Q = boost::math::trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi<Real>(), tol);
BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi<Real>(), 100*tol);
Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi<Real>(), tol);
BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>(), tol);
}
@@ -90,26 +92,28 @@ void test_slowly_converging()
Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
Real error_estimate;
Real Q = boost::math::trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate);
BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi<Real>()/2, 1000*tol);
Real Q = trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate);
BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::half_pi<Real>()/2, 10*tol);
}
template<class Real>
void test_rational_sin()
{
using std::pow;
using std::sin;
using boost::math::constants::two_pi;
using boost::math::constants::half;
std::cout << "Testing that a rational sin function is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
Real a = 5;
auto f = [=](Real x) { Real t = a + sin(x); return 1.0/(t*t); };
Real expected = two_pi<Real>()*a/pow(a*a - 1, 1.5);
Real expected = two_pi<Real>()*a/pow(a*a - 1, 3*half<Real>());
Real tol = 100*std::numeric_limits<Real>::epsilon();
Real Q = boost::math::trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi<Real>(), tol);
BOOST_CHECK_CLOSE(Q, expected, 100*tol);
Real Q = trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi<Real>(), tol);
BOOST_CHECK_CLOSE_FRACTION(Q, expected, tol);
}
BOOST_AUTO_TEST_CASE(trapezoidal)
BOOST_AUTO_TEST_CASE(trapezoidal_quadrature)
{
test_constant<float>();
test_constant<double>();
@@ -123,7 +127,10 @@ BOOST_AUTO_TEST_CASE(trapezoidal)
test_rational_periodic<cpp_bin_float_50>();
test_rational_periodic<cpp_bin_float_100>();
test_bump_function<float>();
test_bump_function<double>();
test_bump_function<long double>();
test_rational_periodic<cpp_bin_float_50>();
test_zero_function<float>();
test_zero_function<double>();
@@ -144,5 +151,6 @@ BOOST_AUTO_TEST_CASE(trapezoidal)
test_rational_sin<float>();
test_rational_sin<double>();
test_rational_sin<long double>();
test_rational_sin<cpp_bin_float_50>();
}