2
$\begingroup$

I am trying to numerically solve the equations associated with the motion of a symmetric top with a fixed point. Details about this system can be found in Goldstein's Classical Mechanics book (3rd edition), pp. 208-223.

Here is the basic code that sets up NDSolve:

a = 6.86;
b = 1.48;
J1 = 1;
J3 = 2;
M = 0.5;
l = 0.05;
g = 9.8;
H = 12.009;
NDSolve[
   {\[Phi]'[t] == (b - a Cos[\[Theta][t]])/Sin[\[Theta][t]]^2,
    \[Psi]'[t] == (J1 a)/J3 - Cos[\[Theta][t]] (b - a Cos[\[Theta][t]])/Sin[\[Theta][t]]^2,
    H == (J1^2 a^2)/(2 J3) + g l M Cos[\[Theta][t]] + 1/2 J1 (\[Theta]'[t]^2 + \[Phi]'[t]^2 Sin[\[Theta][t]]^2),
    \[Phi]'[0.] == 0., \[Theta]'[0.] == 0., \[Theta][0.] == 0.2},
   {\[Phi][t], \[Theta][t], \[Psi][t]},
   {t, 0., 1.}, MaxSteps -> 10^15, AccuracyGoal -> 24
]

Now, I know that this system can be solved through a thorough analysis (like the one presented by Goldstein himself) and that there are Demonstrations of this within Wolfram's community (see here for example).

However, I would like to be able to simply solve the "raw" system of equations (without altering it too much) with numerical methods. Notice that I increased AccuracyGoal and MaxSteps numbers with no effect. What else can be tried?

Improve this question
$\endgroup$
3
  • $\begingroup$ Properly set with start values for all three angles (not derivatives) at t=0 and WorkingPrecision->32, the NIntegrate stops at t=6*10^-7. With a square root , the equations are not necessarily Lipshitz. E.g $$x' =\sqrt(1-x^2)$$ at $x=+-1$ $\endgroup$
    – Roland F
    Sep 8 at 12:05
  • $\begingroup$ Initial conditions should be in agreement with equations. But in your case we have {\[Phi]'[t] == (b - a Cos[\[Theta][t]])/Sin[\[Theta][t]]^2, \[Psi]'[t] == (J1 a)/J3 - Cos[\[Theta][t]] (b - a Cos[\[Theta][t]])/Sin[\[Theta][t]]^2, H == (J1^2 a^2)/(2 J3) + g l M Cos[\[Theta][t]] + 1/2 J1 (\[Theta]'[t]^2 + \[Phi]'[t]^2 Sin[\[Theta][t]]^2)} /. {\[Phi]'[t] ->0., \[Theta]'[t] -> 0., \[Theta][t] -> 0.2} out {False, Derivative[1][\[Psi]][t] == 133.625, False}. $\endgroup$
    – Alex Trounev
    Sep 9 at 2:17
  • $\begingroup$ See also this topic mathematica.stackexchange.com/questions/232976/… $\endgroup$
    – Alex Trounev
    Sep 12 at 6:21

2 Answers 2

Reset to default
2
$\begingroup$

This is Manipulate of Symmetric top gyroscope motion. It has many options to explore the dynamics. I wanted to improve it more, but no time now. But it works and has many options to play with.

enter image description here

enter image description here

Code

Manipulate[
 (*version 9/18/2023*)
 tick;
 Module[{t, phi, theta, zeta, L2, thetas, phis, zetas, thetaW, phiW, 
   zetaW, graph, g, sol, eq1, eq2, eq3, Iz, I0, x, mass1, mass2, cg, 
   coneOff = 10, distToGroud, ic, omegaVector, bodyToInertia, yBody, 
   zBody, L0body, L2body, cgbody, PE, KE, initialSpin, N0, beta, 
   theta2, n, tmp, omegaVectorBody, inertiaToBody, 
   stripLinesSideOfRotor, stripLinesTopOfRotor, plots},
  
  initialSpin = phiD0  (2*Pi);
  mass1 = density*Pi diskr^2*diskThick;(*mass of rotor*)
  mass2 = density*Pi (diskr/10)^2*L0;(*mass of bar*)
  g = 9.81;
  cg = ( mass2*L0/2 + mass1*(L0 + diskThick/2))/(mass1 + mass2);
  I0 = 1/12 mass1 (3*diskr^2 + diskThick^2) + 
    mass1* (L0 + diskThick/2)^2;(*parallel axes*)
  I0 += 1/12 mass2 (3*(diskr/10)^2 + (L0)^2) + mass2*(L0/2)^2;
  Iz = mass1*diskr^2/2;
  Iz += mass2*(diskr/10)^2/2;
  
  eq1 = (mass1 + mass2)*g*(cg) Sin[theta[t]] == 
    I0 theta''[t] + 
     Iz (phi'[t] + zeta'[t] Cos[theta[t]]) zeta'[t] Sin[theta[t]] - 
     I0 (zeta'[t])^2 Sin[theta[t]] Cos[theta[t]];
  (*changed sign*)
  eq2 = 0 == 
    I0 (zeta''[t] Sin[theta[t]] + 
        2 zeta'[t] theta'[t] Cos[theta[t]]) - 
     Iz theta'[t] (phi'[t] + zeta'[t] Cos[theta[t]]);
  eq3 = 0 == 
    Iz (phi''[t] + zeta''[t] Cos[theta[t]] - 
       zeta'[t] theta'[t] Sin[theta[t]]);
  
  ic = {theta[0] == currentTheta, theta'[0] == currentThetaD, 
    zeta[0] == currentZeta, zeta'[0] == currentZetaD, 
    phi[0] == currentPhi, phi'[0] == currentPhiD};
  sol = First@
    NDSolve[Flatten@{eq1, eq2, eq3, ic}, {theta, zeta, phi, theta', 
      zeta', phi'}, {t, 0, delT}, Method -> "StiffnessSwitching"];
  
  thetas = (theta /. sol)[0];
  zetas = (zeta /. sol)[0];
  phis = (phi /. sol)[0];
  thetaW = (theta' /. sol)[0];
  zetaW = (zeta' /. sol)[0];
  phiW = (phi' /. sol)[0];
  
  n = (phiW + zetaW Cos[thetas]);
  N0 = (Iz/I0) n ;
  beta = 2 (mass1 + mass2)*g*cg/I0;
  tmp = N0^2/(2 beta);
  theta2 = 
   ArcCos[tmp - 
     Sqrt[1 - 2 tmp Cos[theta0 Degree] + tmp^2]];(*maximum nutation*)
  
  (*this transformation for only Euler angles zeta and theta (first \
two) *)
  bodyToInertia = {{Cos[zetas], -Cos[thetas] Sin[zetas], 
     Sin[thetas] Sin[zetas]}, {Sin[zetas], 
     Cos[thetas] Cos[zetas], -Cos[zetas] Sin[thetas]}, {0, 
     Sin[thetas], Cos[thetas]}};
  L2 = L0 + diskThick;
  
  L0body = {0, 0, L0};
  L2body = {0, 0, L2};
  cgbody = {0, 0, cg};
  distToGroud = L0 Cos[thetas] - diskr Sin[thetas];
  
  omegaVector = {thetaW , zetaW Sin[thetas], n};
  omegaVectorBody = bodyToInertia . omegaVector;
  stripLinesSideOfRotor = {bodyToInertia . {diskr Cos[#], 
        diskr Sin[#], L2}, 
      bodyToInertia . {diskr Cos[#], diskr Sin[#], L0}} & /@ 
    Range[0, 2 Pi, Pi/4];
  stripLinesTopOfRotor = {bodyToInertia . {diskr Cos[#], diskr Sin[#],
         L2}, bodyToInertia . {0, 0, L2}} & /@ Range[0, 2 Pi, Pi/4];
  
  graph = Graphics3D[
    {
     Rotate[GraphicsGroup[
       {
        {Lighting -> "Neutral", FaceForm[LightGray], Opacity[op], 
         Cylinder[{{bodyToInertia . L0body, 
            bodyToInertia . L2body}}, .98 diskr]},
        {Blue, Line[stripLinesSideOfRotor]},
        {Blue, Line[stripLinesTopOfRotor]},
        
        {Lighting -> "Neutral", FaceForm[LightGray], 
         Cylinder[{bodyToInertia . L0body/coneOff , 
           bodyToInertia . L0body}, diskr/10]}, (*rod*)
        Cone[{{(bodyToInertia . L0body/coneOff), {0, 0, 0}}}, 
         diskr/10],(*botton cone*)
        If[
         showCG, {Red, 
          Sphere[bodyToInertia . {0, 0, cg}, 1.1 (diskr/10)]}]
        }], phis, bodyToInertia . {0, 0, 1}
      ],
     
     If[showTorque,
      {Red, 
       Arrow[{bodyToInertia . {0, 0, cg}, 
         bodyToInertia . {L0, 0, cg}}]}
      ],
     
     If[showH,
      {Black, Arrowheads[Medium], 
       Arrow[{{0, 0, 0}, 
         1.3 Norm[L2body] (omegaVectorBody/Norm[omegaVectorBody])}]}
      ],
     
     If[showBodyAxes,
      {
       Cylinder[{bodyToInertia . {0, 0, -.2}, 
         bodyToInertia . {0, 0, -.1 .5}}, 1.5],
       {Red, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, bodyToInertia . {1, 0, 0}}]},
       Inset[Graphics[
         Text["x"]
         ], 1.1 bodyToInertia . {1, 0, 0}],
       
       {Red, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, bodyToInertia . {0, 1, 0}}]},
       Inset[Graphics[
         Text["y"]
         ], 1.1 bodyToInertia . {0, 1, 0}],
       
       {Red, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, bodyToInertia . {0, 0, 1}}]},
       Inset[Graphics[
         Text["z"]
         ], 1.1 bodyToInertia . {0, 0, 1}]
       }],
     
     If[showPath,
      {Red, Line[tipTimeHistory[[1 ;; timeHistoryIndex]]]}
      ],
     
     If[showInertiaAxes,
      {
       {Arrowheads[Small], Arrow[{{0, 0, 0}, {1, 0, 0}}]},
       Inset[Graphics[
         Text["x"]
         ], {1.1, 0, 0}],
       
       {Arrowheads[Small], Arrow[{{0, 0, 0}, {0, 1, 0}}]},
       Inset[Graphics[
         Text["y"]
         ], {0, 1.1, 0}],
       
       {Arrowheads[Small], Arrow[{{0, 0, 0}, {0, 0, 1}}]},
       
       Inset[Graphics[
         Text["z"]
         ], {0, 0, 1.1}]
       }],
     
     If[showW,
      {Red, Arrowheads[Medium], 
       Arrow[{{0, 0, 0}, 
         1.2 Norm[L2body] (omegaVectorBody/Norm[omegaVectorBody])}]}
      ],
     
     If[showWC,
      {
       {Blue, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, 
          1.2 Norm[
            L2body] ({omegaVectorBody[[1]], 0, 0}/
             Norm[omegaVectorBody])}]},
       {Blue, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, 
          1.2 Norm[
            L2body] ({0, omegaVectorBody[[2]], 0}/
             Norm[omegaVectorBody])}]},
       {Blue, Arrowheads[Small], 
        Arrow[{{0, 0, 0}, 
          1.2 Norm[
            L2body] ({0, 0, omegaVectorBody[[3]]}/
             Norm[omegaVectorBody])}]}
       }
      ],
     
     {Lighting -> "Neutral", FaceForm[LightGray], 
      Cylinder[{{0, 0, -.2}, {0, 0, -.1 .5}}, 1.79 L0]},
     
     If[showSphere,
      {
       {Opacity[.1], Sphere[{0, 0, 0}, 1.79 L0]},
       
       {Red, FaceForm[None], 
        Cylinder[{{0, 0, 1.78 L0 Cos[theta0 Degree]}, {0, 0, 
           1.79 L0 Cos[theta0 Degree]}}, 
         1.79 L0 Sin[theta0 Degree]]},
       
       {Red, FaceForm[None], 
        Cylinder[{{0, 0, 1.78 L0 Cos[theta2]}, {0, 0, 
           1.79 L0 Cos[theta2]}}, 1.79 L0 Sin[theta2]]}
       }
      ],
     
     If[distToGroud <= 0, Text[Style["Crash!", 14, Red], {3, 3, 0}]]
     },
    PlotRange -> {{-zoom, zoom}, {-zoom, zoom}, {-.5, 1.9 L0}},
    SphericalRegion -> True, ImagePadding -> .1, ImageMargins -> 0, 
    If[showPlots, ImageSize -> 250, ImageSize -> 350], 
    ViewPoint -> viewPoint, Boxed -> True
    ];
  
  Which[state == "RUN" || state == "STEP",
   currentTheta = (theta /. sol)[delT];
   currentZeta = (zeta /. sol)[delT];
   currentPhi = (phi /. sol)[delT];
   
   currentThetaD = (theta' /. sol)[delT];
   currentZetaD = (zeta' /. sol)[delT];
   currentPhiD = (phi' /. sol)[delT];
   
   If[currentTime > 0, timeHistoryIndex++];
   thetaTimeHistory[[timeHistoryIndex]] = {currentTime, thetas*180/Pi};
   zetaTimeHistory[[timeHistoryIndex]] = {currentTime, zetas*180/Pi};
   tmp = (bodyToInertia . {0, 0, cg});
   tipTimeHistory[[timeHistoryIndex]] = 1.8 L0 (tmp/Norm[tmp]);
   
   If[timeHistoryIndex == 500,
    timeHistoryIndex = 1;
    currentTime = 0
    ,
    currentTime = currentTime + delT
    ];
   
   If[state == "RUN" && distToGroud > 0,
    tick = Not[tick]
    ]
   ];
  
  PE = (mass1 + mass2)*g*(cg) Cos[thetas];
  KE = 1/2 Iz (initialSpin + zetaW Cos[thetas])^2 + 
    I0 (thetaW^2 + zetaW^2 Sin[thetas]^2);
  If[showPlots, 
   plots = Row[{ListLinePlot[ 
       thetaTimeHistory[[1 ;; timeHistoryIndex]], 
       PlotRange -> {{0, 500*delT}, {0.9 theta0, 1.1 theta2*180/Pi}} ,
        Frame -> True, AspectRatio -> 0.3, 
       FrameLabel -> {{"\[Theta]", None}, {"time (sec)", 
          "Nutation angle (deg) vs. time"}}, ImageSize -> 250, 
       ImagePadding -> {{45, 10}, {30, 45}}, ImageMargins -> 0],
      ListLinePlot[ zetaTimeHistory[[1 ;; timeHistoryIndex]], 
       PlotRange -> {{0, 500*delT}, {0, 360}} , Frame -> True, 
       AspectRatio -> 0.3, 
       FrameLabel -> {{"\[Psi]", None}, {"time (sec)", 
          "precession angle (deg) vs. time"}}, ImageSize -> 250, 
       ImagePadding -> {{45, 10}, {30, 45}}, ImageMargins -> 0]
      }]
   ];
  
  Grid[{
    {Grid[{
       {"c.g.", "w1", "w2", "w3", "P.E.", "K.E.", "total energy"},
       {padIt2[cg, {3, 2}],
        padIt2[omegaVector[[1]], {5, 2}],
        padIt2[omegaVector[[2]], {5, 2}],
        padIt2[omegaVector[[3]], {5, 2}],
        padIt1[PE, {9, 2}],
        padIt2[KE, {9, 2}],
        padIt1[PE + KE, {9, 2}]
        },
       {"spin (hz)", "\!\(\*SubscriptBox[\(\[Theta]\), \(1\)]\)", 
        "\!\(\*SubscriptBox[\(\[Theta]\), \(2\)]\)", "\[Psi]'(hz)", 
        "\[Theta](t)", "n"},
       {padIt1[phiW/(2*Pi), {5, 2}],
        padIt2[theta0, {4, 2}],
        padIt2[theta2*180/Pi, {4, 2}],
        padIt1[zetaW/(2*Pi), {4, 2}],
        padIt1[currentTheta*180/Pi, {4, 2}],
        padIt1[n, {4, 2}]},
       {If[showPlots,
         Grid[{
           {plots},
           {graph}
           }],
         Grid[{
           {graph}
           }]
         ], SpanFromLeft
        }
       }, Frame -> All, FrameStyle -> Gray]
     }
    }]],
 
 Grid[{
   {
    Grid[{
      { Button[
        Text@Style["run", 12], {state = "RUN"; tick = Not[tick]}, 
        ImageSize -> {50, 40}],
       Button[
        Text@Style["step", 12], {state = "STEP"; tick = Not[tick]}, 
        ImageSize -> {50, 40}],
       Button[
        Text@Style["stop", 12], {state = "STOP"; tick = Not[tick]}, 
        ImageSize -> {50, 40}],
       Button[Text@Style["reset", 12], {state = "RESET"; delT = 0.01;
         density = 1;
         diskr = 1;
         diskThick = .6;
         L0 = 2.5;
         showCG = True;
         nStrips = 1;
         currentTheta = 35*Pi/180;
         theta0 = 35;
         currentZeta = 0;
         currentPhi = 0;
         currentThetaD = 0;
         currentZetaD = 0;
         currentPhiD = 2 Pi*15;
         phiD0 = 15;
         showWC = False;
         showW = False;
         showInertiaAxes = True;
         showBodyAxes = True;
         delT = 0.035;
         timeHistoryIndex = 1;
         showPlots = False;
         currentTime = 0;
         showPath = True;
         showH = False;
         showTorque = False;
         tick = Not[tick]}, ImageSize -> {50, 40}](*fix*)
       }
      }, Spacings -> {.5, 0}, Frame -> True, FrameStyle -> Gray
     ]
    },
   {
    Text@Grid[{
       {"Initial disk spin (hz)", 
        Manipulator[
         Dynamic[phiD0, {phiD0 = #, currentPhiD = phiD0 2*Pi; 
            timeHistoryIndex = 1; currentTime = 0;
            currentTheta = theta0*Pi/180;
            currentZeta = 0;
            currentPhi = 0;
            currentThetaD = 0;
            currentZetaD = 0;
            tick = Not[tick]} &], {1, 50, .1}, ImageSize -> Small], 
        Dynamic[padIt1[phiD0, {3, 1}]]},
       
       {"Initial disk angle", 
        Manipulator[
         Dynamic[theta0, {theta0 = #, currentTheta = theta0*Pi/180; 
            timeHistoryIndex = 1; currentTime = 0;
            currentZeta = 0;
            currentPhi = 0;
            currentThetaD = 0;
            currentZetaD = 0;
            currentPhiD = phiD0 2*Pi;
            tick = Not[tick]} &], {1, 45, 1}, ImageSize -> Small], 
        Dynamic[padIt1[theta0, 2]]},
       
       {"simulation time step", 
        Manipulator[
         Dynamic[delT, {delT = #; timeHistoryIndex = 1; 
            currentTime = 0; tick = Not[tick]} &], {.001, 0.05, .001},
          ImageSize -> Small], Dynamic[padIt1[delT, {4, 3}]]},
       
       {"disk density", 
        Manipulator[
         Dynamic[density, {density = #; timeHistoryIndex = 1; 
            currentTime = 0;
            currentZeta = 0;
            currentPhi = 0;
            currentTheta = theta0*Pi/180;
            
            currentZetaD = 0;
            currentPhiD = phiD0 2*Pi;
            currentThetaD = 0;
            
            tick = Not[tick]} &], {1, 100, 1}, ImageSize -> Small], 
        Dynamic[padIt1[density, 3]]},
       
       {"disk radius", 
        Manipulator[
         Dynamic[diskr, {diskr = #; timeHistoryIndex = 1; 
            currentTime = 0;
            currentZeta = 0;
            currentPhi = 0;
            currentTheta = theta0*Pi/180;
            
            currentZetaD = 0;
            currentPhiD = phiD0 2*Pi;
            currentThetaD = 0;
            
            tick = Not[tick]} &], {.1, 1, .1}, ImageSize -> Small], 
        Dynamic[padIt1[diskr, {2, 1}]]},
       
       {"disk thickness", 
        Manipulator[
         Dynamic[diskThick, {diskThick = #; timeHistoryIndex = 1; 
            currentTime = 0;
            currentZeta = 0;
            currentPhi = 0;
            currentTheta = theta0*Pi/180;
            
            currentZetaD = 0;
            currentPhiD = phiD0 2*Pi;
            currentThetaD = 0;
            
            tick = Not[tick]} &], {.01, .6, .01}, ImageSize -> Small],
         Dynamic[padIt1[diskThick, {2, 1}]]},
       
       {"bar length", 
        Manipulator[
         Dynamic[L0, {L0 = #; timeHistoryIndex = 1; currentTime = 0;
            currentZeta = 0;
            currentPhi = 0;
            currentTheta = theta0*Pi/180;
            
            currentZetaD = 0;
            currentPhiD = phiD0 2*Pi;
            currentThetaD = 0;
            
            tick = Not[tick]} &], {1, 2.5, .1}, ImageSize -> Small], 
        Dynamic[padIt1[L0, {2, 1}]]},
       
       {"rotor opacity", 
        Manipulator[
         Dynamic[op, {op = #, tick = Not[tick]} &], {0, 1, .01}, 
         ImageSize -> Small], Dynamic[padIt1[op, {2, 1}]]},
       {"select viewpoint",
        SetterBar[
         Dynamic[viewPoint, {viewPoint = #, tick = Not[tick]} &], {{2,
             0, 3} -> 1, {1, -2, 1} -> 2, {0, -2, 2} -> 
           3, {-2, -2, 0} -> 4, {2, -2, 0} -> 5, {Pi, Pi/2, 2} -> 6}
         ]
        },
       {"zoom",
        Manipulator[
         Dynamic[zoom, {zoom = #, tick = Not[tick]} &], {1, 8, .1}, 
         ImageSize -> Small],
        ""
        },
       {Grid[{
          {"show c.g.", 
           Checkbox[Dynamic[showCG, {showCG = #; tick = Not[tick]} &]],
           "show w components", 
           Checkbox[
            Dynamic[showWC, {showWC = #; tick = Not[tick]} &]]},
          {"show w", 
           Checkbox[Dynamic[showW, {showW = #; tick = Not[tick]} &]],
           "show inertial axes", 
           Checkbox[
            Dynamic[showInertiaAxes, {showInertiaAxes = #; 
               tick = Not[tick]} &]]},
          {"show body axes", 
           Checkbox[
            Dynamic[showBodyAxes, {showBodyAxes = #; 
               tick = Not[tick]} &]],
           "show plots", 
           Checkbox[
            Dynamic[showPlots, {showPlots = #; 
               tick = Not[tick]} &]]},
          {"show space path", 
           Checkbox[
            Dynamic[showPath, {showPath = #; tick = Not[tick]} &]],
           "show sphere", 
           Checkbox[
            Dynamic[showSphere, {showSphere = #; 
               tick = Not[tick]} &]]},
          {"show torque vector", 
           Checkbox[
            Dynamic[showTorque, {showTorque = #; tick = Not[tick]} &]],
           "show angular momentum", 
           Checkbox[Dynamic[showH, {showH = #; tick = Not[tick]} &]]}
          }, Frame -> All, FrameStyle -> Gray], SpanFromLeft
        }
       }, Frame -> True, FrameStyle -> Gray, Alignment -> Left
      ]
    }
   }, Alignment -> Left],
 
 {{showH, False}, None},
 {{showTorque, False}, None},
 {{showPath, True}, None},
 {{showSphere, True}, None},
 {{tick, False}, None},
 {{showPlots, False}, None},
 {{currentTime, 0}, None},
 {{density, 1}, None},
 {{diskr, 1}, None},
 {{diskThick, .6}, None},
 {{L0, 2.5}, None},
 {{state, "STOP"}, None},
 {{showCG, True}, None},
 {{nStrips, 1}, None},
 {{zoom, 4.6}, None},
 {{op, 1}, None},
 {{showWC, False}, None},
 {{showW, False}, None},
 {{showInertiaAxes, True}, None},
 {{showBodyAxes, True}, None},
 
 {{theta0, 35}, None},
 {{currentTheta, 35.0*Pi/180}, None},
 {{currentZeta, 0}, None},
 {{currentPhi, 0}, None},
 {{currentThetaD, 0}, None},
 {{currentZetaD, 0}, None},
 {{currentPhiD, 2 Pi*13}, None},
 {{phiD0, 13}, None},
 
 {{viewPoint, {Pi, Pi/2, 2}}, None},
 {{delT, 0.035}, None},
 {{thetaTimeHistory, Table[{0, 0}, {500}]}, None},
 {{zetaTimeHistory, Table[{0, 0}, {500}]}, None},
 {{tipTimeHistory, Table[{0, 0, 0}, {500}]}, None},
 {{timeHistoryIndex, 1}, None},
 
 SynchronousUpdating -> False,
 ControlPlacement -> Left, Alignment -> Center, ImageMargins -> 0, 
 FrameMargins -> 0,
 TrackedSymbols :> {tick},
 Initialization :>
  (
   mplt = 
    MatrixPlot[Table[Sin[x y/100], {x, -5, 5}, {y, -5, 5}], 
     Frame -> False, ImagePadding -> 0, PlotRangePadding -> 0];
   integerStrictPositive = (IntegerQ[#] && # > 0 &);
   integerPositive = (IntegerQ[#] && # >= 0 &);
   numericStrictPositive = (Element[#, Reals] && # > 0 &);
   numericPositive = (Element[#, Reals] && # >= 0 &);
   numericStrictNegative = (Element[#, Reals] && # < 0 &);
   numericNegative = (Element[#, Reals] && # <= 0 &);
   bool = (Element[#, Booleans] &);
   numeric = (Element[#, Reals] &);
   integer = (Element[#, Integers] &);
   (*--------------------------------------------*)
   padIt1[v_?numeric, f_List] := AccountingForm[v,
     f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   padIt1[v_?numeric, f_Integer] := AccountingForm[Chop[v],
     f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   padIt2[v_?numeric, f_List] := AccountingForm[v,
     f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   padIt2[v_?numeric, f_Integer] := AccountingForm[Chop[v],
     f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True]
   (*--------------------------------------------*)
   )
 (*reference: page 238, applied mechanics, vol 2, dynamics, by \
Housner and Hudon*)
 ]

Link to notebook is here

Improve this answer
$\endgroup$
8
  • $\begingroup$ (+1) I'm always amazed by your manipulators. $\endgroup$
    – E. Chan-López
    Sep 18 at 20:51
  • $\begingroup$ My laptop couldn't handle running for long without getting stuck, what minimum specifications do you think a computer should have to be able to run your manipulator without problems? $\endgroup$
    – E. Chan-López
    Sep 18 at 20:59
  • 1
    $\begingroup$ @E.Chan-López sorry to hear that. I have no problem running this for hrs. I have 128GB RAM and very fast new intel CPU, all on windows 10. I know Mathematica 3D uses lots of RAM and this one does need more RAM I suppose./ $\endgroup$
    – Nasser
    Sep 18 at 21:09
  • 1
    $\begingroup$ @E.Chan-López Once I update to my current PC, yes, I rarely have memory problems. one day, I'd like to get 256 GB RAM PC. If you do lots of 3D graphics, memory is really the most important thing to consider. $\endgroup$
    – Nasser
    Sep 18 at 21:31
  • 1
    $\begingroup$ @E.Chan-López sorry, I have never used JSXGraph. I just googled it, and looks interesting. But I do not use javascript myself. $\endgroup$
    – Nasser
    Sep 24 at 3:59
1
$\begingroup$

The use of the Hamilton function (energy) in terms of the velocities for deriving equations of motion is always a dangerous mistake.

Use the Lagrangian equations (trivial with two constant angular momentum variables corresponding to cyclic $\phi, \psi$) or replace all angular velocities with their conjugate angular momentum in the Hamiltonian and use canonical equations (again trivial with two constant momentum variabled).

In both cases there is a second order ODE for $\theta$ depending an the two constant angular momentums and two integrals for $\phi, \psi$ dependent in the solution for $\theta(t)$.

The direct use of the square root runs into being non-Lipshitz continuous, like e.g. $x'^2 = 1-x^2 $ with solutions $x = \sin (t+c),\ x=1,\ x=-1$ with undecided continuations at common tangents.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ I will do your proposed procedure. However, I am interested in your first comment: "The use of the Hamilton function (energy) in terms of the velocities for deriving equations of motion is always a dangerous mistake." Why? Could you please elaborate on this? Is it related with the energy equation being non-Lipschitz continuous? $\endgroup$
    – Meclassic
    Sep 8 at 23:09
  • 1
    $\begingroup$ The energy equation $v^2 + 2V - 2 E ==0 $ yields the ODE for the inverse function of $t\to x(t)$ by separation of variables $dt == dx/\sqrt(2 (E-V(x)))$ that is ingrable only on an interval of monotony. For the spinning heavy top it produces elliptic integrals with inverses being complicated Jacobi elliptic functions. This is the topic of Klein/Sommerfeld's monumental monography on the Theory of the Spinning Top from about 1900. $\endgroup$
    – Roland F
    Sep 9 at 6:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.