Sympy无法分离变量r与θ：表达式化简求助

我有如下符号表达式：

Delta_T*(-S_xx*r**4 - S_xx*r**2*ri**2 - S_xx*r**2*ro**2 + S_xx*ri**2*ro**2 + S_yy*r**4 + S_yy*r**2*ri**2 + S_yy*r**2*ro**2 - S_yy*ri**2*ro**2)*cos(2*theta)/(4*r**2*(ri**2 + ro**2)*log(ro/ri)) - f(r) + g(theta) + (Delta_T*(S_xx + S_yy)*(ri**2 + ro**2)*log(r) - r**2*(Delta_T*S_xx*sin(theta)**2 - Delta_T*S_xx/2 - Delta_T*S_yy*sin(theta)**2 + Delta_T*S_yy/2) - (-2*Delta_T*S_xx*ri**2*ro**2*sin(theta)**2 + Delta_T*S_xx*ri**2*ro**2 + 2*Delta_T*S_yy*ri**2*ro**2*sin(theta)**2 - Delta_T*S_yy*ri**2*ro**2)/(2*r**2))/(2*ri**2*log(ro/ri) + 2*ro**2*log(ro/ri))

该表达式理论上等价于：

g(theta) - f(r) + (2 * Delta_T * sp.log(r) * (S_xx + S_yy) + Delta_T * sp.cos(2 * theta) * (S_yy - S_xx)) / (4 * sp.log(ro / ri))

（注：可能存在符号差异，但不影响，两者差值应为0）

这意味着可以将变量r与θ分离，或者在Sympy中把r与sin(theta)**2、cos(2*theta)分离，但我尝试后没有进展。

尝试的方法及结果 #

我使用collect()方法进行整理，但效果不佳：

执行代码：

expr_collected = sp.collect(V_r_integral - V_theta_integral, [r, sp.cos(theta)**2, sp.sin(theta)**2])

得到结果：

Delta_T*(S_xx*ri**2*ro**2 - S_yy*ri**2*ro**2 + r**4*(-S_xx + S_yy) + r**2*(-S_xx*ri**2 - S_xx*ro**2 + S_yy*ri**2 + S_yy*ro**2))*cos(2*theta)/(4*r**2*(ri**2 + ro**2)*log(ro/ri)) - f(r) + g(theta) + (Delta_T*(S_xx + S_yy)*(ri**2 + ro**2)*log(r) + r**2*(Delta_T*S_xx/2 - Delta_T*S_yy/2 - (Delta_T*S_xx - Delta_T*S_yy)*sin(theta)**2) + (-Delta_T*S_xx*ri**2*ro**2/2 + Delta_T*S_yy*ri**2*ro**2/2 - (-2*Delta_T*S_xx*ri**2*ro**2 + 2*Delta_T*S_yy*ri**2*ro**2)*sin(theta)**2/2)/r**2)/(2*ri**2*log(ro/ri) + 2*ro**2*log(ro/ri))

可以看到仍存在r²与cos(2*theta)相乘的项。

执行代码：

expr_collected = sp.collect(V_r_integral - V_theta_integral, [r**2, sp.cos(theta)**2, sp.sin(theta)**2, sp.cos(2*theta)])

得到类似结果，未实现变量分离。

可复现代码 #

以下是完整的可运行示例代码，最后一行输出的就是需要改写的表达式：

import sympy as sp
r, theta, H0, ri, ro, rho, T_0, Delta_T, S_xx, S_yy = sp.symbols('r theta H0 ri ro rho T_0 Delta_T S_xx S_yy')
J_r = -2 * H0 * (1/r - 1/(ri**2 + ro**2)*(r + ri**2*ro**2/r**3)) * sp.cos(2*theta)
J_theta = -H0 * (1/(ri**2 + ro**2)*(2*r - 2*ri**2*ro**2/r**3)) * sp.sin(2*theta)
J_r = J_r.subs(H0, -(S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri))
J_theta = J_theta.subs(H0, -(S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri))
T = T_0 - Delta_T * sp.log(r/ro) / sp.log(ro/ri)
grad_T_r = sp.diff(T, r)
grad_T_theta = sp.diff(T, theta)
R = sp.Matrix([[sp.cos(theta), sp.sin(theta)],[-sp.sin(theta), sp.cos(theta)]])
S_cartesian = sp.Matrix([[S_xx, 0], [0, S_yy]])
S_polar = R*S_cartesian*R.T
S_rr = S_polar[0,0]
S_rtheta = S_polar[0,1]
S_thetar= S_polar[1,0]
S_thetatheta= S_polar[1,1]
V = sp.Function('V')(r, theta)
grad_V_r = sp.diff(V, r)
grad_V_theta = 1/r * sp.diff(V, theta)
sigma = 1/rho
eq_r = sp.Eq(J_r, -sigma * grad_V_r -sigma * (S_rr * grad_T_r + S_rtheta * grad_T_theta))
eq_theta = sp.Eq(J_theta, -sigma * grad_V_theta - sigma * (S_thetar * grad_T_r + S_thetatheta * grad_T_theta))
sol_grad_V_r = sp.solve(eq_r, grad_V_r)[0]
sol_grad_V_theta = sp.solve(eq_theta, grad_V_theta)[0]
f = sp.Function('f')(r)
g = sp.Function('g')(theta)
V_r_integral = sp.integrate(sol_grad_V_r, r) + g
V_theta_integral = sp.integrate(sol_grad_V_theta * r, theta) + f
print(V_r_integral - V_theta_integral)

恳请帮助实现该表达式的变量分离与化简。

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