{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6945cd62",
   "metadata": {},
   "source": [
    "## Example Problem - Linstedt Poincare Procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "797e4d60",
   "metadata": {},
   "source": [
    "In this problem, we will apply Lindstedt's method to $\\ddot{x} + \\epsilon x \\dot{x} + x = 0$, for $0 < \\epsilon << 1$, where $x(0) = a_0$, and $\\dot{x}(0) = 0$. We will show that the frequency-amplitude relation for periodic solutions is given by $\\omega = 1 - \\frac{1}{24}a_0^2 \\epsilon^2 + O(\\epsilon^3)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f25ba29d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7f1e21b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Function('x')\n",
    "x0 = Function('x_0')\n",
    "x1 = Function('x_1')\n",
    "x2 = Function('x_2')\n",
    "t = Symbol('t')\n",
    "eps = Symbol('epsilon')\n",
    "w1 = Symbol('w_1')\n",
    "w2 = Symbol('w_2')\n",
    "w = Symbol('w')\n",
    "a0 = Symbol('a_0')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6364f15",
   "metadata": {},
   "source": [
    "Note that to a second order expansion in the frequency due to Lindstedt method, $\\tau = \\omega t$, with $\\omega = 1 + \\epsilon \\omega_1  + \\epsilon^2 \\omega_2$. In addtion, we use the expansion in the solution, $x = x_0(t) + \\epsilon  x_1(t) + \\epsilon^2  x_2(t)$, and substitute both into the differential equation we are trying to solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "989d7206",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\epsilon^{2} x_{2}{\\left(t \\right)} + \\epsilon \\left(\\epsilon^{2} w_{2} + \\epsilon w_{1} + 1\\right) \\left(\\epsilon^{2} x_{2}{\\left(t \\right)} + \\epsilon x_{1}{\\left(t \\right)} + x_{0}{\\left(t \\right)}\\right) \\left(\\epsilon^{2} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\epsilon \\frac{d}{d t} x_{1}{\\left(t \\right)} + \\frac{d}{d t} x_{0}{\\left(t \\right)}\\right) + \\epsilon x_{1}{\\left(t \\right)} + \\left(\\epsilon^{2} w_{2} + \\epsilon w_{1} + 1\\right)^{2} \\left(\\epsilon^{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + \\epsilon \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)}\\right) + x_{0}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "epsilon**2*x_2(t) + epsilon*(epsilon**2*w_2 + epsilon*w_1 + 1)*(epsilon**2*x_2(t) + epsilon*x_1(t) + x_0(t))*(epsilon**2*Derivative(x_2(t), t) + epsilon*Derivative(x_1(t), t) + Derivative(x_0(t), t)) + epsilon*x_1(t) + (epsilon**2*w_2 + epsilon*w_1 + 1)**2*(epsilon**2*Derivative(x_2(t), (t, 2)) + epsilon*Derivative(x_1(t), (t, 2)) + Derivative(x_0(t), (t, 2))) + x_0(t)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = 1 + eps * w1 + eps**2 * w2\n",
    "x = x0(t) + eps * x1(t) + eps**2 * x2(t)\n",
    "expr = w**2 * x.diff(t,t) + w*eps*x*x.diff(t) + x\n",
    "display(expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71535ff6",
   "metadata": {},
   "source": [
    "Next, we expand and collect like terms of $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a1762f64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\epsilon^{7} w_{2} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + \\epsilon^{6} \\left(w_{1} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + w_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + w_{2} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + w_{2} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)}\\right) + \\epsilon^{5} \\cdot \\left(2 w_{1} w_{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + w_{1} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + w_{1} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + w_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + w_{2} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + w_{2} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + w_{2} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)}\\right) + \\epsilon^{4} \\left(w_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + 2 w_{1} w_{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + w_{1} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + w_{1} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + w_{1} x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + w_{2}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + w_{2} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + w_{2} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 w_{2} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)}\\right) + \\epsilon^{3} \\left(w_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 w_{1} w_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + w_{1} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + w_{1} x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 w_{1} \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} + w_{2} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 w_{2} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{2}{\\left(t \\right)} + x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + x_{2}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)}\\right) + \\epsilon^{2} \\left(w_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + w_{1} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 w_{1} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 w_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)}\\right) + \\epsilon \\left(2 w_{1} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\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**7*w_2*x_2(t)*Derivative(x_2(t), t) + epsilon**6*(w_1*x_2(t)*Derivative(x_2(t), t) + w_2**2*Derivative(x_2(t), (t, 2)) + w_2*x_1(t)*Derivative(x_2(t), t) + w_2*x_2(t)*Derivative(x_1(t), t)) + epsilon**5*(2*w_1*w_2*Derivative(x_2(t), (t, 2)) + w_1*x_1(t)*Derivative(x_2(t), t) + w_1*x_2(t)*Derivative(x_1(t), t) + w_2**2*Derivative(x_1(t), (t, 2)) + w_2*x_0(t)*Derivative(x_2(t), t) + w_2*x_1(t)*Derivative(x_1(t), t) + w_2*x_2(t)*Derivative(x_0(t), t) + x_2(t)*Derivative(x_2(t), t)) + epsilon**4*(w_1**2*Derivative(x_2(t), (t, 2)) + 2*w_1*w_2*Derivative(x_1(t), (t, 2)) + w_1*x_0(t)*Derivative(x_2(t), t) + w_1*x_1(t)*Derivative(x_1(t), t) + w_1*x_2(t)*Derivative(x_0(t), t) + w_2**2*Derivative(x_0(t), (t, 2)) + w_2*x_0(t)*Derivative(x_1(t), t) + w_2*x_1(t)*Derivative(x_0(t), t) + 2*w_2*Derivative(x_2(t), (t, 2)) + x_1(t)*Derivative(x_2(t), t) + x_2(t)*Derivative(x_1(t), t)) + epsilon**3*(w_1**2*Derivative(x_1(t), (t, 2)) + 2*w_1*w_2*Derivative(x_0(t), (t, 2)) + w_1*x_0(t)*Derivative(x_1(t), t) + w_1*x_1(t)*Derivative(x_0(t), t) + 2*w_1*Derivative(x_2(t), (t, 2)) + w_2*x_0(t)*Derivative(x_0(t), t) + 2*w_2*Derivative(x_1(t), (t, 2)) + x_0(t)*Derivative(x_2(t), t) + x_1(t)*Derivative(x_1(t), t) + x_2(t)*Derivative(x_0(t), t)) + epsilon**2*(w_1**2*Derivative(x_0(t), (t, 2)) + w_1*x_0(t)*Derivative(x_0(t), t) + 2*w_1*Derivative(x_1(t), (t, 2)) + 2*w_2*Derivative(x_0(t), (t, 2)) + x_0(t)*Derivative(x_1(t), t) + x_1(t)*Derivative(x_0(t), t) + x_2(t) + Derivative(x_2(t), (t, 2))) + epsilon*(2*w_1*Derivative(x_0(t), (t, 2)) + x_0(t)*Derivative(x_0(t), t) + x_1(t) + Derivative(x_1(t), (t, 2))) + x_0(t) + Derivative(x_0(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "expr = expand(expr)\n",
    "expr = collect(expr,eps)\n",
    "display(expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7961634",
   "metadata": {},
   "source": [
    "We can clean this up a bit more. We can create a list of the equations! Let's do this for terms 0, 1, 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1e1f6897",
   "metadata": {},
   "outputs": [],
   "source": [
    "eqLHS = collect(expr,eps,evaluate=False)\n",
    "LHSlist = []\n",
    "for k in [0,1,2]:\n",
    "    LHSlist.append(eqLHS[eps**k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a63234af",
   "metadata": {},
   "outputs": [
    {
     "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"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 w_{1} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "2*w_1*Derivative(x_0(t), (t, 2)) + x_0(t)*Derivative(x_0(t), t) + x_1(t) + Derivative(x_1(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle w_{1}^{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + w_{1} x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + 2 w_{1} \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} + 2 w_{2} \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} + x_{0}{\\left(t \\right)} \\frac{d}{d t} x_{1}{\\left(t \\right)} + x_{1}{\\left(t \\right)} \\frac{d}{d t} x_{0}{\\left(t \\right)} + x_{2}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "w_1**2*Derivative(x_0(t), (t, 2)) + w_1*x_0(t)*Derivative(x_0(t), t) + 2*w_1*Derivative(x_1(t), (t, 2)) + 2*w_2*Derivative(x_0(t), (t, 2)) + x_0(t)*Derivative(x_1(t), t) + x_1(t)*Derivative(x_0(t), t) + x_2(t) + Derivative(x_2(t), (t, 2))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for k in [0,1,2]:\n",
    "    display(LHSlist[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f53aed26",
   "metadata": {},
   "source": [
    "We know that each of these is equal to zero! Let's solve the first one (order $O(\\epsilon^0))$. First we define it as an equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b295ad34",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{0}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{0}{\\left(t \\right)} = 0$"
      ],
      "text/plain": [
       "Eq(x_0(t) + Derivative(x_0(t), (t, 2)), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ODE class  nth_linear_constant_coeff_homogeneous\n"
     ]
    }
   ],
   "source": [
    "eq0 = Eq(LHSlist[0],0)\n",
    "display(eq0)\n",
    "print('ODE class ', classify_ode(eq0)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eac29466",
   "metadata": {},
   "source": [
    "Now we solve, including initial conditions. Note that there is some funny syntax that has to be used to specify the initial condition of the derivative. This is a quirk of SymPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f2ee72c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a_{0} \\cos{\\left(t \\right)}$"
      ],
      "text/plain": [
       "a_0*cos(t)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0sol = dsolve(eq0, ics={x0(0): a0, x0(t).diff(t).subs(t, 0): 0}).rhs\n",
    "display(x0sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77e50b75",
   "metadata": {},
   "source": [
    "Great, so now we have a solution for $x_0$. Let's go look at our order $O(\\epsilon^1)$ equation and substitute in the order $O(\\epsilon^0)$ solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "984909c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - a_{0}^{2} \\sin{\\left(t \\right)} \\cos{\\left(t \\right)} - 2 a_{0} w_{1} \\cos{\\left(t \\right)} + x_{1}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{1}{\\left(t \\right)} = 0$"
      ],
      "text/plain": [
       "Eq(-a_0**2*sin(t)*cos(t) - 2*a_0*w_1*cos(t) + x_1(t) + Derivative(x_1(t), (t, 2)), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ODE class  nth_linear_constant_coeff_variation_of_parameters\n"
     ]
    }
   ],
   "source": [
    "eq1 = Eq(simplify(LHSlist[1].subs(x0(t),x0sol)),0)\n",
    "display(eq1)\n",
    "print('ODE class ', classify_ode(eq1)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7d280ace",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{a_{0}^{2} \\sin{\\left(2 t \\right)}}{6} + \\left(\\frac{a_{0}^{2}}{3} + a_{0} t w_{1}\\right) \\sin{\\left(t \\right)}$"
      ],
      "text/plain": [
       "-a_0**2*sin(2*t)/6 + (a_0**2/3 + a_0*t*w_1)*sin(t)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1sol = dsolve(eq1, ics={x1(0): 0, x1(t).diff(t).subs(t, 0): 0}).rhs\n",
    "display(x1sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73f19ac6",
   "metadata": {},
   "source": [
    "Periodicity requires that $w_1 = 0$ (to remove the secular term)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4d3609b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{a_{0}^{2} \\sin{\\left(t \\right)}}{3} - \\frac{a_{0}^{2} \\sin{\\left(2 t \\right)}}{6}$"
      ],
      "text/plain": [
       "a_0**2*sin(t)/3 - a_0**2*sin(2*t)/6"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1sol_nonsec = x1sol.subs(w1, 0)\n",
    "display(x1sol_nonsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3de59464",
   "metadata": {},
   "source": [
    "Great, so now we have a solution for $x_0$ and $x_1$. Let's go look at our order $O(\\epsilon^2)$ equation and substitute in the order $O(\\epsilon^1)$ and $O(\\epsilon^0)$ solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a183bb78",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{a_{0}^{3} \\cdot \\left(2 \\sin{\\left(t \\right)} - \\sin{\\left(2 t \\right)}\\right) \\sin{\\left(t \\right)}}{6} + \\frac{a_{0}^{3} \\cdot \\left(2 \\cos{\\left(t \\right)} - 2 \\cos{\\left(2 t \\right)}\\right) \\cos{\\left(t \\right)}}{6} - 2 a_{0} w_{2} \\cos{\\left(t \\right)} + x_{2}{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} x_{2}{\\left(t \\right)} = 0$"
      ],
      "text/plain": [
       "Eq(-a_0**3*(2*sin(t) - sin(2*t))*sin(t)/6 + a_0**3*(2*cos(t) - 2*cos(2*t))*cos(t)/6 - 2*a_0*w_2*cos(t) + x_2(t) + Derivative(x_2(t), (t, 2)), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ODE class  nth_linear_constant_coeff_variation_of_parameters\n"
     ]
    }
   ],
   "source": [
    "eq2 = Eq(simplify(LHSlist[2].subs({x0(t):x0sol, x1(t):x1sol_nonsec, w1:0})),0)\n",
    "display(eq2)\n",
    "print('ODE class ', classify_ode(eq1)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d42ed9bd",
   "metadata": {},
   "source": [
    "Now, we solve it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f210aadc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{a_{0}^{3} \\left(- 9 \\cos{\\left(t \\right)} - 8\\right) \\sin^{4}{\\left(t \\right)}}{36} + \\frac{a_{0}^{3} \\cdot \\left(45 \\cos{\\left(t \\right)} - 4\\right) \\sin^{2}{\\left(t \\right)}}{72} + \\frac{a_{0}^{3} \\cos^{5}{\\left(t \\right)}}{4} - \\frac{13 a_{0}^{3} \\cos{\\left(t \\right)}}{36} + \\frac{a_{0}^{3}}{9} + \\left(\\frac{a_{0}^{3} t}{24} - \\frac{a_{0}^{3} \\sin{\\left(3 t \\right)}}{18} + a_{0} t w_{2}\\right) \\sin{\\left(t \\right)}$"
      ],
      "text/plain": [
       "a_0**3*(-9*cos(t) - 8)*sin(t)**4/36 + a_0**3*(45*cos(t) - 4)*sin(t)**2/72 + a_0**3*cos(t)**5/4 - 13*a_0**3*cos(t)/36 + a_0**3/9 + (a_0**3*t/24 - a_0**3*sin(3*t)/18 + a_0*t*w_2)*sin(t)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x2sol = dsolve(eq2, ics={x2(0): 0, x2(t).diff(t).subs(t, 0): 0}).rhs\n",
    "display(x2sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc0a815e",
   "metadata": {},
   "source": [
    "Lots of terms. But! We need to figure out what will give us secular behavior. So we can just collect with respect to $t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cbcb696c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{t: a_0**3*sin(t)/24 + a_0*w_2*sin(t),\n",
       " 1: -a_0**3*sin(t)**4*cos(t)/4 - 2*a_0**3*sin(t)**4/9 + 5*a_0**3*sin(t)**2*cos(t)/8 - a_0**3*sin(t)**2/18 - a_0**3*sin(t)*sin(3*t)/18 + a_0**3*cos(t)**5/4 - 13*a_0**3*cos(t)/36 + a_0**3/9}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x2solterms = collect(expand(x2sol),t,evaluate=False)\n",
    "display(x2solterms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "14144cab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{a_{0}^{3} \\sin{\\left(t \\right)}}{24} + a_{0} w_{2} \\sin{\\left(t \\right)} = 0$"
      ],
      "text/plain": [
       "Eq(a_0**3*sin(t)/24 + a_0*w_2*sin(t), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[-a_0**2/24]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "testeqn=Eq(x2solterms[t],0)\n",
    "display(testeqn)\n",
    "w2sol = solve(Eq(x2solterms[t],0),w2)\n",
    "display(w2sol)\n",
    "w2sol = list(w2sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d2344049",
   "metadata": {},
   "outputs": [],
   "source": [
    "w2val = w2sol[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3087dbb",
   "metadata": {},
   "source": [
    "Therefore, we can get the expression for $w$ that..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a9f106f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{a_{0}^{2} \\epsilon^{2}}{24} + 1$"
      ],
      "text/plain": [
       "-a_0**2*epsilon**2/24 + 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- \\frac{a_{0}^{2} \\epsilon^{2}}{24} + 1\n"
     ]
    }
   ],
   "source": [
    "wsol = Symbol('w_sol')\n",
    "wsol = w.subs(w2,w2val).subs(w1,0)\n",
    "display(wsol)\n",
    "print(latex(wsol))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fc92804",
   "metadata": {},
   "source": [
    "Therefore, our answer is that $\\omega =  - \\frac{a_{0}^{2} \\epsilon^{2}}{24} + 1 + O(\\epsilon^3)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b038975f",
   "metadata": {},
   "source": [
    "Okay, but what if we want to plot our solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "256fe342",
   "metadata": {},
   "outputs": [],
   "source": [
    "xsol = x.subs(x0(t),x0sol).subs(x1(t),x1sol).subs(w1,0).subs(x2(t),x2sol).subs(w2,w2val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9767609f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a_{0} \\cos{\\left(t \\right)} + \\epsilon^{2} \\left(\\frac{a_{0}^{3} \\left(- 9 \\cos{\\left(t \\right)} - 8\\right) \\sin^{4}{\\left(t \\right)}}{36} + \\frac{a_{0}^{3} \\cdot \\left(45 \\cos{\\left(t \\right)} - 4\\right) \\sin^{2}{\\left(t \\right)}}{72} - \\frac{a_{0}^{3} \\sin{\\left(t \\right)} \\sin{\\left(3 t \\right)}}{18} + \\frac{a_{0}^{3} \\cos^{5}{\\left(t \\right)}}{4} - \\frac{13 a_{0}^{3} \\cos{\\left(t \\right)}}{36} + \\frac{a_{0}^{3}}{9}\\right) + \\epsilon \\left(\\frac{a_{0}^{2} \\sin{\\left(t \\right)}}{3} - \\frac{a_{0}^{2} \\sin{\\left(2 t \\right)}}{6}\\right)$"
      ],
      "text/plain": [
       "a_0*cos(t) + epsilon**2*(a_0**3*(-9*cos(t) - 8)*sin(t)**4/36 + a_0**3*(45*cos(t) - 4)*sin(t)**2/72 - a_0**3*sin(t)*sin(3*t)/18 + a_0**3*cos(t)**5/4 - 13*a_0**3*cos(t)/36 + a_0**3/9) + epsilon*(a_0**2*sin(t)/3 - a_0**2*sin(2*t)/6)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a_{0} \\cos{\\left(t \\right)} + \\epsilon^{2} \\left(\\frac{a_{0}^{3} \\left(- 9 \\cos{\\left(t \\right)} - 8\\right) \\sin^{4}{\\left(t \\right)}}{36} + \\frac{a_{0}^{3} \\cdot \\left(45 \\cos{\\left(t \\right)} - 4\\right) \\sin^{2}{\\left(t \\right)}}{72} - \\frac{a_{0}^{3} \\sin{\\left(t \\right)} \\sin{\\left(3 t \\right)}}{18} + \\frac{a_{0}^{3} \\cos^{5}{\\left(t \\right)}}{4} - \\frac{13 a_{0}^{3} \\cos{\\left(t \\right)}}{36} + \\frac{a_{0}^{3}}{9}\\right) + \\epsilon \\left(\\frac{a_{0}^{2} \\sin{\\left(t \\right)}}{3} - \\frac{a_{0}^{2} \\sin{\\left(2 t \\right)}}{6}\\right)\n"
     ]
    }
   ],
   "source": [
    "display(xsol)\n",
    "print(latex(xsol))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8be665e",
   "metadata": {},
   "source": [
    "#### Let $a_0 = 4$. And let's plot for a few values of epsilon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e4707cdc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\epsilon^{2} \\cdot \\left(\\frac{16 \\left(- 9 \\cos{\\left(t \\right)} - 8\\right) \\sin^{4}{\\left(t \\right)}}{9} + \\frac{8 \\cdot \\left(45 \\cos{\\left(t \\right)} - 4\\right) \\sin^{2}{\\left(t \\right)}}{9} - \\frac{32 \\sin{\\left(t \\right)} \\sin{\\left(3 t \\right)}}{9} + 16 \\cos^{5}{\\left(t \\right)} - \\frac{208 \\cos{\\left(t \\right)}}{9} + \\frac{64}{9}\\right) + \\epsilon \\left(\\frac{16 \\sin{\\left(t \\right)}}{3} - \\frac{8 \\sin{\\left(2 t \\right)}}{3}\\right) + 4 \\cos{\\left(t \\right)}$"
      ],
      "text/plain": [
       "epsilon**2*(16*(-9*cos(t) - 8)*sin(t)**4/9 + 8*(45*cos(t) - 4)*sin(t)**2/9 - 32*sin(t)*sin(3*t)/9 + 16*cos(t)**5 - 208*cos(t)/9 + 64/9) + epsilon*(16*sin(t)/3 - 8*sin(2*t)/3) + 4*cos(t)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xsol_sub_a0 = xsol.subs(a0,4)\n",
    "display(xsol_sub_a0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a8b671ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 42324 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "epslist = [0.25,0.1,0.05,0.005]\n",
    "a = plot(xsol_sub_a0.subs(eps,.4),(t,0,20),show=False)\n",
    "a.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "28c9e769",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 98952 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for k in epslist:\n",
    "    p2 = plot(xsol_sub_a0.subs(eps,k),(t,0,20),show=False)\n",
    "    a.extend(p2)\n",
    "a.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b202aa0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}
