- Difficulty: 017/100
- Author: jeremy theler
- Keywords:
`PHASE_SPACE`

,`.=`

,`:=`

,`VAR`

,`INITIAL_CONDITIONS_MODE`

,`FROM_VARIABLES`

,`INCLUDE`

,`OUTPUT_FILE`

,`PRINT`

,`FILE`

,`CONST`

,`end_time`

,`derivative`

,`cos`

,

This is the continuation of a tale about mixing GSL and SUNDIALS that started here. It shows how wasora—which is defined as a *syntactically-sweetened way to ask a computer to perform a certain mathematical calculation*—can be used to solve Hamiltonian mechanical systems no only by just writing down the equations of motion in a text file but also without needing to do so! Just stating the Hamiltonian (or Lagrangian) of a system and them using the `derivative`

functional to have wasora to numerically compute the partial derivatives (using GSL and then to solve the resulting differential equations (using SUNDIALS). Moreover, a scientific video is generated from the computed time series with the wasora plugin besssugo (using SDL), much like those found in Wikipedia.

In this section, we solve the double pendulum that can be mentally built by attaching a massless string of length \(\ell_2\) with a bob of mass \(m_2\) to an existing simple pendulum (see section ref{012-mechanics}) of length \(\ell_1\) and mass \(m_1\). The two angles \(\theta_1\) and \(\theta_2\) are measured from the vertical and the positive direction is taken counter-clockwise:

Very much like the Lorenz system, the double pendulum displays a very complex behavior despite its simple concept. Paying attention only to the Lagrangian and Hamiltonian formulations (as the derivation of the Newton equations for this case are rather cumbersome), we can write the Lagrangian as the difference \(T - V\) of \(m_1\) and \(m_2\) kinetic and potential energies using the angles \(\theta_1\) and \(\theta_2\) as the generalized coordinates:

\[ \mathcal{L}(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) = \frac{1}{2} m_1 \ell_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 \left[\ell_1^2 \dot{\theta}_1^2 + \ell_2^2 \dot{\theta}_2^2 + 2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1-\theta_2)\right] + m_1 g \ell_1\cos \theta_1 + m_2 g \left(\ell_1\cos \theta_1 + \ell_2 \cos \theta_2 \right) \]

and the Hamiltonian as the sum \(T+V\), written as a function of the angles \(\theta_1\) and \(\theta_2\) and of the generalized angular momenta \(p_1 = \partial \mathcal{L}/\partial \dot{\theta}_1\) and \(p_2 = \partial \mathcal{L}/\partial \dot{\theta}_2\):

\[ \mathcal{H}(\theta_1, \theta_2, p_1, p_2) = - \frac{\ell_2^2 m_2 p_1^2 - 2 \ell_1 \ell_2 m_2 \cos(\theta_1-\theta_2) p_1 p_2 + \ell_1^2 (m_1+m_2) p_2^2} {\ell_1^2 \ell_2^2 m_2 \left[-2m_1-m_2+m_2\cos\Big(2(\theta_1-\theta_2)\Big)\right]} - \Big[ m_1 g \ell_1 \cos \theta_1 + m_2 g (\ell_1 \cos \theta_1 + \ell_2 \cos \theta_2) \Big] \]

After performing the tedious task of manually (and carefully) taking the analytical derivatives, we can arrive at the equations of motion. Lagrange says

\[ \begin{align*} 0 &= (m_1+m_2) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos(\theta_1-\theta_2) + m_2 \ell_1 \ell_2 \dot{\theta}_2^2 \sin(\theta_1-\theta_2) + \ell_1 (m_1+m_2) g \sin(\theta_1) \\ 0 &= m_2 \ell_2^2 \ddot{\theta}_2 + m_2 \ell_1 \ell_2 \ddot{\theta}_1 \cos(\theta_1-\theta_2) - m_2 \ell_1 \ell_2 \dot{\theta}_1^2 \sin(\theta_1-\theta_2) + \ell_2 m_2 g \sin(\theta_2) \end{align*} \]

whilst Hamilton states that

\[ \begin{align*} \dot{\theta}_1 &= \frac{p_1 \ell_2 - p_2 \ell_1 \cos(\theta_1-\theta_2)}{\ell_1^2 \ell_2 [m_1 + m_2 \sin^2(\theta_1-\theta_2)]} \\ \dot{\theta}_2 &= \frac{p_2 (m_1+m_2)/m_2 \ell_1 - p_1 \ell_2 \cos(\theta_1-\theta_2)}{\ell_1 \ell_2^2 [m_1 + m_2 \sin^2(\theta_1-\theta_2)]} \\ \dot{p_1} &= -(m_1+m_2) g \ell_1 \sin(\theta_1) - c_1 + c_2 \\ \dot{p_2} &= -m_2 g \ell_2 \sin(\theta_2) + c_1 - c_2 \end{align*} \]

where the expressions \(c_1\) and \(c_2\) are

\[ \begin{align*} c1 &= \frac{p_1 p_2 \sin(\theta_1-\theta_2)}{\ell_1 \ell_2 \Big[m_1+m_2 \sin(\theta_1-\theta_2)^2\Big]} \\ c2 &= \frac{\Big[ p_1^2 m_2 \ell_2^2 - 2 p_1 p_2 m_2 \ell_1 \ell_2 cos(\theta_1-\theta_2) + p_2^2 (m_1+m_2) \ell_1^2)\Big] \sin(2 (\theta_1-\theta_2)}{ 2 \ell_1^2 \ell_2^2 \left[m_1+m_2 \sin^2(\theta_1-\theta_2)\right]^2} \end{align*} \]

As here(real-012-mechanics.html), the problem formulation (i.e. the equations of motions) is selected from the commandline amongst four options:

- Lagrangian formulation with analytical derivatives
- Lagrangian formulation with numerical derivatives
- Hamiltonian formulation with analytical derivatives
- Hamiltonian formulation with numerical derivatives

The files that implement these formulations are:

Lagrange’s equations with the explicit analytical derivatives shown above reduced to a first-order system in `lagrange-analytical.was`

:

```
# the double pendulum solved by the lagrangian formulation
# and analytically computing its derivatives
PHASE_SPACE theta1 theta2 omega1 omega2
# reduction to a first-order system
omega1 .= theta1_dot
omega2 .= theta2_dot
# lagrange equations
0 .= (m1+m2)*l1^2*omega1_dot + \
m2*l1*l2*omega2_dot*cos(theta1-theta2) + \
m2*l1*l2*omega2^2*sin(theta1-theta2) + \
l1*(m1+m2)*g*sin(theta1)
0 .= m2*l2^2*omega2_dot + \
m2*l1*l2*omega1_dot*cos(theta1-theta2) - \
m2*l1*l2*omega1^2*sin(theta1-theta2) + \
l2*m2*g*sin(theta2)
```

Lagrange’s equations written as numerical derivatives of \(\mathcal{L}\) reduced to a first-order system in `lagrange-numerical.was`

:

```
# the double pendulum solved by the lagrangian formulation
# and numerically computing its derivatives
PHASE_SPACE theta1 theta2 dL_dthetadot1 dL_dthetadot2
VAR theta1' theta2' theta1_dot' theta2_dot'
L(theta1,theta2,theta1_dot,theta2_dot) := {
# kinetic energy of m1
1/2*m1*l1^2*theta1_dot^2 +
# kinetic energy of m2
1/2*m2*(l1^2*theta1_dot^2 + l2^2*theta2_dot^2 + 2*l1*l2*theta1_dot*theta2_dot*cos(theta1-theta2))
+ (
# potential energy of m1
m1*g * l1*cos(theta1) +
# potential energy of m2
m2*g * (l1*cos(theta1) + l2*cos(theta2))
) }
dL_dthetadot1 .= derivative(L(theta1, theta2, theta1_dot', theta2_dot), theta1_dot', theta1_dot)
dL_dthetadot2 .= derivative(L(theta1, theta2, theta1_dot, theta2_dot'), theta2_dot', theta2_dot)
dL_dthetadot1_dot .= derivative(L(theta1', theta2,theta1_dot, theta2_dot), theta1', theta1)
dL_dthetadot2_dot .= derivative(L(theta1, theta2',theta1_dot, theta2_dot), theta2', theta2)
```

Hamiltons’s equations with the explicit analytical derivatives shown above in `hamilton-analytical.was`

:

```
# the double pendulum solved by the hamiltonian formulation
# and analytically computing its derivatives
PHASE_SPACE theta1 theta2 p1 p2 c1 c2
theta1_dot .= (p1*l2 - p2*l1*cos(theta1-theta2))/(l1^2*l2*(m1 + m2*sin(theta1-theta2)^2))
theta2_dot .= (p2*(m1+m2)/m2*l1 - p1*l2*cos(theta1-theta2))/(l1*l2^2*(m1 + m2*sin(theta1-theta2)^2))
p1_dot .= -(m1+m2)*g*l1*sin(theta1) - c1 + c2
p2_dot .= -m2*g*l2*sin(theta2) + c1 - c2
c1 .= p1*p2*sin(theta1-theta2)/(l1*l2*(m1+m2*sin(theta1-theta2)^2))
c2 .= { (p1^2*m2*l2^2 - 2*p1*p2*m2*l1*l2*cos(theta1-theta2)
+ p2^2*(m1+m2)*l1^2)*sin(2*(theta1-theta2))/
(2*l1^2*l2^2*(m1+m2*sin(theta1-theta2)^2)^2) }
```

Hamilton’s equations written as numerical derivatives of \(\mathcal{H}\) in `hamilton-numerical.was`

:

```
# the double pendulum solved by the hamiltonian formulation
# and numerically computing its derivatives
PHASE_SPACE theta1 theta2 p1 p2
VAR theta1' theta2' p1' p2'
H(theta1,theta2,p1,p2) := {
- (m1*g*l1*cos(theta1) + m2*g*(l1*cos(theta1) + l2*cos(theta2)))
- (l2^2*m2*p1^2 - 2*l1*l2*m2*cos(theta1-theta2)*p1*p2 + l1^2*(m1+m2)*p2^2)
/(l1^2*l2^2*m2 * (-2*m1-m2+m2*cos(2*(theta1-theta2))))
}
theta1_dot .= +derivative(H(theta1,theta2,p1',p2), p1', p1)
theta2_dot .= +derivative(H(theta1,theta2,p1,p2'), p2', p2)
p1_dot .= -derivative(H(theta1',theta2,p1,p2), theta1', theta1)
p2_dot .= -derivative(H(theta1,theta2',p1,p2), theta2', theta2)
```

The main input is called `double.was`

, which should be called with an argument stating which of the four formulations should be used to solve the double pendulum:

```
# the double pendulum
# solved using different formulations
# parameters
CONST g m1 l1 m2 l2
end_time = 15
m1 = 0.3
m2 = 0.2
l1 = 0.3
l2 = 0.25
g = 9.8
# inital conditions
INITIAL_CONDITIONS_MODE FROM_VARIABLES
theta1_0 = pi/2
theta2_0 = pi
# include the selected formulation
INCLUDE $1.was
# output the results vs. time
PRINT FILE_PATH $1.dat t theta1 theta2 theta1_dot theta2_dot \
l1*sin(theta1) -l1*cos(theta1) \
l1*sin(theta1)+l2*sin(theta2) -l1*cos(theta1)-l2*cos(theta2)
```

```
$ time wasora double.was lagrange-analytical
real 0m0.050s
user 0m0.040s
sys 0m0.008s
$ time wasora double.was lagrange-numerical
real 0m0.357s
user 0m0.352s
sys 0m0.000s
$ time wasora double.was hamilton-analytical
real 0m0.120s
user 0m0.116s
sys 0m0.000s
$ time wasora double.was hamilton-numerical
real 0m0.316s
user 0m0.312s
sys 0m0.000s
$ pyxplot double.ppl
$
```

The script file `double.ppl`

produces the four figures above, which show \(\theta_1\), \(\theta_2\), \(\dot{\theta_1}\) and \(\dot{\theta_2}\) vs. time for the four formulations. The chaotic nature of the double pendulum is illustrated by noting that for small \(t\), all the solutions coincide. However, they start to slightly diverge due to small differences of the numerical integration of the differential equations, which are also slightly different amongst them. These differences start to accumulate and towards the end of the computations, every solution is completely different showing a high sensitivity to perturbations.

In the terminal mimic above a rough measure of the CPU time needed to solve each case is given using the standard UNIX tool `time`

. It can be seen that the Lagrangian formulation with analytical derivatives is the fastest, consuming less than half computational effort than the Hamiltonian formulation with analytical derivatives and eight times less than both completely numerical formulations. However, it should be noted that nowadays (i.e. mid 2010s) CPU time is far less expensive than human expert time, so an eight-fold increase in CPU time is a bargain compared to paying someone to perform analytical derivatives of either the Lagrangian or Hamiltonian, considering that there is no way of telling which of the four solutions above is more accurate.

The reader is encouraged to try different combinations of parameters (masses, lengths and gravity) and initial conditions. Even parametric studies can be performed, as shown here with the period vs. amplitude study.

Once we used wasora to solve the four proposed formulations of the double pendulum, we can use the plugin besssugo to generate a video which better illustrates the chaotic nature of the system. And, of course, we cannot deny that videos are far more sexy than static plots of variables versus time.

On the one hand, we know that videos usually have a fixed frames-per-second value whilst numerical integration of DAE systems have variable time steps. And, on the other hand, we want to have a single video comparing four different solutions that involve functions of time but that are sampled at different times. Therefore, it is best to load the data to be drawn into the video from the output of previous executions instead of solving again the problems at the moment of generating the video. We thus load the positions \(x_1(t)\), \(y_1(t)\), \(x_2(t)\) and \(y_2(t)\) of the centers of both bobs which are located in the sixth, seventh, eighth and ninth column of the output files generated in the previous subsection. As for sure we will have to interpolate the data because the time steps of the video in general will not coincide with the data, we choose the akima interpolation method to obtain smooth functions of time.

We name the cases as

- Lagrange with analytical derivatives
- Lagrange with numerical derivatives
- Hamilton with analytical derivatives
- Hamilton with numerical derivatives

The video is divided into a \(2 \times 2\) grid showing the Lagrangian (cases \(a\) & \(b\)) on the upper half and Hamiltonian (cases \(c\) & \(d\)) on the lower half, with the analytical derivatives (cases \(a\) and \(c\)) on the left half and the numerical derivatives (cases \(b\) and \(d\)) on the right half.

The wasora+besssugo input has `dt=1/200`

to have smooth trajectories. One out of five individual frames are dumped as PNG image fles, giving 40 frames per second. Then the frames are converted into a WebM video using avconv with 20 frames per second, resulting in a video that advances at half the real time speed so as to show with some detail what the pendulums but avoiding boring the reader with an extremely slow motion. Enjoy chaos!

```
# video of four double pendulums
FUNCTION x1a(t) FILE_PATH lagrange-analytical.dat COLUMNS 1 6 INTERPOLATION splines
FUNCTION y1a(t) FILE_PATH lagrange-analytical.dat COLUMNS 1 7 INTERPOLATION splines
FUNCTION x2a(t) FILE_PATH lagrange-analytical.dat COLUMNS 1 8 INTERPOLATION splines
FUNCTION y2a(t) FILE_PATH lagrange-analytical.dat COLUMNS 1 9 INTERPOLATION splines
FUNCTION x1b(t) FILE_PATH lagrange-numerical.dat COLUMNS 1 6 INTERPOLATION splines
FUNCTION y1b(t) FILE_PATH lagrange-numerical.dat COLUMNS 1 7 INTERPOLATION splines
FUNCTION x2b(t) FILE_PATH lagrange-numerical.dat COLUMNS 1 8 INTERPOLATION splines
FUNCTION y2b(t) FILE_PATH lagrange-numerical.dat COLUMNS 1 9 INTERPOLATION splines
FUNCTION x1c(t) FILE_PATH hamilton-analytical.dat COLUMNS 1 6 INTERPOLATION splines
FUNCTION y1c(t) FILE_PATH hamilton-analytical.dat COLUMNS 1 7 INTERPOLATION splines
FUNCTION x2c(t) FILE_PATH hamilton-analytical.dat COLUMNS 1 8 INTERPOLATION splines
FUNCTION y2c(t) FILE_PATH hamilton-analytical.dat COLUMNS 1 9 INTERPOLATION splines
FUNCTION x1d(t) FILE_PATH hamilton-numerical.dat COLUMNS 1 6 INTERPOLATION splines
FUNCTION y1d(t) FILE_PATH hamilton-numerical.dat COLUMNS 1 7 INTERPOLATION splines
FUNCTION x2d(t) FILE_PATH hamilton-numerical.dat COLUMNS 1 8 INTERPOLATION splines
FUNCTION y2d(t) FILE_PATH hamilton-numerical.dat COLUMNS 1 9 INTERPOLATION splines
end_time = x1a_b
skip_frames = 5
dt = 1/200
m1 = 0.3
m2 = 0.2
l1 = 0.3
l2 = 0.25
g = 9.8
# size of each canvas [ px ]
s = 220
# constants
k1 = 0.4
k2 = 1.48
k3 = 0.1
f = 0.27
WINDOW WIDTH 2*s HEIGHT 2*s CAPTION "four almost-exactly-equal double pendulums"
IMAGE FILE_PATH back.png NO_RELOAD
# trajectory canvases
CANVAS traj_a WIDTH 1 HEIGHT 1 TOP 1 LEFT -1 PERSISTENT COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS traj_b WIDTH 1 HEIGHT 1 TOP 1 LEFT 0 PERSISTENT COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS traj_c WIDTH 1 HEIGHT 1 TOP 0 LEFT -1 PERSISTENT COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS traj_d WIDTH 1 HEIGHT 1 TOP 0 LEFT 0 PERSISTENT COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
phase_r = 0.5 + 0.40*sin(4*sqrt(2)*t+1)
phase_g = 0.5 + 0.42*sin(4*sqrt(3)*t+2)
phase_b = 0.5 + 0.44*sin(pi*t+3)
SEGMENT CANVAS traj_a X1 x2a(t-dt) Y1 y2a(t-dt) X2 x2a(t) Y2 y2a(t) COLOR phase_r phase_g phase_b NO_AA
SEGMENT CANVAS traj_b X1 x2b(t-dt) Y1 y2b(t-dt) X2 x2b(t) Y2 y2b(t) COLOR phase_r phase_g phase_b NO_AA
SEGMENT CANVAS traj_c X1 x2c(t-dt) Y1 y2c(t-dt) X2 x2c(t) Y2 y2c(t) COLOR phase_r phase_g phase_b NO_AA
SEGMENT CANVAS traj_d X1 x2d(t-dt) Y1 y2d(t-dt) X2 x2d(t) Y2 y2d(t) COLOR phase_r phase_g phase_b NO_AA
# pendulum canvases
CANVAS pend_a WIDTH 1 HEIGHT 1 TOP 1 LEFT -1 COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS pend_b WIDTH 1 HEIGHT 1 TOP 1 LEFT 0 COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS pend_c WIDTH 1 HEIGHT 1 TOP 0 LEFT -1 COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
CANVAS pend_d WIDTH 1 HEIGHT 1 TOP 0 LEFT 0 COORDINATES 0.5*s k1*s k2*(l1+l2)*s -k2*(l1+l2)*s k2*(l1+l2)*s
# strings 1
SEGMENT CANVAS pend_a X1 0 Y1 0 X2 x1a(t) Y2 y1a(t)
SEGMENT CANVAS pend_b X1 0 Y1 0 X2 x1b(t) Y2 y1b(t)
SEGMENT CANVAS pend_c X1 0 Y1 0 X2 x1c(t) Y2 y1c(t)
SEGMENT CANVAS pend_d X1 0 Y1 0 X2 x1d(t) Y2 y1d(t)
# strings 2
SEGMENT CANVAS pend_a X1 x1a(t) Y1 y1a(t) X2 x2a(t) Y2 y2a(t)
SEGMENT CANVAS pend_b X1 x1b(t) Y1 y1b(t) X2 x2b(t) Y2 y2b(t)
SEGMENT CANVAS pend_c X1 x1c(t) Y1 y1c(t) X2 x2c(t) Y2 y2c(t)
SEGMENT CANVAS pend_d X1 x1d(t) Y1 y1d(t) X2 x2d(t) Y2 y2d(t)
# bobs 1 (and the light reflection)
CIRCLE CANVAS pend_a X x1a(t) Y y1a(t) RAD m1*k3 COLOR black
CIRCLE CANVAS pend_a X x1a(t)+f*m1*k3 Y y1a(t)+f*m1*k3 RAD m1*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_b X x1b(t) Y y1b(t) RAD m1*k3 COLOR black
CIRCLE CANVAS pend_b X x1b(t)+f*m1*k3 Y y1b(t)+f*m1*k3 RAD m1*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_c X x1c(t) Y y1c(t) RAD m1*k3 COLOR black
CIRCLE CANVAS pend_c X x1c(t)+f*m1*k3 Y y1c(t)+f*m1*k3 RAD m1*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_d X x1d(t) Y y1d(t) RAD m1*k3 COLOR black
CIRCLE CANVAS pend_d X x1d(t)+f*m1*k3 Y y1d(t)+f*m1*k3 RAD m1*0.7*k3*f COLOR snow NO_AA
# bobs 2 (and the light reflection)
CIRCLE CANVAS pend_a X x2a(t) Y y2a(t) RAD m2*k3 COLOR black
CIRCLE CANVAS pend_a X x2a(t)+f*m2*k3 Y y2a(t)+f*m1*k3 RAD m2*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_b X x2b(t) Y y2b(t) RAD m2*k3 COLOR black
CIRCLE CANVAS pend_b X x2b(t)+f*m2*k3 Y y2b(t)+f*m1*k3 RAD m2*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_c X x2c(t) Y y2c(t) RAD m2*k3 COLOR black
CIRCLE CANVAS pend_c X x2c(t)+f*m2*k3 Y y2c(t)+f*m1*k3 RAD m2*0.7*k3*f COLOR snow NO_AA
CIRCLE CANVAS pend_d X x2d(t) Y y2d(t) RAD m2*k3 COLOR black
CIRCLE CANVAS pend_d X x2d(t)+f*m2*k3 Y y2d(t)+f*m1*k3 RAD m2*0.7*k3*f COLOR snow NO_AA
BESSSUGO_DRAW
IF mod(step_transient,skip_frames)=0
FILE frame frame-%04g.png round(step_transient/skip_frames) MODE w
BESSSUGO_DUMP FILE frame
ENDIF
```

```
$ rm -f frame-*.png
$ besssugo trajectories.was
$ avconv -v quiet -y -f image2 -framerate 20 -i frame-%04d.png -r 20 -b 800k trajectories.webm
$ rm -f frame-*.png
$
```