scipy.integrate.nquad

scipy.integrate.nquad(func, ranges, args=None, opts=None)

Integration over multiple variables.

Wraps quad to enable integration over multiple variables. Various options allow improved integration of discontinuous functions, as well as the use of weighted integration, and generally finer control of the integration process.

Parameters:

func : callable

The function to be integrated. Has arguments of x0, ... xn, t0, tm, where integration is carried out over x0, ... xn, which must be floats. Function signature should be func(x0, x1, ..., xn, t0, t1, ..., tm). Integration is carried out in order. That is, integration over x0 is the innermost integral, and xn is the outermost. If performance is a concern, this function may be a ctypes function of the form:

f(int n, double args[n])

where n is the number of extra parameters and args is an array of doubles of the additional parameters. This function may then be compiled to a dynamic/shared library then imported through ctypes, setting the function’s argtypes to (c_int, c_double), and the function’s restype to (c_double). Its pointer may then be passed into nquad normally. This allows the underlying Fortran library to evaluate the function in the innermost integration calls without callbacks to Python, and also speeds up the evaluation of the function itself.

ranges : iterable object

Each element of ranges may be either a sequence of 2 numbers, or else a callable that returns such a sequence. ranges[0] corresponds to integration over x0, and so on. If an element of ranges is a callable, then it will be called with all of the integration arguments available. e.g. if func = f(x0, x1, x2), then ranges[0] may be defined as either (a, b) or else as (a, b) = range0(x1, x2).

args : iterable object, optional

Additional arguments t0, ..., tn, required by func.

opts : iterable object or dict, optional

Options to be passed to quad. May be empty, a dict, or a sequence of dicts or functions that return a dict. If empty, the default options from scipy.integrate.quadare used. If a dict, the same options are used for all levels of integraion. If a sequence, then each element of the sequence corresponds to a particular integration. e.g. opts[0] corresponds to integration over x0, and so on. The available options together with their default values are:

• epsabs = 1.49e-08
• epsrel = 1.49e-08
• limit = 50
• points = None
• weight = None
• wvar = None
• wopts = None

The full_output option from quad is unavailable, due to the complexity of handling the large amount of data such an option would return for this kind of nested integration. For more information on these options, see quad and quad_explain.

Returns:

result : float

The result of the integration.

abserr : float

The maximum of the estimates of the absolute error in the various integration results.

See also

quad

1-dimensional numerical integration

dblquad, tplquad

fixed_quad

fixed-order Gaussian quadrature

quadrature

adaptive Gaussian quadrature

Examples

>>> from scipy import integrate
>>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
...                                 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
>>> points = [[lambda (x1,x2,x3) : 0.2*x3 + 0.5 + 0.25*x1], [], [], []]
>>> def opts0(*args, **kwargs):
...     return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
>>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
...                 opts=[opts0,{},{},{}])
(1.5267454070738633, 2.9437360001402324e-14)

>>> scale = .1
>>> def func2(x0, x1, x2, x3, t0, t1):
...     return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0)
>>> def lim0(x1, x2, x3, t0, t1):
...     return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
...             scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
>>> def lim1(x2, x3, t0, t1):
...     return [scale * (t0*x2 + t1*x3) - 1,
...             scale * (t0*x2 + t1*x3) + 1]
>>> def lim2(x3, t0, t1):
...     return [scale * (x3 + t0**2*t1**3) - 1,
...             scale * (x3 + t0**2*t1**3) + 1]
>>> def lim3(t0, t1):
...     return [scale * (t0+t1) - 1, scale * (t0+t1) + 1]
>>> def opts0(x1, x2, x3, t0, t1):
...     return {'points' : [t0 - t1*x1]}
>>> def opts1(x2, x3, t0, t1):
...     return {}
>>> def opts2(x3, t0, t1):
...     return {}
>>> def opts3(t0, t1):
...     return {}
>>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
...                 opts=[opts0, opts1, opts2, opts3])
(25.066666666666666, 2.7829590483937256e-13)