{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 1 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 0 1 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
278 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 
1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 15 "Double Pendulum" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 163 "The double pen
dulum is an interesting example in classical mechanics, as it shows wh
at can happen if one simplifies a problem. Grant Fowles and George Cas
siday in " }{TEXT 257 20 "Analytical Mechanics" }{TEXT -1 381 " (6th e
d., Saunders College Publ. 1999) call it also the problem of two rock \+
climbers dangling on a single rope. Many mechanics texts (and also Dif
ferential Equations with Maple V) discuss only the linearized case whe
re one investigates the normal modes. However, the large-amplitude mot
ion is interesting as it includes chaotic behaviour. The problem is de
alt with in Andre Heck: " }{TEXT 278 21 "Introduction to Maple" }
{TEXT -1 239 " (Springer 1993); he calls the problem the two-link robo
t arm. A perfect classroom demonstration can be made of a more sophist
icated version, namely the double compound pendulum for which one need
s a couple of plates and two ball bearings." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 496 "Fowles and Cassiday almost der
ive the full-blown double pendulum (they could have explained the one \+
missing step easily). We pick the simplified choice of parameters (equ
al mass and equal length), and begin with a graph of the potential ene
rgy. The angles that the two massless bars make with the vertical are \+
called theta and phi respectively. The first term below is the potenti
al energy of the first mass, and the second the potential energy for t
he second mass which is attached to the first." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "Vpot:=(theta,phi)->m*g*l*(1
-cos(theta))+m*g*l*(2-cos(theta)-cos(phi));" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 117 "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]);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "#plot3d(subs(g=10,m=1,l=1,Vp
ot(t,f)),t=0..2*Pi,f=0..2*Pi,axes=boxed,style=patchcontour);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 107 "To derive the kinetic energy expr
ession one starts with the expression that involves the linear velocit
ies " }{TEXT 260 1 "v" }{TEXT -1 6 "1 and " }{TEXT 261 1 "v" }{TEXT 
-1 26 "2 for the two masses. T = " }{TEXT 262 1 "m" }{TEXT -1 4 "/2 (
" }{TEXT 263 1 "v" }{TEXT -1 4 "1^2+" }{TEXT 264 1 "v" }{TEXT -1 78 "2
^2). It is important to realize that the second mass has a velocity gi
ven as " }{TEXT 258 1 "v" }{TEXT -1 32 "1 + its own rotation about mas
s " }{TEXT 259 1 "m" }{TEXT -1 374 "1. The speeds associated with the \+
two rotations are given as length times respective angular velocity. T
he direction of the rotations about the respective support points is g
iven by the theta and phi unit vectors. These have Cartesian component
s [cos(theta), sin(theta)] and [cos(phi), sin(phi)] respectively. Thus
, in Cartesian coordinates in the plane of swinging we have " }{TEXT 
265 1 "v" }{TEXT -1 4 "1 = " }{TEXT 266 1 "l" }{TEXT -1 12 " diff(thet
a(" }{TEXT 267 1 "t" }{TEXT -1 2 ")," }{TEXT 268 1 "t" }{TEXT -1 13 ")
*[cos(theta(" }{TEXT 277 1 "t" }{TEXT -1 13 ")),sin(theta(" }{TEXT 
276 1 "t" }{TEXT -1 9 "))], and " }{TEXT 270 1 "v" }{TEXT -1 3 "2 =" }
{TEXT 269 2 " v" }{TEXT -1 3 "1+ " }{TEXT 271 1 "l" }{TEXT -1 9 "diff(
phi(" }{TEXT 272 1 "t" }{TEXT -1 2 ")," }{TEXT 273 1 "t" }{TEXT -1 11 
")*[cos(phi(" }{TEXT 274 1 "t" }{TEXT -1 10 "),sin(phi(" }{TEXT 275 1 
"t" }{TEXT -1 35 "))]. The kinetic energy reduces to:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 140 "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))));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 85 "In the linearized case the entire bracket with the trig functio
ns is replaced by one." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT -1 95 "We will now carry out the procedure of Euler and Lag
range. The Lagrange function is defined as:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 33 "Lagr:=Tkin-Vpot(theta(t),phi(t));" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 245 "To carry out the variation is not trivia
l in Maple as we have to differentiate with respect to a function. We \+
can help ourselves with substitutions (suggested by A. Heck). An alter
native is to download the calcvar package from the share library." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 147 "It is im
portant to substitute the derivatives first; otherwise the dependent v
ariables in the derivatives are renamed, and then are not recognized.
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "L1:=subs(diff(theta(t),
t)=thdot,theta(t)=th,diff(phi(t),t)=phdot,phi(t)=ph,Lagr);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 166 "Now we can do partial derivatives with r
espect to theta, or w.r.t. diff(theta(t),t), etc., but we cannot diffe
rentiate with respect to time without striking disaster." }}{PARA 0 "
" 0 "" {TEXT -1 25 "So let us proceed slowly:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 20 "dL_dth:=diff(L1,th);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 26 "dL_dthdot:=diff(L1,thdot);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 163 "The latter has to be differentiated with respect to t
ime. We reverse the substitutions (before taking the time derivative o
n the lhs), and collect the expressions:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 176 "EL1:=diff(subs(thdot=diff(theta(t),t),th=theta(t),ph
dot=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;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 38 "We need to get some order into this..." }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "EL1:=collect(EL1/m,[diff(t
heta(t),t$2),diff(phi(t),t$2),diff(theta(t),t),diff(phi(t),t)]);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "EL1:=combine(EL1,trig);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 160 "This is still a formidable equati
on, 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." }}{PARA 0 "
" 0 "" {TEXT -1 51 "We repeat the tricks to obtain the second equation
." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "dL_dph:=diff(L1,ph);" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "dL_dphdot:=diff(L1,phdot)
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 176 "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=ph
i(t),dL_dph)=0:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "EL2:=col
lect(EL2/m,[diff(theta(t),t$2),diff(phi(t),t$2),diff(theta(t),t),diff(
phi(t),t)]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "EL2:=combine(EL2,tr
ig);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 115 "The big question is whet
her Maple's dsolve[numeric] engine can handle this pair of second-orde
r equations directly." }}{PARA 0 "" 0 "" {TEXT -1 86 "We pick initial \+
conditions as one pendulum on top of the other (unstable equilibrium):
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "eta:=10^(-2);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "IC:=theta(0)=Pi,D(theta)(0)=
-eta,phi(0)=Pi,D(phi)(0)=0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "We
 need to pick values for the constants (e.g., in MKSA or SI units):" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "EL1p:=subs(g=10,l=1,EL1);
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "EL2p:=subs(g=10,l=1,EL2
);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "sol:=dsolve(\{EL1p,EL
2p,IC\},\{theta(t),phi(t)\},numeric,method=lsode,output=listprocedure)
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "th:=subs(sol,theta(t))
: ph:=subs(sol,phi(t)):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 123 "We ha
ve to delay evaluation of the variables until actual numbers are passe
d for the time parameter, thus the quotes below:" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 64 "plot(['th(t)','ph(t)'],t=0..30,color=[red,blu
e],numpoints=1000);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "Obviously,
 the pendulum goes wild... " }}{PARA 0 "" 0 "" {TEXT -1 82 "1) Note th
e reluctance to get away from the unstable equilibrium at the beginnin
g." }}{PARA 0 "" 0 "" {TEXT -1 79 "2) Note that the lower pendulum (bl
ue curve) goes into a winding mode at about " }{TEXT 279 1 "t" }{TEXT 
-1 23 "=10, and later unwinds." }}{PARA 0 "" 0 "" {TEXT -1 58 "3) The \+
numpoints option was chosen to avoid weird results." }}{PARA 0 "" 0 "
" {TEXT -1 81 "4) The variables theta(t) and phi(t) keep track of the \+
so-called winding number. " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "Ca
n we animate the solutions? For this we need to reconstruct the Cartes
ian coordinates from the two angles." }}{PARA 0 "" 0 "" {TEXT -1 24 "W
e pick for the lengths " }{TEXT 280 1 "l" }{TEXT -1 3 "=1." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "x1:=t->sin(th(t)):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "y1:=t->-cos(th(t)):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "x2:=t->x1(t)+sin(ph(t)):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "y2:=t->y1(t)-cos(ph(t)):" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 72 "Our animations will be made up fr
om primitive straight-line connections." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 121 "PL:=t->plot([[0,0],[x1(t),y1(t)],[x2(t),y2(t)]],styl
e=line,color=magenta,title=cat(\"time= \",convert(evalf(t,3),string)))
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "PLa:=seq(PL(i/10),i=1.
.300):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "display(PLa,inseq
uence=true,scaling=constrained,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 140 "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." }}{PARA 0 "" 0 "" {TEXT -1 241 "Note that the animat
ion can be stopped and continued later, so that one can go back to the
 upper graph to compare the information. From Maple 6 onward the title
 of the graphs is displayed in the animation (providing the time for e
ach frame)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 255 "Exercise1: Change the intial conditions by a small amount and \+
compare time histories. Change the method of integration (remove metho
d=lsode to use a Runge-Kutta method), and observe whether the results \+
are being reproduced for the same initial conditions." }}{PARA 0 "" 0 
"" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 102 "Exercise2: Extend the
 problem by allowing different masses for the top and bottom pendula (
pendulums)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 13 }{VIEWOPTS 1 1 0 1 1 1803 1 
1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }