Sympy 替代角动量算子

Question

我正在为我的博士学位研究相当复杂和长的方程，我希望我可以使用 sympy 来检查我的方程并更快地获得它们。 这些方程具有球谐函数，我希望可以分解 L^2 运算符以将其替换为 L^2。 但我尝试了很多不同的方法都没有成功，所以我做了几次测试来找出问题出在哪里：

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)
d=sy.Function("d")(zeta)
e=sy.Function("e")(zeta)

#first case without dividing function
test_L=e*a+b*e
Lsquare=a+b

test_L=sy.factor(test_L,Lsquare)
test_L=test_L.replace(Lsquare,Lsqr)
print("test")
print(test_L)
# result : L^2(\theta, \phi)*e(\zeta) => GOOD

#second case with dividing function

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("test")
print(test_L)

# result : (a(\zeta) + b(\zeta)*c(\zeta))*e(\zeta)/c(\zeta) => BAD I would like same form has before 

我想要与第一个测试中相同的形式：L^2*c

我有以下形式的方程的原因：

rac{\left(\sin{\left( heta ight)} rac{\partial^{2}}{\partial heta^{2}} \operatorname{Y^{m}{l}}{\left( heta,\phi 右）} +