ch_11_Systems_of_Nonlinear_Differential_Equations.qmd

---
jupyter: python3
---

# Systems of Nonlinear Differential Equations  {#sec-11}
## Autonomous Systems  {#sec-11-1}

A system of first-order differential equations is called **autonomous** when the system can be written in the form

$$
    \begin{aligned}
    \frac{dx_1}{dt} &=g_1(x_1,x_2,\cdots,x_n) \\ 
    \frac{dx_2}{dt} &=g_2(x_1,x_2,\cdots,x_n) \\ 
        & \;\vdots \\ 
    \frac{dx_n}{dt} &=g_n(x_1,x_2,\cdots,x_n)
    \end{aligned}
$$

  Notice that the independent variable $t$ does not appear explicitly on the right-hand side of each differential equation

* **Second-Order DE as a System**

  Any second-order differential equation $x''=g(x,x')$ can be written as an autonomous system. 
  If we let $y=x'$, the second-order differential equation becomes the system of two first-order equations
  
    $$
    \begin{aligned}
        x' &= y\\ 
        y' &= g(x,y) 
    \end{aligned}
    $$

**Example** $~$ The displacement angle $\theta$ for a pendulum satisfies the nonlinear second-order
  differential equation
  
  $$\frac{d^2 \theta}{dt^2} +\frac{g}{l}\sin\theta=0$$
  
  If we let $x=\theta$ and $y=\dot{\theta}$, $~$ this second-order differential equation can be written as the autonomous system
  
  $$\dot{x} = y,\;\; \dot{y} = -\frac{g}{l}\sin x$$

* **Plane Autonomous System**
  
  When $n=2$, the system is called a <font color="blue">**plane autonomous system**</font>, and we write the system as
  
  $$
  \begin{aligned}
     \frac{dx}{dt} &= P(x,y)\\ 
     \frac{dy}{dt} &= Q(x,y)
  \end{aligned}
  $$
  
  If $P$, $Q$, and the first-order partial derivatives $\partial P/\partial x$,
  $\partial P/\partial y$, $\partial Q/\partial x$, and $\partial Q/\partial y$ are continuous in a
  region $R$ of the plane, then a solution to the plane autonomous system that satisfies 
  $\mathbf{x}(0)=\mathbf{x}_0$ is unique and one of three basic types:

  * A **constant solution**, $\mathbf{x}(t)=\mathbf{x}_0$ for all $t$. A constant solution is called a **critical** or **stationary point**
    
    $$\begin{aligned}
       P(x,y) &= 0 \\ 
       Q(x,y) &= 0
     \end{aligned}$$
    
     Note that since $\mathbf{x}'=\mathbf{0}$, a critical point is a solution of the system of algebraic equations
    
  * An **arc**, $\mathbf{x}(t)$ - a plane curve that does not cross itself
    
    ![](figures/ch11_figure01.png){width="75%" fig-align="center"}
    
  * A **periodic solution** or **cycle**, $\mathbf{x}(t +p)=\mathbf{x}(t)$
    
    ![](figures/ch11_figure02.png){width="25%" fig-align="center"}
    
 **Example** $~$Find all critical points of the following plane autonomous system

$$
\begin{aligned}
 x'&= x^2 +y^2 -6\\ 
 y'&= x^2 -y 
\end{aligned}$$

## Stability of Linear Systems {#sec-11-2}

If $\mathbf{x}_1$ is a critical point of a plane autonomous system and $\mathbf{x}=\mathbf{x}(t)$ is a solution
satisfying $\mathbf{x}(0)=\mathbf{x}_0$, when $\mathbf{x}_0$ is placed near $\mathbf{x}_1$

![](figures/ch11_figure03.png){width="85%" fig-align="center"}

(@) It may return to the critical point 
(@) It may remain close to the critical point without returning
(@) It may move away from the critical point

* **Stability Analysis**

  * A careful geometric analysis of the solutions to the *linear* plane autonomous system
  
    $$
    \begin{aligned}
       x'&= ax +by\\ 
       y'&= cx +dy
    \end{aligned}
     $$
  
    in terms of the eigenvalues and eigenvectors of the coefficient matrix
  
    $$\mathbf{A}=
    \begin{pmatrix}
     a & b\\ 
     c & d
    \end{pmatrix}
     $$
   
    drives the stabilty analysis

  * To ensure that $\mathbf{x}_0=(0,\,0)$ is the only critical point, we will assume that the determinant $\Delta = ad -bc \neq 0$. If $\tau = a + d$ is the trace of matrix $\mathbf{A}$,
  then the characteristic equation, $\mathrm{det}(\mathbf{A} -\lambda\mathbf{I})=0,$ can be rewritten as
  
    $$\lambda^2 -\tau\lambda +\Delta =0$$
  
    Therefore the eigenvalues of $\mathbf{A}$ are 
  
    $$\lambda =\frac{\tau \pm \sqrt{\tau^2 -4\Delta}}{2}$$
  
    and the usual three cases for these roots occur according to whether $\tau^2 -4\Delta$ is positive, negative, or zero

**Example** $~$ Find the eigenvalues of the linear system

$$
    \begin{aligned}
       x'&= -x +y\\ 
       y'&= cx -y
    \end{aligned}
$$
 
in terms of $c$, and use a numerical solver to discover the shapes of solutions corresponding to the cases
$c=\frac{1}{4}$, $4$, $0$, and $-9$. $\mathbf{A}$ has trace $\tau=-2$ and determinant $\Delta=1 -c$, 
and so the eigenvalues are

$$ \lambda =-1 \pm \sqrt{c} $$
 
The nature of the eigenvalues is therefore determined by the value of $c$

```{python}
import numpy as np
import sympy as sp
sp.init_printing(use_unicode=True)

c = sp.Symbol('c')
A = sp.Matrix([[-1, 1], [c, -1]])
A.eigenvects()
```

```{python}
import matplotlib.pyplot as plt

w = 1
xp = np.linspace(-w, w, 50)
yp = np.linspace(-w, w, 50)
x, y = np.meshgrid(xp, yp)

c_ = np.array([1/4, 4, 0, -9])
c_title = [r'$c=\frac{1}{4}$', r'$c=4$', r'$c=0$', r'$c=-9$']

fig = plt.figure(figsize=(6, 6))

for i in range(4):
    
    ax = fig.add_subplot(2, 2, i +1)

    xdot = -x +y
    ydot = c_[i]*x -y
    
    if c_[i] >= 0.0:
        y_1 = -np.sqrt(c_[i]) *xp
        y_2 = np.sqrt(c_[i]) *xp
        ax.plot(xp, y_1, 'r:', xp, y_2, 'g:')    
    ax.streamplot(x, y, xdot, ydot, color='blue')

    ax.set_title(f'{c_title[i]}')

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', 
        labelsize=8)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 5))
    ax.yaxis.set_ticks(np.linspace(-w, w, 5))
    
    if i >= 2: 
        ax.set_xlabel(r'$x$')
    if i == 0 or i == 2: 
        ax.set_ylabel(r'$y$')
```

* <font color="red">Trajectory behaviors in phase portraits can be explained with eigenvalue-eigenvector of $\mathbf{A}$</font>

* **Real Distinct Eigenvalues**, $\tau^2 -4\Delta > 0$

  $$
    \begin{aligned}
        \mathbf{x}(t) &= c_1\mathbf{k}_1 e^{\lambda_1 t} +c_2\mathbf{k}_2 e^{\lambda_2 t}\\ 
        &\;\big\Downarrow \;{\scriptstyle \lambda_1 >\lambda_2}\\ 
        &= e^{\lambda_1 t} \left[c_1\mathbf{k}_1 +c_2\mathbf{k}_2 e^{(\lambda_2 -\lambda_1)t} \right ] \\
        &\;\big\Downarrow \;\,{\scriptstyle t \to \infty}\\
        &\simeq c_1\mathbf{k}_1 e^{\lambda_1 t}
    \end{aligned}
  $$

    * **Both eigenvalues negative**, $\tau^2 -4\Delta > 0$, $\tau<0$, $\Delta>0$
    
        **Stable Node:** Since both eigenvalues are negative, it follows that 
        $\lim_{t \to \infty} \mathbf{x}(t)=\mathbf{0}$ in the direction of $\mathbf{k}_1$ when $c_1 \neq 0$ or
        in the direction of $\mathbf{k}_2$ when $c_1=0$
        
    * **Both eigenvalues positive**, $\tau^2 -4\Delta > 0$, $\tau>0$, $\Delta>0$
    
        **Unstable Node:** $\mathbf{x}(t)$ becomes unbounded in the direction of 
        $\mathbf{k}_1$ when $c_1 \neq 0$ or in the direction of $\mathbf{k}_2$ when $c_1=0$ 
    
    * **Eigenvalues have opposite signs**, $\tau^2 -4\Delta > 0$, $\Delta<0$
    
        **Saddle Point:** When $c_1=0$, $\mathbf{x}(t)$ will approach $\mathbf{0}$ along the line determined by 
        $\mathbf{k}_2$. If $\mathbf{x}(0)$ does not lie on the line determined by $\mathbf{k}_2$, the direction
        determined by $\mathbf{k}_1$ serves as an asymtote for $\mathbf{x}(t)$ 

**Example** $~$ Classify the critical point $(0,0)$ of each of the following linear system 
$\mathbf{x}'=\mathbf{A}\mathbf{x}$ as either a stable node, an unstable node, or a saddle point

<center> 
 $(\text{a})$ $\begin{pmatrix} -2 & -2\\ -2 & -5 \end{pmatrix}$, $\;$
 $(\text{b})$ $\begin{pmatrix} -1 & -2\\ \;\;3 & \;\;4 \end{pmatrix}$, $\;$
 $(\text{c})$ $\begin{pmatrix} 2 & -1\\ 3 & -2 \end{pmatrix}$
 </center>

```{python}
import numpy as np

w = 1
xp = np.linspace(-w, w, 6)
yp = np.linspace(-w, w, 6)
x, y = np.meshgrid(xp, yp)

A = np.array([[[-2, -2], [-2, -5]], 
              [[-1, -2], [3, 4]], 
              [[2, -1], [3, -2]]])
A_title = ['Stable Node', 'Unstable Node', 'Saddle Point']

fig = plt.figure(figsize=(4, 12))

for i in range(3):
    xdot = A[i,0,0]*x +A[i,0,1]*y
    ydot = A[i,1,0]*x +A[i,1,1]*y

    lamda, v = np.linalg.eig(A[i])

    if lamda[0] >= lamda[1]:
        y_1 = v[1,0]/v[0,0]*xp
        y_2 = v[1,1]/v[0,1]*xp
    else:
        y_1 = v[1,1]/v[0,1]*xp
        y_2 = v[1,0]/v[0,0]*xp       

    ax = fig.add_subplot(3, 1, i +1)
    
    ax.plot(xp, y_1, 'r:', xp, y_2, 'g:')    
    ax.streamplot(x, y, xdot, ydot, color='blue')

    ax.set_title(A_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 5))
    ax.set_ylabel(r'$y$')
    
    ax.yaxis.set_ticks(np.linspace(-w, w, 5))
    if i == 2:
        ax.set_xlabel(r'$x$')
```

* **A Repeated Real Eigenvalue**, $\tau^2 -4\Delta = 0$

  The general solution takes on one of two different forms depending on whether **one** or **two** linearly independent eigenvectors can be found for the repeated eigenvalues $\lambda_1$

  * **Two linearly independent eigenvectors**

    If $\mathbf{k}_1$ and $\mathbf{k}_2$ are two linearly independent eigenvectors corresponding to $\lambda_1$,
  then the general solution is given by
  
    $$\mathbf{x}(t)=c_1\mathbf{k}_1 e^{\lambda_1 t} +c_2\mathbf{k}_2 e^{\lambda_1 t}=
    \left(c_1\mathbf{k}_1 + c_2\mathbf{k}_2 \right) e^{\lambda_1 t}$$
   
    If $\lambda_1<0$, the $\mathbf{x}(t)$ approaches $\mathbf{0}$ along the line determined by the vector
    $c_1\mathbf{k}_1 + c_2\mathbf{k}_2$ and the critical point is a **degenerate stable node**. The arrows are reversed when $\lambda_1>0$, and the critical point is a **degenerate unstable node**

  * **A single linearly independent eigenvectors**
  
    When only a single linearly independent eigenvector $\mathbf{k}_{11}$ exists, the general solution is   given by
  
    $$\begin{aligned}
     \mathbf{x}(t)&=c_1\mathbf{k}_{11} e^{\lambda_1 t} +c_2\left(\mathbf{k}_{11} te^{\lambda_1 t} 
            +\mathbf{k}_{12} e^{\lambda_1 t}\right)\\
         &=te^{\lambda_1 t}\left[c_2 \mathbf{k}_{11} +\frac{1}{t} \left(c_1\mathbf{k}_{11} 
            +c_2\mathbf{k}_{12}\right) \right]
    \end{aligned}$$
    
    where $(\mathbf{A} -\lambda_1\mathbf{I})\mathbf{k}_{12}=\mathbf{k}_{11}$. If $\lambda_1<0$, then $\lim_{t \to \infty} te^{\lambda_1 t}=0$ and it follows that $\mathbf{x}(t)$ approaches $\mathbf{0}$ in the line determined by $\mathbf{k}_{11}$. The critical point is again a **degenerate stable node**. When $\lambda_1>0$, $\mathbf{x}(t)$ becomes unbounded as $t$ increases, and the critical point is a **degenerate unstable node**

**Example** $~$ Classify the critical point $(0,0)$ of each of the following linear system $\mathbf{x}'=\mathbf{A}\mathbf{x}$ 

<center>
$(\text{a})$ $\begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}$, $\;$
$(\text{b})$ $\begin{pmatrix} 3 & -18\\ 2 &\; -9 \end{pmatrix}$, $\;$
$(\text{c})$ $\begin{pmatrix} \;\;2 & 4\\ -1 & 6 \end{pmatrix}$
</center>

```{python}
w = 1
xp = np.linspace(-w, w, 6)
yp = np.linspace(-w, w, 6)
x, y = np.meshgrid(xp, yp)

A = np.array([[[1, 0], [0, 1]], 
              [[3, -18], [2, -9]], 
              [[2, 4], [-1, 6]]])
A_title = ['Degenerate Unstable', 
           'Degenerate Stable', 
           'Degenerate Unstable']

fig = plt.figure(figsize=(4, 12))

for i in range(3):
   
    xdot = A[i,0,0]*x +A[i,0,1]*y
    ydot = A[i,1,0]*x +A[i,1,1]*y

    lamda, v = np.linalg.eig(A[i])

    ax = fig.add_subplot(3, 1, i +1)

    if i != 0:
        y_1 = v[1,0] /v[0,0] *xp
        ax.plot(xp, y_1, 'r:')
            
    ax.streamplot(x, y, xdot, ydot, color='blue')

    ax.set_title(A_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 5))
    ax.set_ylabel(r'$y$')
    
    ax.yaxis.set_ticks(np.linspace(-w, w, 5))
    if i == 2: ax.set_xlabel(r'$x$')           
```

* **Complex Eigenvalues**, $\tau^2 -4\Delta < 0$

  If $\lambda_1=\alpha +i\beta$ and $\bar{\lambda}_1$ are the complex eigenvalues and $\mathbf{k}_1=\mathbf{b}_1 +i\mathbf{b}_2$ is a complex eigenvector corresponding to $\lambda_1$, the general solution can be written as $\mathbf{x}=c_1\mathbf{x}_1 +c_2\mathbf{x}_2$

  $$
  \begin{aligned}
    \mathbf{x}_1(t) &=e^{\alpha t}\left(\mathbf{b}_1\cos\beta t -\mathbf{b}_2\sin\beta t\right) \\ 
   \mathbf{x}_2(t) &=e^{\alpha t}\left(\mathbf{b}_2\cos\beta t +\mathbf{b}_1\sin\beta t\right)
   \end{aligned}$$
 
  * **Pure imaginary roots**, $\tau^2 -4\Delta < 0$, $~\tau=0$

    **Center:** When $\alpha=0$, all solutions are ellipses with center at the origin and are periodic   with period $p=2\pi/\beta$. The critical point is called a **center**
  
  * **Nonezero real part**, $\tau^2 -4\Delta < 0$, $~\tau\neq 0$

    **Spiral Point:** When $\alpha<0$, $e^{\alpha t}\to 0$, and the elliptical-like solution spirals closer
  and closer to the origin. The critical point is called a **stable spiral point**. When $\alpha>0$, the 
  effect is the opposite. An elliptical-like solution is driven farther and farther from the origin, and the critical point is now called an **unstable spiral point**

**Example** $~$ Classify the critical point $(0,0)$ of each of the following linear system 
$\mathbf{x}'=\mathbf{A}\mathbf{x}$ 

<center>
$(\text{a})$ $\begin{pmatrix} -1 & 2\\ -1 & 1 \end{pmatrix}$,$\;$
$(\text{b})$ $\begin{pmatrix} -1 & -4\\ \;\;1 & -1 \end{pmatrix}$
</center>

```{python}
w = 1
xp = np.linspace(-w, w, 6)
yp = np.linspace(-w, w, 6)
x, y = np.meshgrid(xp, yp)

A = np.array([[[-1, 2], [-1, 1]],
              [[-1, -4], [1, -1]]])
A_title = ['Center', 'Stable Spiral']

fig = plt.figure(figsize=(4, 8))

for i in range(2):

    ax = fig.add_subplot(2, 1, i +1)
    
    xdot = A[i,0,0]*x +A[i,0,1]*y
    ydot = A[i,1,0]*x +A[i,1,1]*y

    lamda, v = np.linalg.eig(A[i])
            
    ax.streamplot(x, y, xdot, ydot, color='blue')

    ax.set_title(A_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 5))
    ax.set_ylabel(r'$y$')
    
    ax.yaxis.set_ticks(np.linspace(-w, w, 5))
    if i == 1: ax.set_xlabel(r'$x$')
```

![](figures/ch11_figure04.png){width="75%" fig-align="center"}

For a linear plane autonomous system $\mathbf{x}'=\mathbf{A}\mathbf{x}$ with $\mathrm{det}\,\mathbf{A}\neq 0$,
let $\mathbf{x}$ denote the solution that satisfies the initial condition $\mathbf{x}(0)=\mathbf{x}_0$, where
$\mathbf{x}_0\neq\mathbf{0}$

1. $\lim_{t \to \infty}\mathbf{x}(t)=\mathbf{0}$ if and only if the eigenvalues of $\mathbf{A}$ have negative
   real parts. This will occur when $\Delta>0$ and $\tau<0$
1. $\mathbf{x}(t)$ is periodic if and only if the eigenvalues of $\mathbf{A}$ are pure imaginary. This will
   occur when $\Delta>0$ and $\tau=0$
1. In all other cases, given any neighborhood of the origin, there is at least one $\mathbf{x}_0$ in the neighborhood for which $\mathbf{x}(t)$ becomes unbounded as $t$ increases

## Linearization and Local Stability {#sec-11-3}

Here we will use linearization as a means of analyzing nonlinear DEs and nonlinear systems; the idea is to replace them by linear DEs and linear systems. Let $\mathbf{x}_1$ be a critical point of an autonomous system, and let $\mathbf{x}=\mathbf{x}(t)$ denote the solution that satisfies the initial condition $\mathbf{x}(0)=\mathbf{x}_0$, where $\mathbf{x}\neq\mathbf{x}_1$.

* $\mathbf{x}_1$ is a **stable critical point**
  
  when, given any $\rho > 0$, there is a $r>0$ such that if $\mathbf{x}_0$ satisfies $|\mathbf{x}_0 -\mathbf{x}_1|<r$, then $\mathbf{x}(t)$ satisfies $|\mathbf{x}(t) -\mathbf{x}_1|<\rho$ $\,$ for all $t>0$. If, in addition, 
  $\lim_{t \to \infty} \mathbf{x}(t)=\mathbf{x}_1$ whenever $|\mathbf{x}_0 -\mathbf{x}_1|<r$, we call
  $\mathbf{x}_1$ an **asymptotically stable critical point**
  
* $\mathbf{x}_1$ is an **unstable critical point**

  if there is $\rho>0$ with the property that, for any $r>0$, there is at least one $\mathbf{x}_0$ that
  satisfies $|\mathbf{x}_0 -\mathbf{x}_1|<r$, yet the corresponding solution $\mathbf{x}(t)$ satisfies 
  $|\mathbf{x}(t) -\mathbf{x}_1|\geq\rho$ $\,$ for at least one $t>0$

**Example** $~$ Show that $(0,0)$ is a stable critical point of the nonlinear plane autonomous system

$$
x' = -y -x\sqrt{x^2 +y^2}, \;\; y' = x -y\sqrt{x^2 +y^2}$$

From the formulas $r^2=x^2 +y^2$ and $\theta=\tan^{-1}(y/x)$, we obtain

$$
  \frac{dr}{dt}=\frac{1}{r}\left(x\frac{dx}{dt}+y\frac{dy}{dt} \right),\;
  \frac{d\theta}{dt}=\frac{1}{r^2}\left( -y\frac{dx}{dt} +x\frac{dy}{dt}\right)$$
  
Substituting for $dx/dt$ and $dy/dt$ in the expressions for $dr/dt$ and $d\theta/dt$, we obtain

$$\frac{dr}{dt} = -r^2, \; \frac{d\theta}{dt} = 1$$

Using separation of variables, we see that the solution of the system is

$$ r=\frac{1}{t +c_1},\;\;\theta=t +c_2 $$

for $r\neq 0$. If $\mathbf{x}(0)=(r_0,\theta_0)$ is the initial condition in polar coordinates, then

$$ r=\frac{r_0}{r_0 t +1},\;\;\theta=t +\theta_0 $$

Note that $r\leq r_0$ $\,$ for $t\geq 0$, and $r$ approaches $0$ as $t$ increases. Therefore, given $\rho>0$, a solution that starts less than $\rho$ from the origin remains within $\rho$ of the origin for all $t\geq 0$. Hence the critical point $(0,0)$ is asymptotically stable

**Example** $~$ When expressed in polar coordinates, a plane autonomous system takes the form

$$
\begin{aligned}
  \frac{dr}{dt} &= 0.05\,r(3 -r)\\ 
  \frac{d\theta}{dt} &= -1
\end{aligned}
$$

Show that $(x,y)=(0,0)$ is unstable critical point

We see that $dr/dt=0$ when $r=0$ and can conclude that $(x,y)=(0,0)$ is a critical point. $~$ The differential equation can be solved using separation of variables. If $r(0)=r_0$ and $r_0\neq 0,$ then

$$ r=\frac{3}{1 +c_0 e^{-0.15t}} $$

where $c_0=(3 -r_0)/r_0$. Since

$$\lim_{t \to \infty}\frac{3}{1 +c_0 e^{-0.15t}}=3$$

it follows that no matter how close to $(0,0)$ a solution starts, the solution will leave a disk of radius $\epsilon$ about the origin. Therefore $(0,0)$ is an unstable critical point


* **Linearization**

We replace the term $\mathbf{g}(\mathbf{x})$ in the original autonomous system $\mathbf{x}'=\mathbf{g}(\mathbf{x})$ by a linear term $\mathbf{A}(\mathbf{x} -\mathbf{x}_1)$ that most closely approximates $\mathbf{g}(\mathbf{x})$ in a neighborhood of $\mathbf{x}_1$. This replacement process is called **linearization**

When $\mathbf{x}_1$ is a critical point of a plane autonomous system

$$
\begin{aligned}
 x' &= P(x,y)\\ 
 y' &= Q(x,y)
\end{aligned}
$$

$P(x_1,y_1)=Q(x_1,y_1)=0$ and we have

$$
\begin{aligned}
 x' &= P(x,y)\simeq \left.\frac{\partial P}{\partial x}\right|_{(x_1,y_1)}(x -x_1)
         +\left.\frac{\partial P}{\partial y}\right|_{(x_1,y_1)}(y -y_1)\\ 
 y' &= Q(x,y)\simeq \left.\frac{\partial Q}{\partial x}\right|_{(x_1,y_1)}(x -x_1)
         +\left.\frac{\partial Q}{\partial y}\right|_{(x_1,y_1)}(y -y_1)
\end{aligned}
$$

The original system $\mathbf{x}'=\mathbf{g}(\mathbf{x})$ may be approximated in a neighborhood of $\mathbf{x}_1$
by the linear system $\mathbf{x}'=\mathbf{A}(\mathbf{x} -\mathbf{x}_1)$, where

$$\mathbf{A}=
\begin{pmatrix}
 \left.\frac{\partial P}{\partial x}\right|_{(x_1,y_1)} & \left.\frac{\partial P}{\partial y}\right|_{(x_1,y_1)}\\ 
 \left.\frac{\partial Q}{\partial x}\right|_{(x_1,y_1)} & \left.\frac{\partial Q}{\partial y}\right|_{(x_1,y_1)}
\end{pmatrix} = \mathbf{g}'(\mathbf{x}_1)
$$

* **Stability Criteria for Plane Autonomous System**

1. If the eigenvalues of $\mathbf{A}=\mathbf{g}'(\mathbf{x}_1)$ have negative real part, then $\mathbf{x}_1$ is 
   an asymptotically stable critical point
1. If the eigenvalues of $\mathbf{A}=\mathbf{g}'(\mathbf{x}_1)$ have positive real part, then $\mathbf{x}_1$ is 
   an unstable critical point

**Example** $~$ Classify the criticl points of the following plane autonomous system

$$
\begin{aligned}
 x'&= x^2 +y^2 -6\\ 
 y'&= x^2 -y 
\end{aligned}
$$

The critical points are $(\sqrt{2},2)$ and $(-\sqrt{2},2)$ and the Jacobian matrix is

$$
\mathbf{g}'(\mathbf{x})=
\begin{pmatrix}
 2x & \;2y\\ 
 2x &-1 
\end{pmatrix}
$$

and so

$$
\begin{aligned}
    \mathbf{A}_1 &=\mathbf{g}'({\scriptsize(\sqrt{2},2)})=
    {\scriptsize\begin{pmatrix}
    2\sqrt{2} & \;\;4\\ 
    2\sqrt{2} &-1 
    \end{pmatrix}} \\
    \mathbf{A}_2 &=\mathbf{g}'({\scriptsize(-\sqrt{2},2)})={\scriptsize\begin{pmatrix}
    -2\sqrt{2} & \phantom{-}4\\ 
    -2\sqrt{2} &-1 
    \end{pmatrix}}
\end{aligned}
$$

Since the determinant of $\mathbf{A}_1$ is negative, $\mathbf{A}_1$ has a positive real eigenvalue.
Therefore $(\sqrt{2},2)$ is an unstable critical point. $\mathbf{A}_2$ has a positive determinant and 
a negative trace, and so both eigenvalues have negative real parts. It follows that $(-\sqrt{2},2)$ is a
stable critical point

* **Classifying Critical Points**

![](figures/ch11_figure05.png){width="75%" fig-align="center"}

**Example** $~$ The second-order differential equation $mx'' +kx +k_1x^3 =0$, for $k>0$, represents a general
model for the free, undamped oscillations of a mass $m$ attached to a nonlinear spring.
If $k=1$ and $k_1=-1$, the spring is called soft and the plane autonomous system corresponding to the 
nonlinear second-order equation $x'' +x -x^3=0$ is

$$
\begin{aligned}
 x'&= y\\ 
 y'&= x^3 -x
\end{aligned}
$$

Find and classify the critical points

Since $x^3 -x =x(x^2 -1)$, the critical points are $(0,0)$, $(1,0)$, and $(-1,0)$. The corresponding
Jacobian matrices are

$$
\mathbf{A}_1=\mathbf{g}'((0,0))=
\begin{pmatrix}
 \;\;0 & 1\\ 
 -1 & 0  
\end{pmatrix}\;
$$ 

and

$$
\mathbf{A}_2=\mathbf{g}'((1,0))=\mathbf{g}'((-1,0))=
\begin{pmatrix}
  0 & 1\\ 
  2 & 0 
\end{pmatrix}
$$

Since $\mathrm{det}\,\mathbf{A}_2<0$, critical points $(1,0)$ and $(-1,0)$ are both saddle points. The eigenvalues
of $\mathbf{A}_1$ are $\pm i$, and the status of the critical points at $(0,0)$ remains in doubt. It may be
either a stable spiral, an unstable spiral, or a center

* **Phase-Plane Method**

The linearization method can provide useful information on the local behavior of solutions near critical points.
It is of little help if we are interested in solutions whose initial condition $\mathbf{x}(0)=\mathbf{x}_0$ is not
close to a critical point or if we wish to obtain a global view of the family of solution curves. The 
**phase-plane method** is based on the fact that

$$\frac{dy}{dx}=\frac{dy/dt}{dx/dt}=\frac{Q(x,y)}{P(x,y)}$$

and it attempts to find $y$ as a function of $x$ using one of the methods available for solving first-order DEs

**Example**  $~$ Use the phase-plane method to determine the nature of the solutions to $x'' +x -x^3 =0$ in a neighborhood of $(0,0)$

If we let $dx/dt=y$, then $dy/dt=x^3 -x$. From this we obtain the first-order DE

$$\frac{dy}{dx}=\frac{x^3 -x}{y} $$

which can be solved by separation of variables. Integrating gives

$$ y^2=\frac{x^4}{2}-x^2+c $$

If $\mathbf{x}(0)=(x_0,0)$, where $0<x_0<1$, then $c=-\frac{x_0^4}{2} +x_0^2$, and so

$$ y^2=\frac{(2 -x^2 -x_0^2)(x_0^2 -x^2)}{2} $$

Note that $y=0$ when $x=-x_0$. In addition, the right-hand side is positive when $-x_0<x<x_0$, and
so each $x$ has two corresponding values of $y$. The solution $\mathbf{x}(t)$ that satisfies $\mathbf{x}(0)=(x_0,0)$ is therefore periodic, and so $(0,0)$ is a center

```{python}
from scipy.integrate import solve_ivp

def myODE(t, y):
    return [y[1], y[0]**3 -y[0]]

def plotPhasePortrait(ax, g, td, x0, y0, color):

    for i in range(x0.shape[0]):
        sol = solve_ivp(g, [0, td[-1]], 
               [x0[i], y0], t_eval = td)    
        ax.plot(x0[i], y0, color +'x')
        ax.plot(sol.y[0], sol.y[1], color +'-', lw=0.2)    

td = np.linspace(0, 50, 2000)    
        
x0_0, y0_0 = np.linspace( 0.1, 1.0, 10),  0.0
x0_1, y0_1 = np.linspace( 1.6, 2.6, 30), -3.0  
x0_2, y0_2 = np.linspace(-2.6,-1.6, 30),  3.0

fig, ax = plt.subplots(figsize=(6, 6))

plotPhasePortrait(ax, myODE, td, x0_0,  y0_0, 'k')
plotPhasePortrait(ax, myODE, td, x0_1,  y0_1, 'b')
plotPhasePortrait(ax, myODE, td, x0_2,  y0_2, 'r')
    
ax.set_title(r'Phase Portrait: $\times$ - initial conditions')

w = 3
ax.axis((-w, w, -w, w))
ax.set_aspect(aspect='equal')
ax.grid()
   
ax.tick_params(axis='both', direction='in', pad=5)     
ax.xaxis.set_ticks(np.linspace(-w, w, 5))
ax.yaxis.set_ticks(np.linspace(-w, w, 5))
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')

plt.show()
```

## Autonomous Systems as Mathematical Models {#sec-11-4}

Many applications from physics and biology give rise to nonlinear autonomous second-order DEs $x''=g(x,x')$

* Nonlinear Pendulum

  $$ \frac{d^2\theta}{dt^2} +\frac{g}{l}\sin\theta=0 $$

  When we let $x=\theta$ and $y=\theta'$, this second-order differential equation may be rewritten as the dynamical system

  $$
    \begin{aligned}
    x' &= y\\ 
    y' &= -\frac{g}{l}\sin x
    \end{aligned}
  $$

  The critical points are $(\pm k\pi,0)$, and the Jacobian matrix is easily shown to be

  $$\mathbf{g}'((\pm k \pi,0))=
    \begin{pmatrix}
    0 & 1\\ 
    (-1)^{k+1}\frac{g}{l} & 0 
    \end{pmatrix}
  $$

  If $k=2n+1$, $\Delta<0$, and so all critical points $(\pm(2n+1)\pi,0)$ are saddle points. When $k=2n$, the eigenvalues are pure imaginary, and so the nature of these critical points remains in doubt. From

  $$\frac{dy}{dx}=\frac{dy/dt}{dx/dt}=-\frac{g}{l}\frac{\sin x}{y}$$

  it follows that 

  $$\displaystyle y^2=\frac{2g}{l}\cos x +c $$

  If $\mathbf{x}(0)=(x_0,0)$, then 

  $$ y^2=\frac{2g}{l}(\cos x -\cos x_0)$$

  Note that $y=0$ when $x=-x_0$, and that $(2g/l)(\cos x -\cos x_0)>0$ for $|x|<|x_0|<\pi$. Thus each such $x$ has two corresponding values of $y$, and so the solution $\mathbf{x}(t)$ that satisfies $\mathbf{x}(0)=(x_0,0)$ is periodic. We may conclude that $(0,0)$ is a center. In the case of large initial velocities, the pendulum spins in complete circles about the pivot  

```{python}
w1 = 2*np.pi
w2 = np.pi
xp = np.linspace(-w1, w1, 100)
yp = np.linspace(-w2, w2, 100)
x, y = np.meshgrid(xp, yp)

g = 9.8
l = 8
xdot = y
ydot = -g/l*np.sin(x)

fig, ax = plt.subplots(figsize=(6, 3))
      
ax.streamplot(x, y, xdot, ydot, 
             linewidth=0.5, density=2, color='blue')

ax.set_title('Nonlinear Pendulum')
 
ax.axis((-w1, w1, -w2, w2))
ax.set_aspect(aspect='equal')
ax.grid()
    
ax.tick_params(axis='both', direction='in', pad=5)     
ax.xaxis.set_ticks(np.linspace(-w1, w1, 5))
ax.xaxis.set_ticklabels([r'$-4\pi$', r'$-2\pi$', r'$0$',
                         r'$2\pi$', r'$4\pi$'])
ax.set_xlabel(r'$x$')
ax.yaxis.set_ticks(np.linspace(-w2, w2, 3))
ax.yaxis.set_ticklabels([r'$-\pi$', r'$0$', r'$\pi$'])
ax.set_ylabel(r'$y$')

plt.show()
```

* Periodic Solutions of the Pendulum DE

  Suppose a bead with mass $m$ slides along a thin wire whose shape is described by the function $z=f(x)$.
  A wide variety of nonlinear oscillations can be obtained by changing the shape of the wire and by making
  different assumptions about the forces acting on the bead

  ![](figures/ch11_figure06.png){width="55%" fig-align="center"}

  The tangential force $\mathbf{F}$ due to the weight $W=mg$ $\,$ has a magnitude $mg\sin\theta$, and therefore
  the $x$-component of $\mathbf{F}$ is $F_x=-mg\sin\theta\cos\theta$. Since $f'(x)=\tan\theta$, $\,$ we can use the identity $1+\tan^2\theta=\sec^2\theta$ to conclude that

  $$ 
  \begin{aligned}
   F_x &= -mg\sin\theta\cos\theta
   =-mg\tan\theta\cos^2\theta \\
   &=-mg\frac{f'(x)}{1 +\left[f'(x)\right]^2}
  \end{aligned} 
  $$

  A damping force $\mathbf{D}$, in the direction opposite to the motion, is a constant multiple of the velocitcy 
  of the bead. The $x$-component of $\mathbf{D}$ is therefore

  $$ D_x=-\beta\frac{dx}{dt} $$

  It follows from Newton's second law that

  $$mx''=-mg\frac{f'(x)}{1 +\left[f'(x)\right]^2}-\beta x' $$

  and the corresponding plane autonomous system is

  $$
    \begin{aligned}
    x' &= y\\ 
    y' &= -g\frac{f'(x)}{1 +\left[f'(x)\right]^2}-\frac{\beta}{m} y
    \end{aligned}
  $$

  If $\mathbf{x}_1=(x_1,y_1)$ is a critical point of the system, $y_1=0$ and therefore $f'(x_1)=0$. 
  The bead must therefore be at rest at a point on the wire where the tangent line is horizontal.
  When $f$ is twice differentiable, the Jacobian matrix at $\mathbf{x}_1$ is

  $$
    \mathbf{g}'(\mathbf{x}_1)=
    \begin{pmatrix}
    0 & \phantom{-}1\\ 
    -g f''(x_1) & -\frac{\beta}{m} 
    \end{pmatrix}
  $$

  and so 

  $$
    \begin{aligned}
    \tau &= 
    -\beta/m \\ 
    \Delta &= gf''(x_1) \\ 
    \tau^2 -4\Delta &= \beta^2/m^2 -4gf''(x_1) 
    \end{aligned}
  $$

  We can make the following conclusions:

  * $f''(x_1)<0$ : $~$ A relative maximum therefore occurs at $x=x_1$, and since $\Delta<0$, $\,$ an unstable saddle point occurs at $\mathbf{x}_1=(x_1,0)$  

  * $f''(x_1)>0$ and $\beta>0$: $~$ A relative maximum therefore occurs at $x=x_1$, and since $\tau<0$ and $\Delta>0$, $\,\mathbf{x}_1=(x_1,0)$ is a stable critical point. 

    If $\beta^2 > 4gm^2f''(x_1)$, the system is **overdamped** and the critical point is a stable node. If $\beta^2 < 4gm^2f''(x_1)$, the system is **underdamped** and the critical point is a stable spiral point. The exact nature of the stable critical point is still in doubt if $\beta^2 = 4gm^2f''(x_1)$

  * $f''(x_1)>0$ and the system is undamped ($\beta=0$):

    In this case the eigenvalues are pure imaginary, but the phase-plane method can be used to show that the critical point is a center. Therefore solutions with $\mathbf{x}(0)=(x(0),x'(0))$ near $\mathbf{x}_1=(x_1,0)$ are periodic

**Example** $~$ A bead is released from the position $x(0)=x_0$ on the curve $z=\cosh(x)$ with initial velocity $x'(0)=\nu_0$. Use the phase-plane method to show that the resulting solution is periodic when the system is undamped

```{python}
w = 10
xp = np.linspace(-w, w, 200)
yp = np.linspace(-w, w, 200)
x, y = np.meshgrid(xp, yp)

g = 9.8
m = 1
beta = [0, 0.25, 0.5]

sb_title = [r'$\beta=0$', r'$\beta=0.25$', r'$\beta=0.5$']

def df(x):
    return np.sinh(x)

fig = plt.figure(figsize=(5, 15))

for i in range(3):
    
    xdot = y
    ydot = -g*df(x) /(1 +df(x)**2) -beta[i] /m *y

    ax = fig.add_subplot(3, 1, i +1)
    ax.streamplot(x, y, xdot, ydot, density=1.5, 
                  linewidth=0.5, color='blue')

    ax.set_title(sb_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 5))
    ax.set_ylabel(r'$y$')
    ax.yaxis.set_ticks(np.linspace(-w, w, 5))
    if i == 2:
        ax.set_xlabel(r'$x$')
    
plt.show()    
```

* Lotka-Volterra Predator-Prey Model

  There are many predator-prey models that lead to plane autonomous systems with at least one periodic solution.
  The first such model was constructed independently by pioneer biomathematicians A. Lotka (1925) and
  V. Volterra (1926). If $x$ denotes the number of predators and $y$ denotes the number of prey, then
  the Lotka-Volterra model takes the form

  $$
    \begin{aligned}
    x' &= -ax +bxy = x(-a +by)\\ 
    y' &= -cxy +dy = y(-cx +d)
    \end{aligned}
  $$

  where $a$, $b$, $c$, and $d$ are positive constants. The critical points of this plane autonomous system are
  $(0,0)$ and $(d/c,a/b)$, and the corresponding Jocobian matrices are

  $$
  \begin{aligned}
    \mathbf{A}_1 &= \mathbf{g}'\left((0,0)\right)=
      \begin{pmatrix} -a & 0 \\ \phantom{-}0 & d \end{pmatrix} \\
    \mathbf{A}_2 &= \mathbf{g}'\left((d/c,a/b)\right)=
     \begin{pmatrix} 0 & bd/c \\ -ac/b & 0 \end{pmatrix}
  \end{aligned}
  $$

  The critical point $(0,0)$ is a saddle point. Since $\mathbf{A}_2$ has pure imaginary eigenvalues 
  $\lambda=\pm i\sqrt{ad}$, the critical point $(d/c,a/b)$ may be a center. This possibility can be
  investigated using the phase-plane method. Since

  $$ \frac{dy}{dx}=\frac{y(-cx +d)}{x(-a +by)} $$

  we may separate variables and obtain

  $$ \int\frac{-a +by}{y}\,dy=\int\frac{-cx +d}{x}\,dx $$

  so that

  $$ \left(x^d e^{-cx}\right)\left(y^a e^{-by}\right)=F(x)G(y)=c_0 $$

  Typical graphs of the nonnegative functions $F(x)=x^d e^{-cx}$ and $G(y)=y^a e^{-by}$ are shown in

  ![](figures/ch11_figure07.png){width=85%" fig-align="center"}

  It is not hard to show that $F(x)$ has an absolute maximum at $x=d/c$, whereas $G(y)$ has an absolute maximum at $y=a/b$. These graphs can be used to establish the following properties of a solution curve that orginates at a noncritical point $(x_0,y_0)$ in the first quadrant

  * If $y=a/b$, $F(x)G(a/b)=c_0$ has two solutions $x_m$ and $x_M$ that satisfy $x_m<d/c<x_M$

  * If $x_m<x_1<x_M$ and $x=x_1$, $\,$ then $F(x_1)G(y)=c_0$ has exactly two solutions $y_1$ and $y_2$ that
    satisfy $y_1<a/b<y_2$

  * If $x$ is outside the interval $[x_m,x_M]$, $\,$ then $F(x)G(y)=c_0$ has no solutions

**Example** $~$ If we let $a=0.1$, $b=0.002$, $c=0.0025$, and $d=0.2$ in the Lotka-Volterra predator-prey model, the critical point in the first quadrant is $(d/c,a/b)=(80,50)$, and we know that this critical point is a center. The eigenvalues of $\mathbf{g}'((80,50))$ are $\lambda=\pm(\sqrt{2}/10)i$, and so the solutions near the critical point have period $p\simeq 10\sqrt{2}\pi$

```{python}
a = 0.1
b = 0.002
c = 0.0025
d = 0.2

def myODE(t, y):
    return [-a*y[0] +b*y[0]*y[1], -c*y[0]*y[1] +d*y[1]]

fig, ax = plt.subplots(figsize=(6, 6))

td = 80
t_eval = np.linspace(0, td, 300)
x0 = 0
y0 = 50
for i in range(15):
    x0 = x0 +5
    sol = solve_ivp(myODE, [0, td], [x0, y0], 
         t_eval=t_eval)    
    ax.plot(x0, y0, 'ro')
    ax.plot(sol.y[0], sol.y[1], 'b-', linewidth=1)
    
ax.set_title('Lotka-Volterra model')

w = 400
ax.axis((0, w, 0, w))
ax.set_aspect(aspect='equal')
ax.grid()
   
ax.tick_params(axis='both', direction='in', pad=5)     
ax.xaxis.set_ticks(np.linspace(0, w, 5))
ax.set_xlabel(r'$x$, predator')
ax.yaxis.set_ticks(np.linspace(0, w, 5))
ax.set_ylabel(r'$y$, prey')
plt.show()
```

* Lotka-Volterra Competition Model

  A **competitive interaction** occurs when two or more species compete for the food, water, light, and species
  resources of an ecosystem. A number of mathematical models have been constructed that offer insights into conditions that permit *coexistence*. If $x$ denotes the number in species I and $y$ denotes the number in species II, then the Lotka-Volterra model takes the form

  $$
    \begin{aligned}
    x' &= \frac{r_1}{K_1} x(K_1 -x -\alpha_{12}y)\\ 
    y' &= \frac{r_2}{K_2} y(K_2 -y -\alpha_{21}x)
    \end{aligned}
  $$

  Note that in the absence of species II ($y=0$), $x'=(r_1/K_1)x(K_1 -x)$, and so the first population grows logistically and approaches the steady-state population $K_1$. A similar statement holds for species II growing in the absence of speices I. The term $-\alpha_{21}xy$ in the second equation stems from the competitive effect
  of species I on species II. The model therefore assumes that this rate of inhibition is directly proportional to
  the number of possible competitive pairs $xy$ at a particular $t$

  This plane autonomous system has critical points at $(0,0)$, $(K_1,0)$, and $(0,K_2)$. When 
  $\alpha_{12}\alpha_{21}\neq 0$, the lines $K_1 -x -\alpha_{12}y=0$ and $K_2 -y -\alpha_{21}x=0$ intersect to produce a fourth critical point $\hat{\mathbf{x}}=(\hat{x},\hat{y})$

  $$ \hat{x}=\frac{K_1 -\alpha_{12}K_2}{1 -\alpha_{12}\alpha_{21}},\;
  \hat{y}=\frac{K_2 -\alpha_{21}K_1}{1 -\alpha_{12}\alpha_{21}}$$

**Example** $~$ A competitive interaction is described by the Lotka-Volterra competition model

$$
\begin{aligned}
 x' &= 0.004 x\,(50 -x -0.75y)\\ 
 y' &= 0.001 y\,(100 -y -3.0x)
\end{aligned}
$$

Find and classify all critical points of the system

```{python}
w = 120
xp = np.linspace(0, w, 200)
yp = np.linspace(0, w, 200)
x, y = np.meshgrid(xp, yp)

r1_K1 = 0.004
r2_K2 = 0.001
K1 = 50
K2 = 100
a_12 = np.array([0.75, 0.4])
a_21 = np.array([3.00, 1.5])

fig = plt.figure(figsize=(6, 12))

for i in range(2):
    
    ax = fig.add_subplot(2, 1, i +1)
    
    xdot = r1_K1*x*(K1 -x -a_12[i]*y)
    ydot = r2_K2*y*(K2 -y -a_21[i]*x)

    cpt_x = np.array([0, K1, 0, (K1 -a_12[i]*K2)
                      /(1 -a_12[i]*a_21[i])])
    cpt_y = np.array([0, 0, K2, (K2 -a_21[i]*K1)
                      /(1 -a_12[i]*a_21[i])])

    ax.streamplot(x, y, xdot, ydot, density=1.5, 
                  linewidth=0.5, color='blue')
    ax.plot(cpt_x, cpt_y, 'ro')
    ax.set_title(r'$\alpha_{12}=$' +str(a_12[i]) 
                 +r'$, \alpha_{21}=$' +str(a_21[i]))
    ax.axis((0, w, 0, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(0, w, 4))
    ax.set_xlabel(r'$x$')
    ax.yaxis.set_ticks(np.linspace(0, w, 4))
    ax.set_ylabel(r'$y$')
    
plt.show()
```

## Periodic Solutions, Limit Cycles, and Global Stability {#sec-11-5}

In this section, we will investigate the existence of periodic solutions of nonlinear plane autonomous systems
and introduce special periodic solutions called **limit cycles**

An analysis of critical points using linearization can provide valuable information on the behavior of solutions
near critical points and insight into a variety of biological and physical phenomena. However there are some
inherent limitations to this approach. When the eigenvalues of the Jacobian matrix are pure imaginary,
we cannot conclude that there are periodic solutions near the critical point

The first goal of this section is *to determine conditions under which we can either exclude the possibility of periodic solutions or assert their existence*

A second goal is *to determine conditions under which an asymptotically stable critical point is globally 
stable*: $\lim_{t \to \infty}\mathbf{x}(t)=\mathbf{x}_1$ for all initial conditions in 
a simply connected region $\mathbb{R}$

* **Negative Criteria**

**Theorem** If a plane autonomous system has a periodic solution $\mathbf{x}=\mathbf{x}(t)$ in a simply 
connected region $\mathbb{R}$, then the system has at least one critical point inside the corresponding
simple closed curve $C$. If there is a single critical point inside $C$, then that critical point cannot be a 
saddle point

**Corollary** If a simply connected region $\mathbb{R}$ either contains no critical points of a plane autonomous
system or contains a single saddle point, then there are no periodic solutions in $\mathbb{R}$

**Example** $~$ Show that the Lotka-Volterra competition model

$$
\begin{aligned}
 x' &= 0.004 x\,(50 -x -0.75y)\\ 
 y' &= 0.001 y\,(100 -y -3.0x)
\end{aligned}
$$

has no periodic solutions in the first quadrant

Another sometimes useful result can be formulated in terms of the divergence of the vector field 

$$ \mathbf{v}(x,y)=\left(P(x,y), Q(x,y)\right) $$

* **Bendixson Negative Criterion**

If $\displaystyle\mathrm{div}\,\mathbf{v}=\frac{\partial P}{\partial x} +\frac{\partial Q}{\partial y}$ does not
change sign in a simply connected region $\mathbb{R}$, then the plane autonomous system has no periodic solutions in $\mathbb{R}$

**Proof** $~$ Suppose, to the contrary, that there is a periodic solution $\mathbf{x}=\mathbf{x}(t)$ in $\mathbb{R}$, and let $C$ be the resulting simple closed curve and $R_1$ the region bounded by $C$. By using Green's theorem, we obtain

$$ \int_C -Q(x,y)\,dx +P(x,y)\,dy=\iint_{R_1} 
  \left(\frac{\partial P}{\partial x} +\frac{\partial Q}{\partial y}\right)\, dxdy $$

Since $\mathbf{x}=\mathbf{x}(t)$ is a solution with period $p$, $\,$ we have $x'(t)=P(x(t),y(t))$ and
$y'(t)=Q(x(t),y(t))$, and so

$$
\begin{aligned}
 \int_C -&Q(x,y)\,dx +P(x,y)\,dy \\
   &= \int_0^p \left[-Q(x(t),y(t))\,x'(t) +P(x(t),y(t))\,y'(t) \right]\,dt\\ 
   &= \int_0^p -QP +PQ \,dt = 0
\end{aligned}
$$

Since $\displaystyle\mathrm{div}\,\mathbf{v}=\frac{\partial P}{\partial x} +\frac{\partial Q}{\partial y}$ is continuous and does not change sign in $\mathbb{R}$, it follows that either $\mathrm{det}\,\mathbf{v}\geq 0$
or $\mathrm{det}\,\mathbf{v}\leq 0$ in $\mathbb{R}$, and so

$$
  \iint_{R_1} \left(\frac{\partial P}{\partial x} +\frac{\partial Q}{\partial y}\right)\, dxdy \neq 0
$$

This contradiction establishes that there are no periodic solutions in $\mathbb{R}$

**Example** $~$ The sliding bead satisfies the differential equation

$$ mx''=-mg\frac{f'(x)}{1 +\left[f'(x)\right]^2}-\beta x' $$

Show that there are no periodic solutions

**Solution** $~$ The corresponding plane autonomous system is

$$
\begin{aligned}
   x' &= y\\ 
   y' &= -g\frac{f'(x)}{1 +\left[f'(x)\right]^2}-\frac{\beta}{m} y
\end{aligned}
$$

and so $\displaystyle
  \mathrm{div}\,\mathbf{v}=\frac{\partial P}{\partial x} +\frac{\partial Q}{\partial y}=-\frac{\beta}{m}<0$

* **Dulac Negative Criterion**

If $\delta(x,y)$ has continuous first partial derivatives in a simply connected region $\mathbb{R}$ and
$\displaystyle \frac{\partial (\delta P)}{\partial x} +\frac{\partial (\delta Q)}{\partial y}$ does not
change sign in $\mathbb{R}$, $\,$ then the plane autonomous system has no periodic solutions in $\mathbb{R}$

**Example** $~$ Use $\delta(x,y)=1/(xy)$ to show that the Lotka-Volterra competition equations

$$
\begin{aligned}
 x' &= \frac{r_1}{K_1} x(K_1 -x -\alpha_{12}y)\\ 
 y' &= \frac{r_2}{K_2} y(K_2 -y -\alpha_{21}x)
\end{aligned}
$$

have no periodic solutions in the first quadrant

**Solution** $~$ If $\delta(x,y)=1/(xy)$, then

$$ \frac{\partial (\delta P)}{\partial x} +\frac{\partial (\delta Q)}{\partial y}=
  \frac{r_1}{K_1}\left(-\frac{1}{y}\right) +\frac{r_2}{K_2}\left(-\frac{1}{x}\right) <0 $$

in the first quadrant

* **Positive Criteria**

  A region $\mathbb{R}$ is called an **invariant region** for a plane autonomous system if, whenever $\mathbf{x}_0$ is in $\mathbb{R}$, the solution $\mathbf{x}=\mathbf{x}(t)$ remains in $\mathbb{R}$. If $\mathbf{n}(x,y)$ denotes a normal vector on the boundary that points inside the region, then $\mathbb{R}$ will be an invariant region for the plane autonomous system provided $\mathbf{v}(x,y)\cdot\mathbf{n}(x,y)\geq 0$ for all points $(x,y)$ in the boundary

 ![](figures/ch11_figure08.png){width="45%" fig-align="center"}

* **Poincar&eacute;-Benedixon I**

  Let $\mathbb{R}$ be an invariant region for a plane autonomous system and suppose that $\mathbb{R}$ has no
  critical points on its boundary

  $(\text{a})$ $~$ If $\mathbb{R}$ is a type I region that has a single unstable node or an unstable spiral point in its interior, then there is at least one periodic solution in $\mathbb{R}$

  $(\text{b})$ $~$ If $\mathbb{R}$ is a type II region that contains no critical points of the system, 
   then there is at least one periodic solution in $\mathbb{R}$
    
  In either of the two cases, if $\mathbf{x}=\mathbf{x}(0)$ is a nonperiodic solution in $\mathbb{R}$, then
  $\mathbf{x}(t)$ spirals toward a cycle that is a solution to the system. This periodic solution is called a **limit cycle**

**Example** $~$ The Van der Pol equation is a nonlinear second-order differential equation that arises in electronics, and as a plane autonomous system it takes the form

$$\begin{aligned}
 x' &= y\\ 
 y' &= -\mu (x^2 -1)y -x
\end{aligned}$$

Show that the Van der Pol DE has a periodic solution when $\mu>0$

```{python}
w = 3
xp = np.linspace(-w, w, 100)
yp = np.linspace(-w, w, 100)
x, y = np.meshgrid(xp, yp)

fig = plt.figure(figsize=(5, 10))

mu = np.array([1, 3])
vdp_title = [r'$\mu=1$', r'$\mu=3$']

for i in range(2):
    
    ax = fig.add_subplot(2, 1, i +1)
    
    xdot = y
    ydot = -mu[i]*(x*x -1)*y -x

    ax.streamplot(x, y, xdot, ydot, density=1.5,  
       linewidth=0.5, color='blue')
    ax.plot(0, 0, 'ro')
    ax.set_title(vdp_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 4))
    ax.set_xlabel(r'$x$')
    ax.yaxis.set_ticks(np.linspace(-w, w, 4))
    ax.set_ylabel(r'$y$')
    
plt.show()
```

We will assumme that there is a Type I invariant region $\mathbb{R}$ for the corresponding plane autonomous system and that this region contains $(0,0)$ in its interior. The only critical point is $(0,0)$, and the Jacobian matrix is given by

$$
\mathbf{g}'((0,0))=
\begin{pmatrix}
 \;\;0 & 1\\ 
 -1 & \mu
\end{pmatrix}
$$

Therefore, $\tau=\mu$, $\Delta=1$, and $\tau^2 -4\Delta=\mu^2-4$. Since $\mu >0$, the critical point is either an unstable spiral point or an unstable node. By part (a) of **Poincar&eacute;-Benedixon I**, the system has at least one periodic solution in $\mathbb{R}$

```{python}
td = 15
w = 5
mu_ = np.array([1, 2.5])
x0 = np.array([[0.5, 0.5], [3, 3], [-3,-3], [0, 0.1]])

def myODE(t, y):
    return [y[1], -mu*(y[0]*y[0] -1)*y[1] -y[0]]

vdp_title = [r'$\mu=1$', r'$\mu=2.5$']
vdp_line = ['b:', 'r:', 'k:', 'g:']

fig = plt.figure(figsize=(5, 10))

for i in range(2):

    ax = fig.add_subplot(2, 1, i +1)
    
    mu = mu_[i]
    for j in range(4):
        y0 = x0[j]
        t_eval = np.linspace(0, td, 500)
        sol = solve_ivp(myODE, [0, td], y0, t_eval=t_eval)    
        ax.plot(sol.y[0], sol.y[1], vdp_line[j], 
            linewidth=1)
        ax.plot(y0[0], y0[1], 'ro')
    
    ax.set_title(vdp_title[i])

    ax.axis((-w, w, -w, w))
    ax.set_aspect(aspect='equal')
    ax.grid()
    
    ax.tick_params(axis='both', direction='in', pad=5)     
    ax.xaxis.set_ticks(np.linspace(-w, w, 4))
    ax.set_xlabel(r'$x$')
    ax.yaxis.set_ticks(np.linspace(-w, w, 4))
    ax.set_ylabel(r'$y$')
    
plt.show()
```

* **Poincar&eacute;-Benedixon II**

  Let $\mathbb{R}$ be a type I invariant region for a plane autonomous system that has no periodic solutions in $\mathbb{R}$

  $(\text{a})$ $~$ If $\mathbb{R}$ has a finite number of nodes or spiral points,
    then given any initial position $\mathbf{x}_0$ in $\mathbb{R}$, $\lim_{t \to \infty} \mathbf{x}(t)=\mathbf{x}_1$ for some critical point $\mathbf{x}_1$

  $(\text{b})$ $~$ If $\mathbb{R}$ has a single stable node or stable spiral point $\mathbf{x}_1$ in its interior and no critical points on its boundary, then $\lim_{t \to \infty} \mathbf{x}(t)=\mathbf{x}_1$
    for all initial positions $\mathbf{x}_0$ in $\mathbb{R}$: **globally stable**

**Example** $~$ Empirical evidence suggests that the plane autonomous system

$$
\begin{aligned}
 x' &= x^2y -x +1\\ 
 y' &= -x^2y +\frac{1}{2} 
\end{aligned}
$$

has a type I invariant region $\mathbb{R}$ that lies inside the reactangle defined by $0\leq x \leq 2$, $0 \leq y \leq 1$

$(\text{a})$ $~$ Use the Benedixson negative criterion to show that there are no periodic solutions in R

$(\text{b})$ $~$ If $\mathbf{x}_0$ is in $\mathbb{R}$ and $\mathbf{x}=\mathbf{x}(t)$ is the solution satisfying 
$\mathbf{x}(t)=\mathbf{x}_0$, use **Poincar&eacute;-Benedixon II** to find $\lim_{t \to \infty}\mathbf{x}(t)$

## Worked Exercises {.unnumbered}