The Double Pendulum

MATH1902: Numerical Solution of Differential Equations

http://people.sc.fsu.edu/∼jburkardt/classes/math1902 2020/double pendulum/double pendulum.pdf

We track the angles θ_{1} and θ_{2}; masses and lengths are parameters

The double pendulum The double pendulum system can exhibit sensitivity and chaotic solutions.

1 The nonlinear pendulum model

As a warmup, let us begin with the nonlinear pendulum model. The unknown variables are θ and \frac{dθ}{dt}, the deflection angle in radians, and angular velocity in radians/second. For convenience, we will refer to these two variables as y_{1} and y_{2} respectively. Newton’s laws allow us to write the following equations:

y′_{1} = y_{2}
y′_{2} = −\frac{g}{l} sin(y_{1})

where time is measured in seconds, the pendulum length l in meters, and g = 9.8 meters/seconds is the gravitational constant.

A typical initial condition at time t_{0} = 0 might be:

y_{1}(t_{0}) = \frac{π}{6}
y_{2}(t_{0}) = 0

This starting condition means that the maximum deflection of the pendulum will be \frac{π}{6} ≈ 0.5236. In this nonlinear model, the energy is:

E = m g l (1.0 − cos(y_{1})) + 0.5 m y_{2}^{2}

2 The linear pendulum model

Our initial pendulum deflection is relatively small, and the initial angular velocity is zero. This implies that the pendulum will never rise higher than this initial value. Note further that sin(\frac{π}{6} ) = 0.5, in other words, for this problem, the approximation sin(θ) ≈ θ will be reasonably close for all the values of angular deflection that we will encounter. It is this fact that allows us to approximate our nonlinear ODE by the following linear pendulum model, replacing sin(y_{1}) by y_{1}:

y′_{1} = y_{2}
y′_{2} = −\frac{g}{l} y_{1}

The ODE model is linear because both right hand sides are linear in the y variables. This model is less accurate than the nonlinear model, because we’ve made an extra approximation. However, it’s beloved by mathematicians because it’s exactly solvable. For our initial condition, we see that the linear model has the exact solution:

y_{1}(t) = \frac{π}{6} cos(ωt)
y_{2}(t) = −ω \frac{π}{6} sin(ωt)

where, for convenience, we define ω = \sqrt{\frac{g}{l}}. This implies that our solution functions will be periodic with period p = \frac{2π}{ω}. For our problem, ω ≈ 3.1337 and p ≈ 2.0050.

Thus any solution of the linear pendulum problem will involve simple harmonic motion, a combination of sine and cosine functions with period \frac{2π}{ω}, and so a phase plane plot of the solution will be an ellipse (a circle if ω = 1).

In this linear model, the energy is:

E = 0.5 m g l y_{1}^{2} + 0.5 m y_{2}^{2}

3 Comparing two pendulum models

Because they are simplifications of reality, all models are wrong. We know one way in which our linear pendulum model is wrong, namely in the use of the approximation sin(θ) ≈ θ, which gives us a mathematical problem for which we can find an exact solution. We expect that our linearized model is acceptable as long as the oscillations are small. But how much have we given up by this change? It’s easy to answer that question by setting up the MATLAB code for the nonlinear pendulum problem, running both models with a fairly large initial deflection, and comparing.

Our nonlinear pendulum derivative function is simply:

1 function dydt = p e n d u l u m n o n l i n e a r d e r i v ( t , y )

2

3  [ g , l , m, t0 , y0 ] = pendulum parameters ( ) ;

4

5  u = y ( 1 ) ;
6  v = y ( 2 ) ;

7

8  dudt = v ;
9  dvdt = − ( g / l ) ∗ sin ( u ) ; % s i n ( u ) , not u !

10

11  dydt = [ dudt ; dvdt ] ;

12

13  return
14 end

Listing 1: Nonlinear pendulum derivative

and for our initial condition, let’s go \frac{9}{10} of the way to the top, where θ = π:

y_{1}(t_{0}) = 0.9 ∗ π
y_{2}(t_{0}) = 0

If we plot θ and \frac{dθ}{dt}, we can see that the nonlinear pendulum (blue) still has a regular repeating behavior (but it’s not actually a sine or cosine curve!). Instead of the linear period of about p ≈ 2.0050, we see a longer period of about p ≈ 4. Moreover, the shape of the \frac{dθ}{dt} curve is remarkably different for the nonlinear case.

Comparison of solutions for linear and nonlinear pendulums

Moreover, the phase plot shows that the nonlinear pendulum (blue) is not following a circular or elliptical path, but some other shape entirely, and the energy value is very different.

Comparison of solutions for linear and nonlinear pendulums

If we had done this comparison for a problem with a small initial deflection, the two models would have been very close. However, assuming that the nonlinear model is closer to reality, we can see that the linear model is only reliable for small oscillations of the pendulum. It’s easy to perform physical experiments that show that the nonlinear model agrees with reality, aside from friction effects.

4 The double pendulum model

Now we will proceed to a related problem which, at first sight, seems only “slightly” more complicated than our previous cases that we understand reasonably well. We start with one pendulum, and attach another on the end of the first.

To be precise, in the double pendulum problem, a rod of length l_{1} is fixed at one end (0,0), and forms an angle θ_{1} with the downward vertical, so that its endpoint is at (x_{1}, y_{1}) = (l_{1} sin(θ_{1}), l_{1} cos(θ_{1})). (Remember that we measure angles from the downward vertical direction now!) A weight of mass m_{1} is attached to the end of the first rod.

A second rod, of length l_{2} is also attached to this end of the first rod, and can swing freely at the attachment point. This second rod forms an angle θ_{2} with the downward vertical. A weight of mass m_{2} is attached to the free end of the second rod. The position of this weight is (x_{2}, y_{2}) = (x_{1}, y_{1}) + (l_{2} sin(θ_{2}), l_{2} cos(θ_{2})).

Gravity has a force coefficient of g.

The unknowns are the angles θ_{1}, θ_{2} and their time derivatives. Starting from some initial condition, we wish to predict the locations and angular velocities of the two masses for future time.

We store the unknowns in a vector z:

z = [ θ_{1} ; \frac{dθ_{1}}{dt} ; θ_{2} ; \frac{dθ_{2}}{dt} ]

The differential equations become:

z′_{1} = z_{2}
z′_{2} = −\frac{((g(2m_{1} + m_{2}) sin(z_{1}) + m_{2}(g sin(z_{1} − 2z_{3}) + 2(l_{2}z_{4}^{2} + l_{1}z_{2}^{2} cos(z_{1} − z_{3})) sin(z_{1} − z_{3})))}{2l_{1}(m_{1} + m_{2} − m_{2} cos^{2}(z_{1} − z_{3}))}
z′_{3} = z_{4}
z′_{4} = \frac{(((m_{1} + m_{2})(l_{1}z_{2}^{2} + g cos(z_{1})) + l_{2}m_{2}z_{4}^{2} cos(z_{1} − z_{3})) sin(z_{1} − z_{3}))}{l_{2}(m_{1} + m_{2} − m_{2} cos^{2}(z_{1} − z_{3}))}

These derivatives can be evaluated by calling the following function, which is available on the web page.

1 dzdt = double pendulum deriv ( t , z ) ;

Listing 2: double pendulum deriv.m evaluates the ODE.

The system has 5 parameters. To start with, we will assign them the following default values:

• g = 9.81\frac{meters}{second^{2}}

• m_{1} = 1 kilogram;

• m_{2} = 1 kilogram;

• l_{1} = 2 meters;

• l_{2} = 1 meter;

These values of the parameters are available in the following function.

1 [ g , m1, m2, l1 , l 2 ] = double pendulum parameters ( ) ;

Listing 3: double pendulum parameters.m returns parameters.

Time t will be measured in seconds, and angles θ in radians. If you want to run the problem with different values for some of the parameters, it is enough to change the values in double pendulum parameters().

5 The system energy

A physical system has an energy. A model of the physical system can report its estimate of this energy at any time. If the model does not include any friction or other dissipative factors, then the total energy should be constant over time. It’s important to determine whether a model does express this property, and whether an ODE solver for the model accurately preserves the initial energy value.

We know that the midpoint method is one example of an ODE solver that can exactly conserve energy, if that energy can be expressed as a quadratic expression in terms of the variables. As it turns out, the energy of the double pendulum system is “mostly” quadratic, but includes trigonometric terms which we cannot assume the midpoint method will be able to capture perfectly. On the other hand, many other ODE solvers don’t make any claim to energy conservation, so it will be interesting to compute the energy for various solvers and see which ones do best.

Here is the awful expression for the energy E of the double pendulum, for a given set of variables z:

E = 0.5 (m_{1} + m_{2}) l_{1}^{2} z_{2}^{2}
+ 0.5 m_{2} l_{2}^{2} z_{4}^{2}
+ m_{2} l_{1} l_{2} z_{2} z_{4} cos(z_{1} − z_{3})
− (m_{1} + m_{2}) g l_{1} cos(z_{1})
− m_{2} g l_{2} cos(z_{3})

Once you have computed the solution values z over a stretch of time, and want to plot it, you can request the value of the energy from the following function:

1 e = double pendulum energy ( z ) ;

Listing 4: double pendulum energy.m returns the energy

In order for an energy plot to be meaningful, it is useful to also include a plot of the x axis, so that the plot has a scale. This can be done with commands like:

1  h = double pendulum energy ( z ) ; % Compute the energy .
2  hold ( ’ on ’ ) ;  % Allows m u l t i p l e p l o t s in one p i c t u r e
3  plot ( t , h ) ;  % Plot energy
4  plot ( [ t0 , t s t o p ] , [ 0 . 0 , 0 . 0 ] ) ; % Add p l o t o f x a x i s
5  hold ( ’ o f f ’ ) ;  % Finish the p i c t u r e .

Listing 5: Plotting energy + x axis.

6 Simulate the small perturbation problem with Forward Euler

Our first experiment will be timid. We will do a “small perturbation” problem, for which the double pendulum starts at time t_{0} = 0 with the initial condition z(t_{0}) = [0.25, 0, 0, 0]. Thus the first pendulum has been deflected to the small angle of θ = 0.25 radians = 0.25∗ \frac{180}{2π} ≈ 14.3^{∘} and the second pendulum is hanging straight down. We are interested in studying the system over the time period t_{0} = 0 ≤ t ≤ 50 = tstop. We start with our forward Euler code, available on the web page as euler.m.

The data will be numerous, and complicated. In order to have a chance of comprehending what is going on, we will need plots. The solution z is in terms of angles and angular velocities, but we would prefer to see the (x, y) coordinates of the two masses of the pendulum. So once we have computed z, we will convert from angles z_{1} = θ_{1}, z_{3} = θ_{2}, to Cartesian coordinates (x_{1}, y_{1}), (x_{2}, y_{2}). Because we are measuring angles starting from the vertical downward direction (“6 o’clock”) rather than the usual angular origin (“3 o’clock”), these formulas may not look like what you would first guess them to be:

1  x1 =  l 1 ∗ sin ( z1 ) ;
2  y1 =  − l 1 ∗ cos ( z1 ) ;
3  x2 = x1 + l 2 ∗ sin ( z3 ) ;
4  y2 = y1 − l 2 ∗ cos ( z3 ) ;

Listing 6: Converting angles to positions.

Then we will be interested in the following plots:

1. plot (t,x1(t)) and (t,y1(t)), the coordinates of mass m_{1} over time;

2. plot (t,x2(t)) and (t,y2(t)), the coordinates of mass m_{2} over time;

3. plot (x1(t),y1(t)), the position of m_{1};

4. plot (x2(t),y2(t)), the position of m_{2};

5. plot ( t, e(t) ), the energy E(t) over time;

For small perturbations, plots #1 and #2 should seem roughly periodic and roughly regular. For larger perturbations, where one or both pendulums can actually swing through the vertical upright position, we may expect much more complicated behavior. Plot #3 should be a circle, or a portion of a circle for small problems. Plot #4 should be made by circular paths looping along the circle of plot #1. Plot #5 would be a perfectly straight line if our ODE solver was exactly conservative. Instead, we can only hope that the line stays close to its initial value at t = 0. In order to judge whether the variations are significant, it’s helpful to draw a second line for e = 0, so that the plot has an implicit scale.

7 A code to call euler() for the double pendulum

The Forward Euler code uses a fixed stepsize. To get reasonable results over our time interval, it is necessary to use a very large number of steps. Here, we didn’t get believable pictures until trying n = 100, 000 steps. And this is for the small perturbation problem!

1  f = @ double pendulum deriv ;
2  t0 = 0 . 0 ;
3  t s t o p = 5 0 . 0 ;
4  tspan = [ t0 , t s t o p ] ;
5  z0 = [ 0 . 2 5 ; 0 . 0 ; 0 . 0 ; 0 . 0 ] ;
6  n = 100000;

7

8  [ t , z ] = e u l e r ( f , tspan , z0 , n ) ;
9 %
10 % Compute Cartesian c o o r d i n a t e s .
11 %
12  [ g , m1, m2, l1 , l 2 ] = double pendulum data ( ) ;

13

14  x1 = l 1 ∗ sin ( z ( 1 , : ) ) ;
15  y1 = − l 1 ∗ cos ( z ( 1 , : ) ) ;
16  x2 = x1 + l 2 ∗ sin ( z ( 3 , : ) ) ;
17  y2 = y1 − l 2 ∗ cos ( z ( 3 , : ) ) ;
18 %
19 % Plot #1.
20 %
21  figure ( 1 ) ;
22  plot ( t , x1 , ’ g ’ , …
23  t , y1 , ’ r ’ , …
24  ’ LineWidth ’ , 2 ) ;
25  grid ( ’ on ’ ) ;
26  xlabel ( ’ Time ’ ) ;
27  ylabel ( ’ x1 ( t ) , y1 ( t ) ’ ) ;
28  legend ( ’ x1 ’ , ’ y1 ’ ) ;
29  t i t l e ( ’ double pendulum ( x1 , y1 ) ’ )
30  f i l e n a m e = ’ double pendulum x1y1 . png ’ ;
31  print ( ’−dpng ’ , f i l e n a m e ) ;
32  fprintf ( 1 , ’ Graphics saved as ”%s ”\n ’ , f i l e n a m e ) ;
33 %
34 % More p l o t s …
35 %

Listing 7: Solving double pendulum with euler method.

8 Graphic results for the small perturbation problem

The code double pendulum test.m can use the Euler method to solve our small perturbation problem, and produces plots that summarize the results. Because the pendulum is only slightly swinging, the two components of the pendulum seem to be swinging together almost as one.

Time plots of (t, x_{1}), (t, y_{1}) and (t, x_{2}), (t, y_{2}) show regularity and periodicity

The phase plot shows that pendulum #1 stays on a circular arc, and #2 gently loops around it.below it; the energy plot shows conservation.

The graphic results from the small perturbation problem suggest that the solution can exhibit regularity and periodicity for both (x_{1}, y_{1}) and (x_{2}, y_{2}), the positions of the two pendulum masses. We see that the phase portrait is pretty calm, with the first pendulum tracing a thin line, while the thicker line for the second pendulum mass indicates that its motion is somewhat more complicated - it must actually be jiggling up and down just a bit.

9 Experiments

Our small perturbation seems to have allowed the double pendulum to behave almost as though it were still just a single pendulum. We need to shake things up a bit to see what the system is capable of. Let us consider three sets of initial conditions:

1. small: z = [0.25, 0, 0, 0];

2. medium: z = [π/3, 0, 0, 2];

3. large: z = [π, 0, π, 4];

We have already considered the small perturbation case, so it is up to you to run the medium and large perturbations, make the plots, and try to decide what they are telling you.

You may want to check whether the results would change significantly if you increase the number of steps (for a fixed step method) or tighten the error tolerance (for an adaptive method). If the energy conservation improves, for instance, then you might conclude that the first computation was simply not accurate enough. Otherwise, you might wonder whether this ODE system can be too wild to approximate accurately.

The double pendulum test.m file has some commented out portions that make experimentation easy. Simply by uncommenting certain lines, you can switch the ODE solver between euler(), rk4(), midpoint() or ode45(). Especially for the medium and large initial conditions, these solvers may return results that are significantly different.

The function double pendulum parameters() sets the values of g, l1, l2, m1, m2. Changing any value in this file, will cause double pendulum test() to use that new value instead of the ones we talked about initially. In particular, you might investigate what happens if

• you set l2 to 5, or to 1/5;

• you set m2 to 5, or to 1/5;

The initial condition has a huge influence on the solution. What happens if you choose an initial condition so that

• pendulum 1 is not moving: for example, z0 = [ 0, 0, 0, 5 ];

• pendulum 1 and pendulum 2 are at the same angle: for example: z0 = [ pi/3, 0, pi/3, 0 ]

• pendulum 1 and pendulum 2 are at opposite angles: for example: z0 = [ pi/3, 0, -pi/3, 0 ]

• the angular velocity is just large enough that the pendulum will reach the top, and keep going. (Then theta1 and theta2 will not necessarily stay between -pi and pi.)

There is a MATLAB program, written by Alexander Erlich, which can solve the double pendulum problem, display the pendulum system, and make a movie of the results. Watching the animation for various perturbation sizes, it is easy to see how the system can be regular and rhythmic for small perturbation, becoming chaotic when more energy is available, so that the pendulums can both swing full circle.

On the course website, there is a copy of the program, set up to solve the small perturbation problem, in the file double pendulum movie small.m, as well as the resulting animation double pendulum movie small.avi. To change the code to do the medium or large perturbation problems only requires changing the data in the initial segment of the code.

In that initial code, you can also easily see how to choose different values for the gravitational, mass, and length parameters, or to experiment with other initial conditions.

10 HW # 10

Choose some of the double pendulum experiments above, or set up your own. Write a short report on your results.

Send your report to trenchea@pitt.edu.

Professor Trenchea and I have enjoyed presenting this material to you. We hope it has shown you some of the ways that researchers can deal with mathematical problems by taking a computational approach.