Modeling and Simulation in Python

SymPy code for Chapter 16

We can figure out the final temperature of a mixture by setting the total heat flow to zero and then solving for $T$.



In [1]:

    
from sympy import *

init_printing()



In [2]:

    
C1, C2, T1, T2, T = symbols('C1 C2 T1 T2 T')

eq = Eq(C1 * (T - T1) + C2 * (T - T2), 0)
eq









    Out[2]:





$\displaystyle C_{1} \left(T - T_{1}\right) + C_{2} \left(T - T_{2}\right) = 0$



In [3]:

    
solve(eq, T)









    Out[3]:





$\displaystyle \left[ \frac{C_{1} T_{1} + C_{2} T_{2}}{C_{1} + C_{2}}\right]$

We can use SymPy to solve the cooling differential equation.



In [4]:

    
T_init, T_env, r, t = symbols('T_init T_env r t')
T = Function('T')

eqn = Eq(diff(T(t), t), -r * (T(t) - T_env))
eqn









    Out[4]:





$\displaystyle \frac{d}{d t} T{\left(t \right)} = - r \left(- T_{env} + T{\left(t \right)}\right)$

Here's the general solution:



In [5]:

    
solution_eq = dsolve(eqn)
solution_eq









    Out[5]:





$\displaystyle T{\left(t \right)} = C_{1} e^{- r t} + T_{env}$



In [6]:

    
general = solution_eq.rhs
general









    Out[6]:





$\displaystyle C_{1} e^{- r t} + T_{env}$

We can use the initial condition to solve for $C_1$. First we evaluate the general solution at $t=0$



In [7]:

    
at0 = general.subs(t, 0)
at0









    Out[7]:





$\displaystyle C_{1} + T_{env}$

Now we set $T(0) = T_{init}$ and solve for $C_1$



In [8]:

    
solutions = solve(Eq(at0, T_init), C1)
value_of_C1 = solutions[0]
value_of_C1









    Out[8]:





$\displaystyle - T_{env} + T_{init}$

Then we plug the result into the general solution to get the particular solution:



In [9]:

    
particular = general.subs(C1, value_of_C1)
particular









    Out[9]:





$\displaystyle T_{env} + \left(- T_{env} + T_{init}\right) e^{- r t}$

We use a similar process to estimate $r$ based on the observation $T(t_{end}) = T_{end}$



In [10]:

    
t_end, T_end = symbols('t_end T_end')

Here's the particular solution evaluated at $t_{end}$



In [11]:

    
at_end = particular.subs(t, t_end)
at_end









    Out[11]:





$\displaystyle T_{env} + \left(- T_{env} + T_{init}\right) e^{- r t_{end}}$

Now we set $T(t_{end}) = T_{end}$ and solve for $r$



In [12]:

    
solutions = solve(Eq(at_end, T_end), r)
value_of_r = solutions[0]
value_of_r









    Out[12]:





$\displaystyle \frac{\log{\left(\frac{- T_{env} + T_{init}}{T_{end} - T_{env}} \right)}}{t_{end}}$

We can use evalf to plug in numbers for the symbols. The result is a SymPy float, which we have to convert to a Python float.



In [13]:

    
subs = dict(t_end=30, T_end=70, T_init=90, T_env=22)
r_coffee2 = value_of_r.evalf(subs=subs)
type(r_coffee2)









    Out[13]:





sympy.core.numbers.Float



In [14]:

    
r_coffee2 = float(r_coffee2)
r_coffee2









    Out[14]:





$\displaystyle 0.011610223142273859$



In [ ]: