Integration

Tanh-sinh quadrature (quadts)

The function quadts performs numerical integration (quadrature) using the tanh-sinh algorithm. The syntax for integrating a function f between the endpoints a and b is quadts(f, a, b). For example:

>>> from mpmath import *
>>> print quadts(sin, 0, pi)
2.0

Tanh-sinh quadrature is extremely efficient for high-precision integration of analytic functions. Unlike the more well-known Gaussian quadrature algorithm, it is relatively insensitive to integrable singularities at the endpoints of the interval. The quadts function attempts to evaluate the integral to the full working precision; for example, it can calculate 100 digits of pi by integrating the area under the half circle arc x^2 + y^2 = 1 (y > 0):

>>> mp.dps = 100
>>> print quadts(lambda x: 2*sqrt(1 - x**2), -1, 1)
... # doctest:+ELLIPSIS
3.14159265358979323846264338327950288419716939937510582097...

Neat examples

Intervals can be infinite or half-infinite:

>>> mp.dps = 15
>>> print quadts(lambda x: 2/(x**2+1), 0, inf)
3.14159265358979
>>> print quadts(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 - 3.64632996762876e-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 quadts. To do this, simply provide a two-argument function and, instead of two endpoints, provide 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 quadts(f, (0, 1), (0, 1))
0.577215664901532860606512090082
>>> print euler
0.577215664901532860606512090082

>>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
>>> print quadts(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 quadts(f, (0, 1), (0, 1))
1.23370055013616982735431137498
>>> print pi**2 / 8
1.23370055013616982735431137498

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

There is currently no direct support for computing triple or higher dimensional integrals; if desired, this can be done easily by passing a function that calls quadts recursively:

>>> mp.dps = 15
>>> f = lambda x, y: quadts(lambda z: sin(x)/z+y*z, 1, 2)
>>> print quadts(f, (1, 2), (1, 2))
2.91296002641413
>>> print mpf(9)/4 + (cos(1)-cos(2))*log(2)
2.91296002641413

While double integrals are reasonably fast, even a simple triple integral at very low precision is likely to take several seconds to evaluate (harder integrals may take minutes). A quadruple integral will require a whole lot of patience.

Error detection

The tanh-sinh algorithm is not suitable for adaptive quadrature, and does not perform well if there are singularities between the endpoints or if the integrand is oscillatory (such integrals should manually be split into smaller pieces). If the error option is set, quadts 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 quadts the option verbose=True to show detailed progress.

A simple example where the algorithm fails 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, quadts gives a good estimate of the magnitude of the error:

>>> mp.dps = 15
>>> 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 quadts(lambda x: sin(x)/x, 0, inf, error=True)
(mpf('2.3840907358976577'), 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 quadts(lambda x: sin(x**2), 0, inf, error=True)
(mpf('2.5190134849122411e+21'), mpf('1.0'))

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

>>> print quadts(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.38259743528830826'), 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 just 0.7 seconds on the author’s computer.

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

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

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

Performance

The tanh-sinh scheme is efficient enough that analytic 100-digit integrals often can often be evaluated in less than a second. The timings for quadts(lambda x: 2*sqrt(1-x**2), -1, 1) at various precision levels on the author’s computer is:

dps First evaluation Second evaluation
15 0.017 seconds 0.0054 seconds
50 0.085 seconds 0.016 seconds
500 8.9 seconds 0.48 seconds

The second integration at the same precision level is much faster. The reason for this is that the tanh-sinh algorithm must be initalized by computing a set of nodes, and this initalization if often more expensive than actually evaluating the integral. Mpmath automatically caches all computed nodes to make subsequent integrations faster, but the cache is lost when Python shuts down, so if you would frequently like to use mpmath to calculate 1000-digit integrals, you may want to save the nodes to a file. The nodes are stored in a dict TS_cache located in the mpmath.calculus module, which can be pickled if desired.

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. The tanh-sinh quadrature rule is 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 is remarkable because of having 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 the secant root-finder:

>>> def j0zero(n):
...     return secant(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.286529526945860318388688
>>> 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.

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.841416950370044847975007
>>> print sin(t[-1])
0.8414169503700448484253423
>>> print sol[-1][1]
0.5403864502649686951810432
>>> 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.