如何配置Python SciPy的nquad选项，在奇点附近积分以避免除零错误

你遇到的积分无法收敛、触发除零的问题，和limit、epsrel的取值大小没有直接关系，本质是SciPy nquad 默认调用的全局自适应积分逻辑，对紧邻截断边界的窄幅尖峰识别能力不足：

• 哪怕你把c的积分上限截断到0.99，被积函数在c→0.99的区间仍会因为距离c=1的奇点极近，出现宽度极窄、幅值极高的尖峰
• 积分器默认按等长逻辑初始划分积分区间，采样点根本不会落在这个窄尖峰范围内，后续自适应加密也因为初始采样点离尖峰太远，无法识别到尖峰的积分贡献，最终触发除零错误或者返回完全失真的结果
• Mathematica的NIntegrate内置全局奇异性探测逻辑，会自动对尖峰区域加密网格，因此不需要手动干预就能完成计算，但SciPy nquad没有内置的全局奇异性追踪能力，必须手动调整积分结构才能正常运行

不需要把参数调到极端值，核心是手动拆分靠近奇点的积分区间，给尖峰区域单独分配积分精度：

• 对无奇异的a、b维度（积分范围0~30），使用默认精度即可，不需要拉高limit阈值
• 对c维度的积分，手动拆分为两个子区间：远离奇点的[-1, 0.8]、靠近奇点的[0.8, 0.99]，后者单独拉高积分划分上限
• 替换nquad内置的c维度积分逻辑，用手动拆分的quad计算每层c的积分值

可直接运行的修正代码 #

from scipy.integrate import quad

results = []

for d in range(-4, 5):
    f = lambda a, b, c: a**2 * b**2 / (a**2 - 2*a*b*c + b**2 + d)

    # 单独计算单组(a,b)对应的c维度积分，拆分尖峰区间
    def calc_c_integral(a, b):
        # 远离奇点的常规区间
        res_far, _ = quad(lambda c: f(a,b,c), -1, 0.8, limit=100)
        # 靠近c=1奇点的尖峰区间，提高划分精度
        res_peak, _ = quad(lambda c: f(a,b,c), 0.8, 0.99, limit=10000, epsrel=1e-4)
        return res_far + res_peak

    # 嵌套计算b、a维度的外层积分
    def calc_b_integral(a):
        res, _ = quad(lambda b: calc_c_integral(a,b), 0, 30, limit=100)
        return res

    total, _ = quad(calc_b_integral, 0, 30, limit=100)
    results.append([d, 4 * 10**(-8) * total])

print(results)

• 单独给整个c区间[-1, 0.99]设置极大的limit参数没有效果：因为积分器的初始采样点是均匀分布的，窄尖峰可能仅占整个区间长度的1%不到，初始采样点根本不会触碰到尖峰位置，哪怕limit设到100000也无法识别尖峰贡献
• 无限制放宽epsrel阈值只会让积分器更快放弃精度要求，不会解决漏扫尖峰的根本问题

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