Python数值积分时如何避免小r区间下溢问题以提升计算精度

1. 优先使用带权重适配的自适应积分 #

scipy的quad函数内置支持幂函数乘指数形式的权重函数，完全匹配你g(r)的结构，不需要手动计算极小的g(r)值，从根本上避免数值下溢，同时自动在被积函数变化大的区域加密采样：

from scipy.integrate import quad

# weight='alg'对应权重形式为 x^wvar[0] * exp(-wvar[1] * x)
gamma_weighted, error_est = quad(
    lambda r: some_complicated_function(r, params),
    rmin, rmax,
    weight='alg',
    wvar=(alpha, b)
)
# 乘以g(r)的常数因子得到最终积分结果
gamma = a * gamma_weighted

返回的error_est是积分的官方误差估计，可以直接用来验证结果可靠性。如果积分下限rmin为0且α<1存在奇点，quad也会自动适配处理，不需要额外操作。

2. 辛普森法适配优化 #

如果必须使用辛普森法则，可以做两处优化解决下溢和精度问题：

• 替换均匀采样为对数采样，提升小r区域的采样密度，避免漏过小r段的被积函数变化
• 引入缩放因子避免小r区域数值下溢到机器精度以下，积分完成后再还原结果

import numpy as np
from scipy.integrate import simps

# 对数采样，小r区域点更密
r = np.logspace(np.log10(rmin), np.log10(rmax), 5000)
f_val = some_complicated_function(r, params)
# 缩放因子根据g(r)最小值调整，保证缩放后最小值远大于机器epsilon（~1e-16）
SCALE = 1e30
# 用numpy向量化操作替代列表推导，效率和稳定性更高
g_val = a * r**alpha * np.exp(-b*r) * SCALE
gamma = simps(f_val * g_val, r) / SCALE

3. 积分区间拆分（极端小r场景适用） #

如果小r区域g(r)量级极低，但f(r)在小r段的变化规律可近似，可以拆分积分区间处理：

• 选择截断值r_cut，保证r >= r_cut时g(r)不会出现下溢
• [r_cut, rmax]段用上述常规数值积分方法计算
• [rmin, r_cut]段对f(r)做低阶泰勒展开后，结合g(r)的解析积分公式计算，误差可控

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