Chaos in The Double Pendulum
Chaos in The Double Pendulum
Chaos in The Double Pendulum
y x
l1
1 m1
l
2 2 m2
The double pendulum is the simplest example of a time-independent Hamiltonian system that exhibits
chaotic behavior. Its phase space is 4-dimensional.
Chaos can be demonstrated experimentally in the double pendulum[1],[2].
1
Euler-Lagrange Equations
Suppose that the upper pendulum has a massless rod of length `1 and a bob of mass m1 . The lower
pedulum has a massless rod of length `2 and bob of mass m2 . The two rods provide constraints on the
motion of the two masses in the vertical plane:
The 1 equation is
h i
`1 (m1 + m2 )`1 1 + m2 `2 cos(1 2 )2 + m2 `2 sin(1 2 )22 + (m1 + m2 )g sin 1 = 0 , (7)
Small Oscillations
When the amplitudes of oscillation are small, the equations of motion can be linearized to find the normal
modes:
(m1 + m2 )`1 1 + m2 `2 2 + (m1 + m2 )g1 = 0 , (9)
`2 2 + `1 1 + g2 = 0 . (10)
Normal mode solutions have the form:
1 (t)
= 10 eit . (11)
2 (t) 20
The normal mode frequencies are determined by solving the secular equation
m1 (g 2 `1 ) m2 `2 2
=0. (12)
`1 2 (g 2 `2 )
2
The solutions of the secular equation are:
s
2
`1 `2
m1 + m2 (m1 + m2 ) m2 + m1 `1 +`2
2 g(`1 + `2 )
= . (13)
2`1 `2
m1
Hamiltons Equations
To derive the Hamiltonian equations for the double pendulum, start from the Lagrangian
1 1
L=T V = m1 v12 + m2 v22 m1 gy1 m2 gy2 , (19)
2 2
where T is the kinetic energy, and V is the potential energy, and g is the acceleration of gravity.
The generalized momenta can be found from the Lagrangian:
L
p1 = = (m1 + m2 )`21 1 + m2 `1 `2 2 cos(1 2 ) , (20)
1
L
p 2 = = m2 `22 2 + m2 `1 `2 1 cos(1 2 ) . (21)
2
H(1 , 2 , p1 , p2 ) = 1 p1 + 2 p2 L . (22)
3
The velocities can be expressed as functions of the coordinates and momenta:
`2 p1 `1 p2 cos(1 2 )
1 = , (23)
`1 `2 m1 + m2 sin2 (1 2 )
2
The Hamiltonian is
m2 `22 p21 + (m1 + m2 )`21 p22 2m2 `1 `2 p1 p2 cos(1 2 )
H =
2`21 `22 m2 m1 + m2 sin2 (1 2 )
Hamiltons equations for the time rate of change of the generalized momenta are
H K
p1 = = (m1 + m2 )g`1 sin(1 ) , (26)
1 1
H K
p2 = = m2 g`2 sin(2 ) , (27)
2 2
where the generalized kinetic energy
m2 `22 p21 + (m1 + m2 )`21 p22 2m2 `1 `2 p1 p2 cos(1 2 )
K H V = . (28)
2`21 `22 m2 m1 + m2 sin2 (1 2 )
Note that K depends on the generalized coordinates only through the difference 1 2 . Hence
K K m2 `22 p21 + (m1 + m2 )`21 p22 2m2 `1 `2 p1 p2 cos(1 2 )
= = 2 sin[2(1 2 )]
2 1 2`21 `22 m1 + m2 sin2 (1 2 )
p p sin(1 2 )
1 2
`1 `2 m1 + m2 sin2 (1 2 )
m K sin[2(1 2 )] p p sin(1 2 )
= 2 2
1 2 . (29)
m1 + m2 sin (1 2 ) `1 `2 m1 + m2 sin2 (1 2 )
4
Poincare Section
Because the phase space of the system is 4-dimensional and hard to visualize, it is useful to look at a
Poincare section, for example, the intersections of the trajectory with a 2-dimensional slice of the phase
space.
Because the Hamiltonian H is conserved in time, the trajectories actually lie on a 3-dimensional
hypersurface determined by the equation
C++ Program
The following program generates the pendulum trajectory and Poincare section data.
Program 1: http://www.physics.buffalo.edu/phy410-505/topic4/dpendulum.cpp
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <vector>
using namespace std;
#include "linalg.hpp"
#include "odeint.hpp"
using namespace cpl;
double t = x[0], theta_1 = x[1], theta_2 = x[2], p_1 = x[3], p_2 = x[4];
5
cos12 = cos(theta_1 - theta_2),
sin212 = sin(2 * (theta_1 - theta_2)),
denom = 1 + sin12 * sin12;
vector<double> f(5);
f[0] = 1;
f[1] = (p_1 - p_2 * cos12) / (m * L * L * denom);
f[2] = (2 * p_2 - p_1 * cos12) / (m * L * L * denom);
f[3] = - 2 * m * g * sin(theta_1) - C_1 + C_2;
f[4] = - m * g * L * sin(theta_2) + C_1 - C_2;
if (switch_t_and_theta_2) {
double d_theta_2_dt = f[2];
f /= d_theta_2_dt; // change variable from t to theta_2
}
return f;
}
double t = x[0], theta_1 = x[1], theta_2 = x[2], p_1 = x[3], p_2 = x[4];
double sin12 = sin(theta_1 - theta_2),
cos12 = cos(theta_1 - theta_2);
double a = 2 / denom,
b = - 2 * p_1 * cos12 / denom,
c = p_1 * p_1 / denom - E
+ m * g * L * (3 - 2 * cos(theta_1) - cos(theta_2));
6
double discr = b * b - 4 * a * c;
if (discr < 0) {
cerr << " Sorry, E = " << E << " too small, cannot solve for p_theta_2"
<< endl;
exit(EXIT_FAILURE);
}
int main() {
double t = 0;
vector<double> x(5); // extended solution vector
x[0] = t; x[1] = theta_1; x[2] = theta_2; x[3] = p_1; x[4] = p_2;
RK4 rk4;
double dt = 0.001;
double accuracy = 1e-6;
rk4.set_step_size(dt);
rk4.set_accuracy(accuracy);
int crossing = 0;
while (true) {
rk4.adaptive_step(f, x);
if (theta_2_old < 0.0 && x[2] >= 0.0 && x[4] >= 0.0) {
++crossing;
switch_t_and_theta_2 = true;
double dt_save = rk4.get_step_size();
vector<double> x_map = x;
7
double d_theta_2 = -x_map[2];
rk4.set_step_size(d_theta_2);
rk4.step(f, x_map);
switch_t_and_theta_2 = false;
rk4.set_step_size(dt_save);
cout << "Crossing " << crossing << " theta_1 = " << x_map[1]
<< " p_theta_1 = " << x_map[3] << " Energy = "
<< Energy(x_map) << endl;
section_file << x_map[1] << \t << x_map[3] << \n;
}
trajectory_file.close();
section_file.close();
}
Fig. 2 shows two trajectories of the double pendulum with the same energy: one set of initial conditions
results in a multiply periodic trajectory; a different set initial conditions lead to chaos!
2
1
p
-2
-4
-6
-1.5 -1 -0.5 0 0.5 1 1.5
1
Figure 2: Poincare sections with 5000 points for the double pendulum with E = 15. The chaotic trajectory
(red dots) results from choosing initial conditions 1 = 0, 2 = 0, p1 = 0, and the periodic trajectory (green
dots) from 1 = 1.1, 2 = 0, p1 = 0.
8
Homework Problem
First consider small oscillations. Simulate the two normal modes of the double pendulum: measure the
normal mode frequencies and verify the mode behavior by plotting the trajectories 1,2 (t).
Next, explore the full nonlinear behavior by generating Poincare sections for energy E = 1, 5, 10, 15, 40. Try
different initial conditions for each value of the energy E. For example, with E = 15, try the two different
initial values sets (1) 1 = 1.1, and 2 = 0, and (2) 1 = 0, and 2 = 0. Describe qualitatively the types of
motion you observe.
References
[1] T. Shinbrot, C. Grebogi and J. Wisdom, Chaos in a double pendulum, Am. J. Phys. 60, 491 (1992),
http://dx.doi.org/10.1119/1.16860.
[2] R.B. Levien and S.M. Tan, Double pendulum: An experiment in chaos, Am. J. Phys. 61, 1038
(1993), http://dx.doi.org/10.1119/1.17335.
[3] Double pendulum Java applet from a previous course:
http://www.physics.buffalo.edu/gonsalves/phy410-505_fall00/Chapter6/oct27.html.