Compound_Double_Pendulum/Derivation_Compound.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cfe0c1e811986855",
   "metadata": {},
   "source": [
    "# Double Pendulum Equations of Motion  \n",
    "## Lagrangian Formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1148ed921a221a24",
   "metadata": {},
   "source": [
    "![](Image_files/Massive_rods.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8111a08ca5624fc3",
   "metadata": {},
   "source": [
    "#### Compound Double Pendulum System: \n",
    "Figure shows a planar, double compound pendulum hinged at $O$, represented by two mass-carrying uniform rods, connected in series via a frictionless hinge and allowed to move in the $(x,y)$-plane.\n",
    "\n",
    "- **Massive Rods**: The rods $OP_1$ and $P_1P_2$ possess non-zero masses $M_1$ and $M_2$.\n",
    "  - $OP_1$ performs rotation about the fixed axis\n",
    "  - $P_1P_2$ is in plane motion. \n",
    "- **Moment of Inertia**: Signifying resistance to rotational acceleration about the $_{\\text{CM}}$.\n",
    "- **Centres of Mass**: Labelled $C_1$ and $C_2$\n",
    "\n",
    "The system's state is uniquely determined by two rotational angles $\\theta_1$ and $\\theta_2$, offering two degrees of freedom."
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "The parallel axis theorem is used to determine the moment of inertia of a body about any axis, given the moment of inertia of the body about a parallel axis that passes through its center of mass and the perpendicular distance between the two axes.\n",
    "\n",
    "1. Calculate the moment of inertia $I_{\\text{CM}}$ about the axis through the center of mass using the standard formula for the given geometry of the body.\n",
    "2. Measure the distance $d$ from the $\\text{CM}$ to the new axis about which you want to find the moment of inertia.\n",
    "3. Apply the theorem: $I = I_{\\text{CM}} + ML^2$, where:\n",
    "   - $I$ is the moment of inertia about the new axis,\n",
    "   - $I_{\\text{CM}}$ is the moment of inertia about the $\\text{CM}$ axis,\n",
    "   - $M$ is the mass of the body,\n",
    "   - $L$ is the distance between the two parallel axes.\n",
    "\n",
    "This allows the calculation of the rotational kinetic energy $T_{\\text{rot}} = \\frac{1}{2}I\\omega^2$, where $\\omega$ is the angular velocity of the body about the new axis.\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "e24abd269c2f21fe"
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.702730Z",
     "start_time": "2023-12-09T22:36:23.604408Z"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ec116cf8491fab47",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.705399Z",
     "start_time": "2023-12-09T22:36:23.611891Z"
    }
   },
   "outputs": [],
   "source": [
    "# Declare variables & constants\n",
    "t = sp.Symbol(\"t\") \n",
    "l1, l2, M1, M2, g = sp.symbols('l1 l2 M1 M2 g', real=True, positive=True)\n",
    "R1, R2 = sp.symbols('R1 R2', real=True, positive=True)\n",
    "theta1 = sp.Function('theta1')(t)\n",
    "theta2 = sp.Function('theta2')(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "omega1 = sp.diff(theta1, t)\n",
    "omega2 = sp.diff(theta2, t)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.711816Z",
     "start_time": "2023-12-09T22:36:23.616374Z"
    }
   },
   "id": "c270b762820eab30"
  },
  {
   "cell_type": "markdown",
   "id": "5d90fd818fc97833",
   "metadata": {},
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Define Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "29cabcdf55923881",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.735111Z",
     "start_time": "2023-12-09T22:36:23.623346Z"
    }
   },
   "outputs": [],
   "source": [
    "def trans_kinetic_energy(m, dx, dy, dz):\n",
    "    T = sp.Rational(1,2)*m*(dx**2 + dy**2 + dz**2)\n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [],
   "source": [
    "def potential_energy(m, g, h):\n",
    "    V = m*g*h\n",
    "    return V"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.736385Z",
     "start_time": "2023-12-09T22:36:23.629583Z"
    }
   },
   "id": "da27596ecf7ea1c"
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "def moment_of_inertia(M, L, model='uniform', R=None):\n",
    "    \"\"\"\n",
    "    Calculate the moment of inertia for a rod based on the specified model.\n",
    "    Can be 'uniform' for a uniform rod or 'cylindrical' for a cylindrical rod.\n",
    "    The default is 'uniform'.\n",
    "        \n",
    "    Returns:\n",
    "    The moment of inertia of the rod based on the specified model.\n",
    "        \n",
    "    Raises:\n",
    "    ValueError\n",
    "        If the 'cylindrical' model is selected and no radius R is provided.  \n",
    "    \"\"\"\n",
    "    if model == 'uniform':\n",
    "        # The moment of inertia for a thin rod\n",
    "        I_cm = sp.Rational(1, 12) * M * L**2\n",
    "        # Applying parallel axis theorem\n",
    "        I_end = I_cm + M * (L/2)**2\n",
    "        return I_end\n",
    "    elif model == 'cylindrical':\n",
    "        if R is None:\n",
    "            raise ValueError(\"Radius R must be provided for the cylindrical model.\")\n",
    "        # The moment of inertia for a solid cylinder\n",
    "        I_cm = sp.Rational(1, 4) * M * R**2 + sp.Rational(1, 12) * M * L**2\n",
    "        # Applying the parallel axis theorem \n",
    "        I_end = I_cm + M * (L / 2)**2\n",
    "        return I_end"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.737323Z",
     "start_time": "2023-12-09T22:36:23.635115Z"
    }
   },
   "id": "c46e060900bcc161"
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "def rotational_kinetic_energy(M, L, omega, model='uniform', R=None):\n",
    "    I = moment_of_inertia(M, L, model, R) \n",
    "    T_rot = sp.Rational(1, 2) * I * omega**2\n",
    "    return T_rot"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.737402Z",
     "start_time": "2023-12-09T22:36:23.639456Z"
    }
   },
   "id": "f29e41955f863fc9"
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "def form_lagrangian(theta1, theta2, l1, l2, M1, M2, g, R1=None, R2=None, model='uniform'):\n",
    "    \"\"\"\n",
    "    The function computes the Lagrangian by determining the translational and rotational kinetic energies, \n",
    "    as well as the potential energy, of two rods in a pendulum system. \n",
    "    It supports 'uniform' and 'cylindrical' models. (currently)\n",
    "\n",
    "    Returns the symbolic expression for the Lagrangian of the system.\n",
    "    \"\"\"\n",
    "    # Calculate positions of the center of mass\n",
    "    x1 = l1 / 2 * sp.sin(theta1)\n",
    "    y1 = -l1 / 2 * sp.cos(theta1)\n",
    "    x2 = x1 + l2 / 2 * sp.sin(theta2)\n",
    "    y2 = y1 - l2 / 2 * sp.cos(theta2)\n",
    "\n",
    "    # Calculate velocities\n",
    "    xdot1 = sp.diff(x1, t)\n",
    "    ydot1 = sp.diff(y1, t)\n",
    "    xdot2 = sp.diff(x2, t)\n",
    "    ydot2 = sp.diff(y2, t)\n",
    "    \n",
    "    # Define angular velocity\n",
    "    omega_1 = sp.diff(theta1, t)\n",
    "    omega_2 = sp.diff(theta2, t)\n",
    "\n",
    "    # Calculate energies using the previously defined functions\n",
    "    T1_trans = trans_kinetic_energy(M1, xdot1, ydot1, 0)\n",
    "    V1 = potential_energy(M1, g, y1)\n",
    "    T1_rot = rotational_kinetic_energy(M1, l1, omega_1, model, R1)\n",
    "\n",
    "    T2_trans = trans_kinetic_energy(M2, xdot2, ydot2, 0)\n",
    "    V2 = potential_energy(M2, g, y2)\n",
    "    T2_rot = rotational_kinetic_energy(M2, l2, omega_2, model, R2)\n",
    "\n",
    "    # Form the Lagrangian\n",
    "    T = sp.simplify(T1_trans + T1_rot + T2_trans + T2_rot)\n",
    "    V = sp.simplify(V1 + V2)\n",
    "    L = T - V\n",
    "\n",
    "    return L"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:23.749002Z",
     "start_time": "2023-12-09T22:36:23.649086Z"
    }
   },
   "id": "f1997913c2bd0f05"
  },
  {
   "cell_type": "markdown",
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Form Lagrangian"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "fe5d809cc0534789"
  },
  {
   "cell_type": "markdown",
   "source": [
    "`Uniform rod` model"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "aecedf39e29ceffb"
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "7*M1*l1**2*Derivative(theta1(t), t)**2/24 + M2*l1**2*Derivative(theta1(t), t)**2/8 + M2*l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta1(t), t)*Derivative(theta2(t), t)/4 + 7*M2*l2**2*Derivative(theta2(t), t)**2/24 + g*(M1*l1*cos(theta1(t)) + M2*l1*cos(theta1(t)) + M2*l2*cos(theta2(t)))/2",
      "text/latex": "$\\displaystyle \\frac{7 M_{1} l_{1}^{2} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{24} + \\frac{M_{2} l_{1}^{2} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{8} + \\frac{M_{2} l_{1} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta_{1}{\\left(t \\right)} \\frac{d}{d t} \\theta_{2}{\\left(t \\right)}}{4} + \\frac{7 M_{2} l_{2}^{2} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{24} + \\frac{g \\left(M_{1} l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\cos{\\left(\\theta_{2}{\\left(t \\right)} \\right)}\\right)}{2}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L_uniform = form_lagrangian(theta1, theta2, l1, l2, M1, M2, g, model='uniform')\n",
    "display(L_uniform)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:24.898013Z",
     "start_time": "2023-12-09T22:36:23.668979Z"
    }
   },
   "id": "9e0c97e37d81f40f"
  },
  {
   "cell_type": "markdown",
   "source": [
    "`cylindrical_rod` model"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "3e38490403e9fa05"
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d7ee39c49688e33b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:26.099749Z",
     "start_time": "2023-12-09T22:36:24.910749Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "M1*R1**2*Derivative(theta1(t), t)**2/8 + 7*M1*l1**2*Derivative(theta1(t), t)**2/24 + M2*R2**2*Derivative(theta2(t), t)**2/8 + M2*l1**2*Derivative(theta1(t), t)**2/8 + M2*l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta1(t), t)*Derivative(theta2(t), t)/4 + 7*M2*l2**2*Derivative(theta2(t), t)**2/24 + g*(M1*l1*cos(theta1(t)) + M2*l1*cos(theta1(t)) + M2*l2*cos(theta2(t)))/2",
      "text/latex": "$\\displaystyle \\frac{M_{1} R_{1}^{2} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{8} + \\frac{7 M_{1} l_{1}^{2} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{24} + \\frac{M_{2} R_{2}^{2} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{8} + \\frac{M_{2} l_{1}^{2} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{8} + \\frac{M_{2} l_{1} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta_{1}{\\left(t \\right)} \\frac{d}{d t} \\theta_{2}{\\left(t \\right)}}{4} + \\frac{7 M_{2} l_{2}^{2} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{24} + \\frac{g \\left(M_{1} l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\cos{\\left(\\theta_{2}{\\left(t \\right)} \\right)}\\right)}{2}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L_cylindrical = form_lagrangian(theta1, theta2, l1, l2, M1, M2, g, R1, R2, model='cylindrical')\n",
    "display(L_cylindrical)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "462e5a3966d7acf0",
   "metadata": {},
   "source": [
    "----\n",
    "&nbsp;\n",
    "## System of Euler-Lagrange Equations: \n",
    "For each generalised coordinate $\\theta_i$\n",
    "\n",
    "$$\\frac{\\text{d}}{\\text{d}t}\\left(\\frac{\\partial L}{\\partial \\dot{\\theta}_i}\\right)-\\frac{\\partial L}{\\partial \\theta_i}=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [],
   "source": [
    "def euler_lagrange_system(L, q1, q2, model='uniform'):\n",
    "    \"\"\"\n",
    "    Computes the Euler-Lagrange equations for a double pendulum system with two generalized coordinates.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    q1 : sympy.Function\n",
    "        The first generalized coordinate as a function of time.\n",
    "    q2 : sympy.Function\n",
    "        The second generalized coordinate as a function of time.\n",
    "    model : str, optional\n",
    "        The model type of the pendulum system, either 'uniform' or 'cylindrical'. Default is 'uniform'.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    tuple of sympy.Eq\n",
    "        A tuple containing two sympy.Eq objects representing the simplified Euler-Lagrange equations\n",
    "        for the two generalized coordinates.\n",
    "\n",
    "    Raises\n",
    "    ------\n",
    "    ValueError\n",
    "        If an invalid model type is provided.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    The function assumes that the Lagrangians (L_uniform and L_cylindrical) for both model types\n",
    "    are precomputed and accessible in the scope of this function.\n",
    "    \"\"\"\n",
    "    # Derivatives of the generalized coordinates\n",
    "    q1_dot = sp.diff(q1, t)\n",
    "    q2_dot = sp.diff(q2, t)\n",
    "\n",
    "    # Euler-Lagrange equations for q1\n",
    "    partial_L_partial_q1 = sp.diff(L, q1)\n",
    "    partial_L_partial_q1_dot = sp.diff(L, q1_dot)\n",
    "    EL_q1 = sp.diff(partial_L_partial_q1_dot, t) - partial_L_partial_q1\n",
    "\n",
    "    # Euler-Lagrange equations for q2\n",
    "    partial_L_partial_q2 = sp.diff(L, q2)\n",
    "    partial_L_partial_q2_dot = sp.diff(L, q2_dot)\n",
    "    EL_q2 = sp.diff(partial_L_partial_q2_dot, t) - partial_L_partial_q2\n",
    "\n",
    "    # Apply specific simplification factors based on the model\n",
    "    if model == 'uniform':\n",
    "        EL_q1 = EL_q1 * (12 / l1)\n",
    "        EL_q2 = EL_q2 * (12 / (M2 * l2))\n",
    "    elif model == 'cylindrical':\n",
    "        EL_q1 = EL_q1 * 2\n",
    "        EL_q2 = EL_q2 * (4 / M2)\n",
    "\n",
    "    # Simplify the equations\n",
    "    EL_q1_simplified = sp.simplify(EL_q1)\n",
    "    EL_q2_simplified = sp.simplify(EL_q2)\n",
    "    \n",
    "    eqn1 = sp.Eq(EL_q1_simplified, 0)\n",
    "    eqn2 = sp.Eq(EL_q2_simplified, 0)\n",
    "\n",
    "    return eqn1, eqn2"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:26.105597Z",
     "start_time": "2023-12-09T22:36:26.100791Z"
    }
   },
   "id": "cbc9c8ce94067d07"
  },
  {
   "cell_type": "markdown",
   "source": [
    "`uniform_rod` model"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "193efbef3ca7fc46"
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn1\n",
      "<class 'sympy.core.relational.Equality'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq(6*M1*g*sin(theta1(t)) + 7*M1*l1*Derivative(theta1(t), (t, 2)) + 6*M2*g*sin(theta1(t)) + 3*M2*l1*Derivative(theta1(t), (t, 2)) + 3*M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2 + 3*M2*l2*cos(theta1(t) - theta2(t))*Derivative(theta2(t), (t, 2)), 0)",
      "text/latex": "$\\displaystyle 6 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 7 M_{1} l_{1} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)} + 6 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 3 M_{2} l_{1} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)} + 3 M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2} + 3 M_{2} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)} = 0$"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn2\n",
      "<class 'sympy.core.relational.Equality'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq(6*g*sin(theta2(t)) - 3*l1*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2 + 3*l1*cos(theta1(t) - theta2(t))*Derivative(theta1(t), (t, 2)) + 7*l2*Derivative(theta2(t), (t, 2)), 0)",
      "text/latex": "$\\displaystyle 6 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - 3 l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2} + 3 l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)} + 7 l_{2} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)} = 0$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eqn1_u, eqn2_u = euler_lagrange_system(L_uniform, theta1, theta2, model='uniform')\n",
    "uniform_system = [eqn1_u, eqn2_u]\n",
    "for i, eqn in enumerate(uniform_system):\n",
    "    print(f\"eqn{i+1}\\n{type(eqn)}\")\n",
    "    display(eqn)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:26.748487Z",
     "start_time": "2023-12-09T22:36:26.109785Z"
    }
   },
   "id": "fdebce9ebad46cd5"
  },
  {
   "cell_type": "markdown",
   "source": [
    "`cylindrical_rod` model"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "22632ed00ddf8cad"
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn1\n",
      "<class 'sympy.core.relational.Equality'>\n",
      "\n",
      "Eq(M1*R1**2*Derivative(theta1(t), (t, 2))/2 + M1*g*l1*sin(theta1(t)) + 7*M1*l1**2*Derivative(theta1(t), (t, 2))/6 + M2*g*l1*sin(theta1(t)) + M2*l1**2*Derivative(theta1(t), (t, 2))/2 + M2*l1*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2/2 + M2*l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta2(t), (t, 2))/2, 0)\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq(M1*R1**2*Derivative(theta1(t), (t, 2))/2 + M1*g*l1*sin(theta1(t)) + 7*M1*l1**2*Derivative(theta1(t), (t, 2))/6 + M2*g*l1*sin(theta1(t)) + M2*l1**2*Derivative(theta1(t), (t, 2))/2 + M2*l1*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2/2 + M2*l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta2(t), (t, 2))/2, 0)",
      "text/latex": "$\\displaystyle \\frac{M_{1} R_{1}^{2} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)}}{2} + M_{1} g l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + \\frac{7 M_{1} l_{1}^{2} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)}}{6} + M_{2} g l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + \\frac{M_{2} l_{1}^{2} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)}}{2} + \\frac{M_{2} l_{1} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{2} + \\frac{M_{2} l_{1} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)}}{2} = 0$"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn2\n",
      "<class 'sympy.core.relational.Equality'>\n",
      "\n",
      "Eq(R2**2*Derivative(theta2(t), (t, 2)) + 2*g*l2*sin(theta2(t)) - l1*l2*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2 + l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta1(t), (t, 2)) + 7*l2**2*Derivative(theta2(t), (t, 2))/3, 0)\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq(R2**2*Derivative(theta2(t), (t, 2)) + 2*g*l2*sin(theta2(t)) - l1*l2*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2 + l1*l2*cos(theta1(t) - theta2(t))*Derivative(theta1(t), (t, 2)) + 7*l2**2*Derivative(theta2(t), (t, 2))/3, 0)",
      "text/latex": "$\\displaystyle R_{2}^{2} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)} + 2 g l_{2} \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2} + l_{1} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)} + \\frac{7 l_{2}^{2} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)}}{3} = 0$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eqn1_c, eqn2_c = euler_lagrange_system(L_cylindrical, theta1, theta2, model='cylindrical')\n",
    "cylindrical_system = [eqn1_c, eqn2_c]\n",
    "for i, eqn in enumerate(cylindrical_system):\n",
    "    print(f\"eqn{i+1}\\n{type(eqn)}\\n\")\n",
    "    print(eqn)\n",
    "    display(eqn)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:27.664685Z",
     "start_time": "2023-12-09T22:36:26.797681Z"
    }
   },
   "id": "bc3310a17100b8fe"
  },
  {
   "cell_type": "markdown",
   "id": "63780e157ab7d634",
   "metadata": {},
   "source": [
    "#### Which are our coupled 2nd-order differential equations of motion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1740e63f7ad574",
   "metadata": {},
   "source": [
    "----\n",
    "&nbsp;\n",
    "#### Rewrite E-L equations in form;\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\ddot{\\theta}_1+\\alpha_1\\left(\\theta_1, \\theta_2\\right) \\ddot{\\theta}_2=f_1\\left(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2\\right) \\\\\n",
    "& \\ddot{\\theta}_2+\\alpha_2\\left(\\theta_1, \\theta_2\\right) \\ddot{\\theta}_1=f_2\\left(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2\\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "Extract equations from system to avoid name space conflicts"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "23d389cdf9507af8"
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [
    "system = uniform_system\n",
    "eq1, eq2 = system[0], system[1]"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:46.338454Z",
     "start_time": "2023-12-09T22:36:46.328528Z"
    }
   },
   "id": "21ad906ddcfe1df"
  },
  {
   "cell_type": "markdown",
   "id": "94b6ca4fe1870604",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T14:36:40.887790Z",
     "start_time": "2023-12-09T14:36:40.864639Z"
    }
   },
   "source": [
    "#### Isolate the second derivative coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [],
   "source": [
    "def isolate_terms(equation):\n",
    "    \"\"\"\n",
    "    Returns the second derivative terms isolated from the given equation.\n",
    "    \n",
    "    Parameters:\n",
    "        equation (sympy.Eq): The equation from which to isolate the second derivative terms.\n",
    "    \n",
    "    Returns:\n",
    "        list: A list containing the isolated second derivative terms.\n",
    "    \"\"\"\n",
    "    # Define derivative terms\n",
    "    theta1ddot = sp.diff(theta1, t, 2)\n",
    "    theta2ddot = sp.diff(theta2, t, 2)\n",
    "    \n",
    "    terms = []\n",
    "\n",
    "    try:\n",
    "        th1 = sp.Eq(theta1ddot, sp.solve(equation, theta1ddot)[0])\n",
    "        rhs_eq = th1.rhs\n",
    "        alpha1 = sp.together(rhs_eq).as_numer_denom()[1]\n",
    "        eq_new = sp.Eq(th1.lhs * alpha1, th1.rhs * alpha1)\n",
    "        theta1_second = eq_new.lhs\n",
    "        terms.append(theta1_second)\n",
    "    except IndexError:\n",
    "        pass\n",
    "\n",
    "    try:\n",
    "        th2 = sp.Eq(theta2ddot, sp.solve(equation, theta2ddot)[0])\n",
    "        rhs_eq = th2.rhs\n",
    "        alpha2 = sp.together(rhs_eq).as_numer_denom()[1]\n",
    "        eq_new = sp.Eq(th2.lhs * alpha2, th2.rhs * alpha2)\n",
    "        theta2_second = eq_new.lhs\n",
    "        terms.append(theta2_second)\n",
    "    except IndexError:\n",
    "        pass\n",
    "\n",
    "    return terms"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:50.043883Z",
     "start_time": "2023-12-09T22:36:50.037990Z"
    }
   },
   "id": "603798eb2795e9c"
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [],
   "source": [
    "def simplify_system_equations(eq1, eq2, model='uniform'):\n",
    "    \"\"\"\n",
    "    Simplifies the system of equations based on the specified model.\n",
    "    \n",
    "    Parameters:\n",
    "    eq1, eq2 : sympy.Eq\n",
    "        The Euler-Lagrange equations of the system.\n",
    "    model : str, optional\n",
    "        The model type of the pendulum system ('uniform' or 'cylindrical'). Default is 'uniform'.\n",
    "    \n",
    "    Returns:\n",
    "    tuple of sympy.Eq\n",
    "        The simplified equations for the system.\n",
    "    \"\"\"\n",
    "    if model == 'uniform':\n",
    "        # Isolate second derivatives and simplify equations for eq1\n",
    "        second_derivatives_1 = isolate_terms(eq1)\n",
    "        lhs_1 = sum(second_derivatives_1)\n",
    "        rhs_1 = sp.simplify(eq1.lhs - lhs_1) \n",
    "        eqn1 = sp.Eq(lhs_1, rhs_1)\n",
    "\n",
    "        # Isolate second derivatives and simplify equations for eq2\n",
    "        second_derivatives_2 = isolate_terms(eq2)\n",
    "        lhs_2 = sum(second_derivatives_2)\n",
    "        rhs_2 = sp.simplify(eq2.lhs - lhs_2)\n",
    "        eqn2 = sp.Eq(lhs_2, rhs_2)\n",
    "\n",
    "        # Simplify eqn1\n",
    "        eqn1_lhs = sp.simplify(eqn1.lhs / (l1 * (7*M1 + 3*M2)))\n",
    "        eqn1_rhs = sp.simplify(eqn1.rhs / (l1 * (7*M1 + 3*M2)))\n",
    "        eqn1 = sp.Eq(eqn1_lhs, eqn1_rhs)\n",
    "\n",
    "        # Simplify eqn2\n",
    "        eqn2_lhs = sp.simplify(eqn2.lhs / (7 * l2))\n",
    "        eqn2_rhs = sp.simplify(eqn2.rhs / (7 * l2))\n",
    "        eqn2 = sp.Eq(eqn2_lhs, eqn2_rhs)\n",
    "\n",
    "        return eqn1, eqn2\n",
    "    \n",
    "    elif model == 'cylindrical':\n",
    "        # Define second-order derivative terms\n",
    "        theta1_ddot = sp.diff(theta1, t, 2)\n",
    "        theta2_ddot = sp.diff(theta2, t, 2)\n",
    "\n",
    "        # Collect both second derivatives in eq1 and isolate other terms\n",
    "        eq1_lhs_collected = sp.collect(eq1.lhs, (theta1_ddot, theta2_ddot))\n",
    "        eqn1_rhs = sp.sympify(M1*g*l1*sp.sin(theta1) + M2*g*l1*sp.sin(theta1) + M2*l1*l2*sp.sin(theta1 - theta2)*sp.diff(theta2, t)**2/2)\n",
    "        eqn1_lhs = eq1_lhs_collected - eqn1_rhs\n",
    "        eq1_simp = sp.Eq(eqn1_lhs, -eqn1_rhs)\n",
    "\n",
    "        # Collect both second derivatives in eq2 and isolate other terms\n",
    "        eq2_lhs_collected = sp.collect(eq2.lhs, (theta1_ddot, theta2_ddot))\n",
    "        eqn2_rhs = sp.sympify(2*g*l2*sp.sin(theta2) - l1*l2*sp.sin(theta1 - theta2)*sp.diff(theta1, t)**2)\n",
    "        eqn2_lhs = eq2_lhs_collected - eqn2_rhs\n",
    "        eq2_simp = sp.Eq(eqn2_lhs, eqn2_rhs)\n",
    "        \n",
    "        # Simplify further\n",
    "        eqn1_lhs = sp.simplify(eq1_simp.lhs / (M1*R1**2/2 + 7*M1*l1**2/6 + M2*l1**2/2))\n",
    "        eqn1_rhs = sp.simplify(eq1_simp.rhs / (M1*R1**2/2 + 7*M1*l1**2/6 + M2*l1**2/2))\n",
    "        eqn1 = sp.Eq(eqn1_lhs, eqn1_rhs)\n",
    "       \n",
    "        eqn2_lhs = sp.simplify(eq2_simp.lhs / (R2**2 + 7*l2**2/3))\n",
    "        eqn2_rhs = sp.simplify(eq2_simp.rhs / (R2**2 + 7*l2**2/3))\n",
    "        eqn2 = sp.Eq(eqn2_lhs, eqn2_rhs)\n",
    "\n",
    "        return eqn1, eqn2"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:36:56.596913Z",
     "start_time": "2023-12-09T22:36:56.586932Z"
    }
   },
   "id": "d516dc419ef9e17c"
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [],
   "source": [
    "eqn1, eqn2 = simplify_system_equations(eq1, eq2, model='uniform')"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:08.846787Z",
     "start_time": "2023-12-09T22:37:06.392860Z"
    }
   },
   "id": "583af1da90c5d339"
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn1\n",
      "<class 'sympy.core.relational.Equality'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq((3*M2*l2*cos(theta1(t) - theta2(t))*Derivative(theta2(t), (t, 2)) + l1*(7*M1 + 3*M2)*Derivative(theta1(t), (t, 2)))/(l1*(7*M1 + 3*M2)), 3*(2*M1*g*sin(theta1(t)) + 2*M2*g*sin(theta1(t)) + M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)/(l1*(7*M1 + 3*M2)))",
      "text/latex": "$\\displaystyle \\frac{3 M_{2} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)} + l_{1} \\cdot \\left(7 M_{1} + 3 M_{2}\\right) \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)}}{l_{1} \\cdot \\left(7 M_{1} + 3 M_{2}\\right)} = \\frac{3 \\cdot \\left(2 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 2 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}\\right)}{l_{1} \\cdot \\left(7 M_{1} + 3 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eqn2\n",
      "<class 'sympy.core.relational.Equality'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "Eq(3*l1*cos(theta1(t) - theta2(t))*Derivative(theta1(t), (t, 2))/(7*l2) + Derivative(theta2(t), (t, 2)), 3*(2*g*sin(theta2(t)) - l1*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2)/(7*l2))",
      "text/latex": "$\\displaystyle \\frac{3 l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\frac{d^{2}}{d t^{2}} \\theta_{1}{\\left(t \\right)}}{7 l_{2}} + \\frac{d^{2}}{d t^{2}} \\theta_{2}{\\left(t \\right)} = \\frac{3 \\cdot \\left(2 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}\\right)}{7 l_{2}}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "system = [eqn1, eqn2]\n",
    "for i, eqn in enumerate(system):\n",
    "    print(f\"eqn{i+1}\\n{type(eqn)}\")\n",
    "    display(eqn)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:11.448268Z",
     "start_time": "2023-12-09T22:37:11.421212Z"
    }
   },
   "id": "a5f15b7c3f98a6bd"
  },
  {
   "cell_type": "markdown",
   "id": "800441001cc2740f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T14:44:00.372836Z",
     "start_time": "2023-12-09T14:44:00.352456Z"
    }
   },
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Extract coefficients from above equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "af152af6feb7820d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:24.921172Z",
     "start_time": "2023-12-09T22:37:24.912357Z"
    }
   },
   "outputs": [],
   "source": [
    "def extract_coefficient(equation, derivative_term):\n",
    "    \"\"\"\n",
    "    Extracts the coefficient, including the denominator, of the specified derivative term from the left-hand side of the equation.\n",
    "\n",
    "    Parameters:\n",
    "        equation (sympy.Expr): The equation to extract the coefficient from.\n",
    "        derivative_term (sympy.Expr): The derivative term to extract the coefficient of.\n",
    "\n",
    "    Returns:\n",
    "        coeff (sympy.Expr): The coefficient, including the denominator, of the derivative term in the equation.\n",
    "    \"\"\"\n",
    "    lhs_parts = sp.fraction(equation.lhs)\n",
    "    lhs_coeff_term = lhs_parts[0]\n",
    "    lhs_denominator_term = lhs_parts[1]\n",
    "    coeff = lhs_coeff_term.coeff(derivative_term) / lhs_denominator_term\n",
    "\n",
    "    return coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c89312b57bab41d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T14:44:10.376849Z",
     "start_time": "2023-12-09T14:44:10.363324Z"
    }
   },
   "source": [
    "We can now rewrite in the desired form\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\ddot{\\theta}_1+\\alpha_1\\left(\\theta_1, \\theta_2\\right) \\ddot{\\theta}_2=f_1\\left(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2\\right) \\\\\n",
    "& \\ddot{\\theta}_2+\\alpha_2\\left(\\theta_1, \\theta_2\\right) \\ddot{\\theta}_1=f_2\\left(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2\\right)\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "cca59a287e40eeac",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:27.381922Z",
     "start_time": "2023-12-09T22:37:27.347252Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'sympy.core.mul.Mul'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "3*M2*l2*cos(theta1(t) - theta2(t))/(l1*(7*M1 + 3*M2))",
      "text/latex": "$\\displaystyle \\frac{3 M_{2} l_{2} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}}{l_{1} \\cdot \\left(7 M_{1} + 3 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let alpha1 be 2nd derivative coefficient of theta2 from eqn1\n",
    "alpha1 = extract_coefficient(eqn1, sp.diff(theta2, t, 2))\n",
    "print(type(alpha1))\n",
    "display(alpha1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e9fc75eb0a271c5a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:31.574632Z",
     "start_time": "2023-12-09T22:37:31.565866Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'sympy.core.mul.Mul'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "3*l1*cos(theta1(t) - theta2(t))/(7*l2)",
      "text/latex": "$\\displaystyle \\frac{3 l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}}{7 l_{2}}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let alpha2 be 2nd derivative coefficient of theta1 from eqn2\n",
    "alpha2 = extract_coefficient(eqn2, sp.diff(theta1, t, 2))\n",
    "print(type(alpha2))\n",
    "display(alpha2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90ee182009397850",
   "metadata": {},
   "source": [
    "Declare RHS as expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "33d22d562323e1b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:36.690834Z",
     "start_time": "2023-12-09T22:37:36.684573Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'sympy.core.mul.Mul'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "3*(2*M1*g*sin(theta1(t)) + 2*M2*g*sin(theta1(t)) + M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)/(l1*(7*M1 + 3*M2))",
      "text/latex": "$\\displaystyle \\frac{3 \\cdot \\left(2 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 2 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}\\right)}{l_{1} \\cdot \\left(7 M_{1} + 3 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function_1 = eqn1.rhs\n",
    "print(type(function_1))\n",
    "display(function_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f8ac46d10ad943bf",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:40.343593Z",
     "start_time": "2023-12-09T22:37:40.331996Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'sympy.core.mul.Mul'>\n"
     ]
    },
    {
     "data": {
      "text/plain": "3*(2*g*sin(theta2(t)) - l1*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2)/(7*l2)",
      "text/latex": "$\\displaystyle \\frac{3 \\cdot \\left(2 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}\\right)}{7 l_{2}}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function_2 = eqn2.rhs\n",
    "print(type(function_2))\n",
    "display(function_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Create Matrix Equations"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "1a6bfff8ad1d331f"
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "outputs": [],
   "source": [
    "def create_matrix_equation(theta1dd_coeff, theta2dd_coeff, rhs_1, rhs_2):\n",
    "    \"\"\"\n",
    "    This function constructs a matrix A based on symbolic functions alpha1 and alpha2, which represent\n",
    "    coefficients related to the system's inertia, and solves the equation A*[theta1'', theta2''] = -[f1, f2]\n",
    "    for the angular accelerations.\n",
    "\n",
    "    Parameters:\n",
    "    theta1dd_coeff : symbolic expression\n",
    "        The coefficient for theta1's second derivative, as derived from the system's dynamics.\n",
    "    theta2dd_coeff : symbolic expression\n",
    "        The coefficient for theta2's second derivative, as derived from the system's dynamics.\n",
    "    rhs_1 : symbolic expression\n",
    "        The function representing the external force or torque acting on theta1.\n",
    "    rhs_2 : symbolic expression\n",
    "        The function representing the external force or torque acting on theta2.\n",
    "\n",
    "    Returns:\n",
    "    sympy.Eq\n",
    "        The symbolic equation representing the system, with the specified coefficients and forces/torques substituted in.\n",
    "    \"\"\"\n",
    "    # Define the symbolic functions for alpha1 and alpha2\n",
    "    alpha1 = sp.Function('alpha1')(theta1, theta2)\n",
    "    alpha2 = sp.Function('alpha2')(theta1, theta2)\n",
    "    \n",
    "    # Define the symbolic functions for f1 and f2\n",
    "    function1 = sp.Function('f1')(theta1, theta2, sp.diff(theta1, t), sp.diff(theta2, t))\n",
    "    function2 = sp.Function('f2')(theta1, theta2, sp.diff(theta1, t), sp.diff(theta2, t))\n",
    "    \n",
    "    # Left-hand side of the equation\n",
    "    LHS = sp.Matrix([[sp.diff(theta1, t, 2)], [sp.diff(theta2, t, 2)]])\n",
    "    \n",
    "    # A matrix with symbolic functions alpha1 and alpha2\n",
    "    A = sp.Matrix([[1, alpha1], [alpha2, 1]])\n",
    "    \n",
    "    # Right-hand side of the equation before substitution\n",
    "    RHS = sp.Matrix([[function1], [function2]])\n",
    "    \n",
    "    # Invert matrix A and simplify the product with RHS\n",
    "    NewRHS = sp.simplify(A.inv() * RHS)\n",
    "    \n",
    "    # Create the equation\n",
    "    EQS = sp.Eq(LHS, (-1) * NewRHS)\n",
    "    \n",
    "    # Substitute the derived expressions into the equation\n",
    "    EQS_subst = EQS.subs({\n",
    "        alpha1: theta1dd_coeff, \n",
    "        alpha2: theta2dd_coeff, \n",
    "        function1: rhs_1, \n",
    "        function2: rhs_2\n",
    "    })\n",
    "    \n",
    "    RHS_1 = sp.simplify(EQS_subst.rhs[0])\n",
    "    RHS_2 = sp.simplify(EQS_subst.rhs[1])\n",
    "    \n",
    "    return RHS_1, RHS_2"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:49.741454Z",
     "start_time": "2023-12-09T22:37:49.724743Z"
    }
   },
   "id": "442694a891098612"
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "outputs": [],
   "source": [
    "RHS_1, RHS_2 = create_matrix_equation(alpha1, alpha2, function_1, function_2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:55.256090Z",
     "start_time": "2023-12-09T22:37:52.436430Z"
    }
   },
   "id": "14a3cfd8e69d4b09"
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "outputs": [
    {
     "data": {
      "text/plain": "(42*M1*g*sin(theta1(t)) + 9*M2*g*sin(theta1(t) - 2*theta2(t)) + 33*M2*g*sin(theta1(t)) + 9*M2*l1*sin(2*theta1(t) - 2*theta2(t))*Derivative(theta1(t), t)**2/2 + 21*M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)/(l1*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))",
      "text/latex": "$\\displaystyle \\frac{42 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 9 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)} + 33 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + \\frac{9 M_{2} l_{1} \\sin{\\left(2 \\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{2} + 21 M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{l_{1} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(RHS_1)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:37:57.927071Z",
     "start_time": "2023-12-09T22:37:57.905758Z"
    }
   },
   "id": "8886c2278776ef9c"
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "outputs": [
    {
     "data": {
      "text/plain": "3*((7*M1 + 3*M2)*(2*g*sin(theta2(t)) - l1*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2) - 3*(2*M1*g*sin(theta1(t)) + 2*M2*g*sin(theta1(t)) + M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)*cos(theta1(t) - theta2(t)))/(l2*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))",
      "text/latex": "$\\displaystyle \\frac{3 \\left(\\left(7 M_{1} + 3 M_{2}\\right) \\left(2 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}\\right) - 3 \\cdot \\left(2 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 2 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}\\right) \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}\\right)}{l_{2} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(RHS_2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:01.641178Z",
     "start_time": "2023-12-09T22:38:01.627012Z"
    }
   },
   "id": "93e40b4cf2edb71e"
  },
  {
   "cell_type": "markdown",
   "id": "c1ac15d58954cc73",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T14:46:24.411390Z",
     "start_time": "2023-12-09T14:46:24.395612Z"
    }
   },
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Rewrite as first order system\n",
    "\n",
    "Declare new variables $\\omega_1$ and $\\omega_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "outputs": [],
   "source": [
    "def first_order_system(RHS_1, RHS_2): \n",
    "    \"\"\"\n",
    "    Convert a second-order differential equation system to a first-order system in matrix form.\n",
    "\n",
    "    This function takes the expressions for the right-hand side (RHS) of two second-order \n",
    "    differential equations and constructs a first-order system by introducing new functions \n",
    "    omega1 and omega2, which represent the first derivatives of the original functions \n",
    "    theta1 and theta2, respectively.\n",
    "\n",
    "    Parameters:\n",
    "    RHS_1 : sympy.Expr\n",
    "        The right-hand side expression of the differential equation for d^2(theta1)/dt^2.\n",
    "    RHS_2 : sympy.Expr\n",
    "        The right-hand side expression of the differential equation for d^2(theta2)/dt^2.\n",
    "\n",
    "    Returns:\n",
    "    tuple: \n",
    "        - sympy.Eq representing the matrix equation of the first-order system,\n",
    "        - sympy.Expr for the derivative of theta1 with respect to time after substitution,\n",
    "        - sympy.Expr for the derivative of theta2 with respect to time after substitution,\n",
    "        - sympy.Expr for the derivative of omega1 with respect to time after substitution,\n",
    "        - sympy.Expr for the derivative of omega2 with respect to time after substitution.\n",
    "\n",
    "    The matrix equation (MAT_EQ) equates a column matrix of the left-hand sides of the equations \n",
    "    to a column matrix of the right-hand sides, facilitating the representation and manipulation \n",
    "    of the system of equations.\n",
    "\n",
    "    The function returns this matrix equation along with the isolated expressions for the first \n",
    "    derivatives of theta1 and theta2 (eqn1, eqn2), and the expressions for the second derivatives \n",
    "    of omega1 and omega2 (eqn3, eqn4), which correspond to the second derivatives of theta1 and \n",
    "    theta2, respectively.\n",
    "    \"\"\"\n",
    "    # Define new functions for the first derivatives of theta1 and theta2\n",
    "    omega1 = sp.Function('omega1')(t)\n",
    "    omega2 = sp.Function('omega2')(t)\n",
    "\n",
    "    # Cast as first order system\n",
    "    eq1 = sp.Eq(omega1, sp.diff(theta1, t))\n",
    "    eq2 = sp.Eq(omega2, sp.diff(theta2, t))\n",
    "    eq3 = sp.Eq(sp.diff(omega1, t), RHS_1)\n",
    "    eq4 = sp.Eq(sp.diff(omega2, t), RHS_2)\n",
    "\n",
    "    # Assemble the system into matrix form\n",
    "    LHS_FIRST = sp.Matrix([[eq1.lhs], [eq2.lhs], [eq3.lhs], [eq4.lhs]])\n",
    "    RHS_FIRST = sp.Matrix([[eq1.rhs], [eq2.rhs], [eq3.rhs], [eq4.rhs]])\n",
    "\n",
    "    # Create the matrix equation\n",
    "    MAT_EQ = sp.Eq(LHS_FIRST, RHS_FIRST)\n",
    "\n",
    "    # Substitute omega for derivative of theta\n",
    "    derivative_subs = {sp.Derivative(theta1, t): omega1, sp.Derivative(theta2, t): omega2}\n",
    "\n",
    "    # Isolate rhs of equations\n",
    "    eqn1 = eq1.rhs.subs(derivative_subs)\n",
    "    eqn2 = eq2.rhs.subs(derivative_subs)\n",
    "    eqn3 = eq3.rhs.subs(derivative_subs)\n",
    "    eqn4 = eq4.rhs.subs(derivative_subs)\n",
    "\n",
    "    return MAT_EQ, eqn1, eqn2, eqn3, eqn4"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:08.120565Z",
     "start_time": "2023-12-09T22:38:08.104994Z"
    }
   },
   "id": "3bb2b85e96ae13df"
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "outputs": [],
   "source": [
    "MAT_EQ, eqn1, eqn2, eqn3, eqn4 = first_order_system(RHS_1, RHS_2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:10.724815Z",
     "start_time": "2023-12-09T22:38:10.721589Z"
    }
   },
   "id": "b0b5f1cb353df7ad"
  },
  {
   "cell_type": "markdown",
   "source": [
    "----\n",
    "&nbsp;\n",
    "## Derived Quantities"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "97f79e5758add5ea"
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "outputs": [
    {
     "data": {
      "text/plain": "Eq(Matrix([\n[               omega1(t)],\n[               omega2(t)],\n[Derivative(omega1(t), t)],\n[Derivative(omega2(t), t)]]), Matrix([\n[                                                                                                                                                                                                                                                                                   Derivative(theta1(t), t)],\n[                                                                                                                                                                                                                                                                                   Derivative(theta2(t), t)],\n[                (42*M1*g*sin(theta1(t)) + 9*M2*g*sin(theta1(t) - 2*theta2(t)) + 33*M2*g*sin(theta1(t)) + 9*M2*l1*sin(2*theta1(t) - 2*theta2(t))*Derivative(theta1(t), t)**2/2 + 21*M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)/(l1*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))],\n[3*((7*M1 + 3*M2)*(2*g*sin(theta2(t)) - l1*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)**2) - 3*(2*M1*g*sin(theta1(t)) + 2*M2*g*sin(theta1(t)) + M2*l2*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t)**2)*cos(theta1(t) - theta2(t)))/(l2*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))]]))",
      "text/latex": "$\\displaystyle \\left[\\begin{matrix}\\omega_{1}{\\left(t \\right)}\\\\\\omega_{2}{\\left(t \\right)}\\\\\\frac{d}{d t} \\omega_{1}{\\left(t \\right)}\\\\\\frac{d}{d t} \\omega_{2}{\\left(t \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\\\\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\\\\\frac{42 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 9 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)} + 33 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + \\frac{9 M_{2} l_{1} \\sin{\\left(2 \\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}}{2} + 21 M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}}{l_{1} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}\\\\\\frac{3 \\left(\\left(7 M_{1} + 3 M_{2}\\right) \\left(2 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}\\right) - 3 \\cdot \\left(2 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 2 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}\\right) \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}\\right)}{l_{2} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}\\end{matrix}\\right]$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(MAT_EQ)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:13.292654Z",
     "start_time": "2023-12-09T22:38:13.287748Z"
    }
   },
   "id": "7ec29e864868c6f3"
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Right hand sides have been substituted such that;\n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d} t}\\theta_i\\equiv\\omega_i\n",
    "$$"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "3d784d5f14eed864"
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [
    {
     "data": {
      "text/plain": "(42*M1*g*sin(theta1(t)) + 9*M2*g*sin(theta1(t) - 2*theta2(t)) + 33*M2*g*sin(theta1(t)) + 9*M2*l1*omega1(t)**2*sin(2*theta1(t) - 2*theta2(t))/2 + 21*M2*l2*omega2(t)**2*sin(theta1(t) - theta2(t)))/(l1*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))",
      "text/latex": "$\\displaystyle \\frac{42 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 9 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)} + 33 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + \\frac{9 M_{2} l_{1} \\omega_{1}^{2}{\\left(t \\right)} \\sin{\\left(2 \\theta_{1}{\\left(t \\right)} - 2 \\theta_{2}{\\left(t \\right)} \\right)}}{2} + 21 M_{2} l_{2} \\omega_{2}^{2}{\\left(t \\right)} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}}{l_{1} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(eqn3)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:17.462714Z",
     "start_time": "2023-12-09T22:38:17.437264Z"
    }
   },
   "id": "4eb71e437c3983d2"
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "data": {
      "text/plain": "3*((7*M1 + 3*M2)*(2*g*sin(theta2(t)) - l1*omega1(t)**2*sin(theta1(t) - theta2(t))) - 3*(2*M1*g*sin(theta1(t)) + 2*M2*g*sin(theta1(t)) + M2*l2*omega2(t)**2*sin(theta1(t) - theta2(t)))*cos(theta1(t) - theta2(t)))/(l2*(-49*M1 + 9*M2*cos(theta1(t) - theta2(t))**2 - 21*M2))",
      "text/latex": "$\\displaystyle \\frac{3 \\left(\\left(7 M_{1} + 3 M_{2}\\right) \\left(2 g \\sin{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - l_{1} \\omega_{1}^{2}{\\left(t \\right)} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}\\right) - 3 \\cdot \\left(2 M_{1} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + 2 M_{2} g \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + M_{2} l_{2} \\omega_{2}^{2}{\\left(t \\right)} \\sin{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}\\right) \\cos{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)}\\right)}{l_{2} \\left(- 49 M_{1} + 9 M_{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} - \\theta_{2}{\\left(t \\right)} \\right)} - 21 M_{2}\\right)}$"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(eqn4)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-12-09T22:38:21.093846Z",
     "start_time": "2023-12-09T22:38:21.076967Z"
    }
   },
   "id": "894a7c8656f91546"
  },
  {
   "cell_type": "markdown",
   "id": "7a4c42935ce25cc4",
   "metadata": {},
   "source": [
    "----\n",
    "&nbsp;\n",
    "### Numerical Integration in [DoublePendulum Class](DoublePendulum_Compound.py)"
   ]
  }
 ],
 "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.9.12"
  },
  "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": 5
}