cubature#
- scipy.integrate.cubature(f, a, b, *, rule='gk21', rtol=1e-08, atol=0, max_subdivisions=10000, args=(), workers=1)[source]#
Adaptive cubature of multidimensional array-valued function.
Given an arbitrary integration rule, this function returns an estimate of the integral to the requested tolerance over the region defined by the arrays a and b specifying the corners of a hypercube.
Convergence is not guaranteed for all integrals.
- Parameters:
- fcallable
Function to integrate. f must have the signature:
f(x : ndarray, *args) -> ndarray
f should accept arrays
x
of shape:(npoints, ndim)
and output arrays of shape:
(npoints, output_dim_1, ..., output_dim_n)
In this case,
cubature
will return arrays of shape:(output_dim_1, ..., output_dim_n)
- a, barray_like
Lower and upper limits of integration as 1D arrays specifying the left and right endpoints of the intervals being integrated over. Infinite limits are currently not supported.
- rulestr, optional
Rule used to estimate the integral. If passing a string, the options are “gauss-kronrod” (21 node), or “genz-malik” (degree 7). If a rule like “gauss-kronrod” is specified for an
n
-dim integrand, the corresponding Cartesian product rule is used. “gk21”, “gk15” are also supported for compatibility withquad_vec
. See Notes.- rtol, atolfloat, optional
Relative and absolute tolerances. Iterations are performed until the error is estimated to be less than
atol + rtol * abs(est)
. Here rtol controls relative accuracy (number of correct digits), while atol controls absolute accuracy (number of correct decimal places). To achieve the desired rtol, set atol to be smaller than the smallest value that can be expected fromrtol * abs(y)
so that rtol dominates the allowable error. If atol is larger thanrtol * abs(y)
the number of correct digits is not guaranteed. Conversely, to achieve the desired atol, set rtol such thatrtol * abs(y)
is always smaller than atol. Default values are 1e-8 for rtol and 0 for atol.- max_subdivisionsint, optional
Upper bound on the number of subdivisions to perform. Default is 10,000.
- argstuple, optional
Additional positional args passed to f, if any.
- workersint or map-like callable, optional
If workers is an integer, part of the computation is done in parallel subdivided to this many tasks (using
multiprocessing.pool.Pool
). Supply -1 to use all cores available to the Process. Alternatively, supply a map-like callable, such asmultiprocessing.pool.Pool.map
for evaluating the population in parallel. This evaluation is carried out asworkers(func, iterable)
.
- Returns:
- resobject
Object containing the results of the estimation. It has the following attributes:
- estimatendarray
Estimate of the value of the integral over the overall region specified.
- errorndarray
Estimate of the error of the approximation over the overall region specified.
- statusstr
Whether the estimation was successful. Can be either: “converged”, “not_converged”.
- subdivisionsint
Number of subdivisions performed.
- atol, rtolfloat
Requested tolerances for the approximation.
- regions: list of object
List of objects containing the estimates of the integral over smaller regions of the domain.
Each object in
regions
has the following attributes:- a, bndarray
Points describing the corners of the region. If the original integral contained infinite limits or was over a region described by region, then a and b are in the transformed coordinates.
- estimatendarray
Estimate of the value of the integral over this region.
- errorndarray
Estimate of the error of the approximation over this region.
Notes
The algorithm uses a similar algorithm to
quad_vec
, which itself is based on the implementation of QUADPACK’s DQAG* algorithms, implementing global error control and adaptive subdivision.The source of the nodes and weights used for Gauss-Kronrod quadrature can be found in [1], and the algorithm for calculating the nodes and weights in Genz-Malik cubature can be found in [2].
The rules currently supported via the rule argument are:
"gauss-kronrod"
, 21-node Gauss-Kronrod"genz-malik"
, n-node Genz-Malik
If using Gauss-Kronrod for an
n
-dim integrand wheren > 2
, then the corresponding Cartesian product rule will be found by taking the Cartesian product of the nodes in the 1D case. This means that the number of nodes scales exponentially as21^n
in the Gauss-Kronrod case, which may be problematic in a moderate number of dimensions.Genz-Malik is typically less accurate than Gauss-Kronrod but has much fewer nodes, so in this situation using “genz-malik” might be preferable.
References
[1]R. Piessens, E. de Doncker, Quadpack: A Subroutine Package for Automatic Integration, files: dqk21.f, dqk15.f (1983).
[2]A.C. Genz, A.A. Malik, Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region, Journal of Computational and Applied Mathematics, Volume 6, Issue 4, 1980, Pages 295-302, ISSN 0377-0427 DOI:10.1016/0771-050X(80)90039-X
Examples
1D integral with vector output:
\[\int^1_0 \mathbf f(x) \text dx\]Where
f(x) = x^n
andn = np.arange(10)
is a vector. Since no rule is specified, the default “gk21” is used, which corresponds to Gauss-Kronrod integration with 21 nodes.>>> import numpy as np >>> from scipy.integrate import cubature >>> def f(x, n): ... # Make sure x and n are broadcastable ... return x[:, np.newaxis]**n[np.newaxis, :] >>> res = cubature( ... f, ... a=[0], ... b=[1], ... args=(np.arange(10),), ... ) >>> res.estimate array([1. , 0.5 , 0.33333333, 0.25 , 0.2 , 0.16666667, 0.14285714, 0.125 , 0.11111111, 0.1 ])
7D integral with arbitrary-shaped array output:
f(x) = cos(2*pi*r + alphas @ x)
for some
r
andalphas
, and the integral is performed over the unit hybercube, \([0, 1]^7\). Since the integral is in a moderate number of dimensions, “genz-malik” is used rather than the default “gauss-kronrod” to avoid constructing a product rule with \(21^7 \approx 2 \times 10^9\) nodes.>>> import numpy as np >>> from scipy.integrate import cubature >>> def f(x, r, alphas): ... # f(x) = cos(2*pi*r + alphas @ x) ... # Need to allow r and alphas to be arbitrary shape ... npoints, ndim = x.shape[0], x.shape[-1] ... alphas = alphas[np.newaxis, ...] ... x = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim) ... return np.cos(2*np.pi*r + np.sum(alphas * x, axis=-1)) >>> rng = np.random.default_rng() >>> r, alphas = rng.random((2, 3)), rng.random((2, 3, 7)) >>> res = cubature( ... f=f, ... a=np.array([0, 0, 0, 0, 0, 0, 0]), ... b=np.array([1, 1, 1, 1, 1, 1, 1]), ... rtol=1e-5, ... rule="genz-malik", ... args=(r, alphas), ... ) >>> res.estimate array([[-0.79812452, 0.35246913, -0.52273628], [ 0.88392779, 0.59139899, 0.41895111]])
Parallel computation with workers:
>>> from concurrent.futures import ThreadPoolExecutor >>> with ThreadPoolExecutor() as executor: ... res = cubature( ... f=f, ... a=np.array([0, 0, 0, 0, 0, 0, 0]), ... b=np.array([1, 1, 1, 1, 1, 1, 1]), ... rtol=1e-5, ... rule="genz-malik", ... args=(r, alphas), ... workers=executor.map, ... ) >>> res.estimate array([[-0.79812452, 0.35246913, -0.52273628], [ 0.88392779, 0.59139899, 0.41895111]])