are those values of \(x\) and \(y\) such that both \(f(x, y) = 0\) and \(g(x, y) = 0\text{.}\)

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

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

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.

Subsection2.1.1Predator-Prey Systems

Suppose that we have a predator-prey system consisting of a population of foxes (\(F\)) and of rabbits (\(R\)).^{ 1 }

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

we can apply a numerical algorithm to generate a solution for our system.^{ 2 } The graphs of the solutions for \(R(t)\) and \(F(t)\) are given in Figure 2.1.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.1.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.1.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

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.1.4). In fact, the solutions tend towards the equilibrium solution. The phase plane for our modified predator-prey system is given in Figure 2.1.5. The equilibrium solution \((5, 1)\) is an example of a stable equilibrium solution.

where \(x(t)\) is the population of the prey and \(y(t)\) is the population of the predator at time \(t\text{.}\)

(a)

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?

(b)

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

(c)

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?

(d)

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

(e)

Suppose that the predators migrate to a different area if there are more than three times as many prey as predators in that area (\(x > 3y\)) and they move back if there are fewer than three times as many prey as predators. How would you modify the system to take this into account?

(f)

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?

Subsection2.1.2The Spring-Mass Model Revisited

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

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*}
\frac{d^2 y}{dt^2} + \frac{k}{m} y = 0
\end{equation*}

for the motion of a simple harmonic oscillator.

(a)

Rewrite the second-order differential equation, \(d^2 y/dt^2 + (k/m) y = 0\text{,}\) as a system of two first-order differentail equations, where \(v = y'\text{.}\)

(b)

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

(c)

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

(d)

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

(e)

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

Subsection2.1.3Modeling Epidemics

Systems of differential equations are very useful in epidemiology. Differential equations can be used to model various epidemics, including the bubonic plague, influenza, AIDS, the 2015 ebola outbreak in west Africa, and the COVID-19 pandemic. 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:

Susceptible individuals are those 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,

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{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 infected populations interact in an SIR epidemic in Figure 2.1.6.

There are many questions associated 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?

To explore \(SIR\) models in more depth, see Project 2.5.1.

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.

\begin{align*}
\frac{dS}{dt} & = - \alpha SI\\
\frac{dI}{dt} & = \alpha SI - \beta I,
\end{align*}

where \(S\) is the susceptible population and \(I\) is the infected population.

Modify the SIR model to account for the situation where the susceptible population is reproducing at a rate proportional to the current population.

Modify the SIR model to account for the situation where the susceptible population is decreasing at a constant rate such as susceptibles leaving an infected area or city.

Modify the SIR model to account for the situation where the infected population is increasing at a constant rate such as infected individuals entering an infected area or city from the outside.

7.

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.

Subsection2.1.7Using Sage to Solve Systems

We can use Sage to plot the solution of the system

First, let us plot a direction field for our system.

Notice the vectors have different lengths depending on their magnitudes. If we wish all of the vectors to have 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 now have a numerical approximation of the solution to the system (2.1.1)–(2.1.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.