Mitch Richling: Double Pendulum Chaos

Author:	Mitch Richling
Updated:	2025-01-23

1. Introduction
2. Code
3. Making A Movie!
4. References

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{array}{rcl} (1) & \frac{d θ_{1}}{d t} & = & ω_{1} \\ (2) & \frac{d θ_{2}}{d t} & = & ω_{2} \\ (3) & \frac{d ω_{1}}{d t} & = & \frac{- g (2 m_{1} + m_{2}) \sin (θ_{1}) - m_{2} g \sin (θ_{1} - 2 θ_{2}) - 2 m_{2} (ω_{2}^{2} L_{2} + ω_{1}^{2} L_{1} \cos (θ_{1} - θ_{2})) \sin (θ_{1} - θ_{2})}{L_{1} (2 m_{1} + m_{2} - m_{2} \cos (2 θ_{1} - 2 θ_{2}))} \\ (4) & \frac{d ω_{2}}{d t} & = & \frac{2 \sin (θ_{1} - θ_{2}) (ω_{1}^{2} L_{1} (m_{1} + m_{2}) + g (m_{1} + m_{2}) \cos (θ_{1}) + ω_{2}^{2} L_{2} m_{2} \cos (θ_{1} - θ_{2}))}{L_{2} (2 m_{1} + m_{2} - m_{2} \cos (2 θ_{1} - 2 θ_{2}))} \end{array}$

Where

$\begin{array}{lll} L_{1} & = & Length of top bar \\ L_{2} & = & Length of bottom bar \\ m_{1} & = & Mass of top weight \\ m_{2} & = & Mass of bottom weight \\ g & = & Acceleration of gravity \\ θ_{1} & = & Vertical angle of top bar \\ θ_{2} & = & Vertical angle of bottom bar \\ ω_{1} & = & Angular velocity of top bar \\ ω_{2} & = & 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 $θ_{1}$ , and the y coordinate is mapped to $θ_{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 π$ . If it's false, we map from $- π$ to $π$ . 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 ( $θ_{1}, θ_{2}, ω_{1}, and ω_{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):

4. References

Check out the fractals section of my reading list.

All the code used to generate everything on this page may be found on github.

Mitch Richling: Double Pendulum Chaos

Table of Contents

1. Introduction

2. Code

3. Making A Movie!

4. References