scipy.integrate.

求积#

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

多维数组值函数的自适应求积。

给定任意积分规则，此函数返回在由指定超立方体顶点的数组 a 和 b 定义的区域上，达到请求容差的积分估计值。

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

参数:
f可调用对象

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

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

f 应该接受形状为

(npoints, ndim)

并输出形状为

(npoints, output_dim_1, ..., output_dim_n)

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

(output_dim_1, ..., output_dim_n)
a, b类数组

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

rule字符串，可选

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

rtol, atol浮点数，可选

相对和绝对容差。迭代将执行，直到估计误差小于 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_subdivisions整数，可选

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

args元组，可选

传递给 f 的额外位置参数（如果有）。

workers整数或类映射可调用对象，可选

如果 workers 是一个整数，则部分计算将并行完成，细分为指定数量的任务（使用 multiprocessing.pool.Pool）。提供 -1 以使用进程可用的所有核心。或者，提供一个类映射可调用对象，例如 multiprocessing.pool.Pool.map，用于并行评估总体。此评估以 workers(func, iterable) 的形式执行。

points类数组列表，可选

要避免评估 f 的点列表，条件是所使用的规则不在区域边界上评估 f（Genz-Malik 和 Gauss-Kronrod 规则都是如此）。如果 f 在指定点有奇点，这可能很有用。这应该是一个类数组列表，其中每个元素的长度为 ndim。默认为空。参见示例。

返回:
res对象

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

estimatendarray

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

errorndarray

在指定整体区域上的近似误差估计。

status字符串

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

subdivisions整数

执行的细分数量。

atol, rtol浮点数

近似所需的容差。

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 维被积函数（其中 n > 2）使用 Gauss-Kronrod，则通过取一维情况下节点的笛卡尔积来获得相应的笛卡尔积规则。这意味着在 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 维矩形区域数值积分的自适应算法, 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

示例

带矢量输出的一维积分:

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

带无限极限的 2 维积分:

\[\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
在此页上