2005], A symmetric splitting method for rigid body dynamics

It has been known since the time of Jacobi that the solution to the free rigid body Euler equations of motion is given in terms of a certain type of elliptic functions. Using the Arithmetic-Geometric mean algorithm, Abramowitz & Stegun (1992), these functions can be calculated efficiently and accurately. The overall approach yields a faster and more accurate numerical solution to the Euler equations compared to standard numerical ODE and symplectic solvers. In this paper we investigate the possibility of extending this approach to the case of rigid bodies subject to external forces. By using a splitting strategy similar to the one proposed in Reich (1996), we decompose the vector field of our problem in a free rigid body (FRB) problem and another completely integrable vector field. We apply the method to the simulation of the heavy top.


Introduction
The study of numerical discretizations of Hamiltonian systems which preserve one or several of the geometric features of the continuous equations, has been quite popular in the last few decades, [7], [8].Symplectic methods are suited for the numerical solution of Hamiltonian problems.The energy is almost conserved along the numerical flow computed by a symplectic method, with essentially no accumulation of errors in time, [7].Moreover these methods show longterm stability of the integration.Symplectic Runge-Kutta methods are designed to preserve the symplectic structure of canonical Hamiltonian problems, q = ∂H(q, p) ∂p , ṗ = − ∂H(q, p) ∂q .
If the Hamiltonian problem considered is a holonomically constrained Hamiltonian system, (i.e. an algebraic constraint of the type g(q) = 0 is also present), symplectic integration cannot be achieved by the straightforward application of symplectic methods, the use of appropriate partitioned Runge-Kutta methods is required, [7].A efficient way to address constrained Hamiltonian problems is, when possible, to split the constrained Hamiltonian vector field in a sum of two or more constrained Hamiltonian vector fields which can be separately integrated.The flows are then composed together to form the final symplectic numerical solution.
A similar situation arises when the Hamiltonian problem is not given in the canonical form (1), [7].This is the case of the Euler equations describing the motion of a free rigid body (FRB), where I 1 , I 2 and I 3 are the principal moments of inertia.These equations are completely integrable.Energy and angular momentum are preserved along the solution, this means that for all times the two quantities are constant, (here E is the energy and G 2 is the total angular momentum).There is also a non canonical symplectic structure preserved by the flow of (2), [8].By using the two constants of motion it is possible to derive the solution of the equations expressed in terms of Jacobi elliptic functions.Consider if b 3 /a 3 ≤ b 1 /a 1 (we will have a similar situation if b 3 /a 3 ≥ b 1 /a 1 ), the solutions of the Euler equations are where the Jacobi elliptic functions cn, sn and dn, are defined by cn u = cos ϕ, sn u = sin ϕ, Here, the amplitude ϕ is given implicitly by where is an elliptic integral of the first kind with modulus k.
We have and τ is a constant of integration that is used to satisfy the initial conditions.The simulation of rigid body motion is interesting for applications in robotics, structural mechanics, and molecular dynamics, [6], [8].The formulae (5) for the exact solution of the free rigid body can be turned into a numerical method by using efficient numerical approximations of the Jacobi elliptic functions.In this paper we show how this approach is very competitive and can be further applied to problems of rigid bodies subject to external forces.We also refer to [3] and [12] for related literature.Symplectic integration methods for the Euler equations have been constructed by various authors, [5], [13], [9], see also [8] and references therein.In spite of being a source of insight, many of these methods cannot be straightforwardly generalized to the broader class of non canonical Hamiltonian problems, thus their use is limited to the numerical approximation of the FRB equations.However some of these numerical tools have succesfully been applied in the simulation of rigid body dynamics.
Recently McLachlan and Reich have, independently, proposed a splitting method for the rigid body equations [10], [14].In this method the right hand side of (2) is split in the following three terms, Each of the three vector fields is Hamiltonian, and defines a differential equation which is easy to integrate exactly.The appropriate composition of the corresponding flows produces a symplectic approximation of the problem.This symplectic method seems to outperform most of the known and previously proposed strategies of symplectic integration of the FRB problem, [4], [11], and for this reason we will use this method as a comparison method in our numerical experiments.Accurate, symplectic, energy and momentum preserving approximations for the solution of the FRB equations, have been recently addressed in [11].Our numerical tests show that such approximations can be easily achieved by using the exact solution and computing (6) to machine accuracy.
To solve problems modelled by rigid bodies subject to external forces, in [14], a symplectic method based on a splitting of the vector field into a FRB problem and another Hamiltonian vector field has been considered.We will briefly recall this approach in the case of the heavy top in section 4. In the original approach of [14] the author uses the splitting method (7) to solve the FRB problem.In this paper we instead use the formulae (5) and we compute the elliptic functions ( 6) by the Arithmetic-Geometric mean algorithm.We obtain a symmetric second order splitting method, each of the FRB problems is computed more accurately, and in many examples we achieve a lower computational cost.Some technical issues for the implementation of this approach are discussed in section 3.In section 4 we report some numerical experiments comparing the proposed approach to the symplectic splitting of [14], [2].

A splitting method for the heavy top
Efficient integrators for the free rigid body can be used in connection with splitting methods in the numerical approximation of more complex rigid body dynamics.
We study the asymmetric heavy top, subject to a gravitational force equal to one.The method presented here can also be applied to problems of interacting rigid bodies, rigid bodies subject to external forces, rigid body linked by constrains, etc [2].
In the case when the center of mass is located one unit from the origin, the energy for the (asymmetric) heavy top can be formulated in terms of the following Hamil-tonian: where π = (I 1 ω 1 , I 2 ω 2 , I 3 ω 3 ) T is the angular momentum and u = (u 1 , u 2 , u 3 ) is the position of the center of mass.
The Hamiltonian H gives rise to a problem which can be formulated in terms of a system of ordinary differential equations where is the inertia tensor.To derive a symmetric splitting method for the above equations, we start by applying a Strang splitting (Störmer/Verlet splitting).The problem is split into the two systems and respectively, and the numerical scheme is where ϕ represent the exact flows of S 1 and S 2 .One can see that S 1 corresponds to the kinetic part of H , and S 2 to the potential part.
The first equation in (10) is independent from the external force u, and is simply a free rigid body problem.Moreover, the system (10) has one invariant, namely π 2 = constant and hence, the solution of (10) can also be written as where R(t) ∈ SO (3).Here the rotation matrix R(t) satisfies the differential equation Ṙ = skew(T −1 π)R, R(t 0 ) = I, (12) and from this it follows that is a solution for the second equation in (10).
Rewriting the first part of system (10) in terms of the angular velocity, ω = T −1 π, one obtains the Euler equations (2).We can compute π(t), for any t and any initial value π ( j) , to machine accuracy by using the exact solution of the Euler equations, and computing the Jacobi elliptic integrals by the method of descending Landen transformations.Hence, the update of π, on the interval [t j ,t j+1 ], is "exact".However, R(h) is not computed explicitly and therefore must be substituted by an approximation in the update of the center of mass, i.e.
Here R( j) (h) is obtained integrating ( 12) by a second order symmetric Magnus method [7], and taking R(t j ) = I as initial value.This results in the following expression where π ( j+1/2) = π(t j + h/2) is obtained as a byproduct of the update π ( j+1) in (14) with little extra cost (see section 3 for details).The exponential in (15) is computed by Rodrigues formula.Thus, the flow ϕ h is approximated by a second order flow φ h can be calculated exactly and the updating (second order) scheme is finally It is easy to verify that φ −h = I and the overall splitting method has the time-symmetry property,

Implementation issues
For each integration step, as described in the previous section, we have to calculate the exact solution π ( j+1) := π(t j+1 ) and π ( j+1/2) := π(t j+1/2 ) of a free rigid body problem with the initial condition π(t j ) = π ( j) , computed at the previous time step.For this purpose we use the MATLAB standard function ellipj, which returns the values of sn u = sin ϕ, cn u = cos ϕ, dn u = 1 − k 2 sin 2 ϕ, for a given input u = λ(t − τ), with t = (t j+1/2 ,t j+1 ) T and m = k 2 .Here λ and k are given, except for the sign (se section 1), and τ must be determined appropriately to ensure that the initial condition ω(t j ) = ω ( j) is satisfied.To this end, we first find the amplitude ϕ ( j) ∈ [0, 2π], which is uniquely determined from the equations Furthermore, from the sign of ω 3 (t j ), we determine the sign of the constants a 3 and λ.Now we compute τ from the equation where .
The latter integral can be computed to machine accuracy using the descending Landen transformations.We consider the sequence {ϕ n , defined by tan(ϕ where a n , b n are given by the Arithmetic-Geometric Mean series, [1], i.e.
One can show that 2 n a n .
As the Geometric-Arithmetic Mean series converges quadratically, one obtains very accurate approximations of F(ϕ ( j) |k 2 ) in few recursion steps.

Numerical experiments
The splitting method proposed in this paper is compared with the symplectic method of [2] and [14] which we denote in short by MR, and with the classical fourth order Runge-Kutta method (RK4).We also refer to the second order symmetric method, described the previous section, as SEJ.The symplectic method MR is based on a splitting of the the Hamiltonian H into four pieces.
Each of the corresponding Hamiltonian vector fields can be integrated exactly ( H1 , H2 , H3 correspond to the vector fields ( 7)), the symmetric composition of the flows gives rise to the approximation scheme, where Here is the contribution from the kinetic parts, H1 , H2 and H3 .The flows of kinetic parts correspond to elementary rotations in R 3 .E.g. where While the flow for H4 is the same as for the system S 2 (9) of the previous section, i.e. ϕ 4,h = ϕ h .
In the first experiment we consider the integration of the free rigid body without any external forces.In figure 1 we plot on the x axis the number of floating point operations required by the methods to perform the integration on the interval [0, 1], for different choices of the step size.On the y axis we report the corresponding values of the 2-norm of the global error.In all the experiments the reference solution for computing the global error is obtained using the built in function of MATLAB, ode45, setting the absolute and relative tolerance equal to 10e − 14.We compare the MR method, which in this case involves the computation of the three flows corresponding to the Hamiltonians H1 , H2 , H3 only, with the integration performed using the Jacobi elliptic functions.The latter method produces a very accurate solution of the problem (the error is of the size of 10 −14 ) and the error is independent on the step size of integration.With about the same amount of floating point operations the MR method produces a much pourer approximation of the solution, in this case the accuracy depends on the step size and increases with the number of floating point operations.
In this experiment the principal moments of inertial and the initial value for ω are I 1 = 1, I 2 = 2, I 3 = 3 and In the second experiment we perform the numerical integration of the same free rigid body problem as before on the interval [0, 400].In figure 2 we consider the energy error as the difference between the constant exact energy, given by H , and the energy obtained from the numerical methods (with step size h = 0.4).We note that for the Runge-Kutta method there is a visible energy drift, for the MR method we observe a typical behavior of symplectic methods, i.e. the energy is oscillating around a constant value, (the amplitude of the oscillations is about 10 −3 ), the method based on the use of the Jacobi elliptic functions, SEJ, computes a very accurate solution of the free rigid body problem and the energy is conserved to the same accuracy (the energy error is about 10 −14 ).We conclude that computing the analytic solution of the FRB using ellipj of MATLAB, for computing the Jacobi elliptic functions, is a very efficient and accurate numerical approach for FRB.
In figure 3 we report the results of the third experiment.We consider the integration of the heavy top problem with the splitting methods MR and SEJ on the interval [0, 1], where the principal moments of inertia, the initial value for ω and u are I 1 = 1, I 2 = 1.1,I 3 = 2.3, ω 0 = (2, 4, 1) T and u 0 = (0, 0, 1) T , respectively.The small difference between I 1 and I 2 , I 2 −I 1 = 0.1, means that the corresponding rigid body is almost symmetric.Also in this case we plot the number of floating point operations versus the global error for the two splitting methods.The performance of the two methods is quite similar with a slight advantage for SEJ.In figure 4 we repeat the experiment integrating on the interval [0, 400] with step size h = 0.4.Also in this case we report the energy error for the three numerical methods.The energy error has an oscillating behavior for both MR and SEJ.In the next two experiments we consider the integration of a heavy top with principal moments of inertia I 1 = 0.1, I 2 = 1, I 3 = 10, and the initial value for ω, and u taken as ω 0 = (3, 3, 3) T and u 0 = (0, 0, 1) T .We first integrate on the interval [0, 1] and compare the performance of the two splitting methods in terms of floating point operations against global error, figure 5, in this case the advantage of the new splitting method is quite clear.Next we illustrate the qualitative performance of the two methods by comparing the results obtained by using different step sizes.We look at the energy error and the numerical trajectory describing the motion of center of mass, in figures 6 and 7.For step size h = 0.01 the two methods produce similar trajectories (the norm  u SEJ − u MR = O(10 −1 )).The energy error for SEJ is a factor 10 −3 smaller than for MR.Increasing the step size to h = 0.02 and h = 0.04, the amplitude of the oscillations in the energy error increases for both methods.Consistently for all the experiments, SEJ has much smaller energy error than MR.For h = 0.02 the trajectory of the center of mass produced by the MR method is quite different from the one produced with step size h = 0.01, and it becomes completely chaotic for h = 0.04.For the SEJ method the numerical trajectory of the center of mass preserve the same qualitative behavior for h = 0.01, h = 0.02 and h = 0.04, figure 6.Finally, we compare the energy error between the two methods for different values of the inertia tensor.The initial conditions for the angular velocity and the position of the center of mass are chosen ω 0 = (3, 3, 3) and u 0 = (0, 0, 1).The energy error is measured when I 1 and I 2 are held fixed at I 1 = 0.1, I 2 = 0.2 and I 1 = 0.1, I 2 = 5, respectively, while I 3 is varying between I 2 < I 3 < 10.Studying the two graphs in figure 8    appears that the energy error for the SEJ method is generally lower than for the MR method, in particular when the inertia is significantly larger in one direction (I 3 ) and the body is strongly asymmetric.

Conclusions
In this paper we presented a symmetric splitting method for the integration of rigid body problems subject to external forces.The numerical strategy is based on the use of available efficient algorithms for the computation of Jacobi elliptic functions.We show that thanks to its symmetry, the method performs much better then classical integrators, like the Runge-Kutta method of order 4, in the long time integration of the considered problems.We also compared the method with a similar symplectic splitting method which inspired the present approach.In many of the performed experiments the presented symmetric splitting is more efficient then the symplectic splitting, giving smaller global error for the same amount of floating point operations.Moreover the new method presents in many experiments a better energy conservation.This seems to be true especially for problems where the principal moments of inertia are very different in size.

Figure 1 :
Figure 1: Free rigid body.Number of floating point operations against the global error.Integration on the interval [0, 1] with different step sizes.

Figure 2 :
Figure 2: Energy for the free rigid body.Integration on the interval [0, 400].

Figure 3 :
Figure 3: Heavy top.Number of floating point operations against the global error.Integration on the interval [0, 1] with different step sizes.

Figure 5 :
Figure 5: Heavy top.Number of floating point operations against the global error.Integration on the interval [0, 1] with different step sizes.

Figure 6 :
Figure 6: Energy error for the heavy top.Integration on the interval [0, 100], using the method SEJ.
Center of mass u, h = 0.02, MR