{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler Lagrange - pendulum with  oscillating support \n",
    "\n",
    "We define Lagrange function as a difference between kinetic and potential energy:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:Lagangian}\n",
    "L = E_k - E_p\n",
    "\\end{equation}\n",
    "\n",
    "Then equation of motion are given by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:EL}\n",
    "\\frac{d}{dt} \\left (\\frac{\\partial L}{\\partial \\dot{\\varphi}}\\right ) - \\frac{\\partial L}{\\partial \\varphi} = 0\n",
    "\\end{equation}\n",
    "\n",
    "Since this formulation in invariant with respect to the change of system off coordinates, we can use it for many problems in mechanics with constraints. \n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System definition\n",
    "\n",
    "Let us define a system, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load('cas_utils.sage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('l g w0')\n",
    "xy_wsp = ['x','y']\n",
    "uv_wsp = [('phi',r'\\varphi')]\n",
    "\n",
    "to_fun, to_var = make_symbols(xy_wsp, uv_wsp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Horizontal oscillations of a support point\n",
    "\n",
    "We parametrize the system similarily to mathematical pendulum, \n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\label{eq:parametic}\n",
    "x = a \\sin\\left(\\omega t\\right) + l \\sin\\left({\\varphi}\\right) \\\\\n",
    "y =  -l \\cos\\left({\\varphi}\\right)\n",
    "\\end{eqnarray}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# horizontal\n",
    "var('a omega t')\n",
    "x2u = {x:l*sin(phi)+a*sin(omega*t),y:-l*cos(phi)}\n",
    "showmath(x2u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1: Kinetic energy**\n",
    "\n",
    "\n",
    "We have to write kinetic energy in terms of generalized coordinates:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:Ekin}\n",
    "    E_k = \\frac{1}{2}(\\dot x^2 + \\dot y^2)\n",
    "\\end{equation}\n",
    "\n",
    "Using transformation dictionary, to generalized coordinates we have $E_k(\\varphi,,\\dot\\varphi)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ek = 1/2*sum([x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 for x_ in [x,y]])\n",
    "Ek = Ek.trig_simplify()\n",
    "showmath(Ek)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2: Potential energy**\n",
    "\n",
    "Similarily we have to express potential energy:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:Ekin}\n",
    "   E_p = g y\n",
    "\\end{equation}\n",
    "\n",
    "as a function of $E_P(\\varphi)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ep = g*y.subs(x2u)\n",
    "showmath(Ep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3: Langranian**\n",
    "\n",
    "Now we have Lagrangian $L(\\varphi,\\dot\\varphi)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = Ek - Ep\n",
    "showmath(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivation of equations of motion\n",
    "\n",
    "Using Euler-Lagrange formulas \\ref{eq:EL} we write equation of motion in generalized coordinate $\\varphi$. Note, we can differentiate over $\\varphi$ and $\\dot\\varphi$. However to perform time derivative  we first replace symbols representing  variables with functions (i.e. Sage symbolic functions) of time. We have ``to_fun`` dictionary which automatizes this step. Then we use symbolic differenctiation ``diff``. After this operation we bring thhe result back to symbolic variable $\\varphi$ and $\\dot\\varphi$ with ``to_var`` dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EL1 = L.diff(phid).subs(to_fun).diff(t).subs(to_var) - L.diff(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(EL1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "### Small angle approximation\n",
    "\n",
    "Let see what happens if oscillations are small. We can expand in Taylor seried the equations of motion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_lin = EL1.taylor(phi,0,1)\n",
    "showmath(eq_lin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('alpha,omega0')\n",
    "eq_lin2  = (eq_lin/l^2).expand().subs({a:l*alpha,g:l*omega0^2})\n",
    "showmath(eq_lin2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the equations are essentially equivalent to forced harmonic oscillator. The difference might be that the effective amplitude of forcing depends on forcing frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assume(g>0)\n",
    "assume(omega0>0)\n",
    "phi_anal = desolve((eq_lin/l^2).expand()\\\n",
    "                   .subs({a:l*alpha,g:l*omega0^2})\\\n",
    "                   .subs(to_fun).subs({l:1}),\\\n",
    "                   dvar=Phi,ivar=t,contrib_ode=True)\n",
    "showmath(phi_anal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(eq_lin2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical  integration\n",
    "\n",
    "We can numerically compare if the linear approximation works for selected initial conditions and parameters. For this purpose we need to solve Euler-Lagrange equation for $\\ddot\\varphi$, and for following system od 1st order ODEs:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\label{eq:ode}\n",
    "\\frac{d\\varphi}{dt} &=& \\dot\\varphi\\\\\n",
    "\\frac{d\\dot\\varphi}{dt} &=& \\frac{a \\omega^{2} \\cos\\left({\\varphi}\\right) \\sin\\left(\\omega t\\right)}{l} - \\frac{g \\sin\\left({\\varphi}\\right)}{l}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Note that we threat $\\varphi$ and $\\dot\\varphi$ as independent variables. Since in Sage we use formulas where there are represented by different symbolic variables: `phi` and `dphi`, there will be no confusion of \"dot\" and derivative operator. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = EL1.solve(phidd)[0].rhs()\n",
    "showmath(rhs().expand())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear system can be derived in similar way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs_lin = eq_lin.solve(phidd)[0].rhs()\n",
    "showmath(rhs_lin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars = {l:1,g:1,a:.03,omega:1.31}\n",
    "t_end = 60\n",
    "w0 = sqrt(g/l).subs(pars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plug the system of ODE into ``desolve_odeint`` solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode = [phid, rhs.subs(pars)]\n",
    "times = srange(0,t_end,0.1)\n",
    "ics = [0.0, 0.1]\n",
    "sol = desolve_odeint(ode, ics, times, [phi, phid])\n",
    "line( zip(times,sol[::1,0]),figsize=(6,2), )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_lin = [phid, rhs_lin.subs(pars)]\n",
    "times = srange(0,t_end,0.1)\n",
    "ics = [0.0, 0.1]\n",
    "sol_lin = desolve_odeint(ode_lin, ics, times, [phi,phid])\n",
    "\n",
    "line( zip(times[0:],sol[0:,0]),figsize=(6,2) )\\\n",
    " +line( zip(times[0:],sol_lin[0:,0]),color='red')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that for small oscillations the for some time. Then they diverge. One can experiment and simulate both systems for longer times to see that the divergence grows. Also larger amplitudes of driving will make them differ significantly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vertical oscillations\n",
    "\n",
    "In the case of vertical oscillations of a support point, the transformation to generalized coordinates reads:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\label{eq:oscil_vert}\n",
    "x = l \\sin\\left({\\varphi}\\right)\\\\\n",
    "y= -a \\cos\\left(\\omega t\\right) - l \\cos\\left({\\varphi}\\right),\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vertical\n",
    "var('a omega t')\n",
    "x2u = {x:l*sin(phi), y:-l*cos(phi)-a*cos(omega*t)}\n",
    "showmath(x2u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us once again calculate Lagrangian and derive symbolically equations of motion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ek = 1/2*sum([x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 for x_ in [x,y]])\n",
    "Ek = Ek.trig_simplify()\n",
    "Ep = g*y.subs(x2u)\n",
    "L = Ek - Ep\n",
    "EL1 = L.diff(phid).subs(to_fun).diff(t).subs(to_var) - L.diff(phi)\n",
    "showmath(EL1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = EL1.solve(phidd)[0].rhs()\n",
    "showmath( rhs.expand().collect(sin(phi)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to obtain a linear approximation for small $\\varphi$, as in previous case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(EL1.taylor(phi,0,1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that in this case we do not obtain a harmonic oscillator. Forcing term is multiplied by $\\varphi$, i.e. the linearized equation is fundamentally different."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stable inverted pendulum\n",
    "\n",
    "Vertical driving of a support point as a remarkable property - under some conditions the upper steady state can become a stable one. \n",
    "\n",
    "For example for  following parameters and initial condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars = {l:1,a:.2,g:1,omega:10.}\n",
    "w0 = sqrt(g/l).subs(pars)\n",
    "\n",
    "ode = [phid, rhs.subs(pars) ]\n",
    "times = srange(0, 60, 0.1)\n",
    "ics = [pi-1e-1, .0]\n",
    "\n",
    "sol = desolve_odeint(ode, ics, times, [phi,phid])\n",
    "\n",
    "#plt = line( zip(times,sol[::1,0]),figsize=(8,3), ticks=[None,pi/4],\\\n",
    "#     tick_formatter=[None,pi],gridlines=[[],[pi.n()*i for i in range(-100,100,1)]])\n",
    "plt = line( zip(times,sol[::1,0]),figsize=(6,2) )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that for above inital conditions the pendulum oscillates around **$x=\\pi$ state** which is normally unstable fix point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pendulum = [vector([x,y]).subs(x2u).subs(pars).subs({phi:phi_,t:t_})\\\n",
    "            for t_,phi_ in zip(times,sol[:,0])]\n",
    "o_point = [vector([x,y]).subs(x2u).subs(l==0).subs(pars).subs(t==t_)\\\n",
    "           for t_ in times ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@interact\n",
    "def draw_pendulum(ith = slider(0, len(pendulum)-1,1)):\n",
    "    p1,p2 = pendulum[ith], o_point[ith]\n",
    "    plt = line( [p1,p2],xmin=-1, xmax=1, ymin=-1.4, ymax=1.4,\\\n",
    "               aspect_ratio=1,figsize=2,axes=False, title='t=%0.2f'%times[ith])\n",
    "    plt += points([p1,p2],color='red',size=30,gridlines=[None,[0]],\\\n",
    "                 figsize=3,axes=False)\n",
    "    plt += line(pendulum[:ith],thickness=0.9,color='gray',zorder=-10)\n",
    "    return plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#draw_pendulum(1200).save('inverted_pend.png',figsize=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time evolution of the inverted stable pendulum:\n",
    "\n",
    "![Inverted pendulum stalibized by oscillations of a support point](images/inverted_pend.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System with damping\n",
    "\n",
    "Adding damping to the system will make inverted state an stable atractor.\n",
    "\n",
    "$$\\ddot\\varphi  = -2\\gamma \\dot\\varphi + ( -\\omega_0^2 - \\frac{a}{l} \\omega^2 \\cos(\\omega t))\\sin(\\varphi)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('omega,omega0,gama,t,a')\n",
    "pars = {l:1,a:0.152,omega0:1,omega:14.,gama:.1}\n",
    "ode = [phid,\\\n",
    "       (-2*gama*phid+(-omega0^2-a/l*omega^2*cos(omega*t))*sin(phi)).subs(pars)]\n",
    "times = srange(0,30,0.01)\n",
    "ics = [0,2.1]\n",
    "sol = desolve_odeint(ode,ics,times,[phi,phid])\n",
    "line( zip(times,sol[::1,0]),figsize=(8,2), ticks=[None,pi/4],\\\n",
    "     tick_formatter=[None,pi],gridlines=[[],[pi.n()*i for i in range(-100,100,1)]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.7",
   "language": "",
   "name": "sagemath"
  },
  "language": "python",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "nbTranslate": {
   "displayLangs": [
    "en",
    "pl"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "pl",
   "targetLang": "en",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}