The integration range for each dimension may be specified
using a list or tuple. Arguments are interpreted as follows:

quad(f,[x1,x2]) – calculates

quad(f,[x1,x2],[y1,y2]) – calculates

quad(f,[x1,x2],[y1,y2],[z1,z2]) – calculates

Endpoints may be finite or infinite. An interval descriptor
may also contain more than two points. In this
case, the integration is split into subintervals, between
each pair of consecutive points. This is useful for
dealing with mid-interval discontinuities, or integrating
over large intervals where the function is irregular or
oscillates.

If set to true, quad() returns where is the
integral and is the estimated error.

maxdegree

Maximum degree of the quadrature rule to try before
quitting.

verbose

Print details about progress.

Algorithms

Mpmath presently implements two integration algorithms: tanh-sinh
quadrature and Gauss-Legendre quadrature. These can be selected
using method=’tanh-sinh’ or method=’gauss-legendre’ or by
passing the classes method=TanhSinh, method=GaussLegendre.
The functions quadts() and quadgl() are also available
as shortcuts.

Both algorithms have the property that doubling the number of
evaluation points roughly doubles the accuracy, so both are ideal
for high precision quadrature (hundreds or thousands of digits).

At high precision, computing the nodes and weights for the
integration can be expensive (more expensive than computing the
function values). To make repeated integrations fast, nodes
are automatically cached.

The advantages of the tanh-sinh algorithm are that it tends to
handle endpoint singularities well, and that the nodes are cheap
to compute on the first run. For these reasons, it is used by
quad() as the default algorithm.

Gauss-Legendre quadrature often requires fewer function
evaluations, and is therefore often faster for repeated use, but
the algorithm does not handle endpoint singularities as well and
the nodes are more expensive to compute. Gauss-Legendre quadrature
can be a better choice if the integrand is smooth and repeated
integrations are required (e.g. for multiple integrals).

See the documentation for TanhSinh and
GaussLegendre for additional details.

Examples of 1D integrals

Intervals may be infinite or half-infinite. The following two
examples evaluate the limits of the inverse tangent function
(), and the Gaussian integral
:

For nonrectangular areas, one can call quad() recursively.
For example, we can replicate the earlier example of calculating
by integrating over the unit-circle, and actually use double
quadrature to actually measure the area circle:

Both tanh-sinh and Gauss-Legendre quadrature are designed to
integrate smooth (infinitely differentiable) functions. Neither
algorithm copes well with mid-interval singularities (such as
mid-interval discontinuities in or ).
The best solution is to split the integral into parts:

For functions that are smooth (in the sense of being infinitely
differentiable) but contain sharp mid-interval peaks or many
“bumps”, quad() may fail to provide full accuracy. For
example, with default settings, quad() is able to integrate
accurately over an interval of length 100 but not over
length 1000:

where at least one of and is infinite and where
for some slowly
decreasing function . With proper input, quadosc()
can also handle oscillatory integrals where the oscillation
rate is different from a pure sine or cosine wave.

In the standard case when ,
quadosc() works by evaluating the infinite series

where are consecutive zeros (alternatively
some other periodic reference point) of .
Accordingly, quadosc() requires information about the
zeros of . For a periodic function, you can specify
the zeros by either providing the angular frequency
(omega) or the period. In general, you can
specify the -th zero by providing the zeros arguments.
Below is an example of each:

Note that zeros was specified to multiply by the
half-period, not the full period. In theory, it does not matter
whether each partial integral is done over a half period or a full
period. However, if done over half-periods, the infinite series
passed to nsum() becomes an alternating series and this
typically makes the extrapolation much more efficient.

Here is an example of an integration over the entire real line,
and a half-infinite integration starting at :

For an example where zeros becomes necessary, consider the
complete Fresnel integrals

Although the integrands do not decrease in magnitude as
, the integrals are convergent since the oscillation
rate increases (causing consecutive periods to asymptotically
cancel out). These integrals are virtually impossible to calculate
to any kind of accuracy using standard quadrature rules. However,
if one provides the correct asymptotic distribution of zeros
(), quadosc() works:

quadosc() is primarily useful for slowly decaying
integrands. If the integrand decreases exponentially or faster,
quad() will likely handle it without trouble (and generally be
much faster than quadosc()):

Given results from integrations done
with a quadrature of rule of degree , estimate
the error of .

For , we estimate as .

For , we extrapolate
from and under the assumption
that each degree increment roughly doubles the accuracy of
the quadrature rule (this is true for both TanhSinh
and GaussLegendre). The extrapolation formula is given
by Borwein, Bailey & Girgensohn. Although not very conservative,
this method seems to be very robust in practice.

Given a desired precision in bits, estimate the degree
of the quadrature required to accomplish full accuracy for
typical integrals. By default, quad() will perform up
to iterations. The value of should be a slight
overestimate, so that “slightly bad” integrals can be dealt
with automatically using a few extra iterations. On the
other hand, it should not be too big, so quad() can
quit within a reasonable amount of time when it is given
an “unsolvable” integral.

The default formula used by guess_degree() is tuned
for both TanhSinh and GaussLegendre.
The output is roughly as follows:

50

6

100

7

500

10

3000

12

This formula is based purely on a limited amount of
experimentation and will sometimes be wrong.

Main integration function. Computes the 1D integral over
the interval specified by points. For each subinterval,
performs quadrature of degree from 1 up to max_degree
until estimate_error() signals convergence.

summation() transforms each subintegration to
the standard interval and then calls sum_next().

This class implements “tanh-sinh” or “doubly exponential”
quadrature. This quadrature rule is based on the Euler-Maclaurin
integral formula. By performing a change of variables involving
nested exponentials / hyperbolic functions (hence the name), the
derivatives at the endpoints vanish rapidly. Since the error term
in the Euler-Maclaurin formula depends on the derivatives at the
endpoints, a simple step sum becomes extremely accurate. In
practice, this means that doubling the number of evaluation
points roughly doubles the number of accurate digits.

Comparison to Gauss-Legendre:

Initial computation of nodes is usually faster

Handles endpoint singularities better

Handles infinite integration intervals better

Is slower for smooth integrands once nodes have been computed

The implementation of the tanh-sinh algorithm is based on the
description given in Borwein, Bailey & Girgensohn, “Experimentation
in Mathematics - Computational Paths to Discovery”, A K Peters,
2003, pages 312-313. In the present implementation, a few
improvements have been made:

A more efficient scheme is used to compute nodes (exploiting
recurrence for the exponential function)

The nodes are computed successively instead of all at once

Various documents describing the algorithm are available online, e.g.:

Step sum for tanh-sinh quadrature of degree . We exploit the
fact that half of the abscissas at degree are precisely the
abscissas from degree . Thus reusing the result from
the previous level allows a 2x speedup.

This class implements Gauss-Legendre quadrature, which is
exceptionally efficient for polynomials and polynomial-like (i.e.
very smooth) integrands.

The abscissas and weights are given by roots and values of
Legendre polynomials, which are the orthogonal polynomials
on with respect to the unit weight
(see legendre()).

In this implementation, we take the “degree” of the quadrature
to denote a Gauss-Legendre rule of degree (following
Borwein, Bailey & Girgensohn). This way we get quadratic, rather
than linear, convergence as the degree is incremented.