Motion of a Spinning Top

Jan Vrbik

Both approximate and exact solutions for the motion of a spinning top are constructed with the help of quaternions.

Quaternion Algebra

A quaternion is a four-dimensional quantity consisting of a scalar, say A, and a three-dimensional vector a, collectively denoted 𝔸≡(A,a). Addition of two quaternions is component-wise,

(A,a)⊕(B,b)=(A+B,a+b),

(we do not need to add quaternions in this article). Their multiplication follows the rule

(A,a)⊗(B,b)=(A B-a·b ,A b+B a-a×b).

It is important to note that such multiplication is associative (even though noncommutative). This can be verified by the following.

{A_,a_}∘{B_,b_}:={A B-a.b,A b+B a-a\[Cross]b}//TrigReduce//Simplify

({A,{(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3])}}∘{B,{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}})∘{C,{(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])}}-{A,{(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3])}}∘({B,{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}}∘{C,{(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])}})//Expand

When the scalar part of a quaternion is zero, it is called a pure quaternion. Since this constitutes an important special case, we will extend our definition of multiplication accordingly.

{A_,a_}∘b_:={-a.b,A b-a\[Cross]b}//TrigReduce//Simplify

The operation of conjugation simply changes the sign of the vector part of a quaternion.

(OverBar[{A_,a_}]):={A,-a}

Based on the corresponding Taylor expansion, it is possible to evaluate various functions of a quaternion. Of these, the most important is the exponential function, which, in this article, we need only with pure-quaternion arguments.

PQExp[a_]:=Module[{n=(Sqrt[a.a])//PowerExpand},{Cos[n],a /n Sin[n]}]

Rotation

Rotating a 3D vector x with respect to an axis through the origin (with a unit direction of (OverHat[α])) by an angle a can be achieved by

x'=(e^(-α/2))∘x∘(e^(α/2))=(cos (a/2),-(OverHat[α]) sin(( a)/2))∘((Subscript[x, ∥])⊕(Subscript[x, ⊥]))∘(cos(a/2),(OverHat[α]) sin (a/2))=(Subscript[x, ∥])⊕(Subscript[x, ⊥])∘(cos a, (OverHat[α]) sin a)=(Subscript[x, ∥])⊕((Subscript[x, ⊥])cos a+α×(Subscript[x, ⊥]) sin a),

where α=a(OverHat[α]) (a vector of length a and unit direction (OverHat[α])), and (Subscript[x, ∥]) and (Subscript[x, ⊥]) are the parallel and perpendicular (to α) parts of x, respectively []. Note that (Subscript[x, ∥]) and α commute, while (Subscript[x, ⊥]) and α anticommute.

One can show that a quaternion ℝ has the form of (e^(α/2)) if and only if (OverBar[ℝ])∘ℝ≡ℝ∘(OverBar[ℝ])=1 (a pure scalar). Alternately, the same ℝ can be parametrized by the three Euler angles (see the next section).

By taking α=t ω, where t is time, we achieve constant rotation with an uniform angular velocity ω. Note that

ω=2(OverBar[ℝ])∘(OverDot[ℝ]),

where ℝ=(e^((ω/2)t)), and (OverDot[ℝ]) denotes its time derivative.

The previous formula can be generalized to any time-dependent ℝ, as can be seen by differentiating x'=(OverBar[ℝ])∘x∘ℝ, thus

(OverDot[x])'=(OverBar[ℝ])∘x∘(OverDot[ℝ])+(Overscript[ℝ, (OverDot[_])])∘x∘ℝ=(OverBar[ℝ])∘x∘ℝ∘(OverBar[ℝ])∘(OverDot[ℝ])+(Overscript[ℝ, (OverDot[_])])∘ℝ∘(OverBar[ℝ])∘x∘ℝ=(((OverBar[ℝ])∘x∘ℝ∘ω+(OverBar[ω])∘(OverBar[ℝ])∘x∘ℝ)/2)=((x'∘ω-ω∘x')/2)=ω×x',

where ω is now the instantaneous angular velocity at time t, which may change with time. Note that (OverBar[ℝ])∘(OverDot[ℝ]) must be a pure quaternion (take the time derivative of (OverBar[ℝ])∘ℝ=1 to see that). Now, any rotation ℝ can be applied to the axes of the original (inertial) system of coordinates, thus creating the corresponding new, rotating (noninertial) frame. Components of any vector x can then be expressed in either the old, or (often and more conveniently) the new coordinates; the latter will be given by ℝ∘x∘(OverBar[ℝ]). The angular velocity ω, in the new coordinates, will thus have components given by 2(OverDot[ℝ])∘(OverBar[ℝ]). This becomes an important and helpful tool: even though the laws of physics are normally valid only in inertial coordinates, the resulting equations often simplify when expressed in the rotating frame.

Spinning Top

Let us consider an axially symmetrical body (a “top”) of mass M, and parallel and perpendicular (with respect to the body’s symmetry axis) moments of inertia equal to (Subscript[J, ∥]) and (Subscript[J, ⊥]), respectively. We make it spin around its axis, place the bottom tip of its (usually tilted) axis on a nonslip horizontal plane (a desk), and let it continue its motion subject to the vertical gravity field. Let us also assume that the top’s center of mass is at a distance r from its point of contact, that the gravitational acceleration is g, and that there is no friction.

To investigate the top’s motion, we introduce an inertial coordinate system with its origin at the point of contact and its z direction pointing vertically upward. First we position a motionless top vertically (its symmetry axis aligned with the z direction) on the desk’s surface. We then rotate it using a fully general, time-dependent rotation ℝ, parametrized by three Euler angles, thus

ℝ=(e^(𝕜 ϕ/2))∘(e^(𝕚 θ/2))∘(e^(𝕜 ψ/2)),

where ϕ is a rapidly increasing function of time that represents the actual spinning of the top around its symmetry axis (this rotation is applied first), followed by tilting the spinning top by an angle θ (potentially a slowly varying function of time) around the inertial x direction, further followed by slowly rotating the spinning, tilted top around the inertial z direction by an angle ψ, thus creating a so-called precession. Also: 𝕚, 𝕛 (not used in the last formula), and 𝕜 are unit vectors along the x, y, and z directions, respectively, in their pure-quaternion form.

Equations of Motion

We now use the Lagrange technique to find the corresponding equations for the three Euler angles. The kinetic energy of the top equals one-half of the scalar product of its angular velocity with its angular momentum, as its motion is purely rotational. In the rotating frame (this is when it becomes handy), this equals

(((Subscript[J, ∥])((Subscript[ω, z'])^2)+((Subscript[J, ⊥])+M (r^2))(((Subscript[ω, x'])^2)+((Subscript[ω, y'])^2)))/2),

where (Subscript[ω, x']), (Subscript[ω, y']), (Subscript[ω, z']) are the rotating-frame coordinates of ω, and (Subscript[J, ⊥])+M (r^2) (denoted J from now on; similarly, (Subscript[J, ∥]) will be called (Subscript[J, 3])) is the top’s moment of inertia with respect to a line perpendicular to its symmetry axis, and passing through the origin (which is displaced by r from its center of mass). Usually J>(Subscript[J, 3]) (a thin circular disk needs r to be bigger than only one-half of its radius to meet this condition), but our results are valid when J≤(Subscript[J, 3]) as well.

Subtracting the potential energy M g r cos θ≡J m cos θ (by the definition of m) yields the resulting Lagrangian.

𝕚={1,0,0}; 𝕛={0,1,0}; 𝕜={0,0,1};
R=(PQExp[𝕜/2 ϕ[t]]∘PQExp[𝕚/2 θ[t]])∘PQExp[𝕜/2 ψ[t]];
ω=(2D[R,t]∘(OverBar[R]))[[2]];
L={J,J,(Subscript[J, 3])}.(ω^2)/2-J m Cos[θ[t]]//Simplify

Since the result is free of ϕ and ψ (it contains only their time derivatives), the corresponding ((∂L)/(∂(OverDot[ϕ]))) and ((∂L)/(∂(OverDot[ψ]))) must be constants of motion []. The remaining equation is then obtained from

(d/(dθ))(((∂L)/(∂(OverDot[θ]))))=((∂L)/(∂θ)).

This yields the following.

(Subscript[cm, 1])=D[L,ϕ'[t]]//Simplify
(Subscript[cm, 2])=D[L,ψ'[t]]//Simplify
eq=D[D[L,θ'[t]],t]-D[L,θ[t]]//Simplify

The first two expressions are constants of motion, while the last one must be equal to zero.

Introducing a new variable Φ by

(OverDot[ϕ])=((J(OverDot[Φ])+(J-(Subscript[J, 3]))(OverDot[ψ]) cos θ)/(Subscript[J, 3])),

the last three expressions simplify further. It is also important to realize that the transformation () replaces only the ϕ parameter (the speed of spinning); if we are interested in the behavior of ψ and θ only (as is usually the case), we get the same solution regardless of whether we use ϕ (in the old set of equations) or Φ (in the new one).

So, let us see what the new results look like.

(Subscript[cm, 1])=(Subscript[cm, 1])/J/.ϕ'[t]->(J Φ'[t]+(J-(Subscript[J, 3]))ψ'[t]Cos[θ[t]])/(Subscript[J, 3])//Simplify
(Subscript[cm, 2])=(Subscript[cm, 2])/J/.ϕ'[t]->(J Φ'[t]+(J-(Subscript[J, 3]))ψ'[t]Cos[θ[t]])/(Subscript[J, 3])//Simplify
eq=eq/J/.ϕ'[t]->(J Φ'[t]+(J-(Subscript[J, 3]))ψ'[t]Cos[θ[t]])/(Subscript[J, 3])//Simplify

The equations are now much simpler, and the number of parameters has been reduced from three to only one (namely m, whose dimension is (sec^(-2))).

We can now easily infer a possibility of what we call a “steady” solution, with a constant value of θ (say (Subscript[θ, 0])) and

Φ=(Subscript[(OverDot[Φ]), 0])t+(Subscript[Φ, 0]),
ψ=(Subscript[(OverDot[ψ]), 0])t+(Subscript[ψ, 0]),

where (Subscript[(OverDot[Φ]), 0]) and (Subscript[(OverDot[ψ]), 0]) are two constants whose product must be equal to m. To relate the precession speed (Subscript[(OverDot[ψ]), 0]) to the actual spinning speed, say (Subscript[(OverDot[ϕ]), 0]), we have to solve (based on ())

(Subscript[J, 3])(OverDot[(Subscript[ϕ, 0])])=((m J)/(Subscript[(OverDot[ψ]), 0]))+(J-(Subscript[J, 3]))(Subscript[(OverDot[ψ]), 0])cos (Subscript[θ, 0])

for (OverDot[(Subscript[ψ, 0])]), which yields

(OverDot[(Subscript[ψ, 0])])=(((Subscript[J, 3])(OverDot[(Subscript[ϕ, 0])])±(Sqrt[((((Subscript[J, 3])(OverDot[(Subscript[ϕ, 0])])))^2)-4m J (J-(Subscript[J, 3]))cos (Subscript[θ, 0])]))/(2(J-(Subscript[J, 3]))cos (Subscript[θ, 0]))),

where the negative sign corresponds to the usually observed “slow” precession (assuming, from now on, that (OverDot[(Subscript[ϕ, 0])])>0), whereas the positive sign yields a hard-to-achieve “fast” precession. When J=(Subscript[J, 3]) (spherically symmetrical top), we get only the “slow” solution (OverDot[(Subscript[ψ, 0])])=(m/(OverDot[(Subscript[ϕ, 0])])). For a fast-spinning top, the previous formula yields, to a good approximation, (OverDot[(Subscript[ψ, 0])])=((J m)/((Subscript[J, 3])(OverDot[(Subscript[ϕ, 0])]))); the precession speed is practically independent of (Subscript[θ, 0]).

Approximate Solution

The actual motion is usually more complicated than a steady solution of the previous section, as θ may undergo periodic oscillations, called nutation. Assuming that the amplitude of these is small, we can expand the three dependent variables of our equations around the steady-state solution:

θ=(Subscript[θ, s])+ϵ (Subscript[θ, 1]),
(OverDot[ψ])=(Subscript[(OverDot[ψ]), s])+ϵ (Subscript[(OverDot[ψ]), 1]),
(OverDot[Φ])=(m/(Subscript[(OverDot[ψ]), s]))+ϵ (Subscript[(OverDot[Φ]), 1]),

where the first term of each right-hand side is a constant, the second term remains time dependent, and ϵ is small. We substitute these into the three equations, neglecting (ϵ^2) and higher powers of ϵ.

aux={(Subscript[cm, 1]),(Subscript[cm, 2]),eq}/.{θ[t]->(Subscript[θ, s])+ϵ (Subscript[θ, 1])[t],θ''[t]->ϵ (Subscript[θ, 1])''[t],ψ'[t]->(Subscript[(OverDot[ψ]), s])+ϵ (Subscript[ψ, 1])'[t],Φ'[t]->m/(Subscript[(OverDot[ψ]), s])+ϵ (Subscript[Φ, 1])'[t]};
approx=Coefficient[Series[aux,{ϵ,0,1}],ϵ]

Each of the three expressions must be equal to zero.

We solve the first two equations for (Subscript[ψ, 1])' and (Subscript[Φ, 1])' and substitute into the last equation.

aux=Solve[{approx[[1]]==0,approx[[2]]==0},{(Subscript[ψ, 1])'[t],(Subscript[Φ, 1])'[t]}]//Simplify
approx[[3]]/.aux//Simplify

Solving for (Subscript[θ, 1]) (and, subsequently, for (Subscript[ψ, 1]) and (Subscript[Φ, 1])) is now quite easy. Note that the coefficient of (Subscript[θ, 1]) in the last expression is always non-negative.

Also note that (OverDot[ψ]) will be always positive (assuming “slow” precession), whenever the following condition is met:

(( (m-cos (Subscript[θ, s]) (Subscript[(OverDot[ψ]), s]^2))ϵ )/((Subscript[(OverDot[ψ]), s])sin (Subscript[θ, s])))<(Subscript[(OverDot[ψ]), s]).

Reversing the above inequality results in so-called looping orbits (in terms of the path of the axis, when displayed on a unit sphere: see below), with (OverDot[ψ]) periodically changing direction; making the two sides of () equal to each other results in cuspidal orbits (when θ reaches its smallest value, (OverDot[ψ]) becomes zero).

If desired, one can extend the approximate solution to achieve (ϵ^2) (and higher) accuracy.

Exact Solution

One can show that the original differential equation for θ, namely

(Overscript[θ, ..])+(OverDot[Φ])(OverDot[ψ]) sin θ-m sin θ=0

has the following constant of motion.

(Subscript[cm, 3])=(1/2)(((θ'[t])^2)+((ψ'[t])^2)+((Φ'[t])^2))+ψ'[t]Φ'[t]Cos[θ[t]]+m Cos[θ[t]];

This can be easily verified.

D[(Subscript[cm, 3]),t]/θ'[t]/.Solve[{D[(Subscript[cm, 1]),t]==0,D[(Subscript[cm, 2]),t]==0},{ψ''[t],Φ''[t]}][[1]]//Simplify

Making (Subscript[cm, 1]) and (Subscript[cm, 2]) equal to their initial values, we can solve for (OverDot[ψ]) and (OverDot[Φ]) and substitute these into (Subscript[cm, 3]) (minus its initial value, so that the resulting expression must be equal to zero).

(Subscript[cm, 3])=(Subscript[cm, 3])-(1/2)(((Subscript[(OverDot[θ]), 0])^2)+((Subscript[(OverDot[ψ]), 0])^2)+((Subscript[(OverDot[Φ]), 0])^2))-(Subscript[(OverDot[ψ]), 0])(Subscript[(OverDot[Φ]), 0])Cos[(Subscript[θ, 0])]-m Cos[(Subscript[θ, 0])];
eq=(Subscript[cm, 3])/.(aux=Solve[{(Subscript[cm, 1])==(Subscript[(OverDot[Φ]), 0])+(Subscript[(OverDot[ψ]), 0])Cos[(Subscript[θ, 0])],(Subscript[cm, 2])==(Subscript[(OverDot[Φ]), 0])Cos[(Subscript[θ, 0])]+(Subscript[(OverDot[ψ]), 0])},{Φ'[t],ψ'[t]}][[1]])//Simplify;

The resulting differential equation can be simplified by the θ=2arccos λ transformation [].

eqn=TrigExpand[eq/.(Subscript[θ, 0])->2q]/.{(Subscript[(OverDot[θ]), 0])->0,q->ArcCos[(Subscript[λ, 0])]}//Simplify;
eqn=TrigExpand[eqn/.θ[t]->2q]/.{θ'[t]->-2λ'[t]/(Sqrt[1-((λ[t])^2)]),q->ArcCos[λ[t]]}//Simplify;
eqn=Collect[eqn(-1+((λ[t])^2))/2//Simplify,{λ'[t],λ[t]},Simplify]

To simplify the subsequent solution, we have assumed (without loss of generality) that the initial time has been chosen to coincide with θ reaching its maximum value (and, consequently, (Subscript[(OverDot[λ]), 0])=0). The above equation can then be solved, using specific initial conditions, which must satisfy (Subscript[(OverDot[Φ]), 0])(Subscript[(OverDot[ψ]), 0])≥m (due to (), since (Subscript[(Overscript[θ, ..]), 0])≤0); here (Subscript[λ, sq]) provides the solution for (λ^2).

{eqcf,Λ,aux}={Coefficient[eqn+((λ'[t])^2),λ[t],{4,2,0,-2}],((Subscript[λ, 0])^2),aux}/.Thread[{(Subscript[λ, 0]),(Subscript[(OverDot[ψ]), 0]),(Subscript[(OverDot[Φ]), 0]),m}->{0.85,1,3,1}];
δ=x/.Solve[0==eqcf.{(x^2)+3Λ x+3(Λ^2),x+2Λ,1,0},x][[1]];
μ=(Sqrt[-eqcf.{3Λ+δ,1,0,0}]);
k=δ eqcf[[1]]/(μ^2);
(Subscript[λ, sq])=Λ+δ ((JacobiSN[μ t,k])^2)

This can be easily transformed back to θ; one can then find the corresponding Φ and ψ by analytically integrating the expressions for (OverDot[ψ]) and (if desired) (OverDot[Φ]). Finally, the results can be displayed graphically.

(Subscript[(OverDot[ψ]), sol])=ψ'[t]/.aux/.{Cos[θ[t]]->2((λ[t])^2)-1,Cos[(Subscript[θ, 0])]->2Λ-1}/.{((λ[t])^2)->Λ+δ J,((λ[t])^4)->(((Λ+δ J))^2)};
(Subscript[ψ, sol])=3.7+Apart[(Subscript[(OverDot[ψ]), sol]),J]/.{((c_.)/(a_.+b_. J))->((c EllipticPi[-(b/a),JacobiAmplitude[μ t,k],k] JacobiDN[μ t,k])/(a μ (Sqrt[1-k ((JacobiSN[μ t,k])^2)])))};
(Subscript[θ, sol])=2ArcCos[(Sqrt[(Subscript[λ, sq])])];
Show[ParametricPlot3D[{Cos[u]Sin[v],Sin[u]Sin[v],Cos[v]},{v,0,Pi},{u,0,2Pi}],
ParametricPlot3D[Evaluate[{Cos[(Subscript[ψ, sol])]Sin[(Subscript[θ, sol])],Sin[(Subscript[ψ, sol])]Sin[(Subscript[θ, sol])],Cos[(Subscript[θ, sol])]}],{t,0,10}]]

One can now compute and display the top’s motion under various initial conditions, check the accuracy of the approximate solution, and so forth. We will leave it for you to explore. Have fun!

References

D. Hestenes, New Foundations for Classical Mechanics, 2nd ed., Dordrecht/Boston/London: Kluwer Academic Publishers, 1999.

H. Goldstein, Classical Mechanics, 2nd ed., Reading, MA: Addison-Wesley, 1980.

E. T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, 2nd ed., Cambridge: Cambridge University Press, 1917. www.archive.org/details/atreatiseonanal00whitgoog.