Unable to simplify an expression with Sympy

I need to evaluate many integrals, derivatives, etc. symbolically and I have no luck with Sympy. When it is able to perform an integral, the result is gigantic while it should be reducible, but I am totally unable to reach the simplest result. Here's an example:

import sympy as sp

# Define the variables
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')

# Define H0
H0_expr = -(S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri)

# Define the components of the vector field J_r and J_theta
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_x = J_r * sp.cos(theta) - J_theta * sp.sin(theta)
dJ_x_dx = sp.diff(J_x, r) * sp.cos(theta) - 1/r * sp.sin(theta) * sp.diff(J_x, theta) * sp.sin(theta)

T = T_0 + Delta_T * sp.log(r/ro) / sp.log(ro/ri)

# Set up the integral for Q
Q = sp.integrate(r * T * (S_xx - S_yy) * dJ_x_dx, (theta, 0, sp.pi/2), (r, ri, ro))

# Substitute H0 in Q
Q_sub = Q.subs(H0, H0_expr)
Q_sub = sp.simplify(Q_sub)
Q_sub = sp.collect(Q_sub, [ri, ro, sp.log(ro/ri)])
Q_sub = sp.factor(Q_sub)
print(f"This is Q: {Q}")

Running the code yields:

This is Q: -Delta_T*(S_xx - S_yy)**2*(-16*Delta_T*ri**2*log(ri/ro)**2 + 15*pi*Delta_T*ri**2*log(ri/ro)**2 + 80*Delta_T*ri**2*log(ri/ro) + 30*pi*Delta_T*ri**2*log(ri/ro) - 64*Delta_T*ri**2 - 15*pi*Delta_T*ri**2 - 16*Delta_T*ro**2*log(ri/ro)**2 + 15*pi*Delta_T*ro**2*log(ri/ro)**2 + 48*Delta_T*ro**2*log(ri/ro) + 15*pi*Delta_T*ro**2 + 64*Delta_T*ro**2 - 32*T_0*ri**2*log(ri)*log(ro/ri) + 30*pi*T_0*ri**2*log(ri)*log(ro/ri) - 30*pi*T_0*ri**2*log(ro)*log(ro/ri) + 32*T_0*ri**2*log(ro)*log(ro/ri) + 32*T_0*ri**2*log(ro/ri) + 30*pi*T_0*ri**2*log(ro/ri) - 32*T_0*ro**2*log(ri)*log(ro/ri) + 30*pi*T_0*ro**2*log(ri)*log(ro/ri) - 30*pi*T_0*ro**2*log(ro)*log(ro/ri) + 32*T_0*ro**2*log(ro)*log(ro/ri) - 30*pi*T_0*ro**2*log(ro/ri) - 32*T_0*ro**2*log(ro/ri))/(480*rho*(ri**2 + ro**2)*log(ro/ri)**2).

Which is enormous. This should be equal to

$$\frac{\pi\,\Delta T\left(S_{xx}-S_{yy}\right)^2\left[\Delta T\left(r_{o}^{2}-r_{i}^{2}\right)-2\log\left(\frac{r_{o}}{r_{i}}\right)\left(\Delta T\,r_{i}^{2}+T_{0}\left(r_{o}^{2}-r_{i}^{2}\right)\right)\right]}{16\,\rho\left(r_{i}^{2}+r_{o}^{2}\right)\log^{2}\left(\frac{r_{o}}{r_{i}}\right)}$$ 

or

pi*Delta_T*(S_xx - S_yy)**2*(Delta_T*(-ri**2 + ro**2) - 2*(Delta_T*ri**2 + T_0*(-ri**2 + ro**2))*log(ro/ri))/(16*rho*(ri**2 + ro**2)*log(ro/ri)**2)

How can I reach the above "simple" expression using Sympy?

Taking into account the mistake you found into the expression of dJ_x_dx (an extra sin term) and removing the negative sign from H0_expr, I was able to get very very close to the form of your expected solution. Here is how:

import sympy as sp
sp.init_printing()

# Define the variables
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')

# Define H0
H0_expr = (S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri)

# Define the components of the vector field J_r and J_theta
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_x = J_r * sp.cos(theta) - J_theta * sp.sin(theta)
dJ_x_dx = sp.diff(J_x, r) * sp.cos(theta) - 1/r * sp.sin(theta) * sp.diff(J_x, theta)

T = T_0 + Delta_T * sp.log(r/ro) / sp.log(ro/ri)

# Set up the integral for Q
Q = sp.integrate(r * T * (S_xx - S_yy) * dJ_x_dx, (theta, 0, sp.pi/2), (r, ri, ro))

# Substitute H0 in Q
Q_sub = Q.subs(H0, H0_expr)
Q_sub = sp.simplify(Q_sub)
Q_sub = sp.collect(Q_sub, [ri, ro, sp.log(ro/ri)])
Q_sub = sp.factor(Q_sub)
correct_form = sp.pi*Delta_T*(S_xx - S_yy)**2*(Delta_T*(-ri**2 + ro**2) - 2*(Delta_T*ri**2 + T_0*(-ri**2 + ro**2))*sp.log(ro/ri))/(16*rho*(ri**2 + ro**2)*sp.log(ro/ri)**2)
Q_sub

Note that at this point, SymPy is unable to tell if Q_sub is indeed the correct results. The following command should return True, but instead it returns None:

Q_sub.equals(correct_form)

I suspect it has to do with the presence of log(ri/ro) and log(ro/ri). We need to work on the addition at numerator:

num, den = sp.fraction(Q_sub)
a = num.args[-1]
a

After a few manipulations:

b = Delta_T * ro**2 - Delta_T * ri**2
c = a.collect(sp.log(ro / ri)).collect(2*T_0).subs(b, b.simplify()).subs(sp.log(ri/ro), -sp.log(ro/ri)).collect(2 * sp.log(ro/ri))
c

Then:

Q_final = num.subs(a, c) / den
Q_final

You can see it is identical (except for a minus sign that would be non trivial to get it out). You can now confirm the computed result with:

Q_final.equals(correct_form)
# True

NOTE: out of curiosity, I went back and created ro, ri with assumptions: ri, ro = sp.symbols("ri ro", real=True, positive=True). Now, you can directly compare Q_sub to correct_form:

Q_sub.equals(correct_form)
# True

Since you have many integrals to evaluate, I would really consider carefully the application of further manipulation steps: it might get the results more readable, but the risk is to jump into a rabbit hole and waste a lot of time.

---

Original answer: This is not an answer, but another way to verify that your computed expression is different from your theoretical result.

As mentioned above in the comment, you can use spb to visually verify if two expressions evaluates to the same thing. You may want to take the time to insert numbers that represents your physical problem.

Here is a screenshot from Jupyter Notebook.

from spb import *
import ipywidgets
%matplotlib widget

correct = sp.pi*Delta_T*(S_xx - S_yy)**2*(Delta_T*(-ri**2 + ro**2) - 2*(Delta_T*ri**2 + T_0*(-ri**2 + ro**2))*sp.log(ro/ri))/(16*rho*(ri**2 + ro**2)*sp.log(ro/ri)**2)

params = {
    S_xx: (100, 0, 1000),
    S_yy: (200, 0, 1000),
    T_0: (500, 0, 5000),
    rho: (500, 0, 5000),
    (ri, ro): ipywidgets.FloatRangeSlider(value=[20, 30], min=0, max=100),
    # ro: (30, 0, 100),
}
graphics(
    line(Q_sub, (Delta_T, 0, 1000), "Q_sub", params=params), 
    line(correct, (Delta_T, 0, 1000), "computed", params=params),
    ylabel="Q"
)

You can see from the plot that both the sign and magnitude are different.