...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> 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

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 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_{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.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_{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 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_{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.

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.