python - 如何处理结果很小的python中的数值积分？

Question

我遇到了一个问题，即 Python 中的集成返回不正确的积分值和已知解析解。有问题的积分是

LaTex expression for the integral (can't post photos yet)

对于我使用的 sigma 值 (1e-15)，此积分的解的值为 ~ 1.25e-45。然而，当我使用 scipy 集成包来计算这个时，我得到零，我认为这与计算所需的精度有关。

#scipy method
import numpy as np
from scipy.integrate import quad

sigma = 1e-15
f = lambda x: (x**2) * np.exp(-x**2/(2*sigma**2))

#perform the integral and print the result 
solution = quad(f,0,np.inf)[0]
print(solution)
0.0

由于精度是一个问题，我还尝试使用另一个推荐的包 mpmath，它没有返回 0，但与正确答案相差约 7 个数量级。测试较大的 sigma 值会导致解非常接近相应的精确解，但随着 sigma 变小，它似乎变得越来越不正确。

#mpmath method
import mpmath as mp

sigma = 1e-15
f = lambda x: (x**2) * mp.exp(-x**2/(2*sigma**2))

#perform the integral and print the result 
solution = mp.quad(f,[0,np.inf])
print(solution)
2.01359486678988e-52

从这里我可以使用一些建议来获得更准确的答案，因为我希望有信心将 python 集成方法应用于无法解析求解的积分。

Answer 1

您应该为函数添加额外的点作为“中点”，我从 1e-100 到 1 添加了 100 个点以提高准确性。

#mpmath method
import numpy as np
import mpmath as mp

sigma = 1e-15
f = lambda x: (x**2) * mp.exp(-x**2/(2*sigma**2))

#perform the integral and print the result
solution = mp.quad(f,[0,*np.logspace(-100,0,100),np.inf])
print(solution)

1.25286197427129e-45

编辑:原来你需要 10000 个点而不是 100 个点才能得到更准确的结果，1.25331413731554e-45，但计算需要几秒钟。