大数的scipy.integrate.quad精度

大数的scipy.integrate.quad精度

我认为问题是由于增加

np.exp(-x)

很快就变得很小

x

，由于有限的数值精度而导致评估为零。例如，即使

x

小到

x=10**2*

，

np.exp(-x)

其结果也为

3.72007597602e-44

，而

x

order

10**3

或更高的值将得出

0

。

我不知道的实现细节

quad

，但是它可能会在给定的集成范围内对要集成的功能进行某种采样。对于较大的积分上限，大多数

np.exp(-x)

评估样本为零，因此积分值被低估了。（请注意，在这些情况下，所提供的绝对误差

quad

与积分值的阶次相同，这表明后者不可靠。）

避免此问题的一种方法是将积分上限限制为一个数值，高于该数值数值函数会变得非常小（因此对积分值的贡献很小）。从您的代码片段来看，值

10**4

似乎是一个不错的选择，但是，值

10**2

也会导致对积分的准确评估。

避免数值精度问题的另一种方法是使用以 任意
精度算术执行计算的模块，例如

mpmath

。例如，对于

x=10**5

，

mpmath

评估

exp(-x)

如下（使用本机

mpmath

指数函数）

import mpmath as mpprint(mp.exp(-10**5))

3.56294956530937e-43430

请注意此值有多小。在标准硬件数值精度（用于

numpy

）下，该值变为

0

。

mpmath

提供了一个积分函数（

mp.quad

），该函数可以为积分上限的任意值提供准确的积分估计。

import mpmath as mpprint(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))

1.00.9999996504694740.9999999999965160.999999999999997

通过将精度提高到

50

小数点（

15

这是标准精度），我们还可以获得更准确的估算值

mp.mp.dps = 50;print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))

1.00.999999999999999999999999999999999999999998298802620.999999999999999999999999999999999999999999999974630.99999999999999999999999999999999999999999999999998

通常，获得此精度的成本是增加的计算时间。

PS：毋庸置疑，如果您首先能够以分析方式评估您的积分（例如，借助

Sympy

），您会忘记上述所有内容。