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)[source]#

使用双曲正切-正弦积分（tanh-sinh quadrature）对收敛积分进行数值评估。

实际上，双曲正切-正弦积分对于许多被积函数都能实现二次收敛：准确位数大致与函数评估次数呈线性关系[1]。

积分限可以一个或两个都是无穷大，并且端点处的奇点是可接受的。发散积分以及在区间内具有非有限导数或奇点的被积函数超出了其范围，但后者可以通过分别在每个子区间上调用tanhsinh来评估。

参数:

f可调用对象 (callable)

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

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

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

a, b浮点数组类对象 (float array_like)

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

args数组类对象的元组 (tuple of array_like), 可选

传递给f的额外位置数组参数。数组必须彼此以及与a和b的数组可广播。如果所需根的可调用对象需要与x不可广播的参数，则将该可调用对象用f包装，使得f仅接受x和可广播的*args。

log布尔值 (bool), 默认: False

设置为True表示f返回被积函数的对数，并且atol和rtol表示为绝对误差和相对误差的对数。在这种情况下，结果对象将包含积分和误差的对数。这对于数值下溢或溢出可能导致不准确的被积函数很有用。当log=True时，被积函数（f的指数）必须是实数，但它可以是负数，在这种情况下，被积函数的对数是一个虚部为π奇数倍的复数。

maxlevel整型 (int), 默认: 10

算法的最大细化级别。

在第零级，f被调用一次，执行16次函数评估。在每个后续级别，f再被调用一次，函数评估次数大约翻倍。因此，对于许多被积函数，每个连续级别将使结果的准确位数翻倍（直至浮点精度的限制）。

算法将在完成级别maxlevel或满足其他终止条件后终止，以先达到者为准。

minlevel整型 (int), 默认: 2

开始迭代的级别（默认：2）。这不会改变函数评估的总次数或函数被评估的横坐标；它只改变f被调用的次数。如果minlevel=k，则在单次调用中评估级别0到k的所有横坐标处的被积函数。请注意，如果minlevel超过maxlevel，则提供的minlevel将被忽略，并且minlevel将被设置为等于maxlevel。

atol, rtol浮点型 (float), 可选

绝对终止容差（默认：0）和相对终止容差（默认：eps**0.75，其中eps是结果数据类型的精度），分别。如果log为False，则必须是非负且有限的；如果log为True，则必须表示为非负且有限数的对数。当res.error < atol或res.error < res.integral * rtol时，迭代将停止。

preserve_shape布尔值 (bool), 默认: 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可调用对象 (callable), 可选

一个可选的用户提供的函数，在第一次迭代之前和每次迭代之后调用。以callback(res)的形式调用，其中res是一个_RichResult对象，类似于tanhsinh返回的对象（但包含当前迭代的所有变量值）。如果callback引发StopIteration，算法将立即终止，并且tanhsinh将返回一个结果对象。callback不得修改res或其属性。

返回:

res_RichResult

一个类似于scipy.optimize.OptimizeResult实例的对象，具有以下属性。（描述假设值为标量；但是，如果f返回数组，则输出将是相同形状的数组。）

success布尔数组 (bool array)

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

status整型数组 (int array)

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

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

integral浮点数组 (float array)

积分的估计值。

error浮点数组 (float array)

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

maxlevel整型数组 (int array)

使用的最大细化级别。

nfev整型数组 (int array)

f被评估的点的数量。

另请参阅

quad

注意

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

[1]第5节描述了两种误差估计方案：一种试图检测并利用二次收敛；另一种只是简单地比较连续级别的积分估计。虽然两者都没有理论上的严谨性或保守性，但两者在实践中都表现良好。我们的误差估计使用这两种方案中的最小值，其下限为eps * res.integral。

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

参考文献

[1] (1,2,3)

Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. “A comparison of three high-precision quadrature schemes.” Experimental Mathematics 14.3 (2005): 317-329。

[2]

Vanherck, Joren, Bart Sorée, and Wim Magnus. “Tanh-sinh quadrature for single and multiple integration using floating-point arithmetic.” arXiv preprint arXiv:2007.15057 (2020)。

[3]

van Engelen, Robert A. “Improving the Double Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas.” https://www.genivia.com/files/qthsh.pdf

[4]

Takahasi, Hidetosi, and Masatake Mori. “Double exponential formulas for numerical integration.” 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时，函数可以调用参数x，其形状为(4,)或(4, n)，这也是我们观察到的。