{ "cells": [ { "cell_type": "markdown", "id": "b7c0470a-41c2-4728-ba1e-88aceac56cc2", "metadata": {}, "source": [ "## Double pendulum" ] }, { "cell_type": "markdown", "id": "d8d9a01b-e486-46f0-be18-fb8390dc9116", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "a8f2c1cb-e326-4446-8ee9-72d0b6055d62", "metadata": {}, "source": [ "Lagrangian equations of motion for the double pendulum:\n", "\\begin{equation}\n", "\\left( \n", " \\begin{matrix}\n", " \\ddot{\\theta}_1 \\\\ \n", " \\ddot{\\theta}_2 \n", " \\end{matrix} \n", " \\right) =\n", " \\frac{1}{1 - \\alpha_1\\alpha_2}\n", " \\left( \n", " \\begin{matrix} \n", " f_1 - \\alpha_1 f_2\\\\ \n", " -\\alpha_2 f_1 + f_2 \n", " \\end{matrix} \n", " \\right),\n", " \\end{equation}" ] }, { "cell_type": "markdown", "id": "2e4d95cb-fa66-4702-bb28-96fbaca0d0d9", "metadata": {}, "source": [ "where:\n", "\\begin{align}\n", " \\alpha_1(\\theta_1, \\theta_2) \n", " &= \n", " \\frac{l_2}{l_1} \\frac{m_2}{m_1 + m_2} \\cos(\\theta_1 - \\theta_2),\n", " \\\\\n", " \\alpha_2(\\theta_1, \\theta_2) \n", " &= \n", " \\frac{l_1}{l_2} \\cos(\\theta_1-\\theta_2),\n", " \\\\\n", " f_1(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2) \n", " &=\n", " -\\frac{l_2}{l_1} \\frac{m_2}{m_1+m_2} \\, \\dot{\\theta}_2^2 \\sin(\\theta_1 - \\theta_2)\n", " - \\frac{g}{l_1} \\sin\\theta_1,\n", " \\\\\n", " f_2(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2) &\n", " =\n", " \\frac{l_1}{l_2} \\, \\dot{\\theta}_1^2 \\sin(\\theta_1-\\theta_2) \n", " - \\frac{g}{l_2} \\sin\\theta_2.\n", "\\end{align}" ] }, { "cell_type": "markdown", "id": "2b012bbf-8ee7-4613-962f-bc032d566ae9", "metadata": {}, "source": [ "Let $$\\omega_1^2 = \\frac{g}{l_1}, \\quad \\nu = \\frac{l_1}{l_2}, \\quad \\omega_2^2 = \\frac{g}{l_2} = \\nu \\omega_1^2 ,$$\n", "$$\\kappa = \\frac{1}{\\nu} \\frac{m_2}{m_1 + m_2}.$$\n", "Let's introduce dimensionless time $$\\tau = \\omega_1 t.$$" ] }, { "cell_type": "markdown", "id": "5b1cdcb6-be12-40f4-a0d4-eb8af4e18858", "metadata": {}, "source": [ "Then, $$\n", "\\frac{\\mathrm{d}}{\\mathrm{d}t} = \\omega_1 \\frac{\\mathrm{d}}{\\mathrm{d}\\tau}, \\quad\n", "\\frac{\\mathrm{d}^2}{\\mathrm{d}t^2} = \\omega_1^2 \\frac{\\mathrm{d}^2}{\\mathrm{d}\\tau^2}.$$" ] }, { "cell_type": "markdown", "id": "5591511f-5feb-42fe-88c0-98af885a4fed", "metadata": {}, "source": [ "We can write the equation of motion of the pendulum in the following dimensionless form:\n", "\\begin{equation}\n", "\\frac{\\mathrm{d}^2}{\\mathrm{d}\\tau^2} \\left( \n", " \\begin{matrix}\n", " \\theta_1 \\\\ \n", " \\theta_2 \n", " \\end{matrix} \n", " \\right) =\n", " \\frac{1}{1 - \\bar{\\alpha}_1\\bar{\\alpha}_2}\n", " \\left( \n", " \\begin{matrix} \n", " \\bar{f}_1 - \\bar{\\alpha}_1 \\bar{f}_2\\\\ \n", " -\\bar{\\alpha}_2 \\bar{f}_1 + \\bar{f}_2 \n", " \\end{matrix} \n", " \\right),\n", " \\end{equation}" ] }, { "cell_type": "markdown", "id": "66e3ed88-2bf8-4e70-b3c1-d4091f7ae7da", "metadata": {}, "source": [ "where\n", "\\begin{align}\n", " \\bar{\\alpha}_1(\\theta_1, \\theta_2) \n", " &= \n", " \\kappa \\cos(\\theta_1 - \\theta_2),\n", " \\\\\n", " \\bar{\\alpha}_2(\\theta_1, \\theta_2) \n", " &= \n", " \\nu \\cos(\\theta_1-\\theta_2),\n", " \\\\\n", " \\bar{f}_1(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2) \n", " &=\n", " -\\kappa \\, \\dot{\\theta}_2^2 \\sin(\\theta_1 - \\theta_2)\n", " - \\sin\\theta_1,\n", " \\\\\n", " \\bar{f}_2(\\theta_1, \\theta_2, \\dot{\\theta}_1, \\dot{\\theta}_2) &\n", " =\n", " \\nu \\, \\dot{\\theta}_1^2 \\sin(\\theta_1-\\theta_2) \n", " - \\nu \\sin\\theta_2.\n", "\\end{align}" ] }, { "cell_type": "markdown", "id": "83a28d5a-7f9a-4dc8-a45c-6b925b7d6a59", "metadata": {}, "source": [ "Here\n", "$$\\dot{\\theta}_1 = \\frac{\\mathrm{d}\\theta_1}{\\mathrm{d}\\tau}, \\quad \\dot{\\theta}_2 = \\frac{\\mathrm{d}\\theta_2}{\\mathrm{d}\\tau}$$" ] }, { "cell_type": "code", "execution_count": null, "id": "8d0baa5a", "metadata": {}, "outputs": [], "source": [ "using PyPlot" ] }, { "cell_type": "code", "execution_count": null, "id": "f0768159", "metadata": {}, "outputs": [], "source": [ "using OrdinaryDiffEqTsit5" ] }, { "cell_type": "code", "execution_count": null, "id": "78c51fe3", "metadata": {}, "outputs": [], "source": [ "using Printf" ] }, { "cell_type": "code", "execution_count": null, "id": "36da7e57", "metadata": {}, "outputs": [], "source": [ "using Colors" ] }, { "cell_type": "code", "execution_count": null, "id": "93a4be10", "metadata": {}, "outputs": [], "source": [ "function doublependulum!(dudt, u, p, t)\n", " # unpack\n", " theta1 = u[1]\n", " theta2 = u[2]\n", " dtheta1dt = u[3]\n", " dtheta2dt = u[4]\n", " κ, ν = p\n", " \n", " sn, cs = sincos(theta1 - theta2)\n", " a1 = κ*cs\n", " a2 = ν*cs\n", " γ = 1/(1.0 - a1*a2)\n", " \n", " f1 = -κ*sn*dtheta2dt^2 - sin(theta1)\n", " f2 = ν*(sn*dtheta1dt^2 - ν*sin(theta2))\n", "\n", " dudt[1] = u[3]\n", " dudt[2] = u[4]\n", " dudt[3] = γ*(f1 - a1*f2)\n", " dudt[4] = γ*(f2 - a2*f1)\n", " return nothing\n", "end\n" ] }, { "cell_type": "code", "execution_count": null, "id": "13962d21-ccd5-4130-b6ec-0f10e6624eb0", "metadata": {}, "outputs": [], "source": [ "# Parameters of the double pendulum\n", "l1 = 1.0 # length of pendulum1\n", "l2 = 1.0 # length of pendulum2\n", "m1 = 1.0 # mass of pendulum1\n", "m2 = 1.0 # mass of pendulum2" ] }, { "cell_type": "code", "execution_count": null, "id": "5025540b-ab22-47e6-97a7-9d1555841fb3", "metadata": {}, "outputs": [], "source": [ "# Parameters in the equations\n", "ν = l1/l2\n", "κ = 1 / ν * m2 / (m1 + m2)\n", "# Pack the parameters\n", "p = (κ, ν);" ] }, { "cell_type": "code", "execution_count": null, "id": "f48164f3-e4e0-412d-be39-97314faf1cf3", "metadata": {}, "outputs": [], "source": [ "tspan = (0.0, 40.0)\n", "initial = [pi/1.6, pi/1.8, 0.0, 0.0];" ] }, { "cell_type": "code", "execution_count": null, "id": "b50cda3b", "metadata": {}, "outputs": [], "source": [ "\n", "prob = ODEProblem(doublependulum!, initial, tspan, p)\n", "sol = solve(prob, Tsit5(), abstol=1e-7, reltol=1e-7);" ] }, { "cell_type": "code", "execution_count": null, "id": "d8d46f7c", "metadata": {}, "outputs": [], "source": [ "\n", "np = 2001\n", "t = range(tspan[1], tspan[2], np)\n", "sl = sol(t);" ] }, { "cell_type": "code", "execution_count": null, "id": "b5103796", "metadata": {}, "outputs": [], "source": [ "\n", "theta1 = sl[1, :]\n", "theta2 = sl[2, :];" ] }, { "cell_type": "code", "execution_count": null, "id": "549761d8", "metadata": {}, "outputs": [], "source": [ "\n", "# Cartesian coordinates of the masses\n", "x1 = l1 .* sin.(theta1)\n", "y1 = -l1 .* cos.(theta1)\n", "x2 = x1 .+ l2 .* sin.(theta2)\n", "y2 = y1 .- l2 .* cos.(theta2);" ] }, { "cell_type": "code", "execution_count": null, "id": "5d692dfe", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " snapshot_ft(nt, pts_disp, xl, xr, yb, yt, t, x1, y1, x2, y2)\n", "\n", "Display the current positions (x1[nt], y1[nt]) and (x2[nt], y2(nt]), and \n", "the trajectories' fading tails - the most recent pts_disp points\n", "\"\"\"\n", "function snapshot_ft(nt, pts_disp, box, t, x1, y1, x2, y2)\n", " nt_start = max(1, nt-pts_disp)\n", " xl, xr, yb, yt = box\n", " axis(\"square\")\n", " xlim(xl, xr)\n", " ylim(yb, yt)\n", " titl = @sprintf(\"Double pendulum t = %5.2f .. %5.2f\",\n", " t[nt_start], t[nt])\n", " title(titl, family=\"monospace\")\n", " xlabel(L\"$x$\")\n", " ylabel(L\"$y$\")\n", " grid(true)\n", "\n", " # pivot\n", " plot(0.0, 0.0, color=\"black\", marker=\"x\", markersize=5.0)\n", "\n", " # two rods\n", " plot([0.0, x1[nt],x2[nt]], [0.0, y1[nt], y2[nt]], color=\"black\")\n", "\n", " # mass 1\n", " c1 = \"green\"\n", " plot(x1[nt], y1[nt], color=c1, marker=\"o\", markersize=5.0)\n", " \n", " # mass 2\n", " c2 = \"blue\"\n", " plot(x2[nt], y2[nt], color=c2, marker=\"o\", markersize=5.0)\n", "\n", " # fading tail trajectory of mass 2\n", " c2p = parse(RGBA, c2)\n", " for i = nt_start:nt-1\n", " transp = 1.0 - ((nt-i-1)/(pts_disp-1))^2\n", " cc = (c2p.r, c2p.g, c2p.b, transp)\n", " plot([x2[i], x2[i+1]], [y2[i], y2[i+1]], color=cc, linewidth=0.5)\n", " end\n", " return nothing\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "c7f3b558-ce87-47cc-844f-92abfc003f1a", "metadata": {}, "outputs": [], "source": [ "\n", "# Define the plotting limits\n", "axis_lim = (l1 + l2)*1.2\n", "xl = -axis_lim\n", "xr = -xl\n", "yb = -axis_lim\n", "yt = -yb\n", "box = (xl, xr, yb, yt)" ] }, { "cell_type": "code", "execution_count": null, "id": "f3e00df0-87cc-4523-b879-58e487318c47", "metadata": {}, "outputs": [], "source": [ "\n", "# Define the length of the trajectory's tail\n", "f = 0.1 # Fraction of displayed points\n", "pts_disp = round(Int, f*np) # Number of displayed points" ] }, { "cell_type": "code", "execution_count": null, "id": "cc3eb74c-6b27-43f3-9f58-930c55c66bb4", "metadata": {}, "outputs": [], "source": [ "\n", "fig = figure()\n", "for nt = 1:10:np\n", " snapshot_ft(nt, pts_disp, box, t, x1, y1, x2, y2)\n", " display(fig)\n", " IJulia.clear_output(true)\n", " clf()\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "0e1073b7", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.10.5", "language": "julia", "name": "julia-1.10" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.10.5" } }, "nbformat": 4, "nbformat_minor": 5 }