Scipy integrate.quad积分结果不准确，求SciPy/Numpy可靠积分方法

问题重现 #

你遇到的问题是scipy.integrate.quad在处理包含极大上限的积分时，因浮点数精度限制返回错误结果：

from scipy.integrate import quad
xmax=1.0e23
xmin=0.5
def myFUN(x):
    return 1/(1+x)**(3.0)
print( quad( myFUN, xmin,xmax )[0] )

返回结果：

2.3171047961935996e-40

而解析解计算：

def ana(x):
    return -0.5/(1+x)**2.0 
print( ana(xmax)-ana(xmin) )

得到正确结果：

0.2222222222222222

问题原因 #

当x远大于1时（如1e23），双精度浮点数无法区分x和x+1——双精度仅保留约16位有效数字，1e23的量级远大于1，导致1+x被直接近似为x，进而1/(1+x)**3计算为0，积分结果自然趋近于0。

解决方案 #

针对无法手动参数化的复杂积分，可按优先级尝试以下方法：

1. 优化quad的参数设置

通过参数引导积分器聚焦于函数贡献显著的区间：

• 添加采样点：用points参数在函数变化剧烈的区域（如x=1附近）设置采样点，避免积分器在无贡献的大数值区间浪费计算：

  print(quad(myFUN, xmin, xmax, points=[1.0, 10.0, 100.0, 1e6])[0])
  

• 调整迭代次数与误差阈值：增大limit参数提升迭代次数，同时降低epsabs/epsrel设置更严格的误差控制：

  print(quad(myFUN, xmin, xmax, limit=1000, epsabs=1e-12, epsrel=1e-12)[0])
  

2. 使用高精度数值积分

借助mpmath库的任意精度积分功能，避免大数值下的浮点数精度损失：

import mpmath as mp
mp.mp.dps = 50  # 设置50位小数精度
xmax = mp.mpf("1e23")
xmin = mp.mpf("0.5")
def myFUN(x):
    return 1/(1+x)**3
print(mp.quad(myFUN, [xmin, xmax]))

该方法会返回精确结果：0.22222222222222222222222222222222222222222222222222

3. 拆分积分区间

手动将积分拆分为贡献显著区间和可忽略区间：

• 找到临界值X（如双精度下X≈1e16，超过后x+1无法与x区分）
• 仅计算[xmin, X]区间的积分，[X, xmax]区间因函数值趋近于0，贡献可忽略

总结 #

对于复杂积分，优先尝试优化quad的参数设置；若仍存在精度问题，切换到mpmath的高精度积分是最可靠的通用方案。

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