In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *
In [2]:
def update_func(state, t, system):
"""Update the thermal transfer model.
state: State (temp)
t: time
system: System object
returns: State (temp)
"""
r, T_env, dt = system.r, system.T_env, system.dt
T = state.T
T += -r * (T - T_env) * dt
return State(T=T)
In [3]:
def run_simulation(system, update_func):
"""Runs a simulation of the system.
Add a TimeFrame to the System: results
system: System object
update_func: function that updates state
"""
init = system.init
t_0, t_end, dt = system.t_0, system.t_end, system.dt
frame = TimeFrame(columns=init.index)
frame.row[t_0] = init
ts = linrange(t_0, t_end, dt)
for t in ts:
frame.row[t+dt] = update_func(frame.row[t], t, system)
return frame
In [4]:
def make_system(T_init, r, volume, t_end):
"""Makes a System object with the given parameters.
T_init: initial temperature in degC
r: heat transfer rate, in 1/min
volume: volume of liquid in mL
t_end: end time of simulation
returns: System object
"""
init = State(T=T_init)
return System(init=init,
r=r,
volume=volume,
temp=T_init,
t_0=0,
t_end=t_end,
dt=1,
T_env=22)
In [5]:
def func(x):
return (x-1) * (x-2) * (x-3)
modsim.py provides root_bisect, which searches for a root by bisection. The first argument is the function whose roots we want. The second argument is an interval that contains a root.
In [6]:
res = root_bisect(func, [0.5, 1.5])
Out[6]:
The result is an object that contains the root that was found and other information.
In [7]:
res.root
Out[7]:
If we provide a different interval, we find a different root.
In [8]:
res = root_bisect(func, [1.5, 2.5])
Out[8]:
In [9]:
res.root
Out[9]:
If the interval doesn't contain a root, the results explain the error.
In [10]:
res = root_bisect(func, [4, 5])
Out[10]:
We want to find the value of r that makes the final temperature 70, so we define an "error function" that takes r as a parameter and returns the difference between the final temperature and the goal.
In [11]:
def error_func1(r):
"""Runs a simulation and returns the `error`.
r: heat transfer rate, in 1/min
returns: difference between final temp and 70 C
"""
system = make_system(T_init=90, r=r, volume=300, t_end=30)
results = run_simulation(system, update_func)
T_final = get_last_value(results.T)
return T_final - 70
With r=0.01, we end up a little too warm.
In [12]:
error_func1(r=0.01)
Out[12]:
With r=0.02, we end up too cold.
In [13]:
error_func1(r=0.02)
Out[13]:
The return value from root_bisect is an array with a single element, the estimated value of r.
In [14]:
res = root_bisect(error_func1, [0.01, 0.02])
Out[14]:
In [15]:
r_coffee = res.root
Out[15]:
If we run the simulation with the estimated value of r, the final temperature is 70 C, as expected.
In [16]:
coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)
results = run_simulation(coffee, update_func)
T_final = get_last_value(results.T)
Out[16]:
Exercise: When you call root_bisect, it calls error_func1 several times. To see how this works, add a print statement to error_func1 and run root_bisect again.
Exercise: Repeat this process to estimate r_milk, given that it starts at 5 C and reaches 20 C after 15 minutes.
Before you use root_bisect, you might want to try a few values for r_milk and see how close you can get by trial and error. Here's an initial guess to get you started:
In [17]:
r_milk = 0.1
milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)
results = run_simulation(milk, update_func)
T_final = get_last_value(results.T)
Out[17]:
In [18]:
# Solution
def error_func2(r):
"""Runs a simulation and returns the `error`.
r: heat transfer rate, in 1/min
returns: difference between final temp and 20C
"""
system = make_system(T_init=5, r=r, volume=50, t_end=15)
results = run_simulation(system, update_func)
T_final = get_last_value(results.T)
return T_final - 20
In [19]:
# Solution
error_func2(r=0.1)
Out[19]:
In [20]:
# Solution
error_func2(r=0.2)
Out[20]:
In [21]:
# Solution
res = root_bisect(error_func2, [0.1, 0.2])
Out[21]:
In [22]:
# Solution
r_milk = res.root
Out[22]:
In [23]:
# Solution
milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)
results = run_simulation(milk, update_func)
T_final = get_last_value(results.T)
Out[23]:
The following function takes System objects that represent two liquids, computes the temperature of the mixture, and returns a new System object that represents the mixture.
In [24]:
def mix(system1, system2):
"""Simulates the mixture of two liquids.
system1: System representing coffee
system2: System representing milk
returns: System representing the mixture
"""
assert system1.t_end == system2.t_end
V1, V2 = system1.volume, system2.volume
T1, T2 = system1.temp, system2.temp
V_mix = V1 + V2
T_mix = (V1 * T1 + V2 * T2) / V_mix
return make_system(T_init=T_mix,
r=system1.r,
volume=V_mix,
t_end=30)
mix requires the System objects to have temp as a system variable. make_system initializes this variable;
the following function makes sure it gets updated when we run a simulation.
In [25]:
def run_and_set(system):
"""Run a simulation and set the final temperature.
system: System
returns: TimeFrame
"""
results = run_simulation(system, update_func)
system.temp = get_last_value(results.T)
return results
In [26]:
coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)
coffee.temp
Out[26]:
In [27]:
milk = make_system(T_init=5, r=r_milk, volume=50, t_end=30)
milk.temp
Out[27]:
In [28]:
mix_first = mix(coffee, milk)
mix_first.temp
Out[28]:
In [29]:
mix_results = run_and_set(mix_first)
mix_first.temp
Out[29]:
In [30]:
coffee_results = run_and_set(coffee)
coffee.temp
Out[30]:
In [31]:
milk_results = run_and_set(milk)
milk.temp
Out[31]:
Here's what the results look like.
In [32]:
plot(coffee_results.T, label='coffee')
plot(milk_results.T, '--', label='milk')
decorate(xlabel='Time (minutes)',
ylabel='Temperature (C)',
loc='center left')
savefig('figs/chap16-fig01.pdf')
Here's what happens when we mix them.
In [33]:
mix_last = mix(coffee, milk)
mix_last.temp
Out[33]:
In [34]:
mix_last.temp - mix_first.temp
Out[34]:
The following function takes t_add, which is the time when the milk is added, and returns the final temperature.
In [35]:
def run_and_mix(t_add, t_total):
"""Simulates two liquids and them mixes them at t_add.
t_add: time in minutes
t_total: total time to simulate, min
returns: final temperature
"""
coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=t_add)
coffee_results = run_and_set(coffee)
milk = make_system(T_init=5, r=r_milk, volume=50, t_end=t_add)
milk_results = run_and_set(milk)
mixture = mix(coffee, milk)
mixture.t_end = t_total - t_add
results = run_and_set(mixture)
return mixture.temp
We can try it out with a few values.
In [36]:
run_and_mix(t_add=0, t_total=30)
Out[36]:
In [37]:
run_and_mix(t_add=15, t_total=30)
Out[37]:
In [38]:
run_and_mix(t_add=30, t_total=30)
Out[38]:
And then sweep a range of values for t_add
In [39]:
sweep = SweepSeries()
for t_add in linspace(0, 30, 11):
sweep[t_add] = run_and_mix(t_add, 30)
Here's what the result looks like.
In [40]:
plot(sweep, label='final temp', color='C2')
decorate(xlabel='Time added (min)',
ylabel='Final temperature (C)')
savefig('figs/chap16-fig02.pdf')
Now we can use the analytic result to compute temperature as a function of time. The following function is similar to run_simulation.
In [41]:
def run_analysis(system):
"""Computes temperature using the analytic solution.
system: System object
returns: TimeFrame
"""
T_env, r = system.T_env, system.r
T_init = system.init.T
ts = linspace(0, system.t_end)
T_array = T_env + (T_init - T_env) * exp(-r * ts)
# to be consistent with run_simulation,
# we put the array into a TimeFrame
return TimeFrame(T_array, index=ts, columns=['T'])
Here's how we run it. From the analysis (see chap16sympy.ipynb), we have the computed value of r_coffee2
In [42]:
r_coffee2 = 0.011610223142273859
coffee2 = make_system(T_init=90, r=r_coffee2, volume=300, t_end=30)
Out[42]:
In [43]:
results = run_analysis(coffee2)
T_final_analysis = get_last_value(results.T)
Out[43]:
And we can compare to the results from simulation.
In [44]:
coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)
results = run_simulation(coffee, update_func)
T_final_simulation = get_last_value(results.T)
Out[44]:
They are identical except for a small roundoff error.
In [45]:
T_final_analysis - T_final_simulation
Out[45]:
Exercise: Suppose the coffee shop won't let me take milk in a separate container, but I keep a bottle of milk in the refrigerator at my office. In that case is it better to add the milk at the coffee shop, or wait until I get to the office?
Hint: Think about the simplest way to represent the behavior of a refrigerator in this model. The change you make to test this variation of the problem should be very small!
In [46]:
# Solution
## A refrigerator keeps the milk at a constant temperature,
## so it is like a container with r = 0.
## With r_milk = 0, it is best to add the milk at the beginning.
In [ ]: