The symmetrical top is a fundamental problem in classical mechanics, an archetype of rigid-body motion. Numerical methods for its solution are likewise an archetype for simulation and control of machines of all kinds: robots, vehicles, conveyances, etc. While an analytical solution for the symmetrical top exists, it is not of interest here. We are interested in numerical methods that generalize to the infinitely more common cases that can on be solved numerically. In fact, we are interested in solution schemas that can be deploy on real-world problems with little up-front tuning, with acceptable error margins, and with good long-term stability.
The discrete method illustrated in this notebook appears to have these properties. While it introduces oscillatory errors in the fourth decimal place, there is no discernable drift. Contrast to the traditional numerical method, which exhibits long-term instability in the conserved quantities of energy and angular momentum. Furthermore, outside of Mathematica, the traditional solution would be difficult to implement and tune, whereas the discrete solution has only one hyperparameter: the time step.
The discrete method illustrated in this notebook appears to have these properties. While it introduces oscillatory errors in the fourth decimal place, there is no discernable drift. Contrast to the traditional numerical method, which exhibits long-term instability in the conserved quantities of energy and angular momentum. Furthermore, outside of Mathematica, the traditional solution would be difficult to implement and tune, whereas the discrete solution has only one hyperparameter: the time step.
Introduction
Introduction
Traditionally, one solves problems in classical mechanics by writing a continuous Lagrangian function, deriving differential equations of motion from its derivatives, and finally discretizing the equations for numerical solution. We call this late discretization. By contrast, with early discretization, one discretizes the Lagrangian itself at the outset. The resulting equations of motion are algebraic, solved by iterated root-finding.
A nagging problem with late discretization is poor conservation of energy and momenta. In practice, numerical solutions of differential equations have long-term drifts. They might need to be periodically corrected or “rebooted.” By contrast, discrete methods, at least theoretically, conserve energy and momenta by construction, thus being stable in the long run.[1, 2, 3]
Numerical methods for differential equations comprise a vast variety of algorithms and sensitive 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 discrete solution. The two solutions match to three significant figures: within the width of a plotting line. The traditional solution exhibits long-term drift in conserved energy and momenta. While small—in the sixth or seventh decimal place—the poor conservation will eventually become noticeable, requiring a “reboot” to keep numerical errors within required limits. 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 do not diverge.
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. Reading the theory leaves us unsatisfied that the method can be effectively, let alone easily, implemented. Here, we show that it can be.
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 optimizing its implementation.
Outline
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
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 and . Alternatively, we can model the top as an oblate spheroid with semi-major axis and equatorial radius as free variables.
I
1
I
3
l
s
We finesse the issue of consistency between and . We choose numerical values for them and just ride it out. We might consume inconsistencies into a mass adjustment, but such analysis is uninstructive and unnecessary. We reproduce the continuous solution of Reference 4 without worrying over this issue.
l
I
1
The acceleration of gravity and the mass of the top are the other two free variables.
g
m
Continuous Equations of Motion
Continuous Equations of Motion
In[]:=
ClearAll[L,F,V,T,q,ω,ω2,ϕ,θ,ψ,ϕd,θd,ψd,I1,I3,l];
Lagrangian
Lagrangian
Generalized Coordinates and Velocities
Generalized Coordinates and Velocities
Our generalized coordinates are azimuth , co-elevation θ, and nodal angle about the spin axis, ψ, i.e., Euler angles in the -convention (https://mathworld.wolfram.com/EulerAngles.html).
ϕ
x
The corresponding generalized velocities are , , and , written as ϕd, θd, ψd because , , and are not valid Mathematica syntax for variables.
ϕ
θ
ψ
ϕ
θ
ψ
◼
The standard notation is a cruel pun. Consider the dedication of Reference 6, “To all those who, like me, have wondered how in hell you can change without changing .”
∂L∂
ϕ
q
q
◼
We can’t take the derivative of with respect to the time derivative =ϕ'[t] of a function . really means , where is an ordinary parameter variable of the arity-6 Lagrangian function .
L
ϕ
ϕ[t]
∂L∂
ϕ
∂L/∂(ϕd)
ϕd
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
Potential Energy
By inspection of the height of the CG (Center of Gravity) from the contact point:
l
In[]:=
ClearAll[V];V[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=mlgCos[θ];With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},V[X]]
Out[]=
glmCos[θ]
Checked
Checked
Units of measure: .
mass*=energy
2
length
2
time
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[θ]
acting in the direction of increasing θ, that is, vertically downwards, and larger with larger θ.
Kinetic Energy
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 as the state vector . The notation in the kinetic-energy formula merely fetches the scalar value from the one-and-only slot in the matrix formula.
{ϕ,θ,ψ,ϕd,θd,ψd}
X
[[1,1]]
In[]:=
ClearAll[M,L,T];M=
;ω[{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:={{ϕdSin[θ]Sin[ψ]+θdCos[ψ]},{ϕdSin[θ]Cos[ψ]-θdSin[ψ]},{ϕdCos[θ]+ψd}};T[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=ω[X].M.ω[X][[1,1]];L[X:{ϕ_,θ_,ψ_,ϕd_,θd_,ψd_}]:=T[X]-V[X];
I1 | 0 | 0 |
0 | I1 | 0 |
0 | 0 | I3 |
1
2
Checked
Checked
Check T against Reference 4 (Mathematica does not automatically reduce to the right-hand side, and note that === does not work for unknown reasons):
T
In[]:=
With{X={ϕ,θ,ψ,ϕd,θd,ψd}},T[X]==I1(+)+I3//FullSimplify
1
2
2
θd
2
Sin[θ]
2
ϕd
1
2
2
(Cos[θ]ϕd+ψd)
Out[]=
True
In[]:=
1
2
2
ϕd
2
Sin[θ]
2
θd
1
2
2
(ϕdCos[θ]+ψd)
Out[]=
True
Euler-Lagrange Equations
Euler-Lagrange Equations
Writing instead of the pun , expand
∂L/∂v
∂L/∂
q
t
∂L
∂v
∂L
∂q
(
1
)as ranges over and ranges in lock-step over .
v
{ϕd,θd,ψd}
q
{ϕ,θ,ψ}
Here is as a column vector. Note the zeros.
∂L
∂q
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 |
0 |
See dimensions of energy——because , , and ψ, are dimensionless.
mass*
2
length
2
time
ϕ
θ
Here is :
∂L
∂v
In[]:=
ClearAll[dLdqDot];(dLdqDot=(With[{X={ϕ,θ,ψ,ϕd,θd,ψd}},Table[D[L[X],qd],{qd,{ϕd,θd,ψd}}]]//FullSimplify))//MatrixForm
Out[]//MatrixForm=
I3ψdCos[θ]+ 1 2 |
I1θd |
I3(ψd+ϕdCos[θ]) |
Now see dimensions of rather than . The second factor of time in the denominator returns via .
mass*time
2
length
mass*
2
length
2
time
/t
Conjugate Momenta
Conjugate Momenta
For any Lagrangian, (pun intended) are the conjugate momenta . In our case, the first, =, and the third, =, are conserved momenta, or constants of the motion. and are zero, thus so are = and = by the Euler-Lagrange equations; and must be constant because their time derivatives vanish. This is an application of Noether’s Theorem.
∂L
∂
q
p
q
p
ϕ
∂L
∂
ϕ
p
ψ
∂L
∂
ψ
∂L
∂ϕ
∂L
∂ψ
t
∂L
∂
ϕ
p
ϕ
t
t
∂L
∂
ψ
p
ψ
t
p
ϕ
p
ψ
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
Checked
In[]:=
(dLdqDot[[1]]==((I1+I3)ϕd+I3ψdCos[θ]))//FullSimplify
2
Sin[θ]
2
Cos[θ]
Out[]=
True
In[]:=
(dLdqDot[[3]]==(I3(ψd+ϕdCos[θ])))
Conjugate Momenta as Functions of State
Conjugate Momenta as Functions of State
Second-Order Equations of Motion
Second-Order Equations of Motion
Total time derivatives of the conjugate momenta yield second-order terms:
Set up the equations of motion via symbolic rewrites 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
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.
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 various tuning parameters. If we were writing an integrator from scratch, we would have much work to show.
Conservation of Energy
Conservation of Energy
Funny Plotting Business
Funny Plotting Business
We’re now equipped to compute energy and momenta.
Checked
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
Conservation of Momenta
Interactive Exploration
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 Solution
Discrete Lagrangian
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
Discrete Euler-Lagrange Equations
Reference 2 gives us the forced discrete Euler-Lagrange Equations.
Results and Conclusion
Results and Conclusion
The following plot shows angular differences between the continuous and discrete solutions for the first 14 seconds at a time step of 0.01 seconds for the discrete solution. The difference shows a periodic error of amplitude around 0.0025 rad, or about 4.2 arc minutes, and a maximum constant systematic error in the spin angle of around 0.05 rad, or about 1.4 deg. The constant systematic error in the other two angles is much smaller. The systematic error in spin angle is likely due to the bootstrapping procedure and could be mitigated by a more careful initialization of the solution. None of these errors are big enough to be visually discernable in an animation or a line plot of the solutions.
The discrete solution does a qualitatively good job with a coarse time step of 0.01 seconds, quite large in the scheme of things, and presumably much larger than the time step Mathematica chooses internally to NDsolve for the continuous solution. The point of this notebook is not speed of implementation, but validation of both implementation ease and stability of conserved quantities. An optimized implementation could likely be made very fast. The chief suspect is root-finding. We might experiment in the future with a hand-crafted Newton-Raphson method.
Conservation of Energy
Conservation of Energy
Examine the energy of the first 14 seconds. Recall that the continuous solution steadily and unrealistically gains energy.
The energy error of the discrete solution, while oscillating in the fourth figure, does not drift, despite the coarse time step of 0.01 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
Conservation of Momenta
First, compare the drift of the continuously derived azimuthal momentum versus the mere oscillation, albeit much larger, of the error in the discretely derived azimuthal momentum.
Likewise, the continuously derived spin angular momentum drifts while the error in the discretely derived spin momentum merely oscillates in the fourth decimal place.
References
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.
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
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Discrete variational symmetrical top: an archetype for rigid-body dynamics
by Brian Beckman
Wolfram Community, STAFF PICKS, January 1, 2025
https://community.wolfram.com/groups/-/m/t/3348245
by Brian Beckman
Wolfram Community, STAFF PICKS, January 1, 2025
https://community.wolfram.com/groups/-/m/t/3348245