{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $N$ pendula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Consider a system of $N$ pendula suspended on the pendula. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t,g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load('cas_utils.sage')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Triple  pendulum\n",
    "\n",
    "We can easily generalize our symbolic scheme to cases where the number of points is arbitrary.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 3\n",
    "var('t g')\n",
    "for i in range(1,1+N):\n",
    "    var('l%d m%d'%(i,i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_names = [('x%d'%i,'x_%d'%i) for i in range(1,1+N)] \n",
    "xy_names += [('y%d'%i,'y_%d'%i) for i in range(1,1+N)] \n",
    "uv_names = [ ('phi%d'%i,'\\\\varphi_%d'%i) for i in range(1,1+N)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load('cas_utils.sage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_fun, to_var = make_symbols(xy_names,uv_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls = [vars()['l%d'%i] for i in range(1,1+N)]\n",
    "xs = [vars()['x%d'%i] for i in range(1,1+N)]\n",
    "ys = [vars()['y%d'%i] for i in range(1,1+N)]\n",
    "ms = [vars()['m%d'%i] for i in range(1,1+N)]\n",
    "phis = [vars()['phi%d'%i] for i in range(1,1+N)]\n",
    "phids = [vars()['phi%dd'%i] for i in range(1,1+N) ]\n",
    "phidds = [vars()['phi%ddd'%i] for i in range(1,1+N) ]\n",
    "\n",
    "showmath(phis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2u = {x1:l1*sin(phi1),\\\n",
    "       y1:-l1*cos(phi1) }\n",
    "\n",
    "for x_prev,x_,y_prev,y_,l_,phi_ in zip(xs[:-1],xs[1:],ys[:-1],ys[1:],ls[1:],phis[1:]):\n",
    "    x2u[x_] = x2u[x_prev] + l_*sin(phi_)\n",
    "    x2u[y_] = x2u[y_prev] - l_*cos(phi_)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## dAlembert "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform_virtual_displacements(xy_names,uv_names,verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dxs = [vars()['dx%d_polar'%i] for i in range(1,1+N) ]\n",
    "dys = [vars()['dy%d_polar'%i] for i in range(1,1+N) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dAlemb = sum( (m_*x_.subs(x2u).subs(to_fun).diff(t,2))*dx_ for m_,x_,dx_ in zip(ms,xs,dxs) )\n",
    "dAlemb += sum( (m_*x_.subs(x2u).subs(to_fun).diff(t,2) + m_*g)*dx_ for m_,x_,dx_ in zip(ms,ys,dys) )\n",
    "dAlemb = dAlemb.subs(to_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#showmath(dAlemb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dphis = [vars()['dphi%d'%i] for i in range(1,1+N) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eqs = [dAlemb.expand().coefficient(dphi_).trig_simplify() for dphi_ in dphis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(eqs[1].trig_reduce())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol =  solve(eqs,phidds)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(sol[1].rhs().trig_reduce().operands())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars= {m_:1 for m_ in ms}\n",
    "\n",
    "for i,l_ in enumerate(ls):\n",
    "    pars[l_]  = 1/ (i+1)\n",
    "pars[g] = 9.81"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phidds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode =  phids + [sol_.rhs() for sol_ in sol]\n",
    "ode = map(lambda x:x.subs(pars),ode)\n",
    "\n",
    "times = srange(0,5,.01)\n",
    "\n",
    "ics = [0]*(N*2)\n",
    "ics[-1] = 33.01\n",
    "ics[:N]= [0*pi.n()]*N\n",
    "numsol = desolve_odeint(ode,ics,times, phis + phids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_subs = lambda ith: {phi_:numval_ for phi_,numval_ in zip(phis,numsol[ith,:N])}\n",
    "\n",
    "@interact\n",
    "def _(ith=slider(0,numsol.shape[0]-1,step_size=10)):\n",
    "    xnum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in xs]\n",
    "    ynum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in ys]\n",
    "    plt = line(zip(xnum,ynum),\\\n",
    "        xmin=-N,xmax=N,ymin=-N,ymax=N,marker='o')\n",
    "    plt.show(axes=False,figsize=5,aspect_ratio=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euler Lagrange formulation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin = 1/2 * sum( \n",
    "        m_*x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 +\\\n",
    "        m_*y_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 \\\n",
    "                 for m_,x_,y_ in zip(ms,xs,ys))\n",
    "\n",
    "Ekin = Ekin.trig_simplify()\n",
    "Epot = sum(m_*g*y_.subs(x2u) for m_,y_ in zip(ms,ys))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(Ekin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(Epot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = Ekin - Epot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ELs = [L.diff(phid_).subs(to_fun).diff(t).subs(to_var) - L.diff(phi_)\\\n",
    "       for (phi_,phid_) in zip(phis,phids)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ELs = [(EL_/l_).trig_reduce() for EL_,l_ in zip(ELs,ls)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(ELs[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol =  solve(ELs,phidds)[0]\n",
    "#show(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars= {m_:1 for m_ in ms}\n",
    "\n",
    "for i,l_ in enumerate(ls):\n",
    "    pars[l_]  = 0.3+1/ (i+1)\n",
    "pars[g] = 9.81"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode =  phids + [sol_.rhs() for sol_ in sol]\n",
    "ode = map(lambda x:x.subs(pars),ode)\n",
    "\n",
    "times = srange(0,5,.01)\n",
    "\n",
    "ics = [0]*(N*2)\n",
    "ics[-1] = 33.01\n",
    "ics[:N]= [0*pi.n()]*N\n",
    "numsol = desolve_odeint(ode,ics,times, phis + phids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_subs = lambda ith: {phi_:numval_ for phi_,numval_ in zip(phis,numsol[ith,:N])}\n",
    "\n",
    "@interact\n",
    "def _(ith=slider(0,numsol.shape[0]-1,step_size=10)):\n",
    "    xnum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in xs]\n",
    "    ynum = [0]+[x_.subs(x2u).subs(pars).subs(phi_subs(ith)) for x_ in ys]\n",
    "    plt = line(zip(xnum,ynum),\\\n",
    "        xmin=-N,xmax=N,ymin=-N,ymax=N,marker='o')\n",
    "    plt.show(axes=False,figsize=5,aspect_ratio=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.6",
   "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
}