Python中结果极小的数值积分精度问题解决方法咨询

问题分析 #

你要计算的积分是$\int_{0}^{\infty} x^2 e^{-x2/(2\sigma^2)} dx$，解析解为$\frac{\sqrt{\pi}}{2} \sigma^3$。当$\sigma=1e-15$时，解析解约为$1.25e-45$，但：

• scipy.integrate.quad返回0，因为双精度浮点数的精度限制，函数在大部分积分区间内的取值远低于机器可识别的最小非零值，被直接判定为0
• mpmath.quad结果偏差7个数量级，是默认精度不足，加上积分区间的数值特性放大了计算误差

解决方案 #

1. 变量替换简化积分（最优方案） #

做变量替换$t = x/\sigma$，积分可转化为：
$$\int_{0}^{\infty} (\sigma t)^2 e^{-t2/2} \cdot \sigma dt = \sigma^3 \int_{0}^{\infty} t^2 e^{-t2/2} dt$$
后者的积分结果是已知常数$\frac{\sqrt{\pi}}{2} \approx 0.88622692545$，直接计算$\sigma^3$乘以该常数即可得到精确结果：

import numpy as np

sigma = 1e-15
analytical = (np.sqrt(np.pi)/2) * sigma**3
print(analytical)  # 输出约1.2533141373155003e-45

2. 调整mpmath精度适配数值积分 #

如果需要处理无解析解的类似积分，可提高mpmath的工作精度，同时通过变量替换避免直接计算极小值：

import mpmath as mp

sigma = mp.mpf("1e-15")
# 设置100位小数精度，可根据需求调整
mp.mp.dps = 100

# 变量替换后的标准化积分函数
f = lambda t: t**2 * mp.exp(-t**2/2)
# 积分结果乘以sigma^3
solution = sigma**3 * mp.quad(f, [0, mp.inf])
print(solution)  # 输出约1.25331413731550027016490817712e-45

3. scipy的优化思路 #

scipy的quad支持通过epsabs和epsrel参数调整精度，但变量替换仍是这类问题最可靠的处理方式：

import numpy as np
from scipy.integrate import quad

sigma = 1e-15
# 用变量替换后的标准化函数
f = lambda t: t**2 * np.exp(-t**2/2)
# 积分结果乘以sigma^3
solution = sigma**3 * quad(f, 0, np.inf)[0]
print(solution)  # 输出约1.2533141373155003e-45

通用建议 #

对于包含极小参数的积分，优先尝试变量替换将其转化为参数无关的标准形式，从根源上避免数值计算中的下溢或精度丢失。如果必须直接数值积分：

• 使用mpmath时，根据计算需求调高精度（mp.mp.dps）
• 对于scipy，调整精度参数的同时，尽量避免让函数值落入机器精度以下的区间

内容的提问来源于stack exchange，提问作者Adam