Symbolic computation with partial derivatives and a double integral

  • #1
fluidistic
Gold Member
3,955
266
I have an expression of a vector field ##\vec J = J_r \hat r + J_\theta \hat \theta## in polar coordinates, in a region that's a quarter of an annulus, where theta ranges between 0 and pi/2 and r goes from ri to ro.
I need to compute ##-\int _\Omega T(S_{xx}-S_{yy}) \frac{\partial J_x}{\partial x}d\Omega =-(S_{xx}-S_{yy})\int_0^{\pi/2} \int_{ri}^{ro} [T_0 + \Delta T \ln(r/ro) / \ln(ro/ri)] \frac{\partial J_x}{\partial x} rdrd\theta## where all the parameters in the last expressions are known constants.
The thing is that I do not have the ##\hat x## component of ##\vec J##, I have the ##\hat r## and ##\hat \theta## components only.
Therefore I need to obtain ##J_x## first. I get that ##J_x = J_r \cos(\theta)-J_\theta \sin(\theta)##.
Using the chain rule I get that ##\frac{\partial }{\partial x}=\cos(\theta)\frac{\partial}{\partial r}-\frac{\sin(\theta)}{r}\frac{\partial}{\partial \theta}##.
Having derived this, I can now write the Simpy code that should give me the sought integral.

Python:
import simpy 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}")

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).

I wonder what Mathematica (or a better Simpy code) can yield, I strongly suspect the expression can be further simplified. Can someone confirm my result or show me a shoter one, which I could then compare?
 
Physics news on Phys.org
  • #2
It is apparent from the final expression for [itex]Q[/itex] that SymPy is not eliminating one of [itex]\log(r_o/r_i)[/itex] and [itex]\log(r_i/r_o)[/itex] in favour of the other, so you could try explicitly telling it to do that by making the substitution
Python:
Q_sub = Q_sub.subs(sp.log(ri/ro), -sp.log(ro/ri))
before simplifying and collecting terms.

If you explicitly tell SymPy that [itex]r_o[/itex] and [itex]r_i[/itex] are real and positive by declaring them as
Python:
ro = sp.Dummy('ro', real=True, positive=True)
etc. rather than just using sp.symbols, then SymPy might be able to apply further simplifications which don't hold for negative or complex arguments.
 
Last edited:
  • Informative
  • Like
Likes fluidistic and renormalize
  • #3
fluidistic said:
I wonder what Mathematica (or a better Simpy code) can yield, I strongly suspect the expression can be further simplified. Can someone confirm my result or show me a shoter one, which I could then compare?
Here is what Mathematica gives for your heat expression:
1733175460483.png

Cleaning up this result a bit, I repeat it below in LaTeX, along with the Joule heat from your previous thread https://www.physicsforums.com/threads/evaluating-integrals-derivatives-etc-with-simpy.1067101/ :$$Q_{\text{Joule}}=\frac{\pi\,\left(\Delta T\right)^{2}\left(S_{xx}-S_{yy}\right)^2\left[\left(r_{i}^{2}+r_{o}^{2}\right)\log\left(\frac{r_{o}}{r_{i}}\right)+r_{i}^{2}-r_{o}^{2}\right]}{16\,\rho\left(r_{i}^{2}+r_{o}^{2}\right)\log^{2}\left(\frac{r_{o}}{r_{i}}\right)}$$$$Q_{\text{Bridgman}}=\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)}$$Perhaps you can use this to compare to your SymPy expression after applying the suggestions of @pasmith in post #2?
 
Last edited:
  • Informative
Likes fluidistic
  • #4
pasmith said:
It is apparent from the final expression for [itex]Q[/itex] that SymPy is not eliminating one of [itex]\log(r_o/r_i)[/itex] and [itex]\log(r_i/r_o)[/itex] in favour of the other, so you could try explicitly telling it to do that by making the substitution
Python:
Q_sub = Q_sub.subs(sp.log(ri/ro), -sp.log(ro/ri))
before simplifying and collecting terms.

If you explicitly tell SymPy that [itex]r_o[/itex] and [itex]r_i[/itex] are real and positive by declaring them as
Python:
ro = sp.Dummy('ro', real=True, positive=True)
etc. rather than just using sp.symbols, then SymPy might be able to apply further simplifications which don't hold for negative or complex arguments.
Thanks for the suggestion. The expression is still very large after the substitution of the logs.
And the computations is taking forever (I will have to abort) when I specify both ri and ro to be real and positive.
 
  • #5
renormalize said:
Here is what Mathematica gives for your heat expression:
View attachment 354047
Cleaning up this result a bit, I repeat it below in LaTeX, along with the Joule heat from your previous thread https://www.physicsforums.com/threads/evaluating-integrals-derivatives-etc-with-simpy.1067101/ :$$Q_{\text{Joule}}=\frac{\pi\,\left(\Delta T\right)^{2}\left(S_{xx}-S_{yy}\right)\left[\left(r_{i}^{2}+r_{o}^{2}\right)\log\left(\frac{r_{o}}{r_{i}}\right)+r_{i}^{2}-r_{o}^{2}\right]}{16\,\rho\left(r_{i}^{2}+r_{o}^{2}\right)\log^{2}\left(\frac{r_{o}}{r_{i}}\right)}$$$$Q_{\text{Bridgman}}=\frac{\pi\,\Delta T\left(S_{xx}-S_{yy}\right)\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)}$$Perhaps you can use this to compare to your SymPy expression after applying the suggestions of @pasmith in post #2?
That's amazing, thanks a lot! Woah... There's a little typo in your latex (missing square for the S_xx and S_yy term).
This helps a lot.
 
  • Like
Likes renormalize
  • #6
fluidistic said:
There's a little typo in your latex (missing square for the S_xx and S_yy term).
Thanks, fixed it in post #3.
 

Similar threads

Back
Top