Python积分代码执行缓慢且报IntegrationWarning问题求助

问题背景 #

使用scipy.nquad计算多重积分时，出现耗时极久、触发收敛警告后无输出的问题，警告信息如下：

/usr/local/lib/python3.11/dist-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The integral is probably divergent, or slowly convergent.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,
/usr/local/lib/python3.11/dist-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,

优化措施 #

1. 合并冗余函数：将结构完全一致的psi1/psi2、phi1/phi2等函数合并为通用函数，减少代码冗余与重复计算
2. 调整积分区间：避开t1、t2可能导致分母为0的区间，拆分慢收敛区间，降低积分难度
3. 配置积分参数：给nquad传入opts参数，提高细分次数、调整绝对/相对误差阈值，适配慢收敛场景
4. 简化嵌套层级：将F1、F2中的嵌套积分替换为单次一维积分，减少计算开销
5. 增强数值稳定性：对分母项添加极小值避免除零，修正原代码中10*9的笔误为90，防止计算错误

修改后的代码 #

from functools import partial
import numpy as np
import scipy.special as special
from scipy import integrate
from scipy.special import kn, hyp1f1
import matplotlib.pyplot as plt
import math
import time, sys
from scipy.integrate import nquad

# 全局参数
a = 0.01
nab = 10**-9
kz = 0.002
H = 10**-5
start = 100 * H
end = 10**9 * H
step = 10**8 * H
ti = 10 * H
tf = 100 * H
epsilon = 1e-10  # 全局防除零极小值

# 合并通用函数，替代原psi1/psi2、psik1/psik2等重复定义
def psi(m, t):
    return 0.1 + t * (0.1 + t * (-6.4576 + m**2 * (-0.625 - 21.66203 * t) - 100.282 * t - 1.25 * m**4 * t))

def phi(m, t, k):
    return 0.1 + t * (10.30555 + 6.1 * (k - m)**6 * t**2 + 
                      (k - m)**4 * t * (3. + 148.022222 * t) + 
                      t * (224.953549 + 4034.491 * t) + 
                      (k - m)**2 * (1. + t * (50.45555 + 1284.838 * t)))

def psik(m, t, k):
    return 0.1 + t * (0.1 + t * (-6.4576 + (k - m)**2 * (-0.625 - 21.66203 * t) - 
                      100.282 * t - 1.25 * (k - m)**4 * t))

def phim(m, t):
    return 0.1 + t * (10.30555 + 6.1 * m**6 * t**2 + 
                      m**4 * t * (3. + 148.022222 * t) + 
                      t * (224.953549 + 4034.491 * t) + 
                      m**2 * (1. + t * (50.45555 + 1284.838 * t)))

# 重构核心计算函数
def f(m, t, k):
    return -psi(m, t) * phi(m, t, k) - 5.0 * psi(m, t) * psik(m, t, k)

def fm(m, t, k):
    return -psik(m, t, k) * phim(m, t) - 5.0 * psik(m, t, k) * psi(m, t)

def ffmk(t1, t2, m, k):
    f_val = f(m, t1, k)
    return f_val * (fm(m, t2, k) + f_val)

def parf(t1, t2, m, mu, a, k):
    # 确保t1、t2不为零，避免除零错误
    t1_safe = max(abs(t1), epsilon)
    t2_safe = max(abs(t2), epsilon)
    denom = (k**2 - 2 * k * mu * m + m**2)**(3/2)
    denom_safe = denom if denom > epsilon else epsilon
    return (H * a)**2 * k**3 * m**3 * (1/(t1_safe*H)) * (1/(t2_safe*H)) * \
           (4/9)**2 * nab**2 * (1/kz)*0.2 * (1 - 2 * mu**2 + mu**4) * (1/denom_safe)

def ptot(t1, t2, m, mu, a, k):
    return parf(t1, t2, m, mu, a, k) * ffmk(t1, t2, m, k)

def rho(t1, t2, m, mu, a, k):
    # 修正原代码中10*9的笔误为90
    return k**2 * (H * a * 90)**2 * ptot(t1, t2, m, mu, a, k)

# 积分配置：提高细分次数，调整误差阈值适配慢收敛场景
quad_opts = {'limit': 200, 'epsabs': 1e-8, 'epsrel': 1e-8}
nquad_opts = [quad_opts]*4  # 对应4个积分变量的配置

def p21(a, k):
    # 调整t1、t2区间，避开可能的奇点
    ranges = [[-H + epsilon, a], [-H + epsilon, a], [0.1, 1], [-1, 1]]
    integral_result, error = nquad(rho, ranges, args=(a, k), opts=nquad_opts)
    return integral_result

def p22(a, k):
    ranges = [[-ti + epsilon, a], [-ti + epsilon, a], [0.1, 1], [-1, 1]]
    integral_result1, error1 = nquad(rho, ranges, args=(a, k), opts=nquad_opts)
    return integral_result1

delta_t1 = -H + a
delta_t2 = -ti + a

def F1(x):
    res = np.zeros_like(x, dtype=float)
    for i, val in enumerate(x):
        k_val = val/(2*np.pi)
        # 用quad替代嵌套nquad，减少计算层级
        y, err = integrate.quad(lambda a_val: p21(a_val, k_val), H, a)
        res[i] = y / delta_t1
    return res

def F2(x):
    res1 = np.zeros_like(x, dtype=float)
    for i, val in enumerate(x):
        k_val = val/(2*np.pi)
        y, err = integrate.quad(lambda a_val: p22(a_val, k_val), ti + epsilon, a)
        res1[i] = y / delta_t2
    return res

# 生成数据并绘图
X = np.arange(0.1*H, end, step)
F1X = np.nan_to_num(F1(X))
F2X = np.nan_to_num(F2(X))

plt.plot(X, F1X, color='b', label='$\eta_a=-H$')
plt.plot(X, F2X, color='red', label='$\eta_a=-tf$')
plt.xlabel(" $f~ Hz$ ")
plt.ylabel(r"$\rho_h(k,\tau \rightarrow 0)$")
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.savefig('Pi-k.png')

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