如何在python中处理结果较小的数值积分

Question

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

积分的LaTex表达式(还不能发布照片(

对于我使用的西格玛值(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个数量级。测试更大的西格玛值会导致解决方案非常接近相应的精确解决方案，但随着西格玛变小，它似乎变得越来越不正确。

#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结果，但计算需要几秒钟。

Answer 2

由于浮点精度的原因，大多数数值积分器都会遇到数字太小的问题。一种解决方案是在计算之前对积分进行缩放。让q->x/sigma，积分变为：

f = lambda q: sigma**3*(q**2) * np.exp(-q**2/2)
solution = quad(f, 0, np.inf)[0]
# solution: 1.2533156529417088e-45