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?

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?

• python
• sympy
• symbolic-math

• I tried running Q.simplify().combsimp() which reduces it from 450 ops down to 190. This seems like a lot to ask of sympy, I'd suggest trying in Maple if you have access, it might give better results. – Z4-tier Commented yesterday
• Could you use the SymPy form of the latex simplified version in the OP, please. – smichr Commented yesterday
• You print Q .Why not Q_sub? Did you check all of the intermediate sub steps to see if any work? – hpaulj Commented 23 hours ago
• Not sure what you mean @smichr, this seems to be my goal that I am unable to reach, no? I got that simple latex expression from a friend who has Mathematica, but I am unable to reach it with Sympy. – untreated_paramediensis_karnik Commented 22 hours ago
• Yes, sorry @hpaulj , none of the intermediary results look close enough. It's as if Sympy is unable to collect log ro and log ri terms as it should. – untreated_paramediensis_karnik Commented 22 hours ago

2 Answers 2

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.

I inserted the SymPy version of the expected simplified form. If that is correct then your Q_sub is definitely)?) not the same as the simplified form q.

>>> q=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)
>>> n,d = Q_sub.as_numer_denom()
>>> (n*denom(q)-d*numer(q)).expand().subs(log(ri/ro),-log(ro/ri)).expand()==0
False

The log replacement is there to make the log terms appear in the same way.

You can also do a quick numerical check of the difference. If you get a non-zero (and not None) result it can indicate that the expression is not zero, e.g.

>>> dif = Q_sub - q
>>> dif._random()  # supplies random values for variables
-0.0065733029461819 - 0.00737479065426449*I