The Chaotic Motion of A Double Pendulum
The Chaotic Motion of A Double Pendulum
The Chaotic Motion of A Double Pendulum
Carl W. Akerlof
September 26, 2012
The following notes describe the kinematics of the double pendulum. The starting point is a
pendulum consisting of two point masses, m, and m2, suspended by massless wires of length l1
and l2. The treatment of this case can be found at:
http://scienceworld.wolfram.com/physics/DoublePendulum.html
For a real system, the equations of motion depend in a more complicated way on the distribution
of mass that is essential for modeling the physical pendulum used in this experiment.
Figure 1. Point mass double pendulum. Figure 2. Extended mass double pendulum.
x1 l1 sin 1
y1 l1cos 1
x2 l1 sin l2 sin 2
y2 l1cos 1 l2 cos 2
x1 u1 sin 1 v1 cos 1
y1 u1 cos 1 v1 sin 1
x2 l1 sin 1 u2 sin 2 v2 cos 2
y2 l1 cos 1 u2 cos 2 v2 sin 2
KE
1
2 m r
1 1
2
12 m2l12 12 2 m2 u2 l1 cos 2 1 m2 r22 22
PE m1 u1 g cos 1 m2l1 g cos 1 m2 u2 g cos 2
L T V
1
2 m r
1 1
2
12 m2l1212 2 m2 u2 l1 cos 2 1 12 m2 r22 22
m1 u1 g cos 1 m2l1 g cos 1 m2 u2 g cos 2
This leads directly to the equations of motion which we shall investigate shortly.
The aim of this experiment is to compare the actual dynamical behavior of a real physical
pendulum with a mathematical simulation. To this end, we need to characterize the pendulum
properties as accurately as possible and incorporate these into the appropriate equations of
motion.
One approximation will be involved: the motion will be assumed frictionless. This
simplification is driven principally by the lack of any very elegant fundamental theory although
it would actually be fairly trivial to incorporate velocity-dependent damping in the dynamical
modeling.
The vendor for this apparatus, chaoticpendulums.com, has conveniently provided the dimensions
for these parts as shown in Figure 3. Note that all dimensions are in millimeters. The required
spatial mass moments can be calculated analytically for each item. The contribution that each
one provides to the kinetic and potential energies is given below. For future reference, the
distance between the primary and secondary axes is denoted l1 and equals 173 mm.
Bearing: A bearing facilitates low-friction axial rotation shear while constraining the shaft
location. Thus, the inner race of the bearing can be rotating with angular velocity, , while the
outer race is rotating with an angular velocity of . To keep the modeling of this component
simple, it is assumed that the bearing is homogeneous with a rotation rate at radius, r, that is
linearly interpolated by the distances from the inner and outer radii.
KE
1
2 m l
2 2
b 1 1 mb r 2
12 2 mb r 2
12 mb r 2
22
PE mbl1 g cos 1
Primary arm:
1
KE m1 r 2 12
2
PE m1 u g cos 1
The cap screw, spacers and hex nut are all constrained to rotate rigidly with the secondary arm.
KE
1
2 m r
c
2
2 ms r 2 mn r 2 22
PE mc 2 ms mn l1 g cos 1
Secondary arm:
KE
1
2 m l
2 2
2 1 1 2 m2 u l112 cos 2 1 m2 r 2 22
PE m2 l1 g cos 1 m2 u g cos 2
where:
a11 2 m1 r 2 4 mb r 2 m2 2mb mc 2ms mn l12
a12 2 mb r 2
a22 m2 r 2
2 mb r 2 mc r 2 2 ms r 2 mn r 2
b12 m2 u l1
c1 2 m1 u m2 2 mb mc 2 ms mn l1 g
c2 m2 u g
H p11 p22 L
a1112 a12 b12 cos 2 1 12 a2222 c1 cos 1 c2 cos 2
1 1
2 2
By replacing i by pi , one can convert the Hamiltonian into a form that depends only on i and
pi .
a p bp2
1 22 1
a11 a22 b 2
bp1 a11 p2
2
a11 a22 b 2
and
b a12 b12 cos 2 1
Thus:
1 a22 p1 2 a12 b12 cos 2 1 p1 p2 a11 p2
2 2
H c1 cos 1 c2 cos 2
a11 a22 a12 b12 cos 2 1
2
2
H
p1 c1 sin 1 q 1 , 2 , p1 , p2
1
H
p 2 c2 sin 2 q 1 , 2 , p1 , p2
2
q 1 , 2 , p1 , p 2
a
12
b12 cos 2 1 p1 a11 p2 a22 p1 a12 b12 cos 2 1 p2 b12 sin 2 1
In this experiment, we would like to explore how sensitively the dynamical evolution of a system
depends on the initial conditions. Do small perturbations lead to small variations or huge ones?
As a starting point, consider a one-dimensional system whose evolution is governed by:
q f q t
Suppose we want to find out how q will evolve for slightly different initial conditions in the
neighborhood of qo:
f
q f q t f qo t q
q
q f
q(t ) et q 0
q q qo
The value, λ, is the Lyapunov exponent at qo. If it is greater than zero, the evolution of the
system diverges exponentially with time, preventing all attempts to compute the motion at
arbitrary times.
We can take this example to the next level of complexity by considering the simple pendulum.
The Lagrangian for this system is:
1
L T V m r 2 2 m u g cos
2
L
p m r 2
We would like to find out how the motion evolves for small perturbations from the initial
condition specified by θ0 and p0.
2 H 2 H p
2 p
p p m r2
2 H 2 H
p p m u g cos
2 p
1
0
m r 2
p m u g cos p
0
where q (0) a bp is an eigenvector for the matrix with an eigenvalue, λ, and thus
q(t ) et q(0) . To find λ, the following equation must be solved:
m m21
Det 11 0
m12 m22
m u cos
2 g
m r2
m u g cos
: i
2 m r2
1
e
i m u m r2 g cos
m u g cos
:
2 m r2
1
e
m u m r g cos
2
To understand the long term behavior of the system, one must average the Lyapunov exponents
over the dynamical path in the phase space.
It is worth remarking about a generic property of Lyapunov exponents stemming from the
Hamiltonian description. The Jacobian matrix that transforms the phase space volume,
1 p1 2 p2 must have the form:
2 H 2 H 2 H 2 H
p
1 1 p12 p1 2 p1p2
2 H 2 H 2H 2 H
J 2
1 1p1 2 2 1p2
The diagonal elements of this matrix cancel in pairs so that the trace of J is zero. This guarantees
that the sum of eigenvalues, ie. the Lyapunov exponents, is also zero. The physical implication
of this is the invariance of phase space volume for a conservative system. For a dissipative
system, this would not be true.
By defining the following intermediate quantities, the matrix elements can be computed fairly
easily. (See previous parameter definitions.)
j11
s p2 2 a22 u b12 sin 2 1
s2
a22
j12
s
j13
s p2 2 r t b12 sin 2 1
s2
r
j14
s
j21 c1 cos 1 dq
j22 j11
j23 dq
j24 j33
j31 j33
j32 j14
In general, the characteristic equation defining the eigenvalues would be fourth-order in λ. The
traceless nature of J requires the cubic term in λ to vanish and similar other symmetries imposed
by the Hamiltonian origin squelch the linear term as well. Thus, the defining equation for the
eigenvalues takes the form:
2 2
2 b 2 c 0
and thus:
2 b b 2 c
where:
b 1
2
j
2
11 j332 j12 j21 j34 j43 j23 j32 j24 j42
c Det J
Solenoid releases: Two solenoid assemblies are available to hold each arm of the double
pendulum at a fixed predetermined location prior to release. The plungers must be pulled out
manually to hold the arms. Note that the two solenoid assemblies are not quite identical – each
Electronic flash gun: The angular coordinates of the double pendulum arms are captured at
preset times by a short duration burst from a LumoPro LP120 manual flash unit. This device was
chosen because the vendor was reasonably willing to provide information about trigger
requirements. The unit is battery-powered so the appropriate switches must be turned on to
activate. Set intensity switches to NEXT and 1/8 to optimize flash duration and brightness. Please
also make sure to switch unit off before you leave. The LP120 takes about 10 seconds or so to
charge as indicated by a small red LED near the on/off switch. If longer, swap rechargeable
batteries.
Phosphorescent screen: With the bright light from the flash gun, the angular positions of the
double pendulum arms can be recorded by phosphorescent screens, either 12” ä 12” or 6” ä 6” in
area. This will only work in a darkened room. The image decay time is of the order of 30
seconds. This provides a reasonable opportunity to measure the shadows of the pendulum arms.
In case of any need for replacements, the screen material is ZnS vinyl film and was obtained
from Educational Innovations, Inc. and Hanovia Specialty Lighting.
Red flashlight: A Celestron night vision red LED flashlight is available to read the Wixey angle
gauge without bleaching the phosphorescent screen image. A rotating thumbwheel turns the
device on and adjusts the intensity of light. Turn off when not in use.
Leveling table: A small leveling table is available to help set the reference angle for the Wixey
digital gauge. Adjust the three leveling screws so that the 2” diameter steel ball has no tendency
Figure 9. Leveling table with 2” diameter Figure 10. Leveling table with steel block
steel ball. Adjust screws until ball position to set the Wixey digital angle gauge
is neutral. reference direction to vertical.
Validation tests
The physical parameters for the double pendulum used in this experiment are listed in Table I.
The various mass moments were computed analytically from the relevant dimensions. To check
whether these values adequately describe the system, one can constrain the two pendulum arms
in various ways and accurately measure the small amplitude periods. These correspond to the
following configurations:
a. Hold the primary arm fixed and measure the oscillation frequency of the secondary arm.
b. Hold the secondary arm fixed inside the primary arm using a #8 elastic band near the
primary axis.
c1 c2
b2
a11 2 a12 b12 a22
c. Hold the secondary arm fixed to and extended beyond the primary arm using two #8
elastic bands and a small rectangle of box board as a stiffener.
c1 c2
c2
a11 2 a12 b12 a22
The periods and thus the oscillation frequencies can be accurately measured with a Vernier
photogate interfaced to the computer.
Figure 12. Pendulum period Figure 13. Pendulum period Figure 14. Pendulum period
test configuration (a). test configuration (b). test configuration (c).
The previous measurements should validate the pendulum parameters to the accuracy of the
order of 0.1%. The next step is to use these values to predict the pendulum behavior over a finite
range of time of the order of 10 seconds. The most convenient numerical technique is Runge-
I have used the IDL numerical analysis package to compute and plot the behavior of this double
pendulum system as depicted in Figures 15, 16, and 17. This requires two or three dozen lines of
code. I expect that most students are more likely to have experience with MatLab. Program
scripts for both IDL and MatLab are available on the Physics 441/442 Web site. For IDL, you
will need d_chaos.pro, d_pend.pro and lyapunov.pro. MatLab requires d_chaos.m, d_pend.m,
lyapunov.m and rk4.m (The last file performs Runge-Kutta integration. This code is already
embedded in IDL.) In the course of this experiment, you will be required to perform additional
calculations to explore the chaotic motion but these routines should be a useful starting point.
Figure 15. Graphs of θ1(t), p1(t), θ2(t), p2(t) for θ1(0) = 10º, θ2(0) = 10º. The red, green and blue
lines show the behavior for three initial conditions for θ1 and θ2 that vary by 0.01 radians.
The main goal of this experiment is to examine the behavior of the double pendulum as a
function of initial conditions using both computational and experimental techniques. A useful set
of initial conditions to explore is (θ1 = nπ/12, p1 = 0, θ2 = nπ/12, p2 = 0) where n is an integer in
the range from 3 to 9. For small values of n, the experimentally measured pendulum positions
will track the calculations accurately and the dispersion of values at fixed time will be small. As
n gets larger, the differences between the measured and computed displacements will grow,
followed by increases in the dispersions at fixed times. The technique for these measurements is
to set up the electronic flash lamp to trigger at the desired time delay after solenoid release and
then immediately place the Wixey angle gauge and block on the phosphorescent screen to read
off the pendulum arm positions at the time of flash. (See Figures 18 & 19 below.) To determine
dispersion, measure θ1 and θ2 for ten identical releases from the same initial conditions and
compute the standard deviation. You need to be aware that for large values of n, the secondary
pendulum arm can easily have executed one more or one less complete rotation, making
dispersion estimates unreliable under these conditions. For motion in the chaotic regime, plot the
log of dispersion as a function of time.
Figure 18a. Wixey digital angle gauge and retainer block for Figure 18b. Angle gauge
measuring shadow angles on the phosphorescent screen. zero reference position.
The dispersion of phase space trajectories should be simultaneously explored with numerical
techniques. Vary a particular initial condition for (θ1, p1, θ2, p2) by (δθ1, 0, δθ2, 0) where δθ1 and
δθ2 are Gaussian distributed random errors with standard deviations of the order of 0.0005
radians. If you only have access to a uniform random number generator over the interval, [0 ≤ r ≤
1], use rGauss = r1 + r2 + … + r11 + r12 – 6. Generate a hundred or so trajectories to compute the
dispersion. Compare these results qualitatively with the corresponding experimental
measurements.
I noted that the experimentally measured dispersions were not always very stable in value, even
if the mean positions were quite reproducible. That may be an artifact of the strong possibility
that Gaussian initial errors do not propagate to Gaussian trajectory dispersions. This is an open
question for the moment that could be explored by computational methods.
Although the double pendulum does not give rise to interesting geometric behavior, the general
computational procedure often produces surprising results. The figure below was taken from a
recent New York Times Web blog by Steven Strogatz. The image depicts the behavior of
Newton’s method for iteratively computing the solution to z3 = 1. The color of each point in the
figure indicates which of the three roots will be approached from that particular point in the
complex plane.
Figure 20. A map of the complex plane with -1 ≤ Re(z) ≤ 1 and -1 ≤ Im(z) ≤ 1. The three colors
designate which complex root will be found by iteration of Newton’s method for the solution of
z3 = 1. See S. Strogatz.
a b 2 ab
a b 1 2
2 2
x y dx dy
0 0 3
r r 2 x2 r2
0
dx
0
(( x0 x)2 ( y0 y )2 ) dy
24
(16 ( x0 y0 ) r 3 r 2 2 ( x02 y02 ))
a a 2 b2 b
a b i (c d i ) ; d , c
2 2d
References
1. Troy Shinbrot, Celso Grebogi, Jack Wisdom and James A. Yorke, Chaos in a double
pendulum, Am. J. Phys. 60, 491-499 (1992).
2. Tomasz Stachowiak and Toshio Okada, A numerical analysis of chaos in the double
pendulum, Chaos, Solitons and Fractals 29, 417–422 (2006).
3. Steven Strogatz, Finding Your Roots, New York Times, March 7, 2010.