lagrangian double pendulum

lagrangian.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPython console for SymPy 1.0 (Python 3.5.1-64-bit) (ground types: python)\n",
      "\n",
      "These commands were executed:\n",
      ">>> from __future__ import division\n",
      ">>> from sympy import *\n",
      ">>> x, y, z, t = symbols('x y z t')\n",
      ">>> k, m, n = symbols('k m n', integer=True)\n",
      ">>> f, g, h = symbols('f g h', cls=Function)\n",
      ">>> init_printing()\n",
      "\n",
      "Documentation can be found at http://docs.sympy.org/1.0/\n"
     ]
    }
   ],
   "source": [
    "from sympy import *\n",
    "init_session()\n",
    "from sympy.solvers.solveset import linsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create symbols for the variables in the system\n",
    "\n",
    "phi_1, phi_2, l1, l2, m1, m2, g=symbols(\"phi_1 phi_2 l1 l2 m1 m2 g\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# the angle phi is a function of the time\n",
    "\n",
    "phi_1=phi_1(t)\n",
    "phi_2=phi_2(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCcAAAAgCAYAAAAhfxIPAAAABHNCSVQICAgIfAhkiAAAFLdJREFU\neJztnXn0HlV5xz/ZCGERQghQFAgJiUZoIHgQm4aQqD3SFEHKVgq2skjxIBYRK9iCLAVEW+mpgBsc\nf4TKUjhIaytUKyI1Ra2VClSgtuDGEmtAKSiWJf3jO9OZd95Z7+zzPp9zcn55Z+adO/fO9977vPc+\n97lgGIZhGIZhGIZhGIZhGIZhGIZhNMo7gZltP4RhGIZhTBDLgBVJJ6c3+CCGYRiG0RSnA08Di9p+\nEKOTnAd8H3ihofRMj4ZhDK0dGFp+usgQy/he4Ghg37YfxDAMwxgG+wPvRj/wvgisavVp4lkEfK/t\nh3CgD2Xbd9YCHyn4nbLvpY96NC0aRrX0sR1IY2j56SJ9LOM8fccW3rmtm3sswzAMY4hsBVwS+nwU\n8HPg5e08TiKnAZ9s+yEK0pey7TNbAPcAWxb4ThXvpW96NC0aRvX0rR3IYmj56SJ9K+Mifcebgb+s\nItHr0AjOs8DSKm5o9I5dgeNyXvsXJdLZGdg845odgOMjx2Z46W5RIu26sXo0Tl5dTZqmXgncAlyG\nZnvXATuWuN8HgV8D5qL87AmcDbwicl2SRpcBLxG4GL4M2IQ6oDaZCVwInAycCqwHjnC8Vxsag/bL\nti2tNcnZwFkJ56rUfFV6LKNFcNdj21oEeC3S423AfcDVxBu4rv1pn/rhrttdbdJlm2AI/VKYoeWn\naqpoUyat7/g6sLvjM47wYWAjMK2Km3WIWW0/QA/YBriWfGX1ahRwzDWdd8UcXwQ8BuwWOnYYsDJy\n3RJkWLvShBaGWo9cyKurSdPUNsCPGDVK3w/cD2zm+BybIv9eBM5JuDZOo9OQy55/bE/vPssdn6cq\nrkJrMwG2Q7EEtnW4T5saq6ps+6i1JpiBDMf5KddUpfkq9FhGixCvxzgtwrgeq6znLnrcF/gCQZlt\nBdwF/BhYEHO9a3/ah364L3ZXG1Sd37z1A/LZBEPol8IMLT9JlPkNULZNmbS+4xRUZqX5CvBZx++u\nRKMoXeRTdNNlcTpykSnihloXHyV/AJNTgL0c0zkfGSNRTgeeY3Q0bxp6d1HOA97qmH4TWihTj4ZG\nXl1NmqYuQsZ4eEeB7YDngXc4Psf3kYvgrcCfA3unXJtHo9d692mTZcAzwBzv8+uBux3v1bbGwriW\nbR+11gQHoR+8aVSh+ar0WEaLEK/HOC1Cth7L1HMXPf49sEfk2HJk5N4Qc71rf9qHfrhPdlfTVJ3f\novXjPJJtgqH1S0PLTxplfgOUaVMmse+YhwY9/n+TDpfdOuYArwPudPjur6MR76cdvtsEs+im98RL\nwJXI1bZNFiI3pW/lvH4p8O8O6cwHZqMKGuUA5AL0XOjYJuAnjEeyvRx4L5otK0rdWihTj4ZGEV1N\nmqaO9J4tvKPAk8CD3jkXHkHugm8B3gN8O+G6PBo9AXgcONPxWarijcBXgV94n9+AfoTOdbhX2xrz\nKVO2fdNaUxwMfCnlfFWar0qPrlqEZD3GaRHS9Vi2nrvo8UDgDuQ27HMP8FNUvmFc+9M+9MN9s7ua\npur8Fq0faTbB0PqloeUnDdffAGXblEnsOzYCTyBvC8BtcGIFcvG8s+D35qCgHn/lkKYBPwR+APx2\ni89wKhoBy4vvyluUN6E1VnGsRK6dUb4GrI4c+wnwKLDG4RnqxrUeDZEiupokTW0NLEb1PspjwGtq\nTj9Lo2u9v3+EOrIFNT9PGk+hDhA02n84+hF6rMO92tYYNF+2bWutKVaTPgtVlear0qOrFiFZj0la\nhHg9tlXPH0HxTqIeo78kmFX0ce1P+9AP983uapqq81u0fqTZBEPrl4aWnzoo26ZMat9xN/ISAUbd\nN302Q+tMFwP/jWZNfgP4G7S2ag0a5bg3JZE4TmfcFe9QNEq0N/D7yIX0SFSgK4A/Q4GQzkBuH/O9\n5zseuZp2gSbzcLl3r1six/cBzkUzCs+hQaergG965/dHA0Mb0DvfFrgYeCh0j8NQ4LKnkIi2R646\np4auWYs8OJKYBbzNu89maKRvnfccV6D1xnlYAVwa+nw0Kq+5aBZlDXA78DnvviADOm5d2XrgEOAf\nc6ZdBVl1CNzqUZNay9IUNKOrSdaUv64vztPsWbQ8bjYy1ouwGQoKuD1B0KKzgP+IXJem0VXAryDX\n653QTMEGgu2umm7br0eeecegHy6fQR3lv+b4btc0llW2ddCm1prSyhzgVcB3Uq4po/kwrnqsSosw\nqsc8WoRxPbahRZ/XoUGzJ0LHdkYDFndGrnW1S12/B2Z3QXIfWSd15Ne1fvgk2QRD6peGmJ86KNOm\nwOT2HQ+QMoA0BzX6N4aOnY3cPE/xPq9n/MdxHu5ldDBkM4JIov+C3FjeQxBA431o/etHkFsbyG3q\nf5ABUwdTFJsVaCMPdzPqOrMCGZSrQ8fWobW+oG1aHmI0ANhSJJDloc93RNL5HVQePi9HFS6J3bxn\nOw110u8gWF98BcU8Zj5LfHyNPyB+1gQ0EBC35ulgkkdZ05jCbYYoTx2C4vWoSa1laQqa0dWka2oF\n+oF2fsy5a71zLjsp/BewS+jzcajT2ClyXZJGFyJ9RIMd+rGEuti2J9E1jWWVbV6m6IfWmtTKXmig\nJQ1XzVdBlVqEeD2maRFG9VhlnqeoxuPiYmRkr4gcd7VLXb9ndpdI6iProu78FqkfYVxtgiSG9l77\nmJ8p3Nos1zalDEPoO95CaEAnuqzjUhTh8+2hY/cgY+DL3sPuR3F3lVehkeHwWtYDkSEyzcvIEyhg\nhu+G8jyaJbkeeNg79qL3r8zWZlXSRh6+hkak8NL9NPBPjL6TjWhd0Fbe+Y+hGXyfB9As/qe9z8vQ\niNjWoWs+x2inuAAZiXHM89L/DArc9BLwqwRCuxi5JOUNhrYl8TN0a5AB+4uYc3OQwRplI6o0TZFV\nh8CtHjWltSxNQTO6Mk0FI91xrnn+WkiX2BdL0DIxn+vRO3p/6FiaRh9Gncy0yD9/1r0vbXsXNZZV\ntnXRltaa1MrOpJdjGc2XpWotQrwe07QIo3psS4tJLETG9yXAP4eOu9qlrt8zuysgqY+sgybyW6R+\nhKnSJhjaex1aftJwbVPKMJS+42eEApCGPRl2QaMtV0RusBKtf3kIrUOZRfGCX864y/D9yB1uGXLj\njO7Huh/wDVQYPgvR9iauQT58riH+Re2K9tX+35hzJzLuVtNGHjagwR6Q2+ASRmfpAd7t/T0aCffB\nmPs8hLaJeQ3qZHdA28j9rffMN6AZLJ8dkXjiuBRp5nLvsz/r5RuYT3p/l5IvINpGVEbRGYPVaJ/z\nOLYjWKcV5knvXklUpQXIV4f8z0XrUVNay9IUwG9Rv65MU6OGbRR/VNulw426973opXUowdZRLhr1\nabNtL0KXNZaXvmutSa1sQ7pxVEbzZalaixCvx9UkaxG6o8cos4Hr0Mzcn0TOub431++Z3RWQppe9\n0ezztITzUe5B7uNJNJHf1bjVjyyboAhdf69F6Xp+qmyz2uhDhtJ3/IxQHQoPThzufb4t8oUDCQra\nX0tzX8yNt0YjwmcwHlRrR2SAhPEzsQb4OWqYo+lORY4dhNb2faVg+lGS3EGn0LZA38v4vk9Vedgf\nuQlug9YaXUhy0JEnCdxbF3h/f5RwrX/dCzHn/LW6e6BKtj9aC7wWud5+GHUU13vXzUQjclFmoVG5\nsMGwF6OG41Lvb3gGLe19/Sea5QpXjD2RjuLePWgA7Ksxx2eSvk66Ki1AvjoEyfUoTQdV1pe0dBZ4\nf5M0BfXrykVTWXWoj5ragDqYuCjNW6I2tegPxrtQ53VA5PgMtN7ZJ62tz6LNtj0vXW+38tJ3rTXZ\nh25Oeh7KaL4MdWgRxvWYpUXojh6jXIUi1p8bc87VLnXph8HsrjBpevk2wRKWshTNb572oMr6kWUT\n5KUP77UIXbfloNo2y7VNcaWO8oV2+o6nkW0wC3g+vKzDn43/eujY5oy6qKxBmdhE0AADnIRGew8n\nfgeQ2cSPPvn3XB85vwQF0rgzcu1hKBDHMwXTr5syedgKrbW5DFWET3nnk/bXnY2MOFAAEkjeYsY3\n/HaIOTfP+7sBzVr9Eq0p2gWNGN4IfMJLDzTbFZfOdkgn4SBPaxhdS3ks8pzxo6Rnva87UFCXMGtQ\nZ++7c27L6FrmpYzOrvnMJdktsmry1CGIr0d5dVC2vmSlk6UpqF9XRTWVp+z6qKln0WzSLjHn9gD+\nzeGeyxl1I/bZntGYIkltfRHqbtuj6xnz/oPut1tN06bWoJk+9AWCdieOspp31WMdWoRxPWZpEbqj\nxzDnobyHByZ+L/R/V7vUtR82uyugKb0Uye995GsPqqwfYZtgaP1Sk/lp2parGtc2pam+I29f2Ubf\n4ceyGBvU/SDjLo9vQgWwGNjC+9Jp3rlPxNx8E/EBRI5H6++izECxKM6OHD8ZGSnhgBzbeekf4X2+\ngnGS0s/LlMP3y+ZhGUEUc9A6nU3AUQnpfYBgve5MNBr2DzHXHYYE9yzxQc5u8r47E0V4PSMmXz8l\nWNe7hPgfZNPQOuFjQsfC73oZ8vaIBrCC5Pc1k8BFyedGRoX/AWS4AuxO4I4e5TCCWA9FmEp4tjSy\n6hAk16M8OqiivmSlk6WpnZCBVaeuimoqT9n1VVMXoG3Kwq6xi1D+To1cu5jkYEU+NzPqIQH6EbkJ\nOMf7nNbWb+3dY9eMdLrWtkfpertVlin6oTVorg89hPFBEZ8s+yav7l2oQ4swrsc0LUL39AjyJrgg\n5vgnvb+udmmZftjsLlGnXqIUyW/e9qDK+uFqE0QZ2nvtui2XxhTF2qwybYordZQvtNN3rCSYeB8Z\nZbsZFa4/2rsbyuSjwHeRq8V0NAKzH6NbCmbxCPHBq5ajEZg7I8dXo4IIR9ZegBruLyBXz6j7Z1uU\nzcN9aITKD/blj0Z9NyG9hQTrhl4ATvDudWjomvlo68pHkYF3IuNCWoNmH/xRqjMZjSz9CvSuN4Se\n50Xk6hNmE3JL/ENkrE4ncENcg9ZzHc5oAKssXkDuQfuGjs0gcK/aDwVleRxVzqOAjyfcaz9GPRnq\nJKsOQXI9yqODKupLVjpZmnoCuXrVqauimspTdn3V1MfQj7PjQsdOQ9shhqNPr0Kj57dm3O9SNGru\nz8xNQ2ul7wY+5B1L0mgRD7Wut+1db7faoA2tQXN96OMkrwtPs2/q9sysQ4swrsckLUI39bgK6Wch\nijbv/7uBwCPH1S4t0w+b3dW8XorkN297UGX9qMomGNp77botVyVl2hRX6ihfaKfv2IaUmBVvB/4O\nGQ7nocINb0HyPuDzaLZiZvTLJI/czSb+RRzC+BajAF9EgQXDzECjzh9H62umMU4bnhNV5+FaFK08\niQcYn7F6LXov13nfvZhRI+wA79yVKJrr1Siaq89xwB97370I+FM0ahZ127kGBXuK43Dgr5ERe4eX\n1jnI6Ewi632diFyWQAFr1qM9799LYCTOy0jjLsZdwPIwlfFsSWTVIciuRxCvg6q1lpQOZGsK6teV\ni6bS8gT91NQ+aK/oy9Da61sYL8NXIoP2yhz3OwBpYx36gXkhQZn4pGk0TzvbtbY9ia62W2WZoj9a\na6oPnY9+6CUt7XC1b6qiDi1CoMckLUI39biRZHfmC0PXub43137Yx+yu5nHJb5ZNXUX9cLUJkhja\ne+2qLZfGFMXbrLJtiit1lC8023e8jfGtlWPZArlUnlTg5mmV4yYUVKNO2hicqJIT0I/apB+Sy0mP\nllo3K9B7TONdyBUxD3Ube4vIjgKexBTln82lDkG2DqqiqXSyyNJVEU3Vnae2NdUF6q63ZdI8HS2t\n8t0Wl5DPaOxSu1UFU3T/Gesgq/4/yugPxCK4vHcXPQ5Ni9CuHl3LqC/9IwxTM2nkzW8T7zCPTTC0\nfqmJ/LRd/6aovjyL5KnuvqPt8o1yCZoUGWE+mr0IcwSaZYgLkJVEWuV4NfW71pStnPsS7OneNGuR\nWEAjVQtirllHsfdRB7cSVJY4omvB0qi7Mb0CxXxwoagWqqpDeXRQBU2lk5c0XeXVVBN5alJTXaXL\ngxOLGI2ufTXyZsqiS+1WFQxFa0XIU/9vIjk6exYu791Fj0PTIrSrR5cy6lP/CMPUTBp58tvUO8xj\nEwytX6o7P12of1W3WUXzVGff0YXyjXI7odgXvmvGRwnWy4PWyX0IuU2Gt6spw3fQepKqthaqg28R\nbPXUJKtQmX8eBRw8iNG1igBHIreaqt6HK2cil8qk0baujMKtQu8yLmhVHopqoYo6lEcHVdBUOkVI\n01UeTTWRp6Y1ZRRnLYpH4PPGyOckutJuVcWkaS1v/b/du7YpXPQ4NC1Cv/TYt/6RlONDJSu/TdpS\neWyCofVLdeanK/WvyjbLJU919R1dKd8wM9EWp2P5OwY96CVoTcnNwMEFbnwsCqa1CQUremfKtWeR\nHJTKlSLpd42FyHUnuqbyZaFrpqPBia6wEjgl5vhuKOBZFnW/rxloUGBGxfdNo2wdyqODKmgqHRfi\ndJVHU03kqQ1NdZUmZ2iy2oqZyGg/Ge0osR55LK1FA4Y/RH3OypQ0utJuGW4Uqf/zgB+QvB44jTy6\nL6tH02J1uJRR3/pHmDzNZOW3qXeYZhMMrV9qKj9drn+u5M1TE31HV8v3N1GsPsMwDMNwoktuwVeh\ntZkQbDPpB2c6iW7tPmB0g2uANzt8L4/uTY+GYQytHRhafrrIJJfxDRSbzDUMwzAMoHszb8uAZwh2\nL3o92qbSZx3JUe6NyWUx8OUC1+fVvenRMIyhtQNDy08XmeQyXoS2LTUMwzCM3nMGiiHgcxFwPjDX\n+/wwsANaFjev2UczOs4lyGW2SkyPhmEMrR0YWn66yCSX8U1oW+YRpsdcaBiGYRhd5ykUZBlgK7TX\n95fQTPc8tI3vj4G3EsxIGAbAucDvUm0gMNOjYRhDaweGlp8uMqllfAzwTeAb0RMW3M0wDMPoIw+i\nQEpbAvugbbd2RztD3Q+sQMGXHwPuaecRjY7yEooMfjFwm/e5LKZHwzCG1g4MLT9dZBLLeBnwBuCC\nth/EMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAzD6Ab/B/Kf\nu6hF8p8OAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$g l_{1} \\left(m_{1} + m_{2}\\right) \\cos{\\left (\\phi_{1}{\\left (t \\right )} \\right )} + g l_{2} m_{2} \\cos{\\left (\\phi_{2}{\\left (t \\right )} \\right )} + 0.5 l_{1}^{2} m_{1} \\frac{d}{d t} \\phi_{1}{\\left (t \\right )}^{2} + 0.5 m_{2} \\left(l_{1}^{2} \\frac{d}{d t} \\phi_{1}{\\left (t \\right )}^{2} + 2 l_{1} l_{2} \\cos{\\left (\\phi_{1}{\\left (t \\right )} - \\phi_{2}{\\left (t \\right )} \\right )} \\frac{d}{d t} \\phi_{1}{\\left (t \\right )} \\frac{d}{d t} \\phi_{2}{\\left (t \\right )} + l_{2} \\frac{d}{d t} \\phi_{2}{\\left (t \\right )}^{2}\\right)$$"
      ],
      "text/plain": [
       "                                                                       2      \n",
       "                                                       2    ⎛d        ⎞       \n",
       "g⋅l₁⋅(m₁ + m₂)⋅cos(φ₁(t)) + g⋅l₂⋅m₂⋅cos(φ₂(t)) + 0.5⋅l₁ ⋅m₁⋅⎜──(φ₁(t))⎟  + 0.5\n",
       "                                                            ⎝dt       ⎠       \n",
       "\n",
       "    ⎛               2                                                         \n",
       "    ⎜  2 ⎛d        ⎞                               d         d              ⎛d\n",
       "⋅m₂⋅⎜l₁ ⋅⎜──(φ₁(t))⎟  + 2⋅l₁⋅l₂⋅cos(φ₁(t) - φ₂(t))⋅──(φ₁(t))⋅──(φ₂(t)) + l₂⋅⎜─\n",
       "    ⎝    ⎝dt       ⎠                               dt        dt             ⎝d\n",
       "\n",
       "         2⎞\n",
       "        ⎞ ⎟\n",
       "─(φ₂(t))⎟ ⎟\n",
       "t       ⎠ ⎠"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# horrible one-liner for Lagrange\n",
    "# this is kinetic energy minus potential energy\n",
    "\n",
    "L=1/2*m1*l1**2*(phi_1.diff(t))**2+1/2*m2*(l1**2*phi_1.diff(t)**2+l2*phi_2.diff(t)**2+2*l1*l2*phi_1.diff(t)*phi_2.diff(t)*cos(phi_1-phi_2))+m2*l2*cos(phi_2)*g+(m1+m2)*l1*cos(phi_1)*g\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABegAAAAkCAYAAADy1Nr1AAAABHNCSVQICAgIfAhkiAAAG2hJREFU\neJztnXe4HcV1wH+SeBIqNIFEN0IyAoERLTQhxBM9NJvIgBUwvTmUyBAsQWxHQKjG9G4cP0RogWDc\ngEAMmGpMSwIGYRswGDBNVNOL/MfZzd27d8tsn909v+9737tbZnfKmTNndmbOgKIoiqIoiqIoiqIo\niqIoiqIopTOo6ggoiqIAg4HDgOHO8ekVxkVRFEVR2sBQ4OOqIxGC2gWKoiiKoiiKS91tw8HO/88r\njYWiKEoMOwMrO7//E9igwrgoiqIoStM5FZhRdSQiULtAURRFURRFcam7bTgBuBpYNOyGwWEXFEVR\nSmQCMNP5/TQdxasoiqIoSr58B3gb6dzYitoFzeJwYJGqI6FUxmRgStWRUBRFUWpN3W3Dp4GLkI/0\n6s1GURRrGQYs5vy+GVihwrgoiqIoSlPZCrgd+zsGahc0h7nIrDel3ZwDrF91JBRFUZTa0hTb8CRg\nTtWRUBSlPDYGvol0im4DplUam15mAe8go6BeNkOVlaIo9hGms+pK09KjmDEceBaYWNDzs9geahc0\nkx2AMw3vVb1kTh3zagSiFxaLu1FRElDHuqAoihlNtQ2HAU8Ba1QdEUVRimcUcIrneHfgfWDFaqIT\nyATgj75ziwPfLj8qiqIosQTprDrTtPQoZhwHzCvo2VltD7ULmscI4FFgpOH9qpfMqWte7QycW3Uk\nlEZR17qgKEo8TbYNDwJ+WnUkFEUpnsnIztDuSOPiwEKks2wLRwCX+s4dAvQ5f1uXHiNFUZRwgnRW\nnWlaepR4FgUWAOsW9PystofaBc3jWJLNcFO9ZE6d8+oBYNWqI6E0hjrXBUVRommybdgHvAqsnfVB\nVyGjGO8Bk7I+TCmUfuAN4KicnvcFYC/De8/O8J4ViNjZ2GEssJ/v3BDnvSMyvDsrI4D/y/iM1YEb\ngLOQZcHzgGUD7guri4OQZeauf9m1kE7yehnjlYVFgBOBg4HDgHuBr3qu74EsX3odeBP4Usb3VSF/\nZWGDjMVR53ain3bpzaqwvY2I01lJsKFcm5aevDHViSacCmwKLIXI2lrIh8qVskczE3sDj0dc3wjJ\ng5uBx4AfEjz7PS/bo0y7oIkym7W8TMgSdogTdkzEPaqXzGlSXh0KfC/Fu/upxj7Lkl9QD/ny0hb7\nrKnlWnTbUKc+nvbBgqlC9k3reNtsw/OACzLE4//5HjILJ80GU1ORWTU20ld1BEIYjCwbNl0i6jIT\n6RxdkkMclgCuwCyP1gQOz/CeIwPOTwBeAlbxnNsVkScvE8m2fDuLDGwIPIjkeVqWAF6guzE5DulU\nDw2436QuXgF8P0Oc8uAyxIcYwGjgU2DJgt5VpfyZ0EQZCyJLO1ElbdSbZVOHNiIvnWVLuTYtPWGk\n0a9JdWIcC31/nwHfSfGcvLkJOCHk2vrArXRkYhRwFzKzZ1zA/XnYHmXZBbbLbBqKKK8w0obd3olj\nFG3RS3nQpLxa2rl/cML3V2GfZckv9z3+PAvKL1D7rEz7rA7lmsaeKattqEMfT/tgwVQp+yZ1vG22\nYT/wGsnbwx5+Bfw4RbjNMB/FqoIfYJefbi8rk2751heRkaisnIcofRMOJf1o1vFIY+JnFvAh3SNV\ng5Ay8zMX+HrK96eRgUnAL4AB4Ndk+3h6EtKIestsNPAJ8I2A++Pq4v7A6VTbgE4G/oJsTgewJXB/\nge+rWv7iaJqMhZE2nA20UW+WSdVpnUt0G5Gnzqo6rdC89ESRRr8m1YlxPIfYazciH6jXSfGMvBmO\n+IOfHnL9F4je87Ie0tZcE3B/VtujTLvAdplNQ97lFUXasOcDsyOut0kvZaWJefUIstIoKWXbZ1ny\nC4LzLCi/QO0zKM8+q0O5prFnymob6tDH0z5YMFXL/lzC63gbbcNhwAd42sM0X+qHA5sAd6YIdwTw\n7yneWRauLyMb+RPwPPB3CcP9ARl5ysJ45APhI4b3TwJ+m+I9YxAh/UvAtc0Rn4Ufes4tRJa3+Hd1\nPh84BllKk5Q0MvAksCOwLzA/xTu97Iak01tmbzjP3c13b1xd3MH5/y0kX8dljFtatgbuQZQPwFbI\n6P5SBb2vavmLo0kyFkbacLbQRr1ZJlWnNa6NyFNnVZ1WaF56okijX5PoRBOeRZbmfgU4GvjfFM/I\nm42Q/H8g5PoWwO3IclyXR4G36PXvmYftUaZdYLvMpiHP8ooiS9h+ojvWbdJLWWliXt2PfHxJStn2\nWdr8gvA8C8ovUPsMyrPP6lCuaeyZMtqGOvTxtA8WTtWyH1XH22gbfoTU0X73RJoP9FOQJb93Jgw3\ni+CRO8UcV6BNGIb4VJ1K9g3BDkOWCJniLutOynaIn6kgpiJLtPz8Go9AO7wOvEj4TDFbWQxYDRmI\n8fMSsIHvXFRdnAYsjyxpXw5Zarx8XhFNyJvAn53fo4AZwC+BPQt6X9XyZzN5ylgUacNVTZv1ZplU\nnda4NiJPnVV1WqF56cmTpDqxrqyP+Ix9P+T6s4jPfb8rxY/ozGRyycP2KNMuaJrMQn7lFUfasMOB\nNYAnIu5RvWROE/PqSZLtjVWVfZY2vyA8z8LyC9Q+K8s+a2q5ltE21KGPp32wcKqW/ag63lbb8Ek8\nqz2ClogNRXxvrob4w5kPbAP8BPEZNB3xOZV0k8KZ9C4z+TIyUrIOsA+ypHg3JLOmAGcgG1wchfir\nG+PEbz9k6bENlJmGt5z/E4Cnnd+7Iksi3kSMl2WQpRQXIUurtwUuR2beuvefgHSYzkEEZjqiyDdE\nZnv5Z7zsAFwYEa8+5/mbOmlbCpGVh5BNDz4zTN8U4DTP8R5IPi2FjARPB24BfkZnM4WXCPbrdC+w\nC/Dfhu+2Addf1TsB195D9m4YhjSyEF4XxwM/Rz48eFnCd1yW7F6NuLeaiRgHVyJK6uEEz4jCRvmz\nlbxkLI604cqQyTCdeRjS7rVZbxaFjWmNaiOy6Cwb09q09ORJUp1owlBkU9hlgM8Rm20O8DvPPWXb\nv5N87/ezCWIzvOw5twLS0b/Td29W2wOKtQuaLrOQT3mZkDbsBKTOvB5xj+olc5qYVy8g/oi92GCf\n5ZVf0J1nJvkFap9BMfZZW8q1jLYhbbh1ge8i37M+RCYKX4aUgcvGiMeNV5DvlEsCJwNPee6J0hMu\n2gfrYKPsh9XxttqGfyDCDfxwpPJe6zl3LLKU7FDn+F5kZ+gkrIH4qvIylM7OuQ8iyxmOpuOvcjbi\nE/RMxOgHWQrxLtKZKYIBkrkBqSINZ9HZDXgSsozJy9eQdLg85DsG6ey9i/gpO8Bz/lxkhpWXFREl\nHMYqiOFzBKJov0HH3+oFJHNp9GOCN8I9hOCRXxCDLciv006Ej25FMUA2VzADpB+VnOKEPT7g2hXO\ntWU959LURRcb618abJW/KAZovoylCVeGTJroTGiv3iwCW9Oato2Iwta0pqWO6RkgmX5NqhNNeBrZ\nN8hlL6SzuZxzXEX7ezPwo4RhTkY6LFN857PYHkVTR5nNiyLKK23YbenMgsubNpdxUmzOK/cjn4sN\n9lme+QXBeRaVX6D2GeRvn9W1XAfIxz1t3m1DmnBTkIkQ/Z5z85A9e1x2Rj7Ej/Gcm4TU4fU8x3F6\nQvtgHWyV/SL6YFHYXub7A2+7B34XN6chO9oe5Dn3KNIxuIPOaPSd0fHuYT16Z+5sgXRKBiEdkJeR\nkXH3w9MnSIN7NfCMc+4z5y9pZ6koqkjDK8iAB8hGCmPpnq30M7qVUpB/pAXO33jg3zznH0cE2KsY\nxyEdxSCWRmThSmQjjs+BtemMqJ6MLEkx3SBtJMGz1KYjndgPAq4NR4wyPwuQylAn3JG7oI+vrg86\n119X2rro0oT6Z7P82UoZMpY2XBkyaaIzob16M29sTmvebYTNaU1D09ITRhKdaMpEZN8gl6sRnXOc\nc1xF+zuW6NnMfsYjHZlTgPs857PaHkXSFpkNoojyyhJ2CYJXpWSlzWWcFNvz6m26V9dUbZ/lnV8Q\nnGdR+QVqn0G+9lnbytVP3m1DmnCDkAkCd/vCLaCzL84o556LEA8eLk8iXjzcCQYmemIc2gcDu2W/\nzO90dSjzBciK3UWh28XNyshowgV0G1VTkVkQTyE+d/oIrpSLIZXnKHp9eS5Lxz2Ly+POucnI0t6z\nfdc3BH6DJNRlPNKYBzn0j3q/n8sJLoQvIBtpfRxw7QB6l1fklYaNkZHFJZBlHScS7sfpDTozqu5G\nlNQLwE+dd12DzMQy4X/o7pS66R5JRzkui2dEx8dpiKyc7xy7s7/cZ77h/J9EZ5O0qHJagOSB3xjr\nB34YEofRBM/SeYPgZdUueclAnrwWcc0dqXMr91TC66IJeda/JHUvT5LKX1w885Q/aLaMReVlWLg4\nPVeGPs2qM6HZehNEZgc88YrjUTqruvwkTWucjJTZRiTF9nJNSpL0LEq8DWOrfk2iE03xL5P9zHnP\nl4Ejybf9NbUfR9C7eVcYw4CrkNk+3/Zdy2p7FIntdTBP3eolbXnF2URp23IQnVDEBwjbyxiKK+ek\n2Nz+4sRtCCJjn1C9fZa3bEFwnvUTnl+g9pn77rzss7zTCfmXa1H9xSLahjTtwsbIxIVrfWG+6fm9\nI/IhdX5AfJ5CbKcNMNMT2gcTbJb9vPtgUdis61xc+7zHVp/lRHQ7X4C7kcoNcCoyA8cvNAcCc53w\n4wJeOhs4KSRCsxA/n0N9519CRjS8/AMyKjEq4ftNGUgZPksaRiGjmi67I5t4rRjyriPo7sytClyC\nzNZaiBTqTM/1O+ldCgiyXMh/fl9683B3ukdcXfqQdHgV4tp0+/9a33neZs5xXDmd6DzDy1rO/dsG\n3A+i3DcNOL8W6WbvDITELUn4tO5HRiKjeucEXLsd8bXmElYXk5K1/uVV95KSVP5M4pmn/EUxEBEH\n0/BVylhcXgaFS6LnitancToT2qs38yRpWk3KzoY2IoimlWuS9GyBWd22Vb8m0Ykm3IXYzX5eoPcD\nedb2N4lefQb4TlTEPVyB+HQOIi/bI2+aVgeTkKa8TGyiLG35TOCx6Ggnps1lnBTb2183Pgvplq+q\n7LMiZAt68ywuv0DtMzdcHvZZEemE8sp1gGz9xSLahjTtwtec5x0cEddjnXu2Drh2kHNtD+c4Tk9o\nH8x+2c+zDxaF7brOZUsn/PL+Cxc7F5b0nFsUEXq3Qj1Ax+fUeHoJE979kCUrQfwEuNV3biLBibwN\n8fuT9P2mDKQMnyUNk+lsIgayvGEhIhBB/AudZdKTkY0/XFZGRmHfQUZNIbshMx1ZXuRnWefenTzn\njkRGn1y+j4x6+l0phZXTdHqV9+HIzAq3U7ok3b5dLyV4s+OpyIYLSRkIiVuS8Gk/noKMhgf5dXse\ncTPlElUXFwOuR0be48ir/gWV6cKC/iC9/EXpiDzlL4qBiDiYhq9axiA8L4PCJdFzRepTE50J7dWb\neZI0rSYyUmQbUYU+ct87LiB+WdNaVnpM67bN+tVUJ5rwLjKL08+HdG9yBtnb3yR69UmkUxHHXHo/\n5O/t+Z3V9rDNJnDjNC4grrbqVi9zSV9eEG0TZWnLd6Hbt7D3fWn/2lbGZdaHsttfkDb4fc9xlfZZ\nEbLlvsubZ3H5BWqfQX72WRHpDEprUeU6QPr+4lyKaRvStAvTnOPZEfHd17nn7wOuzXGu9WOmJ7QP\nZr/s51XHq9ABRZT59s57RkO38L2FdCy8rmi2QIT9DmTK/QZ0OitRlczPswT7zRyCVFp/B6gfSaR3\n84DRSIZc6RwnWfJWJFnT8BgyouL6GXUL8/ch7xtPZ8nN+siopMufkE0GPqd7oCULLyJLjvy8ivjD\n9/oAm0RHIU5GBmb2c+Jjwt1OOC+bI51d1+fgPyKbFoOMoD7uOfYyhm4fsLawGuGbaIBscLQx3SPT\nExC5uN45jqqLByJyNYPexsVP0fVvUMq/IYjynO38+a+D/fJXJUXLWBRh4Uz1XNH6tAydCfXVm3mS\nNK0mMlJkG5FWX6VJqwlZ01pWekzrtq1yCmY60SVOv/4XvTPA1kNsae8mVHm0v0nsx/dwfFtGsBdi\nN/g/5E91/udhexRhF4CddbBospRXHFnb8j8TvHw9bfkPon1lXGabVHb7CyIf3iX/VdpnRcgW9OZZ\nVH6B2mcuedlnRaQzKK22lKtLUW1D2nbhPmQV4ZYBz9wVWA7xI/8+sHrAPRs44e/BTE9oH8x+2c+r\njsd9M7Jd17m49vn70G1EX49UvLHO8SrIrPcXnUj1Off/DvGP+VDES/zcT+9SAJBOy5L0+rHqR3xv\nvuc5Nw4phFuRgYPfJHh/kWRNw0JkNNKd9TAHOBPxKRXERnTvXv1PdG9+sxJSRq84x30Ej9gEne/z\n/Qcp+8/oHq3Eie8cRACHI7LhKrPpyGjmDIKXGIXxKR3l6zIEmRUBIncfIAbdIGQ06+KQZ21IZ+OR\nMnEr2IiAa9OQEbobI8JfhCy538tz7gjgCTq7P0fVxcuQUXMTbK1/OyJ5dDry0WSDgHtsl78iqVrG\noggLZ6rnytCncTrTTUcb9WaeJE2rSdnZ2kY0rVyTpMe0btsqp2CmE8FMv54GnEVnJtcgZGnr/Uib\n5pJH+5vEfnwFWCYi3tOceI9HBhLcv2vodGzysj3SEGcXNK0OxpG1vOLI2pb/0YnHMPKjbWWcBdvb\nXxA77DnfuarssyJkC3rzLCy/QO0zL3nZZ0WkE+wtVyi2bUjbLnyKfETfAtmLx2UMsA3wMuLP+2DE\nt77XzceqSFntTecDZ5ye0D6Y/bKfVx3P2zasoi0Esc/fIWSvqIOAnzuJnItUPO+Mn9nATchGskGN\nY9Qyr+sQ/zxedkF20PU/6zZkw1ovQ5xnXIxsdDGIXqLeb8JAivB5pmF/JO+D0gbSofNuPrAX8M/I\n0oyTgH9FNkBYGVGCDyJ58jHyUX8kohgf9py/B1GE/4FskrAQWfHwXc97Lqfj98vPDCfsD5x3XIgs\nqYqa7RBXTgfQ+Qi5DjKT7AzgGDqDSkvHvOMu0vn3GoiJWxBjkVlzj9NZUvMaMsK8p+e+1ZHG48KY\n562LzOg7C+n03kDvcqEsddElT9nNWve8zAK+5fw+HfhKxL1J5c8knnnIXxQDBnHwUycZiwsXpeeK\n1qdROhNUbxZBmrTGtYVltRGzEIPJXeo4MSaM7eVadHriyg3s1K9gphNN9evmiG6ah3QcTqR39nre\n9m9c3l8C3BwR5wWELxH2zsLLw/ZIg6ldYHsdzIsyyitLWw4ywStocpYfv16CaN3UljIOoow2qYz2\nF8TP71meYxvssyJkCzp5FpZfoPaZl7ztsyLSCcWX6wDJ29Oi24Ys7cJGTtirkHp+Mr0rrTZ3rl8I\nnId89/K2I3F6wkX7YIKtsp/2O52fIm3DstpCEJn+rcF9jEAaugNNbnaIEt41KX7EqIoP9HmxAyII\nIIUdFI959CqgMpiCdA6jOBJpHE0oqgPnMoH4nc3DGKA6GciTovO4yPcNozPKfzO9I+BBmMpf2fkS\nxIAFcciDNHlpoufyoKz3RFE3vVk0pmkto+xM24gJdGZEgHQUDjIIZ2u5FpkeG+ocNEe/JsEk748k\n3b48SSlKnpPaBbbWQdtIm3YTmbsO2MfgWX69BGa6qY1lXHSbVKYev4Xw/daKJs4+a5ts2SQfRdpn\nNqXThIEK321zH88E7YN1Y5PsZ/lO56co27BsWb4cj7y6X/fHILN5/BEbgszazIMnkKn+6+X0vCI4\nF5nxUTbTkNH+mxA/XNvTu4vvbsjITBV+1e9DljVNiLhnIrK8yAaOorORblKqkgGlw0fIfhibAb8C\nXjIIY5P8xdFWGTPRc3V6Txx105tFY5LWssrOtI3Yge5NPLemd1PPIGwt16LSY0udg/bpV9O8fxhZ\nKh7lP99mktoFttbBJmAqc7c498bh10tgppvaWMZFtkll6vFFEHcEJnEvgjj7rG2yZZN8FGmf2ZRO\nE+pkz9iUb6B9MD82yX6W73R+irANq5DlNYFH/CevQWbLuz6Nl0ec4x9r+NA9ER+eC51nHR5x7xyC\nNw7KQpL328Z4ZImWfxnS4p57BiMf6Kvki8iSo7BlHhcYPKOMcpoGnF3Ac+tGWSO9ScvUdCnz4shS\nflPi5K/OOsI20uSliZ7Lg7LeY0pd9GYZxKW1rLKLaiMWQZYAHwwchgyKfxXpCJ6HDJDPobPRVhi2\nlGsZ6bGtzrWJJHk/FNm8anrBcUpje5i6KkhiF9hSB20lbdqTyNzSwPP0ukEI00uQTDe1oYzLapPK\n1uN/i7jUrZIo+6wNsuXFFvko2j6zJZ02Y3MfLynaB+tgi+wn+U5XhW1YhSy7Xms28V+YiYwSnIL4\nzrke2KnAiCj1ZCpwaMD5VZBN0KpmCFJBh1QdEQuwdSmW6VLmQ5CR7z5kVkQUtsifogRhu94sA1vS\nGtdGXIYYhACjkY2AXN+BB2Lmps+WtELz0qNk40bghILfkcb2MHVVYGoXqMzaw+XAzr5zUXoJzHRT\nW8q4qTr8Guz4zhBkn9mYX0ViS3qLts9sSadSLtoHsyetSb/TtcU23AbZV3Bw3I2KotQT20d6jwAu\n9Z17DlGaLnsgo5evA28CXyonaoqitJjJyAxj1wXIlsD9nuvzCN9wykaalh4lO7tjuAlVCrLYHn67\nwG8TgNoFdWU1ZDN7lzi9BKqbXJqqwycgG7kqiilNrQuKooTTFtvwYmSjY0VRlFLIaymzoihKkRyF\n+Ex2OQk4HljKOX4GGIvMcFi63KilomnpUbLTh+wFtX7F8cjLbYdSD06hY/fF6SVQ3eTSVB1+HbBR\n1ZFQakVT64KiKB3aaBsOQwYX6jiwoChKTcljKbOiKErR7Af8yPk9CpiP+Eo8HOnwzXeu7QOsVHrs\nktO09Cj5MAu4uuI45OG2Q6kPfcANyP5mUXoJVDd5aaIOnwnMrjoSSu1oYl1QFKWbNtqGhyLuJ7tQ\nX92KohTFZOAsZPnRp8jGHmvQvUHHLOBWilt2ryiKYsJ8ZOO6kcC6iN/DVYEngMeBKcgG9y8Bj1YT\nxUQ0LT1KPjyC+KG/A/F5WTZxdoHaBM3jc6RMTwbOBLYjWC/9HtkoTXWT0DQdPhnYiuL3wVCaR9Pq\ngqIo3bTRNhyOuOfaH1hQcVwURWkJupRZURRFUexiY+AuqtmQSl0VKIqiKIqiKC5ttA3PAI4JuqC7\nxSqKUhRvIv5uQZYkzgB+iWwmB6JgPwZeBb5OZ/MfRVEURVGK4QHgWqTzUzZRdoHaBIqiKIqiKO2i\nbbbhdsB45CN9D+riRlGUoohakqhLmRVFURSlGh4E/gZYBXERUBbqqkBRFEVRFEVxaZNtuBpwNLAX\n4s5HURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFUcri\nr0aBdUsG7FM4AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$g l_{1} \\left(m_{1} + m_{2}\\right) \\sin{\\left (\\phi_{1}{\\left (t \\right )} \\right )} + 1.0 l_{1}^{2} m_{1} \\frac{d^{2}}{d t^{2}}  \\phi_{1}{\\left (t \\right )} + 1.0 l_{1} l_{2} m_{2} \\sin{\\left (\\phi_{1}{\\left (t \\right )} - \\phi_{2}{\\left (t \\right )} \\right )} \\frac{d}{d t} \\phi_{1}{\\left (t \\right )} \\frac{d}{d t} \\phi_{2}{\\left (t \\right )} + 0.5 m_{2} \\left(2 l_{1}^{2} \\frac{d^{2}}{d t^{2}}  \\phi_{1}{\\left (t \\right )} - 2 l_{1} l_{2} \\left(\\frac{d}{d t} \\phi_{1}{\\left (t \\right )} - \\frac{d}{d t} \\phi_{2}{\\left (t \\right )}\\right) \\sin{\\left (\\phi_{1}{\\left (t \\right )} - \\phi_{2}{\\left (t \\right )} \\right )} \\frac{d}{d t} \\phi_{2}{\\left (t \\right )} + 2 l_{1} l_{2} \\cos{\\left (\\phi_{1}{\\left (t \\right )} - \\phi_{2}{\\left (t \\right )} \\right )} \\frac{d^{2}}{d t^{2}}  \\phi_{2}{\\left (t \\right )}\\right)$$"
      ],
      "text/plain": [
       "                                         2                                    \n",
       "                                  2     d                                     \n",
       "g⋅l₁⋅(m₁ + m₂)⋅sin(φ₁(t)) + 1.0⋅l₁ ⋅m₁⋅───(φ₁(t)) + 1.0⋅l₁⋅l₂⋅m₂⋅sin(φ₁(t) - φ\n",
       "                                         2                                    \n",
       "                                       dt                                     \n",
       "\n",
       "                                   ⎛        2                                 \n",
       "      d         d                  ⎜    2  d                   ⎛d           d \n",
       "₂(t))⋅──(φ₁(t))⋅──(φ₂(t)) + 0.5⋅m₂⋅⎜2⋅l₁ ⋅───(φ₁(t)) - 2⋅l₁⋅l₂⋅⎜──(φ₁(t)) - ──\n",
       "      dt        dt                 ⎜        2                  ⎝dt          dt\n",
       "                                   ⎝      dt                                  \n",
       "\n",
       "                                                                     2       ⎞\n",
       "       ⎞                    d                                       d        ⎟\n",
       "(φ₂(t))⎟⋅sin(φ₁(t) - φ₂(t))⋅──(φ₂(t)) + 2⋅l₁⋅l₂⋅cos(φ₁(t) - φ₂(t))⋅───(φ₂(t))⎟\n",
       "       ⎠                    dt                                       2       ⎟\n",
       "                                                                   dt        ⎠"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Lagrange says:\n",
    "#  d     d L           d L\n",
    "#  -  x  ---      =   -----\n",
    "#  dt    d phi1'      d phi1\n",
    "\n",
    "da1=L.diff(phi_1.diff(t)).diff(t)    # derive L by phi', derive the result by t\n",
    "\n",
    "b1=L.diff(phi_1)                     # derive L by phi\n",
    "\n",
    "eq1=da1-b1                            # substract the two terms -> this expression is zero\n",
    "eq1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# see above, this code is exactly the same just for phi2\n",
    "\n",
    "da2=L.diff(phi_2.diff(t)).diff(t)\n",
    "\n",
    "b2=L.diff(phi_2)\n",
    "\n",
    "eq2=da2-b2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# constants\n",
    "\n",
    "c_l1=1\n",
    "c_l2=1\n",
    "c_m1=1\n",
    "c_m2=1\n",
    "constants=[(l1, c_l1),(l2,c_l2),(m1,c_m1),(m2,c_m2),(g,9.81)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "phi_2(t)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-29-ed15ada66eb3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdsolve\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0meq1\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mconstants\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0meq2\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mconstants\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mphi_1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mphi_2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:\\Users\\lhk\\Anaconda3\\lib\\site-packages\\sympy\\solvers\\ode.py\u001b[0m in \u001b[0;36mdsolve\u001b[1;34m(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)\u001b[0m\n\u001b[0;32m    582\u001b[0m     \"\"\"\n\u001b[0;32m    583\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0miterable\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0meq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 584\u001b[1;33m         \u001b[0mmatch\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mclassify_sysode\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0meq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    585\u001b[0m         \u001b[0meq\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmatch\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'eq'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    586\u001b[0m         \u001b[0morder\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmatch\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'order'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\lhk\\Anaconda3\\lib\\site-packages\\sympy\\solvers\\ode.py\u001b[0m in \u001b[0;36mclassify_sysode\u001b[1;34m(eq, funcs, **kwargs)\u001b[0m\n\u001b[0;32m   1375\u001b[0m     \u001b[0mfunc_dict\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1376\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mfunc\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mfuncs\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1377\u001b[1;33m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0morder\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1378\u001b[0m             \u001b[0mmax_order\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1379\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0meqs_\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0meq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: phi_2(t)"
     ]
    }
   ],
   "source": [
    "dsolve([eq1.subs(constants),eq2.subs(constants)],[phi_1,phi_2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}