In [2]:

using Interact
using Reactive
using Sundials
import PyPlot



# Problem 1

Draw the phase portrait for the system

\begin{align} \dot{x} & = x \left( x^2 + y^2 - 2 \right)\\ \dot{y} & = - y \left( x^2 + y^2 - 3 x + 1 \right)\\ \end{align}


In [3]:

X = linspace(-3,3,100) * ones(100)';
Y = ones(100) * linspace(-3,3,100)';
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0));






That wasn't particularly enlightening. I think that we should analyze it a bit more.

## A Certain Symmetry

Suppose we substitute $-y$ for $y$ whenever we see it in the equations of motion.

\begin{align} \dot{x} & = x \left( x^2 + (-y)^2 - 2 \right)\\ \dot{-y} & = -( -y) \left( x^2 + (-y)^2 - 3 x + 1 \right)\\ \end{align}

which reduces to:

\begin{align} \dot{x} & = x \left( x^2 + y^2 - 2 \right)\\ \dot{y} & = - y \left( x^2 + y^2 - 3 x + 1 \right)\\ \end{align}

We conclude that the system is symmetric about the x-axis.

## Nullclines

We search for x-nullclines. If $0 = x \left( x^2 + y^2 - 2 \right)$, then either $x = 0$ or $x^2 + y^2 = 2$.



In [4]:

PyPlot.plot(zeros(100),linspace(-3,3,100),color="red"); PyPlot.grid(true);
PyPlot.plot(sqrt(2)*cos(linspace(0,2*π)), sqrt(2)*sin(linspace(0,2*π)),color="red")
PyPlot.title("X-nullclines");






We search for y-nullclines. If $0 = -y(x^2 + y^2 - 3x + 1)$, then either $y=0$ or $x^2 + y^2 - 3x + 1 = 0$. Completing the square, we find that this implies $(x - \frac{3}{2})^2 + y^2 = \frac{5}{4}$. Therefore, we can plot the y-nullclines:



In [5]:

PyPlot.plot(linspace(-3,3,100),zeros(100),color="blue"); PyPlot.grid(true);
PyPlot.plot(sqrt(5.0/4.0)*cos(linspace(0,2*π)) + 1.5, sqrt(5.0/4.0)*sin(linspace(0,2*π)),color="blue")
PyPlot.title("Y-nullclines");






Overlaying the two sets of nullclines onto one graph, we see:



In [6]:

PyPlot.plot(linspace(-3,3,100),zeros(100),color="blue"); PyPlot.grid(true);
PyPlot.plot(sqrt(5.0/4.0)*cos(linspace(0,2*π)) + 1.5, sqrt(5.0/4.0)*sin(linspace(0,2*π)),color="blue")
PyPlot.plot(zeros(100),linspace(-3,3,100),color="red"); PyPlot.grid(true);
PyPlot.plot(sqrt(2)*cos(linspace(0,2*π)), sqrt(2)*sin(linspace(0,2*π)),color="red")
PyPlot.scatter([0.0, 1.0, 1.0, sqrt(2.0), -sqrt(2.0)],[0.0, 1.0, -1.0, 0.0, 0.0])
PyPlot.title("Nullclines");
X = linspace(-1.5,2.2,20) * ones(20)';
Y = ones(20) * linspace(-1.5,1.5,20)';
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0))




Out[6]:

PyObject <matplotlib.quiver.Quiver object at 0xa9aa71ec>



We see intersections at $(0,0), (\pm \sqrt{2},0), (1, \pm 1)$.

## Fixed Points

Recall that our system is:

\begin{align} \dot{x} & = x \left( x^2 + y^2 - 2 \right)\\ \dot{y} & = - y \left( x^2 + y^2 - 3 x + 1 \right)\\ \end{align}

We find that it has jacobian:

$$J(x,y) = \left( \begin{array}{cc} 3 x^2 + 2 y^2 - 2 & 2 x y \\ 3 y - 2 y x & 3x - x^2 - 3 y^2 - 1\\ \end{array} \right)$$


In [7]:

function Jacobian(x,y)
[ 3*x^2 + 2*y^2 - 2  2*x.*y; 3*y - 2*x.*y 3*x - x.^2 - 3*y.^2 - 1]
end




Out[7]:

Jacobian (generic function with 1 method)



### The Fixed Point $(0,0)$

We start by analyzing the fixed point at the origin.



In [8]:

Jacobian(0,0)




Out[8]:

2x2 Array{Int32,2}:
-2   0
0  -1



Thus, the jacobian has eigenvector $(1,0)$ with eigenvalue $-2$ and eigenvector $(0,1)$ with eigenvalue $-1$. The origin is a stable fixed point.



In [9]:

X = linspace(-0.1,0.1,10) * ones(10)'; PyPlot.grid(true);
Y = ones(10) * linspace(-0.1,0.1,10)'; PyPlot.title("Around the fixed point (0,0)"); PyPlot.scatter([0.0],[0.0]);
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0)); PyPlot.axis([-0.1,0.1,-0.1,0.1]);






### The Fixed Points $(\pm \sqrt{2}, 0)$

Next, we analyze the intersections of the circle centered at the origin with the x-axis, the fixed points $(\pm \sqrt{2}, 0)$.

Analytically, we determine that $J(\pm \sqrt{2}, 0) = \left( \begin{array}{cc} 4 & 0 \\ 0 & 3(\pm \sqrt{2} - 1)\\ \end{array} \right)$

The eigenvectors here are also $(1,0)$ and $(0,1)$. For $(- \sqrt{2}, 0)$, the fixed point repels in the x-direction but attracts in the y-direction. For $(+ \sqrt{2}, 0)$, the fixed point repels in both directions.



In [10]:

X = linspace(-1.45,-1.35,10) * ones(10)'; PyPlot.grid(true);
Y = ones(10) * linspace(-0.1,0.1,10)'; PyPlot.title("Around the fixed point (- \sqrt{2},0)"); PyPlot.scatter([-sqrt(2.0)],[0]);
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0)); PyPlot.axis([-1.45,-1.35,-0.1,0.1]);







In [11]:

X = linspace(1.35,1.45,10) * ones(10)'; PyPlot.grid(true);
Y = ones(10) * linspace(-0.1,0.1,10)'; PyPlot.title("Around the fixed point (\sqrt{2},0)"); PyPlot.scatter([sqrt(2.0)],[0]);
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0)); PyPlot.axis([1.35,1.45,-0.1,0.1]);






### The fixed points $(1, \pm 1)$

Finally, we analyze the fixed points $(1,\pm 1)$



In [12]:

Jacobian(1,1)




Out[12]:

2x2 Array{Int32,2}:
3   2
1  -2




In [13]:

Jacobian(1,-1)




Out[13]:

2x2 Array{Int32,2}:
3  -2
-1  -2



Notice that these two matrices have the same determinant and trace. $\Delta = -8$, $\tau = 1$. Therefore, they have the same eigenvalues:

$$\lambda_\pm = \frac{1}{2} \pm \sqrt{ \frac{33}{4} }$$

Becuse of the symmetry, the eigenvectors should be mirror images of each other, if the symmetry holds. Specifically, if $(\alpha, \beta)$ is an eigenvector of $J(1,1)$, then the corresponding eigenvector of $J(1,-1)$ should be $(\alpha, -\beta)$.

We find that the eigenvectors of $J(1,1)$ are $\left( \begin{matrix} 4 \\ -5 \pm \sqrt{33} \\ \end{matrix} \right)$ with eigenvalues $\frac{1 \pm \sqrt{33}}{2}$. Consequently, the eigenvectors of $J(1,-1)$ are $\left( \begin{matrix} 4 \\ 5 \mp \sqrt{33} \\ \end{matrix} \right)$ with the same eigenvalues.



In [14]:

X = linspace(0.9,1.1,10) * ones(10)'; PyPlot.grid(true);
Y = ones(10) * linspace(0.9,1.1,10)'; PyPlot.title("Around the fixed point (1, 1)"); PyPlot.scatter([1],[1]);
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0)); PyPlot.axis([0.9,1.1,0.9,1.1]);
T = linspace(-0.2,0.2)
PyPlot.plot(1 + 4*T, 1 + (-5 + sqrt(33.0))*T);
PyPlot.plot(1 + 4*T, 1 + (-5 - sqrt(33.0))*T);







In [15]:

X = linspace(0.9,1.1,10) * ones(10)'; PyPlot.grid(true);
Y = ones(10) * linspace(-1.1,-0.9,10)'; PyPlot.title("Around the fixed point (1, -1)"); PyPlot.scatter([1],[-1]);
PyPlot.quiver(X,Y,X .* (X .^2 + Y .^2 - 2.0), -Y .* (X.^2 + Y.^2 - 3*X + 1.0)); PyPlot.axis([0.9,1.1,-1.1,-0.9]);
T = linspace(-0.2,0.2)
PyPlot.plot(1 + 4*T, -1 + (5 - sqrt(33.0))*T);
PyPlot.plot(1 + 4*T, -1 + (5 + sqrt(33.0))*T);






Notice that the two phase portraits are mirror images of each other.

## Stable Manifolds

We want to find the stable manifolds of the system. These are trajectories which connect fixed points to other fixed points. We can accomplish this by integrating a trajectory starting at every repelling direction.



In [16]:

function problem1(t, X, Xdot)
Xdot[1] = X[1] .* ( X[1].^2 + X[2].^2 - 2.0)
Xdot[2] = -X[2] .* ( X[1] .^2 + X[2] .^2 - 3.0*X[1] + 1.0)
end