Many situations are best modeled with a system of differential equations rather than a single equation. We have already derived a model that describes how a population of snowshoe hares and and how they interact with one of their primary predators, the lynx (Section 1.1). We denoted the population of hares by \(H(t)\) and the population of lynx by \(L(t)\text{,}\) where \(t\) is the time measured in years and derived the system of differential equations
\begin{align*}
\frac{dH}{dt} & = aH - bHL,\\
\frac{dL}{dt} & = -cL + dHL.
\end{align*}

Just as in first-order differential equations, we can examine the equilibrium solutions of a system. More specifically, we define an equilibrium solution for a system of differential equations
\begin{align*}
\frac{dx}{dt} & = f(x, y)\\
\frac{dy}{dt} & = g(x, y)
\end{align*}
to be those values of \(x\) and \(y\) such that \(f(x, y) = 0\) and \(g(x, y) = 0\text{.}\) That is, an equilibrium solution is a solution where neither \(x(t)\) or \(y(t)\) is changing.

Suppose that we have a predator-prey system consisting of a population of foxes (\(F\)) and of rabbits (\(R\)).^{ 11 }Foxes are omnivores. Their diet consists of small mammals, including rabbits, as well as fruits, berries, and vegetables. They will even eat fish and crabs.
\begin{align*}
\frac{dR}{dt} & = 2R - RF,\\
\frac{dF}{dt} & = -5F + RF,
\end{align*}
we have an equilibrium solution at \(R = 5\) and \(F = 2\text{.}\) That is, the system is in balance and there is just enough prey to support a constant population of predators at the point \((R, F) = (5, 2)\text{.}\)

If the number of rabbits or foxes changes, then the system is no longer in balance. For example, if \(R = 1\) and \(F = 1\text{,}\) then
\begin{align*}
\frac{dR}{dt} & = 1,\\
\frac{dF}{dt} & = -4,
\end{align*}
and the rabbit population will be increasing while the fox population will decrease. If we assume that we have initial conditions
\begin{gather*}
R_0 = R(0) = 1,\\
F_0 = F(0) = 1,
\end{gather*}
we can apply a numerical algorithm to generate a solution for our system.^{ 12 }You will find technology extremely useful when analyzing systems. We will introduce Sage commands for analyzing systems of equations at the end of this section. The graphs of the solutions for \(R(t)\) and \(F(t)\) are given in Figure 2.1. Notice that the solutions are periodic with the same period. Observe that a peak in the rabbit population is followed by a peak in the fox population.

We can graph the solution to our system in a different manner—we can construct a parametric plot of our solution in the \(RF\)-plane. Thus, a point on the graph is given by \((R(t), F(t))\) at time \(t\text{.}\) We can view the solution curve of our system in the \(RF\)-plane in Figure 2.2. The \(RF\)-plane is called the phase plane for our system of differential equations and is analogous to the phase line that we used during our investigation of slope fields for autonomous differential equations. We can plot many solutions to our predator prey-system and even plot direction fields in the phase plane (Figure 2.3).

We will now modify our system by assuming that the rabbit population will grow logistically if there are no predators present. The system can now be written as
\begin{align*}
\frac{dR}{dt} & = aR(1 - R/N) - bRF,\\
\frac{dF}{dt} & = -cF + dRF,
\end{align*}
where \(N\) is the carrying capacity. As a specific example, consider system
\begin{align*}
\frac{dR}{dt} & = 2R\left(1 - \frac{R}{10}\right) - RF,\\
\frac{dF}{dt} & = -5F + RF.
\end{align*}

It is easy to see that we have two equilibrium solutions—one at \((R, F) = (0, 0)\) and one at \((R, F) = (5, 1)\text{.}\) Our solutions now behave very differently from the assumption that the population of the prey grows exponentially. If we have initial values \(R_0 = 1\) and \(F_0 = 1\text{,}\) then our solution is no longer periodic (Figure 2.4). In fact, the solutions tend towards the equilibrium solution. The phase plane for our modified predator-prey system is given in Figure 2.5. The equilibrium solution \((5, 1)\) is an example of a stable equilibrium solution.

Subsection2.1.2The Spring-Mass Model Revisited¶ permalink

Recall the spring-mass model from Section 1.1. We have a mass lying on a flat surface that is attached to one end of a spring with the other end of the spring attached to a wall. The spring displacement is denoted by \(x\text{.}\) If \(x \gt 0\text{,}\) then the spring is stretched. If \(x \lt 0\text{,}\) the spring is compressed. If \(x = 0\text{,}\) then the spring is in a state of equilibrium (Figure 1.4). If the surface is frictionless and we pull on the mass, then the mass will oscillate.

If we use a dashpot to add damping, the second-order differential equation for the spring-mass system is
\begin{equation*}
mx'' = -bx' - kx.
\end{equation*}
or
\begin{equation*}
mx'' + bx' + kx = 0,
\end{equation*}
where \(m\text{,}\) \(b\text{,}\) and \(k\) are all positive constants.

A second-order differential equation can be restated in terms of a system of first order differential equations. Given the differential equation
\begin{equation*}
mx'' + bx' + kx = F(t),
\end{equation*}
with initial position \(x(0) = x_0\) and initial velocity \(x'(0) = x_0'\text{,}\) we can rewrite this equation as a system of first-order differential equations by letting \(v(t) = x'(t)\text{.}\) In this case, the equation becomes
\begin{equation*}
mv' + bv + kx = F(t).
\end{equation*}
We now have a system of first-order differential equations,
\begin{align*}
x' & = v,\\
v' & = \frac{1}{m} F(t) -\frac{b}{m} v -\frac{k}{m} x,
\end{align*}
with initial conditions
\begin{align*}
x(0) & = x_0\\
v(0) & = v_0'.
\end{align*}

Systems of differential equations are very useful in epidemiology. Differential equations can be used to model various epidemics, including the bubonic plague, influenza, AIDS, and the 2015 ebola outbreak in west Africa. To understand how we might model an epidemic, we will consider a very simple situation. We will assume that we have a closed population of size \(N\text{,}\) where immigration, emigration, and birth do not play an important role. We will also ignore any deaths that are not related to our disease.

We will assume that each individual in the population falls into one of the following categories:
\begin{align*}
S(t) & = \text{Susceptible individuals}\\
I(t) & = \text{Infected individuals}\\
R(t) & = \text{Removed individuals}
\end{align*}
Susceptible individuals who do not yet have the disease and can catch the disease from infected individuals. Individuals enter the removed population by either recovering from the disease or dying. If an infected individual recovers, then the individual is immune to the disease. Schematically, we can represent the effect of the disease by the diagram
\begin{equation*}
S \longrightarrow I \longrightarrow R.
\end{equation*}
Since the population is closed, we know that
\begin{equation*}
S(t) + I(t) + R(t) = N.
\end{equation*}
This model is called an SIR model.

We can model how the disease acts with the following system of equations,
\begin{align*}
\frac{dS}{dt} & = - \alpha SI\\
\frac{dI}{dt} & = \alpha SI - \beta I\\
\frac{dR}{dt} & = \beta I.
\end{align*}
We say that \(\alpha\) is the rate of infection and \(\beta\) is the rate at which the infected are removed. That is, an infected individual either dies or recovers after \(1/ \beta\) days. Since
\begin{equation*}
\frac{d}{dt}\left[ S(t) + I(t) + R(t)\right] = \frac{d}{dt} N = 0,
\end{equation*}
we need only solve the system
\begin{align*}
\frac{dS}{dt} & = - \alpha SI\\
\frac{dI}{dt} & = \alpha SI - \beta I.
\end{align*}
Letting \(\alpha = 0.005\) and \(\beta = 0.08\text{,}\) we can see how the susceptible and removed populations interact in an SIR epidemic in Figure 2.6.

There are many questions associate with epidemic models.

Will there be an epidemic?

If there is an epidemic, how many individuals will be infected?

Is there a period of time for which individuals are exposed to the disease but exhibit no symptoms and cannot infect others?

If the disease is endemic, what is the prevalence of the infection?

Can the disease be eradicated or controlled?

What is the effect of the population age structure?

We will later explore how we can adjust the \(SIR\) model to address some of these issues.

A second-order linear equation
\begin{equation*}
a(t) x'' + b(t) x' + c(t) x = g(t)
\end{equation*}
can be written as a system of first-order equations by letting \(v(t) = x'(t)\text{.}\)

An equilibrium solution for a system of differential equations
\begin{align*}
\frac{dx}{dt} & = f(x, y)\\
\frac{dy}{dt} & = g(x, y)
\end{align*}
to be those values of \(x\) and \(y\) such that \(f(x, y) = 0\) and \(g(x, y) = 0\text{.}\)

We can graph solutions to the system
\begin{align*}
\frac{dx}{dt} & = f(x, y),\\
\frac{dy}{dt} & = g(x, y),
\end{align*}
as parametric plots in the \(xy\)-plane. A point on the graph is given by \((x(t), y(t))\) at time \(t\text{.}\) The \(xy\)-plane is called the phase plane for our system of differential equations.

Solution curves in the phase plane of
\begin{align*}
\frac{dx}{dt} & = f(x, y),\\
\frac{dy}{dt} & = g(x, y),
\end{align*}
can act very differently with respect to equilibrium solutions. Consider the two predator-prey systems
\begin{align*}
\frac{dH}{dt} & = aH - bHL,\\
\frac{dL}{dt} & = -cL + dHL,
\end{align*}
and
\begin{align*}
\frac{dH}{dt} & = aH(1 - H/\alpha) - bHL,\\
\frac{dL}{dt} & = -cL + dHL.
\end{align*}
Solutions graph as closed orbits in our example of the first system. In the example of the second system, solutions approach the equilibrium solution as \(t \to \infty\text{.}\)

Consider the following predator-prey system,
\begin{align*}
x' & = 2x(1 - x/2) - 1.2xy\\
y' & = -y + 0.9xy,
\end{align*}
where \(x(t)\) is the population of the prey and \(y(t)\) is the population of the predator at time \(t\text{.}\)

How would you modify this system to include the effect of hunting of the prey at a rate of \(\alpha\) units of prey per unit of time?

How would you modify this system to include the effect of hunting of the predators at a rate proportional to the number of predators?

Suppose the predators discover a second, unlimited source of food, but they still prefer to eat the primary prey when they can catch them. How would you modify the system to include this assumption?

Suppose the predators discover a second source of food that is limited in supply. How would you modify the system to include this assumption?

Suppose that the predators migrate to an area if there is more than five times as many prey in that area (\(x > 5y\)), and that they move away if there are fewer than five times as many prey as predators. How would you modify the system to take this into account?

Suppose that prey move out of an area at a rate proportional to the number of predators in the area. How would you modify the system to take this into account?

2

Consider the equation
\begin{equation*}
\frac{d^2}{dt^2} + \frac{k}{m} y = 0
\end{equation*}
for the motion of a simple harmonic oscillator.

Consider the function \(y(t) = \cos \beta t\text{.}\) Under what conditions on \(\beta\) is \(y(t)\) a solution to the differential equation?

What initial condition (\(t = 0\)) in the \(yv\)-plane corresponds to this solution?

In terms of \(k\) and \(m\text{,}\) what is the period of this solution?

Sketch the solution curve in the \(yv\)-plane that corresponds to this solution.

3

A mass weighing 4 pounds stretches a spring 4 inches.

Consider the function \(y(t) = \cos \beta t\text{.}\) Under what conditions on \(\beta\) is \(y(t)\) a solution to the differential equation?

Formulate an initial-value problem that corresponds to the motion of this undamped mass-spring system if the mass is extended 1 foot from its rest position and released with no initial velocity.

Using the result of the previous exercise, find the solution of this initial-value problem.

4

SEIR model

Subsection2.1.6Project—Using Sage to Solve Systems¶ permalink

We can use Sage to plot the solution of the system
\begin{align}
x' & = -x - y\label{equation-systems01-desolve-rk4-system-1}\tag{2.1}\\
y' & = x + y/5\label{equation-systems01-desolve-rk4-system-2}\tag{2.2}\\
x(0) & = 1\label{equation-systems01-desolve-rk4-system-3}\tag{2.3}\\
y(0) & = -1\label{equation-systems01-desolve-rk4-system-4}\tag{2.4}
\end{align}
First, let us plot a direction field for our system.

Notice the the vectors are different lengths depending on their magnitudes. If we wish all of the vectors to be the same length, we can divide each component by the length of the vector.

We can also add axes labels to our plot.

Now suppose that we wish to plot solutions \(x(t)\) and \(y(t)\) to our system. We can use the Sage command desolve_system_rk4.

We new have a numerical approximation of the solution to the system (2.1)–(2.4). However, our approximation is just a very long list of points. In fact, we get a list of triples, \((t, x, y)\text{.}\) It would be much more useful if we could display a graph of the solution. In the code below, we grab pairs \((t,x)\) and \((t, y)\) and then plot the points using the line command.

We can also add colors, axes labels, and legend colors so that the plot makes more sense.

To plot the solution in the \(xy\)–plane, we will need to select the second two entries in each triple.

We can adjust the thickness of the plot add an arrow to indicate the direction of the solution.