Elastic Pendulum

Physics 499J project
May 5, 1992
Veikko Saar

Student Project for Mathematica for Scientists and Engineers Course
Professor William M. MacDonald
University of Maryland 
College Park MD 20742

Equations of motion

As the motion of any dynamical system can be described by Hamilton's or 
 Lagrange's equations, I define functions which give these equations if kinetic
 and potential energy is known.

Lagrange's equations

Function Lagrange[kinetic energy,potential energy,coordinate,coordinate,...] gives a list
of Lagrange's "expressions" i.e. {expression,expression,...} where expression==0 is
a Lagrange's equation.

Lagrange[T_,U_,x__]:=
                     Module[{l=T-U},
                           Table[Lagreq[l,{x}[[i]]],
                                {i,Length[{x}]}
                           ]
                     ]
Lagreq[L_,x_]:=Simplify[D[D[L,x'[t]],t]-D[L,x[t]]]

Hamilton's equations

Function Hamilton[kinetic energy,potential energy,coordinate,coordinate,...] gives a list
of Hamilton's "eqpressions" i.e. {{exprp,exprq},{exprp,exprq},...} where p'[t]=exprp
and q'[t]=exprq are Hamilton's equations.
   It defines new variables for momenta: for coordinate x[t] the momentum is px[t].

Hamilton[T_,U_,x__]:=
        Module[{moml,h},
              Off[General::spell1];
              moml=Table[ToExpression[
              		                "p"<>ToString[{x}[[i]]]][t],
               	        {i,Length[{x}]}
                   ];
              On[General::spell1];
              h=Sum[moml[[i]] {x}[[i]]'[t],
                   {i,Length[{x}]}]-(T-U)/.
       		  Flatten[
	     	   		 Solve[
					      Table[
					           moml[[i]]==D[T-U,
					           		  	   {x}[[i]]'[t]
					           		  	  ],
			 	   		       {i,Length[{x}]}
			 	          ],
				          Table[{x}[[i]]'[t],
			 	               {i,Length[{x}]}
			 	          ]
	     	         ]
		      ];
			  Table[
		           Flatten[
	   	                  Hameq[h,{x}[[i]],moml[[i]]]
				   ],
				   {i,Length[{x}]}
		 	  ]
		]
Hameq[H_,x_,momx_]:=Simplify[{D[-H,x[t]],D[H,momx]}]

Elastic pendulum:  simplifications

Kinetic and potential energy

The only equations one must define for Mma by himself are those of energies

T=m/2 (q'[t]^2+((r0+q[t]) theta'[t])^2)
U=k q[t]^2/2+m g (r0-(r0+q[t]) Cos[theta[t]])

Here theta is the angle and q the stretch measured from the equilibrium position of the
spring; r0 is the free length of the spring.

Small angle and stretch approximation

For small oscillations, an approximate analytic solution can be found. Thus it doesn't
matter whether one uses Hamilton's or Lagrange's equations.

leq=Lagrange[T,U,q,theta]

These equations are quite complicated, but for small theta and q, one can take 
sin(theta)=theta ,  cos(theta)=1. In this case the velocities theta' and q' are also small. 
Thus,  terms as theta * q and theta' * q' are 1 order smaller than const * q etc and can be 
neglected.

leq1=Expand[Normal[Series[leq,
                         {theta[t],0,1},
                         {theta'[t],0,1}]]]/.
     {theta[t]q[t]->0,theta'[t]q[t]->0,theta'[t]q'[t]->0,
      theta''[t]q[t]->0,q[t]^2->0}

To simplify the equations, I define w0^2=k/m and w1^2=g/r0.

leq2={Expand[leq1[[1]]/m]/.{k/m->w0^2},
	  Expand[leq1[[2]]/(m r0^2)]/.{g/r0->w1^2}}

The solution is

dsol=Flatten[
            DSolve[{leq2[[1]]==0,leq2[[2]]==0},
                   {theta[t],q[t]},t]
     ]

As only the real part of the solution makes sense, I have to get rid of the imaginary one.
I need the following packages.

<<Algebra`ReIm`
<<Algebra`Trigonometry`

Here define the imaginary part of my variables to be zero.

w0/: Im[w0]=0;
w1/: Im[w1]=0;
t/: Im[t]=0;
g/: Im[g]=0;
Table[c/: Im[c[i]]=0,{i,8}];
nsol=dsol/.Table[C[i]->c[2i-1]+I c[2i],{i,4}]

After simplifying and taking the real parts, I get

solq[t_]=Re[ComplexToTrig[q[t]/.nsol]]
soltheta[t_]=Re[ComplexToTrig[theta[t]/.nsol]]

This shows that both coordinates oscillate independently. If w0=integer*w1 then
the movement forms a cycle. For a cycle k=integer^2 mg/r0. This approximation for k will
be used for finding cycles of the real solution.

Reduced coordinates

From now on the solutions will be numerical. To deal with it better, I'll define reduced
cooordinates: the unit mass'll be m kg, unit length r0 meters but I leave unit 
acceleration = 1 m/s^2. Thus the unit time must be sqrt(r0 in meters) seconds,
k has units m/r0. and energy is measured in units m*r0.

unitsub={m->1,r0->1,g->9.8}

The mathematical pendulum

Taking q=constant, I get the equations of the mathematical pendulum

leq3=Simplify[leq]/.{q[t]->0,q'[t]->0,q''[t]->0}

I need to use only one of them

thsubs=Flatten[Solve[leq3[[2]]==0,theta''[t],t]]

This eq can be solved only numerically.  The function that does it is

findsol[angle_]:=
       Flatten[
              NDSolve[{theta''[t]==(theta''[t]/.
                      thsubs/.unitsub),
                      theta'[0]==0,theta[0]==N[angle]},
                      theta[t],{t,0,5}
              ]
       ]

and the solution for theta[0]=90 deg is

val1[t_]=theta[t]/.findsol[Pi/2]

Plot[val1[t],{t,0,5}]

The period of oscillation is

persub=FindMinimum[-val1[t],{t,2}]

matper1=t/.persub[[2]]

And it is measured in units sqrt(r0). In SI, T=2.369*sqrt(r0 [m]) seconds. The period
depends on energy (or on the initial angle). For theta[0]=45 deg, the solution is

val2[t_]=theta[t]/.findsol[Pi/4];
matper2=t/.FindMinimum[-val2[t],{t,2}][[2]]

thus T is smaller for smaller energies.

Chaos and cycles

Equations for finding cycles

The motion of the elastic pendulum is generally chaotic, but in the case of certan initial
conditions the motion forms a cycle. This cycle must be symmetric around the theta 
axis, thus when theta=0 , I require q''=q'=0 too. Since the solution of Lagrange's 
equations requires less memory, I'll define functions for drawing plots both for L and H 
equations.

Solution of Lagrange's equations

To integrate Lagrange's equations, I need to know q and theta' at the point theta=0.
The total energy at the point theta=0 is

Edown=U+T/.{q'[t]->0,theta[t]->0,q[t]->q1,theta'[t]^2->vth12}

where vth12=theta'^2 can be found from the  requirement q''=0

vsol=Flatten[
     Solve[m g+m vth12 ((q1+r0))^2/(q1+r0)-k q1==0,vth12]]

Etot is the total energy of the system. Now I can find q[0] and theta[0]

q1subs=Factor[Solve[Etot==Edown/.vsol,q1]][[1]]/.unitsub;
q1sol=q1/.q1subs

vth1subs=vsol/.q1subs/.unitsub;
vth1sol=vth12/.vth1subs

The solution for q[0] and theta[0] depends only on k and E. I need  a function for finding 
and drawing the plot (dependence of q on theta).

Off[ParametricPlot::ppcom]
Ldraw[E_,K_,T_]:=Module[{su={k->K,Etot->E}},
    ans=Flatten[
               NDSolve[
                 {leq[[1]]==0,leq[[2]]==0,
                 theta[0]==0,theta'[0]==
                 Sqrt[vth1sol],q'[0]==0,
                 q[0]==q1sol}/.su/.unitsub,
                 {theta[t],q[t]},{t,0,T}]];
    ParametricPlot[{theta[t],q[t]}/.ans,{t,0,T 9/10}];
]

Solution of Hamilton's equations

Hamilton's expressions are

heq=Flatten[Hamilton[T,U,q,theta]/.unitsub]

I have to know ptheta[0],pq[0],q[0] and theta[0]=0 in order to solve them.

pth0=D[T-U,theta'[t]]/.
     {q[t]->q1sol,theta'[t]->Sqrt[vth1sol]}/.unitsub
pq0=D[T-U,q'[t]]/.{q'[t]->0}/.unitsub

The function that solves the differential equation is

psolve[su_,T_]:=
      Flatten[
             NDSolve[
                 {heq[[1]]==pq'[t],heq[[2]]==q'[t],
                  heq[[3]]==ptheta'[t],heq[[4]]==theta'[t],
                  theta[0]==0,q[0]==q1sol,
                  ptheta[0]==pth0,
                  pq[0]==0}/.su/.unitsub,
                 {theta[t],q[t],ptheta[t],pq[t]},{t,0,T}]
      ]

Hdraw takes in to parameters one has to plot, for example var1=pq and var2=q.

Hdraw[E_,K_,T_,var1_,var2_]:=Module[{su={k->K,Etot->E}},
     ans=psolve[su,T];
     ParametricPlot[{var2[t],var1[t]}/.ans,{t,0,T 9/10}];
]

I define also a function to draw together plots of q-theta, pq-q and ptheta-theta.

Tdraw[E_,K_,T_]:=Module[{su={k->K,Etot->E}},
     ans=psolve[su,T];
     Show[GraphicsArray[{
        ParametricPlot[{theta[t],q[t]}/.ans,{t,0,T 9/10},
                      DisplayFunction->Identity],
        ParametricPlot[{pq[t],q[t]}/.ans,{t,0,T 9/10},
                      DisplayFunction->Identity],
        ParametricPlot[{ptheta[t],theta[t]}/.ans,{t,0,T 9/10},
                      DisplayFunction->Identity]
        }]
    ]; 
]

N-fold cycles

Cycles usually look like this:

where the number of crossings is 2,4,6,...  but between each two of such cycles (for a 
constant energy, depening on k) there is a back-forth line/cycle:

The q-theta,  pq-q and ptheta-theta plots for this type of cycles are

Tdraw[5,504,10]

Tdraw[1,450,10]

The shape of the cycle depends on the number of crossings N (on the q-theta plot), and on 
the total energy. Using one of the functions
        Ldraw[energy,k,time] for q-theta plot
        Hdraw[energy,k,time,q1,q2] for q1-q2 plot
        Tdraw[energy,k,time] for q-theta, pq-q and ptheta-theta plots,
one can find cycles for different values of E. For example

Hdraw[7,147,5,q,theta]

draws a q-theta graph for which the total energy is 7 (reduced units), k=147 and the
pendulum oscillates for 5 time units. The only requirement is that one must choose
E>=-48.02/k in order not to make the kinetic energy negative. It may seem that the
motion often looks as simple as the graphs above, but if one chooses arbitrary values, 
he usually gets something like this

Tdraw[5,80,10]

Thus one can find the allowed q-theta region that is determined by the principle 
Etot=const.

I calculated numerically the dependence of k=f(N,E) using Ldraw. The result was

tab=
{{61.5,138,220,326,450,595,755,938,1140,1360,1590,1850}, 
 {45.2,123.5,195.5,296,414,550,703,875,1064,1270,1500,1740}, 
 {34.4,110.1,171.4,265.8,376.5,504,648.5,809.5,987.8,1176, 1400, 1630},
 {27.5,97,147.5,236,338,458,592,741,910,1088,1290,1500}, 
 {23,84,124,205,300,410,534,670,825,990,1173,1372}, 
 {19.9,70,100.5,175,260,359,472,597,737,885,1055,1235},
 {17.64,57,78,144,218,305,406,517,641,780,925,1085}};
TableForm[Transpose[tab],TableHeadings->
{{2,3,4,5,6,7,8,9,10,11,12,13},{1,3,5,7,9,11,13}},
TableSpacing->{1}]

where x-axis is the energy and y- the N axis. The values are those of k. Rows, where
N is odd, are those of the back-forth oscillations. The graph k=f(N,E) looks as

ListPlot3D[tab];

The function k=f(E,N=2) is

clien2={61.5,52.5,45,39.3,34.4,30.6,27.5,25,23,21.3,19.9,
        18.7,17.64,16.72};
ListPlot[clien2,PlotStyle->PointSize[0.02],
        AxesLabel->{E,k}];

but its fittings do not give good results.

List of cycles' types

There are several more cycles in addition to the N-fold ones. The following cycles form
a family, just like the N-fold one. One can see that easily, if he looks at the ptheta-theta
plots.  Between each two loop-cycles there is a line-cycle, just as in the case of the 
N-fold cycles. The number of these cycles should therefore be infinity.

SetOptions[NDSolve,MaxSteps->1000];

Tdraw[7,17.9,8]