Abstract

In Discrete variational spherical pendulum: free motion,
[31]
​
we showed that Runge-Kutta integration of quaternion dynamics for the spherical pendulum has unacceptable energy divergence. We thus pursued two alternative numerical approaches: (1) late discretization: continuous equations of motion and Wolfram numerical integrators; and (2) early discretization of the Lagrangian, solution by algebraic root finding—discrete variational integration or DVI.
In that reference, we found that Wolfram’s integrators also have secular—long-term with respect to the time step—energy non-conservation in the sixth significant digit. DVI has no secular energy drift, and that’s the whole point of DVI—it honors the intrinsic abstract geometry of the dynamics! But DVI we found periodic energy error, in the third significant digit. This error is imperceptible in a 20-second animation, so we closed that article declaring the trade-off for energy conservation between small, secular and larger, periodic.
By extending the integration time to 200 seconds, we now find that the Wolfram numerical integrators perform better than expected in energy conservation and that discrete variational integration performs worse than expected in energy conservation. The new “pros” and “cons” are the following:

Late Discretization

con: Wolfram’s numerical integrators suffer secular energy drift in the sixth decimal place.
pro: Wolfram’s numerical integrators occasionally correct the energy drift, keeping it in the sixth decimal place. This sophistication is a surprising, new observation.
pro: Wolfram’s numerical integrators are fast enough for interactive use.

Early Discretization

pro: DVI has no secular energy drift, as expected from its mathematics.
con: The periodic energy noise of DVI not only grows secularly, but accelerates! This is a new observation.
con: DVI, as implemented here, is far too slow for interactive use.
Early discretization and DVI remain of high interest because of their intrinsic secular energy conservation and mathematical beauty, but there is a lot of work to do to make them as fast as standard methods and to understand and mitigate the periodic noise.

Coordinate System

What are good generalized coordinates and velocities for The spherical pendulum? Spherical coordinates have a singularity at the North pole, right where we want (in future notebooks) to control the baton. A better choice are pitch and cone roll (not axis roll!): co-latitude δ and angular displacement η at a clockwise right angle to the great-circle arc of co-latitude. These are spherical coordinate tipped over so the singularities are on the Equator, better than singularities at the North poles because the simulated baton—and its numerical integrators—will, with high likelihood, fly through these singularities, not linger at them.

Global Constants (Initialization Cells)

Rectangular basis vectors,
e
1
,
e
2
,
e
3
, origin
o
, acceleration of Earth’s gravity
g
in
m/
2
sec
, and a default plot range that leaves some padding.
In[]:=
ClearAll[e1,e2,e3,o,g,plotRg];​​e1={1,0,0};e2={0,1,0};e3={0,0,1};​​o={0,0,0};g=-9.81;plotRg={-1.25,1.25};

Dynamical Set-Up

Center of mass (CM)

... in terms of a symbolic constant,
c
z
, the height of the CM in the body frame of the baton:
In[]:=
ClearAll[cm];​​(cm[{δ_,η_}]:=RotationMatrix[-η,e1].RotationMatrix[δ,e2].{0,0,
c
z
});​​cm[{δ,η}]//MatrixForm
Out[]//MatrixForm=
Sin[δ]
c
z
Cos[δ]Sin[η]
c
z
Cos[δ]Cos[η]
c
z

z-Height

The potential energy depends only on the z-height,
zh[δ,η]
, the third component of the CM as the baton rotates around the fixed spherical hinge point:
In[]:=
ClearAll[zh];​​zh[{δ_,η_}]:=-cm[{δ,η}][[3]];​​zh[{δ,η}]
Out[]=
-Cos[δ]Cos[η]
c
z

Symbolical and Numerical Forms

We’ll need symbolical and numerical functions for both potential and kinetic energy. The symbolical functions support Wolfram’s calculus with symbolical functions like
δ[t]
,
η[t]
,
δ'[t]
, and
η'[t]
, yielding the Euler-Lagrange equations of motion, and also Wolfram’s numerical integrators
[15]
.
The numerical energy functions support graphical demonstrations.

numRules

For numerical functions, we replace parametric functions
δ[t],η[t],δ'[t],η'[t]
, with unbound symbols
δ,Dδ,η,Dη
. We also add numerical rule for the mass
m
and CM z-height:
In[]:=
ClearAll[numRules];​​numRules={δ[t]->δ,η[t]->η,δ'[t]->Dδ,η'[t]->Dη,m->1,
c
z
->1.2};

Potential Energy

Here are the symbolical and numerical potential-energy functions, with the acceleration of gravity,
g
, numerical in both cases because we treat it as a global constant:
In[]:=
ClearAll[Vnum,V];​​V[{δ_,η_}]:=mgzh[{δ,η}];V[{δ,η}]​​(*Noticenot:=,SetDelayed,but=,Set*)​​Vnum[{δ_,η_}]=((V[{δ,η}])/.numRules);​​Vnum[{δ[t],η[t]}]
Out[]=
9.81mCos[δ]Cos[η]
c
z
Out[]=
11.772Cos[δ[t]]Cos[η[t]]

Kinetic Energy

Velocity Vector

The kinetic energy depends on the velocity vector, computed via the Wolfram D function:
In[]:=
ClearAll[v];​​(v[{δ_,η_}]:=D[cm[{δ,η}],t]);​​v[{δ[t],η[t]}]//FullSimplify//MatrixForm
Out[]//MatrixForm=
Cos[δ[t]]
c
z
′
δ
[t]
c
z
(-Sin[δ[t]]Sin[η[t]]
′
δ
[t]+Cos[δ[t]]Cos[η[t]]
′
η
[t])
-
c
z
(Cos[η[t]]Sin[δ[t]]
′
δ
[t]+Cos[δ[t]]Sin[η[t]]
′
η
[t])
Here are the symbolical and numerical kinetic-energy functions:
In[]:=
ClearAll[Tnum,T];​​T[{δ_,η_}]:=
1
2
mv[{δ,η}].v[{δ,η}];​​T[{δ[t],η[t]}]//FullSimplify​​(*hand-constructed,forchecking*)​​Tnum[{δ_,η_},{Dδ_,Dη_}]:=
1
2
2
(1.2)
(
2
Dδ
+
2
Cos[δ]
2
Dη
);​​Tnum[{δ_,η_},{Dδ_,Dη_}]:=
1
2
2
(1.2)
2
Dδ
+
1
2
(1+Cos[2δ])
2
Dη
;​​(*Set=,notSetDelayed:=*)​​Tnum[{δ_,η_},{Dδ_,Dη_}]=(T[{δ[t],η[t]}]/.numRules//FullSimplify);​​Tnum[{δ,η},{Dδ,Dη}]
Out[]=
1
2
m
2
c
z
(
2
′
δ
[t]
+
2
Cos[δ[t]]
2
′
η
[t]
)
Out[]=
0.72
2
Dδ
+0.36
2
Dη
+0.36
2
Dη
Cos[2δ]

Lagrangians

L is for the late-discretized solution. Lnum is for the early-discretized Lagrangian,
L
d
.
In[]:=
ClearAll[Lnum,L];​​(L[{δ_,η_}]:=T[{δ,η}]-V[{δ,η}]);L[{δ[t],η[t]}]//FullSimplify​​Lnum[{δ_,η_},{Dδ_,Dη_}]:=Tnum[{δ,η},{Dδ,Dη}]-Vnum[{δ,η}];​​Lnum[{δ,η},{Dδ,Dη}]
Out[]=
0.25m
c
z
(-39.24Cos[δ[t]]Cos[η[t]]+
c
z
(2.
2
′
δ
[t]
+(1.+1.Cos[2δ[t]])
2
′
η
[t]
))
Out[]=
0.72
2
Dδ
+0.36
2
Dη
+0.36
2
Dη
Cos[2δ]-11.772Cos[δ]Cos[η]

Leqns

Wolfram has symbolical magic for the Euler-Lagrange equations. See Reference [15].
In[]:=
ClearAll[Leqns];​​With[{Lcrds={δ[t],η[t]}},​​With[{Lvels=D[Lcrds,t]},​​With[{Lncnf={0,0}},(*non-conservativeforces*)​​Leqns=Flatten@MapThread[​​{q,v,f}
D[D[L[Lcrds],v],t]-D[L[Lcrds],q]==f
,​​{Lcrds,Lvels,Lncnf}]]]]//FullSimplify
Out[]=
{m
c
z
(-9.81Cos[η[t]]Sin[δ[t]]+
c
z
(Cos[δ[t]]Sin[δ[t]]
2
′
η
[t]
+
′′
δ
[t]))0,mCos[δ[t]]
c
z
(-9.81Sin[η[t]]+
c
z
(-2.Sin[δ[t]]
′
δ
[t]
′
η
[t]+Cos[δ[t]]
′′
η
[t]))0}

State-Space Form

Differential equations in state-space form are first-order, thus easier for Wolfram to integrate. We get them by solving Leqns for
δ''[t]
and
η''[t]
:
In[]:=
ClearAll[stEqns];​​stEqns=Equal@@@(Solve[Leqns,{δ''[t],η''[t]}]//FullSimplify//Chop//Part[#,1]&)
Out[]=

′′
δ
[t]Sin[δ[t]]
9.81Cos[η[t]]
c
z
-1.Cos[δ[t]]
2
′
η
[t]
,
′′
η
[t]
9.81Sec[δ[t]]Sin[η[t]]
c
z
+2.Tan[δ[t]]
′
δ
[t]
′
η
[t]

Late Discretization: Numerical Solution of Continuous Equations

Solve the system for
tLim=200
simulated seconds (this time limit is the new feature of this notebook over its ancestor notebook, Reference 31); tLim is a new global constant for the rest of this notebook:
In[]:=
ClearAll[tLim,initialConditions,stSolns];​​tLim=200;​​initialConditions={δ[0]==10°,η[0]==5°,δ'[0]==0,η'[0]==0};​​stSolns=NDSolve[Append[stEqns,initialConditions]/.{
c
z
->1.2},{δ[t],η[t],δ'[t],η'[t]},{t,0,tLim}]
Out[]=
δ[t]InterpolatingFunction
Domain: {{0.,200.}}
Output: scalar
[t],η[t]InterpolatingFunction
Domain: {{0.,200.}}
Output: scalar
[t],
′
δ
[t]InterpolatingFunction
Domain: {{0.,200.}}
Output: scalar
[t],
′
η
[t]InterpolatingFunction
Domain: {{0.,200.}}
Output: scalar
[t]

Plot the Solutions

Animation

This is pretty convincing. Energy is well conserved. We’d only have to reboot it after a long time. If we were satisfied with an integrator that works only inside Wolfram, we’d be just about done. But we want an integrator we can code up in any programming language.

Graphical Support Routines

Early Discretization

Discretize the dynamical variables and the Lagrangian, itself, via Equation 7 of Reference [10]. First, a reminder:

Discrete Euler-Lagrange Equations

The discrete Euler-Lagrange equations, without generalized force terms, are
All we need is FindRoot, which we could replace with an explicit (faster?) Newton-Raphson solver in C, C++, Fortran, or Python, for instance. We bootstrap the integrator by taking the first two positions from the late-discretized solution.
Notice that the angles are qualitatively the same as in late discretization, but show a growth in the envelope of the angles. This is unacceptable and undiagnosed.

Energies of Late Discretization

From the late-discretized solution, we calculate a mean energy of mnE = 11.5490. Energy slowly drifts in the fifth decimal place (sixth significant figure). However, we see that Wolfram occasionally corrects the error drift. That is one of the two primary discoveries of this notebook.

Conclusion

The early-discretized solution accumulates and loses energy. Remarkably, Wolfram’s integrators detect energy divergence and correct themselves, though opaquely and internally. That correction, coupled with their remarkable execution speed, makes the Wolfram integrators suitable for secular energy studies, though not for export to internal applications.
We see, again, that early discretization conserves energy in the mean much better than does late discretization. The periodic error of early discretization, though still acceptably small, is 1000 times greater than the secular error drift of late discretization. That periodic error, along with the slow execution speed of early discretization, renders it unpractical. However, it remains mathematically interesting and practically promising if it can be made fast enough and if the periodic noise can be diagnosed and mitigated.

References

This is a list of general references on discrete variational dynamics and the background mathematics. Not all these references are cited directly in this notebook.
1
.
Jack B Kuipers, Quaternions and Rotation Sequences, Princeton University Press, 1999.
2
.
Jeongseok Lee, Karen Liu, Frank C. Park, Siddartha S. Srinivasa, A Linear-Time Variational Integrator for Multibody Systems, 2018, https://arxiv.org/abs/1609.02898
3
.
Ari Stern, Discrete Geometric Mechanics and Variational Integrators, 2006, http://ddg.cs.columbia.edu/SIGGRAPH06/stern-siggraph-talk.pdf.
4
.
Ethan Eade, Lie Groups for 2D and 3D Transformations, 2017, http://ethaneade.com/lie.pdf.
5
.
NIST: Digital Library of Mathematical Functions, https://dlmf.nist.gov.
6
.
Brian C. Hall, Lie Groups, Lie Algebras, and Representations, Second Edition, 2015, Springer.
7
.
Guangyu Liu, Modeling, stabilizing control and trajectory tracking of a spherical inverted pendulum, PhD thesis, 2007, University of Melbourne,
​https://minerva-access.unimelb.edu.au/handle/11343/37225.
8
.
Wikipedia, Tennis-Racket Theorem, https://en.wikipedia.org/wiki/Tennis_racket_theorem.
9
.
Marin Kobilarov, Keenan Crane, Mathieu Desbrun, Lie Group Integrators for Animation and Control of Vehicles, 2009 (?).
10
.
Ari Stern, Mathieu Desbrun, Discrete Geometric Mechanics for Variational Time Integrators, 2006 (http://www.geometry.caltech.edu/pubs/SD06.pdf).
11
.
Kenth Engø, On The BCH Formula in (3), https://www.researchgate.net/profile/Kenth_Engo-Monsen2/publication/233591614_On_the_BCH-formula_in_so3/links/004635199177f69467000000/On-the-BCH-formula-in-so3.pdf
12
.
Alexander Van-Brunt, Max Visser, Explicit Baker-Campbell-Hausdorff formulae for some specific Lie algebras, https://arxiv.org/pdf/1505.04505.pdf.
13
.
Wikipedia, Baker-Campbell-Hausdorff Formula, https://en.wikipedia.org/wiki/Baker%E2%80%93Campbell%E2%80%93Hausdorff_formula
14
.
Jerrold E. Marsden, Tudor S. Ratiu, Introduction to Mechanics and Symmetry, Second Edition, 2002
15
.
Moylan, Andrew, Stabilized Inverted Pendulum, https://blog.wolfram.com/2011/01/19/stabilized-inverted-pendulum/, and Stabilized n-Link Pendulum, https://blog.wolfram.com/2011/03/01/stabilized-n-link-pendulum/
16
.
Paul R. Halmos, Naive Set Theory, 2015.
17
.
David Lovelock and Hanno Rund, Tensors, Differential Forms, and Variational Principles, 1975.
18
.
Ralph Abraham and Jerrold E. Marsden, Foundations of Mechanics: 2nd Edition, 1980.
19
.
Taeyoung Lee, Melvin Leok, N. Harris McClamroch, Discrete Control Systems, https://arxiv.org/abs/0705.3868.
20
.
Taeyoung Lee, Melvin Leok, N. Harris McClamroch, Global Formulations of Lagrangian and Hamiltonian Mechanics on Manifolds, Springer, https://a.co/d/0lTbL9t.
21
.
Tristan Needham, Visual Differential Geometry and Forms, https://a.co/d/9hoKekz.
22
.
xAct, Efficient tensor computer algebra for the Wolfram Language, http://xact.es/index.html.
23
.
Blanco, Jose-Luis, A tutorial on SE(3) transformation parameterizations and on-manifold optimization, Technical Report #012010, Universidad de Malaga (https://jinyongjeong.github.io/Download/SE3/jlblanco2010geometry3d_techrep.pdf).
24
.
Marsden and West has exhaustive treatment of Noether’s Theorem (http://www.cds.caltech.edu/~marsden/bib/2001/09-MaWe2001/MaWe2001.pdf).
25
.
Maxime Tournier, Notes on Lie Groups, https://maxime-tournier.github.io/notes/lie-groups.html.
26
.
John Baez, Javier P. Muniain, Gauge Fields, Knots and Gravity, (https://a.co/d/jdRnOU0).
27
.
F. C. Park, J. E. Bobrow, S. R. Ploen, A Lie Group Formulation of Robot Dynamics.
28
.
George W. Hart, Multidimensional Analysis, https://a.co/d/8c5RqOX
29
.
Evandro Bernardes, Stephane Viollet, Quaternion to Euler angles conversion: A direct, general and computationally efficient method, PLoS ONE 17(11) e0276302
30
.
Linear Quadratic Regulator, from Underactuated Robotics, https://underactuated.mit.edu/lqr.html
31
.
Brian Beckman, Discrete variational spherical pendulum: free motion, https://community.wolfram.com/groups/-/m/t/3359568

CITE THIS NOTEBOOK

Perspectives on energy conservation: discrete variational spherical pendulum​
by Brian Beckman​
Wolfram Community, STAFF PICKS, April 14, 2025
​https://community.wolfram.com/groups/-/m/t/3444021