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可调用对象

要积分的函数。签名必须是

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

其中 xi 的每个元素都是一个有限实数，argsi 是一个元组，其中可能包含任意数量的可与 xi 广播的数组。f 必须是一个逐元素函数：有关详细信息，请参阅参数 preserve_shape 的文档。它不能修改数组 xi 或 argsi 中的数组。如果在任一端点求值时，f 返回一个具有复数 dtype 的值，则后续参数 x 将具有复数 dtype（但虚部为零）。

a, bfloat 数组式对象

积分的实数下限和上限。必须彼此以及 args 中的数组进行广播。元素可以是无穷大。

args数组式对象的元组，可选

要传递给 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 数组

当算法成功终止时（状态为 0）为 True。否则为 False。

statusint 数组

一个表示算法退出状态的整数。

0 : 算法收敛到指定的容差。-1 : （未使用）-2 : 达到最大迭代次数。-3 : 遇到非有限值。-4 : 迭代由 callback 终止。1 : 算法正常进行（仅在 callback 中）。

integralfloat 数组

积分的估计值。

errorfloat 数组

误差的估计值。仅在完成第二级或更高级别时可用；否则为 NaN。

maxlevel整数数组

使用的最大细化级别。

nfev整数数组

评估 f 的点的数量。

另请参阅

quad

注释

实现 [1] 中描述的算法，并针对有限精度算术进行了一些小的调整，包括 [2] 和 [3] 中描述的一些调整。tanh-sinh 方案最初在 [4] 中引入。

由于横坐标中的浮点误差，在迭代过程中可能会在区间的端点评估函数，但将忽略端点处函数返回的值。

参考文献

[1] (1,2,3)

Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. “三种高精度求积方案的比较。” Experimental Mathematics 14.3 (2005): 317-329.

[2]

Vanherck, Joren, Bart Sorée, and 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, and 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,)，对应于在单个横坐标和所有四个频率下评估被积函数；这用于输入验证并确定存储结果的数组的大小和数据类型。下一个形状对应于在横坐标的初始网格和所有四个频率下评估被积函数。对函数的连续调用使评估函数的横坐标总数翻倍。但是，在稍后的函数评估中，被积函数在较少的频率下评估，因为相应的积分已经收敛到所需的容差。这节省了函数评估以提高性能，但它要求函数接受任何形状的参数。

“矢量值”被积函数，例如那些为与 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 调用该函数，这就是我们观察到的。

返回在新标签页中打开