SIMULATING THE MATHEMATICAL PENDULUM
			  K.6.3.3.a
			Franz Michel, Bernhard Schnizer
										   Institut fuer Theoretische Physik
								Technische Universitaet Graz
								    A-8010  Graz, Austria
								Projekt Multimediale Lehre des
							Zentralen Informatikdienstes der TU Graz
												2002-01-31

Simulation of the pendulum

Parameters

internal parameters

bilderprosek	numbers of pictures per second of the calculated period
gl		thickness of lines
esp		small number

bilderprosek = 4;   
gl = 0.04;              
esp = 
        0.002;                 
 (*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];

Calculations

tzut must be bigger than 1 and a fraction

If[tzut < 1, tzut = 1/tzut]
tzut = Rationalize[tzut]

T	period for δ01

T = 4 EllipticK[Sin[δ01/2]^2];

calculate the second amplitude dependence on tzut

solution = 
        FindRoot[
          4 EllipticK[Sin[x/2]^2] == tzut T, {x, 1}];

ampitude for the second pendulum

δ02 = x /. solution;
δ02 = Mod[δ02, 2 π]

δ0 = {δ01, δ02}

exact solution of the differential equation of the mathematical pendulum

δ[t_] := 
      2 ArcSin[
          Sin[δ0/2] JacobiSN[ t, Sin[δ0/2]^2]]

T2	 periode for δ02
nc	numbers of pictures for the longer periode

T2 = 4 EllipticK[Sin[δ02/2]^2] // N;

    T = {T, T2};

    Tmax = Max[T];

    nc = 
        Round[Tmax bilderprosek] ;

coordinates of the pendulum body as a function of time

r[t_] = {Sin[δ[t]], -Cos[δ[t]]};

winkelzeitplot = 
        Plot[{Evaluate[δ[t][[1]]], 
            Evaluate[δ[t][[2]]]}, {t, 0, 
            Denominator[tzut] Tmax}, 
          PlotStyle ->  {Hue[0.0], Hue[0.6]}, 
          AxesLabel ->  {"Zeit t",  "Ausschlag"}];

zeitlinie = 
        Table[Graphics[{
Hue[0.95], 
              Thickness[
                gl/6], 
Line[{{k/nc Tmax, -π}, {k/
                        nc Tmax,  π}}]}, 
PlotRange ->  {{0, Tmax}, {-π, π}}], 
{k, 0, 
            nc Denominator[tzut] - 1}];

first deviation of δ(t)

dδ[t_] = D[δ[t], t];

curve of motion in the phase-space

phase1 = 
        ParametricPlot[{δ[t][[1]], dδ[
                t][[1]]}, {t, 0,  
            T[[1]]}, 
PlotStyle ->  {Hue[0]},  
          ImageSize  ->  550];
phase2 = 
        ParametricPlot[{δ[t][[2]], dδ[
                t][[2]]}, {t, 0,  
            T[[2]]}, 
PlotStyle ->  {Hue[0.7]}, 
          ImageSize  ->  550];

phase = 
        Graphics[Join[phase1[[1]], phase2[[1]]], 
          Axes ->  True];

coordinate of the pendulum in the phase-space as a function of time

phasenpkt1 = 
        Table[Graphics[{Hue[0], PointSize[gl], 
              Point[{δ[k/nc Tmax][[1]], dδ[
                      k/nc Tmax][[1]]}]}], {k, 0,  
            Denominator[tzut] nc - 1}];
phasenpkt2 = 
        Table[Graphics[{Hue[0.7], PointSize[gl], 
              Point[{δ[k/nc Tmax][[2]], dδ[
                      k/nc Tmax][[2]]}]}], {k, 0,  
            Denominator[tzut] nc - 1}];

put all pictures for the phase-space together

\(\(phase = 
        Table[Graphics[
            Join[phase1[\([1]\)], 
              phase2[\([1]\)], \(phasenpkt1[\([i]\)]\)[\([1]\)], \(phasenpkt2[\([i]\)]\)[\([1]\)]], 
            AxesLabel -> \ {"angle δ(t)", \ \*""(Overscript[δ, ●])(t)""}, Axes -> \ True], {i, 1, 
            Denominator[tzut]\ nc\ }];\)\)

put all pictures for the pendulum together

pendel1 = 
Table[
          Graphics[{Hue[ .0],  
              Thickness[
                gl/3], 
Line[{{0, 
                    0}, Transpose[
                      r[k/nc Tmax]][[1]]}], 
PointSize[gl],  GrayLevel[0], 
              Point[Transpose[
                    r[k/nc Tmax]][[1]]], 
Hue[ .6], 
              Thickness[
                gl/3], 
Line[{{0, 
                    0}, Transpose[
                      r[k/nc Tmax]][[2]]}], 
PointSize[gl],  
              GrayLevel[
                0], 
Point[Transpose[
                    r[k/nc Tmax]][[2]]]}, 
PlotRange -> {{-1, 1}, {-1, 1}}, 
            AspectRatio ->  Automatic], 
{k, 0, 
            Denominator[tzut] nc  - 1}];

simulation of the pendulum

pendel = 
Table[
          Show[GraphicsArray[
{{Graphics[
                    Join[winkelzeitplot[[1]], zeitlinie[[k, 1]]], 
                    winkelzeitplot[[2]]], 
                  pendel1[[k]]}, {phase[[k]]}}], 
AspectRatio  -> Automatic, 
            Prolog  ->  PointSize[ .01], 
            ImageSize ->  600], 
{k, 1, 
            Denominator[tzut] nc }];

    
SelectionMove[EvaluationNotebook[], After, 
        EvaluationCell];
SelectionMove[EvaluationNotebook[], Next, CellGroup];
FrontEndTokenExecute["OpenCloseGroup"];
start /.  zeit ->  0.2  // DisplayForm

Fourieanalysis

Parameters

Calculation

anzahl	numbers of amplitudes in δ0

δ0 = -Sort[-δ0]
anzahl = Dimensions[δ0][[1]]

T	periode

T = 4 EllipticK[Sin[δ0/2]^2] // N

spek[[i]]	fourier coefficients of the (i^th) pendulum
MaxRecursion is set to one to speed up the calculation, since the accuracy is not important in this case.

spek = {};
Do[
      spek = Join[spek, {0}], {i, 1, anzahl}]
Do[spek[[i]]=Table[NIntegrate[δ[t][[i]] Sin[
                      k 2 π/T[[i]] t],{t, 0, T[[i]]/4},Method->Trapezoidal,MaxRecursion ->  1],{k, 1, 7, 2}],
{i, 1, anzahl}]
Do[
        spek[[i]] = spek[[i]]/Max[spek[[i]]] , {i, 1, anzahl}] // 
      Chop

TableForm[spek]

breite	wide of the beam
legendex	starting coordinats for the legend

breite = 0.15;
legendex = 6;
\(\(fanalyse = 
        Graphics[{
Table[{Hue[\((j - 1)\)/
                    anzahl], 
Table[
                  Rectangle[{\((2\ i - 1)\) + \((j - 1)\)\ breite, 
                      spek[\([j, i]\)]}, {\((2\ i - 1)\) + 
                        0.1 + \((j - 1)\)\ breite, 0}], {i, 1, 
                    4}], 

Rectangle[{legendex, 0.88 - \((j - 1)\)\ 0.1}, {legendex + 0.3, 
                    0.92 - \((j - 1)\)\ 0.1}], 
                GrayLevel[0], 
Text[
                  StyleForm[
                    StringJoin["δ0= ", 
                      ToString[δ0[\([j]\)]/\((π)\)\ 180], "°"], \ FontSize -> 16], {legendex + 1.2, 
                    0.9 - 0.1\ \((j - 
                            1)\)}]
}, 
{j, 1, 
                anzahl}]
}];\)\)
Clear[breite, legendex];

Graphics

\(\(Show[fanalyse, Axes -> \ True, 
        AxesLabel -> \ {StyleForm["k", FontSize -> \ 16], 
            StyleForm[\ \*""\!\(a\_k\)/Max[\!\(a\_k\)]"", 
              FontSize -> \ 16]}, \ 
        PlotRange -> \ {{0, 7.8}, {0, 1}}, 
        ImageSize -> 600];\)\)