Action-Angle Variables for Pendulum

First Plots of the phase space orbits

Plot[{-Sqrt[3+Cos[t]],Sqrt[3+Cos[t]]},{t,0,10}]

Plot[{-Sqrt[1/2+Cos[t]],Sqrt[1/2+Cos[t]]},{t,0,10}]

Show[%%,%]

Now consider Pendelum with the Hamiltonian

H = E = (((Subscript[p, θ])^2)/(2m (a^2)))- m g a cosθ

Rescaling the energy, (OverHat[E]) = E/(mga), we have

(Subscript[p, θ]) = ± (Sqrt[2m (a^2)(E+mga cosθ )]) =  ± (Sqrt[2(m^2) (a^3)g ((OverHat[E])+cosθ )])

Let’s leave off the constant prefactor (Sqrt[2(m^2) (a^3)g ]) for simplicity (think of units of angular momentum where this is one).

Calculate J for the rotation (Ee=(OverHat[E]),  t=θ)

Integrate[2Sqrt[Ee+Cos[t]],{t,0,π},Assumptions->{Ee>1}]//PowerExpand

Calculate J for the oscillation (here Ee = -(Subscript[cosθ, 0]),  t0=(Subscript[θ, 0]))

Integrate[4Sqrt[Cos[t]-Cos[t0]],{t,0,t0}]

%[[1]]/. t0->ArcCos[-Ee]//FullSimplify

Jresult[Ee_]:=If[Ee<1,
8Sqrt[2] (Sqrt[1+Ee]) EllipticE[((ArcCos[-Ee])/2),(2/(1+Ee))],
4 Sqrt[2](Sqrt[-1+Ee]) EllipticE[-(2/(-1+Ee))]]

Calculate the derivative and define it as a new function. (The use of Evaluate makes it do the derivative ahead of time.)

DJresult[Eee_]:=If[Eee<1,Evaluate[D[8 Sqrt[2](Sqrt[1+Ee]) EllipticE[((ArcCos[-Ee])/2),(2/(1+Ee))],Ee]/.Ee->Eee],Evaluate[D[4Sqrt[2] (Sqrt[-1+Ee]) EllipticE[-(2/(-1+Ee))],Ee]/.Ee->Eee]]

DJresult[Eee]//Simplify

?EllipticE

?EllipticF

Plot J versus Ee = E / (mga)

Plot[Jresult[Ee],{Ee,-1,5}]

Plot dJ/dEe versus Ee = E / (mga)

Plot[DJresult[Eee],{Eee,-1,3},PlotRange->{0,15}]

The blow up at Ee = +1  is logarithmic and hence integrable (which is why we do not see it in the J versus Ee plot). We can check this by doing a Taylor series about 1:

Series[((2 (Sqrt[2]) EllipticK[-(2/(-1+Eee))])/(Sqrt[-1+Eee])),{Eee,1,2}] /. Arg[-1+Eee]->0

Plot dEe/dJ  versus Ee = E / (mga)
(accounting for the constant that we scaled out earlier, the overall constant for this derivative is exactly ν of the small amplitude solution, (1/(2π))(Sqrt[(g/a)]), which will be shown as 1.0 on the plot)

2 π/DJresult[-0.999999]//Chop

Plot[{1,2 π/DJresult[Eee]},{Eee,-1,3},PlotRange->{0,3}]

Look at the blow up from the right side

Series[1/(((2 (Sqrt[2]) EllipticK[-(2/(-1+Eee))])/(Sqrt[-1+Eee]))),{Eee,1,1}] /. Arg[-1+Eee]->0

Zoom in

Plot[{1,Re[2 π/DJresult[Eee]]},{Eee,0.999,1.001},PlotRange->{0,1}]