Section 2.1 Modeling with Systems
Objectives

To understand and be able to create models using system of differential equations
\begin{align*} \frac{dx}{dt} & = f(x, y)\\ \frac{dy}{dt} & = g(x, y). \end{align*} 
To understand that 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.

To understand that a secondorder linear equation
\begin{equation*} a(t) x'' + b(t) x' + c(t) x = g(t) \end{equation*}can be written as a system of firstorder equations by letting \(v(t) = x'(t)\text{.}\)

To understand that equilibrium solutions for a system of differential equations
\begin{align*} \frac{dx}{dt} & = f(x, y)\\ \frac{dy}{dt} & = g(x, y) \end{align*}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 firstorder 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.
Subsection 2.1.1 PredatorPrey Systems
Suppose that we have a predatorprey 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
and the rabbit population will be increasing while the fox population will decrease. If we assume that we have initial conditions
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 predatorprey 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
where \(N\) is the carrying capacity. As a specific example, consider the system
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 predatorprey system is given in FigureÂ 2.1.5. The equilibrium solution \((5, 1)\) is an example of a stable equilibrium solution.
Activity 2.1.1. Modifying PredatorPrey Systems.
Consider the following predatorprey system,
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?
Subsection 2.1.2 The SpringMass Model Revisited
Recall the springmass 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 secondorder differential equation for the springmass system is
or
where \(m\text{,}\) \(b\text{,}\) and \(k\) are all positive constants.
A secondorder differential equation can be restated in terms of a system of first order differential equations. Given the differential 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 firstorder differential equations by letting \(v(t) = x'(t)\text{.}\) In this case, the equation becomes
We now have a system of firstorder differential equations,
with initial conditions
Activity 2.1.2. A Harmonic Oscillator.
Consider the equation
for the motion of a simple harmonic oscillator.
(a)
Rewrite the secondorder differential equation, \(d^2 y/dt^2 + (k/m) y = 0\text{,}\) as a system of two firstorder 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.
Subsection 2.1.3 Modeling 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 COVID19 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
Since the population is closed, we know that
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
we need only solve the system
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.
Subsection 2.1.4 Important Lessons

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

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{.}\)

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 depending on how close they are to a particular equilibrium solution.
Reading Questions 2.1.5 Reading Questions
1.
Given a system of firstorder differential equations, explain in your own words what it means to be a stable equilibrium solution for the system.
Exercises 2.1.6 Exercises
1.
Verify whether or not the given pair of functions, \((x(t), y(t))\) forms a solution to the system
\(x(t) = e^{2t} + 2e^{5t}\text{,}\) \(y(t) = 3e^{2t} + e^{5t}\)
\(x(t) = 2e^{2t} \text{,}\) \(y(t) = 6 e^{2t} \)
\(x(t) = 2e^{2t} \text{,}\) \(y(t) = 5 e^{2t} \)
\(x(t) = c_1e^{2t} + 2c_2e^{5t}\text{,}\) \(y(t) = 3c_1 e^{2t} + c_2 e^{5t}\)
2.
Verify whether or not the given pair of functions, \((x(t), y(t))\) forms a solution to the system
\(x(t) = 2e^{4t}\cos 3t\text{,}\) \(y(t) = 2e^{4t}\sin 3t\)
\(x(t) = 2e^{4t}\cos 3t  3e^{4t}\sin 3t\text{,}\) \(y(t) = 2e^{4t}\sin 3t + 3e^{4t}\cos 3t\)
\(x(t) = 3e^{4t}\cos 3t  4e^{4t}\sin 3t\text{,}\) \(y(t) = 3e^{4t}\sin 3t + 5e^{4t}\cos 3t\)
\(x(t) =c_1e^{4t}\cos 3t  c_2e^{4t}\sin 3t\text{,}\) \(y(t) = c_1e^{4t}\sin 3t + c_2e^{4t}\cos 3t\)
3.
A system is autonomous if \(dx/dt\) and \(dx/dt\) dop not depend on \(t\text{;}\) that is,
Which of the following systems are autonomous?
 \begin{align*} \frac{dx}{dt} & = 4x 3y\\ \frac{dy}{dt} & = 3x + 4y \end{align*}
 \begin{align*} \frac{dx}{dt} & = 4x 3y + \cos x\\ \frac{dy}{dt} & = 3x + 4y  \sin x \end{align*}
 \begin{align*} \frac{dx}{dt} & = 4x 3y +3tx\\ \frac{dy}{dt} & = 3x + 4y  e^x \end{align*}
 \begin{align*} \frac{dx}{dt} & = 4x \\ \frac{dy}{dt} & = 4y \end{align*}
4.
Change each of the following secondorder initial value problems into a system of equations.
 \begin{align*} x'' + 6x'  2x \amp = 0\\ x(0) \amp = 1\\ x'(0) \amp = 1 \end{align*}
 \begin{align*} y'' + y'  2y \amp = 0\\ y(0) \amp = 2\\ y'(0) \amp = 0 \end{align*}
 \begin{align*} \frac{d^2x}{dt^2} + \frac{dx}{dt}' + 4x \amp = \sin t\\ x(0) \amp = 4\\ x'(0) \amp = 3 \end{align*}
 \begin{align*} x'' + tx'  2 \amp = 0\\ x(0) \amp = 0\\ x'(0) \amp = 1 \end{align*}
5.
Change each of the following secondorder initial value problems into a system of equations.
 \begin{align*} x'' + 6x'  2x \amp = 0\\ x(0) \amp = 1\\ x'(0) \amp = 1 \end{align*}
 \begin{align*} y'' + y'  2y \amp = 0\\ y(0) \amp = 2\\ y'(0) \amp = 0 \end{align*}
 \begin{align*} \frac{d^2x}{dt^2} + \frac{dx}{dt}' + 4x \amp = \sin t\\ x(0) \amp = 4\\ x'(0) \amp = 3 \end{align*}
 \begin{align*} x'' + tx'  2 \amp = 0\\ x(0) \amp = 0\\ x'(0) \amp = 1 \end{align*}
6.
Consider the SIR model
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 massspring 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.
Subsection 2.1.7 Using 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 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 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.
Now let us add a vector field to our plot.
To learn more about how to use Sage to solve systems of differential equations, see www.sagemath.org/doc/reference/calculus/sage/calculus/desolvers.html
. The Sage cell below can be used to make your own computations.