Symbolic computation with partial derivatives and a double integral

  • Thread starter Thread starter fluidistic
  • Start date Start date
  • Tags Tags
    Calculus Python
Click For Summary
SUMMARY

This discussion focuses on the computation of a double integral involving a vector field in polar coordinates, specifically ##\vec J = J_r \hat r + J_\theta \hat \theta##, within a quarter annulus defined by theta from 0 to pi/2 and radius from ri to ro. The integral to compute is ##-\int _\Omega T(S_{xx}-S_{yy}) \frac{\partial J_x}{\partial x}d\Omega##, where ##J_x## is derived from ##J_r## and ##J_\theta## using the transformation ##J_x = J_r \cos(\theta) - J_\theta \sin(\theta)##. The discussion also includes a comparison of results obtained from SymPy and Mathematica, highlighting the need for further simplification in the SymPy output.

PREREQUISITES
  • Understanding of vector fields in polar coordinates
  • Familiarity with partial derivatives and the chain rule
  • Experience with symbolic computation using SymPy
  • Knowledge of integral calculus, specifically double integrals
NEXT STEPS
  • Explore advanced techniques in SymPy for simplifying expressions
  • Learn about Mathematica's capabilities for symbolic computation
  • Investigate the implications of real and positive constraints on variables in SymPy
  • Study the application of logarithmic identities in integral simplification
USEFUL FOR

Mathematicians, physicists, and engineers involved in computational modeling, particularly those working with symbolic computation and integral calculus in polar coordinates.

fluidistic
Gold Member
Messages
3,932
Reaction score
283
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
It is apparent from the final expression for Q that SymPy is not eliminating one of \log(r_o/r_i) and \log(r_i/r_o) 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 r_o and r_i 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   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:
  • Informative
Likes   Reactions: fluidistic
pasmith said:
It is apparent from the final expression for Q that SymPy is not eliminating one of \log(r_o/r_i) and \log(r_i/r_o) 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 r_o and r_i 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.
 
  • Like
Likes   Reactions: renormalize
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

  • · Replies 1 ·
Replies
1
Views
2K