大M值（超5000时项达10^100）下无穷积分的Python求解问题

我需要计算一个从负无穷到正无穷的积分，积分式中包含一个以M为指数的项。当M取5000及以上的较大值时，该项的数值会超过10^100。我编写了Python 3.10代码，在M=1000时可以正常运行，但M更大时无法正常工作。

代码如下：

import numpy as np
from scipy import integrate, special

def integrand(u, M, Z, r):
    sqrt_r = np.sqrt(r)
    sqrt_term = np.sqrt(2 * (1 - r))
    arg1 = (Z + sqrt_r * u) / sqrt_term
    arg2 = (-Z + sqrt_r * u) / sqrt_term
    erf_term = special.erf(arg1) - special.erf(arg2)
    return np.exp(-u**2 / 2) / np.sqrt(2 * np.pi) * erf_term**M

def evaluate_integral(M, Z, r):
    result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))
    return 1 / (2**M) * result

# Example usage:
M = 1000
Z = 5
r = 0.6

integral_value = evaluate_integral(M, Z, r)

运行时出现如下警告：

IntegrationWarning: 检测到舍入误差，无法达到要求的容差。误差可能被低估。
  result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))

且integral_value的结果为nan。请问能否在Python中完成该积分计算，还是必须更换其他工具？

完全可以在Python中完成计算，问题核心是大数溢出和数值不稳定，可以通过以下几种方法解决：

1. 对数转换避免溢出 #

直接计算erf_term**M和2**M会因M过大导致数值溢出，改用对数计算后再取指数，将乘法转加法、除法转减法，从根源上避免溢出：

import numpy as np
from scipy import integrate, special

def integrand_log(u, M, Z, r):
    sqrt_r = np.sqrt(r)
    sqrt_term = np.sqrt(2 * (1 - r))
    arg1 = (Z + sqrt_r * u) / sqrt_term
    arg2 = (-Z + sqrt_r * u) / sqrt_term
    erf_term = special.erf(arg1) - special.erf(arg2)
    # 处理erf_term为0的极端情况，避免对数报错
    log_erf = np.log(erf_term) if erf_term > 1e-300 else -np.inf
    log_gaussian = (-u**2 / 2) - np.log(np.sqrt(2 * np.pi))
    # 转成指数形式返回，参与积分计算
    return np.exp(log_gaussian + M * log_erf)

def evaluate_integral(M, Z, r):
    result, _ = integrate.quad(integrand_log, -np.inf, np.inf, args=(M, Z, r))
    # 1/(2^M)等价于exp(-M*ln2)，同样用对数转换避免溢出
    return result * np.exp(-M * np.log(2))

2. 缩小积分区间，利用函数峰值特性 #

当M很大时，erf_term^M只有在erf_term接近1时才不会趋近于0，结合高斯分布的衰减特性，被积函数的有效贡献区间非常有限。可以把积分区间从(-inf, inf)缩小到有限区间（比如[-10, 10]），区间外的被积函数值可忽略，能大幅减少积分误差：

# 修改integrate.quad的区间参数
result, _ = integrate.quad(integrand_log, -10, 10, args=(M, Z, r))

也可以通过points参数指定函数峰值附近的关键点，帮助积分器更精准地处理核心区域：

result, _ = integrate.quad(integrand_log, -10, 10, args=(M, Z, r), points=[0])

3. 使用高精度计算库 #

如果对数转换后仍存在精度问题，可以用mpmath库进行任意精度计算，它支持高精度的积分和特殊函数运算，彻底解决浮点数精度限制：

import mpmath as mp

mp.mp.dps = 50  # 设置计算精度为50位小数

def integrand_mpmath(u, M, Z, r):
    sqrt_r = mp.sqrt(r)
    sqrt_term = mp.sqrt(2 * (1 - r))
    arg1 = (Z + sqrt_r * u) / sqrt_term
    arg2 = (-Z + sqrt_r * u) / sqrt_term
    erf_term = mp.erf(arg1) - mp.erf(arg2)
    log_erf = mp.log(erf_term)
    log_gaussian = (-u**2 / 2) - mp.log(mp.sqrt(2 * mp.pi))
    return mp.e ** (log_gaussian + M * log_erf)

def evaluate_integral_mpmath(M, Z, r):
    result = mp.quad(lambda u: integrand_mpmath(u, M, Z, r), [-mp.inf, mp.inf])
    return result / (2**M)

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