Skip to main content

The Ordinary Differential Equations Project

Thomas W. Judson

Section 5.1 Linearization

In Chapter 3, we explored how to solve linear systems. However, if we are asked to solve a nonlinear system such as
\begin{align*} \frac{dx}{dt} & = x - 3y + xy^2\\ \frac{dy}{dt} & = 2x - 4y - x^2y, \end{align*}
we are faced with a much more difficult task. In general, it may not be possible to find solutions for a nonlinear system in terms of elementary functions. However, for a given modeling problem, we can ask many questions that may be answered without finding an explicit solution for the associated system of differential equations. For example, what are the equilibrium points? Are the equilibrium points stable? Do closed solution curves exist in the phase plane?

Subsection 5.1.1 Equilibrium Solutions

An equilibrium solution for a first-order system of differential equations
\begin{align*} x'(t) & = f(x, y)\\ y'(t) & = g(x, y) \end{align*}
is a point \((x_0, y_0)\) such that
\begin{equation*} f(x_0, y_0) = g(x_0, y_0) = 0. \end{equation*}
Notice that neither \(x\) nor \(y\) is changing at this point. If we have the initial conditions, \(x(0) = x_0\) and \(y(0) = y_0\text{,}\) then the solution to the system is \(x(t) = x_0\) and \(y(t) = y_0\text{.}\) Of course the interesting question is what happens if our initial conditions are close to an equilibrium solution. Do solutions tend toward the equilibrium solution, away from the equilibrium solution, or is there a combination of the two?
One of the most useful methods of determining the nature of an equilibrium solution for a given nonlinear system is to approximate the nonlinear system with a linear system. More specifically, an equilibrium solution occurs for the linear system
\begin{align*} \frac{dx}{dt} \amp = f(x,y)\\ \frac{dy}{dt} \amp = g(x, y) \end{align*}
when an \(x\)-nullcline intersects a \(y\)-nullcline. That is, when a curve defined by \(dx/dt = f(x,y) = 0\) intersects a curve defined \(dy/dt = g(x,y) = 0\text{,}\) we have an equilibrium solution. Since the \(x\) and \(y\)-nullclines are simply curves in the \(xy\)-plane, we can approximate them locally by intersecting straight lines. Translating the pair of intersection lines to the origin, we obtain a linear system, and we can apply everything that we learned about such systems in Chapter 3.

Example 5.1.1.

Consider the system
\begin{align*} \frac{dx}{dt} & = x - 3y + xy^2\\ \frac{dy}{dt} & = 2x - 4y - x^2y. \end{align*}
From
\begin{align*} \frac{dx}{dt} & = x - 3y + xy^2 = 0\\ \frac{dy}{dt} & = 2x - 4y - x^2y = 0, \end{align*}
we can quickly conclude that the only equilibrium solution to the system is \((0, 0)\text{.}\) The phase portrait for this system is given in Figure 5.1.2. If we have the initial conditions \(x(0) = 1\) and \(y(0) = 1\text{,}\) we can see that the solution tends toward the origin as \(t \to \infty\text{.}\) However, it is unclear from the phase portrait if the solution curves of all initial value problems with initial conditions near the origin tend towards the equilibrium solution as \(t \to \infty\text{.}\)
a direction field of slope arrows with one solution curve approaching the origin and a second solution curve moving away from the origin
Figure 5.1.2. Phase Portrait for \(x' = x - 3y + xy^2\) and \(2x - 4y - x^2y\)
Since the nonlinear terms of the system do not contribute much towards \(dx/dt\) and \(dy/dt\) for values of \(x\) and \(y\) near zero, we can determine the nature of the equilibrium solution by examining the system consisting of only linear terms on the right-hand side of the equation,
\begin{align*} \frac{dx}{dt} & = x - 3y\\ \frac{dy}{dt} & = 2x - 4y. \end{align*}
The matrix for this linear system,
\begin{equation*} A = \begin{pmatrix} 1 & -3 \\ 2 & -4 \end{pmatrix}, \end{equation*}
has eigenvalues \(\lambda = -1\) and \(\mu = -2\) with eigenvectors \(\mathbf u = (3,2)\) and \(\mathbf v = (1,1)\text{,}\) respectively. The general solution for this linear system is
\begin{align*} x(t) & = 3 c_1 e^{-t} + c_2 e^{-2t}\\ y(t) & = 2c_1 e^{-t} + c_2 e^{-2t}. \end{align*}
This indeed suggests that solutions near the origin tend towards the origin as \(t \to \infty\text{.}\) In this case, we say that the equilibrium solution is stable. Of course, if we are given an initial condition such as \(x(0) = -0.5\) and \(y(0) = 1\text{.}\)

Example 5.1.3. A Competing Species Model.

Suppose that \(x\) and \(y\) are the population of two distinct species that compete for the same resources. For example, two species of fish may compete for the same food in a lake or sheep and cattle competing for the same grazing land. Recall from Section 2.2 that we can model two competing species using the following system of first-order differential equations,
\begin{align*} \frac{dx}{dt} & = \alpha x \left(1 - \frac{x}{M} \right) - \beta xy\\ \frac{dy}{dt} & = \gamma y \left(1 - \frac{y}{N} \right) - \delta xy. \end{align*}
The first term in each equation is the logistic growth of each species. The second term in each equation tells how each species is affected by interacting with the competing species.
More specifically, let us consider the following system,
\begin{align} \frac{dx}{dt} & = 2x\left(1-\frac{x}{2}\right) - xy\tag{5.1.1}\\ \frac{dy}{dt} & = 3y\left(1- \frac{y}{3}\right) - 2xy.\tag{5.1.2} \end{align}
It is easy to see that the four equilibrium solutions: \((0,0)\text{,}\) \((0,3)\text{,}\) \((2,0)\text{,}\) and \((1, 1)\text{.}\) We can view the direction field, the phase plane, and some solution curves for this system in Figure 5.1.4.
a direction field of slope arrows and nullclines with two solution curves aproaching an equilibrium solution on the horizontal axis and a third solution curve approaching an equilibrium solution on the vertical axis
Figure 5.1.4. A Competing Species Model
Let us analyze what happens at the equilibrium solution \((1, 1)\text{.}\) If we decide on an appropriate change of variables, we can translate the entire system so that this equilibrium solution is at the origin. If we let
\begin{align*} u & = x - 1\\ v & = y - 1, \end{align*}
then
\begin{gather*} \frac{du}{dt} = \frac{d}{dt}(x - 1) = \frac{dx}{dt},\\ \frac{dv}{dt} = \frac{d}{dt}(y - 1) = \frac{dy}{dt}. \end{gather*}
Equations (5.1.1) and (5.1.2) now become
\begin{align} \frac{du}{dt} & = -u - v - u^2 -uv,\tag{5.1.3}\\ \frac{dv}{dt} & = -2u - v - 2uv - v^2.\tag{5.1.4} \end{align}
As before we consider only the linear part of equations (5.1.4) and (5.1.4) are
\begin{align} \frac{du}{dt} & = -u - v,\tag{5.1.5}\\ \frac{dv}{dt} & = -2u - v.\tag{5.1.6} \end{align}
The idea is that the linear part is a good local approximation to the original equations much like a tangent line is a good local approximation to a smooth function in calculus. We can determine the local nature of the equilibrium solution by examining the eigenvalues of the matrix
\begin{equation*} A = \begin{pmatrix} -1 & -1 \\ -2 & -1 \end{pmatrix}. \end{equation*}
The eigenvalues of \(A\) are \(\lambda = -1 + \sqrt{2} = 0.4142\) and \(\mu = -1 - \sqrt{2} \approx -2.4142\) with eigenvectors \(\mathbf u = (1,\sqrt{2}\,) \approx (1, 1.4142)\) and \(\mathbf v = (1,-\sqrt{2}\,) \approx (1, -1.4142)\text{,}\) respectively. The solution to the linear system is now
\begin{equation*} \mathbf x(t) = c_1 e^{\lambda t} \begin{pmatrix} 1 \\ \sqrt{2} \end{pmatrix} + c_2 e^{\mu t} \begin{pmatrix} 1 \\ -\sqrt{2} \end{pmatrix}. \end{equation*}
or
\begin{align} x(t) \amp = c_1 e^{\lambda t} + c_2 e^{\mu t}\tag{5.1.7}\\ y(t) \amp = c_1 \sqrt{2} e^{\lambda t} - c_2 \sqrt{2} e^{\mu t}.\tag{5.1.8} \end{align}
If we have initial conditions \(x(0) = 3\) and \(y(0) = 4\text{,}\) one can quickly determine that \(c_1 = 3/2 + \sqrt{2} \approx 2.9142\) and \(c_2 = 3/2 - \sqrt{2} \approx 0.0859\text{,}\) and equations (5.1.7) and (5.1.8) become
\begin{align} x(t) \amp = 2.9142 e^{0.4142 t} + 0.0859 e^{-2.4142 t}\tag{5.1.9}\\ y(t) \amp = 4.1213 e^{0.4142 t} - 0.1213 e^{-2.4142 t}.\tag{5.1.10} \end{align}
As \(t \to \infty\) notice that both \(x(t)\) and \(y(t)\) become very large and tend away from the origin.
We can conclude that the equilibrium solution \((1,1)\) is unstable. That is, tend away from the equilibrium solution as \(t \to \infty\) if one population begins with a slight advantage over the other. If neither popluation has an initial advantage over the other, then the solution curve will approach the equilibrium solution as \(t \to \infty\text{.}\)

Example 5.1.5. A Nonpolynomial Example.

If systems such as the following can be approximated by linear systems,
\begin{align*} \frac{dx}{dt} & = y\\ \frac{dy}{dt} & = -y - \sin x. \end{align*}
Certainly, this sytems has an equilibrium solution at \((0,0)\) (Figure 5.1.6). We can expand \(\sin x\) into a power series,
\begin{equation*} \frac{dy}{dt} = -y - \sin x = - y - \left( x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots \right). \end{equation*}
Thus, this system can be approximated by the linear system
\begin{align*} \frac{dx}{dt} & = y\\ \frac{dy}{dt} & = -y - x. \end{align*}
See Figure 5.1.7.
a direction field of slope arrows with one solution curve approaching the origin and a second solution curve moving away from the origin
Figure 5.1.6. Phase portrait of \(x' = y\) and \(y' = -y - \sin x\)
a direction field of slope arrows with two solution curves spiraling into the origin
Figure 5.1.7. Phase portrait of \(x' = y\) and \(y' = -y - x\)

Example 5.1.8.

For example, consider the system
\begin{align} \frac{dx}{dt} \amp = x^2 - 4x - y + 4\tag{5.1.11}\\ \frac{dy}{dt} \amp = x^3 - y.\tag{5.1.12} \end{align}
The \(x\) and \(y\)-nullclines of this system are the curves \(y = x^2 - 4x + 4\) and \(y = x^3\text{,}\) respectively. Since the two nullclines intersect only at \((1, 1)\text{,}\) we have a single equilibrium solution. From the phase plane, it appears that \((1, 1)\) is a stable equilibrium solution. That is, all solution curves starting sufficiently close to \((1, 1)\) will approach the equilibrium solution as \(t \to \infty\) (Figure 5.1.9).
a direction field of slope arrows and nullclines with a solution curve approaching an equilibrium solution
Figure 5.1.9. Phase portrait for the system (5.1.11)–(5.1.12)
Making the substitution \(u = x - 1\) and \(v = y - 1\text{,}\) simply translates the entire system to the origin, and we obtain the system
\begin{align} \frac{du}{dt} \amp = u(u - 2) - v\tag{5.1.13}\\ \frac{dv}{dt} \amp u(u^2 + 3u + 3) - v.\tag{5.1.14} \end{align}
Our new system has \(u\) and \(v\)-nullclines \(v = u(u - 2)\) and \(v = u(u^2 + 3u + 3)\text{,}\) respectively. Notice that we have simply moved the phase portrait of the original system so that our equilibrium solution is now at the origin. Furthermore, we can approximate the \(u\) and \(v\)-nullclines by their tangent lines \(v = -2u\) and \(v = 3u\text{,}\) respectively. From Figure 5.1.9, it appears that we are approximating our original system with a linear system.
a direction field of slope arrows and nullclines with a solution curve approaching the origin
Figure 5.1.10. Phase portrait for the system (5.1.13)–(5.1.14)
In general, our goal is to approximate a nonlinear system
\begin{align} x'(t) \amp = f(x,y)\tag{5.1.15}\\ y'(t) \amp = g(x, y)\tag{5.1.16} \end{align}
near each equilibrium point \((x_0, y_0)\) with a linear system. The idea is to move the equilibrium solution to the origin with a change of coordinates and then approximate the nonlinear system with a linear system. In order to move the equilibrium solution to the origin, we will introduce new variables
\begin{align*} u \amp = x - x_0\\ v \amp = y - y_0. \end{align*}
If \((x,y)\) is close to the equilibrium solution \((x_0, y_0)\text{,}\) then \((u, v)\) will be close to the origin. Under this change of coordinates,
\begin{equation*} \frac{du}{dt} = \frac{d}{dt} (x - x_0) = \frac{dx}{dt} = f(x, y) = f(x_0 + u, y_0 + v). \end{equation*}
Similarly, \(dv/dt = dy/dt = g(x_0 + u, y_0 + v)\text{.}\)
In order to approximate the nonlinear system with a linear system, we will use a Taylor series in two variables. That is, we can write
\begin{align*} f(x_0 + u, y_0 + v) \amp = f(x_0, y_0) + f_x(x_0, y_0)u + f_y(x_0, y_0)v\\ \amp + \frac{1}{2} f_{xx}(x_0, y_0) u^2 + f_{xy}(x_0, y_0) uv + \frac{1}{2} f_{yy}(x_0, y_0) v^2 + \cdots. \end{align*}
If we only use the linear terms of the Taylor series, we obtain a fairly accurate approximation of \(f\) provided we are near the equilibrium solution. Of course, \(f\) and its linear approximation may be quite different for values far away from the equilibrium solution. The linearization of the system of equations (5.1.15)–(5.1.16) now becomes
\begin{gather*} \frac{du}{dt} = f(x_0, y_0) + f_x(x_0, y_0)u + f_y(x_0, y_0)v\\ \frac{dv}{dt} = g(x_0, y_0) + g_x(x_0, y_0)u + g_y(x_0, y_0)v. \end{gather*}
Since \((x_0, y_0)\) is an equilibrium solution, the constant terms vanish in each equation and
\begin{align*} \frac{du}{dt} \amp = f_x(x_0, y_0)u + f_y(x_0, y_0)v\\ \frac{dv}{dt} \amp = g_x(x_0, y_0)u + g_y(x_0, y_0)v. \end{align*}
If we write our linearization in matrix form, we have
\begin{equation*} \begin{pmatrix} du/dt \\ dv/dt \end{pmatrix} = \begin{pmatrix} f_x(x_0, y_0) \amp f_y(x_0, y_0) \\ g_x(x_0, y_0) \amp g_y(x_0, y_0) \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix}. \end{equation*}
The matrix
\begin{equation*} J = \begin{pmatrix} f_x(x_0, y_0)\amp f_y(x_0, y_0) \\ g_x(x_0, y_0) \amp g_y(x_0, y_0) \end{pmatrix} \end{equation*}
is the Jacobian matrix of the system. We can now classify equilibrium solutions of nonlinear systems by examining the eigenvalues of the Jacobian matrix of the system or by using the trace-determinant plane. For example, if \(D =\det(J) \gt 0\) and \(T =\trace(J) \lt 0\text{,}\) we have a sink. We can determine the sink to be spiral or nodal by examining whether \((T, D)\) lies above or below the parabola \(T^2 = 4D\) in the trace-determinant plane.
It is important to note that linearization only tells us the local story. A solution curve might behave quite differently if it is far away from the equilibrium solution.

Example 5.1.11.

In Example 5.1.8, we considered the system
\begin{align*} \frac{dx}{dt} \amp = x^2 - 4x - y + 4\\ \frac{dy}{dt} \amp = x^3 - y. \end{align*}
The Jacobian matrix of this system is
\begin{equation*} J = \begin{pmatrix} f_x(x_0, y_0) \amp f_y(x_0, y_0) \\ g_x(x_0, y_0) \amp g_y(x_0, y_0) \end{pmatrix} = \begin{pmatrix} 2x_0 - 4 \amp -1 \\ 3x_0^2 \amp -1 \end{pmatrix} \end{equation*}
For the equilibrium solution \((1, 1)\text{,}\)
\begin{equation*} J = \begin{pmatrix} -2 \amp -1 \\ 3 \amp -1 \end{pmatrix}. \end{equation*}
Since \(J\) has eigenvalues \(\lambda = (-3 \pm i \sqrt{11})/2\text{,}\) the equilibrium solution will act as a spiral sink for initial values close to \((1,1)\text{.}\)

Activity 5.1.1. Classifying Equilibrium Solutions for Nonlinear Systems.

Find all of the nullclines and equilbrium solutions for each of the following systems. Classify each equilibrium solution as stable or unstable. Plot the nullclines and a direction field for each system.
 1 
You will probably find Sage very useful.
(a)
\begin{align*} \frac{dx}{dt} & = x(2 + x) + 2xy\\ \frac{dy}{dt} & = y(2 + y) + 2xy \end{align*}
(b)
\begin{align*} \frac{dx}{dt} & = x(2 + x) -2xy\\ \frac{dy}{dt} & = y(2 + y) - 2xy \end{align*}
(c)
\begin{align*} \frac{dx}{dt} & = x(2 + x) + 2xy\\ \frac{dy}{dt} & = y(2 + y) - 2xy \end{align*}
(d)
\begin{align*} \frac{dx}{dt} & = x(2 + x) + 2xy\\ \frac{dy}{dt} & = y^2(2 + y) - 2xy \end{align*}

Subsection 5.1.2 When Linearization Fails

There are at least two cases when linearization does not give us the information that we seek. First, it might well be the case that the linear terms vanish in the nonlinear system. For example, the system
\begin{align*} \frac{dx}{dt} \amp = xy\\ \frac{dy}{dt} \amp = -x^2 + xy \end{align*}
has an equilibrium solution at the origin, but this linearization of this system vanishes.
A more subtle example is the system
\begin{align*} \frac{dx}{dt} \amp = y - (x^2 + y^2)x\\ \frac{dy}{dt} \amp = -x - (x^2 + y^2)y. \end{align*}
The origin is an equilibrium solution, and the linearization of this system is
\begin{align*} \frac{dx}{dt} \amp = y\\ \frac{dy}{dt} \amp = -x. \end{align*}
The eigenvalue of the matrix corresponding to the linear system,
\begin{equation*} A = \begin{pmatrix} 0 \amp 1 \\ -1 \amp 0 \end{pmatrix}, \end{equation*}
are \(\lambda = \pm i\text{.}\) Thus, the solution curves of the linear system are simply circles about the origin. However, the nonlinear system acts quite differently. In the nonlinear system, solutions spiral out slowly from the origin (Figure 5.1.12). This system has no periodic solutions.
a direction field of slope arrows with a solution curve spiraling out from the origin towards a circular trajectory
Figure 5.1.12. When linearization fails

Subsection 5.1.3 Important Lessons

  • We can approximate a nonlinear system
    \begin{align*} x'(t) \amp = f(x,y)\\ y'(t) \amp = g(x, y) \end{align*}
    near each equilibrium point \((x_0, y_0)\) with a linear system by using a Taylor series approximation for \(f\) and \(g\text{.}\) The matrix of our linear approximation,
    \begin{equation*} J = \begin{pmatrix} f_x(x_0, y_0) \amp f_y(x_0, y_0) \\ g_x(x_0, y_0) \amp g_y(x_0, y_0) \end{pmatrix}, \end{equation*}
    is the Jacobian matrix of the system. We can classify nonlinear systems by examining the Jacobian matrix of the system and using the trace-determinant plane.
  • Linearization only tells us how solutions behave near the equilibrium point. A solution curve might behave quite differently if it is far away from the equilibrium solution.
  • In some cases, linearization can fail.

Reading Questions 5.1.4 Reading Questions

1.

What does it mean to linearize a nonlinear system?

2.

How might the linearization of a nonlinear system fail?

Exercises 5.1.5 Exercises

Linearization.

For each of the systems in Exercise Group 5.1.5.1–6,
  1. Plot and label the nullclines for equation in the system.
  2. Find all of the equilibrium solutions for the system.
  3. Use the Jacobian to classify each equilibrium solution (spiral source, nodal sink, etc.).
1.
\begin{align*} \frac{dx}{dt} & = -6x - 2x^2 -3xy\\ \frac{dy}{dt} & = y - xy - y^2 \end{align*}
2.
\begin{align*} \frac{dx}{dt} & = x(2 - x - y)\\ \frac{dy}{dt} & = 2y(2 - x - 2y) \end{align*}
3.
\begin{align*} \frac{dx}{dt} & = x(5 - x - y)\\ \frac{dy}{dt} & = y(20 - x - 2y) \end{align*}
4.
\begin{align*} \frac{dx}{dt} & = x(-x^2 - y^2 + 25)\\ \frac{dy}{dt} & = y(x + 2y -10) \end{align*}
5.
\begin{align*} \frac{dx}{dt} & = x(y^2 - x)\\ \frac{dy}{dt} & = y(y - 1) \end{align*}
6.
\begin{align*} \frac{dx}{dt} & = y\\ \frac{dy}{dt} & = - \cos x + y \end{align*}

7.

Consider the following three systems
  1. \begin{align*} x' & = 3 \sin x + y\\ y' & = 4x + \cos y - 1 \end{align*}
  2. \begin{align*} x' & = -3 \sin x + y\\ y' & = 4x + \cos y - 1 \end{align*}
  3. \begin{align*} x' & = -3 \sin x + y\\ y' & = 4x + 3\cos y - 3 \end{align*}
All three systems have an equilibrium solution at \((0, 0)\text{.}\) Which two systems have phase portraits with the same “local picture” near \((0, 0)\text{?}\) Justify your answer.

8.

Let us consider an epidemic model for a city. We make the following additional assumptions about our model.
  • The city has a constant birth rate rate of \(\alpha\) persons per day. All new born babies are susceptible to the disease.
  • Infected people either recover or die after a certain number of days. If an individual recovers, he or she is immune.
If we let \(S(t)\) be the number of susceptible individuals at time \(t\) and \(I(t)\) be the number of infected individuals at time \(t\text{,}\) our assumptions give rise to the following system of differential equations,
\begin{align*} \frac{dS}{dt} & = - \alpha SI + \beta\\ \frac{dI}{dt} & = - \gamma I + \alpha SI. \end{align*}
The constant \(\alpha\) is determined by the probability of an interaction between a susceptible individual and an infected individual, and \(\gamma\) is the rate at which infected individuals are removed from the population. If
\begin{align*} \frac{dS}{dt} & = - \alpha SI + \beta = 0\\ \frac{dI}{dt} & = - \gamma I + \alpha SI = 0, \end{align*}
then both the susceptible and infected populations do not change. This will occur at
\begin{align*} S_0 & = \frac{\gamma}{\alpha}\\ I_0 & = \frac{\beta}{\gamma}. \end{align*}
We are interested in the behavior of solutions near \((S_0, I_0)\text{.}\) If solutions approach this equilbrium point, then the disease will become endemic to the population.

9.

Consider the predator-prey system modeled by the following equations,
\begin{align*} \frac{dx}{dt} & = ax - \alpha xy = x(a - \alpha y)\\ \frac{dy}{dt} & = -by + \beta x y = y(-b + \beta x). \end{align*}
  1. Find the equilibrium solutions of this system.
  2. What are the eigenvalues of the Jacobian for each equilibrium solution?
  3. What, if anything, can be said about the nature of the equilibrium solutions of the system?
PrevTopNext