scipy.nquad与Mathematica数值积分结果不一致问题排查

问题背景 #

我正在用Python的scipy.integrate.nquad实现一个复杂数值积分函数，但计算结果和Mathematica的NIntegrate结果差异明显，希望定位问题并修正Python代码，保留其计算速度快的优势。

Python实现代码 #

import numpy as np
import scipy as sp
import sympy as smp
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.integrate import nquad
import cmath
from scipy.special import gamma # 需要scipy的gamma函数
# 加载scipy.special中的Gegenbauer函数
from scipy.special import eval_gegenbauer as G # gegenbauer

options={'limit':1000, 'epsabs':1e-12, 'epsrel':1e-12}

def IntZ(z,sp,s,l,lp):
    return (1-z**2)**((d-4)/2)/2*G(l,(d-3)/2,z)*G(lp,(d-3)/2,cmath.sqrt(-(3*s*(-1+z**2)+sp*(3+z**2))/((s*(-1+z**2)-sp*(3+z**2)))))
def ker(z,sp,s):
#     return -2+s/(sp-s)+(4*sp*(s+2*sp))/((s+2*sp)**2-s**2*z**2)
    return s*(1/(complex(0,-10**(-20))-s+sp)+(-1+z)/(s+2*sp-s*z)-(1+z)/(s+2*sp+s*z))

# 积分项的所有前置因子
def prefactor(sp,s,sigma,l,lp):
    return ((gamma(l+1)*gamma((d-3)/2)*2**(2*d-3)*(np.pi)**((d-3)/2)*gamma((d-3)/2)*(2*lp+d-3))/(2**(d+1)*(np.pi)**((d-1)/2)*gamma(l+d-3)))*(((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2)))

# 整合所有因子后的最终积分项（sp和z为积分变量）
def integrand(z,sp,s,sigma,l,lp):
    return IntZ(z,sp,s,l,lp)*ker(z,sp,s)*prefactor(sp,s,sigma,l,lp)

def prefactor2(sp,s,sigma,lp):
    return ((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2))

def fell(s,sigma,l,lp):
    return nquad(integrand,[[-1,1],[1,np.inf]],args=(s,sigma,lp,l),opts=[options,options])

Mathematica实现代码 #

d = 5; 
\[Mu] = 0; 
(* z积分采用数值计算 *)
Integral2[s_, \[Sigma]_, l_, lp_, MaxRes_, WorkPres_] := 
 NIntegrate[(l! Gamma[(d - 3)/2])/(
   2^(d + 1) \[Pi]^((d - 1)/2) Gamma[l + d - 3]) (s - \[Mu])^((d - 3)/
    2)/s^(1/2) 2^(2 d - 3) \[Pi]^((d - 3)/2) Gamma[(d - 3)/2] sp^(
   1/2)/(sp - \[Mu])^((d - 3)/2) (2 lp + d - 3) ((sp - 1)/sp)^(
   2 lp + 1/2) (2 Sqrt[(sp - 1) \[Sigma]])/(
   1 - sp - \[Sigma]) (s/(sp - (s + I 10^-20)) + t/(sp - t) + u/(
       sp - u) /. {u -> -s - t} /. {t -> (s (z - 1))/2}) (1 - z^2)^((
   d - 4)/2)/
   2 GegenbauerC[l, (d - 3)/2, 
    z] (GegenbauerC[lp, (d - 3)/2, Sqrt[
        1 + (4 a)/(-a + sp)]] /. {a -> (s t u)/(
         s t + t u + u s)} /. {u -> -s - t} /. {t -> (s (z - 1))/
       2}), {z, -1, 1}, {sp, 1, Infinity}, MaxRecursion -> MaxRes, 
  WorkingPrecision -> WorkPres, Method -> "GlobalAdaptive"]

测试结果对比 #

测试用例1：fell(200,45,4,4) #

• Python的nquad返回（仅实部）：(915.9500431066836, 0.0008916854858398438)
• Mathematica（MaxRecursion=100，WorkingPrecision=100）返回：

  946.243126048692722003932251899650713550717875143021622444266774813501
  3068351836006101088945258421851 - 
    288.75854540972949504131273023372598463063894553182017391668357286673
  70457344996286394039275246094406 I
  

测试用例2：fell(1800,1800,4,4) #

• Python的nquad返回（仅实部）：(5582.867565858222, 90.11617121750788)
• Mathematica（MaxRecursion=100，WorkingPrecision=100）返回：

  6600.16959988860178799109412672538944693500880870136958557262039403043
  5553207340978200772084077156640 - 
    5628.2185973484253452277428879848996442046986842344701373348873276687
  65333223655602231358008282651070 I
  

已尝试调试步骤 #

调整nquad的limit、epsabs、epsrel参数，但结果变化不大，希望找到导致结果差异的根本原因。

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