Chapter 16

Click here to run this chapter on Colab

In the previous chapter we wrote a simulation of a cooling cup of coffee. Given the initial temperature of the coffee, the temperature of the atmosphere, and the rate parameter, r, we can predict how the temperature of the coffee will change over time. Then we used a root finding algorithm to estimate r based on data I invented for the example.

If you did the exercises, you simulated the temperature of the milk as it warmed, and estimated its rate parameter as well.

Now let’s put it together. In this chapter we’ll write a function that simulates mixing the two liquids and use it to answer the question we started with: is it better to mix the coffee and milk at the beginning, the end, or somewhere in the middle?

Mixing Liquids

When we mix two liquids, the temperature of the mixture depends on the temperatures of the ingredients as well as their volumes, densities, and specific heat capacities (as defined in the previous chapter). In this section I’ll explain how.

Assuming there are no chemical reactions that either produce or consume heat, the total thermal energy of the system is the same before and after mixing; in other words, thermal energy is conserved.

If the temperature of the first liquid is \(T_1\), the temperature of the second liquid is \(T_2\), and the final temperature of the mixture is \(T\), the heat transfer into the first liquid is \(C_1 (T - T_1)\) and the heat transfer into the second liquid is \(C_2 (T - T_2)\), where \(C_1\) and \(C_2\) are the thermal masses of the liquids.

In order to conserve energy, these heat transfers must add up to 0:

\[C_1 (T - T_1) + C_2 (T - T_2) = 0\]

We can solve this equation for T:

\[T = \frac{C_1 T_1 + C_2 T_2}{C_1 + C_2}\]

For the coffee cooling problem, we have the volume of each liquid; if we also know the density, \(\rho\), and the specific heat capacity, \(c_p\), we can compute thermal mass:

\[C = \rho V c_p\]

If the liquids have the same density and heat capacity, they drop out of the equation, and we can write:

\[T = \frac{V_1 T_1 + V_2 T_2}{V_1 + V_2}\]

where \(V_1\) and \(V_2\) are the volumes of the liquids.

As an approximation, I’ll assume that milk and coffee have the same density and specific heat. As an exercise, you can look up these quantities and see how good this assumption is.

The following function takes two System objects, representing the coffee and milk, and creates a new System to represent the mixture:

def mix(system1, system2):
    V1, V2 = system1.volume, system2.volume
    T1, T2 = system1.T_final, system2.T_final
    V_mix = V1 + V2
    T_mix = (V1 * T1 + V2 * T2) / V_mix
    return make_system(T_init=T_mix,

The first two lines extract volume and temperature from the System objects. The next two lines compute the volume and temperature of the mixture. Finally, mix makes a new System object and returns it.

This function uses the value of r from system1 as the value of r for the mixture. If system1 represents the coffee, and we are adding the milk to the coffee, this is probably a reasonable choice. On the other hand, when we increase the amount of liquid in the coffee cup, that might change r. So this is an assumption we might want to revisit.

Now we have everything we need to solve the problem.

Mix First Or Last?

First I’ll create objects to represent the coffee and milk. For r_coffee, I’ll use the value we computed in the previous chapter.

r_coffee = 0.0115
coffee = make_system(T_init=90, volume=300, r=r_coffee, t_end=30)

For r_milk, I’ll use the value I estimated in the exercise from the previous chapter.

r_milk = 0.133
milk = make_system(T_init=5, volume=50, r=r_milk, t_end=15)

Now we can mix them and simulate 30 minutes:

mix_first = mix(coffee, milk)
run_simulation(mix_first, change_func)


The final temperature is 61.5 °C which is still warm enough to be enjoyable. Would we do any better if we added the milk last?

I’ll simulate the coffee and milk separately, and then mix them:

run_simulation(coffee, change_func)
run_simulation(milk, change_func)
mix_last = mix(coffee, milk)

After mixing, the temperature is 62.9 °C, so it looks like adding the milk at the end is better. But is that the best we can do?


Adding the milk after 30 minutes is better than adding it immediately, but maybe there’s something in between that’s even better. To find out, I’ll use the following function, which takes the time to add the milk, t_add, as a parameter:

def run_and_mix(t_add, t_total):
    coffee.t_end = t_add
    coffee_results = run_simulation(coffee, change_func)
    milk.t_end = t_add
    milk_results = run_simulation(milk, change_func)
    mixture = mix(coffee, milk)
    mixture.t_end = t_total - t_add
    results = run_simulation(mixture, change_func)

    return mixture.T_final

run_and_mix simulates both systems for the given time, t_add. Then it mixes them and simulates the mixture fir the remaining time, t_total - t_add.

When t_add is0, we add the milk immediately; when t_add is 30, we add it at the end. Now we can sweep the range of values in between:

sweep = SweepSeries()
for t_add in linspace(0, 30, 11):
    sweep[t_add] = run_and_mix(t_add, 30)

Here’s what the results look like:

sweep.plot(label='mixture', color='C2')

decorate(xlabel='Time until mixing (minutes)',
         ylabel='Final emperature (C)')

Note that this is a parameter sweep, not a time series.

The final temperature is maximized when t_add=30, so adding the milk at the end is optimal.


Simulating Newton’s law of cooling isn’t really necessary because we can solve the differential equation analytically. If

\[\frac{dT}{dt} = -r (T - T_{env})\]

the general solution is

\[T{\left (t \right )} = C \exp(-r t) + T_{env}\]

and the particular solution where \(T(0) = T_{init}\) is

\[T_{env} + \left(- T_{env} + T_{init}\right) \exp(-r t)\]

If you would like to see this solution done by hand, you can watch this video:

Now we can use the observed data to estimate the parameter \(r\). If we observe the that temperature at \(t_{end}\) is \(T_{final}\), we can plug these values into the particular solution and solve for \(r\). The result is:

\[r = \frac{1}{t_{end}} \log{\left (\frac{T_{init} - T_{env}}{T_{final} - T_{env}} \right )}\]

The following function takes a System object and computes r:

from numpy import log

def compute_r(system):
    t_end = system.t_end
    T_init = system.T_init
    T_final = system.T_final
    T_env = system.T_env
    r = log((T_init - T_env) / (T_final - T_env)) / t_end
    return r

We can use this function to compute r for the coffee, given the parameters of the problem.

coffee2 = make_system(T_init=90, volume=300, r=0, t_end=30)
coffee2.T_final = 70
r_coffee2 = compute_r(coffee2)

This value is close to the value of r we computed in the previous chapter, 0.115, but not exactly the same. That’s because the simulations use discrete time steps, and the analysis uses continuous time.

Nevertheless, the results of the analysis are consistent with the simulation. To check, we’ll use the following function, which takes a System object and uses the analytic result to compute a time series:

from numpy import exp

def run_analysis(system):
    T_env, T_init, r = system.T_env, system.T_init, system.r
    t_array = linrange(system.t_0, system.t_end, system.dt)    
    T_array = T_env + (T_init - T_env) * exp(-r * t_array)
    system.T_final = T_array[-1]
    return make_series(t_array, T_array)

The first line unpacks the system variables. The next two lines compute t_array, which is a NumPy array of time stamps, and T_array, which is an array of the corresponding temperatures.

The last two lines store the final temperature in the System object and use make_series to return the results in a Pandas Series.

We can run it like this:

coffee2.r = r_coffee2
results2 = run_analysis(coffee2)

The final temperature is 70 °C, as it should be. In fact, the results are identical to what we got by simulation, with a small difference due to rounding.


In this chapter we finished the coffee cooling problem from the previous chapter, and found that it is better to add the milk at the end, at least for the version of the problem I posed.

As an exercise you will have a chance to explore a variation of the problem where the answer might be different.

In the next chapter we’ll move on to a new example, a model of how glucose and insulin interact to control blood sugar.



Use compute_r to compute r_milk according to the analytic solution. Run the analysis with this value of r_milk and confirm that the results are consistent with the simulation.

# Solution

milk2 = make_system(T_init=5, volume=50, r=0, t_end=15)
milk2.T_final = 20
r_milk2 = compute_r(milk2)
# Solution

milk2 = make_system(T_init=5, volume=50, r=r_milk2, t_end=15)
results2 = run_analysis(milk2)
# Solution

milk = make_system(T_init=5, volume=50, r=0.13296, t_end=15)
results = run_simulation(milk, change_func)
# Solution

from numpy import allclose

allclose(results, results2)


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!

# Solution

# A refrigerator keeps the milk at a constant temperature,
# so it is like a container with r = 0.

# With T_init=5 and r_milk = 0, 
# it is best to add the milk at the beginning.