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?

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.

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.

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

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

Example5.3A 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\label{equation-nonlinear01-competing-species-1}\tag{5.1}\\
\frac{dy}{dt} & = 3y\left(1- \frac{y}{3}\right) - 2xy.\label{equation-nonlinear01-competing-species-2}\tag{5.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.4.

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) and (5.2) now become
\begin{align}
\frac{du}{dt} & = -u - v - u^2 -uv,\label{equation-nonlinear01-competing-species-3}\tag{5.3}\\
\frac{dv}{dt} & = -2u - v - 2uv - v^2.\label{equation-nonlinear01-competing-species-4}\tag{5.4}
\end{align}
As before we consider only the linear part of equations (5.4) and (5.4) are
\begin{align}
\frac{du}{dt} & = -u - v,\label{equation-nonlinear01-competing-species-5}\tag{5.5}\\
\frac{dv}{dt} & = -2u - v.\label{equation-nonlinear01-competing-species-6}\tag{5.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}\label{equation-nonlinear01-competing-species-7}\tag{5.7}\\
y(t) \amp = c_1 \sqrt{2} e^{\lambda t} - c_2 \sqrt{2} e^{\mu t}.\label{equation-nonlinear01-competing-species-8}\tag{5.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.7) and (5.8) become
\begin{align}
x(t) \amp = 2.9142 e^{0.4142 t} + 0.0859 e^{-2.4142 t}\label{equation-nonlinear01-competing-species-9}\tag{5.9}\\
y(t) \amp = 4.1213 e^{0.4142 t} - 0.1213 e^{-2.4142 t}.\label{equation-nonlinear01-competing-species-10}\tag{5.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{.}\)

Example5.5A 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)\text{.}\) 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*}

Example5.8

For example, consider the system
\begin{align}
\frac{dx}{dt} \amp = x^2 - 4x - y + 4\label{equation-nonlinear01-example-poly-1}\tag{5.11}\\
\frac{dy}{dt} \amp = x^3 - y.\label{equation-nonlinear01-example-poly-2}\tag{5.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.9).

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\label{equation-nonlinear01-example-poly-3}\tag{5.13}\\
\frac{dv}{dt} \amp u(u^2 + 3u + 3) - v.\label{equation-nonlinear01-example-poly-4}\tag{5.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.9, it appears that we are approximating our original system with a linear system.

To determine the general case, we must approximate a nonlinear system
\begin{align}
x'(t) \amp = f(x,y)\label{equation-nonlinear01-nonlinear-system-1}\tag{5.15}\\
y'(t) \amp = g(x, y)\label{equation-nonlinear01-nonlinear-system-2}\tag{5.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.15)–(5.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.

Example5.11

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

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.12). This system has no periodic solutions.

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.

\begin{align*}
x' & = 3 \sin x + y\\
y' & = 4x + \cos y - 1
\end{align*}

\begin{align*}
x' & = -3 \sin x + y\\
y' & = 4x + \cos y - 1
\end{align*}

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

2

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.

Subsection5.1.5Project—Modeling the Nervous System¶ permalink

Hodgkin and Huxley developed a model for the firing of nerve cells in the 1950s. The model was a system of four differential equations that described the electrochemical transmission of neuronal signals along cell membranes in giant squid. The system itself was very difficult to analyze due to its highly nonlinear nature, but it did earn the two scientists a Nobel Prize for medicine in 1963 (the year after Watson and Crick won their prize).

Fitzhugh and Nagumo managed to simplify the Hodgkin-Huxley model to two equations
\begin{align*}
x' \amp = y + x - x^3/ 3 + I\\
y' \amp = -x + a - by,
\end{align*}
where \(a\) and \(b\) are constants satisfying a certain condition and \(I\) is a parameter. In these equations, \(x\) represents the excitability of the system and \(y\) represents a combination of other forces that return the system to rest. Similar models have been constructed to model the systolic and diastolic movements of the heart.

We will now consider the system
\begin{align*}
\frac{dx}{dt} \amp = x + y - x^3\\
\frac{dy}{dt} \amp = -0.5 x,
\end{align*}
which has been used to study nerve cells and is similar to the Fitzhugh-Nagumo model. The \(x\)-nullcline for this system is \(y = -x + x^3\text{.}\) The \(x\)-component of the direction field above this nullcline is positive. Below the nullcline the \(x\)-component of the direction field is negative.

The \(y\)-nullcline of the system is the \(y\)-axis. To the right of the \(y\)-nullcline, we have \(dy/dt \lt 0\text{,}\) and \(dy/dt \gt 0\) to the left of the \(y\)-nullcline. Thus, the \(x\) and \(y\)-nullclines divide the plane into four regions and all solutions must circulate in a clockwise manner about the origin (Figure 5.13).

The origin is the only equilibrium solution of the system; hence, we can examine the linearization
\begin{align*}
\frac{dx}{dt} \amp = x + y\\
\frac{dy}{dt} \amp = -0.5 x,
\end{align*}
to determine the nature of the equilibrium solution. Since the eigenvalues of this linear system are \((1 \pm i)/2\text{,}\) we can conclude that the origin is a spiral source. A solution that passes close to the origin tends to spiral outward; however, solutions that are far from the origin spiral inward. We shall later see that we must have a closed periodic solution to this system by the Poincaré-Bendixson Theorem.