{ "cells": [ { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "import sympy as smp\n", "from matplotlib import animation\n", "from matplotlib.animation import PillowWriter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A diagram of the problem at hand:\n", "\n", "* The entire system depends on only two angles: $\\theta_1$ and $\\theta_2$\n", "* We will use these variables as our free variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define all symbols we need for this problem using sympy:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "t, g, l1, l2, m1, m2, m3, k, L0 = smp.symbols('t g l_1 l_2 m_1 m_2 m_3 k L_0')\n", "the1, the2 = smp.symbols(r'\\theta_1 \\theta_2', cls=smp.Function)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define $\\theta(t)$ and $\\dot{\\theta}(t)$ and $\\ddot{\\theta}(t)$" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "the1 = the1(t)\n", "the2 = the2(t)\n", "the1_d = smp.diff(the1, t)\n", "the2_d = smp.diff(the2, t)\n", "the1_dd = smp.diff(the1_d, t)\n", "the2_dd = smp.diff(the2_d, t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the $x$ and $y$ coordinates of all three masses" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "x1 = l1*smp.cos(the1)\n", "y1 = -l1*smp.sin(the1)\n", "x2 = 2*x1\n", "y2 = 0\n", "x3 = x2 + l2*smp.sin(the2)\n", "y3 = -l2*smp.cos(the2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define both kinetic and potential energy:\n", "\n", "* Kinetic energy $T$ comes from the motion of the three masses\n", "* Potential energy $V$ comes from both the gravitational potential energy of the masses $mgy$ and the potential energy in the spring $\\frac{1}{2}kx^2$ where $x=x_2 - L_0$" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "T = smp.Rational(1,2) * m1 * (smp.diff(x1,t)**2 + smp.diff(y1,t)**2) \\\n", " +smp.Rational(1,2) * m2 * (smp.diff(x2,t)**2 + smp.diff(y2,t)**2) \\\n", " +smp.Rational(1,2) * m3 * (smp.diff(x3,t)**2 + smp.diff(y3,t)**2)\n", "\n", "V = m1*g*y1 + m2*g*y2 + m3*g*y3 +smp.Rational(1,2) * k * (x2-L0)**2\n", "L =T-V" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can look at the Lagrangian" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle g l_{1} m_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} + g l_{2} m_{3} \\cos{\\left(\\theta_{2}{\\left(t \\right)} \\right)} - \\frac{k \\left(- L_{0} + 2 l_{1} \\cos{\\left(\\theta_{1}{\\left(t \\right)} \\right)}\\right)^{2}}{2} + 2 l_{1}^{2} m_{2} \\sin^{2}{\\left(\\theta_{1}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2} + \\frac{m_{1} \\left(l_{1}^{2} \\sin^{2}{\\left(\\theta_{1}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2} + l_{1}^{2} \\cos^{2}{\\left(\\theta_{1}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{1}{\\left(t \\right)}\\right)^{2}\\right)}{2} + \\frac{m_{3} \\left(l_{2}^{2} \\sin^{2}{\\left(\\theta_{2}{\\left(t \\right)} \\right)} \\left(\\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2} + \\left(- 2 l_{1} \\sin{\\left(\\theta_{1}{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta_{1}{\\left(t \\right)} + l_{2} \\cos{\\left(\\theta_{2}{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta_{2}{\\left(t \\right)}\\right)^{2}\\right)}{2}$" ], "text/plain": [ "g*l_1*m_1*sin(\\theta_1(t)) + g*l_2*m_3*cos(\\theta_2(t)) - k*(-L_0 + 2*l_1*cos(\\theta_1(t)))**2/2 + 2*l_1**2*m_2*sin(\\theta_1(t))**2*Derivative(\\theta_1(t), t)**2 + m_1*(l_1**2*sin(\\theta_1(t))**2*Derivative(\\theta_1(t), t)**2 + l_1**2*cos(\\theta_1(t))**2*Derivative(\\theta_1(t), t)**2)/2 + m_3*(l_2**2*sin(\\theta_2(t))**2*Derivative(\\theta_2(t), t)**2 + (-2*l_1*sin(\\theta_1(t))*Derivative(\\theta_1(t), t) + l_2*cos(\\theta_2(t))*Derivative(\\theta_2(t), t))**2)/2" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get Lagrange's equations\n", "\n", "$$\\frac{\\partial L}{\\partial \\theta_1} - \\frac{d}{dt}\\frac{\\partial L}{\\partial \\dot{\\theta_1}} = 0$$\n", "$$\\frac{\\partial L}{\\partial \\theta_2} - \\frac{d}{dt}\\frac{\\partial L}{\\partial \\dot{\\theta_2}} = 0$$" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "LE1 = smp.diff(L, the1) - smp.diff(smp.diff(L, the1_d), t).simplify()\n", "LE2 = smp.diff(L, the2) - smp.diff(smp.diff(L, the2_d), t).simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve Lagranges equations (this assumes that `LE1` and `LE2` are both equal to zero)\n", "\n", "* We solve these equations (which are **linear** in $\\ddot{\\theta_1}$ and $\\ddot{\\theta_2}$) for $\\ddot{\\theta_1}$ and $\\ddot{\\theta_2}$." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "sols = smp.solve([LE1, LE2], (the1_dd, the2_dd),\n", " simplify=False, rational=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have \n", "\n", "* $\\frac{d^2 \\theta_1}{dt^2} = ...$\n", "* $\\frac{d^2 \\theta_2}{dt^2} = ...$\n", "\n", "These are two second order ODEs! In python we can only solve systems of first order ODEs. Any system of second order ODEs can be converted as follows:\n", "\n", "1. Define $z_1 = d\\theta_1/dt$ and $z_2=d\\theta_2/dt$\n", "2. Then $dz_1/dt = d^2\\theta_1/dt^2$ and $dz_2/dt = d^2\\theta_2/dt^2$\n", "\n", "Now we get a system of 4 first order ODEs (as opposed to 2 second order ones)\n", "\n", "* $d z_1/dt = ...$\n", "* $d\\theta_1/dt = z_1$\n", "* $d z_2/dt = ...$\n", "* $d\\theta_2/dt = z_1$\n", "\n", "We need to convert the **symbolic** expressions above to numerical functions so we can use them in a numerical python solver. For this we use `smp.lambdify`" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "dz1dt_f = smp.lambdify((t,g,k,L0,m1,m2,m3,l1,l2,the1,the2,the1_d,the2_d), sols[the1_dd])\n", "dz2dt_f = smp.lambdify((t,g,k,L0,m1,m2,m3,l1,l2,the1,the2,the1_d,the2_d), sols[the2_dd])\n", "dthe1dt_f = smp.lambdify(the1_d, the1_d)\n", "dthe2dt_f = smp.lambdify(the2_d, the2_d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now define $\\vec{S} = (\\theta_1, z_1, \\theta_2, z_2)$. IF we're going to use an ODE solver in python, we need to write a function that takes in $\\vec{S}$ and $t$ and returns $d\\vec{S}/dt$. In other words, we need to define $d\\vec{S}/dt (\\vec{S}, t)$\n", "\n", "* Our system of ODEs can be fully specified using $d\\vec{S}/dt$ and depends only on $\\vec{S}$ and $t$" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "def dSdt(S, t, g, k, L0, m1, m2, m3, l1, l2):\n", " the1, z1, the2, z2 = S\n", " return [\n", " dthe1dt_f(z1),\n", " dz1dt_f(t,g,k,L0,m1,m2,m3,l1,l2,the1,the2,z1,z2),\n", " dthe2dt_f(z2),\n", " dz2dt_f(t,g,k,L0,m1,m2,m3,l1,l2,the1,the2,z1,z2),\n", " ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the system of ODEs using scipys `odeint` method" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "t = np.linspace(0, 40, 1001) # s\n", "g = 9.81 #m/s^2\n", "k = 30 # N/m\n", "m1=2 # kg\n", "m2=2 # kg\n", "m3=1 # kg\n", "l1 = 1 # m\n", "l2 = 1 # m\n", "L0 = 1.5*l1 # m\n", "ans = odeint(dSdt, y0=[1, -1, -1, 1], t=t, args=(g, k, L0, m1, m2, m3, l1, l2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "25 times per second (number of data points). This will be important for animating later on." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Can obtain $\\theta_1(t)$ and $\\theta_2(t)$ from the answer" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "the1 = ans.T[0]\n", "the2 = ans.T[2]" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[