{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c1e07947-9b5e-4593-8613-1164417ee953",
   "metadata": {},
   "source": [
    "ME 2016 - Nonlinear Dynamical Systems 1\n",
    "\n",
    "Homework 3\n",
    "\n",
    "**Dane Sabo**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "211982e1-a936-47d8-a01e-9a90cfc9de1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sm\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "from math import *"
   ]
  },
  {
   "attachments": {
    "ae0f41aa-f319-4437-bf03-f540055a8e3d.png": {
     "image/png": "[encoded data, 27836 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "a01a46d0-bc05-4b06-a782-431f6f525c9f",
   "metadata": {},
   "source": [
    "![image.png](attachment:ae0f41aa-f319-4437-bf03-f540055a8e3d.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9db2c23d-bf5d-4c4f-a3ea-78c78aa77569",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x78424e30d670>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 224808 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, y = sm.symbols('x,y')\n",
    "xdot = -y + x*sm.cos(sm.pi*x)\n",
    "ydot = x - y**3\n",
    "\n",
    "#Create Functions\n",
    "P = lambda x, y: -y + x*np.cos(np.pi*x)\n",
    "Q = lambda x, y: x - y**3\n",
    "def dX_dt(X, t, args = (P,Q)):\n",
    "    return [P(X[0], X[1]), Q(X[0], X[1])]\n",
    "\n",
    "#Initial Conditions\n",
    "ts = np.linspace(0, 5, 30)\n",
    "ic = np.linspace(-1.2, 1.2, 7)\n",
    "\n",
    "#Make a vectorfield\n",
    "X, Y = np.mgrid[ic[0]:ic[-1]:20j, ic[0]:ic[-1]:20j]\n",
    "u = P(X,Y)\n",
    "v = Q(X,Y)\n",
    "plt.quiver(X, Y, u, v, color = 'red')\n",
    "\n",
    "#plot trajectories\n",
    "plt.title(f'Phase Portrait for $\\\\dot x = $ {xdot} and $\\\\dot y = $ {ydot}')\n",
    "plt.xlabel('x'); plt.ylabel('y'); plt.grid('both')\n",
    "plt.xlim((ic[0],ic[-1])); plt.ylim((ic[0],ic[-1]))\n",
    "\n",
    "for r in ic:\n",
    "    for s in ic:\n",
    "        X0 = [r,s]\n",
    "        Xs = odeint(dX_dt,X0,ts)\n",
    "        plt.plot(Xs[:,0],Xs[:,1], \"royalblue\",linewidth = 0.75)\n",
    "rect = np.array([[-1, -1], [-1, 1], [1, 1], [1, -1], [-1, -1]])\n",
    "plt.plot(rect[:,0], rect[:,1], '-g')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7d4c7594-806b-4711-9f85-24bfd075763a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dV/dt Expression: \n",
      "-pi*x*sin(pi*x) - 3*y**2 + cos(pi*x)\n",
      "\n",
      "dV/dt Simplified:\n",
      "-pi*x*sin(pi*x) - 3*y**2 + cos(pi*x)\n",
      "\n",
      "Sign Changes: \n",
      "[(x, -sqrt(-3*pi*x*sin(pi*x) + 3*cos(pi*x))/3), (x, sqrt(-3*pi*x*sin(pi*x) + 3*cos(pi*x))/3)]\n"
     ]
    }
   ],
   "source": [
    "#Guess a psi for Poincare - Bendixson Theorem\n",
    "psi = 1\n",
    "def sign_change_test(xdot, ydot, psi):\n",
    "    div_psiX = sm.diff(psi*xdot,x) + sm.diff(psi*ydot, y)\n",
    "    print(f'dV/dt Expression: \\n{div_psiX}')\n",
    "    print(f'\\ndV/dt Simplified:\\n{sm.simplify(div_psiX)}\\n')\n",
    "    \n",
    "    sign_changes = sm.solve(sm.Eq(0,div_psiX),x,y)\n",
    "    print(f'Sign Changes: \\n{sign_changes}')\n",
    "    return div_psiX\n",
    "\n",
    "expr = sign_change_test(xdot, ydot, psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5f3e5930-f283-47f2-b6f7-e88131fec12b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 3 y^{2} - 1$"
      ],
      "text/plain": [
       "-3*y**2 - 1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,-1)\n",
    "#This is clearly always < 0. No sign change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "22f92f82-2559-4cf0-9323-46be12b52a4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\pi x \\sin{\\left(\\pi x \\right)} + \\cos{\\left(\\pi x \\right)} - 3$"
      ],
      "text/plain": [
       "-pi*x*sin(pi*x) + cos(pi*x) - 3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "159811ac-97c6-4676-a6e1-dcdf084d32a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 3 y^{2} - 1$"
      ],
      "text/plain": [
       "-3*y**2 - 1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b2493449-3b2b-43d4-8635-2513c65ca16c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\pi x \\sin{\\left(\\pi x \\right)} + \\cos{\\left(\\pi x \\right)} - 3$"
      ],
      "text/plain": [
       "-pi*x*sin(pi*x) + cos(pi*x) - 3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "21246928-c9c0-4c06-833e-8e8dde8b9b33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 33600 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Let's check graphically the range of -pi*x*sin(pi*x) + cos(pi*x) - 3\n",
    "x_range = np.linspace(-1,1,50)\n",
    "y_range = -np.pi*x_range*np.sin(np.pi*x_range) + np.cos(np.pi*x_range) - 3\n",
    "plt.figure()\n",
    "plt.plot(x_range,y_range)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "#Looks like no sign change for our domain of -1<x<1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6ea1e086-12d8-46ad-b738-8a99fcf5577a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 3 y^{2} - 0.1 \\pi \\sin{\\left(0.1 \\pi \\right)} + \\cos{\\left(0.1 \\pi \\right)}$"
      ],
      "text/plain": [
       "-3*y**2 - 0.1*pi*sin(0.1*pi) + cos(0.1*pi)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "424752e0-5840-4193-9045-9530d3c55899",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 3 y^{2} - 0.1 \\pi \\sin{\\left(0.1 \\pi \\right)} + \\cos{\\left(0.1 \\pi \\right)}$"
      ],
      "text/plain": [
       "-3*y**2 - 0.1*pi*sin(0.1*pi) + cos(0.1*pi)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e7396efe-1d3c-4589-95f7-6b6e7c89b3e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\pi x \\sin{\\left(\\pi x \\right)} + \\cos{\\left(\\pi x \\right)} - 0.03$"
      ],
      "text/plain": [
       "-pi*x*sin(pi*x) + cos(pi*x) - 0.03"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1ed10fd6-1bfa-4053-983a-2ae901a069fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\pi x \\sin{\\left(\\pi x \\right)} + \\cos{\\left(\\pi x \\right)} - 0.03$"
      ],
      "text/plain": [
       "-pi*x*sin(pi*x) + cos(pi*x) - 0.03"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "28b6df19-6fd0-4be5-b993-56be66ffdc4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Range (y or x)')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 47912 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Let's check graphically the range of -pi*x*sin(pi*x) + cos(pi*x) - 3\n",
    "the_range = np.linspace(-0.1,0.1,50)\n",
    "y_range = -3*the_range**2 - 0.1*np.pi*np.sin(0.1*np.pi)+np.cos(0.1*np.pi)\n",
    "x_range = -np.pi*the_range*np.sin(np.pi*the_range)+np.cos(np.pi*the_range)-0.03\n",
    "plt.figure()\n",
    "plt.plot(the_range,x_range, the_range, y_range)\n",
    "plt.legend(['X-Range', 'Y-Range'])\n",
    "plt.xlabel('Domain (x or y)')\n",
    "plt.ylabel('Range (y or x)')\n",
    "\n",
    "#Looks like no sign change for our domain of -1<x<1, always positive (towards limit cycle)"
   ]
  },
  {
   "attachments": {
    "656a4fe2-7f79-46cd-a7ae-600dc6cf338f.png": {
     "image/png": "[encoded data, 15928 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "0b628ced-a69f-44b2-b40a-5a8bd986b1f5",
   "metadata": {},
   "source": [
    "![image.png](attachment:656a4fe2-7f79-46cd-a7ae-600dc6cf338f.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1c5a49fe-5309-46e1-be4a-1684edb7f678",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 0),\n",
       " (sqrt(2)*sqrt((alpha - sqrt(alpha**2 - 2*alpha - 3) + 1)/(alpha - 1))*(-alpha - sqrt(alpha**2 - 2*alpha - 3) + 1)/4,\n",
       "  sqrt(alpha/(2*alpha - 2) - sqrt(alpha**2 - 2*alpha - 3)/(2*alpha - 2) + 1/(2*alpha - 2))),\n",
       " (sqrt(2)*sqrt((alpha - sqrt(alpha**2 - 2*alpha - 3) + 1)/(alpha - 1))*(alpha + sqrt(alpha**2 - 2*alpha - 3) - 1)/4,\n",
       "  -sqrt(alpha/(2*alpha - 2) - sqrt(alpha**2 - 2*alpha - 3)/(2*alpha - 2) + 1/(2*alpha - 2))),\n",
       " (sqrt(2)*sqrt((alpha + sqrt(alpha**2 - 2*alpha - 3) + 1)/(alpha - 1))*(-alpha + sqrt(alpha**2 - 2*alpha - 3) + 1)/4,\n",
       "  sqrt(alpha/(2*alpha - 2) + sqrt(alpha**2 - 2*alpha - 3)/(2*alpha - 2) + 1/(2*alpha - 2))),\n",
       " (sqrt(2)*sqrt((alpha + sqrt(alpha**2 - 2*alpha - 3) + 1)/(alpha - 1))*(alpha - sqrt(alpha**2 - 2*alpha - 3) - 1)/4,\n",
       "  -sqrt(alpha/(2*alpha - 2) + sqrt(alpha**2 - 2*alpha - 3)/(2*alpha - 2) + 1/(2*alpha - 2)))]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha, x, y = sm.symbols('alpha, x, y')\n",
    "P = y + x*(alpha - x**2 - y**2)\n",
    "Q = -x + y*(1-x**2-y**2)\n",
    "\n",
    "Peqn = sm.Eq(P, 0)\n",
    "Qeqn = sm.Eq(Q, 0)\n",
    "\n",
    "criticalpoints = sm.solve((Peqn,Qeqn),x,y)\n",
    "display(criticalpoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0b3c6cb5-e13a-4141-a16a-e8eea1f0aa31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 224528 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xdot = P\n",
    "ydot = Q\n",
    "alpha_val = 0.5\n",
    "\n",
    "#Create Functions\n",
    "P = lambda x, y: y + x*(alpha_val - x**2 - y**2)\n",
    "Q = lambda x, y: -x + y*(1 - x**2 - y**2)\n",
    "def dX_dt(X, t, args = (P,Q)):\n",
    "    return [P(X[0], X[1]), Q(X[0], X[1])]\n",
    "\n",
    "#Initial Conditions\n",
    "ts = np.linspace(0, 5, 30)\n",
    "ic = np.linspace(-1.2, 1.2, 7)\n",
    "\n",
    "#Make a vectorfield\n",
    "X, Y = np.mgrid[ic[0]:ic[-1]:20j, ic[0]:ic[-1]:20j]\n",
    "u = P(X,Y)\n",
    "v = Q(X,Y)\n",
    "plt.quiver(X, Y, u, v, color = 'red')\n",
    "\n",
    "#plot trajectories\n",
    "plt.title(f'Phase Portrait for $\\\\dot x = $ {xdot} and $\\\\dot y = $ {ydot}')\n",
    "plt.xlabel('x'); plt.ylabel('y'); plt.grid('both')\n",
    "plt.xlim((ic[0],ic[-1])); plt.ylim((ic[0],ic[-1]))\n",
    "\n",
    "for r in ic:\n",
    "    for s in ic:\n",
    "        X0 = [r,s]\n",
    "        Xs = odeint(dX_dt,X0,ts)\n",
    "        plt.plot(Xs[:,0],Xs[:,1], \"royalblue\",linewidth = 0.75)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "759a57ff-d8a3-49c1-9556-2d41a34d78b8",
   "metadata": {},
   "source": [
    "From calculating the algabraeic expressions for the critical points and looking at the plot, it's pretty clear there is one critical point inside the potential limit cycle and the rest are outside of the potential limit cycle. \n",
    "\n",
    "It also looks like this limit cycle is **stable**. That is because the vectors point in towards the limit cycle. Let me prove this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4b004f14-5c8e-4eb4-88bd-a93666b11602",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dV/dt Expression: \n",
      "alpha - 4*x**2 - 4*y**2 + 1\n",
      "\n",
      "dV/dt Simplified:\n",
      "alpha - 4*x**2 - 4*y**2 + 1\n",
      "\n",
      "Sign Changes: \n",
      "[(-sqrt(alpha - 4*y**2 + 1)/2, y), (sqrt(alpha - 4*y**2 + 1)/2, y)]\n"
     ]
    }
   ],
   "source": [
    "#Guess a psi for Poincare - Bendixson Theorem\n",
    "psi = 1\n",
    "def sign_change_test(xdot, ydot, psi):\n",
    "    div_psiX = sm.diff(psi*xdot,x) + sm.diff(psi*ydot, y)\n",
    "    print(f'dV/dt Expression: \\n{div_psiX}')\n",
    "    print(f'\\ndV/dt Simplified:\\n{sm.simplify(div_psiX)}\\n')\n",
    "    \n",
    "    sign_changes = sm.solve(sm.Eq(0,div_psiX),x,y)\n",
    "    print(f'Sign Changes: \\n{sign_changes}')\n",
    "    return div_psiX\n",
    "\n",
    "expr = sign_change_test(xdot, ydot, psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "063e52d9-875c-4dc9-83b6-6dca251ffa43",
   "metadata": {},
   "source": [
    "Now let me draw a rectangle around this limit cycle the same as the last problem. A rectangle at points (-1, -1), (-1, 1), (1, -1), (1, 1), and a smaller internal rectangle at (-0.1, -0.1), (-0.1, 0.1), (0.1, 0.1), (0.1, -0.1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a854e3d8-dee6-4084-be05-89e0686557b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 y^{2} - 3$"
      ],
      "text/plain": [
       "alpha - 4*y**2 - 3"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,-1)\n",
    "#This is clearly always < 0. No sign change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "488b307d-1f03-4f33-8a0c-76dffc194d2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 x^{2} - 3$"
      ],
      "text/plain": [
       "alpha - 4*x**2 - 3"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9b011142-9718-447e-ae43-9d30cff58568",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 y^{2} - 3$"
      ],
      "text/plain": [
       "alpha - 4*y**2 - 3"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "8b82bc08-3e2d-4db8-8424-0678fc762a68",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 x^{2} - 3$"
      ],
      "text/plain": [
       "alpha - 4*x**2 - 3"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee28ae38-576e-4bc5-8c53-45e8dea3de1a",
   "metadata": {},
   "source": [
    "For $0 < \\alpha < 1$ these are all clearly always negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "31752100-ae0a-4a84-8ba9-ef4d7147a5e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 y^{2} + 0.96$"
      ],
      "text/plain": [
       "alpha - 4*y**2 + 0.96"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7dc2aee7-f0f7-40d6-bad3-aee37f7b550a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 y^{2} + 0.96$"
      ],
      "text/plain": [
       "alpha - 4*y**2 + 0.96"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "3c9c2d41-f3b5-4026-a8e7-d782388d7686",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 x^{2} + 0.96$"
      ],
      "text/plain": [
       "alpha - 4*x**2 + 0.96"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "126023fa-4546-4fa4-a4a8-70de4ef70cfa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\alpha - 4 x^{2} + 0.96$"
      ],
      "text/plain": [
       "alpha - 4*x**2 + 0.96"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "519cbcda-bcf9-4e05-8d99-4356d16ffd5c",
   "metadata": {},
   "source": [
    "For x or y = 0.1, it's clear this is always positive for $ 0 < \\alpha < 1 $.\n",
    "As a result, there must be a limit cycle between these two rectangles that is stable, as flow enters the region on all boundaries."
   ]
  },
  {
   "attachments": {
    "a49d2138-fa02-44c3-bad9-c574795b9cc7.png": {
     "image/png": "[encoded data, 25624 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "0bc31580-8784-4e92-bbc8-93ee06cd9bf2",
   "metadata": {},
   "source": [
    "![image.png](attachment:a49d2138-fa02-44c3-bad9-c574795b9cc7.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44598ad7-a1a8-486a-8d27-839c657966c3",
   "metadata": {},
   "source": [
    "**Problem A**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "72a7604d-a938-45f6-9de4-4c2aa6acda05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 200712 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xdot = y\n",
    "ydot = -x - (1+x**2+x**4)*y\n",
    "\n",
    "#Create Functions\n",
    "P = lambda x, y: y\n",
    "Q = lambda x, y: -x - (1+x**2+x**4)*y\n",
    "def dX_dt(X, t, args = (P,Q)):\n",
    "    return [P(X[0], X[1]), Q(X[0], X[1])]\n",
    "\n",
    "#Initial Conditions\n",
    "ts = np.linspace(0, 5, 30)\n",
    "ic = np.linspace(-1.2, 1.2, 7)\n",
    "\n",
    "#Make a vectorfield\n",
    "X, Y = np.mgrid[ic[0]:ic[-1]:20j, ic[0]:ic[-1]:20j]\n",
    "u = P(X,Y)\n",
    "v = Q(X,Y)\n",
    "plt.quiver(X, Y, u, v, color = 'red')\n",
    "\n",
    "#plot trajectories\n",
    "plt.title(f'Phase Portrait for $\\\\dot x = $ {xdot} and $\\\\dot y = $ {ydot}')\n",
    "plt.xlabel('x'); plt.ylabel('y'); plt.grid('both')\n",
    "plt.xlim((ic[0],ic[-1])); plt.ylim((ic[0],ic[-1]))\n",
    "\n",
    "for r in ic:\n",
    "    for s in ic:\n",
    "        X0 = [r,s]\n",
    "        Xs = odeint(dX_dt,X0,ts)\n",
    "        plt.plot(Xs[:,0],Xs[:,1], \"royalblue\",linewidth = 0.75)"
   ]
  },
  {
   "attachments": {
    "7767ef13-07dc-41ce-a54e-beef95d98016.png": {
     "image/png": "[encoded data, 41956 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "fe2f9eb3-0261-4dd9-bd11-b9f235651b41",
   "metadata": {},
   "source": [
    "![image.png](attachment:7767ef13-07dc-41ce-a54e-beef95d98016.png)"
   ]
  },
  {
   "attachments": {
    "f74938ea-d577-425c-b7b1-70a42814e493.png": {
     "image/png": "[encoded data, 46540 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "80d5238b-e994-4493-82b9-66a90fe9251e",
   "metadata": {},
   "source": [
    "![image.png](attachment:f74938ea-d577-425c-b7b1-70a42814e493.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
