Sympy替换角动量算子L²遇分式场景失效，求解决方案

问题场景 #

在博士研究的方程推导中，需用Sympy验证含球面谐波的复杂方程，核心需求是将展开的L²算子形式替换为L²(θ,φ)符号，但含分式的表达式无法按预期完成替换。

测试代码与问题复现 #

以下测试代码展示了两种场景的结果差异：

import sympy as sy

zeta,theta,phi=sy.symbols(r"\zeta \theta \phi")
Lsqr=sy.Function(r"L^2")(theta,phi)

a=sy.Function("a")(zeta)
b=sy.Function("b")(zeta)
c=sy.Function("c")(zeta)
e=sy.Function("e")(zeta)

# 无分式情况：替换符合预期
test_L=e*a+b*e
Lsquare=a+b
test_L=sy.factor(test_L,Lsquare)
test_L=test_L.replace(Lsquare,Lsqr)
print("无分式测试结果：")
print(test_L)
# 输出: L^2(\theta, \phi)*e(\zeta)

# 含分式情况：替换不符合预期
test_L=e*a/c+b*e
Lsquare=a/c+b
test_L=sy.factor(test_L,Lsquare)
test_L=test_L.replace(Lsquare,Lsqr)
print("\n含分式测试结果：")
print(test_L)
# 输出: (a(\zeta) + b(\zeta)*c(\zeta))*e(\zeta)/c(\zeta)
# 期望输出: L^2(\theta, \phi)*e(\zeta)

实际待处理方程 #

原始方程：
$$\frac{\left(\sin{\left(\theta \right)} \frac{\partial^{2}}{\partial \theta^{2}} \operatorname{Y^{m}{l}}{\left(\theta,\phi \right)} + \cos{\left(\theta \right)} \frac{\partial}{\partial \theta} \operatorname{Y^{m}{l}}{\left(\theta,\phi \right)}\right) \rho_{0}{\left(\zeta,\theta \right)}}{\sin{\left(\theta \right)}} + \frac{\rho_{0}{\left(\zeta,\theta \right)} \frac{\partial^{2}}{\partial \phi^{2}} \operatorname{Y^{m}_{l}}{\left(\theta,\phi \right)}}{\sin^{2}{\left(\theta \right)}}$$

期望转换为：
$$\rho_0(\zeta,\theta) L^2 Y_l^m(\theta,\phi)$$

解决方案 #

问题根源是factor()会合并分式改变表达式结构，导致无法匹配目标替换式。改用collect()提取公共因子后直接替换即可：

修改后的测试代码 #

import sympy as sy

zeta,theta,phi=sy.symbols(r"\zeta \theta \phi")
Lsqr=sy.Function(r"L^2")(theta,phi)

a=sy.Function("a")(zeta)
b=sy.Function("b")(zeta)
c=sy.Function("c")(zeta)
e=sy.Function("e")(zeta)

# 含分式情况：使用collect提取公共因子后替换
test_L=e*a/c+b*e
Lsquare=a/c+b
# 提取公共因子e
test_L_collected = sy.collect(test_L, e)
# 替换目标表达式
test_L_replaced = test_L_collected.replace(Lsquare, Lsqr)
print("修改后的含分式测试结果：")
print(test_L_replaced)
# 输出: L^2(\theta, \phi)*e(\zeta)

实际方程的处理代码 #

import sympy as sy

zeta,theta,phi=sy.symbols(r"\zeta \theta \phi")
Y = sy.Function(r"Y_l^m")(theta, phi)
rho0 = sy.Function(r"\rho_0")(zeta, theta)
Lsqr = sy.Function(r"L^2")(theta, phi)

# 定义原始方程
term1 = (sy.sin(theta)*sy.diff(Y, theta, 2) + sy.cos(theta)*sy.diff(Y, theta)) * rho0 / sy.sin(theta)
term2 = rho0 * sy.diff(Y, phi, 2) / (sy.sin(theta)**2)
original_eq = term1 + term2

# 提取公共因子rho0
collected_eq = sy.collect(original_eq, rho0)
# 定义L²算子的展开形式
Lsquare_expanded = (sy.sin(theta)*sy.diff(Y, theta, 2) + sy.cos(theta)*sy.diff(Y, theta))/sy.sin(theta) + sy.diff(Y, phi, 2)/(sy.sin(theta)**2)
# 替换得到目标形式
target_eq = collected_eq.replace(Lsquare_expanded, Lsqr)
print("转换后的实际方程：")
sy.pprint(target_eq)
# 输出: \rho_0(\zeta, \theta) L^2(\theta, \phi)

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