Preserving first integrals with symmetric Lie group methods

The discrete gradient approach is generalized to yield integral preserving methods for differential equations in Lie groups.


Introduction
Our point of departure is the system of differential equationṡ where the unknown x = x(t) is a curve on some Lie group x(t) ∈ G and the dot over x signifies differentiation with respect to t. The map f ∶ G → g where g ≡ T e G is the Lie algebra corresponding to G, and the dot should be interpreted as the derivative of right (resp. left) multiplication, e.g. f (x) ⋅ x ∶= T e R x f (x). Such equations occur for instance in mechanical systems where the Lie group could be the special orthogonal group SO(3) or the special Euclidean group SE(3). We shall consider in particular the case where the system (1) possesses one or more first integrals.
Here we define a first integral to be any function H ∶ G → R which is invariant on solutions d dt H(x(t)) = ⟨dH, F ⟩ = 0.
Any differential equation (1) on a Lie group having H as first integral can be formulated via a bivector (dual two-form) ω and the differential of H aṡ A Riemannian metric on G defines an inner product (⋅, ⋅) x on every tangent space T x G which is also varying smoothly with x. With such a metric one may define the Riemannian gradient vector field as the unique vector field gradH satisfying ⟨dH, v⟩ x = (gradH x , v) x for every v ∈ T x G. An example of a bivector ω to be used in (2) is then provided by means of the wedge product Note that ω is not uniquely defined by F and H. For a choice x = (x 1 , . . . , x d ) of local coordinates on the group, we may write In this way, we find the coordinate version of (2) x = S ω ∇H (4) for the skew-symmetric d × d-matrix S ω = ω T − ω and ω is the strictly upper triangular matrix with entries ω ij (x), 1 ≤ i < j ≤ d. The formulation (4) has been the starting point for energy preserving integrators devised by Gonzalez [8]. In fact, McLachlan et al. [13] showed that under relatively general circumstances, any vector field F on R d with a first integral H can be written in the form F (x) = S(x)∇H (5) for some skew-symmetric matrix S(x). Note that the bivector is not required to be nondegenerate. For Hamiltonian systems, the Hamiltonian itself is a first integral and an accompanying bivector can be inferred from the symplectic two-form Ω by inversion. In coordinates one can represent the symplectic two-form by a skew-symmetric matrix S Ω in a similar way as in (3) and (4) with respect to the basis {dx i ∧ dx j , i < j}. By definition, this matrix will be invertible, and its inverse is precisely (2) is easily generalised to the case with k independent first integrals H 1 , . . . , H k . We may now replace the bivector with a k + 1-vector ω and write (1) aṡ Also in this case we can find, by means of a Riemannian structure, an example of a feasible Earlier work on energy-preserving methods on Lie-groups was presented in [12], see also the references therein. In this paper we shall generalise the notion of discrete gradient methods to the situation where the phase space is a Lie group. The paper extends results presented in [4] by also considering high order methods, a larger class of manifolds and further examples.
2 Discrete differentials in Lie groups

A review of the situation in Euclidean space
The idea of discrete gradient methods, is to consider some approximation to the exact gradient in (5) ∇H ∶ R d × R d → R d satisfying the following two conditions Many such discrete gradients have been proposed in the literature, and we give here a few examples. The averaged vector field discrete gradient is discussed for instance in [13] and more recently in the PDE setting [7,3] Another discrete gradient which was proposed in [8] is the midpoint gradient An integral preserving integrator for (5) is readily given as whereS is a skew-symmetric matrix which approximates S. It is required to satisfy the consistency conditionS(u, u) = S(u). An easy calculation, using (6) and (10) now shows that the last identity follows sinceS is skew-symmetric. The skew-symmetric matrixS used in the integrator is not unique, and the freedom can be used to improve the approximation (10) in various ways, for instance to increase its order of convergence.

The Lie group setting
For Lie groups, none of the definitions (6) or (10) can be used since they both involve vector space operations, not generally defined on Lie groups. In a coordinate free setting, we also find it more convenient to replace the gradient and its discrete counterpart by dual quantities, the differential as indicated in (2). Below, we also use the notion of a discrete differential rather than a discrete gradient. As is the tradition for Lie group integrators [10,5], approximations are introduced through some finite dimensional action and a trivialisation principle is applied. By this, we mean that tangent vectors at some x ∈ G, i.e. v x ∈ T x G can be represented via either left or right translation of a vector ξ ∈ T e G ≅ g, where g is the Lie algebra of G. In the present paper we will just for convenience choose right translation. Defining the right multiplication operator and similarly, any vector p x ∈ T * x G is identified by some µ ∈ g * through where ⟨⋅, ⋅⟩ is some duality pairing between T * x G and T x G, as well as between g * and g. For example, if the Lie group G and its Lie algebra g both are realized as m × m-matrices, the dual elements can also be represented as matrices and the duality pairing could be given as ⟨p, v⟩ = trace(p T v). In this case one simply has v x = R x * ξ = ξ ⋅ x, and µ = R * x p x = p x ⋅ x T . We also note that in Euclidean space (G = R d and group operation is +), left and right translation are both realised by the identity map, e.g. R x * ξ = ξ.
We shall introduce the trivialised discrete differential of a function H as a mapdH ∶ G×G → g * satisfying the following identities generalised from (6) and (7) By identifying vectors and co-vectors in Euclidean space through the standard inner product, and noting that the logarithmic map log ∶ G → g in this setting is simply the identity map, thus log(v ⋅ u −1 ) = v − u, we recover the standard discrete gradient conditions (6) and (7) when Euclidean space is chosen as our Lie group. We also need to introduce a trivialised approximation to the bivector ω in (2). For this purpose we define, for any pair of points (u, v) ∈ G × G, an exterior 2-form on the linear space g * which we denote byω(u, v) thusω ∶ G × G → Λ 2 (g * ). We impose the consistency condition Of course, in practice,ω need only be defined in some suitable neighborhood of the diagonal subset {(x, x), x ∈ G}. Introducing coordinates, the formω plays a similar role as the skewsymmetric matrixS(x n , x n+1 ) in (10). We may further define our numerical method as follows For the reader who is unfamiliar with the notation, the definition of F (x n , x n+1 ) ∈ g using coordinates as in (4) is a skew-symmetric matrix approximating S ω (x n ) and where we have also expresseddH(x n , x n+1 ) in coordinates. From the defining relation (11), it follows immediately that the method preserves H since

Examples of trivialised discrete differentials
An example of a trivialised discrete differential, generalising (8) is Here, we have introduced a straight line in the Lie algebra between the points 0 and a = log(v⋅u −1 ). By applying exp to each point on the curve and multiplying by u, we obtain a curve on the Lie group between the points u and v, this is (ξ). Finally the (trivialisation of the) differential dH is averaged along this curve to obtain the AVF type of trivialised discrete differential. Considering (ξ) with u and v interchanged yields (1 − ξ), and from this it easily follows thatdH(u, v) = dH(v, u). The Gonzalez midpoint gradient can be generalised as well, for instance by introducing an inner product on the Lie algebra, we denote it (⋅, ⋅). We apply "index lowering" to any element η ∈ g by defining η ♭ ∈ g * to be the unique element satisfying ⟨η ♭ , ζ⟩ = (η, ζ) for all ζ ∈ g. We can then introduce a generalisation of (9) as where c ∈ G, is some point typically near u and v. One may for instance choose c = exp(η 2) ⋅ u, which implies symmetry, i.e.dH(u, v) =dH(v, u).
We have now presented two examples of trivialised discrete differentials which are both symmetric in the two arguments u and v. One observes from the definition (13)

More general manifolds
More interesting examples of mechanical systems can be found for instance in the larger class of homogeneous manifolds. There are various ways to devise discrete gradient methods in this setting, or for even more general classes of manifolds. From now on, we assume that M is a smooth manifold for which a retraction map is available, retracting the tangent bundle T M into M . This is a very basic and straightforward approach, certainly there are other ways to devise integral preserving numerical schemes. A retraction is a map Denote by φ p the restriction of φ to T p M and let 0 p be the zero-vector in T p M . Following [1], we impose the following conditions on φ This implies in particular that φ p is a diffeomorphism from some neighbourhood U of 0 p to its In what follows, we shall always assume that the step size used in the integration is sufficiently small such that both the initial and terminal point of the step are contained in such a set W . Furthermore, we shall assume that there is a given map c defined on some open subset of M ×M containing all diagonal points (p, p), for which c(p, q) ∈ M . Typically c(p, q) will be either p or q or some kind of centre point between p and q to be defined later. We shall always require that c(p, p) = p for any p ∈ M .
We assume as before the existence of a bivector ω on M and a first integral H such that the differential equation can be written in the form (2). We introduce, for any pair of points p and q on M , an approximate bivectorω(p, q) such that The discrete differential of a function H can now be defined for any pair of points (p, q) ∈ M × M as a covectordH(p, q) ∈ T * c(p,q) M satisfying the relations where c = c(p, q) is the map referred to above. We define the integrator as It follows immediately that the method is symmetric if the following three conditions are satisfied: 1. The map c is symmetric, i.e. c(p, q) = c(q, p) for all p and q.
2. The discrete differential is symmetric in the sense thatdH(p, q) =dH(q, p).
The condition 1) can be achieved by solving the equation with respect to c.
We can now write down a version of the AVF type discrete differential. Let Similarly, assuming that M is Riemannian, we can define the following counterpart to the Gonzalez midpoint discrete gradient where we may require that c(p, q) satisfies (17) for the method to be symmetric.
Example 1. We consider the sphere M = S n−1 where we represent its points as vectors in R n of unit length, p 2 = 1. The tangent space at p is then identified with the set of vectors in R n orthogonal to p with respect to the Euclidean inner product (⋅, ⋅). A retraction is its inverse is defined in the cone {q ∶ (p, q) > 0} where the geodesic midpoint between p and q in terms of the standard Riemannian metric on S n−1 . We compute the tangent map of the retraction to be As a toy problem, let us consider a mechanical system on S 2 . Since the angular momentum in body coordinates for the free rigid body is of constant length, we may assume (p, p) = 1 for all p and we can model the problem as a dynamical system on the sphere. But in addition to this, the energy of the body i preserved, i.e.
which we may take as the first integral to be preserved. Here the inertia tensor is I = diag(I 1 , I 2 , I 3 ). The system of differential equations can be written as follows where the righthand side in both equations refer to the representation in R 3 . A symmetric consistent approximation to ω would bē ω(p, q)(α, β) = ( p + q 2 , α × β) We write ξ = c + γ ξ with the notation in (18), this is a linear function of the scalar argument ξ. and thus, φ c (γ ξ ) = ξ ξ from (20). We therefore derive for the AVF discrete gradient This integral is somewhat complicated to solve analytically. Instead, we may consider the discrete gradient (19) where we take as Riemannian metric the standard Euclidean inner product restricted to the tangent bundle of S 2 . We obtain the following version of the discrete differential in the chosen representation The corresponding method is symmetric, thus of second order, and in Figure 1 we used this method to draw the trajectories of the free rigid body problem.

Methods of higher order
One viable way to obtain higher order variants of the proposed methods is by following the collocation strategy proposed by Hairer in [9] and Cohen and Hairer in [6], see also [2], and generalise it to the Lie group setting. Let c 1 , . . . , c s be distinct real numbers such that 0 ≤ c i ≤ 1 and b i ≠ 0 for all i. We consider the polynomial σ(ξh) of degree s such that whered with X j ∶= R x0 exp(σ j ). The numerical solution after one step is The collocation polynomial is obtained by integrating where j j = 1, . . . , s are the Lagrange basis functions and so

Energy preservation
To prove energy preservation of the proposed method we consider the path X(ξh) = R x0 exp(σ(ξh)) ∈ G such that X(0) = x 0 and X(h) = x 1 and integrate along this path obtaining and (26), we obtain and further so that finally Order Considering the change of variables x(ξh) = exp(Ω(ξh))x 0 ξ ∈ [0, 1], and x the solution of (2), by differentiation we get the following differential equation for Ω: Depending on the choice of quadrature points and weights, the collocation method (22), (23) and (24) approximates Ω at a certain order p, this suffices to guarantee that the overall method attains the same order, see also [6] and [14].

Examples and numerical experiments
In the numerical experiments we consider the equations of the attitude rotation of a free rigid body in unit quaternions and a problem of elasticity the equations for pseudo-rigid bodies on the cotangent bundle of GL(3) (or SL (3)).

Attitude of a free rigid body
The set with the quaternion product is a Lie group. The corresponding Lie algebra is and can be identified with R 3 . We consider the representation of the Euler attitude equations for the free rigid body in unit and I is the inertia tensor (a 3 × 3 fixed, diagonal matrix). The Euler- where I 3 is the 3 × 3 identity matrix andq is defined as follows by means of the components of q,q The energy function preserved along q(t) is We can choose the Euclidean inner product on to play the role of the Riemannian metric on S 3 , and the corresponding Riemannian gradient is with I 4 the 4 × 4 identity. It follows that the bivector ω = grad H∧F grad H 2 can be expressed by the 4 × 4 rank-2 matrix where ξ = f (q) and γ = gradH q ⋅ q c . As a discrete bivector to construct the method we choosē Identifying s 3 with its dual, and using the Gonzales trivialised discrete differential we get The energy-preserving symmetric Lie group method is Alternatively we can use the averaged trivialised discrete differential obtained averaging the Riemannian gradient as and q(ξ) = exp(ξ log(q ′ ⋅ q c )) ⋅ q.
In figure 2 we apply a symmetric energy-preserving method of order 2 and 4 to the free rigid body problem in the formulation presented in this section, and we compare them with an explicit Lie group method of order 2 based on Heun's Runge-Kutta formula. We have used the discrete trivialized differential (27) in the Lie group method of order two and the formulation outlined in section 4 for the 4-th order variant of the method. The collocation points for this method are c 1,2 = 1 2 ∓ √ 3 6 . The integrals are approximated with accurate quadrature formulae.

Pseudo-rigid bodies
We consider the Hamiltonian equations describing rotating homogeneous elastic rigid bodies [12], [11]. A pseudo-rigid body is a three dimensional elastic body whose deformation gradient F = F (t) ∈ GL + (3) is assumed to be constant throughout the body B: for all X ∈ B ⊂ R 3 , the deformation is always given by the matrix vector product F X and F is not depending on X.
In the incompressible case F ∈ SL(3). The configuration space of a pseudo-rigid body is the group G equal to GL + (3) or SL(3), and the flow of the corresponding Hamiltonian equations evolves on T * G ≈ G × g * . With the semidirect-product group structure induced by the group multiplication in G, G × g * is also a Lie group. In our particular example GL + (3) × gl(3) * , we identify gl(3) * with its dual, and we use coordinates (F, P ) where F ∈ GL + (3) and P ∈ gl(3) respectively. The corresponding Lie algebra is gl(3) × gl(3) * . We denote with W (C) the stored energy function depending on the Cauchy-Green tensor C = F T F , and with E the inertia tensor. The Hamiltonian function is where ⟨P, V ⟩ = tr(P T V ) is the standard matrix inner product, i.e. the duality pairing between between T G and T * G. The canonical Hamilton's equations take the forṁ where ( δH δF , δH δP ) = (−2F ∇W (F T F ), P E −1 ). In our experiments we consider a St Venant-Kirchhoff material leading to a stored energy function W (C) = 1 2 λ(tr(C − I)) 2 + µ tr((C − I) 2 ), with tr denoting the trace operator and λ and µ the Lamé constants (such that µ > 0 and 3λ + 2µ > 0, µ = 1 and λ = 1 3 in the experiments).
The bivector ω is represented in this case simply by the 6 × 6 inverse Darboux matrix The semi-direct product group multiplication in G × g * is and as matrix operation Ad * , we denote with Exp and Log the exponential and logarithm between G × g * and g × g * respectively, these are defined by: Exp((ξ, µ)) = (exp(ξ), dexp ξ (µ)), Log((g, σ)) = (log(g), (dexp * − log(g) ) −1 (σ)), where exp and log are the corresponding maps for G. We consider and denote with (γ 1 , γ 2 ) ∶= R * (F ,P ) dH (F ,P ) ∈ gl(3) * × gl(3) the trivialized differential in the point (F ,P ) = Exp(η F 2, η P )(F 0 , P 0 ), we have and in matrix form ad * γ2 (P ) = γ T 2P −P γ T 2 . The trivialised discrete differential (15) becomes in coordinates and where the metric is deduced by the standard matrix inner product, The energy preserving method can then be formulated as In figure 3 we report the results of a simulation for this problem. We compare an explicit Lie group method (Heun's method), a symmetric Lie group method and the symmetric energypreserving Lie group method presented in this section. The symmetric Lie group method is obtained by setting α = 0 in (30). All methods are Lie group methods of order 2. We have halved the step-size for the explicit second order Lie group method, to avoid instability. All Lie group methods have the property that det(F n ) > 0 for all n. The explicit Lie group method fails to preserve the energy of the problem ( figure 3 top left), the symmetric Lie group methods have both a much smaller energy error (figure 3 top right), the symmetric energy preserving Lie group method preserves the energy to very high precision. We report here experiments with diagonal initial values for F and P . We have performed experiments also with non diagonal initial values obtaining similar results. The performance of the methods is relying on the accurate computation of matrix functions, and in particular matrix logarithms. In figure 3 (bottom right) we show the difference in the determinants of F for the two symmetric methods: the symmetric one denoted (sym) and the symmetric and energy preserving one denoted (EP). We plot det(F sym n )−det(F EP n ) for n = 1, . . . , 8000, this measures how the two numerical solutions depart from each other with time.