Mitch Richling: Double Pendulum Chaos
| Author: | Mitch Richling |
| Updated: | 2025-01-23 |
Table of Contents
1. Introduction
This page was inspired by a reddit!
The equations of motion for the double pendulum may be written as a system of first degree differential equations:
\[\begin{eqnarray} \frac{\mathrm{d}\theta_1}{\mathrm{d}t} & = & \omega_1 \\ \frac{\mathrm{d}\theta_2}{\mathrm{d}t} & = & \omega_2 \\ \frac{\mathrm{d}\omega_1}{\mathrm{d}t} & = & \frac{-g (2 m_1+m_2) \sin(\theta_1)-m_2 g \sin(\theta_1-2 \theta_2)-2 m_2 (\omega_2^2 L_2+\omega_1^2 L_1 \cos(\theta_1-\theta_2))\sin(\theta_1-\theta_2)}{L_1 (2 m_1+m_2-m_2 \cos(2 \theta_1-2 \theta_2))} \\ \frac{\mathrm{d}\omega_2}{\mathrm{d}t} & = & \frac{2 \sin(\theta_1-\theta_2) (\omega_1^2 L_1 (m_1+m_2)+g (m_1+m_2) \cos(\theta_1)+\omega_2^2 L_2 m_2 \cos(\theta_1-\theta_2))}{L_2 (2 m_1+m_2-m_2 \cos(2 \theta_1-2 \theta_2))} \end{eqnarray}\]
Where
\[\begin{array}{lll} L_1 & = & \text{ Length of top bar} \\ L_2 & = & \text{ Length of bottom bar} \\ m_1 & = & \text{ Mass of top weight} \\ m_2 & = & \text{ Mass of bottom weight} \\ g & = & \text{ Acceleration of gravity} \\ \theta_1 & = & \text{ Vertical angle of top bar} \\ \theta_2 & = & \text{ Vertical angle of bottom bar} \\ \omega_1 & = & \text{ Angular velocity of top bar} \\ \omega_2 & = & \text{ Angular velocity of bottom bar} \end{array}\]
The idea is to run a double pendulum simulation for each pixel. Each pixel is mapped to a pair of angles – the x coordinate is mapped to \( \theta_1 \), and the y coordinate is mapped to \( \theta_2 \). These angles are used for the initial conditions for the double pendulum equations. We then solve these equations over a time span of 2 seconds using 200 steps of Euler's method. We dump an image of the final state of the simulation at the end.
Two ways are provided to map pixels to angles controlled by the boolean CENTER. If it's true, then the angles are mapped left to right & top to bottom
from \(0\) to \(2\pi\). If it's false, we map from \(-\pi\) to \(\pi\). The effect is that when CENTER is true the larger angles are at the
center of the image, and they are at the corners otherwise.
2. Code
#include "ramCanvas.hpp" #include "MRMathCPP.hpp" int main(void) { std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now(); const bool CENTER = true; const int IMXSIZ = 7680/2; const int IMYSIZ = 7680/2; const double m1 = 1.0; const double m2 = 1.0; const double L1 = 1.0; const double L2 = 1.0; const int steps = 200; const double h = 0.01; const double g = 9.80665; const double p = std::numbers::pi; mjr::ramCanvas3c8b theRamCanvas(IMXSIZ, IMYSIZ); if (CENTER) theRamCanvas.newRealCoords(0, 2*p, 0, 2*p); else theRamCanvas.newRealCoords(-p, p, -p, p); theRamCanvas.clrCanvasToBlack(); # pragma omp parallel for schedule(static,1) for(int y=0;y<theRamCanvas.getNumPixY();y++) { if ((y%100)==0) std::cout << y << std::endl; for(int x=0;x<theRamCanvas.getNumPixX();x++) { double theta1 = theRamCanvas.int2realX(x); double omega1 = 0.0; double theta2 = theRamCanvas.int2realY(y); double omega2 = 0.0; for(int i=0; i<steps; i++) { double cd = cos(theta1-theta2); double sd = sin(theta1-theta2); double o12 = omega1*omega1; double o22 = omega2*omega2; double denom = 2*m1+m2-m2*cos(2*theta1-2*theta2); double new_omega1 = (-g*(2*m1+m2)*sin(theta1)-m2*g*sin(theta1-2*theta2)-2*sd*m2*(o22*L2+o12*L1*cd))/(L1*denom); double new_omega2 = (2*sd*(o12*L1*(m1+m2)+g*(m1+m2)*cos(theta1)+o22*L2*m2*cd))/(L2*denom); theta1 += h*omega1; omega1 += h*new_omega1; theta2 += h*omega2; omega2 += h*new_omega2; } double r = mjr::math::ivl::wrapCO(theta2, 2*p)/(2*p); double b = mjr::math::ivl::wrapCO(theta1, 2*p)/(2*p); double g = mjr::math::ivl::wrapCC((omega1+omega2)/(4*p), 1.0); theRamCanvas.getPxColorRefNC(x, y).setChansRGB_dbl(r, g, b); } } theRamCanvas.writeTIFFfile((CENTER ? "doublePendulum_center.tiff" : "doublePendulum_corner.tiff")); std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime; std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl; return 0; }
Source Link: doublePendulum.cpp
The above program will generate the following images (the first is with CENTER=true):
3. Making A Movie!
We can animate the above process, and drop a movie frame for every time step. The program below demonstrates one of the neat features of MRaster – we can create images that have floating point doubles as pixels. I'm using that feature below to store the four values (\(\theta_1, \theta_2, \omega_1, \text{ and } \omega_2\)) necessary to save the state of the equations at each pixel. In a way, I'm doing the Euler method on equations being evaluated on images. This way we only have to do one step of the Euler method for each frame of the movie. How cool is that?
#include "ramCanvas.hpp" #include "MRMathCPP.hpp" int main(void) { std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now(); const bool CENTER = true; const int NUMFRM = 1080; const int IMXSIZ = 7680/8; const int IMYSIZ = 7680/8; const double m1 = 1.0; const double m2 = 1.0; const double L1 = 1.0; const double L2 = 1.0; const double h = 0.01; const double g = 9.80665; const double p = std::numbers::pi; mjr::ramCanvas1c64F theta1canvas(IMXSIZ, IMYSIZ); mjr::ramCanvas1c64F omega1canvas(IMXSIZ, IMYSIZ); mjr::ramCanvas1c64F theta2canvas(IMXSIZ, IMYSIZ); mjr::ramCanvas1c64F omega2canvas(IMXSIZ, IMYSIZ); mjr::ramCanvas3c8b pcolorCanvas(IMXSIZ, IMYSIZ); if (CENTER) pcolorCanvas.newRealCoords(0, 2*p, 0, 2*p); else pcolorCanvas.newRealCoords(-p, p, -p, p); pcolorCanvas.clrCanvasToBlack(); for(int frame=0; frame<NUMFRM; frame++) { # pragma omp parallel for schedule(static,1) for(int y=0;y<pcolorCanvas.getNumPixY();y++) { for(int x=0;x<pcolorCanvas.getNumPixX();x++) { double theta1, omega1, theta2, omega2; if (frame == 0) { theta1 = pcolorCanvas.int2realX(x); omega1 = 0.0; theta2 = pcolorCanvas.int2realY(y); omega2 = 0.0; } else { theta1 = theta1canvas.getPxColorRefNC(x, y).getC0(); omega1 = omega1canvas.getPxColorRefNC(x, y).getC0(); theta2 = theta2canvas.getPxColorRefNC(x, y).getC0(); omega2 = omega2canvas.getPxColorRefNC(x, y).getC0(); double cd = cos(theta1-theta2); double sd = sin(theta1-theta2); double o12 = omega1*omega1; double o22 = omega2*omega2; double denom = 2*m1+m2-m2*cos(2*theta1-2*theta2); double new_omega1 = (-g*(2*m1+m2)*sin(theta1)-m2*g*sin(theta1-2*theta2)-2*sd*m2*(o22*L2+o12*L1*cd))/(L1*denom); double new_omega2 = (2*sd*(o12*L1*(m1+m2)+g*(m1+m2)*cos(theta1)+o22*L2*m2*cd))/(L2*denom); theta1 += h*omega1; omega1 += h*new_omega1; theta2 += h*omega2; omega2 += h*new_omega2; } theta1canvas.drawPoint(x, y, theta1); omega1canvas.drawPoint(x, y, omega1); theta2canvas.drawPoint(x, y, theta2); omega2canvas.drawPoint(x, y, omega2); double r = mjr::math::ivl::wrapCO(theta2, 2*p)/(2*p); double b = mjr::math::ivl::wrapCO(theta1, 2*p)/(2*p); double g = mjr::math::ivl::wrapCC((omega1+omega2)/(4*p), 1.0); pcolorCanvas.getPxColorRefNC(x, y).setChansRGB_dbl(r, g, b); } } pcolorCanvas.writeTIFFfile((CENTER ? "doublePendulumM_center_" : "doublePendulumM_corner_") + mjr::math::str::fmt_int(frame, 4, '0') + ".tiff"); # pragma omp critical std::cout << "FRAME(" << frame << "): " << "DONE" << std::endl; } std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime; std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl; return 0; }
Source Link: doublePendulumM.cpp
The above program will generate the following movies (the first is with CENTER=true):