使用SymPy计算拉普拉斯算子时出现除零错误的问题

问题背景 #

计算函数 $f(x,y)=r^{2/3}(1-x2)(1-y^2)\sin(\frac{2}{3}\theta)$ 的拉普拉斯算子（$r,\theta$ 为极坐标，$0<\theta<2\pi$），函数在 $(-1,0)$ 附近光滑。使用SymPy实现的代码如下：

import sympy
x,y = sympy.symbols('x y')
r = sympy.sqrt(x**2+y**2)
A1 = sympy.acos(x/r)
A2 = 2*sympy.pi-sympy.acos(x/r)  
theta = sympy.Piecewise((A1,y>=0),(A2,y<0)) #compute polar coordinate with range [0,2pi)

expr = r**(2/3)*(1-x**2)*(1-y**2)*sympy.sin((2/3)*theta)
uxx = sympy.diff(expr,x,x)
uyy = sympy.diff(expr,y,y)
lapl=-(uxx+uyy)
fxy = sympy.lambdify([x,y],lapl,'math')
u_ = sympy.lambdify([x,y],expr,'math')

遇到的问题：

• 直接计算 fxy(-1,0) 触发除零错误；
• y取-0.01、-0.001等趋近于0的负数时结果收敛，但y取-1e-8时再次触发除零错误。

错误原因 #

1. 自定义theta的分段导数问题：手动用Piecewise定义theta时，对其求导后会引入含y的分母项，代入y=0时直接触发除零；当y取极小值时，浮点数精度误差会导致部分项的计算出现0/0或分母溢出的情况。
2. acos(x/r)的数值稳定性问题：当y趋近于0且x=-1时，x/r = -1/sqrt(1+y²)，浮点数计算中可能出现参数触及acos定义域边缘的情况，间接引发除零或非法操作。

解决方法 #

1. 改用内置atan2定义极角theta #

SymPy的atan2(y, x)直接返回[0, 2π)范围的极角，符号求导和数值计算的稳定性远优于手动分段定义，避免导数中的除零项。修改后的代码：

import sympy
x,y = sympy.symbols('x y')
r = sympy.sqrt(x**2+y**2)
theta = sympy.atan2(y, x)  # 替换自定义的Piecewise theta

expr = r**(2/3)*(1-x**2)*(1-y**2)*sympy.sin((2/3)*theta)
uxx = sympy.diff(expr,x,x)
uyy = sympy.diff(expr,y,y)
lapl=-(uxx+uyy)
# 先化简表达式，提升数值稳定性
lapl_simplified = sympy.simplify(lapl)
fxy = sympy.lambdify([x,y], lapl_simplified,'math')
u_ = sympy.lambdify([x,y], expr,'math')

2. 符号计算后先化简再数值化 #

对拉普拉斯算子的符号表达式调用sympy.simplify()，可以消除冗余的分母项，减少数值计算时的精度问题。

3. 直接计算符号极限获取(-1,0)处的值 #

对于(-1,0)点，可以通过符号极限计算准确值，避免数值代入的问题：

# 计算y→0⁻时的极限
limit_val = sympy.limit(sympy.limit(lapl_simplified, y, 0, dir='-'), x, -1)
print(limit_val.evalf())  # 输出数值结果

验证 #

修改后的代码可以直接计算fxy(-1,0)，不会触发除零错误；y取极小负数时的数值计算也会保持稳定，结果收敛于符号极限值。

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