chap15soln


Modeling and Simulation in Python

Chapter 15

Copyright 2017 Allen Downey

License: Creative Commons Attribution 4.0 International


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 *

The coffee cooling problem

I'll use a State object to store the initial temperature.


In [2]:
init = State(T=90)


Out[2]:
values
T 90

And a System object to contain the system parameters.


In [3]:
coffee = System(init=init,
                volume=300,
                r=0.01,
                T_env=22,
                t_0=0,
                t_end=30,
                dt=1)


Out[3]:
values
init T 90 dtype: int64
volume 300
r 0.01
T_env 22
t_0 0
t_end 30
dt 1

The update function implements Newton's law of cooling.


In [4]:
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)

Here's how it works.


In [5]:
update_func(init, 0, coffee)


Out[5]:
values
T 89.32

Here's a version of run_simulation that uses linrange to make an array of time steps.


In [6]:
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

And here's how it works.


In [7]:
results = run_simulation(coffee, update_func)


Out[7]:
T
0 90
1 89.32
2 88.6468
3 87.9803
4 87.3205
5 86.6673
6 86.0207
7 85.3804
8 84.7466
9 84.1192
10 83.498
11 82.883
12 82.2742
13 81.6714
14 81.0747
15 80.484
16 79.8991
17 79.3201
18 78.7469
19 78.1795
20 77.6177
21 77.0615
22 76.5109
23 75.9658
24 75.4261
25 74.8919
26 74.3629
27 73.8393
28 73.3209
29 72.8077
30 72.2996

Here's what the results look like.


In [8]:
plot(results.T, label='coffee')
decorate(xlabel='Time (minutes)',
         ylabel='Temperature (C)')


And here's the final temperature:


In [9]:
coffee.T_final = get_last_value(results.T)
T_final = get_last_value(results.T)


Out[9]:
72.2996253904031

Encapsulation

Before we go on, let's define a function to initialize System objects with relevant parameters:


In [10]:
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)

Here's how we use it:


In [11]:
coffee = make_system(T_init=90, r=0.01, volume=300, t_end=30)
results = run_simulation(coffee, update_func)
T_final = get_last_value(results.T)


Out[11]:
72.2996253904031

Exercises

Exercise: Simulate the temperature of 50 mL of milk with a starting temperature of 5 degC, in a vessel with the same insulation, for 15 minutes, and plot the results.

By trial and error, find a value for r that makes the final temperature close to 20 C.


In [12]:
# Solution

milk = make_system(T_init=5, t_end=15, r=0.133, volume=50)
results = run_simulation(milk, update_func)
T_final = get_last_value(results.T)


Out[12]:
20.00135627897414

In [13]:
plot(results.T, label='milk')
decorate(xlabel='Time (minutes)',
         ylabel='Temperature (C)')



In [ ]: