# How to numerically solve symmetric top with fixed point motion equations with Mathematica?

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?

• 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$ Commented Sep 8, 2023 at 12:05
• 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}. Commented Sep 9, 2023 at 2:17
• See also this topic mathematica.stackexchange.com/questions/232976/… Commented Sep 12, 2023 at 6:21

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.

## 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,
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],
Arrow[{{0, 0, 0}, bodyToInertia . {1, 0, 0}}]},
Inset[Graphics[
Text["x"]
], 1.1 bodyToInertia . {1, 0, 0}],

Arrow[{{0, 0, 0}, bodyToInertia . {0, 1, 0}}]},
Inset[Graphics[
Text["y"]
], 1.1 bodyToInertia . {0, 1, 0}],

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,
Arrow[{{0, 0, 0},
1.2 Norm[L2body] (omegaVectorBody/Norm[omegaVectorBody])}]}
],

If[showWC,
{
Arrow[{{0, 0, 0},
1.2 Norm[
L2body] ({omegaVectorBody[[1]], 0, 0}/
Norm[omegaVectorBody])}]},
Arrow[{{0, 0, 0},
1.2 Norm[
L2body] ({0, omegaVectorBody[[2]], 0}/
Norm[omegaVectorBody])}]},
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];

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"},
},
{"spin (hz)", "\!$$\*SubscriptBox[\(\[Theta]$$, $$1$$]\)",
"\!$$\*SubscriptBox[\(\[Theta]$$, $$2$$]\)", "\[Psi]'(hz)",
"\[Theta](t)", "n"},
{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;
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;
tick = Not[tick]} &], {1, 50, .1}, ImageSize -> Small],

{"Initial disk angle",
Manipulator[
Dynamic[theta0, {theta0 = #, currentTheta = theta0*Pi/180;
timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentPhiD = phiD0 2*Pi;
tick = Not[tick]} &], {1, 45, 1}, ImageSize -> Small],

{"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;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {1, 100, 1}, ImageSize -> Small],

Manipulator[
Dynamic[diskr, {diskr = #; timeHistoryIndex = 1;
currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {.1, 1, .1}, ImageSize -> Small],

{"disk thickness",
Manipulator[
Dynamic[diskThick, {diskThick = #; timeHistoryIndex = 1;
currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {.01, .6, .01}, ImageSize -> Small],

{"bar length",
Manipulator[
Dynamic[L0, {L0 = #; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {1, 2.5, .1}, ImageSize -> Small],

{"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},
{{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}],
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] &);
(*--------------------------------------------*)
f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"},
(*--------------------------------------------*)
f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"},
(*--------------------------------------------*)
f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"},
(*--------------------------------------------*)
f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"},
(*--------------------------------------------*)
)
(*reference: page 238, applied mechanics, vol 2, dynamics, by \
Housner and Hudon*)
]


• (+1) I'm always amazed by your manipulators. Commented Sep 18, 2023 at 20:51
• 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? Commented Sep 18, 2023 at 20:59
• @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./ Commented Sep 18, 2023 at 21:09
• @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. Commented Sep 18, 2023 at 21:31
• @E.Chan-López sorry, I have never used JSXGraph. I just googled it, and looks interesting. But I do not use javascript myself. Commented Sep 24, 2023 at 3:59

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.

• 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? Commented Sep 8, 2023 at 23:09
• 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. Commented Sep 9, 2023 at 6:02

For completeness, the heavy symmetric spinning top system can be solved using the "raw" system of equations deriving from its Lagrangian.

Without going into too much detail about deriving these equations, it suffices to properly define the Lagrangian function (which can be found in Goldstein's book and also some other books), and use the EulerEquations from the VariationalMethods package in Mathematica. Here is the MWE along with a graphical visualization:

Needs["VariationalMethods"]
motionEq =
EulerEquations[
I1/2 (\[Theta]'[t]^2 + \[Phi]'[t]^2 Sin[\[Theta][t]]^2) +
I3/2 (\[Psi]'[t] + \[Phi]'[t] Cos[\[Theta][t]])^2 -
m*g*l*Cos[\[Theta][t]], {\[Phi][t], \[Theta][t], \[Psi][t]}, t];
ics = {\[Theta][0.] == 43 Degree, \[Theta]'[0.] == 0., \[Phi][0.] ==
0, \[Phi]'[0.] == 1., \[Psi][0.] == 0., \[Psi]'[0.] == 5.};
{m, g, l, I1, I3, T} = {2., 9.81, 0.1, 2., 1., 50.};
sol = NDSolve[
Join[motionEq, ics], {\[Phi][t], \[Theta][t], \[Psi][t]}, {t, 0.,
T}];
\[Theta][t_] = \[Theta][t] /. sol[[1]];
\[Phi][t_] = \[Phi][t] /. sol[[1]];
\[Psi][t_] = \[Psi][t] /. sol[[1]];
X[t_] = {Sin[\[Theta][t]] Cos[\[Phi][t]],
Sin[\[Theta][t]] Sin[\[Phi][t]], Cos[\[Theta][t]]};
p1 = ParametricPlot3D[X[t], {t, 0, T}, PlotStyle -> {{Red, Thick}},
Boxed -> False, Axes -> False];
p2 = Graphics3D[{Opacity[0.5], Sphere[{0, 0, 0}, 1]}];
p3 = Graphics3D[{
Blue, Thick, Arrow[{{0, 0, 0}, {1.1, 0, 0}}],
Arrow[{{0, 0, 0}, {0, 1.1, 0}}], Arrow[{{0, 0, 0}, {0, 0, 1.1}}],
Text[
Style[
"\[Theta](0)=" <> ToString[\[Theta][0]/Degree] <> "\[Degree]",
Black, 15], {0, 0, 1.3}]
}];
Show[p1, p2, p3, PlotRange -> All]
`