...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

#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 𝑶(h^{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.

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 C_{0}^{∞}(∝) 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 L_{1} 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 L_{1} norm of the integral. This is computed along
with the integral, and is an essential component of the algorithm. First,
the L_{1} norm establishes a scale against which the error can be measured. Second,
the L_{1} 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

κ(S_{n}) := Σ_{i}^{n} |x_{i}|/|Σ_{i}^{n} x_{i}|

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

J_{n}(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.