Chapter 4

Building and Exploring Models

In the last two chapters, we learned how to find solutions to differential equations and how to understand their behavior without finding a formula. But where do the equations come from? If we begin with a question about the world, how do we decide which equation to write?

This is the task of mathematical modeling. We begin by identifying the quantities that seem relevant to our question. Then we decide which changes we need to follow and which quantities we can reasonably treat as constant. Using these ingredients and their derivatives, we try to turn an ordinary description of what is happening into a precise mathematical relationship.

There is an art to these choices. The same description may suggest several equations, and a quantity we treat as constant for one question may need to vary for another. Our aim is to build a model simple enough that we can understand it and make predictions with it, but detailed enough to capture the behavior that matters. Once we have an equation, mathematics tells us what it implies. We can then compare those consequences with observations and use what we learn to revise our choices.

We will begin by returning to Newton's law of cooling, this time slowing down to examine how we build the model. Then we will investigate population growth, letting measurements guide us from a simple model to a more complicated one, before adding harvesting. Elsewhere in the chapter, models of interacting populations, epidemics, and moving springs will ask us to follow more quantities and relate them in new ways. The chapter closes with a climate appendix, where we apply the same modeling approach to Earth's temperature. Ecology, climate, and mechanics will keep returning throughout the book as we develop the mathematics their models lead us to.

4.1Newton's Law of Cooling

Suppose we leave a hot cup of coffee on the table and want to know how warm it will be twenty minutes from now. What quantities might matter? The coffee's current temperature matters, of course. But so does the temperature of its surroundings: the same coffee would cool differently in a warm kitchen and a cold garage.

Let us begin with those two quantities: the temperature of the coffee and the temperature of the room. Which should we treat as variables, and which as constants?

Over twenty minutes, the coffee's temperature may change substantially, while the room's temperature changes very little. So a reasonable first choice is to let the coffee's temperature be a function $𝑇 (𝑡)$ , while treating the room's temperature as a constant $𝑇 r o o m$ . We are also simplifying by describing the coffee with just one temperature, ignoring differences between the liquid near the surface and the liquid deeper in the cup.

These are choices about our model. Over a single second, we might reasonably treat both temperatures as nearly constant. Over a whole day, the room may warm and cool enough that we need to follow its temperature too. We choose what to hold fixed by considering the time scale and the question we want to answer.

Now we have our ingredients. What can we say about how they are related? A cup hotter than the room cools down, a cup colder than the room warms up, and a larger temperature difference produces a faster change. The quantity

𝑇 - 𝑇 r o o m

therefore seems useful: its sign tells us which way the temperature should change, and its magnitude tells us how far the coffee is from matching its surroundings.

But “a larger difference produces a faster change” does not tell us exactly how much faster. We have to make a choice. One simple possibility is proportionality: doubling the temperature difference doubles the rate of change. With the sign chosen so that hot coffee cools, this gives

𝑇' = - 𝑘 (𝑇 - 𝑇 r o o m), 𝑘 > 0 .

Why do we need the proportionality constant $𝑘$ ? Even before considering any details of heat transfer, the units tell us that something is missing without it. The temperature difference is measured in degrees, while $𝑇'$ is measured in degrees per unit time. Suppose we change our time unit from seconds to hours. At the same physical instant, the temperature difference has the same numerical value, but the rate expressed in degrees per hour is $3600$ times the rate expressed in degrees per second! The numerical value of $𝑘$ must change by the same factor.

Thus $𝑘$ has units of inverse time. Its value also describes how quickly this particular object responds to its surroundings: a larger $𝑘$ gives a faster change for the same temperature difference. Our verbal observations have not determined that value; we would need measurements or further physical information. Treating $𝑘$ as constant throughout the cooling is another assumption of our simple model.

Let us check the equation against the observations that motivated it. When $𝑇 > 𝑇 r o o m$ , the right side is negative, so the coffee cools. When $𝑇 < 𝑇 r o o m$ , it is positive, so the coffee warms. And when the temperatures agree, the rate is zero. The equation says what we intended!

We have recovered Newton's law of cooling, which we already solved in chapter 2:

𝑇 (𝑡) = 𝑇 r o o m + (𝑇 (0) - 𝑇 r o o m) 𝑒 - 𝑘 𝑡 .

Once we have reason to trust this model---for example, because we have seen it agree with measurements in many situations---we can put it to work. But our formula still contains $𝑘$ . The cooling law itself does not tell us its value: that depends on the thermal properties of the object and its surroundings, as well as our choice of time unit.

Fortunately, we do not need to derive $𝑘$ from those physical details. We can measure our way to it! Suppose we know the constant room temperature and measure the hot coffee's temperature at time $0$ . Our formula then has only one undetermined constant left. Measuring the coffee's temperature once more, at a later time, gives an equation we can solve for $𝑘$ . With that value in hand, the formula predicts the temperature at times we have not measured.

This is a very useful feature of modeling. We began by wanting an entire function---the coffee's future temperature---but the model reduced what we needed to learn to finitely many numbers. A few measurements determine those numbers, and the equation supplies the rest of the history. In the problems, you will use two measurements of a cooling cup of coffee to predict a third.

Try this with the measurements in Figure figure 4.1. Here we have more than two observations, so we can also see how well one choice of $𝑘$ describes the whole cooling history.

Figure 4.1 The dots are reported coffee-temperature measurements from *Active Prelude to Calculus*. The initial temperature is $186 \circ F$ and the room temperature is $71 \circ F$ ; keep both fixed and adjust $𝑘$ to fit the curve. With $𝑘 = 0.030 m i n - 1$ , the model follows the main cooling trend, although it does not pass through every measurement.

Of course, determining the constant assumes that we have chosen a suitable model. Let us now turn to population growth, where measurements will help us both choose the constants and investigate whether the equation itself needs to change.

4.2Population Growth

How can we predict the size of a growing population? Let us begin with a population living where food and space are plentiful, and try to describe its growth with a single quantity: the amount of population present, $𝑁 (𝑡)$ .

4.2.1From Exponential to Logistic Growth

What determines how quickly this quantity changes? New individuals are born, while others die. For a first model, suppose each individual makes approximately the same average contribution to the population's net growth. Then twice as many individuals should produce twice as much net growth over the same short interval. As with cooling, a simple way to express this assumption is proportionality:

𝑁' = 𝑟 𝑁 .

The coefficient $𝑟$ measures the net growth per unit of population per unit time.

Over time, the conditions in which a population lives can change, and its growth rate can change with them. In the United States, the annual birth rate is estimated to have been about 50 births per 1,000 people around 1800; in 2024 it was 10.7. Over those centuries, family life, economic conditions, and expectations about having children changed enormously as the country shifted from a predominantly rural society toward an urban one. Births are only one contribution to population growth---deaths and migration matter too---but even this one contribution is clearly something we cannot treat as constant forever. (Historical estimates; 2024 measurements.)

But we may not need to describe centuries of change. Over a shorter window---perhaps a few hours for a yeast culture, or a decade for a human population---the conditions controlling growth may change relatively little while the population itself changes appreciably. This gives us a reason to try treating $𝑟$ as constant in our first simplified model. As with the coffee and its surroundings, we are choosing to follow one quantity while holding another approximately fixed.

The units work just as they did for cooling. Since $𝑁'$ has units of population per time and $𝑁$ has units of population, $𝑟$ must have units of inverse time. Its numerical value depends both on our time unit and on how quickly this particular population grows.

We already know the solution. If $𝑁 (𝑡 0) = 𝑁 0$ , then

𝑁 (𝑡) = 𝑁 0 𝑒 𝑟 (𝑡 - 𝑡 0) .

For $𝑟 > 0$ , our assumption predicts exponential growth. Once we measure the initial population, the only undetermined constant is $𝑟$ . We can therefore try what we did with coffee: use measurements to choose the coefficient, and see how well the resulting curve describes the population.

In 1913, Tor Carlson measured the amount of yeast growing in brewer's wort at hourly intervals. Let us imagine watching his experiment as it happens. We have seen the first six measurements, taken an hour apart, and want to use them to build a prediction.

The first measurement gives $𝑁 (1) = 9.6$ , so our exponential model becomes

𝑁 (𝑡) = 9.6 𝑒 𝑟 (𝑡 - 1) .

Just as with the coffee, the initial measurement fixes one constant, leaving a single parameter for us to determine. Here it is $𝑟$ , the growth coefficient. Changing $𝑟$ changes how quickly the curve rises, while every choice passes through the same measured starting point.

Figure 4.2 The dots show Carlson's first six measurements. Adjust $𝑟$ to fit the exponential curve to them, keeping $𝑁 (1) = 9.6$ fixed. With $𝑟 = 0.50 h - 1$ , the model predicts $117.0$ at hour 6; Carlson measured $119.1$ .

This looks promising! One constant describes the rapid growth quite well. We have done more than draw a curve through measurements: we have chosen a differential equation, used the measurements to determine its remaining parameter, and obtained a function we can evaluate at future times.

But those future values are predictions. To find out whether they are useful, we have to let the experiment continue.

Figure 4.3 Now reveal the remaining measurements. The amount of yeast levels off near $660$ , while our exponential curve continues upward. Try changing $𝑟$ : a value that matches the later measurements misses the early growth, and a value that matches the early growth misses the later measurements.

At hour $18$ , our model predicts

𝑁 (18) = 9.6 𝑒 0.50 (18 - 1) \approx 47, 182,

while Carlson measured $659.6$ . The prediction is about seventy-two times too large!

The same model that worked well over a few hours has failed over a longer interval. And adjusting its one constant cannot fix the problem. We need to reconsider the assumption that each unit of population keeps contributing the same net growth per unit time.

Can we use the measurements to check our assumption? We said that each unit of population contributes the same net growth per unit time. To investigate that contribution, divide the total growth rate by the amount of population present:

𝑁' 𝑁 .

This is called the per-capita growth rate. Our exponential model says that it equals the constant $𝑟$ . But what do the measurements say?

We do not have measurements of $𝑁'$ directly. We have the amount of yeast at particular times. Between two consecutive measurements, we can estimate the derivative using a difference quotient:

𝑁' \approx 𝑁 𝑛 + 1 - 𝑁 𝑛 Δ 𝑡 .

To estimate the growth per unit of population, we also need the amount of yeast present during that interval. It changes from $𝑁 𝑛$ to $𝑁 𝑛 + 1$ , so let us use their average. This gives

e s t i m a t e d p e r - c a p i t a g r o w t h r a t e = 𝑁 𝑛 + 1 - 𝑁 𝑛 (𝑁 𝑛 + 1 + 𝑁 𝑛 2) Δ 𝑡 .

We can compute this number for each pair of consecutive measurements, then plot it against the average amount of yeast present.

Figure 4.4 Each dot compares the average amount of yeast during an interval with its estimated per-capita growth rate. If our exponential model were accurate, the dots would lie near one horizontal level. Instead, the rate decreases as the amount of yeast increases. Try fitting a straight line to this relationship.

The data suggest a way forward! The per-capita rate is not constant, but its decrease looks approximately linear. We could try many curves through these points; a line is a simple next choice.

Let $𝑟$ denote the rate when the population is very small, and let $𝐾$ denote the population at which the line reaches zero. The line with these two intercepts is

𝑁' 𝑁 = 𝑟 (1 - 𝑁 𝐾) .

Multiplying by $𝑁$ , we obtain the logistic equation:

𝑁' = 𝑟 𝑁 (1 - 𝑁 𝐾) .

Notice what we have changed. The growth per unit of population now varies as the population grows. We still have constant parameters, $𝑟$ and $𝐾$ , but they describe a changing rate. The measurements have helped us replace one assumption with another. The parameter $𝑟$ still has units of inverse time, while $𝐾$ has the same units as $𝑁$ .

We have not established that the decrease must be linear. We have chosen a simple relationship suggested by the data. Now we need to investigate what this new equation predicts---and whether those predictions describe the experiment.

Before solving anything, let us see what the equation tells us:

𝑁' = 𝑟 𝑁 (1 - 𝑁 𝐾) .

When $𝑁$ is much smaller than $𝐾$ , the factor $1 - 𝑁 / 𝐾$ is nearly $1$ , so

𝑁' \approx 𝑟 𝑁 .

Our revised model therefore retains the early exponential growth that worked so well! The correction becomes important as the population grows. At $𝑁 = 𝐾$ , it makes the growth rate zero; above $𝐾$ , it makes the rate negative.

We can summarize these observations with a phase line:

0 ⟶ 𝐾 ⟵ .

Both $𝑁 = 0$ and $𝑁 = 𝐾$ are equilibrium solutions. A population between them grows, but uniqueness prevents its solution from crossing $𝐾$ . By the reasoning from chapter 3, it approaches $𝐾$ as time increases. A population starting above $𝐾$ decreases toward the same equilibrium. This is why we call $𝐾$ the carrying capacity of the model: it is the population level toward which every positive starting population tends.

So we already know something useful. Our new equation predicts growth that eventually levels off, just as Carlson's measurements do. But does it level off at the right height, and does it get there at the right pace?

The calculus from chapter 2 lets us make this comparison precisely. Away from the equilibria, we can separate variables:

𝑑 𝑁 𝑁 (1 - 𝑁 / 𝐾) = 𝑟 𝑑 𝑡 .

The partial-fraction identity

1 𝑁 (1 - 𝑁 / 𝐾) = 1 𝑁 + 1 𝐾 - 𝑁

lets us integrate the left side. For an initial population $𝑁 (𝑡 0) = 𝑁 0$ with $0 < 𝑁 0 < 𝐾$ , solving the resulting equation for $𝑁$ gives

𝑁 (𝑡) = 𝐾 1 + (𝐾 - 𝑁 0 𝑁 0) 𝑒 - 𝑟 (𝑡 - 𝑡 0) .

Now we can plot this solution against the measurements. As before, the first measurement fixes the initial condition. The two parameters left to choose are $𝑟$ and $𝐾$ , and the per-capita plot has already suggested values for both.

Figure 4.5 The logistic solution starts at Carlson's first measurement, $𝑁 (1) = 9.6$ . With $𝑟 = 0.54 h - 1$ and $𝐾 = 663$ , it follows the measurements through both the early growth and the later plateau. Adjust $𝑟$ and $𝐾$ yourself. The approximate values suggested by the per-capita plot, $𝑟 \approx 0.53 h - 1$ and $𝐾 \approx 666$ , also give a close fit.

This is worth celebrating! We began with a simple assumption, watched its prediction fail, and used the measurements to suggest a revision. Solving the revised equation then produced a curve that follows the whole experiment.

The logistic model is only a little more complicated than the exponential model: we added one parameter and let the growth per unit of population decrease linearly. Yet that change captures both the early rapid growth and the later slowdown. We have found a model that is still simple enough to understand, but now describes the important behavior over a much longer interval.

We built the logistic model by examining one yeast culture. How much of its success carries over to other populations?

In the 1930s, Georgy Gause measured populations of the single-celled organism Paramecium growing in laboratory cultures. His measurements fluctuate more than Carlson's yeast record, but we can recognize the same broad pattern: rapid growth at first, followed by a slowdown toward a persistent population level.

Try fitting the logistic model to these measurements. We need different values of $𝑟$ and $𝐾$ , but the equation itself stays the same. It does not reproduce every fluctuation. It does capture the larger pattern around which those fluctuations occur.

Now make a much larger leap. Instead of a laboratory culture observed over days, consider the population of the United States, recorded by the census every ten years. These measurements span more than two centuries, during which the country's territory, economy, medicine, and ways of life changed enormously. Can the same simple equation describe anything useful here?

Figure 4.6 Choose among Carlson's yeast, Gause's two *Paramecium* cultures, and the U.S. census from 1790 through 2020. For each record, the first measurement fixes the initial population; adjust $𝑟$ and $𝐾$ to fit the remaining observations. Notice the different population and time units. How much of each record can you describe with just these two parameters?

The census comparison is striking. A logistic curve captures much of the population's long sweep, despite the enormous complexity of the history behind those measurements. This is the kind of success we are looking for in a simple model: a small amount of mathematics organizes a great deal of observed behavior.

But interpreting its parameters requires judgment. In the laboratory cultures, the measurements show a recognizable plateau. The census record has not settled at a population ceiling, so its fitted $𝐾$ is an extrapolation beyond what we have observed.

And remember our earlier discussion of changing birth rates. The logistic model allows growth per unit of population to change, but it describes that change using only the population itself. A successful fit does not tell us that crowding explains changes in family size, or that immigration and technology can be ignored when predicting the future.

We can learn something from the broad pattern without claiming to have explained every cause. Whether that is enough depends on the question we want the model to answer.

4.2.2Harvesting

So far, we have asked how a population grows. Now suppose we want to use that population: we have a trout pond, and we would like to catch fish. How much can we harvest without eventually emptying the pond?

Let $𝑃 (𝑡)$ denote the fish population. Suppose its natural growth is described by the logistic model,

𝑃' = 𝑟 𝑃 (1 - 𝑃 𝐾) .

Fishing introduces another process that changes $𝑃$ . To include it, we need to decide how the harvesting works. Let us suppose fish are removed at a constant rate of $𝐻$ fish per year, regardless of how many fish are currently in the pond. This contributes a negative term:

𝑃' = 𝑟 𝑃 (1 - 𝑃 𝐾) - 𝐻 .

For a concrete example, suppose the pond has a carrying capacity of 200 fish and a growth coefficient of 1 per year. Including harvesting gives

𝑃' = 𝑃 (1 - 𝑃 200) - 𝐻,

where time is measured in years and $𝐻$ is the number of fish removed per year.

The population stays constant when natural growth exactly replaces the fish we catch:

𝑃 (1 - 𝑃 200) = 𝐻 .

We can find these equilibria graphically. Plot the natural growth rate against the population: it is a parabola. The constant harvesting rate is a horizontal line, and their intersections are the populations where growth and harvesting balance.

Where the parabola lies above the line, the population grows; where it lies below, the population declines. So this one picture also tells us how to draw the phase-line arrows!

Figure 4.7 Compare the pond's natural growth rate with the constant harvesting rate. Their intersections give the equilibria, and the signs of their difference determine the phase-line arrows and solution histories. Increase $𝐻$ and watch how the possible futures change.

For $𝐻 = 32$ fish per year, the two rates balance at

𝑃 = 40 a n d 𝑃 = 160 .

If the population falls below $40$ , its natural growth cannot replace the fish we remove, and it keeps declining. Between $40$ and $160$ , growth exceeds harvesting, so the population increases. Above $160$ , harvesting again exceeds natural growth, and the population decreases. Thus $40$ is an unstable threshold, while $160$ is an attracting equilibrium.

The same harvest can therefore have very different consequences depending on how many fish we begin with!

What happens if we increase the quota? The horizontal line rises, and the two intersections approach one another. The growth curve reaches its maximum at $𝑃 = 100$ , where natural growth adds 50 fish per year.

At $𝐻 = 50$ , the two equilibria meet. Our equation becomes

𝑃' = - (𝑃 - 100) 2 200 .

A population above $100$ decreases toward the equilibrium, but one below $100$ keeps declining. If we raise the harvest above $50$ , the pond cannot replace the catch at any population: there are no positive equilibria.

The model has helped us identify a maximum sustainable harvesting rate, but it has also shown us why operating at that maximum is precarious. At $𝐻 = 50$ , even a small drop below the equilibrium leads to continued decline. Errors in our population estimate or growth model, or an unusually poor year, could matter enormously. A harvesting policy needs room for that uncertainty.

There is another limitation we can read directly from the equation. At $𝑃 = 0$ , it still says $𝑃' = - 𝐻$ : our constant quota keeps removing fish from an empty pond! We must stop interpreting the solution as a population when it reaches zero.

We could instead model harvesting with a term such as $- 𝑞 𝑃$ , removing fish at a rate proportional to the population. That rate would automatically approach zero as the pond empties. But it would describe a different harvesting policy and give us a different equation to investigate.

Once again, choosing the term is part of choosing the model. One constant subtraction created a population threshold and changed which long-term outcomes were possible.

4.3Modeling with Systems

So far, we have built models around a single changing quantity: the temperature of an object or the size of a population. But sometimes predicting one quantity requires us to follow another. An object changes the temperature of the room around it, which in turn affects how the object cools. Predators change the population of their prey, while the available prey affects how the predators grow.

In these situations, we need several unknown functions and equations describing how they change together. We are naturally led to a system of differential equations. Let us return to temperature and population to see how this happens.

4.3.1The Coffee and the Room

When we modeled a cooling cup of coffee, we treated the room's temperature as constant. But the hot object also warms its surroundings. What if that change is large enough to matter?

Let $𝑇 (𝑡)$ be the object's temperature and $𝑅 (𝑡)$ the room's temperature. We will still describe each with a single temperature, but now both are unknown functions of time. Our rule for the object stays the same:

𝑇' = - 𝑘 (𝑇 - 𝑅), 𝑘 > 0 .

We need a new rule for $𝑅'$ . For this first model, suppose the room's temperature changes only because of its interaction with the object. When the object is hotter than the room, the room warms; when the object is colder, the room cools. A larger temperature difference should produce a faster change. Just as before, let us choose proportionality as a simple way to make that description precise:

𝑅' = 𝑐 (𝑇 - 𝑅), 𝑐 > 0 .

The new constant $𝑐$ describes how strongly the room responds. There is no reason for it to equal $𝑘$ : the object and the room need not change temperature equally quickly. In fact, for a cup of coffee in a large room, we would expect $𝑐$ to be much smaller than $𝑘$ : the coffee can cool substantially while the room barely warms.

Putting the two rules together gives our system:

𝑇' = - 𝑘 (𝑇 - 𝑅), 𝑅' = 𝑐 (𝑇 - 𝑅) .

If we set $𝑐 = 0$ , then $𝑅' = 0$ : the room stays at its initial temperature, and we recover our original cooling model exactly!

For $𝑐 > 0$ , the equations must be studied together. To know how quickly the object cools, we need the room's current temperature. But that temperature is itself changing in response to the object! A solution now consists of two functions, $𝑇 (𝑡)$ and $𝑅 (𝑡)$ , satisfying their equations together. To begin a prediction, we need both initial temperatures.

We have not yet developed systematic methods for finding formulas for solutions of systems like this. Learning how to solve and understand systems will be one of the main topics of the rest of the book. But we can already explore them numerically. Euler's method asks us to compute both rates from the current temperatures, then use those rates to advance both temperatures by one short time step. We must compute both rates before updating either temperature, so they refer to the same moment.

Figure 4.8 The object begins at $80 \circ C$ and the room at $20 \circ C$ . The solid curve follows the object and the dashed curve follows the room. Begin with $𝑐 = 0$ , then increase $𝑐 / 𝑘$ to make the room respond more strongly. Compare $𝑐 / 𝑘 = 0$ , $0.05$ , $0.25$ , and $0.5$ , keeping $𝑘 = 0.05 m i n - 1$ . These are illustrative numerical solutions, not fitted temperature measurements.

The two temperatures approach agreement. When $𝑐$ is small compared with $𝑘$ , they meet near the room's initial temperature. Increasing $𝑐 / 𝑘$ makes the room warm more, and the two temperatures approach a higher common value. Try changing $𝑘$ while keeping $𝑐 / 𝑘$ fixed: the histories unfold at a different speed, but approach the same temperature.

Allowing one more quantity to change gave us a different prediction and two equations to study together. Let us now see how interactions between populations lead to the same mathematical need.

4.3.2Interacting Populations

Suppose we want to predict how a population of hares will change. Earlier, we tried to do this using only the number of hares. But if lynx are hunting them, that number alone is no longer enough. The same hare population might grow when there are few lynx and decline when there are many. And the lynx population depends on how much food the hares provide. As with the coffee and the room, predicting either quantity requires us to follow both.

Let $𝐻 (𝑡)$ and $𝐿 (𝑡)$ be the hare and lynx populations. Before putting them together, let us consider how each might change on its own.

In the absence of lynx, suppose the hares have plenty of food and space. Our earlier assumption of a constant average contribution to net growth gives

𝐻' = 𝑎 𝐻, 𝑎 > 0 .

We recognize exponential growth. This is our first approximation to what the hares would do without predators.

What should happen to the lynx without hares to eat? We expect their population to decline, but that does not mean every lynx dies at the same instant. For a simple model, suppose that over each short interval, approximately the same fraction of the surviving lynx dies. Twice as many remaining lynx then means twice as many deaths per unit time. This gives

𝐿' = - 𝑐 𝐿, 𝑐 > 0,

and hence exponential decline, $𝐿 (𝑡) = 𝐿 (0) 𝑒 - 𝑐 𝑡$ . A constant death rate per individual is a modeling choice; it does not describe every detail of starvation. Like our other continuous population models, it also stops being a useful description when individual animals matter.

Now let the two populations interact.

What about encounters? If the populations mix uniformly, doubling the number of hares should double the opportunities for any lynx to find one, and doubling the number of lynx should double the total opportunities again. The simplest expression with both properties is the product $𝐻 𝐿$ . Let $𝑏 𝐻 𝐿$ be the rate at which encounters remove hares, and let $𝑑 𝐻 𝐿$ be the rate at which food from those encounters contributes to lynx growth. Then

𝐻' = 𝑎 𝐻 - 𝑏 𝐻 𝐿, 𝐿' = - 𝑐 𝐿 + 𝑑 𝐻 𝐿 .

This is the Lotka--Volterra predator--prey model. The coefficients $𝑏$ and $𝑑$ need not be equal: one successful encounter removes one hare, but its contribution to predator reproduction depends on how food is converted into new lynx.

The limiting cases check the bookkeeping. If $𝐿 = 0$ , then $𝐻' = 𝑎 𝐻$ and hares grow exponentially. If $𝐻 = 0$ , then $𝐿' = - 𝑐 𝐿$ and lynx decline exponentially. If either population is zero, the interaction term vanishes. The signs of all four terms match the verbal story.

We have built a system because the populations affect one another. A solution consists of both histories, $𝐻 (𝑡)$ and $𝐿 (𝑡)$ , and a prediction needs both starting populations.

We do not need to fit a particular ecosystem in order to explore one consequence of the model. Choose simple population and time units, take $𝑎 = 𝑏 = 1$ and $𝑐 = 𝑑 = 𝛾$ . This gives the illustrative family

𝐻' = 𝐻 (1 - 𝐿), 𝐿' = 𝛾 𝐿 (𝐻 - 1),

Here $𝛾$ controls how quickly the predator population responds relative to the prey population. Let us compute both histories and see what the interaction produces.

Figure 4.9 For the illustrative family $𝐻' = 𝐻 (1 - 𝐿)$ , $𝐿' = 𝛾 𝐿 (𝐻 - 1)$ , follow the solid hare history and the dashed lynx history. Both populations repeatedly rise and fall! Change the relative predator response $𝛾$ and compare the timing of their peaks.

The populations oscillate! When hares are abundant, the lynx have enough food for their population to grow. More lynx then put more pressure on the hares, whose population declines. With less food available, the lynx begin to decline too. Fewer predators give the hares a chance to recover, and the cycle begins again. In the computed histories, the hare peak comes before the lynx peak: food becomes abundant first, and the predator population responds afterward.

Remember the conclusion from chapter 3: a nonconstant solution of a single smooth autonomous equation $𝑁' = 𝑓 (𝑁)$ must be increasing or decreasing. It cannot keep rising and falling periodically. But knowing $𝐻$ alone no longer determines $𝐻'$ : the same number of hares can be growing or declining, depending on the lynx population. By following two interacting quantities, we have opened up a new kind of behavior.

We can explore these oscillations numerically now. Explaining their structure mathematically will be part of our later study of systems.

This is a prediction worth taking back to evidence. The next figure compares the Hudson Bay hare--lynx fur record, a laboratory Paramecium--Didinium experiment, and a prey--predatory-mite experiment. Each population has been divided by its own observed maximum, so we are comparing relative timing, not fitting one scaled model to all three records.

Figure 4.10 Does the prey peak first outside our computed example? Compare the Hudson Bay hare--lynx record, a Paramecium--Didinium chemostat, and Huffaker's mite experiment. Each series is divided by its own observed maximum, so we are comparing timing rather than abundance or fitted parameters. The Hudson Bay values are pelts rather than population counts, and the collection method changed in 1903. The Huffaker series is a secondary transcription included only for qualitative comparison.

The resemblance is striking. Yet fur returns depend on trapping as well as abundance, and normalization reveals timing while hiding scale. The evidence supports a qualitative mechanism; it does not establish one four-parameter model for every predator--prey system.

In Chapter 3, we sometimes compressed a solution picture onto a phase line by keeping the current value and leaving out time. We can do something similar here, but the current state consists of two populations.

Return to our illustrative model

𝐻' = 𝐻 (1 - 𝐿), 𝐿' = 𝛾 𝐿 (𝐻 - 1) .

At any time $𝑡$ , the two heights on the solution graphs give us one point:

(𝐻 (𝑡), 𝐿 (𝑡)) .

Put the hare population on the horizontal axis and the lynx population on the vertical axis. As time passes, this point moves through the plane, tracing a parametric curve. We call this plane the phase plane and the curve a trajectory.

The differential equations tell us how the point moves. Its horizontal velocity is $𝐻'$ , and its vertical velocity is $𝐿'$ . For example, at $(𝐻, 𝐿) = (1.5, 0.7)$ ,

𝐻' = 1.5 (1 - 0.7) = 0.45, 𝐿' = 𝛾 (0.7) (1.5 - 1) = 0.35 𝛾 .

Both populations are increasing, so the point moves upward and to the right. We can repeat this calculation at other points to draw a field of arrows.

Our equations are autonomous: their rates depend on the current populations, without explicitly depending on time. This is what lets us assign one fixed arrow to each point in the plane. Whenever the populations have those values, the equations prescribe the same velocity.

Figure 4.11 The histories on the left and the trajectory on the right show the same numerical solution. The horizontal coordinate of the moving point is the current hare population; its vertical coordinate is the current lynx population. Follow the point through the population cycle, or drag either picture to inspect a particular instant. Changing $𝛾$ changes both views.

The arrows and trajectories together form a phase portrait. Here the numerical trajectory appears to make a closed loop: the populations return together to their earlier values as the cycle repeats. We have observed this in a computation; we have not yet proved that the trajectory closes.

This curve is not generally the graph of lynx population as a function of hare population. The same number of hares can occur with different numbers of lynx, at different stages of the cycle. Nor does the curve alone tell us when the populations reach a particular point. The time histories retain that information; the phase-plane picture makes their motion together easier to see.

Limited food for the prey

We can now see how each assumption entered the equations: hare reproduction contributed $𝑎 𝐻$ , predator decline contributed $- 𝑐 𝐿$ , and encounters contributed the two $𝐻 𝐿$ terms. That also tells us where to revise the model. Let us try two extensions, beginning with the hares' own food supply.

Without lynx, our model says the hares grow exponentially forever. But hares need grass and other plants to eat, and the available land cannot support unlimited growth. We have already built a simple model for this kind of limitation! Replace the exponential growth term with logistic growth:

𝐻' = 𝑎 𝐻 (1 - 𝐻 𝐾) - 𝑏 𝐻 𝐿, 𝐿' = - 𝑐 𝐿 + 𝑑 𝐻 𝐿 .

Here $𝐾$ describes the carrying capacity of the hare population in the absence of lynx. The terms describing encounters stay the same. We have changed one assumption about how the hares grow on their own.

How much difference does this make? If $𝐻$ stays much smaller than $𝐾$ , the factor $1 - 𝐻 / 𝐾$ stays close to $1$ , so the correction is small. We might expect our original model to give much the same prediction. Let us compare them using the same initial populations and the same coefficients. Using the scaled populations from our earlier example gives

𝐻' = 𝐻 (1 - 𝐻 𝐾) - 𝐻 𝐿, 𝐿' = 𝛾 𝐿 (𝐻 - 1) .

Figure 4.12 Compare the original model with the logistic-prey revision, both starting at $𝐻 (0) = 1.5$ and $𝐿 (0) = 0.7$ , with $𝛾 = 1$ . Both populations share the same axes: blue for hares and red for lynx. The dashed curves show the original model; the solid curves include limited food for the hares. Start with $𝐾 = 100$ and the first $20$ time units, then reduce $𝐾$ or extend the view to $300$ . These are illustrative model comparisons, not fits to new measurements.

With a large carrying capacity, the hares are far from using up all the grass available to them. Over the first few cycles, the original model does a pretty good job! The extra term makes little difference to the behavior we see. For smaller carrying capacities, the peaks in the revised model visibly shrink as time passes. The cycles no longer repeat in the same way.

Let us follow that change in the phase plane. A repeating population cycle would trace the same closed loop each time. What happens to the trajectory when the successive peaks get smaller?

Figure 4.13 Follow the revised model over many cycles, starting from the same populations as above. The dashed loop is the original model; the blue trajectory includes limited food, and the arrows show the revised rates. Gold highlights the last $20$ time units and the endpoint. Begin with $𝐾 = 20$ , $𝛾 = 1$ , and $300$ time units. Adjust $𝐾$ , the predator rate $𝛾$ , or the elapsed time to see how the trajectory develops.

The trajectory winds inward: each trip around the center corresponds to a smaller swing in the populations. The computation suggests that both populations eventually settle at constant values. Try increasing $𝐾$ . The trajectory stays close to the original loop for more turns before the inward drift becomes noticeable.

We did not need to account for grass to explain why predators and prey can oscillate. But including it changes our prediction about whether those oscillations persist. If food is plentiful and we want to follow just a few cycles, the simpler model may be enough. If we want to predict what happens over many cycles, even a small limitation on the hares' food supply can matter. How far ahead we want to predict is part of choosing the model!

Hunting interacting populations

For our second extension, return to the original predator–prey model and add hunting. The trout pond taught us that a harvesting term is a modeling decision, and that removing a fixed fraction of the population is the policy which shuts itself off as the population disappears. So suppose hunters remove hares at a per-capita rate $𝑒$ and lynx at a per-capita rate $𝑓$ . We will consider $0 \leq 𝑒 < 1$ and $𝑓 \geq 0$ , so the hares can still grow when lynx are absent. Taking $𝛾 = 1$ , the two new terms are $- 𝑒 𝐻$ and $- 𝑓 𝐿$ :

𝐻' = 𝐻 (1 - 𝐿) - 𝑒 𝐻 = 𝐻 (1 - 𝑒 - 𝐿), 𝐿' = 𝐿 (𝐻 - 1) - 𝑓 𝐿 = 𝐿 (𝐻 - 1 - 𝑓) .

Before computing anything, guess. Hunting hares should mean fewer hares, and hunting lynx should mean fewer lynx. Now look for the balance point, the state at which neither population changes. When both populations are positive, the hare rate vanishes only when $𝐿 = 1 - 𝑒$ , and the lynx rate vanishes only when $𝐻 = 1 + 𝑓$ . So the balance sits at

(𝐻, 𝐿) = (1 + 𝑓, 1 - 𝑒) .

Read this carefully, because it contradicts the guess. The hare balance depends only on how hard we hunt lynx, and the lynx balance depends only on how hard we hunt hares. Hunting hares leaves the hare balance exactly where it was and lowers the lynx balance instead.

Figure 4.14 Hunt both populations at fixed per-capita rates $𝑒$ and $𝑓$ . The histories still rise and fall, but now about the dashed balance levels $𝐻 = 1 + 𝑓$ and $𝐿 = 1 - 𝑒$ , and the trajectory loops around the new balance point rather than the hollow unhunted balance at $(1, 1)$ . Raise the hare hunting rate and watch which dashed level moves.

The figure shows why the guess was wrong. Hunting hares does remove hares, but every hare we take is one the lynx do not eat. The hares can only hold their own once the lynx have declined to $1 - 𝑒$ , so it is the lynx that have the lower balance. Hunting lynx has the mirror effect: a hunted lynx population needs more food to hold steady, so the lynx are balanced only when the hares reach $1 + 𝑓$ , and it is the hares that have the higher balance. The cycles persist throughout, with each population spending part of every cycle above its balance and part below. What the hunting rates control is the level the cycles surround.

Nothing about the individual terms announced this. The term $- 𝑒 𝐻$ is as simple as a term can be, and on the hares alone its effect would be just what we guessed. Its effect on the coupled system became visible only when we balanced both rates at once. This is essentially the observation that first led Volterra to the model: when fishing in the Adriatic fell during the First World War, the share of predatory fish in the catch rose. The model captures that direction of change. Easing the harvest raises the predator balance and lowers the prey balance.

Learning rates from an outbreak

An epidemic gives us another reason to follow several populations. The number of infectious people matters, but it is not enough: the same number can spread an illness very differently depending on how many people are still able to catch it. Let $𝑆 (𝑡)$ count these susceptible people, $𝐼 (𝑡)$ those who are infectious and can transmit the illness, and $𝑅 (𝑡)$ those who have recovered from it. Here we use numbers of people rather than the fractions plotted in Chapter 1. For a first model, suppose recovered people cannot catch the illness again during this outbreak, and neglect births, deaths, and people entering or leaving the population we are following.

How quickly do new infections occur? We can reuse our reasoning about predator--prey encounters. If people mix uniformly, doubling the susceptible population should double an infectious person's opportunities to pass on the illness. Doubling the infectious population should double the total opportunities again. This suggests an infection rate $𝑏 𝑆 𝐼$ , where the positive constant $𝑏$ accounts for how often people meet and how readily those encounters transmit the illness.

For recovery, suppose approximately the same fraction of the infectious population recovers over each short interval of a given length. Twice as many infectious people then means twice as many recoveries per unit time, giving a recovery rate $𝛾 𝐼$ for another positive constant $𝛾$ .

Now follow where people go. Each new infection removes one person from $𝑆$ and adds one to $𝐼$ . Each recovery removes one from $𝐼$ and adds one to $𝑅$ . So the susceptible population only loses people, the infectious population both gains and loses them, and the recovered population only gains:

𝑆' = - 𝑏 𝑆 𝐼, 𝐼' = 𝑏 𝑆 𝐼 - 𝛾 𝐼, 𝑅' = 𝛾 𝐼 .

These are the SIR equations. The same $𝑏 𝑆 𝐼$ appears with opposite signs in the first two equations because an infection transfers a person between the groups; the recovery terms work the same way. We have three equations for three unknown histories.

As with the cooling constant, our assumptions have given us a model without determining its constants. We can learn about $𝑏$ and $𝛾$ by choosing values, computing the resulting histories, and comparing them with measurements.

In October 1967, a common-cold outbreak on Tristan da Cunha was recorded through daily reports of illness and recovery. A published table of 21 daily observations gives us two curves to compare with our model. We use the number currently ill as an approximation to $𝐼$ , and the number recovered from this outbreak as $𝑅$ . This assumes that the period with symptoms approximates the period when someone can transmit the illness.

There is another unknown here: how many people were susceptible at the start? We must estimate $𝑆 (0)$ as well, keeping the observed $𝐼 (0) = 1$ and $𝑅 (0) = 0$ . Take $𝑡 = 0$ at the first observation and measure time in days. For each choice of the three unknown values, a numerical solution gives us both histories to compare with the data.

Figure 4.15 Dots show the recorded illness and recovery counts; curves show numerical solutions of the SIR equations. Adjust the transmission rate $𝑏$ , recovery rate $𝛾$ , and initial susceptible count $𝑆 (0)$ to compare both histories. The fitted values $𝑏 \approx 0.0204 p e r s o n - 1 d a y - 1$ , $𝛾 \approx 0.267 d a y - 1$ , and $𝑆 (0) \approx 40.3$ were chosen by numerically minimizing the sum of squared discrepancies in both curves.

The model captures the rise and fall of illness along with the accumulating recoveries. Try changing one parameter at a time. Increasing $𝑏$ makes new infections occur faster at the same values of $𝑆$ and $𝐼$ , while increasing $𝛾$ moves people out of the infectious group faster. Changing $𝑆 (0)$ changes how many people can catch the illness in the first place. Each choice affects both curves, so matching the illness history alone is not enough: the recovery measurements give us another way to check it.

The equations can also suggest a more direct way to estimate a constant. In the exercises, we use $𝑅' = 𝛾 𝐼$ and the area under the illness curve to obtain $𝛾 \approx 0.28$ per day. This is close to the value found by fitting both curves, though the two procedures need not give exactly the same answer.

There is still room to improve the model. The fit misses the sharp illness peak and predicts a slower decline. Its estimate of about forty initially susceptible people is also much smaller than the island's population: the fit does not tell us whether the others were immune or whether uniform mixing was a poor assumption. These are estimates within our model, not direct measurements of its mechanisms.

We could revise those assumptions too. An SEIR model adds an exposed group between susceptible and infectious, allowing time between catching an illness and being able to transmit it. The published comparison considers this extension using the same data.

From the coffee and the room to interacting populations, the need for systems has come from asking what else we must keep track of. Once those quantities affect one another, their equations belong together. We can already compute approximate solutions and compare them with observations; learning how to solve and understand these systems will be one of the main tasks of the rest of the book.

4.4Forces and Second-Order Equations

Pull a mass attached to a horizontal spring away from its resting position and let go. It moves back and forth. How could we build a model of that motion?

Let $𝑥 (𝑡)$ be the displacement from the resting position, positive to the right and negative to the left. The mass $𝑚$ and the spring's stiffness will matter too; we will treat both as constant. For a first model, suppose the mass slides without friction and ignore air resistance. We want to predict $𝑥 (𝑡)$ from how we start the motion.

4.4.1From a restoring force to an equation

The spring pulls the mass back toward its resting position. When $𝑥 > 0$ , the force points left; when $𝑥 < 0$ , it points right. A greater displacement also produces a stronger pull. As with cooling, we need to choose a precise relationship. Let us try proportionality: doubling the displacement doubles the force. Then

𝐹 s p r i n g = - 𝑘 𝑥, 𝑘 > 0 .

The constant $𝑘$ measures the stiffness, and the minus sign makes the force oppose the displacement. This is Hooke's law, a useful approximation for many springs over a limited range of stretching and compression.

What does this force tell us about $𝑥$ ? Recall that $𝑥'$ is velocity: it says how quickly the position changes, including the direction of motion. Its derivative $𝑥 ″$ is acceleration, the rate at which velocity changes. Newton's second law says that total force equals mass times acceleration. Since the spring supplies the only horizontal force in our model,

𝑚 𝑥 ″ = - 𝑘 𝑥 .

Our description has led us to a second-order differential equation! The force gives us a rule for the second derivative. In particular, a mass to the right of its resting position has negative acceleration. It might still be moving right, but the spring is slowing that motion and will pull it back. The same position can occur on the outward trip and on the return trip, with different velocities.

This is why knowing the initial position alone cannot determine the motion. We must also specify how we launch the mass:

𝑥 (0) = 𝑥 0, 𝑥' (0) = 𝑣 0 .

These are our two initial conditions. As with interacting populations, we have found that predicting one quantity requires more information than its current value.

4.4.2What motion does the model predict?

Dividing by $𝑚$ gives $𝑥 ″ = - (𝑘 / 𝑚) 𝑥$ . We have seen this kind of equation in chapter 2: taking two derivatives of a sine or cosine gives a negative multiple of the original function. In particular,

𝑑 2 𝑑 𝑡 2 c o s (𝜔 𝑡) = - 𝜔 2 c o s (𝜔 𝑡),

so choosing $𝜔 = \sqrt 𝑘 / 𝑚$ gives the multiplier we need. The sine works the same way, and we can check that

𝑥 (𝑡) = 𝑥 0 c o s (𝜔 𝑡) + 𝑣 0 𝜔 s i n (𝜔 𝑡), 𝜔 = \sqrt 𝑘 𝑚,

satisfies both the equation and our initial conditions: at $𝑡 = 0$ , its value is $𝑥 0$ and its derivative is $𝑣 0$ .

Figure 4.16 For $𝑥 ″ = - 𝑥$ , both masses begin at the same position, but one has velocity $+ 𝑣 0$ and the other $- 𝑣 0$ . Choose the common position and velocity magnitude, then play the motions. Their histories leave the same point with opposite slopes and follow different paths.

The model predicts oscillation. For a nonzero motion, the time to complete one cycle is

T = 2 𝜋 𝜔 = 2 𝜋 \sqrt 𝑚 𝑘 .

A heavier mass oscillates more slowly; a stiffer spring makes it oscillate more quickly. We can also work backward: measuring the period for a known mass lets us determine $𝑘$ , much as temperature measurements let us find the cooling constant. You will do this in the exercises.

4.4.3Revising the model: resistance to motion

Our formula predicts that the oscillations continue forever, with the same amplitude. A spring we release on a table usually settles down. What did we leave out? We deliberately ignored resistance to motion. Let us revisit that choice.

A resisting force should point opposite the velocity: left when the mass moves right, and right when it moves left. For a simple model, suppose its magnitude is proportional to speed. With a new positive constant $𝑏$ , this gives

𝐹 d r a g = - 𝑏 𝑥' .

Notice the difference between the two forces. The spring responds to position; drag responds to velocity. They can point in the same direction or in opposite directions, depending on how the mass is moving. Adding their contributions to the total force gives

𝑚 𝑥 ″ = - 𝑘 𝑥 - 𝑏 𝑥', o r 𝑚 𝑥 ″ + 𝑏 𝑥' + 𝑘 𝑥 = 0 .

We have changed the model by adding one term. Setting $𝑏 = 0$ recovers our first equation. For $𝑏 > 0$ , we do not yet have a systematic solution method, but we can compute motions and see what the new term does.

Figure 4.17 Compare the models with and without drag, starting from the same position and velocity. Here $𝑚 = 1 k g$ , $𝑘 = 1 N m - 1$ , $𝑥 (0) = 1 m$ , and $𝑥' (0) = - 0.2 m s - 1$ . Initially, the solid curve uses $𝑏 = 0.4 N s m - 1$ ; the dashed curve uses $𝑏 = 0$ . Vary $𝑏$ from $0$ to $1$ and compare the motions. These are computed solutions for illustrative parameters.

With these positive choices of $𝑏$ , the successive swings shrink and the mass approaches its resting position. The added force has captured something our first model missed. Whether this particular resistance law describes a real spring well, and which value of $𝑏$ to use, would require measurements.

How can we compute these motions with the first-order methods we already know? Let $𝑣 = 𝑥'$ be velocity. Then the second-order equation becomes

𝑥' = 𝑣, 𝑣' = - 𝑘 𝑚 𝑥 - 𝑏 𝑚 𝑣 .

The first equation says how velocity changes position; the second says how the forces change velocity. Starting from $(𝑥 0, 𝑣 0)$ , Euler's method can advance both together, just as it advanced the coffee and room temperatures. Conversely, substituting $𝑣 = 𝑥'$ into the second equation recovers our original second-order equation.

Position and velocity together form the state of this model. We can describe its motion with one second-order equation or two first-order equations. So the two extensions we have met in this chapter are connected: learning to work with systems will also help us understand equations with higher derivatives.

4.5Appendix: Modeling Earth's Temperature

Sunlight warms Earth, while radiation carries energy away into space. Can we turn this description into a model of Earth's temperature? Let us try the same process we used for coffee and populations: choose what to hold constant, build an equation, and see what its predictions tell us.

4.5.1A planet with one temperature

Our ingredients are the surface temperature, the sunlight Earth absorbs, and the radiation it emits. Let $𝑇 (𝑡)$ be the mean surface temperature, setting aside differences between places. Average the absorbed sunlight over the whole planet and call its rate $𝑄$ . We will treat $𝑄$ as constant, ignoring daily and seasonal changes. For our first attempt, suppose all radiation emitted by the surface escapes directly to space.

How should these quantities determine $𝑇'$ ? Absorbing energy warms the planet; emitting energy cools it. Just as with growth and harvesting, we need a positive contribution and a negative one. The temperature should rise when more energy arrives than leaves, and fall when more leaves than arrives.

To make this precise, suppose a fixed amount of energy is needed to raise the temperature by one degree. Call that amount $𝐶$ , measured per square meter of surface. This constant is the model's heat capacity: a larger $𝐶$ means the same energy gain produces a smaller temperature change. Then $𝐶 𝑇'$ is the rate at which energy accumulates per square meter, so

𝐶 𝑇' = 𝑄 - r a t e o f e n e r g y e m i t t e d p e r s q u a r e m e t e r .

We need a rule for the outgoing radiation. Here physics supplies one: an ideal emitter, called a blackbody, radiates at rate $𝜎 𝑇 4$ per square meter. This is the Stefan--Boltzmann law, with $𝜎 \approx 5.67 \times 10 - 8 W m - 2 K - 4$ . It uses absolute temperature, so we must measure $𝑇$ in kelvin. Choosing to approximate the surface's radiation by this law gives

𝐶 𝑇' = 𝑄 - 𝜎 𝑇 4 .

We now have a model! The temperature varies, while $𝑄$ , $𝐶$ , and $𝜎$ are constant. The absorbed sunlight is approximately $𝑄 = 238 W m - 2$ : about 238 joules arrive each second per square meter, after reflection and averaging over the planet. This value comes from observations of Earth's energy budget. With a choice of $𝐶$ and a starting temperature, we can compute a solution.

Figure 4.18 Start at $288 K$ and watch the model cool. Try a colder initial temperature and a larger heat capacity. The default $𝐶 = 2 \times 108 J m - 2 K - 1$ is illustrative; the curves show the model's predictions, not measurements of Earth's temperature history.

At a constant temperature, $𝑇' = 0$ , so $𝑄 = 𝜎 𝑇 4$ : incoming sunlight and outgoing radiation balance. Solving gives

𝑇 = (𝑄 𝜎) 1 / 4 \approx 255 K .

Above this temperature, radiation removes more energy than sunlight supplies, so the planet cools; below it, the planet warms. That explains where the plotted solutions are heading. Changing $𝐶$ changes how quickly they get there, but not the temperature where the two rates balance.

The prediction is about $- 18 \circ C$ ( $0 \circ F$ ). Earth's average surface temperature is closer to $15 \circ C$ ( $59 \circ F$ ), or $288 K$ . Our model is too cold, and adjusting $𝐶$ cannot fix it! What did we leave out?

4.5.2Adding an atmosphere

We let radiation escape straight from the surface to space. But Earth's atmosphere absorbs some of that radiation and emits radiation of its own, including radiation back toward the surface. To include this interaction, we need to follow the atmosphere's temperature too, just as we followed both the coffee and the room.

Call the surface temperature $𝑇 𝑠 (𝑡)$ and the atmospheric temperature $𝑇 𝑎 (𝑡)$ . We will describe the atmosphere as one layer with one temperature, ignoring differences with height.

Sunlight is concentrated at shorter wavelengths than the radiation emitted by Earth's surface, and the atmosphere absorbs them differently. For a simple model, suppose sunlight passes through the atmosphere and warms the surface, while all infrared radiation from the surface is absorbed by the atmosphere. Keep the same absorbed sunlight rate $𝑄$ as before.

Absorbing radiation warms the atmosphere. The atmosphere also emits radiation of its own: downward toward the surface and upward toward space. Let us apply the same blackbody law to each side of our layer, giving an emission rate $𝜎 𝑇 4 𝑎$ in each direction. For now, these exchanges of radiation are the only ways our model gains and loses heat.

We can now write the surface's equation. It gains energy from sunlight at rate $𝑄$ and from the atmosphere at rate $𝜎 𝑇 4 𝑎$ , and loses energy by emitting radiation at rate $𝜎 𝑇 4 𝑠$ . Choosing a constant surface heat capacity $𝐶 𝑠$ gives

𝐶 𝑠 𝑇' 𝑠 = 𝑄 + 𝜎 𝑇 4 𝑎 - 𝜎 𝑇 4 𝑠 .

The atmosphere gains energy from the surface at rate $𝜎 𝑇 4 𝑠$ . Its upward and downward emissions each remove energy at rate $𝜎 𝑇 4 𝑎$ . With its own constant heat capacity $𝐶 𝑎$ , its equation is

𝐶 𝑎 𝑇' 𝑎 = 𝜎 𝑇 4 𝑠 - 2 𝜎 𝑇 4 𝑎 .

The factor of two counts the two directions. Each temperature now affects the other's rate of change, giving us a system. Choosing both initial temperatures lets us explore it numerically as we did with the coffee and the room.

Figure 4.19 Begin with the surface at $288 K$ and the atmosphere at $255 K$ . Change either initial temperature or the surface heat capacity and follow both temperatures. The illustrative defaults are $𝐶 𝑠 = 2 \times 108$ and $𝐶 𝑎 = 2 \times 107 J m - 2 K - 1$ ; these values have not been fitted to an observed warming history.

The temperatures settle at different values. Setting both derivatives to zero tells us that incoming and outgoing radiation balance for each:

𝑄 + 𝜎 𝑇 4 𝑎 = 𝜎 𝑇 4 𝑠, 𝜎 𝑇 4 𝑠 = 2 𝜎 𝑇 4 𝑎 .

Substituting the second equation into the first gives $𝜎 𝑇 4 𝑎 = 𝑄$ , and then $𝜎 𝑇 4 𝑠 = 2 𝑄$ . Thus

𝑇 𝑎 = (𝑄 𝜎) 1 / 4 \approx 255 K, 𝑇 𝑠 = (2 𝑄 𝜎) 1 / 4 \approx 303 K .

Adding the atmosphere has made the surface warmer: it now receives radiation from the atmosphere as well as sunlight. Unlike the coffee and room, these temperatures need not become equal, because energy keeps entering from the Sun and leaving for space. This simple greenhouse model shows how an atmosphere can maintain a warmer surface.

But $303 K$ is about $30 \circ C$ ( $86 \circ F$ ), so our revised prediction is too warm. One assumption to revisit is that the atmosphere absorbs all radiation from the surface. In the exercises, you will let some pass through and see how the prediction changes.

We have repeated the same process: a choice of quantities and simple rules gave us an equation; its prediction led us to revise those choices. Adding a second temperature captured an important interaction, and the new prediction suggests where to look next.

Problems

Modeling begins before an equation is solved. In every problem, name the state, say what each term means, and check signs or units before trusting a calculation. The explorations ask you to make new modeling choices; more than one answer may be defensible when the assumptions are stated clearly.

Check Your Understanding

1Read the assumptions from the equation

For each model, identify the state variables and parameters, describe the mechanism represented by every term, and give one limiting case which checks your interpretation.

(a)

$𝑁' = 0.4 𝑁 (1 - 𝑁 / 800)$ .

(b)

$𝑉' = 12 - 0.03 𝑉$ .

(c)

$𝐻' = 𝑎 𝐻 - 𝑏 𝐻 𝐿$ , $𝐿' = - 𝑐 𝐿 + 𝑑 𝐻 𝐿$ .

(d)

$𝑚 𝑥 ″ = - 𝑘 𝑥 - 𝑏 𝑥'$ .

2Predicting the coffee's temperature

A cup of coffee sits in a room held at $20 \circ C$ . Its temperature is initially $80 \circ C$ , and ten minutes later it is $50 \circ C$ . Assume that Newton's law of cooling describes its temperature throughout this interval and afterward.

(a)

Use the two measurements to determine $𝑘$ , taking time in minutes. State its units and write the resulting temperature function with no undetermined constants.

(b)

Predict the coffee's temperature twenty minutes after the first measurement.

3Exponential first, logistic later

A bacterial culture begins with mass $𝑁 (0) = 12$ mg and initially has a per-capita growth rate of $0.7 h - 1$ .

(a)

Write the exponential model and predict the mass after six hours.

(b)

Suppose the environment has carrying capacity $𝐾 = 500$ mg. Write the logistic revision and use its exact solution to predict the mass after six hours.

(c)

At what population does the logistic model's total growth rate reach its maximum? At what population has its per-capita rate fallen to half of the low-density value? Explain why these are the same population in this model.

4Harvesting a pond

For $𝑃' = 𝑃 (1 - 𝑃 / 200) - 𝐻$ , answer the following without first using the quadratic formula.

(a)

Complete the square in the natural-growth term to show that its maximum is $50$ at $𝑃 = 100$ .

(b)

For $𝐻 = 42$ , find both equilibrium populations and classify them from the sign of $𝑃'$ .

(c)

A manager wants an attracting equilibrium of at least $140$ fish. Determine the largest constant quota the model permits. Then give one reason to choose a smaller quota in practice.

5A tank as a balance law

A tank initially holds $300$ L. Water enters at $18$ L/min, and a drain removes $4 %$ of the current volume per minute.

(a)

Derive the differential equation and check the units of both terms.

(b)

Find and classify the equilibrium volume without solving the equation.

(c)

Solve the initial value problem. Use the formula to answer when the volume is within $5$ L of equilibrium.

(d)

The physical tank holds only $400$ L. Does the mathematical solution ever overflow it? Explain how this feasibility check changes if the initial volume is $500$ L.

6What changes a climate equilibrium?

The absorbed sunlight in the climate appendix can be written as $𝑄 = 𝑆 0 4 (1 - 𝛼)$ , where $𝑆 0$ is the incoming solar intensity and $𝛼$ is the fraction reflected. The factor $1 / 4$ averages sunlight intercepted by Earth's disk over its spherical surface.

Suppose absorbed sunlight increases by a constant amount $𝐹 0$ . The equation becomes $𝐶 𝑇' = 𝑆 0 4 (1 - 𝛼) - 𝜎 𝑇 4 + 𝐹 0$ .

(a)

Derive the equilibrium temperature and show directly which parameters can move it.

(b)

Without solving for $𝑇 (𝑡)$ , explain what doubling $𝐶$ changes and what it cannot change.

(c)

Compute the new equilibrium when $𝑆 0 = 1361 W m - 2$ , $𝛼 = 0.30$ , and $𝐹 0 = 4 W m - 2$ . Compare it with the unforced equilibrium. Use $𝜎 = 5.67 \times 10 - 8 W m - 2 K - 4$ .

7Two histories, one moving point

Consider the chapter's predator–prey model with $𝛾 = 1$ :

𝐻' = 𝐻 (1 - 𝐿), 𝐿' = 𝐿 (𝐻 - 1) .

(a)

Calculate $(𝐻', 𝐿')$ at each state

(2, 12), (2, 2), (12, 2), (12, 12) .

Draw the four points in the phase plane and attach an arrow showing the velocity at each. Interpret each arrow in terms of growing or declining populations.

(b)

In the chapter's figure showing time histories beside the phase-plane trajectory, locate a maximum of the hare population and the next maximum of the lynx population. Mark the corresponding points on a sketch of the trajectory. At which point is the motion vertical, and at which is it horizontal? Explain using the derivatives.

(c)

At the states $(2, 12)$ and $(2, 2)$ , the hare population is the same. Can an equation $𝐻' = 𝑓 (𝐻)$ describe the hare's rate of change at both states? Explain what information is missing when we keep track of hares alone.

8Transfers must cancel

Here we use fractions of the population rather than counts. If the constant total population is $𝑁$ , divide each count by $𝑁$ . The count-based infection term $𝑏 𝑆 𝐼$ then becomes a fraction-based term with coefficient $𝛽 = 𝑏 𝑁$ ; the recovery coefficient $𝛾$ stays the same. Below, $𝑆, 𝐼, 𝑅$ denote the fractions.

𝑆' = - 𝛽 𝑆 𝐼, 𝐼' = 𝛽 𝑆 𝐼 - 𝛾 𝐼, 𝑅' = 𝛾 𝐼 .

A quantity computed from a solution is conserved if it remains constant along that solution.

(a)

Differentiate $𝑆 + 𝐼 + 𝑅$ and explain why every transfer cancels. Deduce that $𝑆 + 𝐼 + 𝑅 = 1$ for all times of existence if the initial fractions sum to $1$ .

(b)

Explain why knowing $𝑆 (𝑡)$ and $𝐼 (𝑡)$ determines $𝑅 (𝑡)$ . Write the two equations we could solve first, and the formula recovering $𝑅$ afterward.

(c)

At the state $(𝑆, 𝐼, 𝑅) = (0.80, 0.15, 0.05)$ with $𝛽 = 0.6$ and $𝛾 = 0.2$ , compute $(𝑆', 𝐼', 𝑅')$ and interpret every sign.

(d)

Perform one simultaneous Euler step with $ℎ = 0.1$ . Verify that the three new fractions still sum to $1$ .

(e)

Add the three Euler update formulas at an arbitrary state. Show that their sum is unchanged by every simultaneous step, apart from rounding errors in a computer calculation.

(f)

A student updates $𝑆$ , then uses that new value of $𝑆$ while updating $𝐼$ . Explain why this is not the simultaneous Euler update and why exact cancellation may be lost.

9Spring data determine parameters

A $0.50$ kg mass on an undamped spring completes one oscillation every $1.2$ s.

(a)

Use $T = 2 𝜋 \sqrt 𝑚 / 𝑘$ to determine the spring stiffness $𝑘$ , including its units.

(b)

The mass begins $4$ cm to the right of equilibrium with velocity $- 0.10$ m/s. Find the constants in $𝑥 (𝑡) = 𝐴 c o s (𝜔 𝑡) + 𝐵 s i n (𝜔 𝑡)$ .

(c)

Find the first time at which the mass passes through equilibrium. Check that the sign of its velocity then agrees with the initial motion.

(d)

Let $𝑣 = 𝑥'$ . Write the spring equation as a first-order system for $(𝑥, 𝑣)$ , using the parameters you found, and give its initial state. Explain why specifying the initial position alone would not determine the motion.

Explorations

10A threshold built one factor at a time

Suppose a population has trouble finding mates when its numbers are small. We want a model in which populations below $75$ decline toward extinction, populations between $75$ and $600$ grow, and populations above $600$ decline toward $600$ .

(a)

Draw the desired phase line. Where must the rate vanish?

(b)

Start with the logistic rate $𝑟 𝑁 (1 - 𝑁 / 600)$ , where $𝑟 > 0$ . On which interval does it have the wrong sign?

(c)

Find a dimensionless factor, linear in $𝑁$ , which is negative below $75$ , zero at $75$ , and positive above $75$ . Multiply the logistic rate by your factor and check all the required signs.

(d)

Describe the futures of populations starting at $40$ , $200$ , and $800$ . Explain why they cannot cross the equilibrium populations.

(e)

If time is measured in years, what units must $𝑟$ have?

11Constant harvest or proportional harvest?

Replace the fixed quota in the trout model by removal proportional to the population: $𝑃' = 𝑃 (1 - 𝑃 / 200) - 𝑞 𝑃$ .

(a)

Find and classify the equilibria as $𝑞$ varies.

(b)

Compare the behavior near $𝑃 = 0$ with constant harvesting. Which policy prevents the model from demanding removal from an empty pond?

(c)

Compare the constant quota $𝐻 = 32$ with a proportional policy which removes the same number of fish when $𝑃 = 160$ . Choose $𝑞$ and compare their predictions when $𝑃 = 40$ . What assumption about policy does each term encode?

12The vanishing snowball

A spherical snowball melts at a rate proportional to its surface area. Its radius is $10$ cm initially and $8$ cm five minutes later. Assume its density remains constant.

(a)

Start with volume $𝑉 = 43 𝜋 𝑟 3$ and surface area $𝐴 = 4 𝜋 𝑟 2$ . Translate the melting assumption into an equation for $𝑉'$ , then use the chain rule to obtain an equation for $𝑟'$ .

(b)

Determine when the snowball disappears. Why is the radius linear in time even though both volume and area are nonlinear functions of radius?

(c)

Name one physical change near the end of melting which could make the model fail.

13An oven that is still warming

An object is placed in an oven just as the oven is switched on. Both initially have temperature $20 \circ C$ . Suppose the oven's temperature is prescribed by

𝐴 (𝑡) = 200 - 180 𝑒 - 𝑡 / 10, 𝑡 \geq 0,

where time is measured in minutes and temperature in degrees Celsius. Model the object by a single temperature $𝑇 (𝑡)$ , and assume it obeys Newton's law with $𝑘 = 1 / 20 m i n - 1$ . Assume the object does not appreciably affect the oven's temperature.

(a)

Write the initial value problem for $𝑇$ . What does the model predict for $𝑇' (0)$ ? Explain why that does not mean the object will remain at its initial temperature.

(b)

In Chapter 2 we subtracted the constant surrounding temperature. Try the corresponding substitution $𝑢 = 𝑇 - 𝐴 (𝑡)$ here. Derive

𝑢' = - 𝑘 𝑢 - 𝐴' (𝑡) .

What does the extra term represent?

(c)

Solve for $𝑇 (𝑡)$ using an integrating factor. Check both the differential equation and the initial condition.

(d)

Show that $20 < 𝑇 (𝑡) < 𝐴 (𝑡) < 200$ for every $𝑡 > 0$ , and find the limiting temperature of the object. Sketch the two temperature histories together.

(e)

The temperature difference $𝐴 (𝑡) - 𝑇 (𝑡)$ starts at zero and eventually returns toward zero. When is it largest? Use the differential equation to explain why this is also when the object warms most rapidly.

14A partially transparent greenhouse

Modify the climate appendix's greenhouse model by assuming the atmosphere absorbs a fraction $𝜀$ of the surface's infrared radiation and transmits the rest directly to space, where $0 \leq 𝜀 \leq 1$ . The layer emits $𝜀 𝜎 𝑇 4 𝑎$ upward and the same amount downward. Sunlight still passes through the atmosphere, supplying energy to the surface at rate $𝑄$ per unit area.

Let $𝐶 𝑠, 𝐶 𝑎 > 0$ be the constant heat capacities per unit area of the surface and atmosphere.

(a)

Write differential equations for $𝑇 𝑠$ and $𝑇 𝑎$ . Account for each incoming and outgoing contribution, and check that $𝜀 = 1$ recovers the chapter's equations.

(b)

Add the equations to find

𝑑 𝑑 𝑡 (𝐶 𝑠 𝑇 𝑠 + 𝐶 𝑎 𝑇 𝑎) .

Identify the two ways radiation escapes to space. Explain why energy exchanged between the surface and atmosphere disappears from this combined balance.

(c)

Suppose $0 < 𝜀 \leq 1$ . Set both derivatives equal to zero and solve for the equilibrium surface temperature in terms of $𝑄, 𝜎, 𝜀$ . Do the heat capacities affect this temperature?

(d)

Now set $𝜀 = 0$ in the original differential equations. What happens to the atmospheric temperature? Find the surface equilibrium and compare it with the limit of your formula from part (c) as $𝜀 \to 0$ . Why does the equilibrium calculation no longer determine an atmospheric temperature?

(e)

Find the value of $𝜀$ which makes the surface equilibrium $288 K$ when the bare-planet equilibrium is $255 K$ . Explain what matching this one number does—and does not—establish about the model.

15Vaccination as a transfer

Suppose susceptible individuals in an SIR model are vaccinated at a per-capita rate $𝜈$ and move directly to the removed compartment.

(a)

Modify all three equations. Make the new transfer appear once with a minus sign and once with a plus sign.

(b)

Show that $𝑆 + 𝐼 + 𝑅$ remains conserved.

(c)

At what instantaneous rate does vaccination reduce the susceptible population when $𝑆 = 0.7$ and $𝜈 = 0.03$ per day? Explain why a constant term $- 0.03$ would represent a different program.

Python: Models as Experiments

16Recover a growth law from Carlson's data

The arrays below contain the eighteen hourly yeast measurements used in the chapter.

import numpy as np
import matplotlib.pyplot as plt

t = np.arange(1, 19, dtype=float)
N = np.array([
    9.6, 18.3, 29.0, 47.2, 71.1, 119.1, 174.6, 257.3, 350.7,
    441.0, 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6
])

(a)

Compute the midpoint population and midpoint estimate of per-capita growth for every consecutive pair. Plot the estimated rate against the midpoint population. Label the axes with units.

(b)

Use np.polyfit to fit a line to this rate plot. Read estimates of $𝑟$ and $𝐾$ from the line's two intercepts.

(c)

Plot the data together with the logistic solution using your estimates. Describe one feature the curve captures and one feature the two-parameter model leaves unexplained.

17Force a climate model

Suppose the absorbed sunlight in the climate appendix changes by a prescribed amount $𝐹 (𝑡)$ , while reflection and the radiation law retain their chosen values. These are hypothetical experiments with the model; the inputs below are not measured solar histories.

Adapt your Euler function from Chapter 2 to the model $𝑇' = (𝑄 - 𝜎 𝑇 4 + 𝐹 (𝑡)) / 𝐶$ .

Figure 4.20 In this illustrative experiment, the absorbed sunlight increases by $4 W m - 2$ at year 1. The three curves use $0.5 𝐶 0$ , $𝐶 0$ , and $2 𝐶 0$ , with $𝐶 0 = 2 \times 108 J m - 2 K - 1$ . They approach the same new equilibrium, about $1.06 K$ warmer, at different speeds. The following calculations let you explore other prescribed inputs.

Use $𝑄 = 238 W m - 2$ and start from the corresponding unforced equilibrium. Measure time in years, take $𝐶 = 8$ and $40 W y r m - 2 K - 1$ , use $ℎ = 0.02$ year, and simulate fifty years. The right-hand side then returns kelvins per year.

(a)

Write three forcing functions: a step from $0$ to $4 W m - 2$ at year $1$ , a ramp increasing by $0.1 W m - 2$ per year, and a zero-mean sinusoid of amplitude $4 W m - 2$ and period ten years.

(b)

Plot the response to each forcing for two values of $𝐶$ . Keep every other parameter and the time grid fixed. Which conclusions are robust when $𝐶$ changes?

(c)

For the periodic forcing, measure the time difference between one forcing peak and the following temperature peak. How does this lag change with $𝐶$ ?

18Does one recovery rate describe the outbreak?

These are the daily illness and recovery counts from the Tristan da Cunha record used in the chapter. Time is measured in days from the first observation. As in the text, use the number currently ill as an approximation to the infectious population $𝐼$ .

I = [1,1,3,7,6,10,13,13,14,14,17,10,6,6,4,3,1,1,1,1,0]
R = [0,0,0,0,5,7,8,13,13,16,16,24,30,31,33,34,36,36,36,36,37]

(a)

Starting from $𝑅' = 𝛾 𝐼$ , show that a constant recovery rate must satisfy

𝛾 = 𝑅 (𝑏) - 𝑅 (𝑎) \int 𝑏 𝑎 𝐼 (𝑡) 𝑑 𝑡

on any interval with a nonzero denominator. Check the units.

(b)

Estimate the integral over $[0, 20]$ by adding trapezoid areas between consecutive observations. Use it to estimate $𝛾$ and compare your answer with the value obtained by fitting both curves in the chapter.

(c)

Repeat the calculation separately for $[0, 10]$ and $[10, 20]$ . If the observations followed the model exactly and we knew the integrals exactly, how would these estimates compare? How closely do your estimates agree?

(d)

Use your whole-record estimate $̂ 𝛾$ to calculate

̂ 𝑅 (𝑛) = 𝑅 (0) + ̂ 𝛾 𝑛 - 1 \sum 𝑗 = 0 𝐼 (𝑗) + 𝐼 (𝑗 + 1) 2, 𝑛 = 1, \dots, 20 .

Plot these reconstructed recovery counts alongside the observations. Why must the final counts agree? Do the intermediate counts agree?

(e)

Discuss what this comparison tells us about using one constant recovery rate. Consider the discrete observations, the integral approximation, and the assumption that being ill coincides with being infectious. Does disagreement by itself establish that the biological recovery rate changed?

19Simulate transfers and audit conservation

Write a function returning the three SIR rates and use simultaneous vector Euler updates

state = state + h * rates(state)

to simulate an epidemic from $(𝑆, 𝐼, 𝑅) = (0.99, 0.01, 0)$ . Use $𝛽 = 0.6 d a y - 1$ , $𝛾 = 0.2 d a y - 1$ , $ℎ = 0.05$ day, and run for sixty days.

(a)

Plot all three compartments against time and plot S + I + R on a second axis. Report the largest deviation of the total from $1$ .

(b)

Add the vaccination transfer from the previous exploration. Compare the peak infected fraction and the time of that peak for $𝜈 = 0, 0.01, 0.03 d a y - 1$ .

(c)

Deliberately implement the sequential-update mistake from the routine problem. Plot the total again. Explain how a conservation law became a test of your code rather than merely a fact about the equations.