Ctrl+K
scipy.integrate.

cubature#

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

多維陣列值函數的自適應 cubature 積分。

給定任意積分規則，此函數返回在由陣列 a 和 b 定義的超立方體角落區域上，達到請求容差的積分估計值。

並非所有積分都保證收斂。

參數:
f可呼叫物件 (callable)

要積分的函數。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, b類陣列 (array_like)

積分的下限和上限，以 1 維陣列形式指定要積分區間的左右端點。 限制可以是無限的。

rule字串 (str)，選用

用於估計積分的規則。 如果傳遞字串，選項為 “gauss-kronrod”（21 節點）或 “genz-malik”（7 度）。 如果為 n 維被積函數指定了像 “gauss-kronrod” 這樣的規則，則會使用對應的笛卡爾積規則。“gk21”、“gk15” 也被支援，以與 quad_vec 相容。 請參閱「Notes」部分。

rtol, atol浮點數 (float)，選用

相對和絕對容差。 迭代會執行直到誤差估計小於 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整數 (int)，選用

要執行的細分次數上限。 預設值為 10,000。

args元組 (tuple)，選用

要傳遞給 f 的額外位置引數（如果有的話）。

workers整數 (int) 或類 map 可呼叫物件 (callable)，選用

如果 workers 是整數，則部分計算會並行完成，並細分為這麼多任務（使用 multiprocessing.pool.Pool）。 提供 -1 以使用 Process 可用的所有核心。 或者，提供類似 map 的可呼叫物件，例如 multiprocessing.pool.Pool.map，以並行評估族群。 此評估會以 workers(func, iterable) 的形式執行。

points類陣列 (array_like) 的列表，選用

避免在某些點評估 f 的點列表，條件是所使用的規則不會在區域邊界上評估 f（所有 Genz-Malik 和 Gauss-Kronrod 規則都是這種情況）。 如果 f 在指定點具有奇異點，則這可能很有用。 這應該是類陣列的列表，其中每個元素的長度為 ndim。 預設為空。 請參閱「Examples」部分。

返回:
res物件 (object)

包含估計結果的物件。 它具有以下屬性

estimatendarray

在指定的整體區域上，積分值的估計。

errorndarray

在指定的整體區域上，近似誤差的估計。

status字串 (str)

估計是否成功。 可以是：“converged”、“not_converged”。

subdivisions整數 (int)

執行的細分次數。

atol, rtol浮點數 (float)

請求的近似容差。

regions: 物件 (object) 列表

包含網域較小區域上積分估計值的物件列表。

regions 中的每個物件都具有以下屬性

a, bndarray

描述區域角落的點。 如果原始積分包含無限限制，或是在由 region 描述的區域上進行積分，則 a 和 b 將會是轉換後的座標。

estimatendarray

在此區域上，積分值的估計。

errorndarray

在此區域上，近似誤差的估計。

註解

此演算法使用類似於 quad_vec 的演算法，後者本身是基於 QUADPACK 的 DQAG* 演算法的實作，實作了全域誤差控制和自適應細分。

用於 Gauss-Kronrod 正交積分的節點和權重來源可以在 [1] 中找到，而用於計算 Genz-Malik cubature 積分中節點和權重的演算法可以在 [2] 中找到。

目前透過 rule 引數支援的規則為

  • "gauss-kronrod"，21 節點 Gauss-Kronrod

  • "genz-malik"，n 節點 Genz-Malik

如果對於 n 維被積函數（其中 n > 2）使用 Gauss-Kronrod，則將透過取 1 維情況下節點的笛卡爾積來找到對應的笛卡爾積規則。 這表示在 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, Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region，《Journal of Computational and Applied Mathematics》，第 6 卷，第 4 期，1980 年，第 295-302 頁，ISSN 0377-0427 DOI:10.1016/0771-050X(80)90039-X

範例

具有向量輸出的 1 維積分:

\[\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 避免奇異點的 1 維積分

\[\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
本頁內容