scipy.integrate.

tanhsinh

scipy.integrate.tanhsinh(f, a, b, *, args=(), log=False, maxlevel=None, minlevel=2, atol=None, rtol=None, preserve_shape=False, callback=None)

使用 tanh-sinh 正交數值評估收斂積分。

實際上，對於許多被積函數，tanh-sinh 正交可以實現二次收斂：準確位數 的數量大致與函數評估的次數成線性比例 [1]。

積分的上限或下限可以是無限的，端點的奇異點也是可以接受的。發散積分和在區間內具有非有限導數或奇異點的被積函數超出範圍，但後者可以透過在每個子區間上單獨呼叫 tanhsinh 來評估。

參數:

f可呼叫物件 (callable)

要積分的函數。簽名必須是

f(xi: ndarray, *argsi) -> ndarray

其中 xi 的每個元素都是有限實數，而 argsi 是一個元組，可以包含任意數量的可與 xi 廣播的陣列。f 必須是元素級函數：請參閱參數 preserve_shape 的文件以瞭解詳細資訊。它不得變更陣列 xi 或 argsi 中的陣列。如果 f 在端點評估時傳回具有複數 dtype 的值，則後續引數 x 將具有複數 dtype（但虛部為零）。

a, bfloat array_like

積分的實數下限和上限。必須可與彼此以及 args 中的陣列廣播。元素可以是無限的。

argsarray_like 元組，選用

要傳遞給 f 的其他位置陣列引數。陣列必須可與彼此以及 a 和 b 的陣列廣播。如果需要尋找根的可呼叫物件需要無法與 x 廣播的引數，請使用 f 包裝該可呼叫物件，以便 f 僅接受 x 和可廣播的 *args。

logbool，預設值：False

設定為 True 表示 f 傳回被積函數的對數，並且 atol 和 rtol 表示為絕對誤差和相對誤差的對數。在這種情況下，結果物件將包含積分和誤差的對數。這對於數值下溢或溢位會導致不準確的被積函數很有用。當 log=True 時，被積函數（f 的指數）必須是實數，但它可以是負數，在這種情況下，被積函數的對數是一個複數，其虛部是 π 的奇數倍。

maxlevelint，預設值：10

演算法的最大精細化層級。

在第零層級，f 被呼叫一次，執行 16 次函數評估。在每個後續層級，f 再被呼叫一次，大約使已執行的函數評估次數加倍。因此，對於許多被積函數，每個連續層級都會使結果中的準確位數加倍（達到浮點精度的限制）。

演算法將在完成層級 maxlevel 或滿足另一個終止條件後終止，以先到者為準。

minlevelint，預設值：2

開始迭代的層級（預設值：2）。這不會變更函數評估的總次數或函數評估的橫座標；它僅變更呼叫 f 的次數。如果 minlevel=k，則在單次呼叫中評估層級 0 到 k 的所有橫座標處的被積函數。請注意，如果 minlevel 超過 maxlevel，則提供的 minlevel 會被忽略，並且 minlevel 設定為等於 maxlevel。

atol, rtolfloat，選用

絕對終止容差（預設值：0）和相對終止容差（預設值：eps**0.75，其中 eps 是結果 dtype 的精度）。當 res.error < atol + rtol * abs(res.df) 時，迭代將停止。誤差估計如 [1] 第 5 節中所述。雖然在理論上不嚴謹或保守，但據說在實務中效果良好。如果 log 為 False，則必須為非負數且有限，如果 log 為 True，則必須表示為非負數且有限數的對數。

preserve_shapebool，預設值：False

在下文中，「f 的引數」指的是陣列 xi 和 argsi 中的任何陣列。讓 shape 為 a、b 和 args 的所有元素（在概念上與傳遞到 f 的 xi` 和 ``argsi 不同）的廣播形狀。

• 當 preserve_shape=False（預設值）時，f 必須接受任何 可廣播形狀的引數。

• 當 preserve_shape=True 時，f 必須接受形狀為 shape 或 shape + (n,) 的引數，其中 (n,) 是正在評估函數的橫座標數。

在任一種情況下，對於 xi 中的每個純量元素 xi[j]，f 傳回的陣列都必須在相同索引處包含純量 f(xi[j])。因此，輸出的形狀始終是輸入 xi 的形狀。

請參閱範例。

callback可呼叫物件，選用

要在第一次迭代之前和每次迭代之後呼叫的可選使用者提供函數。呼叫為 callback(res)，其中 res 是類似於 _differentiate 傳回的 _RichResult（但包含所有變數的目前迭代值）。如果 callback 引發 StopIteration，則演算法將立即終止，並且 tanhsinh 將傳回結果物件。callback 不得變更 res 或其屬性。

傳回值:

res_RichResult

類似於 scipy.optimize.OptimizeResult 實例的物件，具有以下屬性。（描述的撰寫方式就好像值將是純量；但是，如果 f 傳回陣列，則輸出將是相同形狀的陣列。）

successbool 陣列

當演算法成功終止時為 True（狀態 0）。否則為 False。

statusint 陣列

表示演算法結束狀態的整數。

0 : 演算法收斂到指定的容差。-1 : (未使用) -2 : 已達到最大迭代次數。-3 : 遇到非有限值。-4 : 迭代被 callback 終止。1 : 演算法正常進行中（僅在 callback 中）。

integralfloat 陣列

積分的估計值。

errorfloat 陣列

誤差的估計值。僅在完成第二層或更高層級時可用；否則為 NaN。

maxlevelint 陣列

使用的最大精細化層級。

nfevint 陣列

評估 f 的點數。

另請參閱

quad

註解

實作 [1] 中描述的演算法，並針對有限精度算術進行了細微調整，包括 [2] 和 [3] 中描述的一些調整。tanh-sinh 方案最初在 [4] 中引入。

由於橫座標中的浮點誤差，函數可能會在迭代期間在區間的端點進行評估，但函數在端點傳回的值將被忽略。

參考文獻

[1] (1,2,3)

Bailey, David H.、Karthik Jeyabalan 和 Xiaoye S. Li。「三種高精度正交方案的比較。」Experimental Mathematics 14.3 (2005): 317-329。

[2]

Vanherck, Joren、Bart Sorée 和 Wim Magnus。「使用浮點算術進行單次和多次積分的 Tanh-sinh 正交。」arXiv preprint arXiv:2007.15057 (2020)。

[3]

van Engelen, Robert A. 「改進雙指數正交 Tanh-Sinh、Sinh-Sinh 和 Exp-Sinh 公式。」 https://www.genivia.com/files/qthsh.pdf

[4]

Takahasi, Hidetosi 和 Masatake Mori。「數值積分的雙指數公式。」Publications of the Research Institute for Mathematical Sciences 9.3 (1974): 721-741。

範例

在您的瀏覽器中試試看！

評估高斯積分

>>> import numpy as np
>>> from scipy.integrate import tanhsinh
>>> def f(x):
...     return np.exp(-x**2)
>>> res = tanhsinh(f, -np.inf, np.inf)
>>> res.integral  # true value is np.sqrt(np.pi), 1.7724538509055159
1.7724538509055159
>>> res.error  # actual error is 0
4.0007963937534104e-16

對於遠離零的引數，高斯函數（鐘形曲線）的值幾乎為零，因此有限區間上的積分值幾乎相同。

>>> tanhsinh(f, -20, 20).integral
1.772453850905518

但是，在不利的積分限制下，積分方案可能無法找到重要區域。

>>> tanhsinh(f, -np.inf, 1000).integral
4.500490856616431

在這種情況下，或者當區間內存在奇異點時，將積分分成幾個部分，端點位於重要點。

>>> tanhsinh(f, -np.inf, 0).integral + tanhsinh(f, 0, 1000).integral
1.772453850905404

對於涉及非常大或非常小量級的積分，請使用對數積分。（為了說明目的，以下範例顯示了常規積分和對數積分都有效的情況，但對於更極端的積分限制，對數積分將避免在正常評估積分時遇到的下溢。）

>>> res = tanhsinh(f, 20, 30, rtol=1e-10)
>>> res.integral, res.error
(4.7819613911309014e-176, 4.670364401645202e-187)
>>> def log_f(x):
...     return -x**2
>>> res = tanhsinh(log_f, 20, 30, log=True, rtol=np.log(1e-10))
>>> np.exp(res.integral), np.exp(res.error)
(4.7819613911306924e-176, 4.670364401645093e-187)

積分的限制和 args 的元素可以是可廣播的陣列，並且逐元素執行積分。

>>> from scipy import stats
>>> dist = stats.gausshyper(13.8, 3.12, 2.51, 5.18)
>>> a, b = dist.support()
>>> x = np.linspace(a, b, 100)
>>> res = tanhsinh(dist.pdf, a, x)
>>> ref = dist.cdf(x)
>>> np.allclose(res.integral, ref)
True

預設情況下，preserve_shape 為 False，因此可呼叫物件 f 可以使用任何可廣播形狀的陣列呼叫。例如

>>> shapes = []
>>> def f(x, c):
...    shape = np.broadcast_shapes(x.shape, c.shape)
...    shapes.append(shape)
...    return np.sin(c*x)
>>>
>>> c = [1, 10, 30, 100]
>>> res = tanhsinh(f, 0, 1, args=(c,), minlevel=1)
>>> shapes
[(4,), (4, 34), (4, 32), (3, 64), (2, 128), (1, 256)]

為了瞭解這些形狀的來源 - 並更好地理解 tanhsinh 如何計算準確的結果 - 請注意，c 的值越高，對應的正弦波頻率越高。頻率越高的正弦波使被積函數更複雜，因此需要更多的函數評估才能達到目標精度

>>> res.nfev
array([ 67, 131, 259, 515], dtype=int32)

初始 shape (4,) 對應於在單個橫座標和所有四個頻率下評估被積函數；這用於輸入驗證並確定儲存結果的陣列的大小和 dtype。下一個形狀對應於在橫座標的初始網格和所有四個頻率下評估被積函數。對函數的後續呼叫使已評估函數的橫座標總數加倍。但是，在後續函數評估中，被積函數在較少的頻率下評估，因為對應的積分已收斂到所需的容差。這節省了函數評估以提高效能，但它要求函數接受任何形狀的引數。

「向量值」被積函數，例如那些為與 scipy.integrate.quad_vec 一起使用而編寫的被積函數，不太可能滿足此要求。例如，考慮

>>> def f(x):
...    return [x, np.sin(10*x), np.cos(30*x), x*np.sin(100*x)**2]

此被積函數與 tanhsinh 不相容；例如，輸出的形狀將與 x 的形狀不同。這樣的函數可以透過引入其他參數來轉換為相容形式，但這會很不方便。在這種情況下，更簡單的解決方案是使用 preserve_shape。

>>> shapes = []
>>> def f(x):
...     shapes.append(x.shape)
...     x0, x1, x2, x3 = x
...     return [x0, np.sin(10*x1), np.cos(30*x2), x3*np.sin(100*x3)]
>>>
>>> a = np.zeros(4)
>>> res = tanhsinh(f, a, 1, preserve_shape=True)
>>> shapes
[(4,), (4, 66), (4, 64), (4, 128), (4, 256)]

在這裡，a 和 b 的廣播形狀為 (4,)。使用 preserve_shape=True，可以使用形狀為 (4,) 或 (4, n) 的引數 x 呼叫函數，這就是我們觀察到的。

返回在新分頁中開啟