|
Home | Libraries | People | FAQ | More |
#include <boost/math/quadrature/trapezoidal.hpp> 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); template<class F, class Real, class Policy> Real trapezoidal(F f, Real a, Real b, Real tol, size_t max_refinements, Real* error_estimate, Real* L1, const Policy& pol); }}} // namespaces
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 is 𝑶(h2). 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. This can be seen by examination of the 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. Hence the trapezoidal rule is essentially optimal for periodic integrands.
Other classes of integrands which are integrated efficiently by this method are the C0∞(∝) bump functions and bell-shaped integrals over the infinite interval. For details, see Trefethen's SIAM review.
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 = 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
double I = trapezoidal(f, 0.0, two_pi<double>(), 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 good choice for the tolerance is √ε, which is the default. If the
integrand is periodic, then the number of correct digits should double on
each interval halving. Hence, once the integration routine has estimated
that the error is √ε, then the actual error should be ~ε. If the integrand is
not periodic, then reducing the error to
√ε takes much longer, but is nonetheless possible without becoming a major performance
bug.
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:
size_t max_refinements = 15; 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. This is achieved by passing additional
pointers into the routine:
double error_estimate; double L1; double I = adaptive_trapezoidal(f, 0, two_pi<double>(), tolerance, max_refinements, &error_estimate, &L1); if (error_estimate > tolerance*L1) { double I = some_other_quadrature_method(f, 0, two_pi<double>()); }
The final argument is the L1 norm of the integral. This is computed along with the integral, and is an essential component of the algorithm. First, the L1 norm establishes a scale against which the error can be measured. Second, the L1 norm can be used to evaluate the stability of the computation. This can be formulated in a rigorous manner by defining the condition number of summation. The condition number of summation is defined by
κ(Sn) := Σin |xi|/|Σin xi|
If this number of ~10k, 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 ∫ |f| dt/|∫ 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 and reliability.
As an example, we consider evaluation of Bessel functions by trapezoidal quadrature. The Bessel function of the first kind is defined via
Jn(x) = 1/2Π ∫-ΠΠ 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. In double precision, x = 17 and n = 25 gives a sum which is so poorly conditioned that zero correct digits are obtained.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
References:
Trefethen, Lloyd N., Weideman, J.A.C., The Exponentially Convergent Trapezoidal Rule, SIAM Review, Vol. 56, No. 3, 2014.
Stoer, Josef, and Roland Bulirsch. Introduction to numerical analysis. Vol. 12., Springer Science & Business Media, 2013.
Higham, Nicholas J. Accuracy and stability of numerical algorithms. Society for industrial and applied mathematics, 2002.