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

And a System object to contain the system parameters.

In [3]:
coffee = System(init=init,

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)

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)

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)


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,

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)


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 goes here

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

In [ ]: