Discrete Variational Symmetrical Top: an Archetype for Rigid-Body Dynamics

Brian Beckman, December 2024, Revised January 2025

Abstract

We compare two numerical solutions for the free motion of a spinning, symmetrical top. One solution is by the standard method of integrating differential equations, discretizing the independent time variable as late as possible in the solution process. The other solution begins by discretizing the dynamical variables and the Lagrangian as early as possible in the solution process, then solving derived algebraic recurrences by root-finding.
Our contribution is to find that the early-discretized solution has no secular, long-term drift in conserved quantities (energy and momenta), unlike the standard, late-discretized solution, which exhibits small but steady drifts in the sixth or seventh significant figure.
While the early-discretized solution has no secular drifts, its conserved quantities have oscillatory noise in the fourth or fifth significant figure, about 2000 times larger in amplitude than the drift in the standard solution. Considering its other benefits—ease of implementation and tuning—we find that the discrete method is worth further study and optimization.

Introduction

Traditionally, one solves problems in classical mechanics via differential equations of motion, discretized for numerical solution as last as possible in the solution process. We call this tactic late discretization. By contrast, one can discretize the dynamical variables and the Lagrangian at the outset and avoid differential equations and their numerical methods altogether. We call this tactic early discretization. The resulting equations of motion are algebraic, solved by iterated root-finding. Theoretically, at least, the early-discretized methods conserve energy and momenta by construction.[1, 2, 3]
Numerical methods for differential equations comprise a vast variety of algorithms, sensitive to multiple tuning parameters. Bringing a solution to ground can be an expensive process of trial-and-error and compromises. By contrast, discrete methods have only one tuning parameter—the time step.
In this notebook, we contrast a traditional numerical solution for the symmetrical top with a new discrete solution. The two solutions match to three significant figures: within the width of a plotting line. The traditional solution exhibits secular drift in conserved energy and momenta. While small—in the sixth or seventh significant figure—the poor conservation will eventually become noticeable, requiring a “reboot.”. By contrast, the discrete solution exhibits only stable oscillatory error in the conserved quantities. While larger—in the fourth or fifth decimal place—these errors are not noticeable, mostly because they are stable in the long run.
The contribution of this notebook is applying discrete dynamics to this problem from scratch. While we have seen theoretical discrete dynamics for the top,[9] we have not seen a convincing numerical demonstration. The theory leaves us unsatisfied that the method can be effectively implemented, let alone easily. Here, we show that it can be both.
Despite these demonstrable benefits, the discrete implementation in this notebook is too slow. However, with this and our prior successes with discrete harmonic oscillators,[7, 8] we are convinced that the natural benefits of discrete dynamics justify further study and optimization.

References

1
.
Marsden and West, 2001, Discrete Mechanics and Variational Integrators, http://www.cds.caltech.edu/~marsden/bib/2001/09-MaWe2001/MaWe2001.pdf
2
.
Matthew West, Variational Integrators, Caltech Thesis, May 2004, https://thesis.library.caltech.edu/2492/1/west_thesis.pdf
3
.
Marin Kobilarov and Jerrold E. Marsden, Discrete Geometric Optimal Control on Lie Groups, IEEE Transactions on Robotics, Vol. 27, No. 4, August 2011.
4
.
https://hepweb.ucsd.edu/ph110b/110b_notes/node36.html
5
.
https://bingweb.binghamton.edu/~suzuki/Math-Physics/LN-13S_Symmetric_top.pdf (vetted against Reference 4 in a separate notebook)
6
.
William L. Burke, Applied Differential Geometry, Cambridge University Press; 1st edition (January 12, 2008)
7
.
Brian C. Beckman, Harmonic Oscillator via Discrete Lagrangian, https://community.wolfram.com/groups/-/m/t/3161109, 2024
8
.
Brian C. Beckman, Discrete Damped Harmonic Oscillator, https://community.wolfram.com/groups/-/m/t/3174169, 2024
9
.
A. I. Bobenko, Yu. B. Suris, Discrete time Lagrangian mechanics on Lie groups, with an application to the Lagrange top, https://arxiv.org/abs/solv-int/9810018v1, 1998

Outline

1
.
Set up the problem.
2
.
Numerically solve the continuous Euler-Lagrange equations, vetting Reference 4.
3
.
Explore constants of the motion. Observe numerical drift in energy and momenta.
4
.
Formulate and solve the discrete Lagrangian. Compare results with the continuous solution. Observe only stable oscillatory noise.

Problem Set-Up

There are four parameters or free variables for the problem. Because the top is symmetrical about its spinning axis, its moment-of-inertia tensor has only two independent components, free variables
I
1
and
I
3
. Alternatively, we can model the top as an oblate spheroid with semi-major axis
l
and equatorial radius
s
as free variables.
The acceleration of gravity
g
and the mass
m
of the top are the other two free variables.

Continuous Equations of Motion

In[]:=
ClearAll[L,F,V,T,q,ω,ω2,ϕ,θ,ψ,ϕd,θd,ψd,I1,I3,l];

Lagrangian

Generalized Coordinates and Velocities

Choose generalized coordinates azimuth,
ϕ
, co-elevation, θ, and nodal angle about the spin axis, ψ, i.e., Euler angles in the
x
-convention (https://mathworld.wolfram.com/EulerAngles.html).
The corresponding generalized velocities are
•
ϕ
,
•
θ
, and
•
ψ
, written as ϕd, θd, ψd because
•
ϕ
,
•
θ
, and
•
ψ
are not valid Mathematica syntax for variables.
◼
  • The standard notation
    ∂L∂
    •
    ϕ
    is a cruel pun. Consider the dedication of Reference 6, “To all those who, like me, have wondered how in hell you can change
    •
    q
    without changing
    q
    .”
    • ◼
  • We can’t take the derivative of
    L
    with respect to the time derivative
    •
    ϕ
    =ϕ'[t]
    of a function
    ϕ[t]
    .
    ∂L∂
    •
    ϕ
    really means
    ∂L/∂(ϕd)
    , where
    ϕd
    is an ordinary parameter variable of the arity-6 Lagrangian function
    L[ϕ,θ,ψ,ϕd,θd,ψd]
    .
    • ◼
  • Reference 6 develops the differential-geometry background, from a physicist’s perspective, necessary for the theoretical developments in References 1, 2, 3, and 9. That’s all we have to say about theory in this notebook.

    • Potential Energy

    By inspection of the height
    l
    of the CG (Center of Gravity) from the contact point:
    In[]:=
    ClearAll[V];​​V[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=mlgCos[θ];​​With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},​​V[X]]
    Out[]=
    glmCos[θ]

    Checked

    Units of measure:
    mass*
    2
    length
    
    2
    time
    =energy
    .
    Intuition check: when θ is small, the top is upright and potential energy is high, just like cosine.
    The force on the top is
    In[]:=
    With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},​​-D[V[X],θ]]
    Out[]=
    glmSin[θ]
    zero at the pole, then acting in the direction of increasing θ, that is, vertically downwards, and larger with larger θ.

    Kinetic Energy

    Consider an inertial reference frame instantaneously oriented with the current Euler angles of the top, but not rotating with the top. Measure the overall angular-velocity vector of the top in this frame in terms of the time derivatives of the angles. Write kinetic energy in terms of that angular-velocity vector in that frame. Denote the 6-tuple
    {ϕ,θ,ψ,ϕd,θd,ψd}
    as the state vector
    X
    . The
    [[1,1]]
    notation in the kinetic-energy formula merely fetches the scalar value from the one-and-only slot in the matrix formula.
    In[]:=
    ClearAll[M,L,T];​​M=
    I1
    0
    0
    0
    I1
    0
    0
    0
    I3
    ;​​ω[{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:={​​{ϕdSin[θ]Sin[ψ]+θdCos[ψ]},​​{ϕdSin[θ]Cos[ψ]-θdSin[ψ]},​​{ϕdCos[θ]+ψd}};​​T[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=
    1
    2
    ω[X].M.ω[X][[1,1]];​​L[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=T[X]-V[X];

    Checked

    Check T against Reference 4 (Mathematica does not automatically reduce
    T
    to the right-hand side, and note that === does not work for unknown reasons):
    In[]:=
    With{X={ϕ,θ,ψ,ϕd,θd,ψd}},​​T[X]==
    1
    2
    I1(
    2
    θd
    +
    2
    Sin[θ]
    2
    ϕd
    )+
    1
    2
    I3
    2
    (Cos[θ]ϕd+ψd)
    //FullSimplify
    Out[]=
    True
    Check our
    L
    against
    L
    hand-copied from Reference 4.
    In[]:=
    1
    2
    I1(
    2
    ϕd
    2
    Sin[θ]
    +
    2
    θd
    )+
    1
    2
    I3
    2
    (ϕdCos[θ]+ψd)
    -mlgCos[θ]​​==With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},L[X]]//FullSimplify
    Out[]=
    True

    Euler-Lagrange Equations

    Writing
    ∂L/∂v
    instead of the pun
    ∂L∂
    •
    q
    , expand
    
    t
    ∂L
    ∂v
    -
    ∂L
    ∂q
    =0
    (
    1
    )
    as
    v
    ranges over
    {ϕd,θd,ψd}
    and
    q
    ranges in lock-step over
    {ϕ,θ,ψ}
    .
    Here is
    ∂L
    ∂q
    as a column vector. Note the zeros.
    In[]:=
    ClearAll[dLdq];​​(dLdq=(With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},​​Table[D[L[X],q],{q,{ϕ,θ,ψ}}]]//FullSimplify))//MatrixForm
    Out[]//MatrixForm=
    0
    (glm-I3ϕdψd+(I1-I3)
    2
    ϕd
    Cos[θ])Sin[θ]
    0
    See dimensions of energy—
    mass*
    2
    length
    
    2
    time
    —because
    ϕ
    ,
    θ
    , and ψ, are dimensionless.
    Here is
    ∂L
    ∂v
    :
    In[]:=
    ClearAll[dLdqDot];​​(dLdqDot=(With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},​​Table[D[L[X],qd],{qd,{ϕd,θd,ψd}}]]//FullSimplify))//MatrixForm

    Conjugate Momenta

    Reference 5 exhibits alternative, first-order equations of motion via these constants. We have, in a separate notebook, shown that References 4 and 5 are mutually consistent.

    Checked

    Manually check against Reference 4:

    Conjugate Momenta as Functions of State

    Second-Order Equations of Motion

    Total time derivatives of the conjugate momenta yield second-order terms:
    Set up symbolic rewrites of the equations of motion for numerical solution via NDSolve:
    It’s not easy to interpret, or even believe, these equations straight-up, but the numerical solutions deliver the goods.

    Numerical Solution

    Pick numerical values for the free variables. The “4” reminds us we’re working with Reference 4. In a later notebook, we reconcile these values with Reference 5.
    Change tN if you want a longer integration.
    Mathematica has no trouble, numerically, with these equations, given plausible initial conditions. There is a great deal of automation hidden in this call: automatic (and opaque) choice of algorithm, time step, and values for various tuning parameters. If we were writing an integrator from scratch, we would have much work to do.

    Conservation of Energy

    Funny Plotting Business

    We’re now equipped to compute energy and momenta.

    Checked

    Energy, about 14 Joule, is conserved to five figures. It is, however, trending upwards due to numerical effects, which we will not see in the discrete solution.

    Conservation of Momenta

    Interactive Exploration

    NDSolve is so fast that we can put it in a Manipulate to play around with initial conditions and integration time.
    Graphics helpers:
    A rendering function:

    Discrete Solution

    Discrete Lagrangian

    Recall the Lagrangian from above.
    The third, unnumbered equation on page 643 of Reference 3 has a useful definition of the discrete Lagrangian, with h as time step.
    It is beautifully simple and easy to memorize, though it has units of action rather than of energy.

    Discrete Euler-Lagrange Equations

    Reference 2 gives us the forced discrete Euler-Lagrange Equations.

    Results and Conclusion

    The discrete solution does a qualitatively good job with a coarse time step of 0.001 seconds, quite large in the scheme of things, and presumably much larger than the time step Mathematica chooses for the continuous solution internally to NDsolve.
    The following plot shows angular differences between the continuous and discrete solutions for the first 14 seconds at a time step of 0.001 seconds for the discrete solution. The angles show both secular-trending and periodic differences in the sixth significant figure. The plot does not reveal the origin of the secular differences: are they from the discrete or the continuous solution, or a combination? This question deserves further work.

    Conservation of Energy

    Examine the energy of the first 14 seconds. Here, we compare the energy fluctuations, each scaled to its standard deviation. Notice no discernible secular energy growth in the discrete solution, despite the fact that the standard deviation of the discrete solution is 2300 times greater.
    The energy error of the discrete solution, while oscillating in the fourth figure, does not drift, despite the coarse time step of 0.001 seconds. In perspective, the initial condition for spin rate is 5 radians per second (~48 RPM), so a hundredth of a second is a little less than a hundredth of a revolution, or around 3 or 4 degrees, maybe just enough time to track precession and nutation. The fact that the discrete solution tracks energy in the mean without drift with such a coarse time step is remarkable, and is one of the great selling points of discrete dynamics.

    Conservation of Momenta

    First, compare the drift of the continuously derived azimuthal momentum versus the mere oscillation, albeit 2300 times larger, of the error in the discretely derived azimuthal momentum.
    The spin angular momentum shows a similar lack of secular drift in the discrete solution, scaled by its standard deviation, which is 1000 times larger than the standard deviation of the continuous solution.