如何用tplquad/nquad积分不连续函数？解决积分收敛报错

我在对函数f(t1,t2,m)执行三重积分时，遇到多条IntegrationWarning错误：积分可能发散或收敛缓慢、已达最大细分次数(50)、算法不收敛且检测到舍入误差。想了解如何在tplquad或nquad中，像quad那样通过设置epsrel、epsabs、points参数，或使用if-else条件判断来完成积分。

我的代码如下：

from scipy import integrate 
import math
import scipy

k=3;a=0.1;nab=10**-9;n=0.94;kz= 0.002;
# 定义被积函数
def f(m,t1,t2):
    return (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + m*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*m**2)*t1**2 + 8*(3 + 2*m**2)*t1**3 + (6 + 4*m**2)*t1**4 + 6*t1**5 + 3*t1**6))) * ((50 * t1**2) +5* (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + (k - m)*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*(k - m)**2)*t1**2 + 8*(3 + 2*(k - m)**2)*t1**3 + (6 + 4*(k - m)**2)*t1**4 + 6*t1**5 + 3*t1**6)))) 


# 定义积分区间
a = 0.1
b = 5


# 执行积分
result, error = integrate.tplquad(f, 0.1, 2, lambda x: 0.1, lambda x: 2, lambda x,y: 0.1, lambda x,y: 2 )

print("result:", result)

解决方法 #

1. 传递精度控制参数（epsrel、epsabs） #

tplquad底层调用quad，可直接在函数末尾添加epsrel（相对误差阈值）和epsabs（绝对误差阈值）参数，调整收敛要求：

# 修改积分调用，设置精度参数
result, error = integrate.tplquad(
    f, 
    0.1, 2, 
    lambda x: 0.1, lambda x: 2, 
    lambda x,y: 0.1, lambda x,y: 2,
    epsrel=1e-4,  # 相对误差，可按需调整
    epsabs=1e-6   # 绝对误差
)

2. 指定奇点/断点（points参数） #

若积分区间内存在函数突变、奇点，可给内层积分传递points参数，引导积分器在这些点附近细分。nquad的参数更灵活，支持给每个维度单独设置：

# 用nquad实现，给每个维度指定参数
def limits():
    # 顺序：m, t1, t2的积分区间
    return [
        [0.1, 2],
        [0.1, 2],
        [0.1, 2]
    ]

result, error = integrate.nquad(
    f, 
    limits,
    opts=[
        {'points': [1.0]},  # m的奇点/断点
        {'epsrel': 1e-4},   # t1的精度参数
        {}                  # t2使用默认参数
    ]
)

3. 函数内添加条件判断处理数值异常 #

检查分母是否接近0，避免除以0导致的数值爆炸，返回极小值跳过异常点：

def f(m,t1,t2):
    # 先计算两个分母
    denom1 = 4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*m**2)*t1**2 + 8*(3 + 2*m**2)*t1**3 + (6 + 4*m**2)*t1**4 + 6*t1**5 + 3*t1**6)
    denom2 = 4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*(k - m)**2)*t1**2 + 8*(3 + 2*(k - m)**2)*t1**3 + (6 + 4*(k - m)**2)*t1**4 + 6*t1**5 + 3*t1**6)
    
    # 分母接近0时返回0，避免数值错误
    if abs(denom1) < 1e-10 or abs(denom2) < 1e-10:
        return 0.0
    
    # 原函数计算逻辑
    part1 = (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + m*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/denom1)
    part2 = ((50 * t1**2) +5* (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + (k - m)*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/denom2)
    return part1 * part2

4. 增加最大细分次数（limit参数） #

通过limit参数提高积分器的最大细分次数，给积分器更多收敛机会：

result, error = integrate.tplquad(
    f, 
    0.1, 2, 
    lambda x: 0.1, lambda x: 2, 
    lambda x,y: 0.1, lambda x,y: 2,
    limit=100  # 把最大细分次数从50提高到100
)

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