Integration

Basic numerical integration (quad)

The function quad performs numerical integration (quadrature). The syntax for integrating a function f between the endpoints a and b is quad(f, [a, b]). For example:

>>> from mpmath import *
>>> mp.dps = 15
>>> print quadts(sin, [0, pi])
2.0

Integrals can typically be evaluated to arbitrary precision. For example, quad can easily calculate 50 digits of pi by integrating the area under the half circle arc x^2 + y^2 = 1 (y > 0):

>>> mp.dps = 50
>>> print quad(lambda x: 2*sqrt(1 - x**2), [-1, 1])
3.1415926535897932384626433832795028841971693993751

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'. 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, mpmath caches computed nodes.

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 well and the nodes are more expensive to compute. Gauss-Legendre quadrature might be a better choice if the integrand is smooth and repeated integrations are required.

Neat examples

Intervals can be infinite or half-infinite:

>>> mp.dps = 15
>>> print quad(lambda x: 2/(x**2+1), [0, inf])
3.14159265358979
>>> print quad(lambda x: exp(-x**2), [-inf, inf])**2
3.14159265358979

Complex integrals are also supported. The next example computes Euler’s constant gamma by integrating around the pole of the Riemann zeta function at z = 1:

>>> print 1/(2*pi)*quadts(lambda x: zeta(exp(j*x)+1), [0, 2*pi])
(0.577215664901533 - 1.04534339656676e-23j)

Functions with integral representations, such as the gamma function, can be implemented directly from the definition:

>>> def Gamma(z):
...     return quadts(lambda t: exp(-t)*t**(z-1), [0, inf])
...
>>> print Gamma(1)
1.0
>>> print Gamma(10)
362880.0
>>> print Gamma(1+1j)
(0.498015668118356 - 0.154949828301811j)

Multiple integrals

It is possible to calculate double integrals with quad. To do this, simply provide a two-argument function and two intervals. The first interval specifies the range for the x variable and the second interval specifies the range of the y variable:

>>> f = lambda x, y: cos(x+y/2)
>>> print quadts(f, [-pi/2, pi/2], [0, pi])
4.0

Here are some more difficult examples taken from MathWorld (all except the second contain corner singularities):

>>> mp.dps = 30
>>> f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
>>> print quad(f, [0, 1], [0, 1])
0.577215664901532860606512090082
>>> print euler
0.577215664901532860606512090082

>>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
>>> print quad(f, [-1, 1], [-1, 1])
3.17343648530607134219175646705
>>> print 4*log(2+sqrt(3))-2*pi/3
3.17343648530607134219175646705

>>> f = lambda x, y: 1/(1-x**2 * y**2)
>>> print quad(f, [0, 1], [0, 1])
1.23370055013616982735431137498
>>> print pi**2 / 8
1.23370055013616982735431137498

>>> print quad(lambda x, y: 1/(1-x*y), [0, 1], [0, 1])
1.64493406684822643647241516665
>>> print pi**2 / 6
1.64493406684822643647241516665

Triple integrals are also supported:

>>> mp.dps = 15
>>> print quadgl(lambda x,y,z: x*y/(1+z), [0, 1], [0, 1], [1, 2])
0.101366277027041
>>> print (log(3)-log(2))/4
0.101366277027041

There is currently no direct support for computing quadruple or higher dimensional integrals; if desired, this can be done easily by passing a function that calls quad recursively. While double integrals are reasonably fast, even a simple triple integral at very low precision is likely to take several seconds to evaluate. (This is one application where Gaussian quadrature is often superior.) A quadruple integral will require a whole lot of patience.

Error detection

Neither the tanh-sinh nor the Gauss-Legendre rules are applied adaptively, and do not perform well if there are singularities present or if the integrand is highly oscillatory. (The tanh-sinh algorithm copes well with singularities at the endpoints, but not if they are located in between.)

Problematic integrals can sometimes be evaluated by manually splitting them into smaller pieces.

If the error option is set, quad will return an error estimate along with the result; although this estimate is not always correct, it can be useful for debugging. You can also pass quad the option verbose=True to show detailed progress.

A simple example where both algorithms fail is the function f(x) = abs(sin(x)), which is not smooth at x = pi. In this case, a close value is calculated, but the result is nowhere near the target accuracy; however, quad gives a good estimate of the magnitude of the error:

>>> mp.dps = 15
>>> quadgl(lambda x: abs(sin(x)), [0, 2*pi], error=True)
(mpf('4.0008720436680933'), mpf('0.001'))
>>> quadts(lambda x: abs(sin(x)), [0, 2*pi], error=True)
(mpf('3.9990089417677899'), mpf('0.001'))

This highly oscillatory integral should be pi/2 = 1.57:

>>> print quad(lambda x: sin(x)/x, [0, inf], error=True)
(mpf('0.90019013929034331'), mpf('1.0'))

The next integral should be approximately 0.627 but quadts generates complete nonsense both in the result and the error estimate (the error estimate is somewhat arbitrarily capped at 1.0):

>>> print quad(lambda x: sin(x**2), [0, inf], error=True)
(mpf('7.0886401859297505e+18'), mpf('1.0'))

However, oscillation is not a problem if suppressed by sufficiently fast (preferrably exponential) decay. This integral is exactly 1/2:

>>> print quad(lambda x: exp(-x)*sin(x), [0, inf])
0.5

Another illustrative example is the following double integral, which quadts will process for several seconds before returning a value with very low accuracy:

>>> mpf.dps = 15
>>> f = lambda x, y: sqrt((x-0.5)**2+(y-0.5)**2)
>>> quadts(f, [0, 1], [0, 1], error=1)
(mpf('0.38259743528830825'), mpf('1.0e-6'))

The problem is due to the non-analytic behavior of the function at the midpoint (1/2, 1/2). We can do much better by splitting the area into four pieces (because of the symmetry, we only need to evaluate one of them):

>>> f = lambda x, y: 4*sqrt((x-0.5)**2 + (y-0.5)**2)
>>> print quadts(f, [0.5, 1], [0.5, 1])
0.382597858232106
>>> print (sqrt(2) + asinh(1))/6
0.382597858232106

The value agrees with the known answer and the running time in this case is less than a second on the author’s computer.

Even for analytic integrals on finite intervals, there is no guarantee that quad will be successful. A few examples of integrals for which quad currently fails to reach full accuracy are:

quad(lambda x: sqrt(tan(x)), [0, pi/2])
quad(lambda x: atan(x)/(x*sqrt(1-x**2)), [0, 1])
quad(lambda x: log(1+x**2)/x**2, [0, 1])
quad(lambda x: x**2/((1+x**4)*sqrt(1-x**4)), [0, 1])

(It is possible that future improvements to the quad implementation will make these particular examples work.)

Oscillatory quadrature (quadosc)

It is common to encounter integrals with an integrand of the form f(x) = g(x)*sin(a*x+b), where g(x) is a smoothly decaying function and the range of integration is the entire real line or half the real line. Ordinary quadrature rules are practically useless for integrals of this kind unless g(x) decays exponentially.

Various algorithms more appropriate for oscillatory integration are known, with varying degrees of sophistication. The function quadosc in mpmath uses a very simple strategy: it treats the integral as an infinite alternating series, where each term is the integral over a single half-period. This series will in general converge slowly, but can be accelerated using generic summation algorithms. An advantage of this method is that the oscillation need not be strictly periodic; it is sufficient that the zeros follow a “smooth” distribution.

Specifying the rate of oscillation

The user must supply quadosc with details about the oscillation rate of f(x) via the period or zeros parameter.

If f(x) behaves like a scaled sine or cosine, with a constant distance between each zero, it is sufficient to provide the period. For example, if f sin(3*x), the period should be set to 6*pi. Note that the period is defined as the distance between two consecutive zeros (one period of a sine wave consisting of both one positive and one negative “lobe”).

The parameter zeros is a function that, given an integer n, returns the location where f crosses the axis the nth time, counting from the starting point of integration. To elaborate, the convention is as follows:

  • If integrating from a to inf, zeros(1), zeros(2), ... give the 1st zero to the right of a, the 2nd zero, etc
  • If integrating from -inf to b, zeros(-1), zeros(-2), ... give the 1st zero to the left of b, the 2nd zero, etc
  • If integrating from -inf to inf, zeros(1), zeros(2), ... give the positive zeros (away from 0) and zeros(-1), zeros(-2), ... give the negative zeros

The precise indexing of the zeros is not too significant; it is fine if the first positive zero occurs at zeros(-1) or zeros(2), as long as the zeros are properly ordered.

For a periodic function with period p, it is equivalent to pass period=p or zeros = lambda n: p*n/2 as a parameter.

Examples with simple periodic integrands

The integral of the sinc function sin(x)/x is a perfect match for quadosc. In this case, both the alternating and nonalternating versions work well:

>>> mp.dps = 30
>>> print quadosc(lambda x: sin(x)/x, [0, inf], period=2*pi)
1.57079632679489661923132169164
>>> print pi/2
1.57079632679489661923132169164

Here with zeros instead of period:

>>> print quadosc(lambda x: sin(x)/x, [0, inf], zeros=lambda n: n*pi)
1.57079632679489661923132169164

In the following integral, the period is different, and the first zero is not located at zero (this does not pose a problem):

>>> print quadosc(lambda x: sin(3*x+4)/x**2, [1, inf], period=6*pi)
0.260725538961809481605366600383

The next example involves one of the author’s favorite integrals. The example demonstrates that the integration can be done over the entire real line:

>>> print quadosc(lambda x: cos(x)/(1+x**2), [-inf, inf], period=2*pi)
1.15572734979092171791009318331

The integral has the following neat value:

>>> print pi/e
1.15572734979092171791009318331

Nonperiodic integrands

Important examples of oscillatory integrals with a nonconstant rate of oscillation include those involving Bessel functions. We will illustrate with the J0 function, whose positive zeros occur roughly at 2.40, 5.52, 8.65, and so on. The zeros can be calculated to high accuracy by starting from a simple asymptotic approximation and calling findroot:

>>> def j0zero(n):
...     return findroot(j0, pi*(n-0.25))
...

By definition, the J0 function is normalized so that its integral from 0 to infinity is 1. This should now be easy to verify:

>>> mp.dps = 15
>>> print quadosc(j0, [0, inf], zeros=j0zero)
1.0

We can also try a non-pure Bessel integral:

>>> print quadosc(lambda x: j0(x)/(x+1), [0, inf], zeros=j0zero)
0.754610025770972

Finally, let us consider the integral of sin(x^2) from 0 to infinity (the Fresnel sine integral). Trying to calculate this integral to high accuracy with a generic adaptive quadrature algorithm is a waste of time, and numerical integration based on Fourier analysis does not work either since the integrand is not periodic. But it is easy to locate the zeros, and this is sufficient information for quadosc:

>>> mp.dps = 30
>>> print quadosc(lambda x: sin(x**2), [0, inf], zeros=lambda n: sqrt(n*pi))
0.626657068657750125603941321203

The analytical value of the Fresnel integral is:

>>> print sqrt(pi/8)
0.626657068657750125603941321203

Nonalternating integrands

By default, quadosc rewrites the integral as an alternating series. If the option alt=0 is passed, it instead rewrites the integral as a nonalternating series, with each term being the integral over a full period. This works better if the integrand is oscillatory but does not change signs:

>>> mp.dps = 25
>>> print quadts(lambda x: 1/x**2+sin(x)/x**4, [1, inf])
1.286529507455656027135662
>>> print quadosc(lambda x: 1/x**2+sin(x)/x**4, [1, inf], period=2*pi)
1.280615829835256783092592
>>> print quadosc(lambda x: 1/x**2+sin(x)/x**4, [1, inf], period=2*pi, alt=0)
1.286529535596167393119348

This integral can be evaluated exactly in terms of the cosine integral function, which shows that the last of these three values is correct:

>>> print 1 + (cos(1) + ci(1) + sin(1))/6
1.286529535596167393119348

Solving ODEs (odeint)

Systems of ordinary differential equations can be integrated using the odeint function, which takes as input:

  • A function that computes the vector of derivatives
  • A vector of initial values
  • A vector of positions/time values to step between

It returns the list of function values at each step. As a simple example, we can consider the system x’’ + x = 0. By introducing the variables x = x, y = x’, this second-order differential equation can be rewritten as the first-order system x’ = y, y’ = -x. The analytic solutions with initial values x = 0, y = 1 are of course sin(t) and cos(t):

>>> mp.dps = 25
>>> t = arange(0, 1, 1e-4)
>>> sol = odeint(lambda (x,y),t: (y,-x), (0, 1), t)
>>> print sol[-1][0]
0.8414169503700448479750065
>>> print sin(t[-1])
0.8414169503700448484253423
>>> print sol[-1][1]
0.5403864502649686951810426
>>> print cos(t[-1])
0.5403864502649686944799696

When integrating from 0 to 1, the result with a step size of 1e-4 is in this case accurate to about 18 digits. Currently, odeint uses a simple implementation the fourth-order Runge-Kutta algorithm, with which one must unfortunately set an impractically small step size to obtain much higher precision.