Time dependent
solutions of the pendulum equations

by Robert C. Pickett, Gordon E. Lindgren, and Ronald K Anderson, Department of Physics, Bemidji State University, Bemidji, MN 56601.  Please direct questions and comments to this address or email them to rpickett@vax1.bemidji.msus.edu.  Publications using results from these programs must contain a reference to the authors' papers below.  This note must not be removed.

Theory

This  notebook  discusses methods of solving the simple pendulum equation.    It  is  the first in a series of notebooks showing how Mathematica can be used  to find exact power series solutions (to arbitrary order) of several differential equations.

A point mass at the end of a massless rigid rod is constrained to move without friction on a vertical circle about a fixed point.  If the polar angle (phi) is measured counter-clockwise from vertically downward, the dimensionless form of the simple pendulum equation for the polar angle  (phi) is:

                         	phi'' +  sin(phi)  =  0  where  phi(0) =  phi00, and  phi'(0) = phi10.
	
To reduce the pendulum equation to a finite polynomial format, the new variables  omega =  phi', x = sin(phi),  and  y = - cos(phi) are introduced.   This  yields the system of first order equations :

                                          phi' = omega  with  phi(0) =  phi00,
                                         	omega' = - x with  omega(0) =  phi10 = omega00,
                                         	x'  =  - omega y  with  x(0) = x00 =  sin(phi00), and
                                         	y'  = omega x  with  y(0) = y00  =  -  cos(phi00).

To solve the system of first order differential equations, we define the following power series :

                                          phi = Sum[Phi[n] t^n,{n,0,Infinity}],
                                          omega = Sum[Omega[n] t^n,{n,0,Infinity}],
	                                         x = Sum[X[n] t^n,{n,0,Infinity}],
                                         	y = Sum[Y[n] t^n,{n,0,Infinty}],
                                         	omega y = Sum[OmegaY[n] t^n,{n,0,Infinity}], and
                                         	omega x = Sum[OmegaX[n] t^n,{n,0,Infinity}].

One can extract the recursion equations and the initial values:

                                          Phi[n_]:=Omega[n-1}/n,	 Phi[0]=phi00,
                                         	Omega[n_]:=-X[n-1]/n,	 Omega[0]=omega00,
	                                         X[n_]:=-OmegaY[n-1]/n,	 X[0]=x00,
                                          Y[n_]:=OmegaX[n-1]/n, and Y[0]=y00 where
                                          OmegaY[n_]:=Sum[Omega[m]Y[n-m],{m,0,n}], and
                                          OmegaX[n_]:=Sum[Omega[m]X[n-m],{m,0,n}].

Since  phi00,  omega00,  x00, and y00 are known, all the coefficients can be computed.   If the potential energy is  defined as  y and the  kinetic energy as  omega^2 /2 , then the  total energy  is   e = omega^2 /2 + y .  Of course, we could have chosen to use only one dependent variable, e.g., phi, x, or y.

The following Mathematica   function is used throughout the notebook for both symbolic and numerical solutions.

Clear[PendulumPowerSeries];

PendulumPowerSeries[]:={
Clear[Phi,Omega,X,Y,OmegaY,OmegaX,
	phi00,omega00,x00,y00,e];
Phi[n_]:=Omega[n-1]/n;	Phi[0]=phi00;
Omega[n_]:=Omega[n]=-X[n-1]/n;	Omega[0]=omega00;
X[n_]:=X[n]=-OmegaY[n-1]/n;	X[0]=x00;
Y[n_]:=Y[n]=OmegaX[n-1]/n;	Y[0]=y00;
OmegaY[n_]:=Together[Sum[Omega[m]Y[n-m],{m,0,n}]];
OmegaX[n_]:=Together[Sum[Omega[m]X[n-m],{m,0,n}]];
};

Numerical calculations

The potential energy y of the  pendulum is defined to be +1 at the top, -1 at the bottom andthus 0 at a point level with the support point.  Depending on the value of the total energy e, the pendulum will move in one of three possible modes:  rotation, libration, or stagnation.  The
system rotates if e>1, librates (oscillates) if e<1, and stagnates (ultimately stops at the top) when e=1.  The latter case defines the separatrix which separates the rotation and libration regions.

The initial conditions are defined for motion with the least possible kinetic energy in the rotation and the libration regions .  The pendulum released from rest (omega00=0) at the angle phi=phi00 yields libration.  Rotation results when the pendulum is given an angular velocity omega00=Sqrt[2(e-1)] (with e>1) at the top (phi=Pi).  Finally, stagnation results when the pendulum is given just enough energy to make it to the top (in infinite time).  This will happen if omega00=2 when phi00=0.  Since the pendulum bob moves along a unit circle, the initial values of x and of y are as shown.

The time t1 is defined for the three cases as follows:  If e=1 then the theoretical time t1 to reach the top (period/2) is infinity.  When the pendulum is released from rest, the theoretical time to reach the bottom (period/4) is EllipticK[(e+1)/2].  The theoretical half-period for the rotation case is EllipticK[2/(e+1)]/Sqrt[(e+1)/2].

The following Mathematica   function is used throughout the notebook for numerical solutions.

Clear[NumericalPendulumPowerSeries];

NumericalPendulumPowerSeries[]:={
PendulumPowerSeries[];
omega00=If[e<1,0,If[e>1,Sqrt[2(e-1)],2]]//N;
x00=If[e<1,Sqrt[1-e^2],0]//N;
y00=If[e==1,-1,If[e<1,e,1]]//N;
phi00=If[e==1,0,If[e>1,Pi,ArcSin[x00]]]//N;
t1=If[e==1,Infinity,If[e<1,EllipticK[(e+1)/2],
	EllipticK[2/(e+1)]/Sqrt[(e+1)/2]]]//N;
};

Angular frequency and the order of approximation needed to obtain 0.0001% accuracy

A plot of the angular frequency vs energy is shown.  The solid line is found from the theoretical times while the dots result from the power series solutions.  The number of terms required to achieve 0.0001% accuracy is also shown.  Though as yet unexplained, the fine structure appears to be real.

The following code generates the angular frequency plots

Data tables

lst0={{-0.99,38},{-0.9,50},{-0.8,58},{-0.7,60},
  {-0.6, 58},{-0.5,60},{-0.4,62},{-0.3,64},
  {-0.2,58},{-0.1,60},{0,62},{0.05,72},{0.1,88},
  {0.2,92},{0.3,96},{0.4,98},{0.45,110},
  {0.5,122},{0.6,128},{0.7,124},{0.75,128},
  {0.8,144},{0.85,152},{0.9,162},{1.1,162},
  {1.2,146},{1.3,128},{1.4,122},{1.5,118},
  {1.6, 114},{1.7,112},{1.8,112},{1.9,106},
  {2.,104},{2.1,104},{2.2,102},{2.3,100},
  {2.4,98},{2.5,98},{2.6,96},{2.7,88},{2.8,88}, 
  {2.9,86},{3.,86}};

lst1={{-0.99, 0.998748041485949},{-0.9,0.987299130653834},
  {-0.8,0.974173519145861},{-0.7,0.96058179576567},
  {-0.6,0.946478398387498},{-0.5,0.931809286065574},
  {-0.4,0.916510570535938},{-0.3,0.900505758446479},
  {-0.2,0.883702216599437},{-0.1,0.865985854627746},
  {0,0.847213878360317},{0.1,0.827203744544505},
  {0.2,0.80571569883907},{0.3,0.7824260614623946},
  {0.4,0.7568782701771136},{0.5,0.7283958731642441},
  {0.6,0.6959034894812896},{0.7,0.6575078436052007},
  {0.8,0.6092865924585951},{0.9,0.5401013608825669},
  {1.1,1.0980013736494},{1.2,1.256156857315814},
  {1.3,1.372974409013424},{1.4,1.470390224803143},
  {1.5,1.556087568392301},{1.6,1.633761611337674},
  {1.7,1.705507036414285},{1.8,1.772640225002212},
  {1.9,1.836049561492536},{2.,1.896367895636579},
  {2.1,1.954061426676085},{2.2, 2.009490157721578},
  {2.3,2.062934632560696},{2.4,2.114620664714678},
  {2.5,2.164732301038067},{2.6,2.213422877734657},
  {2.7,2.260820966671498},{2.8,2.307036114921357},
  {2.9,2.352163006019662},{3.,2.396282110028069}};

Show[Plot[Pi/2/EllipticK[(e+1)/2],{e,-1,1},
	DisplayFunction->Identity],
	ListPlot[lst1,DisplayFunction->Identity],
	Plot[Pi Sqrt[(e+1)/2]/EllipticK[2/(e+1)],
	{e,1,3},DisplayFunction->Identity],
	DisplayFunction->$DisplayFunction,
	Frame->True,Axes->False,FrameLabel->{"e","nu"},
	PlotLabel->"Angular Frequency vs Energy",
	Prolog->AbsolutePointSize[3]];
ListPlot[lst0,Frame->True,Axes->False,FrameLabel->{
	"e","Order"},Prolog->AbsolutePointSize[3],
	PlotLabel->"Order of Approximation vs Energy"];

The following code can be used to calculate the data tables

lst0={};
lst1={};

NumericalPendulumPowerSeries[];
e=3//N;
OrderOf=88;
Do[Omega[n]; X[n]; Y[n],{n,0,OrderOf}];

OrderOf=86;
Do[Omega[n]; X[n]; Y[n],{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];
answer=FindRoot[y+1,{t,t1}];

Abs[((t/.answer)/t1-1)]

Abs[((t/.answer)/t1-1)]<10^-6

{e,OrderOf}

{e,If[e<1,1/2,1]N[Pi]/t/.answer}

AppendTo[lst0,{e,OrderOf}];
AppendTo[lst1,{e,If[e<1,1/2,1]N[Pi]/t/.answer}];

Libration (-1 < e < 1)

The following code generates the libration region data, with e=0

NumericalPendulumPowerSeries[];
e=0;	(*e=0*)
OrderOf=62;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

Time plots

Plot[phi,{t,-t1,t1},Frame->True,FrameLabel->
	{"t","phi"},PlotLabel->"phi vs t"];
Plot[omega,{t,-t1,t1},Frame->True,FrameLabel->
	{"t","omega"},PlotLabel->"omega vs t"];
Plot[x,{t,-t1,t1},PlotRange->All,Frame->True,
	FrameLabel->{"t","x"},PlotLabel->"x vs t"];
Plot[y,{t,-t1,t1},Frame->True,FrameLabel->{"t","y"},
	PlotLabel->"y vs t"];

Phase space plots

ParametricPlot[{{phi,omega},{-phi,-omega}},{t,-t1,t1},
	Frame->True,FrameLabel->{"phi","omega"},
	PlotLabel->"omega vs phi"];
Vx=Sum[X[n+1](n+1)t^n,{n,0,OrderOf}];
ParametricPlot[{{x,Vx},{-x,-Vx}},{t,-t1,t1},Frame->True,
	FrameLabel->{"x","Vx"},PlotLabel->"Vx vs x"];
Vy=Sum[Y[n+1](n+1)t^n,{n,0,OrderOf}];
ParametricPlot[{y,Vy},{t,-t1,t1},Frame->True,
	FrameLabel->{"y","Vy"},PlotLabel->"Vy vs y"];

Separatrix ( e == 1 )

At the separatrix, the time needed for the pendulum to reach unstable equilibrium (stopping at the top) is infinite.  The exact series solution for a dependent variable, e.g., phi, will theoretically contain an infinite number of terms.  In practice, of course, the series is truncated.  The worst case absolute truncation error is bounded by the absolute value of the first non-zero term beyond the truncation.  The plots below each consists of a theoretical curve and power series solution of order 400.  In each case, the power series solution deviates sharply from theory at t2+.1, about 1.59 sec.  Remember, the bob started at the bottom (phi=0) at time t=0 with just enough energy (e=1) to reach the top in infinite time.  The power series solution in y does not deviate from theory by more than 0.0001%, until about 1.49 sec.  The series solution is plotted a small time beyond time t2 (0.1 sec for phi).  The additional time is chosen just large enough to display a detectable deviation without disturbing the scale of the plot.

The following code generates the separatrix data

NumericalPendulumPowerSeries[];
e=1;
OrderOf=400;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf+2}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];
next=If[EvenQ[OrderOf],OrderOf+2,OrderOf+1];
t2=t/.FindRoot[Y[next]t^next/(
	1-2Sech[t]^2)==10^-6,{t,1.49176}];

Time plots

Show[Plot[phi,{t,-(t2+.1),(t2+.1)},
	DisplayFunction->Identity],
	Plot[2ArcTan[Sinh[t]],{t,-5,5},DisplayFunction->
	Identity],DisplayFunction->$DisplayFunction,
	Frame->True,FrameLabel->{"t","phi"},PlotLabel->
	"phi vs t"];
Show[Plot[omega,{t,-(t2+.074),(t2+.074)},
	DisplayFunction->Identity],
	Plot[2Sech[t],{t,-5,5},DisplayFunction->Identity],
	DisplayFunction->$DisplayFunction,Frame->True,
	FrameLabel->{"t","omega"},PlotLabel->"omega vs t"];
Show[Plot[x,{t,-(t2+.052),(t2+.052)},
	DisplayFunction->Identity],
	Plot[2Tanh[t]Sech[t],{t,-5,5},DisplayFunction->
	Identity],DisplayFunction->$DisplayFunction,
	Frame->True,FrameLabel->{"t","x"},PlotLabel->
	"x vs t"];
Show[Plot[y,{t,-(t2+.052),(t2+.052)},
	DisplayFunction->Identity],
	Plot[1-2Sech[t]^2,{t,-5,5},DisplayFunction->
	Identity],DisplayFunction->$DisplayFunction,
	Frame->True,FrameLabel->{"t","y"},
	PlotLabel->"y vs t"];

Phase space plots

The separatrix ( e = 1) phase space plots below demonstrate that the power series solutions are very accurate for small values of time but differ from theory at some point .  Again, this point is at a time of about 1.49 sec, when phi is somewhat more than 2.25 radians.

Show[ParametricPlot[{phi,omega},{t,-(t2+.074),(t2+.074)},
	DisplayFunction->Identity],ParametricPlot[
	{2ArcTan[Sinh[t]],2Sech[t]},{t,-5,5},
	DisplayFunction->Identity],Frame->True,
	FrameLabel->{"phi","omega"},PlotLabel->
	"omega vs phi",DisplayFunction->$DisplayFunction];
Vx=Sum[X[n+1](n+1)t^n,{n,0,OrderOf}];
Show[ParametricPlot[{x,Vx},{t,-(t2+.03),(t2+.03)},
	DisplayFunction->Identity],ParametricPlot[
	{2Tanh[t]Sech[t],(3-Cosh[2t])Sech[t]^3},{t,-5,5},
	DisplayFunction->Identity],Frame->True,
	FrameLabel->{"x","Vx"},PlotLabel->"Vx vs x",
	DisplayFunction->$DisplayFunction];
Vy=Sum[Y[n+1](n+1)t^n,{n,0,OrderOf}];
Show[ParametricPlot[{y,Vy},{t,-(t2+.031),(t2+.031)},
	DisplayFunction->Identity],
	ParametricPlot[{1-2Sech[t]^2,4Sech[t]^2Tanh[t]},
	{t,-5,5},DisplayFunction->Identity],
	Frame->True,FrameLabel->{"y","Vy"},PlotLabel->
	"Vy vs y",DisplayFunction->$DisplayFunction];

Rotation (1 < e < Infinity)

When e>1, the pendulum will never come to rest.  In the plots below, omega, x, and y are cyclic in time whereas phi continues to increase as time increases.

The following code generates the rotation region data, with e=2

NumericalPendulumPowerSeries[];
e=2; (* e=2 *)
OrderOf=104;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

Time plots

Plot[phi,{t,-t1,t1},Frame->True,FrameLabel->
	{"t","phi"},PlotLabel->"phi vs t"];
Plot[omega,{t,-t1,t1},Frame->True,FrameLabel->
	{"t","omega"},PlotLabel->"omega vs t"];
Plot[x,{t,-t1,t1},PlotRange->All,Frame->True,
	FrameLabel->{"t","x"},PlotLabel->"x vs t"];
Plot[y,{t,-t1,t1},Frame->True,FrameLabel->{"t","y"},
	PlotLabel->"y vs t"];

Phase space plots

ParametricPlot[{phi,omega},{t,-t1,t1},
	Frame->True,FrameLabel->{"phi","omega"},
	PlotLabel->"omega vs phi"];
Vx=Sum[X[n+1](n+1)t^n,{n,0,OrderOf}];
	ParametricPlot[{x,Vx},{t,-t1,t1},Frame->True,
	FrameLabel->{"x","Vx"},PlotLabel->"Vx vs x"];
Vy=Sum[Y[n+1](n+1)t^n,{n,0,OrderOf}];
	ParametricPlot[{y,Vy},{t,-t1,t1},Frame->True,
	FrameLabel->{"y","Vy"},PlotLabel->"Vy vs y"];

Symbolic calculations

To this point, all solutions have been evaluated numerically.   Now, symbolic programming is used to demonstrate that, to the order of the series, the power series solutions are exact.

Constants of the  motion

The constants of the motion are :	x^2 + y^2 == 1, and
                                        e == omega^2 /2 + y.
The power series solutions of the two equations are demonstrated to be exact to the order of the approximation.

OrderOf=10;
PendulumPowerSeries[];
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

Simplify[Normal[Series[omega^2/2+y,{t,0,OrderOf}]]]==
	omega00^2/2+y00
Simplify[Normal[Series[Together[x^2+y^2],
	{t,0,OrderOf}]]/.{x00->Sin[phi00],y00->
		-Cos[phi00]}]==
	1

Differential equations

The differential equations are :	phi'' + Sin[phi] == 0,	
	                                    omega''+omega^3 /2-e  omega == 0,
	                                    x'' +2 e x - 3 Sqrt[1-x^2] x == 0, and
	                                    y''+2e  y-3 y^2+1 == 0.
The power series solutions of the four differential equations are demonstrated to be exact to the order of the approximation.

OrderOf=6;
PendulumPowerSeries[];
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf+2}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf+2}];
x=Sum[X[n]t^n,{n,0,OrderOf+2}];
y=Sum[Y[n]t^n,{n,0,OrderOf+2}];

Simplify[Together[Normal[Series[Together[
	D[phi,{t,2}]+Sin[phi]],{t,0,OrderOf}]
		]/.{e->omega00^2/2+y00}]/.{
			x00->Sin[phi00],y00->-Cos[phi00]}]==
	0
Simplify[Together[Normal[Series[Together[
	D[omega,{t,2}]+omega^3 /2-e omega],
		{t,0,OrderOf}]/.{e->omega00^2/2+y00}]
			]/.{x00->Sin[phi00],y00->-Cos[phi00]}]==
	0
Simplify[Together[Normal[Series[Together[
	D[x,{t,2}]+ 2 e x -3 Sqrt[1-x^2] x],
		{t,0,OrderOf}]]/.{e->omega00^2/2+y00
			}]/.{y00->Sqrt[1-x00^2]}]==
	0
Simplify[Together[Normal[Series[Together[
	D[y,{t,2}]+2e y-3y^2+1],
		{t,0,OrderOf}]]/.{e->omega00^2/2+y00}
			]/.{x00->Sin[phi00],y00->-Cos[phi00]}]==
	0

Series approximation for phi, omega, x, and y.

OrderOf=2;
PendulumPowerSeries[];
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Simplify[Sum[Phi[n]t^n,{n,0,OrderOf}]/.{
	x00->Sin[phi00],y00->-Cos[phi00]}];
omega=Simplify[Sum[Omega[n]t^n,{n,0,OrderOf}]/.{
	x00->Sin[phi00],y00->-Cos[phi00]}];
x=Simplify[Sum[X[n]t^n,{n,0,OrderOf}]/.{
	x00->Sin[phi00],y00->-Cos[phi00]}];
y=Simplify[Sum[Y[n]t^n,{n,0,OrderOf}]/.{
	x00->Sin[phi00],y00->-Cos[phi00]}];
Print["phi = ",phi];
Print["omega = ",omega];	Print["\n"];
Print["x = ",x];	Print["\n"];
Print["y = ",y];

Libration	( 	- 1  <  e	 <  1 	)

The solutions approach the theoretical solutions as e->-1

The theoretical solutions are : Limit[phi/phi00,e->-1]==Cos[t],
	                                              Limit[omega/phi00,e->- 1]==-Sin[t],
	                                              Limit[x/x00,e ->-1]==Cos[t], and
                                               Limit[(y+1)/(e+1),e->-1]==Cos[t]^2.
The four power series solutions are demonstrated to be exact to the order of the approximation.

PendulumPowerSeries[];
y00=e;
x00=Sqrt[1-y00^2];
phi00=ArcSin[x00];
omega00=0;
OrderOf=10;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

(Simplify[phi/phi00/.ArcSin[Sqrt[1-e^2]]->
	Sqrt[1-e^2]]/.{e->-1})==
	Normal[Series[Cos[t],{t,0,OrderOf}]]
(Simplify[omega/phi00/.ArcSin[Sqrt[1-e^2]]->
	Sqrt[1-e^2]]/.{e->-1})==
	Normal[Series[-Sin[t],{t,0,OrderOf}]]
Limit[x/x00,e->-1]==
	Normal[Series[Cos[t],{t,0,OrderOf}]]
Limit[(y+1)/(e+1),e->-1]==
	Normal[Series[Cos[t]^2,{t,0,OrderOf}]]

The solutions approach the theoretical solutions as e->1

The theoretical solutions are : Limit[phi,e->1]==2ArcTan[Sinh[t]],
                                               Limit[omega,e->1]==2Sech[t],
                                               Limit[x,e->1]==2Tanh[t]Sech[t], and
                                               Limit[y,e->1]==1-2Sech[t]^2 .
The four power series solutions are demonstrated to be exact to the order of the approximation.

PendulumPowerSeries[];
phi00=0;
phi10=Sqrt[2(e+1)];
x00=0;
y00=-1;
OrderOf=10;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

(phi/.e->1)==
	Normal[Series[2ArcTan[Sinh[t]],{t,0,OrderOf}]]
(omega/.e->1)==
	Normal[Series[2Sech[t],{t,0,OrderOf}]]
(x/.e->1)==
	Normal[Series[2Tanh[t]Sech[t],{t,0,OrderOf}]]
(y/.e->1)==
	Normal[Series[1-2Sech[t]^2,{t,0,OrderOf}]]

Separatrix	( 	e = 1	 )

The exact solution are : phi ==2 ArcTan[ Sinh[ t ] ],
	                                    omega == 2 Sech[ t ],
	                                    x == 2 Tanh[ t ] Sech[ t ], and
	                                    y == 1 - 2 Sechh[ t ]^2.
The four power series solutions are demonstrated to be exact to the order of the approximation.

PendulumPowerSeries[];
phi00=0;
phi10=2;
x00=0;
y00=-1;
OrderOf=10;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];

phi=Sum[Phi[n]t^n,{n,0,OrderOf}]==
	Normal[Series[2ArcTan[Sinh[t]],{t,0,OrderOf}]]
omega=Sum[Omega[n]t^n,{n,0,OrderOf}]==
	Normal[Series[2Sech[t],{t,0,OrderOf}]]
x=Sum[X[n]t^n,{n,0,OrderOf}]==
	Normal[Series[2Tanh[t]Sech[t],{t,0,OrderOf}]]
y=Sum[Y[n]t^n,{n,0,OrderOf}]==
	Normal[Series[1-2Sech[t]^2,{t,0,OrderOf}]]

Rotation	(	1 < e < Infinity		)

The solutions approach the theoretical solutions as e->1

The theoretical solutions are : Limit[phi,e->1]==2ArcTan[Sinh[t]],
                                               Limit[omega,e->1]==2Sech[t],
                                               Limit[x,e->1]==2Tanh[t]Sech[t], and
                                               Limit[y,e->1]==1-2Sech[t]^2 .
The four power series solutions are demonstrated to be exact to the order of the approximation.

PendulumPowerSeries[];
phi00=0;
omega00=Sqrt[2(e+1)];
x00=0;
y00=-1;
OrderOf=10;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

(phi/.e->1)==
	Normal[Series[2ArcTan[Sinh[t]],{t,0,OrderOf}]]
(omega/.e->1)==
	Normal[Series[2Sech[t],{t,0,OrderOf}]]
(x/.e->1)==
	Normal[Series[2Tanh[t]Sech[t],{t,0,OrderOf}]]
(y/.e->1)==
	Normal[Series[1-2Sech[t]^2,{t,0,OrderOf}]]

The solutions approach the theoretical solutions as e->Infinity

The theoretical solutions are : Limit[omega/Sqrt[2e],e->Infinity]==1,
                                               Limit[x/Sin[Sqrt[2e]t],e->Infinity]==-1, and
                                               Limit[y/Cos[Sqrt[2e]t],e->Infinity]==1.
The three power series solutions are demonstrated to be exact to the order of the approximation.

PendulumPowerSeries[];
phi00=Pi;
omega00=Sqrt[2(e-1)];
x00=0;
y00=1;
OrderOf=4;
Do[Omega[n];X[n];Y[n],{n,0,OrderOf}];
phi=Sum[Phi[n]t^n,{n,0,OrderOf}];
omega=Sum[Omega[n]t^n,{n,0,OrderOf}];
x=Sum[X[n]t^n,{n,0,OrderOf}];
y=Sum[Y[n]t^n,{n,0,OrderOf}];

Limit[Together[Normal[Series[omega/Sqrt[2(e-y)],
	{t,0,OrderOf}]]],e->Infinity]==1
Limit[Together[x/Normal[Series[Sin[Sqrt[2e]t],
	{t,0,OrderOf}]]],e->Infinity]==-1
Limit[y/Normal[Series[Cos[Sqrt[2e]t],{t,0,OrderOf}]],
	e->Infinity]==1