Two Mass-Spring-Damper with SymPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 style=\"text-align:center\">Mass-Spring-Damper Example</h1>\n",
    "<h3 style=\"text-align:center\">MCHE 513: Intermediate Dynamics</h3> \n",
    "<p style=\"text-align:center\">Dr. Joshua Vaughan <br>\n",
    "<a href=\"mailto:joshua.vaughan@louisiana.edu\">joshua.vaughan@louisiana.edu</a><br>\n",
    "<a href=\"http://www.ucs.louisiana.edu/~jev9637/\">http://www.ucs.louisiana.edu/~jev9637/</a></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we'll determine the equations of motion for the simple mass-spring-damper system shown in Figure 1. The system consists of two masses, $m_1$ and $m_2$, connected by spring, $k$, and damper, $c$. There is a force on the first mass, $f$. The position of first mass is represented by $x_1$, and the second mass by $x_2$.\n",
    "\n",
    "<p style=\"text-align:center\">\n",
    "\t<img src=\"http://shared.crawlab.org/TwoMassSpringDamper_Forced.png\" alt=\"Forces Two-Mass-Spring-Damper\" width=50%><br>\n",
    "    <strong>Figure 1: Forced Two-Mass-Spring-Damper</strong>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the SymPy Module\n",
    "import sympy\n",
    "\n",
    "# Import the necessary sub-modules for dynamics\n",
    "from sympy.physics.mechanics import LagrangesMethod, Lagrangian\n",
    "from sympy.physics.mechanics import ReferenceFrame, Particle, Point\n",
    "from sympy.physics.mechanics import dynamicsymbols\n",
    "\n",
    "# initiate better printing of SymPy results\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to define the generalized coordinates. We'll treat the two masses as point masses, and identify them as such to SymPy. We also need to identify the other symbols in the problem, the two masses, $m_1$ and $m_2$, the spring constant, $k$, damping coefficient, $c$, and force, $f$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, f = dynamicsymbols('x_1 x_2 f')\n",
    "\n",
    "# Also define the first derivatives of each\n",
    "x1_dot, x2_dot, f_dot = dynamicsymbols('x_1 x_2 f', 1)\n",
    "\n",
    "# Define the symbols for the other paramters\n",
    "mass1, mass2, k, c = sympy.symbols('m_1 m_2 k c', isPostive = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the reference frame\n",
    "N = ReferenceFrame('N')\n",
    "\n",
    "#Define a world coordinate origin\n",
    "O = Point('O')\n",
    "O.set_vel(N, 0)\n",
    "\n",
    "# Define the two points and set their velocities\n",
    "# m1p = Point('m1p')\n",
    "m1p = O.locatenew('m1p', x1 * N.x)\n",
    "m1p.set_vel(N, x1_dot * N.x)\n",
    "\n",
    "# m2p = Point('m2p')\n",
    "m2p = O.locatenew('m2p', x2 * N.x)\n",
    "m2p.set_vel(N, x2_dot * N.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAA2BAMAAAC//xeaAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMpmJZlQi77ur\nRHZNUE1LAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIoUlEQVRoBdVafYhUVRQ/b75nd2eclP6oP5ph\nMzCw2swko49HIEUQu2VqoOikGARim2BSBA7VH1HkDoSYRLhFRBHV9sVimi1WUAjLQNFfgmMklJis\nYh+Ssf3Ove++efPefW8+vE7ugXc/zjn3d+4777777jkzRHOMNn9eMDrj1NoHjeL1BCxVHSgZNfQ4\nHTaK1xOweDVx3qih5TRcNQrYC7B4LfuvUTt30r6aUcDegGXOGLYzZnajufjZbbBbY8TrrXWURjt4\ndJfSvlzqYhsT2dyGjlJpBy81orT/9/oBOYPnIidSZ2m6rUm3j0fbI232UpiYkNain/NRVvqJfmw9\nsw7wcuWc3RqwJxoxsUFdeXqRxtoe5qWrKFKjRIlFL36sUfKxovCyQCF6Ghfj0T0vvlxgxmVA63kO\n1tKEZvvPToj5fc/lzUR9s7MzghFZRODR62LkQFXi0Uezs5FQvRNmxd0nRxLjQZsxTHdBleaz5ERQ\nrOVE4dH18PdyskoY2S4eVHtAcpvuq+SngsZ48xrD6q2j3hcUazlReLy7ZrCctnaApzVinNlXYchi\ntf/JQgB7OTgP4fGNo46PBsRaRhRevEqULBFd3QGe1ohx5rC4+T57+N4AdHYIrGW4bsCVnAjItYwI\nPHoVI+JlolitfTwMMU+5LXtPT+e2fukiPyNauelda0Rj2+zo8FlbSvMjlDvyzzTRKvTz45LpKzvB\n473g6FcH65QcCcXzwV+ibqr/EBV/p5fw3SsJE1jxTbShELPxDWNe/6h8jWkbOjnNNwTsBp6TfNDi\nOYmESQzAPkPpUigepD2g++Il2lemnZRYNyTMHfAZ7S/9QrTia+bGK/I1pufRSZxjVoBcPJV80OI5\niYTbMJy9lIHpELyAAQOMYDKnUKwRArkjciowsdhnxnofK0DcN/XV5GtMw+CEBNUunko+aPGcRMIH\nwOF9JjETigepcdIkc7AiVxJ96HrhVhidNytJ2B/DvTe8AKeR9MLfQhooFJ5KPmjxnEQCvGBxhkJ4\nIQQvYODiGZpkDlbkMsqdt8SyhAXn2SXdBb/+TXD5HRBvxDGyiH5AJ2Qt8AqXeCSTD3o8Xn9EeCNy\nM1Sg3FA4Hhs2TJpkzjv8OJKl13xesGYoZrP1VOETlMILvDtupBdI7I5hyTYXj+J1Hu94wYcnEwnY\nHVOljC12xzA8xjBO4hl4UPFI07jhcll5ATsEU6ZERX5cdJqKI44X+Bi4zqoTrQY/5BvRwKPNPJx3\nHCYfnjxRQiM5dQ2KciieGGy88CdzsBhxVs6ssZUXTkqTmKZ4/pvO0s6zddnmhfv6DsgX4kqPS0Vf\n2cBzkg96vO1iGN6s7OQoTk12KJ7QM13IZ6BF5c8V6FgBhbX+58ruJQ8LhiiER4hP0KAsa7Y8OzrJ\nBy2ek0jgEzQTnyFb4glNGB91Gk19Ds4jaA/LOCMgsgJRyZzMuIDpq6FaYOPg600w4FQFkk9QRlOx\numCFFir5oMVzEgnOchHRVBDPWv2UzfCx44lHyPpiSJqSwbhscyn7A7jDcMpOCBlnBPiKSOZkV16o\nsy6/ozgi3090nruSdt+NczMOj9LWfG7LXYNbelLJBy2eSiTsF2OtEqogXqyCqBsEpDJ8McJt4mC8\nidBXwXkT39vxZAR47m0kc0Re4V06RblxL5BoOz4Vc382INUzovCcBymcG8RDRLORMWNbfkM5r8Lt\nQKrTG5wLBV3BaxjnmFRdXuoZUAaMEPoOzj1Pk1fkp64KaLjvF9HSZuGu5q6nF4XXyLj58QBwEFMv\noB4QYGJ/wZlFLkcX3xucu0x/g5cUzjHWuLxccYQXeM28nTg3mqyXXfVgI+0ThnuhOzxYvICtuoJa\nemETWiDeSL3E/XgZC6aGRgiJ3VxlBDgroCjCCxm8gEd/PTGam7aVtqb2b1LhXugOD0/tD3ihDtMD\nO07YRNfRrm/QOYEzmyfYbwrOIVaUvumvxCx/0Jny3ozAKskTZYQX6BWPXmjzQZ8k3Avd4eGT+CcC\nljLM9BdwFLP+zr/xFzqTuDzBvuiPgZcuofBQ8gyfamR0z2fezAyE25xL6UV5Qel0VEd4oSOchnIW\na0F4AaxPEZQct3AzHHqQJ9gXfbzy6syHlkOf5G1yovumjIA8+Uily9kL829hWui+EZjxhmoSPmFC\nGEqeYF/0VXBO1rU8conNildj1TjRPR9Z4mXwhp0LFb4Wg4PXvjU4OMVtJ4buvuKVtnFw8PbBwcXG\n8BhI7I5jFdSHiZ6oxacPMU/cNTWCfe67wblQUEUGIbsT3bMX3IwAe0LR5bwW1BwPEO0soPMW7wTF\nmvxrA78R1Aj2m4JzNVDUyZM8VkT3/Ea4GQHELi7NBS/g1MRvPI0QLUFKsFgpo8O7oyfY574bnEPU\noEf7S+jE6yh4d3QzArxDKpoLXohXrLepWKLdlH4PAfo8u4zZIxj3Bvvcd4NzdXNc5+u8r8jonkNI\nNyPAWQFFc8EL1tZTNsWnKLv62wKeeX49L3EsaG+wz303OFc3hzp1dyk9e6Nz3s4NeTIC6hDBynPB\nC56bajQbJ2j5yVN9PkMGyYnuRVAGd8Eb4hypFNv1wu61n6kh0XWb54W28cKsqWCc+CcUkOpv1Q1Q\n0f12KVTRlKtqVd1mVMOawkcqSsGVBSMvV+RptI/nGdTc3O90nWCfZN8qNWvJnoruPRmB+Tq9FryB\nAvWNt9DpRGwAzx/AyL5znyFzyU4IATtMOTFEU8vuG6KBc1pJd0wDeDIYb5iX/acbDF1rDzMbGTed\nSgQvdsasF0zjRUzdrCiu/1W2ayOm8bqeSEcDj5U7Um+pbBqvpUEjCncYQWmAmMZrIF/CVnLELLhp\nPLOzC0NDhssomcYzOrkwMCQN94bJuuGbxutmDp2PWUHEPwYYI9N4xiYWBZRd9pjzR58orfZlpvHa\nt3wxmjFk4oYuBsA3th28/wAowMmoYtW5MwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- 0.5 k \\left(\\operatorname{x_{1}}{\\left (t \\right )} - \\operatorname{x_{2}}{\\left (t \\right )}\\right)^{2} + \\frac{m_{1} \\left(\\frac{d}{d t} \\operatorname{x_{1}}{\\left (t \\right )}\\right)^{2}}{2} + \\frac{m_{2} \\left(\\frac{d}{d t} \\operatorname{x_{2}}{\\left (t \\right )}\\right)^{2}}{2}$$"
      ],
      "text/plain": [
       "                                         2                 2\n",
       "                              ⎛d        ⎞       ⎛d        ⎞ \n",
       "                           m₁⋅⎜──(x₁(t))⎟    m₂⋅⎜──(x₂(t))⎟ \n",
       "                       2      ⎝dt       ⎠       ⎝dt       ⎠ \n",
       "- 0.5⋅k⋅(x₁(t) - x₂(t))  + ─────────────── + ───────────────\n",
       "                                  2                 2       "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m1 = Particle('m1', m1p, mass1)\n",
    "m2 = Particle('m2', m2p, mass2)\n",
    "\n",
    "# Set up the force list - each item follows the form:\n",
    "#    (the location where the force is applied, its magnitude and direction)\n",
    "forces = [(m1p, f * N.x + c * (x2_dot - x1_dot) * N.x), # damping and external force on m1\n",
    "          (m2p, -c * (x2_dot - x1_dot) * N.x)]          # damping on m2\n",
    "\n",
    "# Define the potential energy for the spring\n",
    "m1.potential_energy = 1 / 2 * k * (x1 - x2)**2\n",
    "\n",
    "# Form the Lagrangian, then simplify and printd\n",
    "L = Lagrangian(N, m1, m2)\n",
    "L.simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "LM = LagrangesMethod(L, [x1, x2], forcelist = forces, frame = N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAA/BAMAAAD3QRp/AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARM1UEIl2MiLd76uZ\nu2Z9/d+fAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAL9UlEQVR4Ae1cX2wcRxn/7nbPd3bOl7z0IeqD\nj0SUVK3UI5aAFpAXQyuSBmwhWqgUohWggJJUuJV4QQjfA0KK1DQmBZoqhRjKH/FQ1eJPCBUVp6iE\nojTY6lOkKuFUQalK1bpx0gQBPb4/M7czu7N7e4rPjZqM5Nn5Zn7ft7/5dmZvb+8nw1inU4frJS0D\no53OCoyNj29IA1zvh8L4OObo3Z6Ir3+heWVTfOLdn6Olwtz1HPXIQN2/0APRY9i9jhZnerjFhvvF\nx9wHaSI1b3XX0QNCd6pP1tn4oT6jrSocqU02ryyivY7KSxJta59Bs/GFhT7DrSZ8K/i/ucJ4do6K\ngYTb2WfUHvhdfYZbTfhO2PHt6SsLaOfoJAfbsXuzI2ilTZ2vOEZ64ocDhxd3PRwbYLvWjPXaZkTE\nRcbGErWjHdf9KIpieyjLJGLlqMIfAP6jZdcHwQH2XtdMhOyNH2rFnD71S+moqL2th5V9XNvOY0TE\nQcb2SKFGoCiK7QKwY18TLCJWjkotwldb5Xk6xsoWAP8w+PVYdx58ZTnudKt0FJv2ANnbmjBp98as\niIiDjI1NmQqBoii2C8BvJxpgEbFyNNIg/EhjaCHuB1BrAXg413sTQznwP1JOXqgaK3J8TZn6QPYE\nQCnUHY6jSUST8VIc3NQoqBmF7Kj4b34MwCJi5Wg2IOhUc/Tv3CCjW4abuGLqANu7PbqRA39aYb1Q\nGgW1mQ+rfn0g+wSu13nd4TiaRDQZL3QAsctNjbBmFNuXFgJYRKwcvcrokZnZb9luZO3Hv+FpgOKG\n+FgO/B7l44XSKM3xsSIH6cOa7e9j4z3drmTDJKLJeGESRz1uajRiRiE7KsTNJqJyVNr9FUThJcRS\nOPaBv3LjpU579vwMNwH2Adz4txdDvFupju4hB34sELQXyrHY8vZi4KEW1G57u9x5n/SSXXjmf8cA\nnlcd6pBKRJPxQkHecPw7e8Mb7t+g/KKpUEdqFI3mY+3pIy0YsohIjiqPwCZEnGVUVC0GRZzJg9xx\nDmu8VUCtzqazYnzp9LM0aOFnVaa9kIZwDzS+OtYGGMW/6gXAFJXr2Ms2L/WXyDJKWmBNxgsFPH02\n9N/a4M1r0tKr63zTWR/GiUiORqbh1xjpFh1NHUfruLq++Xu2HsGa1lnWlx/Gfx5+RQ4W/qC6sl5I\nQ5jsu4JP42G4gdULQzNQPjVH3WTTPQ9uJ8soaYE1GS9kcLlxCCoXobasSRshsJlvOgdn4kQkR4sB\nR3sc6/UdKdThP9nGWgj/GFt0qygvY+Xf/BiWQxgtiT8Ms03st/AHcfJQ2rTp5p9s2rSAzcUWVni7\n2IDVdlpxMlmy6Z4Hs/gHk3SOx3gfMhFHYCYTxfX9N8G7AOuWhLSeCk1oGQPmmY5sFpuI5OhJ4pRc\nRzCBpKMcIQOVIwK7CuG/C7xqMEcRPr6OjtwfkjvnyKOoUY6m6JScI0LokhKYc4QYL8QKC27U6hyM\nhJo0d0L5njPSyDWd+xBLxCIinCP/osRQe616SUyAkz/Flqwj3DuFZQjA8cIqhoeJAL0sPK8s7PRC\nrCh5FUqM7LXqnkDniPbaGPgAG2nULEwkGViT8UIB40bFdTjWDPRuZWqjOl6u6byIoWwiso7eRpY4\n9IycyV8Gulfj7ghewFpyhDuiVPdmXPfsGB5+QL4Wfiygrm6Oasv+SglNumfDZ0frOkdkfwY+iB9B\n2G8WIZIMHL9nj0zT5T9RaegcMbUiTAUULd90fo5Im4jk6BT4f8ahPRQJCddV0N0w1dI5wq/21YW7\nsJpmjFnF8PKNxsK/quBeyA28YVyiHBFyKPT/q3NE9ik/BPgTDppFiCQDazJeKGjc028AvBdjy4XV\nUxEr13T4BmETkRyVdr1MJ1GXGxEc9L7z8Mb5UJ1uIz5ZnWvjMySvMKGk6hj+Ne628KcV0gu5gZdp\n35exhVuldLRe69yk7ke0dQ68jgN8o2YoV4pIMrAm44UCxpPuBfjIP7s5Emr+FhrONx3+CmATkRzJ\nGeQe6p/8YuOhQ7SsVLmdj8NNMekB1SpxfGG6QGm08DcpDy9UDXWgR34u6kNc2fHnb8E4AmsyXiiQ\nqGbSmtrdQTSgLrlFzxiF0hJZFhErR7x0t83MBrDZcLuD2/QlkIr+GikW1nH8g+MfDbDbxHdfuvtN\nHDGKrAzs8Oa5V9mlkK1Y5QisyfjNGBaYtKLmP3enMZw9nQ8X24S1iFg5Ks/j8PPwDYDLBJTy0FH8\ncoDlDNd+nQ9GFcc/1enwqIGvLRh4sznaFKvyl/+E1FL2pPTG6mTgJBntIqQVteHOv3U/QI/p/HAb\nQy0iVo74YfsJ3NKF+Siobh0wnXUnHnPgh2nzuUqFF3Y0omxJb9RttyIiOsX2eGSlUSNEFCXCY+vl\n59i0iNg5wgvoX4Zznxha+KTlSUalTfUrVJklD/4fpoPVftiyANiuNWO9thkRSZCxganUCBZFsZ2U\nZRKxc+S1AH5RvtSuhtNOT1dnb3x8tbiiDKSvN7Vcp7VzBO/HNyBf2tcuHEvbHcmgvfHVIOm1Jj29\nqOVUAsRytCbUr5qT5FQCXNM5yqkEuKZzlPkyLFrtOXK0mP/exHH7xUdk1raFPCebeU6ZlaMHJMBU\nnjgGJhsfGsh3tjmVVwmQkaOyesDLFjwk55mNvzHp8A71bM2rBMjIUTEQ8viWo6+SjS+1+wo2QPDO\nFCVA4pQZOTrJYLfgwXwMtWP2xH/PxqdZfAZjMP2MESh6du7xBI4ubp4UK4qiIqfnKEsgoZ6cj6sg\nxsGtQjDx+wy01dyx63MAV6NaIj1H+HoKi1tVUGymaRd64w9SWKvcEbC5meQQcKsM0RnMkn5GA7UF\nv26KbCP9jYCCu3nyYBRFYdNzNNIgiFtVQK9XJvBda0gQq/TGD7ctB/D/8MeAemrz+MIcYIXa3Rc4\nYijbfcYuxNI56DdLXhiNmy03T0LUWvg6a9l8UZaeI3zThsWtKjiMIyfwos0TxCq98VX1cRm5fS2g\n9ugC/TZ2Naol0nMk7+mdqgJ+l+rWLvTGD81TRswiOVq/AIUVKM3xCJ/BwGScMULtxyb9hElKidVU\nS7hydO/uGTwbrhMsTsHDUMutXciF1yuFokuRHE1h0LfgalRLOHK0v1GsI/mzagb6YCgK8HcN3rKJ\n38E0lo6GjsHCl7u/cGq05OjgNJQvDl4toSQb+txyNLjaYg4eduToZ7CugWMZAgn6GbNaR4z8ZMKB\nEpWhY7Dw8hOtCVc5alGOBq6WUJIN8/x4LyTxh0vMwbAoR1rngGIXLo9jrVUF1GEoCkZww9PGl9/l\ntaMLLzoGC1/hTy5D8ADmXltsYZhBqiWYkZ6XpZZwaC6YC76xH5NGtx5ZkmZ8HRkCCZpzJBnoesYa\nkY7Bwqeto/X4ubYCA1dLKMmGcLXUEi4xB8OSOZqiBYJF5SgmeODdRXvHrV3APdi93UQ6Bgtfvizx\no1rW0egS1C6sgVpClBXC01ZLJDUXwjGZo/VtAMxBlkCC7sFO7QK6uQQSFj7tcw2fIav1NVBLiLJC\neNpqiaTmAueDJZmj6gLsmMGRPTyuVQUAhqKAvqY4tQvo4tmCCv5GY+ExF7GC68g/T2qc7TNroJZg\nRl2e8qmTprkQoskcwT2vhzQ2FlBN2g6OYyoKbMkAw7pVDM8/Clv46lIXK41TT93ShiMAd+7+HYla\nBq2WYEZ6XpZaQkY2JsQfjhypKdCNFrSqoDstyTt9F8ESfxpO4pWOwcQXQ/btUSkPvN5zjDQjJFz7\nUUsIVs/r7iAK5ogig+k5qk4jIi54EK2BLRmIzpLEKx2DvnYhYqcC0yGtrTwwR/MMMSMkfPpRSwjW\npZZwRJEzpeeoTOTiggelKBhtsvck12YVxysdg4n/l4lPbSsPWHW1hDByqSUU16T4Iz1HWQIJ9c7s\nTGKKKSoEE/9owsnVoTy6Q2aEbme8cYA7KL06xXGEtlN48nAURaEzcoSrJFVVkPLmNAe+Rls4Rxno\nu9pUnkSsj3e1MBiBhFylHEkaKKQvtUTGOhqMQOLZgc49b/BeagkrTlaOLOA1bGCOrv9Plqzrz/+T\n5ePX/ydLVo7wf7J86P9m7/eIRAcs1gAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}c \\left(\\frac{d}{d t} \\operatorname{x_{1}}{\\left (t \\right )} - \\frac{d}{d t} \\operatorname{x_{2}}{\\left (t \\right )}\\right) + k \\left(\\operatorname{x_{1}}{\\left (t \\right )} - \\operatorname{x_{2}}{\\left (t \\right )}\\right) + m_{1} \\frac{d^{2}}{d t^{2}} \\operatorname{x_{1}}{\\left (t \\right )} - f{\\left (t \\right )}\\\\- c \\left(\\frac{d}{d t} \\operatorname{x_{1}}{\\left (t \\right )} - \\frac{d}{d t} \\operatorname{x_{2}}{\\left (t \\right )}\\right) - 1.0 k \\left(\\operatorname{x_{1}}{\\left (t \\right )} - \\operatorname{x_{2}}{\\left (t \\right )}\\right) + m_{2} \\frac{d^{2}}{d t^{2}} \\operatorname{x_{2}}{\\left (t \\right )}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡                                                     2              ⎤\n",
       "⎢  ⎛d           d        ⎞                           d               ⎥\n",
       "⎢c⋅⎜──(x₁(t)) - ──(x₂(t))⎟ + k⋅(x₁(t) - x₂(t)) + m₁⋅───(x₁(t)) - f(t)⎥\n",
       "⎢  ⎝dt          dt       ⎠                            2              ⎥\n",
       "⎢                                                   dt               ⎥\n",
       "⎢                                                                    ⎥\n",
       "⎢                                                           2        ⎥\n",
       "⎢    ⎛d           d        ⎞                               d         ⎥\n",
       "⎢- c⋅⎜──(x₁(t)) - ──(x₂(t))⎟ - 1.0⋅k⋅(x₁(t) - x₂(t)) + m₂⋅───(x₂(t)) ⎥\n",
       "⎢    ⎝dt          dt       ⎠                                2        ⎥\n",
       "⎣                                                         dt         ⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(LM.form_lagranges_equations())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The <tt>LagrangesMethod</tt> class gives us lots of information about the system. For example, we can output the mass/inertia matrix and the forcing terms. Note that the forcing terms include what might be conservative forces and would therefore normally appear in a stiffness matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFgAAAAyBAMAAAAq1AE7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhBEqyJ2u93N\nZplQnf8bAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAB60lEQVRIDe2Wzy8DQRTHv7o/VLdFcHGyB4lI\nqDb+gK6osxI3iTbizsGFUxMXNz2R3lwQDiJcHDfpX8A/oHFxQFJEcVGz/TH7ZneSbns2l31v3idv\n35uZ/e5grPaOQGOm9obh9EIgFivpJEaCoQ411xbWHu+tZr72cDSh3QaG74CpwPAZkMo06PZlfAGb\nCQLrE/vFHX3ytPk24aH9MrhA4FDkHNlDJAHVFkjmhNmWLZqN2XoZ84aNAxMPUJ+ufDDLLMCZbM5p\nYhdQfLCvDKSAZ2BUBoM1mKIN4gW4hv6tSTLjAnjIkJoxBK2KfntPBrNNYbnqo95guIpYBVHTlMFG\nQmudtTqsX6FvC8qrJYO1ySOLZm7a7OFfDTfmO6LKIA16bfFshJ+/Cl6C+CJMAjLzH6ar0uVqrE2P\n0ywy281cwmpeRpA5DscGEbJJQGZyOMLO6KeMIHMcHtiC/kECMpPD2TL0Hw/hFR8OM91Qvb8AV3zi\nl84368JlP8zFJ5TvtSksK4OLj5FXvynMGlR8DbbEx8gxASBlRI4R8y0dF5/GsvKa2ab0O3UJg4sP\nDOfb5DBusGQJJJNbLj6IOyEXXi5ue1i44hMri7CXFPx1bDDfzSwEPY5amj0JDPfUapXAcPNFwcro\nDu7oCtHJ5eQPvcqDCnUKW+cAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}m_{1} & 0\\\\0 & m_{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡m₁  0 ⎤\n",
       "⎢      ⎥\n",
       "⎣0   m₂⎦"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Output the inertia/mass matrix of the system\n",
    "LM.mass_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAA1BAMAAAAnolWHAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhB2It3Nq5m7\nZkQincJbAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMmUlEQVRoBe1bb4xcVRU/O2/ezr+d2VWqQkJg\nQrQCDWEIhoAa+4j4gVTdwUSgVt1HhBJbP8w3mxTckaKF+oeJVWAlgUGxYMQ4qVChxXRMRARpXEiE\nVEMYjR9Mv7ARKhal4+/cf+/e9968t/uB9osn2XvPPed3zrlz3r133ntnliiDDgcZyhTVWvEpLk6P\n6NzxP6PAn45Y5hbcYe4oGx/m2p8OQGE8nqMzrvqEiV0cGlYwV7jD3FE2/liu/ekA+Fd9fI7WWZEr\nbWsAdoM7zB1l4wvLuQ5OC6DoJuGAM4ktW/c6YzkoiY/ycoomF/+NFKMUkYyQUNzKknonIU8RRHNM\nm6cwsLy5SSgdtx36txedsdLdIvpqcjL5+O22e+a3XPR+Ibq7valFnz9TqRFh51O/VQPTlYaCfdgI\nFJMCpWiOKfOkLds7ZHtzk1AY2BFqA2iTdA+Rv5v8XkKTj1+M2+yl60Uy5eF8n1LfQ36fLh6pke4q\nQAK9WY9V76dAKZpjyjzpRxtbZHtzkzDVsiNMtab79ljy9QGRt0K0Xqm8UDGUjy8va6zs63NU6DH7\nvX0dtK8r6YCqbZqakyPT/gXcRqJCaCSCSYPac0zO0z/5RSLbm5uE+bYdYKEz8ydHIJRlTLfWI7pG\nQb1QMZSPrw01VvYzffLEjuvzuKE2HyJMdan6hsSYdje4J7EK54xEMGlQe47JefIlJNubm4QXHP9T\nwfynHIEY7EBbbhJVRmJIXih7onz8dOwDzPapIS5/n30UutwSIULleCIJpS5Ut+HvfvxZlAJlD2aO\nyXlyHMebTsL6rQEMkWiLGvsu+7UYHh0vz58IlAaH27E/PhdSbSAFXqgUlI/X11pbLAyo8S8e/Pyl\nx/DJB962AAN5fJaxLuzA08A++t99RIfYwKE4lD2YOSbmWX9wz4AcbyoJO1qVHvw+7/iOBofblYDo\naiE4ghZbk+psAPJCNDES+MJTj7PYxhdjS3yxSUVxu/o0LY5ooXXza8vKgui1Jlgr8AxUYiEfhdyl\nCJoSMznP2ZDI8aaS8G6qtuD4h65zM5rp3Ui06REx3oWWV4zXFUPyQtnbrcB/kH7KMhtfOmmjiBYH\nKgmYVZc2Xtv+grYg2sOsFbiMCfJZRJezwqEImhIzOc/FAJvF9iaTUFNH0rdw7tx7B+jOgGbHkjic\nf/YyWhn9HHC8NYsraApLS/eetbTUB5vA76b5DuQ2viQOgM0c4Y7vQ2e2Az7fm3R4ABGILdRuswJP\njeQ+p3lWm2m60JSY8XnKdex4k0mYGrLjySuBNo6gNUnw+YIK5+i9EE2cGP9dXuLiIxl8fCXgYPQ4\nL9UeFU7Snq+E4FUSsK2ZosA87QX8ySSgN2RBU2Im53k+LB1vMgkLTekxth1qZgsfeA8AlwsQlndj\nhdrU6EobL5Q92hieNrYhtPHFNw1WMDNDqvMinOpR7bh/UuWIN5A3oK8yJArMC/g18on+wXKLHGgy\nZnKez8HY8SaTMLtMhBj0qOUbrL9ClYBFhfazaGUScNAVel5A9R5rnJUQw9M3WW/j498OuFkS27yG\nWfXqK/7rBWVBm4huBm8F5qPsAvoQvjIYY5EDTcZMzvO9MHa8qTOhT1sCqG6wfIP1erTQZtFW3r0q\nCRtwxfvXommyxklCDC9vwm08PrRLu+iawD9BpQG90KkO6Q2RBFiUbjvvYBdQKzB7O+iHRPtdFy40\nGTMxT7E5HW8yCXTeK/CO5SY+sgkCpLj655+gV0+EKglYjaUjy/hWDyTOC2WPNobnO1Oxeg2+NhSi\nqPnM1ofE98DOCx/ni7P9JlYhQgWHcpfIDszL+pZXoOYT1SIXmoyZmKdYjo43lQTllI8LQ/6BL7d2\n3vmYGask8C0p0w7ZkRcqJo5vNBsBVDa+orHKJL3TFkYr9+FuOS51jTzJpMRMzLMwZDvbm5sEvXSE\n9+uD+TbtFaxsLhEdP5wwrZcd+SopFMdffelH28DYeLm5lOHETlsYgAwsL3LiAcqgmEmJGZ/nhyvL\njLS9uUnAKKJDdCWRdZzvvGufUD4j2uQjahx/7ngcx/858p7FyQgGoQLPyGxvNvIUJhkzMc9vXy/s\nbG9uEujHluMzaBs17KwonXxhoZysDX+7Bc9gZYQ4oDQUkliG4ig5juaYmOdLvxIQ21ssCVaa/Tfp\nyGen+9clopTEcno5Ll8Fvt6MG6WPZYSE7laW1OV6SChdQTTHxDw10PIWSwJuPAytK76xXAubZpzH\n5OPTr3Ce33deH0sCfSwKeexL25cb+4JIkMPl48VTZY6X06GOJyF7DtF9cTZOa9eK13anuM9JQnHk\nzAf3xWuibPz0mny9g+CcJHzSDe313HHeKBvf6OfZnyJ9ThJ+4E4D98Vrohz8RWty9s6Bs5NQC+3I\nfF9sjyUffRvFdbn4cjtukjqe8I2psZMnoBFRb74XnW/a7CQsOrMU98WRQ8XJ771qJ6GQ99EJsYWf\nHsS0/v4XAxZVbir+kvzfdJnHY5PsuM2uNkU4cKVDzweOgGVDIXkYLf9pyk4CHmotEvfF1liy9+C9\nQ2pBCq+Fr0ygySoOUWklpq+04Ao0NR43+e2zVCOCIj+72qRhov8cv6yKUaVjiljWbSFlJ+EMxwnf\nRyeIn3Y8fBj9nEKXGUg+/jsGKxnUPS5grnLhzWhnW0JqPU9VUwpT8Ql4obCinxC9S3JRy09NqohV\nCCNxZhLceqy4L44sFcfPvbVeVJCKkrAK/FMxd89hjm3IqkKu3m5YT9ZTKYWp+AS8UBjTmUS/CCRr\nWl5mT2Llzsk/Lc9MQuxFEN8XazvT8+N6uYmLN1KiaCXk429QNrp7C+91WhjIJJwvxfqFAEZ51SY2\n8EI0oAfaiSSIVxG6iHW/hHE7IQmFrTdCyVfYIr4v5mH9gbeLY/1+Z7tV7BHYKAn5eHWtdQz/P0hC\niFH1le0B0X102UMYIIJFmdUmxnkhGkl3td254iCOiliHNGpSEkq7aAmYmWEEtDmUKZCDokjRESjs\nglS0HSILid/whzZEDn4+iEDgSqhGzTfBzLTx4tl/ffprb2PAFhFF1abciliDi1t2bH67yueXeFd7\nNPKZvhKmmuLNQqUfAR3u2emAige7LNuFP2wzU5BKSwIxvtCpctYc/OIIkohKWAkiCRD9HvXqm3zM\nWlhEmD3MWoUp1yF0XohGULnPnRW73FKr+3LI+U+RmwRdGzrcFurygEiLSNajxiusuYavjtdl9hz8\nmYIUXbC09PWlJVG+iOPLHVF1cPCLmBVFxSSzHSA93KkhJUywiCA1TAnjVVTEdOWW56piT43k+SXq\nN6KQxc4mbYezpXLiSvBOAmCSYApMwio6E6QTbhlfHonaCn8kthaFodhKIByM+BEJcRHz4lF53+/Y\nVqRZMNyoA2LjCLy8lo5DSL0QDZO6FbNicxIW8LeqJPiiUozDuA8Diw4HalC7oa2TgOXdsAtS9nZw\n8fL3GA5eFCutAM8TvQrPdBbXoxdGMke8gTQ51SaZBMchcF6IhukjVOqgE3OVsXk7mCIWXuxrcreD\nluI88sBXh1og+wX90v8DMz2dBCw2pyBlJ8HFUzmEGwcf+3Yg3Czx+UIDojvpVdTqm9KCZYI2JQpT\nrkOAvBANqIgfIXASxFxlbD4YTREr92A8SD7XG6bn2FtEV5DPlUOaDnn3yu2wAbm2C1J2Ely8/FGk\ng8eHdqjc8tfRQo92Uv19qAnOBk2oYaHJrTbJleA4BNAL0YB2nLfhEXRyrtIJP9SaIpZVyEpfCYWL\nXmI/8dLhBvnBC3f16uO7VRKwqkyBiW3sJLh4efvr4ON3jP76bQGV+1Ta/7M2LvH0gTYcwkKTW22S\nSXAcAuiFaEAPjMf/xiqVc5WxG12riKVvdABNTwI7YcKhExF+qll64r7ACORK0De1O7TcHIxx/N/o\n78A4+Lu1UVavLRIYmQSt1hPwwgSQSMa2y06lbgTLToK1ZEj8VLPSikw9sVlkjq0HqOsUIo4v7r30\naahsvNxckcd0TlsktJcIiVbrJzhTEbPwKrZddlr1A5TzKCx+qjnbMa5LT7wV8uAZIfF7orOaOB6P\nxyustvD1voWfzEqLhD6vImYZ6NiqErMZKv7TlL0Sqi2NwzM+/7Tzr9FYc/KVR6LQsxp8OdBOMnsZ\nYSJk0gRSDEpDIeSs2plFEuyf+scM7esrfqq5LbnYJr3dWgX+xVi4CUMZYYIS5/IyqyZWmhy7W3lU\n7ziFLPFTf+efPhwTvuEwJH6qeUD8isLIsph8vLouWU5OhU7800dWIG9ktOKnmsduNOM8Jh9fa+f5\n+L/+VGbgf03lK9UAMrrjAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}c \\left(- \\frac{d}{d t} \\operatorname{x_{1}}{\\left (t \\right )} + \\frac{d}{d t} \\operatorname{x_{2}}{\\left (t \\right )}\\right) - 0.5 k \\left(2 \\operatorname{x_{1}}{\\left (t \\right )} - 2 \\operatorname{x_{2}}{\\left (t \\right )}\\right) + f{\\left (t \\right )}\\\\- c \\left(- \\frac{d}{d t} \\operatorname{x_{1}}{\\left (t \\right )} + \\frac{d}{d t} \\operatorname{x_{2}}{\\left (t \\right )}\\right) - 0.5 k \\left(- 2 \\operatorname{x_{1}}{\\left (t \\right )} + 2 \\operatorname{x_{2}}{\\left (t \\right )}\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡  ⎛  d           d        ⎞                                   ⎤\n",
       "⎢c⋅⎜- ──(x₁(t)) + ──(x₂(t))⎟ - 0.5⋅k⋅(2⋅x₁(t) - 2⋅x₂(t)) + f(t)⎥\n",
       "⎢  ⎝  dt          dt       ⎠                                   ⎥\n",
       "⎢                                                              ⎥\n",
       "⎢      ⎛  d           d        ⎞                               ⎥\n",
       "⎢  - c⋅⎜- ──(x₁(t)) + ──(x₂(t))⎟ - 0.5⋅k⋅(-2⋅x₁(t) + 2⋅x₂(t))  ⎥\n",
       "⎣      ⎝  dt          dt       ⎠                               ⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Output the forcing terms of the system\n",
    "LM.forcing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border: 0px;\n",
    "        height: 1px;\n",
    "        text-align: center;\n",
    "        background: #333;\n",
    "        background-image: -webkit-linear-gradient(left, #ccc, #333, #ccc); \n",
    "        background-image:    -moz-linear-gradient(left, #ccc, #333, #ccc); \n",
    "        background-image:     -ms-linear-gradient(left, #ccc, #333, #ccc); \n",
    "        background-image:      -o-linear-gradient(left, #ccc, #333, #ccc);\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Licenses\n",
    "Code is licensed under a 3-clause BSD style license. See the licenses/LICENSE.md file.\n",
    "\n",
    "Other content is provided under a [Creative Commons Attribution-NonCommercial 4.0 International License](http://creativecommons.org/licenses/by-nc/4.0/), CC-BY-NC 4.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n",
       "\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsx.otf');\n",
       "        font-weight: bold;\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsi.otf');\n",
       "        font-style: oblique;\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunso.otf');\n",
       "        font-weight: bold;\n",
       "        font-style: oblique;\n",
       "    }\n",
       "\n",
       "    div.cell {\n",
       "        max-width: 1100px;\n",
       "    }\n",
       "    \n",
       "    h1 {\n",
       "        font-family: Computer Modern;\n",
       "    }\n",
       "    \n",
       "    h4 {\n",
       "        margin-top: 12px;\n",
       "        margin-bottom: 3px;\n",
       "    }\n",
       "\n",
       "    div.text_cell_render {\n",
       "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 145%;\n",
       "        font-size: 130%;\n",
       "        width: 100%;\n",
       "        max-width: 1100px;\n",
       "    }\n",
       "    \n",
       "    .CodeMirror {\n",
       "        font-family: \"Source Code Pro\", source-code-pro, Consolas, monospace;\n",
       "    }\n",
       "    \n",
       "    .warning {\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    }  \n",
       "  \n",
       "   \n",
       "    hr.style-end {\n",
       "        border: 0px !important;\n",
       "        height: 1px !important;\n",
       "        text-align: center !important;\n",
       "        background: #333 !important;\n",
       "        background-image: -webkit-linear-gradient(left, #ccc, #333, #ccc) !important; \n",
       "        background-image:    -moz-linear-gradient(left, #ccc, #333, #ccc) !important; \n",
       "        background-image:     -ms-linear-gradient(left, #ccc, #333, #ccc) !important; \n",
       "        background-image:      -o-linear-gradient(left, #ccc, #333, #ccc) !important; \n",
       "    }\n",
       "\n",
       "    hr.style-end:after {\n",
       "        content: &#x269C !important;\n",
       "        left: 50% !important;\n",
       "        position: absolute !important;\n",
       "        /* Controls the whitespace around the symbol */\n",
       "        padding: 0px !important;\n",
       "        background: #fff !important;\n",
       "    }\n",
       "    \n",
       "/*  Center figures, etc\n",
       "    .ui-wrapper {\n",
       "        margin-left: auto !important;\n",
       "        margin-right: auto !important;\n",
       "    }\n",
       "*/\n",
       "    \n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This cell will just improve the styling of the notebook\n",
    "# You can ignore it, if you are okay with the default sytling\n",
    "from IPython.core.display import HTML\n",
    "import urllib.request\n",
    "response = urllib.request.urlopen(\"https://cl.ly/1B1y452Z1d35\")\n",
    "HTML(response.read().decode(\"utf-8\"))"
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}