Section 3.1 “Solving” Differential Equations Using Graphs
Our goal in this section is to learn to use Python to create graphs similar to what we saw in
Figure 2.12. Such a graph starts with a system of differential equations and draws one possible
graphical solution to that system. The graphical solution drawn depends on the equations themselves, the value(s) of any parameter(s), and the starting values for each population in the system, that is, the values when
\(t=0\) of
\(S(t)\text{,}\) \(I(t)\text{,}\) and any other populations being modeled. For now, we specify that this type of solution is a
graphical solution, but as the book continues, we will often refer to a graphical solution as simply a
solution.
One thing we will
almost never do is solve an entire system of differential equations in the way you may have seen in a calculus course, where we wrote the solution as a formula. To remind you how this often looks in calculus class—that is, how it looks to solve by writing the solution of a differential equation as a formula—we show an example in
Activity 3.2. Then
Activity 3.3 begins the discussion about why we use a different approach in this text.
Activity 3.2.
As we mentioned back in
Chapter 2, a differential equation is just an equation in which at least one term is a derivative. This means that an equation like
\(dy/dt=\cos t + 3\) is a differential equation.
Solving a differential equation often means finding a way to write an equivalent relationship, without including a derivative. Using this interpretation of “solving”, solve
\begin{equation*}
\frac{dy}{dt}=\cos t + 3\text{.}
\end{equation*}
Answer.
The general solution is
\(y=\sin t + 3t + C\text{,}\) where
\(C\) can be any constant. In this case, we solved by computing the antiderivatives of both sides of the differential equation.
Activity 3.3.
Now consider the SI model:
\begin{align*}
\frac{dS}{dt} \amp = -\beta S I \\
\frac{dI}{dt} \amp = \beta S I \text{.}
\end{align*}
What are the challenges to solving this system of equations, where “solving” means determining formulas that do not involve derivatives for
\(S(t)\) and
\(I(t)\text{?}\)
Hint.
There are many ways to answer this. One consideration is that each of the differential equations involves both
\(S\) and
\(I\text{.}\) This means the differential equations do not stand alone: each depends on the other. This interrelationship means we cannot just compute two separate antiderivative formulas.
Another consideration is that the right-hand side of the differential equation in
Activity 3.2 involves only
\(t\text{,}\) that is, only the independent variable. The dependent variable
\(y\) does not appear except in the derivative term
\(dy/dt\text{.}\) However, in the SI model system of equations, the dependent variables
\(S\) and
\(I\) appear in the derivatives
and elsewhere in the equations. Most antiderivative methods learned in a typical introductory Calculus course do not apply in this situation.
You may have thought of yet other challenges as well.
As discussed in
Activity 3.3, solving a system of differential equations is challenging. As we include more than two compartments, finding solution formulas for
\(S(t)\text{,}\) \(I(t)\text{,}\) and additional populations becomes yet more challenging, and is often impossible. Even when we can determine solution formulas, they are typically quite complicated to work with.
Instead, our approach is to compute
numerical solutions, which are solutions in the form of numbers. The main idea is that Python starts at the initial population values, then uses the system of differential equations to “step through” the full range of time, computing population values step-by-step at later times by using information about how the population values change throughout. The system of differential equations is made up of derivatives, which are rates of change, so they exactly describe how the populations change at all times
\(t\text{.}\) The population values computed in the numerical solution can be printed, or more often, graphed, to see the results in the form of a graphical solution.
We noted earlier that graphical solutions may be called simply
solutions. Similarly, numerical solutions are sometimes called simply
solutions. And we have seen that formula-based solutions can be called simply
solutions. It is important to use context to determine which kind of solution is being discussed, in this text or in other reading.
Here is an outline of the numerical solution process.
-
Tell Python the formulas for the system of differential equations.
-
Tell Python which numbers to use for any parameters in the model. We select the numbers for a particular numerical solution, and we can use different numbers to create a different numerical solution.
-
Tell Python the
initial conditions for each population, that is, the size of each population at the starting time of the model (which is usually
\(t=0\)).
-
Tell Python when to start and stop the model, that is, the starting and ending values for time
\(t\text{.}\) Also tell Python what time steps to use as Python “steps through” from the starting to the ending time.
-
Ask Python to use all this information to numerically solve the model and plot graphical solution curves.
Activity 3.4.
Read the Python code in the block below, and click the
Evaluate (Sage) button to view the resulting graph. Think carefully about steps 1-5 in the outline above. Which lines of code address which steps of the outline? (Notice that you may need to scroll down within the code block. Alternatively, click the rectangular icon in the top right of the code block. This opens a full-screen version of the code. To return to the regular screen, click again the rectangle in the top right.)
Answer.
-
The formulas for the system of differential equations appear most directly in line 8, in order. The first formula is
-beta * S * I, and the second formula is
beta * S * I. Collectively, these are named
dy/dt. More generally, the lines 6–9 define the system of differential equations.
-
In line 11, we give
\(\beta\) the value of 0.02. This number can be changed, and the resulting numerical solution then changes accordingly.
-
In line 12, the initial conditions are set to
\(S(0)=49\) and
\(I(0)=1\text{.}\) These initial conditions are collectively named
y0. The initial conditions appear in square brackets and are separated by a comma. This syntax should look familiar from line 8, when we entered the right-hand sides of the differential equations all in one command, using the same syntax.
-
The time range for this numerical solution appears in line 13, using the
arange command from the
numpy package downloaded in line 2. The code says
t = np.arange(0.0, 8.0, 0.1). This means that time starts at
\(t=0.0\) and continues till the step just before
\(t=8.0\text{.}\) Therefore
\(t=7.9\) is the last time step for which values of
\(S(t)\) and
\(I(t)\) are computed. In addition, the
0.1 tells Python to use a step size of 0.1 in the “step through” process for computing
\(S(t)\) and
\(I(t)\) values.
-
We ask in line 14 for Python to compute a solution, called
sol. To compute this solution, Python uses the command
odeint and the inputs
hz (the differential equations in the model),
y0 (the initial conditions for this solution), and
t_range (the range of time across which to compute the solution). These commands produce a collection of numbers that describe the values of
\(S(t)\) and
\(I(t)\) at all the times specified in
t_range. These numbers are arranged in two columns: one column for all the
\(S(t)\) values, and one column for all the
\(I(t)\) values.
The commands in lines 16-18 work together to create the plot we see after clicking
Evaluate (Sage). Line 16 plots the blue
\(S(t)\) curve, using the time inputs
t_range and the values of
\(S(t)\) throughout those times. The values of
\(S(t)\) are described as
sol[:, 0]. The
: in
sol[:, 0] indicates that we want all rows of the solution, which means the solution at all times
\(t\text{.}\) The
0 in
sol[:, 0] means that we want the
\(S(t)\) values for each of those times. Python typically starts counting with 0 rather than 1, so the column of
\(S(t)\) values is numbered as column
0 rather than column
1.
Line 17 plots the orange
\(I(t)\) curve, using the time inputs
t_range and the values of
\(I(t)\) throughout those times, which are described as
sol[:, 1]. This command is very similar to
sol[:, 0], but whereas the
\(S(t)\) values were in column
\(0\) of
sol, the
\(I(t)\) values are in column
1 of
sol.
Here is a summary of the Python code above, with some of the commands explained in more detail.
In the lines
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
we import online packages that make it possible for Python to do some of the things we are about to do. Saying we are importing “as” gives us a name, such as plt or np, by which to call each package and use the tools within that package.
The lines
def hz(y, t):
S, I = y
dydt = [-beta * S * I, beta * S * I]
return dydt
define the system of differential equations, using the command def. The system is named hz, for “human-zombie”, and the (y, t) means that the model depends on variables \(y\) and \(t\text{.}\) The line S, I = y defines \(y\) as comprising both \(S\) and \(I\text{,}\) in that order. This means that \(y\) collectively represents all our dependent variables. Taking the derivative of \(y\text{,}\) then, means taking each of the derivatives \(dS/dt\) and \(dI/dt\text{.}\) The formulas for these derivatives, which are on the right-hand sides of the SI model, appear in order in the line dydt = [-beta * S * I, beta * S * I]. Finally, the line return dydt ends the definition and says that whenever we later use the definition, we will be using the derivatives defined in dydt.
A quick way to think about these lines is that they define our system of differential equations and give them the name
hz.
The lines
beta = 0.02
y0 = [49, 1]
t_range = np.arange(0.0, 8.0, 0.1)
provide values for \(\beta\text{,}\) the initial conditions \(S(0)\) and \(I(0)\text{,}\) and \(t\text{.}\) Any of these values can be changed, and the code can be run again, to produce a different numerical simulation of the same system of differential equations. Notice that the initial conditions are listed in the same order as how we defined the differential equations: \(S(0)=49\) is first, and \(I(0)=1\) is second.
The solution of the system of differential equations
sol = odeint(hz, y0, t_range)
is a numerical solution, and it includes multiple parts. The solution is named sol. The Python command odeint computes the solution. For odeint to run, it needs three things: (1) a system of differential equations hz, (2) initial conditions y0 for the differential equations, and (3) a time range t_range across which to create the solution. Though odeint computes the numerical solution, it does not show us the graphical solution: the plots we create later show us these results. Instead, the numerical solution created by odeint consists of lists of numbers for \(S\) and for \(I\text{,}\) at the times specified in t_range. (If you would like to see the numerical values computed by sol, include the line print(sol) in the block of Python code above, sometime after the line sol = odeint(hz, y0, t_range). When you then evaluate the code, the list of values will appear.)
The lines
plt.plot(t_range, sol[:, 0], color='blue')
plt.plot(t_range, sol[:, 1], color='orange')
create, but do not show, the plots of graphical solutions to the system of differential equations. For each plot, we provide the input values t_range and the output values from the numerical solution sol. As described above, sol is a set of numbers which are the values of the populations \(S(t)\) and \(I(t)\) at all the times in t_range. The numbers in sol are structured into two lists, list 0 and list 1. (Python often counts things, such as lists, starting with 0 rather than starting with 1.) The syntax sol[:,0] describes all the \(S(t)\) data. The colon : means “all” the elements, and the 0 means list 0, which is the list containing values of \(S(t)\text{.}\) The syntax sol[:,1] describes all the \(I(t)\) data, with : again representing “all” and 1 representing list 1, which contains \(I(t)\) values.
The lines above include colors for each graph,
'blue' for the
\(S(t)\) graph and
'orange' for the
\(I(t)\) graph. Stating colors is optional.
Finally,
plt.show()
shows the two plots, on a shared set of axes.
To build further understanding of how the Python code works, it helps to make several small changes to the code. Simultaneously, we continue to build intuition for how the output graph is affected by our selected parameter values, initial conditions, time range, and more. Try this in
Activity 3.5.
Activity 3.5.
Start with the same Python code as in
Activity 3.4 and try various small changes to the code. Identify where and how in the code to make the changes suggested below, and discuss the effects of each change on the resulting graph. Here and elsewhere, whenever we work with Python, feel free to try your own adaptations to the code.
-
Change the value of
\(\beta\text{.}\) Identify where in the Python code to make this change, and comment on how the output graph changes when
\(\beta\) increases or decreases. (While you can change
\(\beta\) to any value you like, recall that a realistic outbreak model requires
\(\beta > 0\text{.}\))
-
Change the initial conditions so that
\(S(0) \neq 49\text{,}\) \(I(0) \neq 1\text{,}\) or both. Comment on how these changes affect the resulting graph. (To stay in the realm of biologically feasible outbreaks, keep
\(S(0) > 0\) and
\(I(0) > 0\text{.}\))
-
Change the time range. Continue to begin at
\(t=0\text{,}\) but change the end time to something different than
\(t=8\text{.}\) Discuss reasons why a modeler may want to change the time range, and identify ways in which changing the time range causes the output graph to change.
-
Change the step size for time
\(t\) so that Python steps forward in increments of 1.0 instead of 0.1. Identify where in the Python code to make this change, and observe any ways in which the output graph looks different after the change.
These are all the basics for numerically solving a system of differential equations. In
Section 3.2 we create a graph in which data points and solution curves appear all on the same set of axes. Combining these plots provides a visual way to change the value of
\(\beta\) to improve the model’s fit to data points. Then, also in
Section 3.2, we discuss helpful ways to visualize information on graphs, so readers can best understand the results that emerge from models.
