scipy.integrate.

cubature

scipy.integrate.cubature(f, a, b, *, rule='gk21', rtol=1e-08, atol=0, max_subdivisions=10000, args=(), workers=1, points=None)

多维数组值函数的自适应立方体积分。

给定任意积分规则，此函数返回在由数组 a 和 b 定义的超立方体角点所定义的区域上的积分估计值，该估计值满足请求的公差。

并非所有积分都能保证收敛。

参数:

fcallable

要积分的函数。f 必须具有以下签名

f(x : ndarray, *args) -> ndarray

f 应该接受形状为 x 的数组

(npoints, ndim)

并输出形状为

(npoints, output_dim_1, ..., output_dim_n)

在这种情况下，cubature 将返回形状为

(output_dim_1, ..., output_dim_n)

a, barray_like

积分的下限和上限，以 1D 数组形式指定要积分的区间左右端点。极限可以是无限的。

rulestr, 可选

用于估计积分的规则。如果传递字符串，则选项为“gauss-kronrod”（21 个节点）或“genz-malik”（7 度）。如果为 n 维被积函数指定了类似“gauss-kronrod”的规则，则使用相应的笛卡尔积规则。“gk21”、“gk15”也支持与 quad_vec 兼容。请参阅注释。

rtol, atolfloat, 可选

相对和绝对公差。执行迭代，直到误差估计值小于 atol + rtol * abs(est)。此处，rtol 控制相对精度（正确位数），而 atol 控制绝对精度（正确的小数位数）。要实现所需的 rtol，请将 atol 设置为小于预期 rtol * abs(y) 的最小值，以便 rtol 控制允许的误差。如果 atol 大于 rtol * abs(y)，则不保证正确的位数。相反，要实现所需的 atol，请设置 rtol，使 rtol * abs(y) 始终小于 atol。rtol 的默认值为 1e-8，atol 的默认值为 0。

max_subdivisionsint, 可选

执行的最大细分次数的上限。默认值为 10,000。

argstuple, 可选

传递给 f 的其他位置参数（如果有）。

workersint 或类似 map 的可调用对象, 可选

如果 workers 是一个整数，则部分计算将并行完成，并细分为这么多任务（使用 multiprocessing.pool.Pool）。提供 -1 以使用进程可用的所有核心。或者，提供类似 map 的可调用对象，例如 multiprocessing.pool.Pool.map，以并行评估种群。此评估将作为 workers(func, iterable) 执行。

pointsarray_like 的列表, 可选

在以下条件下，避免在 f 处评估的点的列表：所使用的规则不会在区域边界上评估 f（所有 Genz-Malik 和 Gauss-Kronrod 规则都是这种情况）。如果 f 在指定点处存在奇点，这将很有用。这应该是一个 array-likes 的列表，其中每个元素的长度为 ndim。默认值为空。请参阅示例。

返回:

resobject

包含估计结果的对象。它具有以下属性

estimatendarray

指定区域上的积分值的估计值。

errorndarray

指定区域上近似值的误差估计值。

statusstr

估计是否成功。可以是：“converged”（已收敛）、“not_converged”（未收敛）。

subdivisionsint

执行的细分次数。

atol, rtolfloat

请求的近似值公差。

regions: 对象的列表

包含域较小区域上积分估计值的对象列表。

regions 中的每个对象都具有以下属性

a, bndarray

描述区域角点的点。如果原始积分包含无限限制或位于由 region 描述的区域上，则 a 和 b 位于转换后的坐标中。

estimatendarray

此区域上的积分值的估计值。

errorndarray

此区域上近似值的误差估计值。

注释

该算法使用与 quad_vec 类似的算法，该算法本身基于 QUADPACK 的 DQAG* 算法的实现，实现了全局误差控制和自适应细分。

用于 Gauss-Kronrod 正交的节点和权重的来源可以在 [1] 中找到，而用于计算 Genz-Malik 立方体中的节点和权重的算法可以在 [2] 中找到。

当前通过 rule 参数支持的规则是

• "gauss-kronrod", 21 节点 Gauss-Kronrod

• "genz-malik"，n 节点 Genz-Malik

如果对 n 维被积函数使用 Gauss-Kronrod，其中 n > 2，则通过取 1D 情况下的节点笛卡尔积来找到相应的笛卡尔积规则。这意味着在 Gauss-Kronrod 情况下的节点数会呈指数级增长，即 21^n，这在维度数量适中的情况下可能会出现问题。

Genz-Malik 通常不如 Gauss-Kronrod 精确，但节点少得多，因此在这种情况下，使用“genz-malik”可能更可取。

无限限值通过适当的变量变换来处理。假设 a = [a_1, ..., a_n] 和 b = [b_1, ..., b_n]

如果 \(a_i = -\infty\) 且 \(b_i = \infty\)，则第 i 个积分变量将使用变换 \(x = \frac{1-|t|}{t}\)，且 \(t \in (-1, 1)\)。

如果 \(a_i \ne \pm\infty\) 且 \(b_i = \infty\)，则第 i 个积分变量将使用变换 \(x = a_i + \frac{1-t}{t}\)，且 \(t \in (0, 1)\)。

如果 \(a_i = -\infty\) 且 \(b_i \ne \pm\infty\)，则第 i 个积分变量将使用变换 \(x = b_i - \frac{1-t}{t}\)，且 \(t \in (0, 1)\)。

参考文献

[1]

R. Piessens, E. de Doncker, Quadpack: 一个自动积分的子程序包, 文件: dqk21.f, dqk15.f (1983).

[2]

A.C. Genz, A.A. Malik, 关于算法 006 的评论：一个用于 N 维矩形区域数值积分的自适应算法, 计算与应用数学杂志, 第 6 卷, 第 4 期, 1980, 第 295-302 页, ISSN 0377-0427 DOI:10.1016/0771-050X(80)90039-X

示例

在您的浏览器中尝试一下！

具有向量输出的一维积分:

\[\int^1_0 \mathbf f(x) \text dx\]

其中 f(x) = x^n 且 n = np.arange(10) 是一个向量。由于没有指定规则，将使用默认的“gk21”，它对应于具有 21 个节点的 Gauss-Kronrod 积分。

>>> 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       ])

具有任意形状数组输出的 7 维积分:

f(x) = cos(2*pi*r + alphas @ x)

对于某些 r 和 alphas，积分在单位超立方体 \([0, 1]^7\) 上执行。由于积分处于中等维度，因此使用“genz-malik”而不是默认的“gauss-kronrod”，以避免构建具有 \(21^7 \approx 2 \times 10^9\) 个节点的乘积规则。

>>> 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]])

使用 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]])

具有无限限值的二维积分:

\[\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

使用 points 避免奇点的一维积分

\[\int^{ 1 }_{ -1 } \frac{\sin(x)}{x} \text dx\]

有必要使用 points 参数来避免在原点评估 f。

>>> def sinc(x):
...     return np.sin(x)/x
>>> res = cubature(sinc, [-1], [1], points=[[0]])
>>> res.estimate
 1.8921661

返回在新标签页中打开