Double Pendulum

The double pendulum is an interesting example in classical mechanics, as it shows what can happen if one simplifies a problem. Grant Fowles and George Cassiday in Analytical Mechanics (6th ed., Saunders College Publ. 1999) call it also the problem of two rock climbers dangling on a single rope. Many mechanics texts (and also Differential Equations with Maple V) discuss only the linearized case where one investigates the normal modes. However, the large-amplitude motion is interesting as it includes chaotic behaviour. The problem is dealt with in Andre Heck: Introduction to Maple (Springer 1993); he calls the problem the two-link robot arm. A perfect classroom demonstration can be made of a more sophisticated version, namely the double compound pendulum for which one needs a couple of plates and two ball bearings.

Fowles and Cassiday almost derive the full-blown double pendulum (they could have explained the one missing step easily). We pick the simplified choice of parameters (equal mass and equal length), and begin with a graph of the potential energy. The angles that the two massless bars make with the vertical are called theta and phi respectively. The first term below is the potential energy of the first mass, and the second the potential energy for the second mass which is attached to the first.

> restart;

> Vpot:=(theta,phi)->m*g*l*(1-cos(theta))+m*g*l*(2-cos(theta)-cos(phi));

> with(plots):

Warning, the name changecoords has been redefined

> contourplot(subs(g=10,m=1,l=1,Vpot(t,f)),t=0..2*Pi,f=0..2*Pi,filled=true,contours=[5,10,15,20,25,30,35,40,45,50,55]);

> #plot3d(subs(g=10,m=1,l=1,Vpot(t,f)),t=0..2*Pi,f=0..2*Pi,axes=boxed,style=patchcontour);

To derive the kinetic energy expression one starts with the expression that involves the linear velocities v 1 and v 2 for the two masses. T = m /2 ( v 1^2+ v 2^2). It is important to realize that the second mass has a velocity given as v 1 + its own rotation about mass m 1. The speeds associated with the two rotations are given as length times respective angular velocity. The direction of the rotations about the respective support points is given by the theta and phi unit vectors. These have Cartesian components [cos(theta), sin(theta)] and [cos(phi), sin(phi)] respectively. Thus, in Cartesian coordinates in the plane of swinging we have v 1 = l diff(theta( t ), t )*[cos(theta( t )),sin(theta( t ))], and v 2 = v 1+ l diff(phi( t ), t )*[cos(phi( t ),sin(phi( t ))]. The kinetic energy reduces to:

> Tkin:=m*l^2*(diff(theta(t),t)^2+1/2*diff(phi(t),t)^2+diff(theta(t),t)*diff(phi(t),t)*(cos(theta(t))*cos(phi(t))+sin(theta(t))*sin(phi(t))));

In the linearized case the entire bracket with the trig functions is replaced by one.

We will now carry out the procedure of Euler and Lagrange. The Lagrange function is defined as:

> Lagr:=Tkin-Vpot(theta(t),phi(t));

To carry out the variation is not trivial in Maple as we have to differentiate with respect to a function. We can help ourselves with substitutions (suggested by A. Heck). An alternative is to download the calcvar package from the share library.

It is important to substitute the derivatives first; otherwise the dependent variables in the derivatives are renamed, and then are not recognized.

> L1:=subs(diff(theta(t),t)=thdot,theta(t)=th,diff(phi(t),t)=phdot,phi(t)=ph,Lagr);

Now we can do partial derivatives with respect to theta, or w.r.t. diff(theta(t),t), etc., but we cannot differentiate with respect to time without striking disaster.

So let us proceed slowly:

> dL_dth:=diff(L1,th);

> dL_dthdot:=diff(L1,thdot);

The latter has to be differentiated with respect to time. We reverse the substitutions (before taking the time derivative on the lhs), and collect the expressions:

> EL1:=diff(subs(thdot=diff(theta(t),t),th=theta(t),phdot=diff(phi(t),t),ph=phi(t),dL_dthdot),t)-subs(thdot=diff(theta(t),t),th=theta(t),phdot=diff(phi(t),t),ph=phi(t),dL_dth)=0;

We need to get some order into this...

> EL1:=collect(EL1/m,[diff(theta(t),t$2),diff(phi(t),t$2),diff(theta(t),t),diff(phi(t),t)]);

> EL1:=combine(EL1,trig);

This is still a formidable equation, but it is readable at least. It is hard to keep the combined trig functions, and to divide out the length at the same time.

We repeat the tricks to obtain the second equation.

> dL_dph:=diff(L1,ph);

> dL_dphdot:=diff(L1,phdot);

> EL2:=diff(subs(thdot=diff(theta(t),t),th=theta(t),phdot=diff(phi(t),t),ph=phi(t),dL_dphdot),t)-subs(thdot=diff(theta(t),t),th=theta(t),phdot=diff(phi(t),t),ph=phi(t),dL_dph)=0:

> EL2:=collect(EL2/m,[diff(theta(t),t$2),diff(phi(t),t$2),diff(theta(t),t),diff(phi(t),t)]):

> EL2:=combine(EL2,trig);

The big question is whether Maple's dsolve[numeric] engine can handle this pair of second-order equations directly.

We pick initial conditions as one pendulum on top of the other (unstable equilibrium):

> eta:=10^(-2);

> IC:=theta(0)=Pi,D(theta)(0)=-eta,phi(0)=Pi,D(phi)(0)=0;

We need to pick values for the constants (e.g., in MKSA or SI units):

> EL1p:=subs(g=10,l=1,EL1);

> EL2p:=subs(g=10,l=1,EL2);

> sol:=dsolve({EL1p,EL2p,IC},{theta(t),phi(t)},numeric,method=lsode,output=listprocedure):

> th:=subs(sol,theta(t)): ph:=subs(sol,phi(t)):

We have to delay evaluation of the variables until actual numbers are passed for the time parameter, thus the quotes below:

> plot(['th(t)','ph(t)'],t=0..30,color=[red,blue],numpoints=1000);

Obviously, the pendulum goes wild...

1) Note the reluctance to get away from the unstable equilibrium at the beginning.

2) Note that the lower pendulum (blue curve) goes into a winding mode at about t =10, and later unwinds [for initial conditions theta=Pi, phi=0].

3) The numpoints option was chosen to avoid weird results.

4) The variables theta(t) and phi(t) keep track of the so-called winding number.

5) After t=15 the lower pendulum seems to be gaining energy (steepening slope) at the expense of the upper one [for initial conditions theta=Pi, phi=0].

Can we animate the solutions? For this we need to reconstruct the Cartesian coordinates from the two angles.

We pick for the lengths l =1.

> x1:=t->sin(th(t)):

> y1:=t->-cos(th(t)):

> x2:=t->x1(t)+sin(ph(t)):

> y2:=t->y1(t)-cos(ph(t)):

Our animations will be made up from primitive straight-line connections.

> PL:=t->plot([[0,0],[x1(t),y1(t)],[x2(t),y2(t)]],style=line,color=magenta,title=cat("time= ",convert(evalf(t,3),string))):

> PLa:=seq(PL(i/10),i=1..300):

> display(PLa,insequence=true,scaling=constrained,axes=boxed);

Observe how there are moments when the upper pendulum is almost stationary at the bottom, while the lower one is spinning at a maximum rate.

Note that the animation can be stopped and continued later, so that one can go back to the upper graph to compare the information.

Exercise1: Change the intial conditions by a small amount and compare time histories. Change the method of integration (remove method=lsode to use a Runge-Kutta method), and observe whether the results are being reproduced for the same initial conditions.

Exercise2: Extend the problem by allowing different masses for the top and bottom pendula (pendulums).

>