使用SciPy quad积分时遭遇ZeroDivisionError问题求助

原本用SymPy做积分计算速度太慢，切换到SciPy的quad后性能提升，但运行代码时触发ZeroDivisionError。

原代码如下：

from scipy.integrate import quad
import sympy as S

x = S.Symbol('x')
y = S.Symbol('y')
xi = 0.75
a = 10; b = 10; f1 = 0.5; f2 = 0.5; f0 = f1+f2; al = -f1/f0; be = -f2/f0 
F0 = f0*(al*(x**2/a**2)*xi+be*(y*2/b**2)**xi+1)
j2 = F0.diff(y,2)
jj2 = S.lambdify([y],j2,'scipy')
J2_ = quad(jj2,-a,a)
J2 = (J2_[0]*a**2)/f0

错误信息：

File "C:\Users\Mikhail\Desktop\robpy\cyc.py", line 60, in raschet
J2_ = quad(jj2,-a,a)
File "C:\Users\Mikhail\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\integrate\_quadpack_py.py", line 463, in quad
retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
File "C:\Users\Mikhail\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\integrate\_quadpack_py.py", line 575, in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
File "<lambdifygenerated-2>", line 2, in _lambdifygenerated
ZeroDivisionError: float division by zero

F0对y的二阶导数表达式为：
-0.0118585412256314*(y**2)**0.75/y**2

推测错误源于分母y²，但用SymPy直接积分S.integrate(j2,(y,-a,a))能正常得到结果-0.15。需求是支持用户输入不同初始数据，无法直接用浮点数重构F0，需要解决除零问题。

方法1：用SymPy化简表达式后再lambdify #

原二阶导数表达式可以通过SymPy化简消除分母的y²，因为(y²)^0.75 = |y|^1.5，化简后表达式变为-0.011858... / |y|^0.5，避免除零。

修改代码如下：

from scipy.integrate import quad
import sympy as S

x = S.Symbol('x')
y = S.Symbol('y')
xi = 0.75
a = 10; b = 10; f1 = 0.5; f2 = 0.5; f0 = f1+f2; al = -f1/f0; be = -f2/f0 
F0 = f0*(al*(x**2/a**2)*xi+be*(y*2/b**2)**xi+1)
j2 = F0.diff(y,2)
# 化简表达式
j2_simplified = S.simplify(j2)
jj2 = S.lambdify([y], j2_simplified, 'scipy')
J2_ = quad(jj2, -a, a)
J2 = (J2_[0]*a**2)/f0

方法2：自定义包装函数处理奇点 #

由于原函数在y=0处是可积奇点（SymPy能正常积分），可以写一个包装函数，当y接近0时用极限值代替，避免除零。

修改代码如下：

from scipy.integrate import quad
import sympy as S
import numpy as np

x = S.Symbol('x')
y = S.Symbol('y')
xi = 0.75
a = 10; b = 10; f1 = 0.5; f2 = 0.5; f0 = f1+f2; al = -f1/f0; be = -f2/f0 
F0 = f0*(al*(x**2/a**2)*xi+be*(y*2/b**2)**xi+1)
j2 = F0.diff(y,2)
# 计算y→0时的极限值
limit_val = float(S.limit(j2, y, 0))
# 生成原函数的lambda
jj2_raw = S.lambdify([y], j2, 'scipy')
# 包装函数处理奇点
def jj2(y_val):
    if np.abs(y_val) < 1e-10:
        return limit_val
    return jj2_raw(y_val)

J2_ = quad(jj2, -a, a)
J2 = (J2_[0]*a**2)/f0

方法3：利用SciPy quad的奇点处理参数 #

quad支持通过points参数指定积分区间内的奇点，让积分器专门处理该点，避免除零错误。

修改代码如下：

from scipy.integrate import quad
import sympy as S

x = S.Symbol('x')
y = S.Symbol('y')
xi = 0.75
a = 10; b = 10; f1 = 0.5; f2 = 0.5; f0 = f1+f2; al = -f1/f0; be = -f2/f0 
F0 = f0*(al*(x**2/a**2)*xi+be*(y*2/b**2)**xi+1)
j2 = F0.diff(y,2)
jj2 = S.lambdify([y], j2, 'scipy')
# 指定奇点位置为0
J2_ = quad(jj2, -a, a, points=[0])
J2 = (J2_[0]*a**2)/f0

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