Evaluating integrals, derivatives, etc. with Simpy...

[fluidistic]

I would like to evaluate expressions with Simpy, but unfortunately I am unable to get a simple answer, the one I would get by hand if I had the time to perform all the computations. As far as I understand, Mathematica does it and yields 4 times the Simpy result, which is a big worry since I wish to trust the result I get.

I have the expression of a vector field $\vec J$ in polar coordinates, in a region of a quarter of an annulus (2D), theta goes from 0 to pi/2 and r goes from ri to ro. I need to compute $\int _\Omega \rho |\vec J|^2d\Omega =\int_0^{\pi/2} \int_{ri}^{ro}\vec J \cdot \vec J rdrd\theta$.

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

# Calculate the magnitude squared of J (J^2)
J_squared = J_r**2 + J_theta**2

# Set up the integral for Q_total
Q_total = sp.integrate(r * rho * J_squared, (theta, 0, sp.pi/2), (r, ri, ro))

# Substitute H0 in Q_total
Q_total_sub = Q_total.subs(H0, H0_expr)
Q_total_sub = sp.simplify(Q_total_sub)
print(f"This is Q_total, the total Q: {Q_total_sub}.")

This yields:

This is Q_total, the total Q: pi*Delta_T**2*(S_xx - S_yy)**2*(-ri**4 + 2*ri**2*(ri**2 + ro**2) + ri**2*(ri**2 + 2*ro**2) + ro**4 - 2*ro**2*(ri**2 + ro**2) - ro**2*(2*ri**2 + ro**2) - 2*(ri**2 + ro**2)**2*log(ri) + 2*(ri**2 + ro**2)**2*log(ro))/(32*rho*(ri**4 + 2*ri**2*ro**2 + ro**4)*log(ro/ri)**2)

If I take it from here, by hand, I reach $Q=\pi (\Delta T)^2(S_\text{xx}-S_\text{yy})^2\frac{\left[ri^2-ro^2+(ri^2+ro^2)\ln\left( \frac{ro}{ri} \right)\right]}{16\rho(ri^2+ro^2)\ln^2(\frac{ro}{ri})}$. How can I reach this result using Simpy? Can someone else confirm that Mathematica or other software yields 4 times this result?

 

[renormalize]

I reach $Q=\pi (\Delta T)^2(S_\text{xx}-S_\text{yy})^2\frac{\left[ri^2-ro^2+(ri^2+ro^2)\ln\left( \frac{ro}{ri} \right)\right]}{16\rho(ri^2+ro^2)\ln^2(\frac{ro}{ri})}$. How can I reach this result using Simpy? Can someone else confirm that Mathematica or other software yields 4 times this result?

Actually, Mathematica gives precisely the same answer (no factor of 4 difference!):