scipy.integrate 的意外行为答案

【问题标题】：Unexpected behaviour of scipy.integratescipy.integrate 的意外行为
【发布时间】：2019-08-13 12:26:22
【问题描述】：

我想使用 scipy.integrate 进行一些数值计算。我只是跑了一个小例子来尝试它并遇到了一些意想不到的行为。

我编写了一些干净的代码来演示这个问题。我使用一个非常简单的指数分布来测试。

这是我的代码：

import numpy as np
import sys
import scipy as sc
from scipy import integrate


print(sys.version)
print(np.version.version)
print(sc.version.version)
print()

r1 = integrate.quad(lambda x: sc.exp(-x), 0, 10)
r2 = integrate.quad(lambda x: sc.exp(-x), 0, 100000)
r3 = integrate.quad(lambda x: sc.exp(-x), 0, np.inf)

print(r1)
print(r2)
print(r3)

print()
r4 = integrate.quad(lambda x: sc.exp(-x), 0, 10000)
print(r4)

输出是

3.7.2 (default, Jan  2 2019, 17:07:39) [MSC v.1915 64 bit (AMD64)]
1.15.4
1.1.0

r1 (0.9999546000702375, 2.8326146575791917e-14)
r2 (2.0614532085314573e-45, 4.098798466247153e-45)
r3 (1.0000000000000002, 5.842606996763696e-11)

r4 (1.0, 1.6059202674761255e-14)

我希望所有输出总是大约为 1。但是在 r2 中，我得到了一个非常小的值。奇怪的是，当积分到无穷大 (r3) 或非常小的边界 (r1) 时，问题并没有出现。此外，通过将限制降低一个数量级 (r4)，我也得到了完美的结果。

有谁知道为什么这个问题会出现在 scipy 中？我会称之为错误，但也许我违反了一些限制？我如何提前知道以防止在我应用的问题中出现错误的结果？

提前谢谢你

full_output 的输出：

r2 (2.0614532085314573e-45, 4.098798466247153e-45, {'neval': 63, 'last': 2, 'iord': array([      1,       2,       3,       4,       5, 6357060, 6357108,
       4259932, 6357102, 7274595, 6553710, 3342433, 7077980, 6422633,
       7536732, 7602281, 2949221, 6357104, 7012451, 6750305, 7536741,
       7536732, 6881379, 7929968, 7274588, 7602288, 7143529, 7995497,
       6029413, 7209055, 7077998, 3014771, 7340131, 3604531, 7798829,
       7209065, 6357087, 6553709, 3407926, 7340078, 6553721, 3276846,
       5046318, 7209057, 6684777, 7536741,     116, 6619136, 7602291,
             0], dtype=int32), 'alist': array([0.00000000e+000, 5.00000000e+004, 0.00000000e+000, 0.00000000e+000,
       6.88436472e-272, 3.80218509e-136, 2.65902947e-068, 2.20016853e-034,
       1.04474528e-019, 3.09734336e-016, 9.03970673e-019, 8.23342652e-316,
       8.23342968e-316, 8.23343284e-316, 8.23343601e-316, 8.23343917e-316,
       8.23344233e-316, 8.23344549e-316, 8.23344865e-316, 8.23345182e-316,
       8.23345498e-316, 8.23345814e-316, 8.23346130e-316, 8.23346446e-316,
       8.23346763e-316, 8.23347079e-316, 8.23347395e-316, 8.23347711e-316,
       8.23348027e-316, 8.23348344e-316, 8.23348660e-316, 8.23348976e-316,
       8.23349292e-316, 8.23349608e-316, 8.23349925e-316, 8.23350241e-316,
       8.23350557e-316, 8.23350873e-316, 8.23351189e-316, 8.23351506e-316,
       8.23351822e-316, 8.23352138e-316, 8.23352454e-316, 8.23352770e-316,
       8.23353087e-316, 8.23353403e-316, 8.23353719e-316, 8.23354035e-316,
       8.23354351e-316, 8.23354668e-316]), 'blist': array([5.00000000e+004, 1.00000000e+005, 0.00000000e+000, 0.00000000e+000,
       6.88436472e-272, 3.80218509e-136, 2.65902947e-068, 2.20016853e-034,
       1.04474528e-019, 3.09734336e-016, 9.03970673e-019, 1.20736675e+285,
       1.05117823e-153, 1.05132391e-153, 1.05146958e-153, 3.79823888e-258,
       1.61465766e+184, 3.11517960e+161, 4.26137323e+257, 6.01346953e-154,
       6.01366349e-154, 1.19632546e-153, 3.64465882e-086, 1.31100174e-259,
       1.20679441e-153, 1.20679327e-153, 3.24245662e-086, 3.64465882e-086,
       6.01357764e-154, 1.20679441e-153, 5.75105581e+072, 2.20791354e+214,
       1.27734658e-152, 5.29444423e+160, 6.19633416e+223, 2.25563599e-153,
       8.21947530e+223, 6.09892510e-013, 1.06097757e-153, 2.86747940e-110,
       6.06154135e-154, 6.06445477e-154, 6.96312298e-077, 3.00226946e-067,
       6.03810921e-154, 1.30421760e-076, 1.21438942e-067, 4.61448322e-072,
       8.51221910e-053, 3.73237334e+069]), 'rlist': array([2.06145321e-045, 0.00000000e+000, 6.73898103e+149, 3.51023756e+151,
       4.50937881e-292, 9.43293441e-314, 4.65203811e+151, 6.99386802e-283,
       3.53886392e-308, 1.33360313e+241, 1.15420781e+171, 9.30281767e+242,
       1.17364463e+214, 3.12671297e+185, 2.85341794e-313, 8.18432962e-085,
       6.45840689e+170, 4.42638830e-239, 9.78681729e+199, 3.38460675e+125,
       3.11732880e+150, 9.78747303e+199, 2.27948172e-191, 1.04972250e+214,
       4.77402433e+180, 1.12985581e+277, 3.16464606e-307, 1.33360315e+241,
       1.76252970e-310, 1.02318154e-012, 1.15549302e-313, 1.03539814e-308,
       1.33360293e+241, 5.67421675e-311, 5.00120719e-162, 6.46048250e-313,
       1.68400738e-019, 1.10811151e-302, 1.66468912e-312, 1.09403545e-303,
       1.27613271e-303, 7.10020498e-270, 4.99875566e-111, 9.11927054e-304,
       9.11571045e-304, 9.11749048e-304, 9.11571042e-304, 9.60205653e+303,
       5.43239349e-312, 1.79972786e-304]), 'elist': array([4.09879847e-045, 0.00000000e+000, 6.47287707e+170, 5.98178835e-154,
       1.69375668e+190, 4.44389806e+252, 1.12297399e+219, 1.87673453e-152,
       7.20706153e+159, 1.27826731e-152, 2.43812981e-152, 5.52716101e+228,
       6.01346953e-154, 1.57761457e+214, 7.19938459e+252, 3.94357072e+180,
       3.44210870e+175, 3.62478142e+228, 1.23732543e-259, 3.53810655e+155,
       4.81222029e+233, 1.06843264e-258, 9.15000112e+199, 4.26614628e+180,
       3.53387914e+246, 2.35509149e+251, 1.69375944e+190, 1.57762309e+214,
       6.19634286e+223, 8.95533289e-106, 5.98148090e-154, 1.17914189e+195,
       5.42869734e+213, 6.72794695e+199, 5.30383390e+180, 1.02188594e-152,
       2.16452413e+233, 7.50052033e+247, 6.98907523e+096, 7.69843824e+218,
       3.23097122e+174, 9.84214185e-154, 1.36723829e+161, 1.19346501e+243,
       1.94670285e+227, 2.21366476e+214, 8.95533289e-106, 8.75378213e+247,
       1.87673453e-152, 2.50722129e-310])})

´´´

【问题讨论】：

您能否尝试将关键字参数full_output=1 用于r2。这将在输出中添加一个字典，其中包含一些有关计算的额外信息，请参阅here。我的第一个猜测是您在计算中遇到了“反甜点”。
@Chgad 我在上面添加了输出
分区子区间的结果在键rlist的数组中。积分算法选择的第一个区间是0.0到5.0e+4，结果是2.06e-45和一个两倍的绝对误差（参见elist 的第一个元素）。这根本不正确。进一步比较区间会发现更多这样不切实际的结果。我建议玩弄积分的上限。也许选择了 99900 之类的东西？
仅供参考：这是关于quad 的常见问题解答。 SO上有很多这样的相关问题；这里只有三个我很容易找到，因为我提供了答案：stackoverflow.com/questions/29179778/…、stackoverflow.com/questions/47193045/…、stackoverflow.com/questions/40536467/…

标签： python numpy scipy numeric

【解决方案1】：

这不是一个错误，它与积分的数值精度有关，并且您正在积分一个在大部分时间间隔内（几乎）为 0 的函数。来自docs：

请注意，与使用积分区间的大小可能无法正确积分这个方法。

根据您的输出，该函数仅使用两个 (last=2) 间隔，在每个间隔上评估 rlist=(2.06145321e-045, 0.00000000e+000, ..) 的值（有关输出的更多详细信息，请参阅文档）

您可以在区间中添加点以强制例程使用更接近左边界的点。

a = quad(lambda x: np.exp(-x), 0, 1e9, points=np.logspace(-10,3,10))
print(a)
(0.9999999999999997, 2.247900608926337e-09)

补充说明（感谢@norok2）：请注意，points 是有界积分区间中的一系列断点，其中可能出现被积函数的局部困难（例如，奇异点、不连续点）。在这种情况下，我不是用它来指出不连续性，而是强制quad 在左边界附近执行更多积分步骤，使用对数间隔，因为我有一个指数函数（这当然是任意和对于这个函数，因为我知道它的形状）。

【讨论】：

非常感谢您详细解释的回答。我想知道为什么 +infinity 集成效果更好。在这种情况下，它会自动进行对数区间评估吗？
10000 的限制给出了很好的结果。是否有任何经验法则可以使用什么限制来获得良好的精度和发生上述事情的风险很小？
@Mstain 您可以更好地解释points 参数。 points 是有界积分区间中的一系列断点，其中可能出现被积函数的局部困难（例如，奇异点、不连续点）。序列不必排序。因此，您不需要那么多间隔，选择它们的方法是确保捕获其中的困难。例如，在这种情况下，只要np.exp(-b) 足够小足够，points=(b, np.inf) 就会同样有效。对于任意函数，没有简单的方法可以预先知道这一点。
（points 中的值越多，计算负担就越大）。另一种方法是自己执行采样，然后使用scipy.integrate.trapz() 或类似方法进行数值积分。
没错。另见this post

【解决方案2】：

无需在（非常大的）区间内将积分转换为积分。形式的积分有一个特定的积分方案

即Gauss-Laguerre quadrature。它也包含在quadpy（我的一个项目）中。试试看

import numpy
import quadpy

scheme = quadpy.e1r.gauss_laguerre(1)

val = scheme.integrate(lambda x: numpy.ones(x.shape[1:]))

print(val)

1.0

【讨论】：