If we consider the phase portrait of a system
\begin{align*}
\frac{dx}{dt} & = f(x, y)\\
\frac{dy}{dt} & = g(x, y),
\end{align*}
we may see closed orbits and various types of equilibrium solutions. If we make a slight change in the system, then we might reasonably expect the phase portrait to change slightly. For example, a closed orbit might expand or contract, an equilibrium solution might shift, or a spiral source might tighten up, but we would not expect new equilibrium solutions to appear or a spiral source to turn into a spiral sink. However, this may be exactly what happens. If such a change occurs in the system through an adjustment of a parameter, we say a bifurcation has occurred.

The simplest bifurcations occur when an equilibrium solution appears, disappears, or splits into two or more equilibrium solutions. For example, the system
\begin{align*}
\frac{dx}{dt} & = \alpha + x^2\\
\frac{dy}{dt} & = -y
\end{align*}
has no equilibrium solution if \(\alpha \gt 0\) (Figure 5.33). However, we have a single equilibrium solution at the origin if \(\alpha = 0\) and two equilibrium solutions, \((\pm \sqrt{-\alpha}, 0)\text{,}\) if \(\alpha \lt 0\text{.}\)

The Jacobian of the system
\begin{align*}
\frac{dx}{dt} & = \alpha + x^2\\
\frac{dy}{dt} & = -y
\end{align*}
is
\begin{equation*}
J(x,y)
=
\begin{pmatrix}
2x & 0 \\
0 & -1
\end{pmatrix}.
\end{equation*}
If \(\alpha = 0\text{,}\) then \(J\) has eigenvalues 0 and \(-1\text{.}\) This is not an elementary equilibrium point such a saddle, a sink, or a spiral source. We call this type of equilibrium solution a saddle-node. On the right half of the \(xy\)-plane, the equilibrium solution resembles a saddle, while on the left half it resembles a sink (Figure 5.34).

If \(\alpha \lt 0\text{,}\) then
\begin{equation*}
J
=
\begin{pmatrix}
2 \sqrt{\alpha} & 0 \\
0 & -1
\end{pmatrix}
\text{ or }
\begin{pmatrix}
-2 \sqrt{\alpha} & 0 \\
0 & -1
\end{pmatrix}.
\end{equation*}
In the first case, we have a saddle since the eigenvalues are real and of opposite sign. In the second case, we have a nodal sink, since both eigenvalues are negative. The two equilibrium solutions move in opposite directions as \(\alpha\) decreases. We can summarize what happens for various values of \(\alpha\) with a bifurcation diagram such as the one in Figure 5.35.

Now let us consider an entirely different type of bifurcation by examining the system
\begin{align*}
\frac{dx}{dt} & = \alpha x + 5y - x(x^2 + y^2)\\
\frac{dy}{dt} & = -5x + \alpha y - y(x^2 + y^2).
\end{align*}
The origin is an equilibrium solution for all values of \(\alpha\text{.}\) The linearization of our system is
\begin{align*}
\frac{dx}{dt} & = \alpha x + 5y\\
\frac{dy}{dt} & = -5x + \alpha y.
\end{align*}
Since the eigenvalues of this system are \(\alpha \pm 5i\text{,}\) we can easily determine the nature of of our linearization. If \(\alpha \lt 0\text{,}\) we have a spiral sink. As \(\alpha\) moves from negative values to positive values, the origin changes to a center (\(\alpha = 0\)), and then to a spiral source (\(\alpha \gt 0\)).

In the case of the nonlinear system, the origin is still a spiral sink for \(\alpha \gt 0\) (Figure 5.36). If \(\alpha \gt 0\text{,}\) the origin is still a spiral source (Figure 5.37). However, something quite different happens at \(\alpha = 0\text{.}\) The origin is stable equilibrium solution with solutions spiraling into the origin very slowly. As \(\alpha\) increases past zero, the equilibrium solution destabilizes and becomes a source. In addition, a closed orbit of radius \(\sqrt{\alpha}\) develops. Solutions inside of this closed orbit spiral out towards the orbit, while solutions outside of the orbit spiral inward (Figure 5.38).

To see exactly what happens as \(\alpha\) passes zero and becomes positive, we will rewrite our system in polar coordinates. If we make the substitution \(x = r \cos \theta\) and \(y = r \sin \theta\text{,}\) our nonlinear system can be rewritten as
\begin{align*}
\frac{dr}{dt} & = \alpha r - r^3\\
\frac{d\theta}{dt} & = -5.
\end{align*}
If \(\alpha \lt 0\text{,}\) the origin is a sink since \(\alpha r - r^3 \lt 0\) for all \(r \gt 0\text{.}\) In this case all solutions tend towards the origin as \(t \to \infty\text{.}\) When \(\alpha \gt 0\text{,}\) the origin is still an equilibrium solution. Moreover, \(r' = 0\text{,}\) when \(r = \sqrt{\alpha}\text{.}\) We also know that \(r' \gt 0\) for \(0 \lt r \lt \sqrt{\alpha}\) and \(r' \lt 0\) if \(r \gt \sqrt{\alpha}\text{.}\) So the circle of radius \(\sqrt{\alpha}\) is a periodic solution with the trajectory moving clockwise since \(\theta' = -5 \lt 0\text{.}\) All nonzero solutions spiral towards this orbit as \(t \to \infty\text{.}\) This type of bifurcation is called a Hopf bifurcation. No new equilibrium solutions arise, but a periodic solution develops as \(\alpha\) passes through the bifurcation value.

Theorem5.39Hopf Bifurcation

Suppose that \((x_0(\alpha), y_0(\alpha))\) is a equilibrium solution for the family of systems
\begin{align*}
\frac{dx}{dt} & = f_\alpha(x, y)\\
\frac{dy}{dt} & = g_\alpha(x, y)
\end{align*}
parameterized by \(\alpha\text{,}\) and the Jacobian matrix for the system evaluated at this equilibrium has eigenvalues \(a(\alpha) \pm i b(\alpha)\text{.}\) Assume that at some \(\alpha = \alpha_0\text{,}\) we have \(a(\alpha_0) = 0\text{,}\) \(a'(\alpha_0) \gt 0\text{,}\) \(b(\alpha_0) \neq 0\text{,}\) and that the equilibrium point is asymptotically stable. Then there exists an \(\alpha_1\) such that the system has a periodic solution encircling the equilibrium solution for \(\alpha_0 \lt \alpha \lt \alpha_1\text{.}\)^{ 26 }For a proof and description of the Hopf Bifurcation Theorem see C. Chicone. Ordinary Differential Equations with Applications. Springer-Verlag, New York, 1999.

Recall that our example system
\begin{align*}
\frac{dx}{dt} & = \alpha x + 5y - x(x^2 + y^2)\\
\frac{dy}{dt} & = -5x + \alpha y - y(x^2 + y^2)
\end{align*}
has an equilibrium at \((0, 0)\) with eigenvalues \(\alpha \pm 5i\text{.}\) If we let \(\alpha_0 = 0\text{,}\) then \(a(0) = 0\text{,}\) \(a'(0) = 1\text{,}\) and \(b(0) = 5\text{.}\) Thus, the hypothesis of the Hopf Bifurcation Theorem are satisfied and we are guaranteed a period solution surrounding our equilibrium solution.

Example5.40Van der Pol's equation

Van der Pol's equation, a simple nonlinear equation having applications in electrical engineering and mathematical biology, is \(x'' - x(1 - x^2)x' + x = 0\text{.}\) This equation can be written as the system
\begin{align}
\frac{dx}{dt} \amp = y\label{equation-nonlinear04-van-der-pol-a}\tag{5.30}\\
\frac{dy}{dt} \amp = -x + x(1 - x^2)y.\label{equation-nonlinear04-van-der-pol-b}\tag{5.31}
\end{align}
The phase portrait for Van der Pol's equation is given in Figure 5.41. The origin is the only equilibrium solution to Van der Pol's equation, and one might guess that this solution acts like a spiral source by examining the phase portrait. If we examine the system consisting of only linear terms on the right-hand side of the equation,
\begin{align*}
\frac{dx}{dt} \amp = y\\
\frac{dy}{dt} \amp = - x + y
\end{align*}
we might get a better sense of what is happening. The matrix for this linear system,
\begin{equation*}
A =
\begin{pmatrix}
0 \amp 1 \\
-1 \amp 1
\end{pmatrix},
\end{equation*}
has eigenvalues \(\lambda = (1 \pm i\sqrt{3})/2\text{.}\) This suggests that the origin looks like a spiral source at least locally. Indeed, the Jacobian matrix for (5.30)–(5.30) is
\begin{equation*}
J
=
\begin{pmatrix}
0 \amp 1 \\
-1 - 2 x_0 y_0 \amp 1 - x_0^2
\end{pmatrix},
\end{equation*}
which agrees with \(A\) for \((x_0, y_0) = (0,0)\text{.}\)

Let us examine how a bifurcation might occur in a predator-prey model. In our model, we will assume that a predator's appetite is satiated when food is abundant. If this is the case, an increase in the prey population has little effect on the interaction terms in our model. We might model this satiable predator-prey scenario by the equations
\begin{align*}
\frac{dx}{dt} & = (a -bx)x - \frac{cxy}{d + \alpha x}\\
\frac{dy}{dt} & = -ey + \frac{fxy}{d + \alpha x},
\end{align*}
where \(x(t)\) is the prey population at time \(t\) and \(y(t)\) is the predator population at time \(t\text{.}\) Here, \(a, \ldots, f\) and \(\alpha\) are constants. Let us choose the constants \(a, \ldots, f\) so that our system becomes
\begin{align*}
\frac{dx}{dt} & = (1 - x)x - \frac{xy}{0.3 + \alpha x}\\
\frac{dy}{dt} & = -0.5y + \frac{xy}{0.3 + \alpha x}.
\end{align*}
The constant \(\alpha\) is know as the satiation constant. The larger the value of \(\alpha\text{,}\) the more quickly a predator's appetite is satiated as the prey population increases.

Let us see examine the phase portrait for several values of \(\alpha\text{.}\) If \(\alpha = 0\text{,}\) then there is no satiation effect. In this case, all trajectories inside the population quadrant spiral asymptotically towards a single stable equilibrium point (Figure 5.42). If \(\alpha = 1.35\text{,}\) we have a very similar phase portrait with solutions spiraling in towards a stable equilibrium solution (Figure 5.43). However, if \(\alpha = 0.9\text{,}\) we have a periodic trajectory. All other trajectories are pulled toward this attracting periodic orbit (Figure 5.44).

As \(\alpha\) continues to increase a Hopf bifurcation occurs at \(\alpha \approx 0.5\text{.}\) At this point the equilibrium solution destabilizes and spawns an attracting periodic orbit. As \(\alpha\) continues to increase, the amplitude of this periodic orbit increases. However, at \(\alpha \approx 0.85\text{,}\) the amplitude of the periodic solution begins to decrease. At \(\alpha \approx 1.2\text{,}\) the equilibrium point re-stabilizes and the periodic solution is absorbed. For \(\alpha > 1.2\text{,}\) we only have a stable equilibrium solution in the first quadrant. The \(x\)-coordinate of our equilibrium solution is
\begin{equation*}
x = \frac{0.15}{1 - 0.5 \alpha},
\end{equation*}
and we can summarize our findings with a bifurcation diagram (Figure 5.45).

If we make a slight change in the system, then we might reasonably expect the phase portrait to change slightly. For example, a closed orbit might expand or contract, an equilibrium solution might shift, or a spiral source might tighten up, but we would not expect new equilibrium solutions to appear or a spiral source to turn into a spiral sink. However, this may be exactly what happens. If such a change occurs in the system through an adjustment of a parameter, we say a bifurcation has occurred.

A Hopf bifurcation if an attracting periodic solution encircling an equilibrium solution develops as \(\alpha\) passes through the bifurcation value.

Suppose that \((x_0(\alpha), y_0(\alpha))\) is a equilibrium solution for the family of systems
\begin{align*}
\frac{dx}{dt} & = f_\alpha(x, y)\\
\frac{dy}{dt} & = g_\alpha(x, y)
\end{align*}
parameterized by \(\alpha\text{,}\) and the Jacobian matrix for the system evaluated at this equilibrium has eigenvalues \(a(\alpha) \pm i b(\alpha)\text{.}\) Assume that at some \(\alpha = \alpha_0\text{,}\) we have \(a(\alpha_0) = 0\text{,}\) \(a'(\alpha_0) \gt 0\text{,}\) \(b(\alpha_0) \neq 0\text{,}\) and that the equilibrium point is asymptotically stable. Then there exists an \(\alpha_1\) such that the system has a periodic solution encircling the equilibrium solution for \(\alpha_0 \lt \alpha \lt \alpha_1\text{.}\)

A satiable predator-prey scenario can be modeled by the equations
\begin{align*}
\frac{dx}{dt} & = (a -bx)x - \frac{cxy}{d + \alpha y}\\
\frac{dy}{dt} & = -ey + \frac{fxy}{d + \alpha y},
\end{align*}
where \(x(t)\) is the prey population at time \(t\) and \(y(t)\) is the predator population at time \(t\text{.}\) Here, \(a, \ldots, f\) and \(\alpha\) are constants. The constant \(\alpha\) is know as the satiation constant. The larger the value of \(\alpha\text{,}\) the more quickly a predator's appetite is satiated as the prey population increases.

Consider the system
\begin{align*}
x' & = y\\
y' & = -x - \frac{y}{4} + x^2.
\end{align*}

Verify that the function
\begin{equation*}
L(x, y) = \frac{x^2}{2} + \frac{y^2}{2} - \frac{x^3}{3}
\end{equation*}
is a Lyapunov function for the system.

Sketch the level sets of \(L\text{.}\)

What can you conclude about the phase portrait of the system from (1) and (2)?

2

Let \(G(x, y) = x^3 - 3xy^2\text{.}\)

What is the gradient system with vector field given by the gradient of \(G\text{?}\)

Sketch the graph of \(G\) and the level sets of \(G\text{.}\)

Sketch the phase portrait of the gradient system in (1).

3

Suppose that the smell of a bunch of dead fish in the region \(-2 \leq x \leq 2\text{,}\) \(-2 \leq y \leq 2\) is given by the function
\begin{equation*}
S(x, y) = x^2 + y^2 - \frac{x^4 + y^4}{4} - 3x^2 y^2 + 100.
\end{equation*}

What is the gradient system who vector field is the gradient of \(S\text{?}\)

Use Sage to graph the phase portrait of the system.

How many dead fish are there and where are they?

Using the results from part (2), sketch the level sets of \(S\text{.}\)

Why is this model not realistic for large values of \(x\) or \(y\text{?}\)

4

Consider the system of differential equations
\begin{align*}
x' & = x(-x -y + 1)\\
y' & = y(-ax -y + b),
\end{align*}
where \(a\) and \(b\) are parameters with \(a, b \gt 0\text{.}\) Suppose that this system is only defined for \(x, y \geq 0\text{.}\)

Use nullclines to sketch the phase portrait for this system for several values of \(a\) and \(b\text{.}\)

Determine the values of \(a\) and \(b\) at which a bifurcation occurs.

5

Show that the quadratic polynomial
\begin{equation*}
V(x, y) = ax^2 + 2bxy + cy^2
\end{equation*}
is positive definite with minimum at \((x, y) = (0, 0)\) if and only if \(a \gt 0\) and \(ac - b^2 \gt 0\text{.}\)

6

Use the positive definite function \(V(x, y) = x^2 + y^2\) to argue that the system
\begin{align*}
x' & = -x - 2y^2\\
y' & = 2xy - y^3
\end{align*}
has an asymptotically stable equilibrium solution at \((x, y) = (0, 0)\text{.}\)

In 1963, an MIT professor, Edward N. Lorenz, published a paper on his research in meteorology. Using differential equations, Lorenz had developed a simplified system to model certain weather-related phenomena. When he analyzed the system, however, he found that the trajectories of the solutions were incredibly convoluted and effectively unpredictable for certain parameters. The Lorenz equations can be written as
\begin{align*}
\frac{dx}{dt} & = -\sigma x + \sigma y\\
\frac{dy}{dt} & = \rho x - y - xz\\
\frac{dz}{dt} & = -\beta z + xy,
\end{align*}
where \(\sigma\text{,}\) \(\rho\text{,}\) and \(\beta\) are constants. The Existence and Uniqueness Theorem for systems of differential equations guarantees a unique solution for each set of initial conditions,
\begin{align*}
x(0) & = x_0\\
y(0) & = y_0\\
z(0) & = z_0.
\end{align*}
However, Lorenz discovered that the trajectories of the solutions were incredibly convoluted and effectively unpredictable for certain parameters. For certain values of \(\sigma\text{,}\) \(\rho\text{,}\) and \(\beta\text{,}\) the trajectories are extremely sensitive to initial conditions. Since real data always has some inherent uncertainty, initial values are never precisely known, and we may have little success modeling real world phenomena. In addition, solutions can stay in a bounded region of the three dimensional version of the phase plane and wind through the region along an incredibly convoluted path. There is much more freedom to move around in three dimensions than there is in two.

Lorenz noticed that the system behaved strangely, when he let \(\sigma =10\text{,}\) \(\rho = 28\text{,}\) and \(\beta = 8/3\text{.}\) Thus, our system
\begin{align*}
\frac{dx}{dt} & = -10 x + 10 y\\
\frac{dy}{dt} & = 28 x - y - xz\\
\frac{dz}{dt} & = -\frac{8}{3} z + xy,
\end{align*}
defines a vector field in \({\mathbb R}^3\text{,}\)
\begin{equation*}
{\mathbf F}(x, y, z)
=
\left(
-10 x + 10 y,
28 x - y - xz,
-\frac{8}{3} z + xy
\right),
\end{equation*}
and the equilibrium solutions occur exactly when this vector field is zero. That is, \((x,y,z)\) is an equilibrium solution if
\begin{equation*}
{\mathbf F}(x, y, z)
=
\left(
-10 x + 10 y,
28 x - y - xz,
-\frac{8}{3} z + xy
\right)
=
(0, 0, 0).
\end{equation*}
It is easy to check that \((0, 0, 0)\text{,}\) \((6\sqrt{2}, 6\sqrt{2}, 27)\text{,}\) and \((-6\sqrt{2}, -6\sqrt{2}, 27)\) are the equilibrium solutions to the Lorenz system for the parameters that we have chosen. With the exception of these equilibrium solutions it is probably impossible to find solutions for the system, so we will have to use a numerical algorithm to approximate solutions (Figure 5.46).

SubsubsectionThe Equilibrium Solutions

If we wish to understand the nature of the equilibrium solutions of the Lorenz system, it makes sense to linearize the system. The Jacobian matrix of the Lorenz system is
\begin{equation*}
J(x, y, z)
=
\begin{pmatrix}
-10 & 10 & 0 \\
28 - z & -1 & -x \\
y & x & -8/3
\end{pmatrix}.
\end{equation*}
At the origin, the Jacobian is
\begin{equation*}
J
=
\begin{pmatrix}
-10 & 10 & 0 \\
28 & -1 & 0 \\
0 & 0 & -8/3
\end{pmatrix}.
\end{equation*}
Thus, the linear system that approximates the Lorenz system at the origin is
\begin{align*}
\frac{dx}{dt} & = - 10x + 10 y\\
\frac{dy}{dt} & = 28x - y\\
\frac{dz}{dt} & = -\frac{8}{3} z.
\end{align*}
This system decouples into two systems,
\begin{align*}
\frac{dx}{dt} & = - 10x + 10 y\\
\frac{dy}{dt} & = 28x - y,
\end{align*}
and
\begin{equation*}
\frac{dz}{dt} = -\frac{8}{3} z.
\end{equation*}
Consequently, the origin is a sink in the \(z\)-direction and a saddle in the \(xy\)-plane.

Now let us consider what happens at the other two equilibrium solutions. For the equilibrium point \((6\sqrt{2}, 6\sqrt{2}, 27)\text{,}\) we have the Jacobian matrix
\begin{equation*}
J
=
\begin{pmatrix}
-10 & 10 & 0 \\
1 & -1 & -6 \sqrt{2} \\
6 \sqrt{2} & 6\sqrt{2} & - 8/3
\end{pmatrix}.
\end{equation*}
The eigenvalues of this matrix are
\begin{align*}
\lambda_1 & \approx -13.8\\
\lambda_2 & \approx 0.094 + 10.2i\\
\lambda_3 & \approx 0.094 - 10.2i.
\end{align*}
Therefore, this equilibrium point is a spiral saddle that has a line of solutions that tend toward the equilibrium point as \(t\) increases and a plane of solutions that spiral towards the equilibrium point as \(t\) decreases. That is, solutions approach this equilibrium point along the direction of the eigenvalue of the negative eigenvector, then spiral slowly away along the plane corresponding to the complex eigenvalues. Similarly, the equilibrium solution \((-6\sqrt{2}, -6\sqrt{2}, 27)\) is a spiral saddle.

Now that we have some idea of how the solutions of the Lorenz system behave near the equilibrium solutions, we can start to understand global behavior of the system. A trajectory approaches one of the two equilibrium solutions \((\pm 6\sqrt{2}, \pm 6\sqrt{2}, 27)\) along the straight line corresponding to the negative eigenvalue. When the trajectory gets close enough to the equilibrium solutions, it begins to spiral away. When the spiral is large enough, the trajectory come under the influence of the saddle at the origin and either returns to repeat it previous pattern or moves over to the other equilibrium solution (Figure 5.47).

Consider the Lorenz system
\begin{align*}
x' & = -a x + ay\\
y' & = rx -y -xz\\
z' & = -bz + xy.
\end{align*}
In this exercise, we shall show that the solutions of the system are bounded.

If
\begin{equation*}
V(x, y, z) = rx^2 +ay^2 + a(z - 2r)^2,
\end{equation*}
Show that
\begin{equation*}
\dot{V} = -2a [rx^2 + y^2 + b(z - r)^2 - br^2].
\end{equation*}

It is true that \((z - r)^2 \geq (z - 2r)^2 - r^2\text{.}\) Assuming this, show that
\begin{equation*}
rx^2 + y^2 + b(z - r)^2 \geq mV - br^2,
\end{equation*}
where \(m\) is the smallest of the three numbers 1, \(1/a\text{,}\) and \(b/2a\text{.}\)

Use parts (1) and (2) to show that
\begin{equation*}
\dot{V} \leq -2maV + 4abr^2.
\end{equation*}

Use part (3) to show that \(\dot{V} \leq -2abr^2 \lt 0\) everywhere outside of the ellipsoid
\begin{equation*}
R = \left\{ (x, y, z) | V(x, y, z) \leq 3br^2/m \}\right\}.
\end{equation*}

Use part (4) to show that the ellipsoid is invariant, and that every solution curve ends up inside of \(R\text{.}\)