Synopsis

#include <boost/math/quadrature/trapezoidal.hpp>
namespace boost{ namespace math{ namespace quadrature {

template<class F, class Real>
auto trapezoidal(F f, Real a, Real b,
                 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
                 size_t max_refinements = 12,
                 Real* error_estimate = nullptr,
                 Real* L1 = nullptr);

template<class F, class Real, class Policy>
auto trapezoidal(F f, Real a, Real b, Real tol, size_t max_refinements,
                 Real* error_estimate, Real* L1, const Policy& pol);

}}} // namespaces

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 is 𝑶(h²). 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 and the Bernoulli numbers. 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 C₀^∞(∝) 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.0, boost::math::constants::two_pi<double>());

The integrand must accept a real number argument, but can return a complex number. This is useful for contour integrals (which are manifestly periodic) and high-order numerical differentiation of analytic functions. An example using the integral definition of the complex Bessel function is shown here:

auto bessel_integrand = [&n, &z](double theta)->std::complex<double>
{
    std::complex<double> z{2, 3};
    using std::cos;
    using std::sin;
    return cos(z*sin(theta) - 2*theta)/pi<double>();
};

using boost::math::quadrature::trapezoidal;
std::complex<double> Jnz = trapezoidal(bessel_integrand, 0.0, pi<Real>());

Other special functions which are efficiently evaluated in the complex plane by trapezoidal quadrature are modified Bessel functions and the complementary error function. Another application of complex-valued trapezoidal quadrature is computation of high-order numerical derivatives; see Lyness and Moler for details.

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 L₁ 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 repeatedly, generating an exponential explosion in function evaluations. 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 12, 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.0L, 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 = trapezoidal(f, 0.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 L₁ norm of the integral. This is computed along with the integral, and is an essential component of the algorithm. First, the L₁ norm establishes a scale against which the error can be measured. Second, the L₁ 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

If this number of ~10^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 ∫ |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

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.

Lyness, James N., and Cleve B. Moler. Numerical differentiation of analytic functions. SIAM Journal on Numerical Analysis 4.2 (1967): 202-210.

Gil, Amparo, Javier Segura, and Nico M. Temme. Computing special functions by using quadrature rules. Numerical Algorithms 33.1-4 (2003): 265-275.