<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Expression Manipulation example Ch4
Start with an expression for the volume of a hemi-spherical capped cylinder, and its mass.  Determine length to diameter ratio that gives the minimum mass for a given volume.
Step1:  Do the same with the expression for M
Step2:  Substitute and solve
<ASSISTANT_TASK:>
Python Code:
r, L, D, Vs, Ms = sp.symbols('r L D V M', positive=True)  # NB symbols for V and M here, expressions  below
V = (sp.Rational(4,3)*sp.pi*r**3 + sp.pi*r**2*L).subs({r: D/2})
V
A = (4*sp.pi*r**2 + 2*sp.pi*r*L).subs({r:D/2})
A
rho, p, sigma = sp.symbols('rho p sigma', positive=True)
M = sp.pi*D**2*rho*p/2/sigma*(L+D/2)
M
# Step 1: copy original expresions to check that our manipulations don't change them
Vcopy, Mcopy = V, M
Vcopy.equals(V), Mcopy.equals(M)
# Factor the expression for V
V.factor()
# Divide by D**3 to get D in denominator
V.factor()/D**3
# Expand to get constant term plus one in L/D (quanity of interest)
(V.factor()/D**3).expand()
# collect to get common factor of pi out front (in general use collect_const(expr, c) by pi is special const)
(V.factor()/D**3).expand().collect(sp.pi)
# Multiply by D**3 to get back to original expression
V = D**3*(V.factor()/D**3).expand().collect(sp.pi)
V
Vcopy.equals(V)
# Alternative - replace factor-divide-expand-factor with simplify 
(Vcopy / D**3).simplify().collect(sp.pi)*D**3
Explanation: Expression Manipulation example Ch4
Start with an expression for the volume of a hemi-spherical capped cylinder, and its mass.  Determine length to diameter ratio that gives the minimum mass for a given volume.
End of explanation
M = (M/D**3).expand().collect(sp.pi * p * rho / sigma)*D**3
M
Mcopy.equals(M)
Explanation: Do the same with the expression for M
End of explanation
Veq = sp.Eq(Vs, V)
Veq
t = sp.symbols('t', positive=True)
# Not what we want
Veq.subs({D**3: t})
Veq.xreplace({D**3: t})
# Doesn't work: doesn't keep the L/D ratio we want
sp.solve(Veq.xreplace({D**3: t}), t, dict=True)[0]
# Solution: introduce new variable for L/D
xi = sp.symbols('xi', positive=True)
# Veq.xreplace({D**3: t, L/D: xi}) - doesn't work because of 4 in denominator
# Veq.xreplace({D**3: t, L/D/4: xi/4}) - works because applied to 
Veq.xreplace({D**3: t}).subs({L/D: xi})  # also works
# We can now find the solution for D**3 in terms of V and L/D
D3 = sp.solve(Veq.xreplace({D**3: t, L/D/4: xi/4}), t)[0].subs(xi, L/D)
D3
# And finally get the expression for M in terms of L/D constrained by the volume specification
# (text replaces M for this but I believe that then Mcopy.equals(M) would be False)
Mcons = M.xreplace({D**3: D3})
Mcons
# Minimize by plotting.  First scale by constants and input V
M2 = Mcons/(Vs*p*rho/sigma)
M2
_ = sp.plot(M2.subs(L/D, xi), (xi, 0, 10), xlabel='L/D', ylabel='Dimensionless Mass')
# sp.sympify() aka sp.S() converts non-SymPy to SymPy
sp.S(8).is_real
sp.srepr(D3)
sp.S(sp.srepr(D3))
import inspect
# Message Resolution Order of an object
inspect.getmro(sp.Add)
expr1 = x**2 + x*sp.cos(x) - 1/x
expr1
expr1.atoms()
expr1.free_symbols
expr1.args
expr1.func
m, A, B, C, D = sp.symbols('m A B C D')
t = m/2*(1+1/A)
expr = C**2*t + D**2*t + B
expr
expr.func(*(term/t for term in expr.args))
x, y, z, a, b = sp.symbols('x:z a:b')
expr = x**2 + 2*y*z
expr
integ = sp.Integral(expr, (x, a, b))
integ
integ.bound_symbols
t = sp.symbols('t')
integ.replace(x, t)
sp.srepr(integ)
n = sp.symbols('n', positive=True, integer=True)
sp.Sum(1/2**n, (n, 0, sp.oo)).doit()
n = sp.symbols('n', positive=True, integer=True)
sp.summation(1/2**n, (n,0,sp.oo))
from sympy.stats import Normal
mu0, mu1 = sp.symbols('mu:2', real=True)
mu0, mu1
sigma0, sigma1 = sp.symbols('sigma:2', positive=True)
sigma0, sigma1
N0 = Normal('N_0', mu0, sigma0)
N1 = Normal('N_1', mu1, sigma1)
N0, N1
x = sp.symbols('x', real=True)
p = sp.stats.density(N0)(x) # .subs({mu0:0, sigma0: 1})
q = sp.stats.density(N1)(x)
sp.integrate(p*(sp.log(p) - sp.log(q)), (x, -sp.oo, sp.oo)).simplify()
Explanation: Substitute and solve: ~start by writing volume equation as LHS - RHS~ actually I prefer using sp.Eq
End of explanation
<END_TASK>