Skip to main content

The Ordinary Differential Equations Project

Thomas W. Judson

Section 2.1 Modeling with Systems

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

Subsection 2.1.1 Predator-Prey Systems

Suppose that we have a predator-prey system consisting of a population of foxes (\(F\)) and of rabbits (\(R\)).
 1 
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.
 2 
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.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.
two oscillating curves with the curve for the foxes reaching a lower peak after the higher peak for the rabbit population
Figure 2.1.1. A simple predator-prey system
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).
a closed loop where the horizontal axis is the rabbit population and the vertical axis is the fox population
Figure 2.1.2. A solution curve in the \(RF\)-plane
a closed loop on a direction field of slope arrows with the rabbit population plotted against the fox population
Figure 2.1.3. The phase plane for a predator-prey system
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 the 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.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.
two oscillating curves that approach a constant value with the higher valued curve representing the rabbit population
Figure 2.1.4. Solutions for a modified predator-prey system
three curves on a direction field of slope arrows with the curves spiraling in towards an equilibrium solution
Figure 2.1.5. Phase plane for a modified predator-prey system

Activity 2.1.1. Modifying Predator-Prey Systems.

Consider the following predator-prey system,
\begin{align*} x' & = ax(1 - x/N) - bxy\\ y' & = -cy + dxy, \end{align*}
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 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
\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*}

Activity 2.1.2. A Harmonic Oscillator.

Consider the equation
\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.

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 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:
\begin{align*} s(t) & = \text{Susceptible individuals}\\ i(t) & = \text{Infected individuals}\\ r(t) & = \text{Removed individuals} \end{align*}
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,
\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 infected populations interact in an SIR epidemic in Figure 2.1.6.
susceptible and infected populations plotted against time with the susceptible population decreasing to zero and the infected population increasing to a peak and then decreasing to zero
Figure 2.1.6. The susceptible and removed populations for an SIR epidemic
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 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{.}\)
  • 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 first-order 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
\begin{align*} \frac{dx}{dt} & = 4x + 2y\\ \frac{dy}{dt} & = 3x - y \end{align*}
  1. \(x(t) = e^{-2t} + 2e^{5t}\text{,}\) \(y(t) = -3e^{-2t} + e^{5t}\)
  2. \(x(t) = 2e^{-2t} \text{,}\) \(y(t) = -6 e^{-2t} \)
  3. \(x(t) = 2e^{-2t} \text{,}\) \(y(t) = -5 e^{-2t} \)
  4. \(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
\begin{align*} \frac{dx}{dt} & = 4x -3y\\ \frac{dy}{dt} & = 3x + 4y \end{align*}
  1. \(x(t) = 2e^{4t}\cos 3t\text{,}\) \(y(t) = 2e^{4t}\sin 3t\)
  2. \(x(t) = 2e^{4t}\cos 3t - 3e^{4t}\sin 3t\text{,}\) \(y(t) = 2e^{4t}\sin 3t + 3e^{4t}\cos 3t\)
  3. \(x(t) = 3e^{4t}\cos 3t - 4e^{4t}\sin 3t\text{,}\) \(y(t) = 3e^{4t}\sin 3t + 5e^{4t}\cos 3t\)
  4. \(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\) do not depend on \(t\text{;}\) that is,
\begin{align*} \frac{dx}{dt} & = f(x,y)\\ \frac{dy}{dt} & = g(x,y). \end{align*}
Which of the following systems are autonomous?
  1. \begin{align*} \frac{dx}{dt} & = 4x -3y\\ \frac{dy}{dt} & = 3x + 4y \end{align*}
  2. \begin{align*} \frac{dx}{dt} & = 4x -3y + \cos x\\ \frac{dy}{dt} & = 3x + 4y - \sin x \end{align*}
  3. \begin{align*} \frac{dx}{dt} & = 4x -3y +3tx\\ \frac{dy}{dt} & = 3x + 4y - e^x \end{align*}
  4. \begin{align*} \frac{dx}{dt} & = 4x \\ \frac{dy}{dt} & = 4y \end{align*}

4.

Change each of the following second-order initial value problems into a system of equations.
  1. \begin{align*} x'' + 6x' - 2x \amp = 0\\ x(0) \amp = 1\\ x'(0) \amp = 1 \end{align*}
  2. \begin{align*} y'' + y' - 2y \amp = 0\\ y(0) \amp = 2\\ y'(0) \amp = 0 \end{align*}
  3. \begin{align*} \frac{d^2x}{dt^2} + \frac{dx}{dt}' + 4x \amp = \sin t\\ x(0) \amp = 4\\ x'(0) \amp = -3 \end{align*}
  4. \begin{align*} x'' + tx' - 2 \amp = 0\\ x(0) \amp = 0\\ x'(0) \amp = 1 \end{align*}

5.

Change each of the following second-order initial value problems into a system of equations.
  1. \begin{align*} x'' + 6x' - 2x \amp = 0\\ x(0) \amp = 1\\ x'(0) \amp = 1 \end{align*}
  2. \begin{align*} y'' + y' - 2y \amp = 0\\ y(0) \amp = 2\\ y'(0) \amp = 0 \end{align*}
  3. \begin{align*} \frac{d^2x}{dt^2} + \frac{dx}{dt}' + 4x \amp = \sin t\\ x(0) \amp = 4\\ x'(0) \amp = -3 \end{align*}
  4. \begin{align*} x'' + tx' - 2 \amp = 0\\ x(0) \amp = 0\\ x'(0) \amp = 1 \end{align*}

6.

Consider the SIR model
\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.
  1. Modify the SIR model to account for the situation where the susceptible population is reproducing at a rate proportional to the current population.
  2. 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.
  3. 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.
  1. Consider the function \(y(t) = \cos \beta t\text{.}\) Under what conditions on \(\beta\) is \(y(t)\) a solution to the differential equation?
  2. 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.
  3. 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
\begin{align} x' & = -x - y\tag{2.1.1}\\ y' & = x + y/5\tag{2.1.2}\\ x(0) & = 1\tag{2.1.3}\\ y(0) & = -1\tag{2.1.4} \end{align}
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.
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.
PrevTopNext