Symbolic computation with partial derivatives and a double integral

  • Thread starter Thread starter fluidistic
  • Start date Start date
  • Tags Tags
    Calculus Python
Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
5 replies · 3K views
fluidistic
Gold Member
Messages
3,934
Reaction score
286
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?
 
on Phys.org
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:
Reply
  • Informative
  • Like
Likes   Reactions: fluidistic and renormalize
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:
Reply
  • Informative
Likes   Reactions: fluidistic
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.
 
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.
 
Reply
  • Like
Likes   Reactions: renormalize