On some aspects of the discretization of the Suslov problem

In this paper we explore the discretization of Euler-Poincar\'e-Suslov equations on $SO(3)$, i.e. of the Suslov problem. We show that the consistency order corresponding to the unreduced and reduced setups, when the discrete reconstruction equation is given by a Cayley retraction map, are related to each other in a nontrivial way. We give precise conditions under which general and variational integrators generate a discrete flow preserving the constraint distribution. We establish general consistency bounds and illustrate the performance of several discretizations by some plots. Moreover, along the lines of [14] we show that any constraints-preserving discretization may be understood as being generated by the exact evolution map of a time-periodic non-autonomous perturbation of the original continuous-time nonholonomic system.


Introduction
The Lagrangian formulation of mechanical systems with nonholonomic constraints has been extensively studied in recent years (see [3,19,26] for a complete description and extensive bibliographies). In short, this kind of systems are characterized by so-called nonholonomic constraints, i.e. constraints involving both configuration as well as velocity variables, and which can not be integrated to purely configuration-dependent constraints (in this case the constraints are called holonomic). Moreover, the preservation of the nonholonomic constraints by suitable integrators is a delicate issue upon which a lot of attention has been put by the Geometric Mechanics comunity (see the recent works, such as [10,14,15,16,22,23,34], which have introduced numerical integrators for nonholonomic systems with very good energy behavior and preservation properties). The approaches in most of these references are based on the ideas of [33,38]. In these works, the continuous variational principles are replaced by discrete ones aiming to obtain proper integrators approximating the continuous dynamics; we will call the integrators related to this framework variational integrators. Analogously, in the case of nonholonomic mechanics, the continuous Lagrange-d'Alembert principle, which provides the actual dynamics, is replaced by a discrete Lagrange-d'Alembert principle on the discrete phase space. Of special interest are the seminal works on nonholonomic integration [10,34], where a discrete version of the Lagrange-d'Alembert principle is proposed by introducing a proper discretization of the nonholonomic distribution.
Reduction theory is one of the fundamental tools in the study of mechanical systems with symmetry and it concerns the removal of symmetries and the associated conservation laws (see [32] for more details). Roughly speaking, we will say that a Lagrangian system is invariant under a Lie group if the Lagrangian function is invariant under the tangent lift of the action of the Lie group on the configuration manifold. Furthermore, for a symmetric mechanical system the process of reduction eliminates the directions along the group variables and thus provides a system with fewer degrees of freedom. When the configuration manifold is the Lie group itself, the case this work is concerned about, the process of reduction for unconstrained Lagrangian mechanics leads from the Euler-Lagrange equations to the Euler-Poincaré ones. Needless to say, nonholonomic systems may also possess symmetries, and the geometrical treatment of such situations was studied in [5,8,9,24] among other references. On the other hand, variational integrators for reduced systems were carefully studied from the theoretical point of view in [28,29]. The combination of these two elements, say symmetric nonholonomic Lagrangian systems and variational integrators, was studied in [13] for a ll system, i.e. a nonholonomic system whose configuration manifold is the Lie group G, and where both the Lagrangian and the constraint distribution are invariant with respect to the induced left action of G on T G. The dynamics of this kind of system is described by the Euler-Poincaré-Suslov equations [25,41,42]. Other interesting discretizations of symmetric nonholonomic systems can be found in [11].
The Euler-Poincaré-Suslov equations on SO(3) are the object of study in this paper, more concretely from the point of view of the results presented in [20]. As already pointed out, the Euler-Poincaré-Suslov equations are obtained by reduction from the nonholonomic equations produced by the Lagrange-d'Alembert principle. This particular reduction process leads to fewer degrees of freedom and also turns the dynamical equations from second-order differential algebraic to firstorder. Moreover, the configuration manifold is a Lie algebra, concretely so(3), a fact which allows the direct application of many results on dynamical systems defined on linear spaces. Particularly, we focus on their discretizations, both direct ones and those coming from discretizations of the Lagrange-d'Alembert principle, and meticulously study the conditions under which the constraints are preserved. In terms of consistency of integrators, we wonder about the relationship between these in unreduced and reduced settings. We find out that this relationship is nontrivial, and can be completely worked out when we choose a particular retraction map relating SO(3) and so (3), say the Cayley map. The study of consistency also requires the introduction of a suitable metric on SO (3). The Suslov problem on SO(3), i.e. a rigid body with some of its body angular velocity components set equal to zero, is treated. General considerations and analysis of this system (and its generalizations) may be found in [4,12,44]. Although we will not consider integrability issues (in that regard we refer to [42,44]), we shall carefully study the consistency order of general and variational discretizations, and show their performance by several simulations (for extensive simulations by other reduced variational integrators we refer to [27]). Finally, we address the discrepancies between continuous and discrete dynamics from a non-autonomous perspective, thereby taking into account previous results in [17,20].
The paper is structured as follows: In §2, §3 and §4 we provide a comprehensive introduction to nonholonomic mechanics, variational nonholonomic integrators, symmetries in nonholonomic Lagrangian systems, the Euler-Poincaré-Suslov equations and, finally, variational integrators for reduced systems. We prove the equivalence of the reaction forces, say Lagrange multipliers, coming from the unreduced and reduced settings, see proposition 4.2. Moreover, we provide a new version of the discrete Euler-Poincaré-Suslov equations, defined on the Lie algebra through the incorporation of a retraction map, in corollary 4.4. §5 establishes the link between the results in [20] and the reduced framework. First, §5.3 accounts for the relationship between the consistency order of the unreduced and reduced setups, which is nontrivial and can be completely worked out when the discrete reconstruction equation is given by a Cayley retraction map. Then, in subsection §5.4, we study general discretizations of the Euler-Poincaré-Suslov problem which preserve the nonholonomic constraints, as well as sufficient conditions in the case of variational integrators guaranteeing such a preservation (proposition 5.2). The section §6 accounts for an exhaustive study of the Suslov problem on SO(3). We give consistency results both for general and variational discretizations, and illustrate their performance through some plots. Finally, §7 is devoted to the study of distribution-preserving discretizations of the Euler-Poincaré-Suslov problem viewed as perturbation of the continuous dynamics. For this, we apply the result by Fiedler and Scheurle [17] to the system under consideration, and establish a bound of the perturbation produced in the unreduced dynamics by a reduced integrator, and conversely.
Throughout the paper, we use Einstein's convention for the summation over repeated indices unless the opposite is stated.

Preliminaries
2.1. The nonholonomic setting. We define a nonholonomic Lagrangian system as a triple (G, L, D), where G is a finite dimensional manifold, L : T G → R is the Lagrangian function and D ⊂ T G is a constraint distribution, which we assume to be a constant rank linear vector subbundle of T G and non-integrable in the Frobenious sense. In this work we are interested in the particular instance G = SO(3); however, the results introduced in the preliminary sections refer to any Lie group, and therefore we stick to the notation G. Locally, the linear constraints are written as follows: where (g i ,ġ i ), i = 1, ..., n, are coordinates of T G and rank (D) = n − m, with m < n. The annihilator D • is locally given by where the one-forms µ α are linearly independent.
In addition to the distribution, we need to specify the dynamical evolution of the system through the Lagrangian function. In nonholonomic mechanics, the procedure leading from the Newtonian point of view to the Lagrangian one is given by the Lagrange-d'Alembert principle. This principle says that a curve g : I ⊂ R → Q describes an admissible motion of the system if δ t2 t1 L (g (t) ,ġ (t)) dt = 0 with respect to all variations such that δg (t) ∈ D g(t) , t 1 ≤ t ≤ t 2 and the fixed endpoint condition is satisfied, and if the velocity of the motion itself satisfies the constraints. It is to be noted that the Lagrange-d'Alembert principle is not variational since we are imposing the constraints on the curve after the extremization. Thus, one may consider the intrinsic data defining the Lagrangian nonholonomic problem to be given by the triple (G, L, D). Using Lagrange-d'Alembert's principle, we arrive at the nonholonomic equations, which in coordinates read where λ α , α = 1, ..., m are "Lagrange multipliers". The right-hand side of equation (2a) represents the reaction forces due to the constraints, and equations (2b) represent the constraints themselves. For more details we refer to [2,3]. As mentioned above, the concern of this paper is the discretization of this kind of systems for G = SO(3) and, moreover, when they are left invariant under the same group. One of the advantages of this group, and all finite-dimensional matrix groups in general, is that they can be embedded in R n , for n big enough. In [20], the discretization of nonholonomic systems in R n was treated, more concretely discretizations that preserve the principal geometrical object, which is the distribution manifold D. Using previous results in [39,40], it was proven the existence and uniqueness of solutions for this problem, even in the non-autonomous case, by reducing the basic equations of motion (2), say the nonholonomic DAE, to a set of ordinary differential equations, say the nonholonomic ODE, on D. Furthermore, upon a main result in [17] it is shown that any integrator that preserves the distribution may be understood as being generated by the exact evolution map of a time-periodic non-autonomous perturbation of the original continuous system. Thus, these results apply in the case of the considered matrix group, although some subtle issues appear when dealing with the group structure as it will be shown. In regards of the discretization, the following definition will be of main importance: Definition 2.1. By a (p, s) order discretization of the nonholonomic problem we understand a sequence of points By | · | we understand the usual Euclidean metric. We note that the notation |g(t k + ) − g k+1 | might be easily misunderstood, since in general G is not a linear space. However, the particular meaning of this metric on SO(3) will be clarified when considering the consistency of reduced systems. On the other hand, the dynamics of (g(t), v g (t), λ(t)) is determined by (2), after introducing the second order conditionġ = v g and projecting the resulting differential algebraic equations onto D, which provides an ODE evolving on D and a unique λ α in terms of the other variables. We note that the evolution prescribed by (2) is local.
2.2. Discrete mechanics and variational integrators. The variational integrators are a kind of geometric integrators for the Euler-Lagrange equations which retain their variational character and also, as a consequence, some of the main geometric properties of the continuous system, such as symplecticity and momentum conservation (see [18,33,38]). In the following we will summarize the main features of this type of geometric integrators. A discrete Lagrangian is a map L d : G × G → R, which may be considered as an approximation of the action integral defined by a continuous Lagrangian L : T G → R, that is where g(t) is a solution of the Euler-Lagrange equations, say joining g(t 0 ) = g 0 and g(t 0 + ) = g 1 for small enough > 0, where will be considered as the step size of the integrator.
Define the action sum S d : G N +1 → R corresponding to the Lagrangian L d by where g k ∈ G for 0 ≤ k ≤ N , N is the number of discretization steps. The discrete variational principle states that the solutions of the discrete system determined by L d must extremize the action sum given fixed endpoints g 0 and g N . By extremizing S d over g k , 1 ≤ k ≤ N − 1, we obtain the system of difference equations or, in coordinates, , and x and y represent the n first and n last variables, respectively, of the function L d .
These equations are usually called the discrete Euler-Lagrange equations. Under some regularity hypotheses (the matrix D 12 L d (g k , g k+1 ) is supposed to be regular), it is possible to define from (6) We will refer to the F L d flow, and also (with some abuse of notation) to the equations (6), as a variational integrator. This kind of integrators are symplectic-momentum preserving in some sense, which makes them interesting from the Geometric Integration perspective (see [18,33] for more details.) 2.3. Nonholonomic integrators. Discretizations of the Lagrange-d'Alembert principle for Lagrangian systems with nonholonomic constraints have been introduced in [10,34] as a nonholonomic extension of variational integrators. To define a discrete nonholonomic system providing a discrete flow on a submanifold of G×G and, as well, a corresponding version of discrete Euler-Lagrange equations (6), one needs three ingredients: a discrete Lagrangian, a constraint distribution D ⊂ T G and a discrete constraint space D d ⊂ G × G. Definition 2.2. A discrete nonholonomic system is given by the triple (G×G, L d , D d ): (1) D d is a submanifold of G × G of dimension 2n − m with the additional property that (2) L d : G × G → R is the discrete Lagrangian.
We define the discrete Lagrange-d'Alembert principle (DLA) to be the extremization of the action sum among all sequences of points {g k } 0:N with given fixed end points g 0 , g N , where the variations must satisfy δg k ∈ D g k (in other words δg k ∈ ker µ α ) and (g k , g k+1 ) ∈ D d for all k ∈ {0, ..., N − 1}. This leads to the following set of discrete nonholonomic equations For the sake of clarity, the condition (8b) may be rewritten as φ α d (g k , g k+1 ) = 0, where φ α d : G × G → R is the set of m functions whose annihilation defines D d (in other words, a suitable discretization of the nonholonomic constraints (1)).
Equations (8), where λ α = (λ k+1 ) α is chosen appropriately by projecting onto D d , define a local discrete nonholonomic flow map F nh where g k+1 satisfyes (8) provided that (g k−1 , g k ) ∈ D d , if and only if a regularity condition, say the matrix   To obtain distribution preserving integrators within this setting is not straightforward, since the equations (8) are defined on the discrete space G × G and moreover determine a discrete flow on the discrete distribution D d . Therefore, to locally relate (g k , g k+1 ) ∈ G × G and (g k , v g k ) ∈ T G is in order, a task accomplished in [20] by means of discrete difference maps ρ [34].
is a neighborhood of the diagonal I d in G × G, and V (Z) denotes a neighborhood of the zero section of T G, i.e Z : G → T G s.t. Z(g) = 0 g ∈ T g G, which satisfies the following conditions: Using them, the so-called velocity nonholonomic integrators are: defining a flowF L d : Tg k G → Tg k+1 G; (g k , vg k ) → (g k+1 , vg k+1 ). In general, this flow is not D−presrving, but sufficient conditions for that to hold are given, which in general involve an appropriate redefinition of the discrete nodes {g k } → {g k } and the particular discrete constraints given by φ α d . The result is stated in the following proposition (see [20]): Proposition 2.1. Assume that (g k , vg k ) ∈ Dg k . If D d is defined by φ α d := µ α • ρ : G×G → R, then (g k+1 , vg k+1 ) defined by the velocity nonholonomic integrator (10), i.e.F L d , belongs to Dg k+1 . In other words,F L d : Dg k → Dg k+1 .

Symmetries
3.1. Symmetries of Lagrangian systems. We say that the Lagrangian is invariant under a group action Φ : Such a symmetry allows to define a reduced Lagrangian system on the reduced phase space T G/G = (G × g)/G ∼ = g, say l : g → R (where we have employed the left trivialization mapping tr : is the left action and g is the Lie algebra of G). More concretely, the reduced Lagrangian is obtained using the symmetry by L(g,ġ) = L(g −1 g, g −1ġ ) = L(e, ξ) =: l(ξ), where we are employing the shorthand notation T g hġ := hġ (we shall equally employ the notation T g h ≡ ( h ) * ), e ∈ G is the identity element and g −1ġ := ξ ∈ g is called the reconstruction equation. Oviously, employing the right-invariance leads to equivalent results.
The reduction process allows to diminish the degrees of freedom of the given system and to obtain a first-order ODE describing the dynamics, instead of a second-order one. The paradigmatic example in the context of this work is the Euler-Poincaré equations. The Euler-Lagrange equations (4) are second-order, and obtained by the extremization of the action functional t2 t1 L(g(t),ġ(t))dt with fixed endpints δg(t 1 ) = 0, δg(t 2 ) = 0 and arbitrary variations δg ∈ T g G. Now, we assume a left-symmetric Lagrangian. Thus, we consider the left-translated variations δg ∈ T g G, namely g −1 δg = η, where η ∈ g. The endpoint condition δg(t 1 ) = δg(t 2 ) = 0 implies directly that η(t 1 ) = η(t 2 ) = 0. Furthermore, the variations δξ are not free anymore but determined by the equations ξ = g −1ġ and g −1 δg = η. After an easy calculation we arrive at where the adjoint action on the algebra stands for ad ξ η = [ξ, η]; [·, ·] : g × g → g represents the Lie algebra bracket. Now, the extremization of the action functional for an arbitrary curve η(t) in the Lie algebra with fixed endpoints η(t 1 ) = η(t 2 ) = 0, provides the so-called Euler-Poincaré equations Under regularity conditions, these equations provide the solution curve ξ(t) ⊂ g. Finally, the reconstruction equationġ = g ξ will provide the configuration curve g(t). Thus, in case of a symmetric Lagrangian, the second-order equations (4) may be replaced by the two first-order ones d dt ∂l ∂ξ = ad * ξ ∂l ∂ξ andġ = g ξ, which can be solved consecutively. However, they are not "Lagrangian" in the strict sense, since the variations (12) are not free anymore.
3.2. The nonholonomic case: the Euler-Poincaré-Suslov equations. In the nonholonomic case, besides a left-invariant L we shall consider as well a leftinvariant distribution D (which we will call a left-left or ll system). D ⊂ T G is left-invariant if and only if there exists a subspace g D ⊂ g such that D g = ( g ) * g D ⊂ T g G for any g ∈ G. Let a α ∈ g * , α = 1, ..., m, be independent annihilators of the subspace g D , i.e. g D = {ξ ∈ g | a α , ξ = 0 , α = 1, ..., m} .
Consequently, the left-invariant constraints on T G are defined by the equations a α , g −1ġ = 0. This last equation establishes the relationship between the nonholonomic constraints of a left-left system and the usual linear ones (1). More concretely, the nonholonomic constraints are locally given by φ α (g,ġ) = µ α (g)ġ, and therefore Thus, µ α (g) = * g −1 a α . The left invariance is easy to prove through φ α (g g,gġ) = µ α (g g),gġ = * (g g) −1 a α , ( g ) * ġ = * The dynamics of a ll system is determined by the so-called Euler-Poincaré-Suslov equations according to [25]: Furthermore, as in the unconstrained case, the dynamics of the group variables g is obtained through the reconstruction equationġ = g ξ. We note that the equations in (17) are the analogue to those in (2) in the case of ll reduced systems. For more details on the Euler-Poincaré-Suslov equations and their integrability we refer to [13,42].

Variational integrators and reduced systems
When one intends to discretize a Lagrangian problem defined on a Lie group, the most natural way is to directly discretize the equations (4). Through the discretization of variational principles, Marsden and collaborators obained in [28,29] the discrete Euler-Poincaré equations, which define an algorithm on the reduced space that is shown to be equivalent to the discrete Euler-Lagrange equations [33,38] in the sense of reconstruction. This algorithm is obtained by assuming left-invariance of the discrete Lagrangian function L d introduced in §2.2. More concretely, with respect to the induced group action on G × G we assume This projection map defines the reduced discrete Lagrangian l d : G → R for any G−invariant L d by l d • π = L d ; furthermore if we define the points in the quotient space as W k := g −1 k g k+1 , we obtain and the reduced action sum Applying calculus of variations to this action functional we obtain the so-called discrete Euler-Poincaré equations ( [28]), for the sequences where W k and r W k+1 are the left-and right-translation of the Lie group and denotes the derivative. Roughly speaking, the discrete Euler-Poincaré equations (21) generate, under regularity conditions, a discrete flow W k → W k+1 on G, generating the sequence {W k } 0:N −1 , which, after reconstruction g k+1 = g k W k , approximates the original flow g(t) defined by the original Euler-Lagrange equations (4). To see the relationship between the continuous and discrete Euler-Poincaré equations, (14) and (21) respectively, is not straightforward (one might argue, for instance, that they are defined in different spaces, g and G respectively) and needs more elaboration as shown next. 4.1. Discretization using natural charts. The above mentioned relationship between G and g is achieved by means of the so-called retraction map τ : g → G, which is an analytic local diffeomorphism near the identity such that τ (ξ)τ (−ξ) = e, where ξ ∈ g. In this sense, τ provides a natural chart on G, and W k are regarded as small displacements on the Lie group, linking g k and g k+1 . Thus, it is possible to express each term through a Lie algebra element that can be regarded as the averaged velocity of this displacement. Two standard choices for τ are the exponential map and the Cayley map. Regarding ξ as the "velocity" of a given displacement in the group, it can be expressed as which in general is an element of a nonlinear space, can now be represented by the vector ξ k . Moreover, the discrete Euler-Poincaré equations (21) can be expressed in terms of the velocities ξ k , the retraction map and its derivatives [7], which are defined by: Definition 4.1. Given a retraction map τ : g → G, its right trivialized tangent dτ ξ : g → g and its inverse dτ −1 ξ : g → g, are such that for g = τ (ξ) ∈ G and η, ξ ∈ g, the following holds ). Using these definitions, variations δξ and δg are constrained by where g η k := g −1 k δg k ; this expression is obtained by straightforward differentiation of ξ k = τ −1 (g −1 k g k+1 )/ . Note that η k forms the sequence {η k } 0:N where η 0 = η N = 0 since δg 0 = δg N = 0. Now, let us define a Lagrangian functioñ l d : g → R as a suitable approximation of the reduced action functional (13).
Remark 4.2. The naming ofl d : g → R is controversial. We are allocating the name discrete for L d : G × G → R and l d : G → R, the latter in the reduced case. These are considered as the discrete counterparts of the continuous Lagrangians L : T G → R and l : g → R. We note that T G and g are the Lie algebroids associated to the Lie groupoids G × G and G, respectively. Thus, in the general framework introduced in [43], we associate the discrete Lagrangian L d : G → R to the continuous one L : AG → R, where G and AG are the linked Lie groupoid and Lie algebroid, respectively. Therefore, to calll d : g → R a "discrete" Lagrangian function could be misleading.
After the introduction of these ingredients, the following result establishes an alternate version of the Euler Poincaré equations on the Lie algebra.
(2) The discrete Euler-Poincaré equations defined on the Lie algebra hold The last equation may be rewritten as where it has been taken into account that Ad * [16] for the proof). This equations, which can be undersood as a variational integrator for a reduced Lagrangian system, define, under regularity conditions, a discrete flow on the Lie algebra ξ k → ξ k+1 , generating a sequence {ξ k } 0:N −2 , which approximates the continuous solution of the Euler-Poincaré equations (14). The discrete configuration variables are recovered by the discrete reconstruction equation, namely 4.2. The nonholonomic setting. The extension of the previous results, and also of those in [6], to the nonholonomic left-left setup defined on Lie groups has been developed in [13].
whose annihilation defines completely D d . The main result is stated in the next theorem [13].
G → R be the reduced discrete Lagrangian, and D ⊂ T G, D d ⊂ G × G be the constraint distribution (also left-invariant) and the discrete constraint submanifold respectively. Then the following assertions are equivalent: (20), with respect to variations δW k induced by the constraint variations δg k ∈ D g k . (4) {W k } 0:N −1 satisfies the reduced nonholonomic equations or discrete Euler-Poincaré-Suslov equations: for k = 0, ..., N − 2, where again λ α = (λ k+1 ) α are chosen appropriately.
The constraints ϕ α d : G → R determine the so called constrained displacement subvariety S [13], say With regard to the Lagrange multipliers in equations (28) and (29), to determine if they are equal needs more ellaboration. The right hand sides of (28a) and (29a) represent the reaction forces provoked by the nonholonomic constraints, forces that involve the Lagrange multipliers. Thus, it is interesting to study whether in the unreduced and reduced setups these forces are equal. For that we need to check both if, in the two setups, the discrete dynamical equations and the constraint subspaces are equivalent. The first item is treated in the next result. Proof. We prove this in the unreduced-to-reduced way and by direct computations. By the left invariance of L d we have Taking variations of the right hand side we obtain where in the last line we have taken into account that W k = g −1 k g k+1 and l d (W k ) = L d (e, W k ). Rearranging the summation index after setting δg 0 = δg N = 0 we arrive at and consequently (29a) just by shifting k → k + 1. Taking variations of the left hand side of (31) we arrive directly at (28a), proving the claim.
On the other hand, we note that the subspaces determined by (28b) and (29b), say D d ⊂ G × G and S ⊂ G respectively, despite related to each other, are different in principle. Consequently we conclude that the multipliers resulting from the projection of the equations (28a) and (29a) onto these spaces will be different in general.
Considering a retraction map and equation (24), the following result is a corollary of theorem (4.3).
Proof. By direct computations we obtain where in the last line we have rearranged the summation index taking into ac- where we have used the shift k → k + 1, and therefore the claim holds.
Note that the second equation in (32) defines a subset, which we will denote by g D ⊂ g, g D := kerφ α d , such that in general g D = g D . Note as well, that the suitable Lagrange multipliers (λ k+1 ) α must be chosen appropriately in this case by projecting onto g D . The equations (32) define a discrete flow F nh , only under some regularity conditions, conditions which may be obtained using the implicit function theorem and locally ammounts to the invertibility of the following matrix Using coordinates we can write the upper-left term as Here, we are defining locally the pull-back of the trivialized inverse tangent retraction map as a linear operator by we shall see that these definitions are useful in the case of matrix groups.

Application to SO(3)
The group SO(3) is formed by 3 × 3 real orthogonal matrices such that their determinant is equal to one (with some abuse of naming, we will sometimes call these matrices orthonormal). It can be also interpreted as the +1-determinant connected component of O (3), which is the group of rigid transformation in R 3 . From the physical point of view this is relevant since, among other reasons, it is the configuration manifold of the the Lagrangian and Hamiltonian formulations of the rigid body problem.
The corresponding Lie algebra so (3), which is the algebra of real 3 × 3 skewsymmetric matrices, is isomorphic to the Euclidean space R 3 through the isomor-phism· : R 3 → so (3), ω →ω, given bŷ for ω = (ω 1 , ω 2 , ω 3 ) ∈ R 3 . In this representation, the antisymmetric bracket operation is the standard vector product in Furthermore, in following sections the definition of a distance between group and algebra elements will be relevant, respectively. In the latter case, it is common to pick the Killing form, since it is invariant under all the automorphisms of so(3). Namely, whereξ,η ∈ so(3) andξη represents the usual matrix product. As it is wellknown, this bilinear form corresponds to the usual Euclidean product on R 3 , say for ξ, η ∈ R 3 , we have (ξ,η) e = ξ · η = ξ 1 η 1 + ξ 2 η 2 + ξ 3 η 3 . Thus, the distance betweenξ,η ∈ so(3) may be defined by Concerning the Lie group, to pick a distance is a subtle issue, even more if one wants to stick close to Euclidean metrics. In the reduced setting, the update of the Lie group elements by the developed variational integrators is given by Ω k+1 = Ω k τ ( ω) (25), where Ω k ∈ SO(3) andω ∈ so(3). In case of τ =exp, it is clear that Ω k+1 would also belong to SO(3). Nevertheless, for later applications, we will pick a different retraction map, for instance a suitable truncation of the exponential. This choice in general prevents Ω k+1 to lie in SO(3). Following this idea, we establish the self-distance in SO(3) as dist(Ω, Ω) = |I − Ω T Ω|, where we use the Euclidean metric of matrices, induced by |Ω| 2 = ij |Ω ij | 2 . As it is natural, if Ω ∈ SO(3), then dist(Ω, Ω) = 0. In other words, what (37) measures is how far is a matrix from being orthonormal in terms of an Euclidean metric in a matrix space. This generalizes to measuring how far is Ω 1 from being Ω 2 within SO(3). It is easy to prove that dist: SO(3) × SO(3) → R is indeed a metric, and moreover allows to establish the following lemma, which is very useful for our purposes 1 . Then where H is a 3 × 3 real matrix defined by 9 real entries H ij = O( r+1 )h ij ( ), where at least one of the nine h ij is independent of . 1 In case of SO(3) matrices it might be more natural to consider the Frobenius matrix norm ||Ω 1 −Ω 2 || 2 F = 2tr(I−Ω T 1 Ω 2 ). However, this metric only allows to establish bounds for the diagonal entries of the difference matrix I − Ω T 1 Ω 2 , which is not enough to establish the consistency bounds in the context of this work.
Proof. Let us define I −Ω T 1 Ω 2 = H, then it follows directly that Ω 1 −Ω 2 = Ω 1 H. On the other hand, the condition dist(Ω 1 , Thus, all the matrix entries may be factorized as H ij = O( r+1 )h ij ( ), where at least one of the entries h ij must be independent of , since otherwise the square root would be of order O( r+2 ) as → 0. (3). Concerning the retraction map τ , we employ the Cayley map since it is simple and also computational efficient. Say, it is suitable for iterative integration and optimization problems since its derivatives do not have any singularities that might otherwise cause difficulties for gradient-based methods [18]. The Cayley map cay : g → G is defined by cay(ξ) = (e − ξ 2 ) −1 (e + ξ 2 ) and is valid for a general class of quadratic groups. The quadratic Lie groups are those defined as

The Cayley map in SO
where P ∈ GL(n, R) is a given matrix (here, GL(n, R) denotes the general linear group of degree n; O(n) or SO(n) are examples of quadratic Lie groups). The corresponding Lie algebra is g = Ω ∈ gl(n, R) | P Ω + Ω T P = 0 .
The right trivialized derivative and inverse of the Cayley map are defined by for ξ, η ∈ g. Using the identification (34), its particular form in the case of SO (3) is One of the most interesting properties of the Cayley map is that its image belongs SO(3), i.e. cay(ω) cay(ω) T = cay(ω) T cay(ω) = I. Furthermore, the linear maps dτ ξ and dτ −1 ξ are represented by the 3 × 3 matrices We point out to the reader that the same notation | · | will be used for the Euclidean norms on R 3 , so(3) and matrices. Taking into account (37), Ω k+1 = Ω k cay(ω) and (39), after a straighforward calculation we arrive at for Ω k ∈ SO(3). Therefore, the update Ω k+1 defined by (25) in the case τ =cay remains in SO(3).

The Suslov problem on SO(3)
. The Suslov problem is a well-known example of a left-left system, introduced in [41]. It describes the motion of a rigid body suspended at its center of mass in the presence of a constraint that forces the projection of the body angular velocity along a direction fixed in the body to vanish. The generalization to the group SO(n), called the generalized Suslov problem, has been described, both from continuous and discrete points of view, in [12,13,21,44].
Let I = (I ij ) be the inertia tensor I : so(3) → so(3) * of the body, I −1 = (I ij ) its inverse, and ω ∈ R 3 be the body angular velocity vector. The dynamics is determined through the reduced Lagrangian (11), which in this case reads l : Let a be the direction, fixed in the body, of the vanishing component of the angular velocity. Thus, the nonholonomic constraints read a, ω = 0, where we are using the bracket as the pairing between (R 3 ) * and R 3 . Taking into account (17) and (41), the Euler-Poincaré-Suslov equations, written in R 3 , are where λ ∈ R is the Lagrange multiplier. Componentwise, we have and a 1 ω 1 + a 2 ω 2 + a 3 ω 3 = 0. Without loss of generality, we may choose a as the third component of the body frame {e 1 , e 2 , e 3 } in R 3 , say a = (0, 0, 1). Then, the constraint becomes ω 3 = 0 and Moreover, (42) becomes where, from now on, i = {1, 2}. After a straightforward computation we arrive aṫ where I m = I 11 I 12 I 21 I 22 is considered to be non-degenerate; moreover (45) In other words, what happens to the Euler-Poincaré-Suslov equations (42) in the presence of the nonholonomic constraint ω 3 = 0 is that they decouple into a differential part (44), this is a system of nonlinear ODEs which we will denote bẏ ω = f (ω) for simplicity, and an algebraic part (45). 5.3. Consistency order. Now, we extend the consistency results from [20] to the present case, when G = SO(3). We consider a general Euler-Poincaré-Suslov problem in SO(3) and therefore we may have two scalar constraints. Thus, we keep the a α notation instead of the single a as in the previous subsection, where we introduced the specific example of the Suslov problem. In the definition 2.1 we have introduced the consistency order of a general integrator of the nonholonomic DAE (2). However, the reduction process leads to the Euler-Poincaré-Suslov equations (17), which are first-order instead of second-order equations. We take into account all these facts to establish a relationship between the (p, s) consistency of a discretization of the nonholonomic DAE (2) (for instance velocity nonholonomic integrators (10)) and a general discretization of (17), and conversely (always in the left-left framework). For that purpose, we employ the notions of distance (36) and (37). Namely, definition 2.1 implies where the continuous dynamics (Ω(t), v Ω (t), λ(t)) is determined by (2) and the discrete one (Ω k+1 , v Ω k+1 , λ k+1 ) by the prescribed integrator. On the other hand, the notion of (p, s) order integrator may redefined for the Euler-Poincaré-Suslov equations to be an integrator such that |ω(t k + )−ω k+1 | ∼ O( p+1 ) and |λ(t k + )−λ k+1 | ∼ O( s+1 ), where (ω(t), λ(t)) is given by (17) and (ω k+1 , λ k+1 ) by the prescribed integrator. The relationship between both the continuous dynamics is established by the reconstruction equation, i.e. v Ω (t) =Ω(t) = Ω(t)ω(t), while the discrete ones are linked through Ω k+1 = Ω k cay( ω k ). Furthermore, we note that the left translation ofω ∈ so(3) by a group element Ω ∈ SO(3), i.e. Ωω, is nothing but the matrix product Ωω, where Ω ∈ SO(3) and ω ∈ R 3 as established by the hat isomorphism (34).
The link between consistency orders at the group and algebra level, which is not straightforward, is discussed in the following result.  (3), and such that Ω k+1 = Ω k cay( ω k ). Then the following assertions hold: 1. the (p, s) integrator (DAE) generates a (p, s) integrator of (17), 2. the (p, s) integrator (EPS) generates at least a (min(p, 1), s) integrator of (2).
Proof. In first place we take into account the dynamical variables, say Ω(t), v Ω (t) andω(t).
1) By construction, the (p, s) integrator (definition 2.1) implies |I − Ω k+1 Ω(t k + ) T | ∼ O( r+1 ) and therefore it follows from lemma 5.1 that where we identify H directly with O( r+1 ). On the other hand, regarding the velocities we have v Ω (t k + ) − v Ω k+1 = O( l+1 ). (This last expression could be misleading from the geometric point of view, since we are adding vectors belonging to two different vector spaces, say T Ω(t k + ) SO(3) and T Ω k+1 SO(3); however, from the reconstruction equation these are nothing but rotated vectors in R 3 , i.e. Ω(t k + ) ω(t k + ) and Ω k+1 ω k+1 and therefore we can understand the sum as embedded in R 3 .) Thus, taking into account the reconstruction equation we arrive at v , where the right hand side is identified with a vector in R 3 . Furthermore and Ω(t k ) = Ω k , ω(t k ) = ω k , from the last expression we obtain . Now, taking the norm on both sides and noting that p = min(r, l), we finally obtain 2) First we consider the consistency order in the group, namely |I − Ω k+1 Ω(t k + ) T |. For this, it will be enough to set up the Taylor expansion Ω(t k + ) = Ω(t k ) + Ω (t k ) + 2 2Ω (t k ) to second order. On the other hand, Ω k+1 = Ω k cay( ω k ) = . Taking this into account and considering thatω T = −ω forω ∈ so(3) we have that Moreover, sinceΩ = Ωω, we haveΩ(t k ) = Ω(t k )ω 2 (t k ) + Ω kω (t k ). With all these ingredients, recalling that Ω(t k ) = Ω k and after some computations we arrive at In general,ω(t k ) = 0 and consequently, However, ifω(t k ) = 0 we get |I − Ω k+1 Ω(t k + ) T | ∼ O( 3 ) and then the consistency order is at least min(p, 2). Thus, we can say that the minimum consistency order is given by min(p, 1), and higher orders depend on nontrivial cancelations. Regarding the velocity, since we are considering a p−th order integrator in the algebra, i.e.ω(t k + ) −ω k+1 = O( p+1 ), the following chain of equalities holds and furthermore, employing (47) and that Ω k+1 = Ω k cay( ω k ), we obtain Taking norms, we arrive at |v Ω (t k + ) − v Ω k+1 | ∼ O( min(1,p)+1 ). (The same argument presented above regarding the sum of vectors at T Ω(t k + ) SO(3) and T Ω k+1 SO(3) holds.) Finally, we consider the algebraic variable, i.e. the Lagrange multiplier λ. Since we consider the left-left framework for SO(3), equation (16) reads a α ,ω = a α , Ω TΩ = a α , ( Ω T ) * Ω = * Ω −1 a α ,Ω = µ α (Ω),Ω , which in other words means that the constraint distributions in both the unreduced and reduced setups are related by the group lift, determining the same subspace of so(3) D ⊂ so (3), a subset prescribing a particular subspace so * (3) D ⊂ so * (3) through the Killing metric. Thus, the reaction forces involving the continuous dynamics of λ(t) will also be equal if the continuous dynamical equations are the same. This can be proved in the case l(ω) = 1 2 Iω, ω noting that the left hand side of the reduced and unreduced equations, both elements of so * (3), say are equal, a fact that follows from straightforward computation setting L(Ω,Ω) = 1 2 IΩ TΩ , Ω TΩ . Consequently, the consistency bound |λ(t k + ) − λ k+1 | ∼ O( s+1 ) remains the same. This assures the complete claim.
Remark 5.3. We observe that the dynamical part of theorem 5.2, namely the p → p, p → min(p, 1) consistency relationships, is general for any kind of Lagrangian functions, since in any case the unreduced and reduced continuous setups are related to each other by the reconstruction equation. In the case of the Lagrange multipliers, the consistency order depends on the particular Lagrangian function, and therefore the generality of the result is left as subject of further study.
Remark 5.4. We note that the result in theorem 5.2 is generalizable to any group which can be embedded in R n , for n big enough. However, when dealing with a general Lie group G the situation is more subtle, since to define a suitable distance is not easy in principle. In terms of Riemannian geometry, nevertheless, the distance is always prescribed once we have defined a positive definite bilinear form in the algebra (·, ·) e : g × g → R + . Furthermore, there is a bijective correspondence between left-invariant metrics on a Lie group G and the inner products on the Lie algebra g (see [36,37] for further details). This leads to the definition of a leftinvariant Riemannian metric on G (·, ·) g : T g G×T g G → R + through left translation, say for any u, v ∈ T g G, and furthermore to the Riemannian distance where γ : [t 1 , t 2 ] → G is the geodesic with γ(t 1 ) = g, γ(t 2 ) =g, such that g,g are close enough in terms of natural charts ( §4.1). This fact allows to relate a (p, s) integrator for (2) and a (q, s) integrator for (17) in both directions, say unreduced to reduced and conversely, by the following lemma.
However, this is not enough to establish the full relationship between the consistency of integrators, since the bound of |v g (t k + ) − v g k+1 | needs to be established also in both directions. This is left as subject of further study, but necessarily this bound has to depend on the structure of the bilinear form (·, ·) e and its link to the particular Riemannian metric (·, ·) g .
Second, it is interesting to note that the variational procedure described in the- As an interesting feature, we also point out that there is no need in this case to redefine the nodes {ω k } 0:N −1 , as happens in general for nonholonomic systems defined in Q = R n (see [20] for further details).
The preservation of the constraints in the variational setting, when we restric ourselves to the discrete Euler-Poincaré-Suslov equations (theorem 4.3), is less obvious. Locally, these equations generate a flow inside the submanifold S ⊂ SO(3) (30) if a regularity condition is satisfied. Again, this condition is obtained employing the implicit function theorem and is equivalent to the nondegeneracy of the matrix where (∇r W k+1 ) * makes sense locally, in analogy with the case of the Euler-Poincaré-Suslov equations defined on the Lie algebra and their regularity condition (33). Moreover, the relationship between S and D d ⊂ SO ( Taking into account all these elements, we can establish at this point a result concerning the D−preservation of the discrete Euler-Poincaré-Suslov equations (29). Proof. The discrete Euler-Poincaré-Suslov equations generate the sequence {W k } 0:N −2 evolving on S, which by "unreduction" generate {(Ω k , Ω k+1 )} 0:N −1 evolving on D d according to W k = Ω T k Ω k+1 . Thus, φ α (Ω k , Ω k+1 ) = 0. Furthermore, (53) implies that φ α (Ω k , Ω k+1 ) = µ α (Ω k+1 )Ω k+1 = 0, and the claim follows.
Note that (53) gives us some freedom concerning the choice ofΩ k , that is the nodes for discrete integration (a fact that also showed up in the case of R n treated in [20]). Nevertheless, the natural choice for ρ is ). It is easy to understand that this choice of a finite difference map establishes a local isomorphism S ∼ = so(3) D through the retraction map τ , since

Discretizations of the Suslov Problem on SO(3)
In this section we carefully study some particular discretizations of the Suslov problem in SO(3) introduced in §5.2.

General discretizations DREPS.
It is interesting to note that any discretization DREPS(ω k , λ k+1 ; ω k+1 ) = 0 2 (50) (which implies ω k+1 3 = ω k 3 = 0), respecting the local regularity condition (51) and applied to (43), leads to particular discretizations of (44) and (45). Conversely, a general discretization of (44) prescribes a DREPS of the Euler-Poincaré-Suslov equation on SO(3), for which we can derive an interesting result. Prior to its statement, we recall thatω = f (ω) is a system of nonlinear ODE defined in a vector space, and consequently we can apply a numerical method of arbitrary consistency order (Euler, midpoint rule, Runge-Kutta, etc.). Furthermore, the decoupling of (43) into differential and algebraic parts allows some freedom in the choice of its DREPS. In particular, we choose where we pick λ k+1 according to (45) as a natural choice. We assume l 1 , l 2 to be smooth functions chosen in such a way that they generate a p−th order consistent one-step discretization of (44), i.e. |ω(t k + ) − ω k+1 | ∼ O( p+1 ), with p ≥ 1. Namely Now we have all the necessary ingredients to establish the result.
Proof. The dynamical part is obvious, namely the order |ω(t k + )−ω k+1 | ∼ O( p+1 ) is given by assumption. On the other hand, given the choice λ k+1 = λ(ω k+1 ) according to (45) and the ω consistency bound, using the Taylor expansion of λ we have Once the (p, p) integrator is confirmed for the reduced system, the second statement, this is the (min(p, 1), 1) consistency for the unreduced problem, follows directly from Theorem 5.2.
Remark 6.1. The previous result may be refined concerning the algebraic part, generating a (p, s) integrator with s > p. For this purpose, we set in (54), with l λ chosen as shown next. In the new scenario, we have (where we omit the l λ variables for the sake of simplicity) where we take the Taylor expansion of λ(ω k+1 + O( p+1 )) up to the term s+1 (s = n p + n − 1, with n an integer). Thus, it is obvious that choosing l λ such that we obtain the (p, s) integrator as claimed.
6.2. Variational nonholonomic integrators. As presented in §4, the variational integrators setting provides us with a framework for the numerical integration of reduced systems. In particular, corollary 4.4 prescribes a particular DREPS, namely, Regarding the discrete constraintsφ d , although we have some freedom for their choice,φ d (·) = a, · is the most suitable, since it implies the reduced constraints preservation.
For the case concerning this work, we will assume furthermore that ω 3 = 0, which yields the following form of the right trivialized inverse derivative of the Cayley map in SO (3): As pointed out in the discussion below the equation (33), we see that the right trivialized derivative of some retraction maps, in this case the Cayley map in SO(3) is given by a matrix, and therefore in local coordinates it has two indices: dcay −1 Settingl d (ω) = l(ω) as a first order approximation of the action s(ω) = t1+ t1 l(ω) dt (13), where l : so(3) → R is given by (41), and applying (56) with ϕ d (·) = a, · (note that the constraint is single in the Suslov problem), we obtain the following algorithm which accounts for the discretization of the dynamical part of the Euler-Poincaré-Suslov equations (we recall that I m = I 11 I 12 where the rescaling λ k+1 → −λ k+1 / 2 has been introduced, accounts for the algebraic part (note that in this case we do not have freedom to fix the algebraic part since it is directly given by the variational integrators setup). This rescaling may be understood in the context of the construction of the variational integrator from a discretization map ρ : . More concretely, as can be shown in view of proposition 5.2, any discretization φ α ∝ µ α • ρ preserves the constraint for the discrete flow (Ω k ,Ω k ) → (Ω k+1 ,Ω k+1 ). Thus, setting φ α = −1 2 µ α • ρ accounts for the mentioned rescaling.
Once we have derived the discrete equations, we may establish the following result regarding their consistency order with respect the continuous dynamics.
Proposition 6.2. The numerical method (58), (59), is a (2, ) consistent method with respect to the Euler-Poincaré-Suslov problem (17) defined on SO(3), and consequently (1, ) with respect to (2); where by we mean that the integrator is not neccessarily consistent with respect to the multipliers.
Proof. To prove this result, we just consider the Taylor expansion ω(t k + ) = ω(t k ) + ω(t k ) + 2 2ω (t k ) + O( 3 ), whereω is given in (44), and compare, order by order, with ω k+1 provided by (58). In first place, we calculate the second order time derivatives, namelÿ On the other hand, the first two terms in (58) imply that ω k+1 = ω k +O( ) (ensuring consistency) and, moreover, that which implies first order consistency in comparison to (44). Futhermore, ω k+1 = ω k + O( ) implies that the third and fourth terms in (58) cancel out at O( 2 ) order; thus, we realize that the relevant term in (58) at this order is just .
Using (60), we obtain Plugging these terms into (61), we obtain Now, substraction from the expressions ofω 1 ,ω 2 presented above, we obviously have . Furthermore, the 1 4 factors in (58) prevent the O( 3 ) terms on the discrete and continuous sides to coincide. Regarding the multipliers, it is straightforward to see from (59) that λ k+1 = ω k which is different from zero, in general, and consequently the discrete multiplier is not consistent with the continuous one. This shows the first claim. Concerning the second, it is enough to apply theorem 5.2.
Therefore, the variational integrators setting generates a second-order consistent method on the dynamical side (a fact which is interesting, since we are considering a first-order consistent discrete Lagrangianl d , and therefore we would expect a firstorder consistent numerical method in the spirit of [33]) while it is not neccessarily consistent in the algebraic one. Needless to say, this is a drawback in the numerical scheme. However, due to the decoupling between the two parts mentioned above (we can obtain the ω k+1 values independently of the λ k+1 's), we may perform, besides the 2 reescaling, a discrete shift as described in remark 6.1, generating in consequence a (2, s) method for (17), and therefore a (1, s) method for (2). Unfortunately, such a shift cannot be understood in general as an alternate choice of the discretization map ρ, as long as the particular map is not just a rescaling anymore but involves nonlinear terms depending on the ω variables.
We illustrate these developments by the following discussion and plots. Before going into details, we point out that the energy is preserved in nonholonomic dynamics along the solutions. This can be easily shown by direct calculation when we define the Lagrangian energy by E L = ∂L ∂ġ iġ i − L(g,ġ), where we consider again the original setting expressed in equations (2). Taking the time derivative of E L , we arrive at according to (2). Thus, as long as g(t) is a solution of the Lagrange-d'Alembert equations, we have that µ α (g)ġ = 0, and consequently E L is conserved. In the reduced case, using similar arguments it is proven in [13] that the Euler-Poincaré-Suslov equations preserve the reduced energy where we consider the reduced Lagrangian l(ξ) (11). In our case, the reduced energy along the solutions reads E l (ω) = 1 2 Iω,ω = ω 1 (I 1i ω i )+ω 2 (I 2i ω i ) (with i = {1, 2}). The preservation of this energy should be taken into account as a favorable property in the performance of nonholonomic integrators, as shown below.
We consider an homogeneous rigid body with inertia matrix I = and initial values ω 1 (0) = 0.4 and ω 2 (0) = 0.5 (with proper unities). We shall display the performance of an order (2,1) DREPS, in terms of (54), and the variational nonholonomic integrator (58) and (59), which is also order 2 in the dynamical variables as proved in proposition 6.2. More concretely, the integrator DREPS corresponds to the midpoint rule in (54), namely l 1 (ω k , ω k+1 ) = − (note as well that in this case there is no dependence in l 1 , l 2 ). In order to achieve order 2 consistency with respect to (45), we determine λ k+1 according to (54).
On the other hand, the variational nonholonomic integrator (58) and (59) corresponds to the setupl d (ω) = l(ω),φ d (ω) = a,ω and the retraction map τ = cay. What we observe is that, for small time steps , the behavior of both integrators is indistinguishable except with respect to the multipliers. As to be expected due to proposition 6.2, the nonholonomic variational integrator produces an inconsistent discretization of the Lagrange multipliers. However, we display also the performance of both integrators over a big time interval, noticing that the variational integrator's performance is much better in all respects, mainly with respect to the preservation of energy, where we observe a fast decay of the midpoint rule. Moreover, we observe that the variational integrator even works quite well if the step size of the time discretization is taken to be relatively big. Figure 1. In this figure we display the performance of the midpoint rule (DREPS(ω k , λ k+1 ; ω k+1 ) = 0, with inertia matrix I and initial values ω 1 (0) and ω 2 (0) introduced above) for the nonholonomic rigid body with time step = 10 −3 . The solid red line is obtained through a RK4 integrator (which we consider an accurete approximation of the continous nonlinear dynamics in the short term), while the blue dots represent the performance of the midpoint rule. The plots (a) and (b) correspond to the dynamical variables ω 1 , ω 2 , while (c) displays the Lagrange multipliers λ. On the other hand (d) shows the inconsistent multipliers generated by the nonholonomic variational integrator. Finally, (e) and (f ) shows the preservation of the constraints and the energy E l (ω) up through round off errors, respectively.  Figure1) and the variational integrator (58), (59), for a time step = 10 0 = 1 (we recall that this integrator is also order 2 consistent in the dynamical variables). The former is represented by the green points and the latter by the blue ones, while the solid red line still represents the performance of a RK4. Variables ω 1 (a), ω 2 (b), λ (c) and E l (d) are displayed, while (e) represents the preservation of the constraints by the variational integrator up through round off errors. We observe a better performance of the variational setup, mostly with respect to the preservation of energy, a fact which, considering the bigger time step, leads to the conclusion that its convergence to the actual solution is much faster and its long-term behavior will be more accurate.

Discretization of the Euler-Poincaré-Suslov problem as perturbation
So far, we have focused on the discretization of the Euler-Poincaré-Suslov equations on SO (3), paying attention to its consistency order and D−preservation. In this section we return to any Lie group G (as long as it may be embedded in R n ) and its Lie algebra g, since the relevant results from above apply in this general case as well; in the last part we will particularize G = SO(3).
As it is well-known, by discretizing the dynamics we introduce some discrepancies with respect to the continuous system even when the discretization is performed in some kind of structure-preserving fashion. It is needless to mention that this is a central and fundamental question for all kinds of numerical investigations, especially concerning the long-term evolution of dynamical systems. Regarding this issue, we refer to [17] where a positive answer to the following question is given: Is it possible to embed a numerical scheme approximating the continuous-time flow of a set of autonomous ordinary differential equations (ODE) into the time evolution corresponding to a non-autonomous perturbation of the original autonomous ODE? The positive result may be phrased as: Any p−th order discretization of an autonomous ODE can equivalently be viewed as the time− period map of a suitable −periodic non-autonomous perturbation of the original ODE (where is the fixed discretization lenght); while the precise statement is: Theorem 7.1. Suposse that h ∈ C r (R n , R n ), r ≥ 1, and consider the autonomous ODEẋ Let F (t, x) be the fundamental solution of (62) satisfying F (0, x) = 0, and assume that there are an integer ρ ≥ 1, a continuous function C : [0, ∞) → [0, ∞) and a one-step difference approximation of step size which is consistent of order p, i.e.
Then, there exists a function d( , t/ , x), as smooth as h and periodic in t of period , such that if G(t, s; , x) 3 , G(s, s; , x) = x, is the fundamental solution of the non-autonomous, −periodic ODĖ then where G( , 0; , ·) : R n → R n is the Poincaré map (period map) for (63), corresponding to initial time s = 0.
Following these ideas, in [20] the case of nonholonomic dynamics when Q = R n was explored, and the following theorem was proved. It is stated here for later use. Then, any p−th order discretization of the nonholonomic ODE (evolving on D) may be embedded into the time evolution of a non-autonomous perturbation of the original nonholonomic DAE (2) of the forṁ Similarly, an analogous result can be obtained in the reduced setting. There are some differences though; first, we note that the Euler-Poincaré-Suslov equations (17) are first-order rather than second, as (2) are; second, the constraints a α , ξ = 0 determine a subset g D ⊂ g of a linear space instead of a regular submanifold. The procedure to obtain the perturbation of the nonholonomic dynamics produced by a given discretization can be split into three steps, say: (1) Define an ODE evolving on g D from the Euler-Poincaré-Suslov DAE (17) by projection, which we will call Euler-Poincaré-Suslov ODE; (2) Choose appropriate coordinates in such equations to apply the theorem 7.1; (3) Undo the projection process to recover the perturbed Euler-Poincaré-Suslov DAE. Thus, we obtain: where we have used a local form of the equations where C e ab are the structure constants of g.
Proof. Before going into details, we introduce some notation. We denote m ≡ (m ab ) := ∂ 2 l ∂ξ a ∂ξ b , while (m ab ) ≡ m −1 denotes its inverse (recall that (m ab ) is regular and positive-definite). With this, the first equation in (17) may be rewritten (1) Since µ α (g) = * g −1 a α represent a set of independent one-forms spanning D • g ⊂ T * g G, it is easy to see that the a α spanning (g D ) • ⊂ g * are also linear independent. Therefore, we can decompose the algebra as with respect to the metric (m ab ). Thus, from (66) we can derive an ODE evolving on g D by eliminating the Lagramge multiplier. More particularly, if we take the time derivative of the nonholonomic constraints we obtain a α ,ξ = 0 (note that a α are constant), which, after replacingξ by the right hand side of the equation (66), yields λ α (ξ) = −C αβ a β , f (ξ) , (68) where C αβ := a α , m −1 a β and (C αβ ) = C −1 . It is possible to prove the invertibility of C using geometric arguments (see for instance [20], lemma 2.4), or just by arguing that m is full-rank and a α constant-rank. Finally, we obtain the Euler-Poincaré-Suslov ODE evolving on g D given bẏ as long as a α , ξ = 0, where h(ξ) := f (ξ) + λ α (ξ) m −1 a α and λ α (ξ) is defined in (68).
Roughly speaking,d( , t/ , ξ) is the "component" ofd( , t/ , ξ) along g D , as we show next. By a direct calculation, we see that the first equation in (72) leads to (71) by eliminating the Lagrange multipliers, which gives where λ α (ξ) is defined in (68). Plugging this into (72) we finally obtain (71), including as well the identificationd = m −1d − C αβ m −1 a α a β , m −1d . This finishes the proof.
Of course, it would be ideal for a further analysis, if the perturbed Euler-Poincaré-Suslov equations in (72) admitted a Lagrangian structure of some kind, i.e. if it resulted from a suitable reduction of the Lagrange-d'Alembert principle. The question, when this is true and how this is related to specific properties of both, the underlying discretization of the unperturbed problem as well as the interpolation procedure, is left as a subject of further study.
Finally, when we restric to G = SO (3), it follows directly from theorems 5.2, 7.2 and 7.3 the next corollary.

Conslusions
We follow the ideas presented in [20] to study the discretization of the Euler-Poincaré-Suslov equations on SO(3). First, concerning the consistency order of the discretizations corresponding to the unreduced and reduced settings, we find that when the discrete reconstruction equation is given by a Cayley retraction map both consistency orders are related to each other, remaining the same in the unreduced to reduced direction, and spoiled to first order (when the original discretization has order bigger than one) in the reduced to unreduced one. Furthermore, we study distribution preserving integrators, giving precise conditions under which this property is ensured from the discretization of the Lagrange-d'Alembert principle. On the other hand, we also show that this preservation may be shown for more general numerical schemes, and give a precise algorithm to achieve it based on the general reduced framework. These definitions and proved results are illustrated through the example of the Suslov problem on SO(3). We consider two different integrators, one of them based on the midpoint rule while the other is based on a variational scheme. As proved previously, both are order-2 consistent in the dynamical variables, but it is numerically shown that the variational one converges faster to the actual solution and shows a better long-term behaviour. Finally, concerning the discretization understood as perturbation, we prove that any distribution preserving integrator of the Euler-Poincaré-Suslov equations may be understood as a non-autonomous perturbation of the continuous dynamics itself. In the case of SO(3) and considering a discrete reconstruction based on the Cayley map, the last result combines with the consistency bounds in order to provide the order of the perturbation produced in the unreduced dynamics by a reduced integrator, and conversely.