如何消除Python三重积分运行时的IntegrationWarning警告？

1. 消除RuntimeWarning（除零操作） #

问题根源 #

• k=0的情况：X = np.arange(0, 50, 0.1)包含k=0，而函数f中有2/k**2项，直接触发除零错误。
• t1/t2=0的情况：t1、t2的积分区间是[-4, 0.1]，包含t=0点，1/(H*t1)和1/(H*t2)在t=0处存在奇点，导致积分计算中出现除零。

修复方案 #

• 调整k的起始值：将X的起始点改为极小正数，避开k=0：

  X = np.arange(0.01, 50, 0.1)
  

• 拆分t的积分区间处理奇点：将t1、t2的积分区间拆分为[-4, -1e-6]和[1e-6, 0.1]，避开t=0的奇点（假设原积分在t=0处收敛），拆分后分别计算积分再求和。

2. 解决IntegrationWarning（无法达到容差） #

问题根源 #

• 积分上限为无穷大，无穷积分的数值计算容易收敛缓慢，难以达到默认精度要求。
• 默认的积分精度参数（epsabs=1.49e-08，epsrel=1.49e-08）过于严格，导致无法满足容差标准。

修复方案 #

• 调整积分精度参数：在tplquad中设置更大的epsabs和epsrel，降低精度要求，比如：

  y, err = integrate.tplquad(..., epsabs=1e-4, epsrel=1e-4)
  

• 变量替换处理无穷积分：将m的积分从[0, inf)转换为[0, 1)，令m = u/(1-u)，则dm = du/(1-u)**2，替换后积分区间变为有限区间，更容易收敛。修改后的函数和积分限如下：

  f_new = lambda t1, t2, u, k: (k**3 * (u/(1-u))**3) * (1/(H * t1)) * (1/(H * t2)) * (2/k)**2 * (1/(1-u)**2)
  y, err = integrate.tplquad(f_new, ti, a, lambda x: ti, lambda x: a, lambda x, y: 0, lambda x, y: 1, args=(val,))
  

完整修复后的代码 #

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

H = 4
ti = -H
a = 0.1

f = lambda t1, t2, m, k: (k**3 * m**3) * (1/(H * t1)) * (1/(H * t2)) * (2/k)**2

# 避开k=0，调整起始值
X = np.arange(0.01, 50, 0.1)

def F(x):
    res = np.zeros_like(x)
    for i, val in enumerate(x):
        # 拆分四个子区间计算积分，避开t=0奇点
        y1, err1 = integrate.tplquad(f, ti, -1e-6, lambda x: ti, lambda x: -1e-6, 
                                     lambda x, y: 0, lambda x, y: float('inf'), 
                                     args=(val,), epsabs=1e-4, epsrel=1e-4)
        y2, err2 = integrate.tplquad(f, ti, -1e-6, lambda x: 1e-6, lambda x: a, 
                                     lambda x, y: 0, lambda x, y: float('inf'), 
                                     args=(val,), epsabs=1e-4, epsrel=1e-4)
        y3, err3 = integrate.tplquad(f, 1e-6, a, lambda x: ti, lambda x: -1e-6, 
                                     lambda x, y: 0, lambda x, y: float('inf'), 
                                     args=(val,), epsabs=1e-4, epsrel=1e-4)
        y4, err4 = integrate.tplquad(f, 1e-6, a, lambda x: 1e-6, lambda x: a, 
                                     lambda x, y: 0, lambda x, y: float('inf'), 
                                     args=(val,), epsabs=1e-4, epsrel=1e-4)
        res[i] = y1 + y2 + y3 + y4
    return res

plt.plot(X, F(X))
plt.title("P(k), H=4")  # 修正原代码标题与H值不一致的问题
plt.xlabel("k")
plt.ylabel("P")
plt.savefig('P(k).png')
plt.show()

额外说明 #

• 如果原始积分在t=0处实际发散，需要重新确认积分区间的合理性，结合物理或数学背景调整区间。
• 变量替换处理无穷积分的方法稳定性更强，建议优先尝试。

内容的提问来源于stack exchange，提问作者Dr. phy