THE SPHERICAL PENDULUM WITH NON-ZERO ANGULAR MOMENTUM 2
K.6.3.3.b
Franz Michel, Bernhard Schnizer

Institut fuer Theoretische Physik
Technische Universitaet Graz
A-8010  Graz, Austria
Projekt Multimedial Lehre des
Zentralen Informatikdienstes der TU Graz

2002-02-11

Introduction

Parameters

Internal Parameters

bilder		numbers of pictures per second of calculation time

bilder = 3 Round[tt];

start		Button to start animations

start = 
        ButtonBox[
          StyleBox["Run Animation", "Title", 
            FontFamily -> "Helvetica"], 
          ButtonFunction :> 
            FrontEndExecute[{SelectionMove[EvaluationNotebook[], Before, 
                    ButtonCell]
SelectionMove[
                    EvaluationNotebook[], Previous, 
                    CellGroup]
SelectionAnimate[
                    SelectedNotebook[], AnimationDisplayTime -> zeit]}], 
          ButtonEvaluator -> Automatic, 
          Active -> True];

Needs["Graphics`Colors`"];

lt1  = 
        Thickness[  .001] ;                       (*  thickness of lines  *)

First Animation

Calculations

pg		polynomial in the equation of motion for the spherical pendulum :  \(TraditionalForm\`u\&● = \ \(\@\(\(2\)\(\ \)\)\) \(ω\_0\) \@\(pg(u)\)\)
pu		decomposition of pg. u1, u2 and u3 are the roots of pg
sod		argument u for the maximum and the minimum of pg

Clear[u1, u2, u3,  ϵ,  pu, 
        pg,  λ, λ1, λ2];
pu  =  (u  -  u1) (u  -  u2) (u  -  
              u3);
pg  =  
        u^3  -  ϵ u^2  -  
          u  +  ϵ  - λ^2/2;

dpg  =  D[pg, u];

eqd  =  dpg  ==  0;
sod  =  Solve[eqd, u]
Remove[eqd];

The second derivative of pg  is calculated to find which value of sod gives the maximum; and which is the minimum

D[dpg, u] /. sod // Together    
Remove[sod];

u3 is expressed as a function of u1 and u2:

eq3  = 
        Coefficient[pu, u]  == 
          Coefficient[pg, u];
so3  =  Solve[eq3, 
          u3]//Flatten 
Remove[
        eq3];

Energy ϵ is given as a function of  u1, u2  and  u3:

eq2 = 
        Coefficient[pu, u^2]  == 
          Coefficient[pg, u^2];
so2  =  Solve[eq2, ϵ] // Flatten
Remove[eq2];

Angular momentum  λ  is given as a function of ϵ, λ, u1, u2 and u3 :

eq0  =  (pu  /.  
              u  ->  
                0)  ==  (pg  /.  {λ  ->  
                  Sqrt[ll], u  ->  0});
so0  =  Solve[eq0, ll] // Flatten;
λ1  =  Sqrt[
          ll]/.so0 
Remove[eq0];

\(k2\  = \ \(k\^2\  = \ \ \(square\ of\ modulus\  = \ 
          parameter\ of\ elliptic\ integrals . \ \ su\ is\ an\ \(\(abbreviation\)\(.\)\)\)\)\)

k2  =  (u2 - u1)/(u3 - 
              u1);                 
su  =  Sqrt[(u3  -  u1)/2];

tt			period
ψ[t_]			parameter of the trajectory as a function of time  t;
up[s_]			latitude as a function of the parameter  ψ
ph[s_]			azimuth as a function of the parameter  ψ
Δph			shift of azimuth after one period

tt =  4  EllipticK[k2]/su;
ψ[t_] = 
        JacobiAmplitude[su t, 
          k2];                              
up[ψ_]  =  
        u1  +  (u2  -  
                u1) Sin[ψ]^2;
ph[ψ_]  = λ1/2/su ( 
            1/(1  -  
                      u1) EllipticPi[(u2 - u1)/(1  -  
                        u1), ψ, 
                  k2] + 
1/(1  +  
                      u1) EllipticPi[-(u2 - u1)/(1  +  
                        u1), ψ, 
                  k2]);         
Δph = ph[2  π] - ph[0] - 2  π;

drehm			matrix of rotation
Since the evaluation of Elliptic Functions in Mathematica takes a long time it is faster to calculate just one period; 
thereafter  the matrix drehm is used to find the coordinates of further trajectory points..

drehm = {{Cos[Δph], -Sin[Δph], 
            0}, {Sin[Δph], Cos[Δph], 0}, {0, 0, 
            1}};

numerical values are assigned to the variables

u1  =   u1N ;   
u2  =   u2N  ;
u3  =  
        u3  /.  
          so3;                            
ϵ  =  ϵ  /. so2;
drehmN = drehm // N ;      
phN[t_] = ph[t] // N;
upN[t_] = up[t] // N;
ΔphN = Δph // N;

cz1		lower latitude circle
cz2.		upper latitude cirle
periode1	coordinates of the pendulum bob as a function of time
xyzRange	PlotRange

Clear[farbe]
SetOptions[ParametricPlot3D, 
        DisplayFunction -> Identity];
cz1  =  
        ParametricPlot3D[{Sqrt[1  -  u1^2] Cos[v], 
            Sqrt[1  -  u1^2] Sin[v],  u1, lt1}, {v, 0, 
            2  Pi}];
cz2  =  
        ParametricPlot3D[{Sqrt[1  -  u2^2] Cos[v], 
            Sqrt[1  -  u2^2] Sin[v],  u2,  lt1}, {v, 0, 
            2  Pi}];
periode1  =  
ParametricPlot3D[{Sqrt[
                1  -  upN[ψ[t]]^2] Cos[phN[ψ[t]]], 
            Sqrt[1  -  upN[ψ[t]]^2] Sin[phN[ψ[t]]], 
            upN[ψ[t]], Hue[farbe]}, {t, 0, tt},  
          AspectRatio  ->  Automatic];

xyzRange = 
        If[u2 > 0, {{-1, 1}, {-1, 1}, {u1, 
              u2}}, {{-Sqrt[1 - u2^2], 
              Sqrt[1 - u2^2]}, {-Sqrt[1 - u2^2], Sqrt[1 - u2^2]}, {u1, 
              u2}}];

Animation

the curve of motions turnes Δθ after each period

periodeviel = 
        Table[Show[
            Join[{cz1, cz2}, 
              Table[
periode1 /. {Line[
                        l : {_, __}] :> ({Line[{MatrixPower[drehmN, 
                                      i] . #1[[1]], 
                                  MatrixPower[drehmN, 
                                      i] . #1[[2]]}]} &) /@ 
                        Partition[l, 2, 1], 
                    farbe ->  (1  - i)/
                        perzahl}, 
{i, 0, j}]], 
            PlotRange -> xyzRange, 
            DisplayFunction ->  $DisplayFunction, 
            ImageSize -> 500, Background ->  RGBColor[1, 1, 1], 
            Axes  ->  True], 
{j, 0, 
            perzahl}];
SelectionMove[EvaluationNotebook[], After, 
        EvaluationCell];
SelectionMove[EvaluationNotebook[], Next, 
        CellGroup];
FrontEndTokenExecute["OpenCloseGroup"];
start /. zeit ->  0.5 // DisplayForm

A Periodic Spherical Pendulum

Input Parameters

Calculations

u3, ϵ and λ are now a function of  u1

Clear[u1, u3, ϵ, λ1]  
u2 = u22N;
u3  =  
        u3  /.  
          so3 ;                     
ϵ  =  ϵ  /. so2;
λ1  =  Sqrt[ll] /. so0;
tt;

The ratio   2π/Δph  is plotted as a function of u1. If 2  π/Δph  is an integer the motion of the pendlum is periodic.

subperU1 = 
        Plot[(2  π)/Δph , {u1, -0.999999999, u2}, 
          PlotLabel ->  "Number of oscillations per period as a \
function of lower latitude cosine u1.\n\n", 
          AxesLabel -> {"cos(θ1)", "Δ\
φ"}, ImageSize ->  600];

Find u1 according to subper :

u1 = 
      Re[u1 /. FindRoot[Δph == 
                2 π/ subper, {u1, -0.9999}] // Chop]

Both values must be identical:

Δph
2 π/ subper // N

If not, you get an error message:

If[Or[u1 < -1,u1 > u2,Round[1000 Δph] ≠ 
              Round[1000 2 π/ subper]],
NotebookPut[Notebook[{Cell["
With the two values of  subper and u1 is no motion of the pendulum possible. \
Chose a new value for u1 or subper . The graph in the subsection " 
                  Calculations " is very helpful to find a proper subper \
for your u1.
", "Text"]}, 
Background -> Hue[0], 
            WindowSize -> {Fit, Fit}, 
            WindowTitle -> "Errormessage"]]

        Abort]

The functions below are  evaluated numerically to speed up calculations:

drehmN = drehm  // N ;
phN[t_] = ph[t] // N;
upN[t_] = up[t] // N;

The trajectory is computed for one period:

periode1 =  
ParametricPlot3D[
{Sqrt[
                1  -  upN[ψ[t]]^2] Cos[phN[ψ[t]]], 
            Sqrt[1  -  upN[ψ[t]]^2] Sin[phN[ψ[t]]], 
            upN[ψ[t]], Hue[farbe]}, 
{t, 0, tt},  
          AspectRatio  ->  Automatic, 
          PlotRange ->  {{-1, 1}, {-1, 1}, {-1, 
                1}}];

cz1		lower latitude circle
cz2		upper latidute cirle
fix		stick where the pendulum is fixed
periodeviel	howl curve of motion calculated with the rotation-matrix drehmN

cz1  =  
        ParametricPlot3D[{Sqrt[1  -  u1^2] Cos[v], 
            Sqrt[1  -  u1^2] Sin[v],  u1, lt1}, {v, 0, 
            2  Pi}];
cz2  =  
        ParametricPlot3D[{Sqrt[1  -  u2^2] Cos[v], 
            Sqrt[1  -  u2^2] Sin[v],  u2,  lt1}, {v, 0, 
            2  Pi}];
fix = 
        Graphics3D[{Thickness[0.01], 
            Line[{{0, 0, 0}, {0, 0, 1}}]}];
perzahl = subper - 1
periodeviel = 
        Table[
periode1 /. 
            Line[l : {_, __}] :> ({Hue[ .2 (i  -  1)], 
                      Line[{MatrixPower[drehmN, i] . #1[[1]], 
                          MatrixPower[drehmN, i] . #1[[2]]}]} &) /@ 
                Partition[l, 2, 1], 
{i, 0, 
            perzahl}];
Show[periodeviel, cz1,  cz2, 
        PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, 
        DisplayFunction ->  $DisplayFunction, ImageSize -> 400, 
        Background ->  RGBColor[1, 1, 1], 
        ViewPoint -> {0.893,  -1.927,  2.634}];

linie = string of the pendulum

linie  =  
        Table[
Graphics3D[{PointSize[
                0.02], 
Line[{{0, 0, 
                    0}, 
{Sqrt[
                        1  -  upN[ψ[t/bilder tt]]^2] Cos[
                        phN[ψ[t/bilder tt]]], 
                    Sqrt[1  -  upN[ψ[t/bilder tt]]^2] Sin[
                        phN[ψ[t/bilder tt]]], 
                    upN[ψ[
                        t/bilder tt]]}
}], 
Point[{Sqrt[1  -  upN[ψ[t/bilder tt]]^2] Cos[
                      phN[ψ[t/bilder tt]]], 
                  Sqrt[1  -  upN[ψ[t/bilder tt]]^2] Sin[
                      phN[ψ[t/bilder tt]]], 
                  upN[ψ[
                      t/bilder tt]]}]}, 
AspectRatio  ->  Automatic, 
            PlotRange ->  {{-1, 1}, {-1, 1}, {-1, 
                  1}}], 
{t, 0, 
            bilder subper}];

Simulation

Table[
        Show[{periodeviel, fix, cz1, cz2, linie[[i]]}, 
          ImageSize ->  600, 
          ViewPoint -> {1.300,  -2.400,  2.000}, Boxed ->  False, 
          Axes -> False, Background ->  RGBColor[1, 1, 1], 
          DisplayFunction ->  $DisplayFunction], {i, 1, 
          bilder subper}];
SelectionMove[EvaluationNotebook[], After, 
        EvaluationCell];
SelectionMove[EvaluationNotebook[], Next, 
        CellGroup];
FrontEndTokenExecute["OpenCloseGroup"];
start /.  zeit ->  ditime // DisplayForm