{
 "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 0x7fdf6bb15880>]"
      ]
     },
     "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": 21,
   "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": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(x,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y,-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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": 24,
     "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": 25,
   "id": "5b6610b7-3ae1-474e-bdcd-99f29fa603bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - x^{4} - x^{2} - 1$"
      ],
      "text/plain": [
       "-x**4 - x**2 - 1"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdot = y\n",
    "ydot = -x - (1+x**2+x**4)*y\n",
    "zeta = 1\n",
    "fun = sm.simplify(sm.diff(zeta*xdot,x) + sm.diff(zeta*ydot,y))\n",
    "fun"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce9a8a1a-4790-4825-8ddb-89145a4260f9",
   "metadata": {},
   "source": [
    "This function cannot change sign for all real numbers. Therefore, there is no limit cycle in the 2d plane.\n",
    "\n",
    "**Problem B**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "03ed9400-4f8a-4d93-95c4-8dc567a802dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 - x$"
      ],
      "text/plain": [
       "2 - x"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdot = x - x**2 + 2*y**2\n",
    "ydot = y*(x+1)\n",
    "\n",
    "zeta = 1\n",
    "\n",
    "fun = sm.simplify(sm.diff(zeta*xdot,x) + sm.diff(zeta*ydot,y))\n",
    "fun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ef02629e-77eb-4876-b1bf-a14c4167818e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 77916 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Create Functions\n",
    "P = lambda x, y: x - x**2 + 2*y**2\n",
    "Q = lambda x, y: y*(x+1)\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, 0.1, 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",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eca4d569-15b4-4167-a4ae-5360acd6e256",
   "metadata": {},
   "source": [
    "Now we know we can't just use Bendixson's criterion. But! The x axis looks invariant. Let's take a closer look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "98e0ce12-6e96-4b4c-825f-56d39aad8478",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dV/dt Expression: \n",
      "2 - x\n",
      "\n",
      "dV/dt Simplified:\n",
      "2 - x\n",
      "\n",
      "Sign Changes: \n",
      "[(2, y)]\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2$"
      ],
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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)\n",
    "\n",
    "expr.subs(x,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "826ecb2a-a9aa-4617-b1fc-9fb207b84622",
   "metadata": {},
   "source": [
    "So we know the x axis is invariant. Now to know if we have a limit cycle or not, we just need to be sure there are no critical points on the x axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "899989af-6470-4837-9674-2784f372be09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-1, -1), (-1, 1), (0, 0), (1, 0)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Peqn = sm.Eq(xdot, 0)\n",
    "Qeqn = sm.Eq(ydot, 0)\n",
    "\n",
    "criticalpoints = sm.solve((Peqn,Qeqn),x,y)\n",
    "display(criticalpoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fed0883c-3a2e-457b-8fa7-1a837bbafc2d",
   "metadata": {},
   "source": [
    "We have one troublesome critical point. Let's check if this is the right sign for us to say there's no limit cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "a27622ff-ab8e-4d8a-ab98-03631d28ff24",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 y^{2}$"
      ],
      "text/plain": [
       "2*y**2"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdot.subs(x, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "770cfd39-1620-4284-8115-202953c1c9b3",
   "metadata": {},
   "source": [
    "This is always positive, just like our other information about the x axis. Therefore there is no limit cycle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e023c7ab-7b27-4897-a378-5268d3d2c38d",
   "metadata": {},
   "source": [
    "**Problem C**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "e800d1c1-5823-47b7-8425-c73fc5a07bf5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 2 x^{2} - 6$"
      ],
      "text/plain": [
       "-2*x**2 - 6"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdot = y**2 - 2*x\n",
    "ydot = 3 - 4*y - 2*x**2*y\n",
    "\n",
    "zeta = 1 \n",
    "\n",
    "fun = sm.simplify(sm.diff(zeta*xdot,x) + sm.diff(zeta*ydot,y))\n",
    "fun"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92bc04df-6400-470f-9c08-f39b992ce675",
   "metadata": {},
   "source": [
    "This function does not change sign for all real numbers. Therefore there is no limit cycle thanks to Bendixson's Criterion."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77be06d4-3205-4c4f-b989-22ae7f906db4",
   "metadata": {},
   "source": [
    "**Problem D**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "b7cd9c6d-1552-4295-b68e-d8025d393e5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - x^{2} - 3$"
      ],
      "text/plain": [
       "-x**2 - 3"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xdot = -x + y**3 - y**4\n",
    "ydot = 1 - 2*y - x**2*y + x**4\n",
    "\n",
    "zeta = 1\n",
    "\n",
    "fun = sm.simplify(sm.diff(zeta*xdot,x) + sm.diff(zeta*ydot,y))\n",
    "fun"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c0bb802-bc3d-475d-ba8e-15576cda7926",
   "metadata": {},
   "source": [
    "This function also does not change sign for all real numbers. There is no limit cycle thanks to Bendixson's Criterion."
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "658651a7-69e8-4852-bb1f-8f767c3715b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\epsilon^{3} \\left(\\omega^{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + x_{1}^{2}{\\left(t \\right)}\\right) + \\epsilon^{2} \\left(\\omega^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + 2 \\omega \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)}\\right) + \\epsilon \\left(2 \\omega \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} + x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)}\\right) + x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "epsilon**3*(omega**2*Derivative(x_1(t), (t, 2)) + x_1(t)**2) + epsilon**2*(omega**2*Derivative(x_0(t), (t, 2)) + 2*omega*Derivative(x_1(t), (t, 2)) + 2*x_0(t)*x_1(t)) + epsilon*(2*omega*Derivative(x_0(t), (t, 2)) + x_0(t)**2 + x_1(t) + Derivative(x_1(t), (t, 2))) + x_0(t) + Derivative(x_0(t), (t, 2))"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sm.Function('x')\n",
    "x0 = sm.Function('x_0')\n",
    "x1 = sm.Function('x_1')\n",
    "\n",
    "t = sm.Symbol('t')\n",
    "eps = sm.Symbol('epsilon')\n",
    "omega = sm.Symbol('omega')\n",
    "\n",
    "#Set up series expansion\n",
    "x = x0(t) + eps*x1(t)\n",
    "dt = (1+eps*omega)\n",
    "\n",
    "expr = dt**2*x.diff(t,t) + x + eps*x**2 \n",
    "sm.collect(sm.expand(expr),eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "965b5c7f-cd15-4d86-9b17-cecbfa77ec0a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{0}{\\left(t \\right)} = C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)}$"
      ],
      "text/plain": [
       "Eq(x_0(t), C1*sin(t) + C2*cos(t))"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solve for O(1)\n",
    "O_1 = sm.dsolve(x0(t) + x0(t).diff(t,t))\n",
    "O_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "b507187a-0c37-4f67-a10c-2237fa3538fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{1}{\\left(t \\right)} = \\frac{- 2 C_{1} \\omega \\left(\\epsilon \\omega + 1\\right)^{2} \\sin{\\left(t \\right)} - 2 C_{2} \\omega \\left(\\epsilon \\omega + 1\\right)^{2} \\cos{\\left(t \\right)} - \\left(C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)}\\right)^{2}}{\\left(\\epsilon \\omega + 1\\right)^{2} + 1}$"
      ],
      "text/plain": [
       "Eq(x_1(t), (-2*C1*omega*(epsilon*omega + 1)**2*sin(t) - 2*C2*omega*(epsilon*omega + 1)**2*cos(t) - (C1*sin(t) + C2*cos(t))**2)/((epsilon*omega + 1)**2 + 1))"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solve for O(eps)\n",
    "x0 = O_1.rhs #extract x0 solution\n",
    "\n",
    "#Solve for x1\n",
    "O_eps = sm.dsolve(sm.simplify(2*omega*dt**2*x0 + x0**2 + x1(t) + dt**2*x1(t)), x1(t))\n",
    "O_eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "102fc8da-f0c6-43a9-99c2-64580071b3ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{1}{\\left(t \\right)} = - \\frac{\\left(C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)}\\right)^{2}}{2}$"
      ],
      "text/plain": [
       "Eq(x_1(t), -(C1*sin(t) + C2*cos(t))**2/2)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Now a hiccup. We need to make sure we have no secular terms. By inspection, we can do this by setting $\\omega$ = 0.\n",
    "O_eps = O_eps.subs(omega, 0)\n",
    "O_eps.simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2be9ea31-fb5a-49cc-888f-45b79d8fa22e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0 = C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)} - \\frac{\\epsilon \\left(C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)}\\right)^{2}}{2}$"
      ],
      "text/plain": [
       "Eq(0, C1*sin(t) + C2*cos(t) - epsilon*(C1*sin(t) + C2*cos(t))**2/2)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Now plug everything back into our expression for x:\n",
    "xsoln = sm.Eq(0, x0 + eps*O_eps.rhs)\n",
    "xsoln"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "e6c41e17-6b70-425a-8b9d-10c83f4b3e94",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\epsilon^{9} \\omega_{2} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\epsilon^{8} \\left(\\omega_{1} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{2} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\omega_{2} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)}\\right) + \\epsilon^{7} \\left(2 \\omega_{1} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\omega_{1} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\omega_{2} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{2} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\omega_{2} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)}\\right) + \\epsilon^{6} \\left(2 \\omega_{1} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\omega_{1} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{1} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\omega_{1} x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\omega_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\omega_{2} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{2} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)}\\right) + \\epsilon^{5} \\left(2 \\omega_{1} \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + 2 \\omega_{1} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{1} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\omega_{1} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{1} x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\omega_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + \\omega_{2} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\omega_{2} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{2} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + x_{2}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)}\\right) + \\epsilon^{4} \\left(\\omega_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + 2 \\omega_{1} \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + \\omega_{1} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 \\omega_{1} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{1} x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\omega_{1} x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{1} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\omega_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + \\omega_{2} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{2} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{2} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 x_{1}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)}\\right) + \\epsilon^{3} \\left(\\omega_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 \\omega_{1} \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + \\omega_{1} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{1} x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{1} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 \\omega_{1} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + \\omega_{2} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{2} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{1}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\frac{d}{d t} x_{2}{\\left(t \\right)}\\right) + \\epsilon^{2} \\left(\\omega_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + \\omega_{1} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{1} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 \\omega_{1} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}{\\left(t \\right)} - \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)}\\right) + \\epsilon \\left(2 \\omega_{1} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{1}{\\left(t \\right)} - \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)}\\right) + x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "epsilon**9*omega_2*x_2(t)**2*Derivative(x_2(t), t) + epsilon**8*(omega_1*x_2(t)**2*Derivative(x_2(t), t) + 2*omega_2*x_1(t)*x_2(t)*Derivative(x_2(t), t) + omega_2*x_2(t)**2*Derivative(x_1(t), t)) + epsilon**7*(2*omega_1*x_1(t)*x_2(t)*Derivative(x_2(t), t) + omega_1*x_2(t)**2*Derivative(x_1(t), t) + 2*omega_2*x_0(t)*x_2(t)*Derivative(x_2(t), t) + omega_2*x_1(t)**2*Derivative(x_2(t), t) + 2*omega_2*x_1(t)*x_2(t)*Derivative(x_1(t), t) + omega_2*x_2(t)**2*Derivative(x_0(t), t) + x_2(t)**2*Derivative(x_2(t), t)) + epsilon**6*(2*omega_1*x_0(t)*x_2(t)*Derivative(x_2(t), t) + omega_1*x_1(t)**2*Derivative(x_2(t), t) + 2*omega_1*x_1(t)*x_2(t)*Derivative(x_1(t), t) + omega_1*x_2(t)**2*Derivative(x_0(t), t) + omega_2**2*Derivative(x_2(t), (t, 2)) + 2*omega_2*x_0(t)*x_1(t)*Derivative(x_2(t), t) + 2*omega_2*x_0(t)*x_2(t)*Derivative(x_1(t), t) + omega_2*x_1(t)**2*Derivative(x_1(t), t) + 2*omega_2*x_1(t)*x_2(t)*Derivative(x_0(t), t) + 2*x_1(t)*x_2(t)*Derivative(x_2(t), t) + x_2(t)**2*Derivative(x_1(t), t)) + epsilon**5*(2*omega_1*omega_2*Derivative(x_2(t), (t, 2)) + 2*omega_1*x_0(t)*x_1(t)*Derivative(x_2(t), t) + 2*omega_1*x_0(t)*x_2(t)*Derivative(x_1(t), t) + omega_1*x_1(t)**2*Derivative(x_1(t), t) + 2*omega_1*x_1(t)*x_2(t)*Derivative(x_0(t), t) + omega_2**2*Derivative(x_1(t), (t, 2)) + omega_2*x_0(t)**2*Derivative(x_2(t), t) + 2*omega_2*x_0(t)*x_1(t)*Derivative(x_1(t), t) + 2*omega_2*x_0(t)*x_2(t)*Derivative(x_0(t), t) + omega_2*x_1(t)**2*Derivative(x_0(t), t) - omega_2*Derivative(x_2(t), t) + 2*x_0(t)*x_2(t)*Derivative(x_2(t), t) + x_1(t)**2*Derivative(x_2(t), t) + 2*x_1(t)*x_2(t)*Derivative(x_1(t), t) + x_2(t)**2*Derivative(x_0(t), t)) + epsilon**4*(omega_1**2*Derivative(x_2(t), (t, 2)) + 2*omega_1*omega_2*Derivative(x_1(t), (t, 2)) + omega_1*x_0(t)**2*Derivative(x_2(t), t) + 2*omega_1*x_0(t)*x_1(t)*Derivative(x_1(t), t) + 2*omega_1*x_0(t)*x_2(t)*Derivative(x_0(t), t) + omega_1*x_1(t)**2*Derivative(x_0(t), t) - omega_1*Derivative(x_2(t), t) + omega_2**2*Derivative(x_0(t), (t, 2)) + omega_2*x_0(t)**2*Derivative(x_1(t), t) + 2*omega_2*x_0(t)*x_1(t)*Derivative(x_0(t), t) - omega_2*Derivative(x_1(t), t) + 2*omega_2*Derivative(x_2(t), (t, 2)) + 2*x_0(t)*x_1(t)*Derivative(x_2(t), t) + 2*x_0(t)*x_2(t)*Derivative(x_1(t), t) + x_1(t)**2*Derivative(x_1(t), t) + 2*x_1(t)*x_2(t)*Derivative(x_0(t), t)) + epsilon**3*(omega_1**2*Derivative(x_1(t), (t, 2)) + 2*omega_1*omega_2*Derivative(x_0(t), (t, 2)) + omega_1*x_0(t)**2*Derivative(x_1(t), t) + 2*omega_1*x_0(t)*x_1(t)*Derivative(x_0(t), t) - omega_1*Derivative(x_1(t), t) + 2*omega_1*Derivative(x_2(t), (t, 2)) + omega_2*x_0(t)**2*Derivative(x_0(t), t) - omega_2*Derivative(x_0(t), t) + 2*omega_2*Derivative(x_1(t), (t, 2)) + x_0(t)**2*Derivative(x_2(t), t) + 2*x_0(t)*x_1(t)*Derivative(x_1(t), t) + 2*x_0(t)*x_2(t)*Derivative(x_0(t), t) + x_1(t)**2*Derivative(x_0(t), t) - Derivative(x_2(t), t)) + epsilon**2*(omega_1**2*Derivative(x_0(t), (t, 2)) + omega_1*x_0(t)**2*Derivative(x_0(t), t) - omega_1*Derivative(x_0(t), t) + 2*omega_1*Derivative(x_1(t), (t, 2)) + 2*omega_2*Derivative(x_0(t), (t, 2)) + x_0(t)**2*Derivative(x_1(t), t) + 2*x_0(t)*x_1(t)*Derivative(x_0(t), t) + x_2(t) - Derivative(x_1(t), t) + Derivative(x_2(t), (t, 2))) + epsilon*(2*omega_1*Derivative(x_0(t), (t, 2)) + x_0(t)**2*Derivative(x_0(t), t) + x_1(t) - Derivative(x_0(t), t) + Derivative(x_1(t), (t, 2))) + x_0(t) + Derivative(x_0(t), (t, 2))"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sm.Function('x')\n",
    "x0 = sm.Function('x_0')\n",
    "x1 = sm.Function('x_1')\n",
    "x2 = sm.Function('x_2')\n",
    "\n",
    "t = sm.Symbol('t')\n",
    "eps = sm.Symbol('epsilon')\n",
    "omega_1 = sm.Symbol('omega_1')\n",
    "omega_2 = sm.Symbol('omega_2')\n",
    "\n",
    "#Set up series expansion\n",
    "x_expr = x0(t) + eps*x1(t) + eps**2*x2(t)\n",
    "dt = (1+eps*omega_1+eps**2*omega_2)\n",
    "\n",
    "expr = dt**2*x_expr.diff(t,t) + eps*(x_expr**2 - 1) * dt * x_expr.diff(t) + x_expr \n",
    "sm.collect(sm.expand(expr),eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "c70d7149-3628-4a60-b261-5eedc052e72b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Order 0: \n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "x_0(t) + Derivative(x_0(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Order 1: \n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 \\omega_{1} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{1}{\\left(t \\right)} - \\frac{d}{d t} x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "2*omega_1*Derivative(x_0(t), (t, 2)) + x_0(t)**2*Derivative(x_0(t), t) + x_1(t) - Derivative(x_0(t), t) + Derivative(x_1(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Order 2: \n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\omega_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + \\omega_{1} x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} - \\omega_{1} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 \\omega_{1} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 \\omega_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}^{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + 2 x_{0}{\\left(t \\right)} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}{\\left(t \\right)} - \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "omega_1**2*Derivative(x_0(t), (t, 2)) + omega_1*x_0(t)**2*Derivative(x_0(t), t) - omega_1*Derivative(x_0(t), t) + 2*omega_1*Derivative(x_1(t), (t, 2)) + 2*omega_2*Derivative(x_0(t), (t, 2)) + x_0(t)**2*Derivative(x_1(t), t) + 2*x_0(t)*x_1(t)*Derivative(x_0(t), t) + x_2(t) - Derivative(x_1(t), t) + Derivative(x_2(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#getting the first 3 orders:\n",
    "eps_ords = sm.collect(sm.expand(expr),eps,evaluate=False)\n",
    "eps_list= []\n",
    "for order in range(3):\n",
    "    eps_list.append(eps_ords[eps**order])\n",
    "    print(f'\\nOrder {order}: ')\n",
    "    display(eps_list[order])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "d071eaab-7b56-4a69-bc0d-9d231fb780cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{0}{\\left(t \\right)} = C_{1} \\sin{\\left(t \\right)} + C_{2} \\cos{\\left(t \\right)}$"
      ],
      "text/plain": [
       "Eq(x_0(t), C1*sin(t) + C2*cos(t))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Solve 0th order\n",
    "ord_0 = sm.dsolve(eps_list[0],x0(t))\n",
    "display(ord_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "7d7fd828-c8d1-4773-9f43-8f4cb9ba51f4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{1}{\\left(t \\right)} = - \\frac{C_{1} \\left(C_{1}^{2} - 3 C_{2}^{2}\\right) \\cos^{3}{\\left(t \\right)}}{8} + \\frac{\\left(3 C_{1}^{3} - C_{1}^{2} C_{2} t - 9 C_{1} C_{2}^{2} - 8 C_{1} \\omega_{1} t - C_{2}^{3} t + 4 C_{2} t + 8 C_{3}\\right) \\cos{\\left(t \\right)}}{8} + \\frac{\\left(- C_{1}^{3} t + 4 C_{1}^{2} C_{2} - C_{1} C_{2}^{2} t + 8 C_{1} \\omega_{1} + 4 C_{1} t + 8 C_{2} \\omega_{1} t + C_{2} \\left(3 C_{1}^{2} - C_{2}^{2}\\right) \\cos^{2}{\\left(t \\right)} - 4 C_{2} + 8 C_{4}\\right) \\sin{\\left(t \\right)}}{8}$"
      ],
      "text/plain": [
       "Eq(x_1(t), -C1*(C1**2 - 3*C2**2)*cos(t)**3/8 + (3*C1**3 - C1**2*C2*t - 9*C1*C2**2 - 8*C1*omega_1*t - C2**3*t + 4*C2*t + 8*C3)*cos(t)/8 + (-C1**3*t + 4*C1**2*C2 - C1*C2**2*t + 8*C1*omega_1 + 4*C1*t + 8*C2*omega_1*t + C2*(3*C1**2 - C2**2)*cos(t)**2 - 4*C2 + 8*C4)*sin(t)/8)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Solve 1st Order\n",
    "ord_1 = sm.dsolve(eps_list[1].subs(x0(t),ord_0.rhs),x1(t))\n",
    "display(ord_1.simplify())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "84909933-3519-43c8-90a2-94dc13511795",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-C1**3*sin(t) - C1**2*C2*cos(t) - C1*C2**2*sin(t) + 4*C1*sin(t) - C2**3*cos(t) + 4*C2*cos(t))/(8*(C1*cos(t) - C2*sin(t)))]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Collect t terms to solve for omega\n",
    "ord_1_terms = sm.collect(sm.expand(ord_1.rhs), t, evaluate=False)\n",
    "\n",
    "#Solve for omega that makes that t term equal to 0.\n",
    "omega_1_soln = sm.solve(sm.Eq(0,ord_1_terms[t]), omega_1)\n",
    "omega_1_soln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "5f46858f-7448-47d4-8f65-588728554a40",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eq(x_1(t), (-C1**4*cos(2*t)**2 + 6*C1**4*cos(2*t) + 3*C1**4 + 4*C1**3*C2*sin(2*t) + 2*C1**3*C2*sin(4*t) + 6*C1**2*C2**2*cos(2*t)**2 - 2*C1**2*C2**2*cos(2*t) - 28*C1**2*C2**2 - 8*C1**2*cos(2*t) + 8*C1**2 + 12*C1*C2**3*sin(2*t) - 2*C1*C2**3*sin(4*t) + 16*C1*C3*cos(2*t) + 16*C1*C3 + 16*C1*C4*sin(2*t) - C2**4*cos(2*t)**2 + C2**4 - 8*C2**2*cos(2*t) + 8*C2**2 - 16*C2*C3*sin(2*t) + 16*C2*C4*cos(2*t) - 16*C2*C4)/(32*C1*cos(t) - 32*C2*sin(t)))\n"
     ]
    }
   ],
   "source": [
    "#Plug this solved omega_1 back into our expression for the order 1 terms\n",
    "ord_1 = sm.simplify(ord_1.subs(omega_1, list(omega_1_soln)[0]))\n",
    "ord_1\n",
    "print(ord_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "4eb17f76-b2d2-4a56-a6df-78555443aa41",
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "closing parenthesis ')' does not match opening parenthesis '{' (3522621134.py, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[44], line 2\u001b[0;36m\u001b[0m\n\u001b[0;31m    ord_2 = sm.simplify(sm.dsolve(eps_list[2].subs({x1(t): ord_1.rhs), x0(t): ord_0.rhs},x2(t)))\u001b[0m\n\u001b[0m                                                                    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m closing parenthesis ')' does not match opening parenthesis '{'\n"
     ]
    }
   ],
   "source": [
    "#Again for order 2\n",
    "ord_2 = sm.simplify(sm.dsolve(eps_list[2].subs({x1(t): ord_1.rhs), x0(t): ord_0.rhs},x2(t)))\n",
    "ord_2\n",
    "print(ord_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64f88151-64b0-4a6f-90df-6d33be83aa65",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Collect t terms to solve for omega_2\n",
    "ord_2_temrs = sm.collect(sm.expand(ord_w.rhs), t, evaluate=False)\n",
    "\n",
    "#Solve for omega_2 that makes that t term equal to 0.\n",
    "omega_2_soln = sm.solve(sm.Eq(0,ord_2_terms[t]), omega_2)\n",
    "omega_2_soln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b215845a-9b51-4aea-b947-c6d530a05a27",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plug this solved omega_2 back into the expression for the order 2 terms:\n",
    "ord_2 = sm.simplify(ord_2.subs(omega_2, list(omega_2_soln)[0]))\n",
    "ord_2\n",
    "print(ord_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32cbfb5c-7e61-431c-b7e7-e06c2ee05bb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_soln = sm.Eq(x(t), ord_0.rhs + eps*ord_1.rhs + eps**2*ord_2.rhs"
   ]
  }
 ],
 "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
}
