Bug
Links
def cubature(f, a, b, *, rule='gk21', rtol=1e-08, atol=0, max_subdivisions=10000, args=(), workers=1, points=None)
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
----------
f : callable
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, b : array_like
Lower and upper limits of integration as 1D arrays specifying the left and right
endpoints of the intervals being integrated over. Limits can be infinite.
rule : str, 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 with `quad_vec`. See Notes.
rtol, atol : float, 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 from
``rtol * abs(y)`` so that `rtol` dominates the allowable error. If `atol` is
larger than ``rtol * abs(y)`` the number of correct digits is not guaranteed.
Conversely, to achieve the desired `atol`, set `rtol` such that
``rtol * abs(y)`` is always smaller than `atol`. Default values are 1e-8 for
`rtol` and 0 for `atol`.
max_subdivisions : int, optional
Upper bound on the number of subdivisions to perform. Default is 10,000.
args : tuple, optional
Additional positional args passed to `f`, if any.
workers : int or map-like callable, optional
If `workers` is an integer, part of the computation is done in parallel
subdivided to this many tasks (using :class:`python:multiprocessing.pool.Pool`).
Supply `-1` to use all cores available to the Process. Alternatively, supply a
map-like callable, such as :meth:`python:multiprocessing.pool.Pool.map` for
evaluating the population in parallel. This evaluation is carried out as
``workers(func, iterable)``.
points : list of array_like, optional
List of points to avoid evaluating `f` at, under the condition that the rule
being used does not evaluate `f` on the boundary of a region (which is the
case for all Genz-Malik and Gauss-Kronrod rules). This can be useful if `f` has
a singularity at the specified point. This should be a list of array-likes where
each element has length ``ndim``. Default is empty. See Examples.
Returns
-------
res : object
Object containing the results of the estimation. It has the following
attributes:
estimate : ndarray
Estimate of the value of the integral over the overall region specified.
error : ndarray
Estimate of the error of the approximation over the overall region
specified.
status : str
Whether the estimation was successful. Can be either: "converged",
"not_converged".
subdivisions : int
Number of subdivisions performed.
atol, rtol : float
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, b : ndarray
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.
estimate : ndarray
Estimate of the value of the integral over this region.
error : ndarray
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 where ``n > 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 as ``21^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.
Infinite limits are handled with an appropriate variable transformation. Assuming
``a = [a_1, ..., a_n]`` and ``b = [b_1, ..., b_n]``:
If :math:`a_i = -\infty` and :math:`b_i = \infty`, the i-th integration variable
will use the transformation :math:`x = \frac{1-|t|}{t}` and :math:`t \in (-1, 1)`.
If :math:`a_i \ne \pm\infty` and :math:`b_i = \infty`, the i-th integration variable
will use the transformation :math:`x = a_i + \frac{1-t}{t}` and
:math:`t \in (0, 1)`.
If :math:`a_i = -\infty` and :math:`b_i \ne \pm\infty`, the i-th integration
variable will use the transformation :math:`x = b_i - \frac{1-t}{t}` and
:math:`t \in (0, 1)`.
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**:
.. math::
\int^1_0 \mathbf f(x) \text dx
Where ``f(x) = x^n`` and ``n = 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`` and ``alphas``, and the integral is performed over the unit
hybercube, :math:`[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 :math:`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]])
**2D integral with infinite limits**:
.. math::
\int^{ \infty }_{ -\infty }
\int^{ \infty }_{ -\infty }
e^{-x^2-y^2}
\text dy
\text dx
>>> def gaussian(x):
... return np.exp(-np.sum(x**2, axis=-1))
>>> res = cubature(gaussian, [-np.inf, -np.inf], [np.inf, np.inf])
>>> res.estimate
3.1415926
**1D integral with singularities avoided using** `points`:
.. math::
\int^{ 1 }_{ -1 }
\frac{\sin(x)}{x}
\text dx
It is necessary to use the `points` parameter to avoid evaluating `f` at the origin.
>>> def sinc(x):
... return np.sin(x)/x
>>> res = cubature(sinc, [-1], [1], points=[[0]])
>>> res.estimate
1.8921661
Améliorations / Corrections
Vous avez des améliorations (ou des corrections) à proposer pour ce document : je vous remerçie par avance de m'en faire part, cela m'aide à améliorer le site.
Emplacement :
Description des améliorations :