Exercise 1.44: The double pendulum

This namespace explores Exercise 1.44 from Sussman and Wisdom's Structure and Interpretation of Classical Mechanics, using the SICMUtils Clojure library and the Clerk rendering environment.

Lagrangian

Start with a coordinate transformation from theta1, theta2 to rectangular coordinates. We'll generate our Lagrangian by composing this with an rectangular Lagrangian with the familiar form of T - V.

T describes the sum of the kinetic energy of two particles in rectangular coordinates.

V describes a uniform gravitational potential with coefficient g, acting on two particles with masses of, respectively, m1 and m2. Again, this is written in rectangular coordinates.

Form the rectangular Lagrangian L by subtracting (V m1 m2 g) from (T m1 m2):

Form the final Langrangian in generalized coordinates (the angles of each segment) by composing L-rect with a properly transformed angles->rect coordinate transform!

The Lagrangian is big and hairy:

Let's simplify that:

Better yet, let's render it as LaTeX, and create a helper function, render-eq to make it easier to render simplified equations:

And here are the equations of motion for the system:

What do these mean?

• the system has two degrees of freedom: $\theta_1$ and $\theta_2$.
• at any point t in time, the two equations above, full of first and second
• order derivatives of the position functions, will stay true
• the system can use these equations to simulate the system, one tick at a time.

Simulation

Next, let's run a simulation using those equations of motion and collect data on each coordinate's evolution.

Here are the constants specified in exercise 1.44:

masses in kg:

lengths in meters:

g in units of m/s^2:

And two sets of initial pairs of theta1, theta2 angles corresponding to chaotic and regular initial conditions:

Composing Lagrangian->state-derivative with L-double-pendulum produces a state derivative that we can use with our ODE solver:

Finally, two default parameters for our simulation. We'll record data in steps of 0.01 seconds, and simulate to a horizon of 50 seconds.

run! will return a sequence of 5001 states, one for each measured point in the simulation. The smaller-arity version simply passes in default masses and lengths, but you can override those with the larger arity version if you like.

(The interface here could use some work: integrate-state-derivative tidies this up a bit, but I want it out in the open for now.)

Run the simulation for each set of initial conditions and show the final state. Chaotic first:

Looks good:

Next, the regular initial condition:

Peek at the final state:

Measurements, Data Transformation

Next we'll chart the measurements trapped in those sequences of state tuples.

The exercise asks us to graph the energy of the system as a function of time. Composing Lagrangian->energy with L-double-pendulum gives us a new function (of a state tuple!) that will return the current energy in the system.:

energy-monitor returns a function of state that stores an initial energy value in its closure, and returns the delta of each new state's energy as compared to the original.

Finally, the large transform-data function. The charting library we'll use likes Clojure dictionaries; transform-data turns our raw data into a sequence of dictionary with all values we might care to explore. This includes:

• the generalized coordinate angles
• the generalized velocities of those angles
• the rectangular coordinates of each pendulum bob
• :d-energy, the error in energy at each point in the simulation

Here is transform-data:

The following forms transform the raw data for each initial condition and bind the results to chaotic-data and regular-data for exploration.

Here's the final, transformed chaotic state:

And the similar regular state:

Data Visualization Utilities

Vega-Lite allows us to visualizing the system.

I am not a pro here, but this does the trick for now.

First, a function to transform the dictionaries we generated above into a sequence of x, y coordinates tagged with distinct IDs for each pendulum bob's points:

segments-data generates endpoints of the pendulum segments, bob-to-bob:

Visualizations

a helper function that should be in clojure.core

With those tools in hand, let's make some charts. I'll call this first chart the system-inspector; this function will return a chart that will let us evolve the system with a time slider and monitor both angles, the energy error, and the pendulum bob itself as it evolves through time.

Here's a system monitor for the chaotic initial condition:

And again for the regular initial condition:

Generalized coordinates, velocities

angles-plot generates a plot, with no animation, showing both theta angles and their associated velocities:

Here are the angles for the chaotic initial condition:

And the regular initial condition:

Double Double Pendulum!

No visualizations here yet, but the code works well.

Looks good:

Next, the regular initial condition:

Peek at the final state: