Integration (scipy.integrate)


>>> help(integrate)
 Methods for Integrating Functions given function object.

   quad          -- General purpose integration.
   dblquad       -- General purpose double integration.
   tplquad       -- General purpose triple integration.
   quadrature    -- Integrate with given tolerance using Gaussian quadrature.
   romberg       -- Integrate func using Romberg integration.

 Methods for Integrating Functions given fixed samples.

   trapz         -- Use trapezoidal rule to compute integral from samples.
   cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
   simps         -- Use Simpson's rule to compute integral from samples.
   romb          -- Use Romberg Integration to compute integral from

   See the special module's orthogonal polynomials (special) for Gaussian
      quadrature roots and weights for other weighting factors and regions.

 Interface to numerical integrators of ODE systems.

   odeint        -- General integration of ordinary differential equations.
   ode           -- Integrate ODE using VODE and ZVODE routines.

General integration (quad)

The function quad is provided to integrate a function of one variable between two points. The points can be \(\pm\infty\) (\(\pm\) inf) to indicate infinite limits. For example, suppose you wish to integrate a bessel function jv(2.5,x) along the interval \([0,4.5].\)


This could be computed using quad:

>>> print result
(1.1178179380783249, 7.8663172481899801e-09)
>>> print I





is the Fresnel sine integral. Note that the numerically-computed integral is within \(1.04\times10^{-11}\) of the exact result — well below the reported error bound.


\[I(a,b)=\int_{0}^{1} ax^2+b \, dx.\]

This integral can be evaluated by using the following code:

>>> from scipy.integrate import quad
...     return a * x + b
>>> a = 2
>>> b = 1
>>> I = (2.0, 2.220446049250313e-14)

Infinite inputs are also allowed in quad by using \(\pm\) inf as one of the arguments. For example, suppose that a numerical value for the exponential integral:



>>> from scipy.integrate import quad
>>> def expint(n, x):
>>> vec_expint = vectorize(expint)
>>> vec_expint(3,arange(1.0,4.0,0.5))
>>> special.expn(3,arange(1.0,4.0,0.5))


>>> print result
(0.33333333324560266, 2.8548934485373678e-09)
>>> I3 = 1.0/3.0
>>> print I3
>>> print I3 - result[0]


General multiple integration (dblquad, tplquad, nquad)


An example of using double integration to compute several values of \(I_{n}\) is shown below:

>>> from scipy.integrate import quad, dblquad
>>> def I(n):
>>> print I(4)
(0.25000000000435768, 1.0518245707751597e-09)
>>> print I(3)
(0.33333333325010883, 2.8604069919261191e-09)
>>> print I(2)
(0.49999999999857514, 1.8855523253868967e-09)

As example for non-constant limits consider the integral



>>> from scipy.integrate import dblquad
>>> area
(0.010416666666666668, 1.1564823173178715e-16)


The integral from above


can be calculated as

>>> from scipy import integrate
>>> N = 5
>>> def f(t, x):
(0.20000000000002294, 1.2239614263187945e-08)

Note that the order of arguments for f must match the order of the integration bounds; i.e. the inner integral with respect to \(t\) is on the interval \([1, \infty]\) and the outer integral with respect to \(x\) is on the interval \([0, \infty]\).



can be evaluated by means of

>>> from scipy import integrate
>>> def f(x, y):
>>>     return x*y
>>> def bounds_y():
>>>     return [0, 0.5]
>>> def bounds_x(y):
(0.010416666666666668, 4.101620128472366e-16)

which is the same result as before.

Gaussian quadrature


Romberg Integration


Integrating using Samples

If the samples are equally-spaced and the number of samples available is \(2^{k}+1\) for some integer \(k\), then Romberg romb integration can be used to obtain high-precision estimates of the integral using the available samples. Romberg integration uses the trapezoid rule at step-sizes related by a power of two and then performs Richardson extrapolation on these estimates to approximate the integral with a higher-degree of accuracy.



>>> from scipy.integrate import simps
>>> import numpy as np
>>> def f(x):
...    return x**2
>>> def f2(x):
...    return x**3
>>> x = np.array([1,3,4])
>>> y1 = f1(x)
>>> I1 = integrate.simps(y1, x)
>>> print(I1)

This corresponds exactly to

\[\int_{1}^{4} x^2 \, dx = 21,\]

whereas integrating the second function

>>> y2 = f2(x)
>>> I2 = integrate.simps(y2, x)
>>> print(I2)

does not correspond to

\[\int_{1}^{4} x^3 \, dx = 63.75\]


Faster integration using Ctypes

A user desiring reduced integration times may pass a C function pointer through ctypes to quad, dblquad, tplquad or nquad and it will be integrated and return a result in Python. The performance increase here arises from two factors. The primary improvement is faster function evaluation, which is provided by compilation. This can also be achieved using a library like Cython or F2Py that compiles Python. Additionally we have a speedup provided by the removal of function calls between C and Python in quad - this cannot be achieved through Cython or F2Py. This method will provide a speed increase of ~2x for trivial functions such as sine but can produce a much more noticeable increase (10x+) for more complex functions. This feature then, is geared towards a user with numerically intensive integrations willing to write a little C to reduce computation time significantly.

ctypes integration can be done in a few simple steps:

1.) Write an integrand function in C with the function signature double f(int n, double args[n]), where args is an array containing the arguments of the function f.

double f(int n, double args[n]){
    return args[0] - args[1] * args[2]; //corresponds to x0 - x1 * x2


$ gcc -shared -o -fPIC testlib.c



>>> import ctypes
>>> from scipy import integrate
>>> lib = ctypes.CDLL('/**/') # Use absolute path to testlib
>>> func = lib.f # Assign specific function to name func (for simplicity)
>>> func.restype = ctypes.c_double
>>> func.argtypes = (ctypes.c_int, ctypes.c_double)

Note that the argtypes will always be (ctypes.c_int, ctypes.c_double) regardless of the number of parameters, and restype will always be ctypes.c_double.


>>> integrate.nquad(func, [[0,10],[-10,0],[-1,1]])
(1000.0, 1.1102230246251565e-11)

And the Python tuple is returned as expected in a reduced amount of time. All optional parameters can be used with this method including specifying singularities, infinite bounds, etc.

Ordinary differential equations (odeint)



given initial conditions \(\mathbf{y}\left(0\right)=y_{0}\), where \(\mathbf{y}\) is a length \(N\) vector and \(\mathbf{f}\) is a mapping from \(\mathcal{R}^{N}\) to \(\mathcal{R}^{N}.\) A higher-order ordinary differential equation can always be reduced to a differential equation of this type by introducing intermediate derivatives into the \(\mathbf{y}\) vector.



with initial conditions \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) and \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.\) It is known that the solution to this differential equation with these boundary conditions is the Airy function


which gives a means to check the integrator using special.airy.

First, convert this ODE into standard form by setting \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) and \(t=z\). Thus, the differential equation becomes


In other words,


As an interesting reminder, if \(\mathbf{A}\left(t\right)\) commutes with \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) under matrix multiplication, then this linear differential equation has an exact solution using the matrix exponential:


However, in this case, \(\mathbf{A}\left(t\right)\) and its integral do not commute.


The following example illustrates the use of odeint including the usage of the Dfun option which allows the user to specify a gradient (with respect to \(\mathbf{y}\) ) of the function, \(\mathbf{f}\left(\mathbf{y},t\right)\).

>>> from scipy.integrate import odeint
>>> from scipy.special import gamma, airy
>>> y0 = [y0_0, y1_0]
>>> def func(y, t):
>>> def gradient(y, t):
>>> x = arange(0, 4.0, 0.01)
>>> t = x
>>> ychk = airy(x)[0]
>>> print ychk[:36:6]
[ 0.355028  0.339511  0.324068  0.308763  0.293658  0.278806]
>>> print y[:36:6,1]
[ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806]
>>> print y2[:36:6,1]
[ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806]