Energy preserving methods on Riemannian manifolds
Elena Celledoni, Sølve Eidnes, Brynjulf Owren, Torbjørn Ringholm
EEnergy preserving methods on Riemannian manifolds
Elena Celledoni ∗ , Sølve Eidnes ∗ , Brynjulf Owren ∗ , Torbjørn Ringholm * May 22, 2018
Abstract
The energy preserving discrete gradient methods are generalized to finite-dimensionalRiemannian manifolds by definition of a discrete approximation to the Riemannian gradient,a retraction, and a coordinate center function. The resulting schemes are intrinsic anddo not depend on a particular choice of coordinates, nor on embedding of the manifoldin a Euclidean space. Generalizations of well-known discrete gradient methods, such asthe average vector field method and the Itoh–Abe method are obtained. It is shown howmethods of higher order can be constructed via a collocation-like approach. Local and globalerror bounds are derived in terms of the Riemannian distance function and the Levi-Civitaconnection. Some numerical results on spin system problems are presented.
Keywords : Geometric integration, discrete gradients, Riemannian manifolds, numericalanalysis.
Classification : 37K05, 53B99, 65L05, 82-08
A first integral of an ordinary differential equation (ODE) is a scalar-valued function on the phasespace of the ODE that is preserved along solutions. The potential benefit of using numericalmethods that preserve one or more such invariants is well-documented, and several energy-preserving methods have been developed in recent years. Among these are the discrete gradientmethods, which were introduced for use in Euclidean spaces in [1], see also [2]. These methodsare based on the idea of expressing the ODE using a skew-symmetric operator and the gradient ofthe first integral, and then creating a discrete counterpart to this in such a way that the numericalscheme preserves the energy.For manifolds in general, one can use the same schemes expressed in local coordinates. Adrawback is that the numerical approximation will typically depend on the particular choiceof coordinates and also on the strategy used for transition between coordinate charts. Another * Department of Mathematical Sciences, Norwegian University of Science and Technology, N–7491 TrondheimElena Celledoni: [email protected] ; Sølve Eidnes: [email protected] ; Brynjulf Owren: [email protected] ; Torbjørn Ringholm: [email protected]
This work was supported by the European Union’s Horizon 2020 research and innovation programme under the MarieSkłodowska-Curie grant agreement No. 691070. a r X i v : . [ m a t h . NA ] M a y lternative is to use a global embedding of the manifold into a larger Euclidean space, but then ittypically happens that the numerical solution deviates from the manifold. Even if the situationcan be amended by using projection, it may not be desirable that the computed approximationdepends on the particular embedding chosen. Crouch and Grossmann [3] and Munthe-Kaas[4, 5] introduced different ways of extending existing Runge–Kutta methods to a large class ofdifferentiable manifolds. Both these approaches are generally classified as Lie group integrators,see [6] or the more recent [7] for a survey of this class of methods. They can also both beformulated abstractly by means of a post-Lie structure which consists of a Lie algebra with aflat connection of constant torsion, see e.g. [8]. In the present paper we shall state the methodsin a slightly different context, using the notion of a Riemannian manifold. It is then natural tomake use of the Levi-Civita connection, which in contrast to the post-Lie setting is torsion-free,and which in general has a non-zero curvature. For our purposes it is also an advantage that theRiemannian metric provides an intrinsic definition of the gradient. Taking an approach morein line with this, Leimkuhler and Patrick [9] considered mechanical systems on the cotangentbundle of a Riemannian manifold and succeeded in generalising the classical leap-frog schemeto a symplectic integrator on Riemannian manifolds.Some classical numerical methods in Euclidean spaces preserve certain classes of invariants;for instance, symplectic Runge–Kutta methods preserve all quadratic invariants. This can beuseful when there is a natural way of embedding a manifold into a linear space by using con-straints that are expressed by means of such invariants. An example is the 2-sphere which canbe embedded in (cid:82) by adding the constraint that these vectors should have unit length. Theclassical midpoint rule will automatically ensure that the numerical approximations remain onthe sphere as it preserves all quadratic invariants. In general, however, the invariants preservedby these methods are expressed in terms of coordinates. Hence the preservation property of themethod may be lost under coordinate changes if the invariant is no longer quadratic. In [10], ageneralization of the discrete gradient method to differential equations on Lie groups and a broadclass of manifolds was presented. Here we develop this further by introducing a Riemannianstructure that can be used to provide an intrinsic definition of the gradient as well as a means tomeasure numerical errors.The structure of this paper is as follows: In section 2, we formulate the problem to be solvedand introduce discrete Riemannian gradient methods, as well as presenting some particularexamples with special attention to a generalization of the Itoh–Abe discrete gradient. We alsobriefly discuss the Euclidean setting as a special choice of manifold and show how the standarddiscrete gradient methods are recovered in this case. In the third section, we consider higherorder energy preserving methods based on generalization of a collocation strategy introduced byHairer [11] to Riemannian manifolds. We present some error analysis in section 4, and shownumerical results in section 5, where the methods are applied to spin system problems. Consider an initial value problem on the finite-dimensional Riemannian manifold ( M , g ) , ˙ u = F ( u ), u (0) = u ∈ M . (2.1)2e denote by F ( M ) the space of smooth functions on M . The set of smooth vector fields anddifferential one-forms are denoted Γ ( T M ) and Γ ( T ∗ M ) respectively, and for the duality pairingbetween these two spaces we use the angle brackets 〈· , ·〉 .A first integral associated to a vector field F ∈ Γ ( T M ) is a function H ∈ F ( M ) such that 〈 d H , F 〉 vanishes identically on M . First integrals are preserved along solutions of (2.1), dd t H ( u ( t )) = 〈 d H ( u ( t )), ˙ u ( t ) 〉 = 〈 d H ( u ( t )), F ( u ( t )) 〉 = The fact that a vector field F has a first integral H is closely related to the existence of a tensorfield Ω ∈ Γ ( T M ⊗ T ∗ M ) = : Γ ( T M ) , skew-symmetric with respect to the metric g , such that F ( u ) = Ω ( u ) grad H ( u ), (2.2)where grad H ∈ Γ ( T M ) is the Riemannian gradient, the unique vector field satisfying 〈 d H , ·〉 = g (grad H , · ) . Any ODE (2.1) where F is of this form preserves H , since dd t H ( u ) = 〈 d H ( u ), ˙ u 〉 = (cid:173) d H ( u ), Ω grad H ( u ) (cid:174) = g (grad H ( u ), Ω grad H ( u )) = A converse result is detailed in the following proposition.
Proposition 1.
Any system (2.1) with a first integral H can be written with an F of the form (2.2).The skew tensor field Ω can be chosen so as to be bounded near every nondegenerate criticalpoint of H .Proof. Similar to the proof of Proposition in [2], we can write an explicit expression for apossible choice of Ω , Ω y = g (grad H , y ) F − g ( F , y ) grad Hg (grad H , grad H ) . (2.3)Clearly, g ( y , Ω y ) = for all y . Since H is a first integral, g ( F , grad H ) = 〈 d H , F 〉 = , so Ω grad H = F . For a proof that Ω is bounded near nondegenerate critical points, see [2].In fact, such a tensor field Ω often arises naturally from a two-form ω through Ω y = ω ( · , y ) (cid:93) . Awell-known example is when ω is a symplectic two-form. Note that Ω is not necessarily unique.Retractions, viewed as maps from T M to M , will play an important role in the methods wediscuss here. Their formal definition can be found e.g. in [12]: Definition 1.
Let φ be a smooth map φ : T M → M and let φ p denote the restriction of φ to T p M ,with p being the zero-vector in T p M . Then φ is a retraction if it satisfies the conditions1. φ p is defined in an open ball B r p (0 p ) ⊂ T p M of radius r p about p ,2. φ p ( x ) = p if and only if x = p ,3. T φ p | p = Id T p M .A generic example of a retraction on ( M , g ) is obtained via the Riemannian exponential,setting φ p ( x ) = exp p ( x ) , i.e. following along the geodesic emanating from p in the direction x .3 .2 The discrete Riemannian gradient method We adapt the discrete gradients in Euclidean space to discrete Riemannian gradients (DRG) on ( M , g ) by means of a retraction map φ and a center point function c . Definition 2.
A discrete Riemannian gradient is a triple (grad, φ , c ) where1. c : M × M → M is a continuous map such that c ( u , u ) = u for all u ∈ M ,2. grad : F ( M ) → Γ ( c ∗ T M ) ,3. φ : T M → M is a retraction,such that for all H ∈ F ( M ) , u ∈ M , v ∈ M , c = c ( u , v ) ∈ M , H ( v ) − H ( u ) = g (grad H ( u , v ), φ − c ( v ) − φ − c ( u )), (2.4) grad H ( u , u ) = grad H ( u ). (2.5)The DRG grad H is a continuous section of the pullback bundle c ∗ T M , meaning that π ◦ grad H = c , where π : T M → M is the natural projection. We also need to define an ap-proximation to be used for the tensor field Ω ∈ Γ ( T M ) . To this end we let Ω ∈ Γ ( c ∗ T M ) be acontinuous skew-symmetric tensor field such that Ω ( u , u ) = Ω ( u ) ∀ u ∈ M . Inspired by [10, 13], we propose the scheme u k + = φ c k ( W ( u k , u k + )), c k = c ( u k , u k + ) (2.6) W ( u k , u k + ) = φ − c k ( u k ) + h Ω ( u k , u k + ) grad H ( u k , u k + ), (2.7)where h is the step size. The scheme (2.6)–(2.7) preserves the invariant H , since H ( u k + ) − H ( u k ) = g (grad H ( u k , u k + ), φ − c k ( u k + ) − φ − c k ( u k )) = g (grad H ( u k , u k + ), h Ω ( u k , u k + ) grad H ( u k , u k + )) = Here and in the following we adopt the shorthand notation c = c ( u , v ) as long as it is obviouswhat the arguments of c are.The Average Vector Field (AVF) method has been studied extensively in the literature; someearly references are [14, 2, 15]. This is a discrete gradient method, and we propose a correspond-ing DRG satisfying (2.4)-(2.5) as follows: grad AVF H ( u , v ) = (cid:90) ( T γ ξ φ c ) T grad H ( φ c ( γ ξ )) d ξ , γ ξ = (1 − ξ ) φ − c ( u ) + ξφ − c ( v ), (2.8)where ( T x φ c ) T : T φ c ( x ) M → T x M is the unique operator satisfying g (( T x φ c ) T a , b ) = g ( a , T x φ c b ), ∀ x , b ∈ T c M , a ∈ T φ c ( x ) M . To avoid cluttered notation we will just write grad for the triple (grad, φ , c ) in the sequel. grad MP H ( u , v ) = grad H ( c ( u , v )) + H ( v ) − H ( u ) − g (grad H ( c ( u , v )), η ) g ( η , η ) η (2.9)where η = φ − c ( v ) − φ − c ( u ) .Note that both these DRGs involve the gradient of the first integral. This may be a disad-vantage if H is non-smooth or if its gradient is expensive to compute. Also, the implicit natureof the schemes requires the solution of an n -dimensional nonlinear system of equations at eachtime step. An alternative is to consider the Itoh–Abe discrete gradient [16], also called the co-ordinate increment discrete gradient [2], which in certain cases requires only the solution of n decoupled scalar equations. We now present a generalization of the Itoh–Abe discrete gradientto finite-dimensional Riemannian manifolds. Definition 3.
For any tangent space T c M one can choose a basis { E , ..., E n } composed of tangentvectors E i , i =
1, ..., n , orthonormal with respect to the Riemannian metric g . Then, given u , v ∈ M , there exists a unique { α i } ni = so that φ − c ( v ) − φ − c ( u ) = n (cid:88) i = α i E i . The
Itoh–Abe DRG of the first integral H is then given by grad IA H ( u , v ) = n (cid:88) j = a j E j , (2.10)where a j = H ( w j ) − H ( w j − ) α j if α j (cid:54)= g (grad H ( w j − ), T φ c ( η j − ) E j ) if α j = w j = φ c ( η j ), η j = φ − c ( u ) + j (cid:88) i = α i E i . We refer to [13] for proof that this is indeed a DRG satisfying (2.4)-(2.5).
Let M = V be an (cid:82) -linear space, and let g be the Euclidean inner product, g ( x , y ) = x T y . Theoperator Ω is a solution dependent skew-symmetric n × n matrix Ω ( u ) . For any u ∈ V , we have T u V ≡ V . The retraction φ : V → V is defined as φ p ( x ) = p + x , the Riemannian exponential on V , so that φ − c ( v ) − φ − c ( u ) = v − u . The gradient grad H is an n -vector whose i th component is5 H ∂ u i , and the definition of the discrete Riemannian gradient coincides with the standard discretegradient, since (2.4) now reads H ( v ) − H ( u ) = grad( u , v ) T ( v − u ). Furthermore, (2.6)-(2.7) simply becomes the discrete gradient method introduced in [1], givenby the scheme u k + − u k = h Ω ( u k , u k + ) grad H ( u k , u k + ), (2.11)where Ω is a skew-symmetric matrix approximating Ω . Typical choices are Ω ( u k , u k + ) = Ω ( u k ) ,or Ω ( u k , u k + ) = Ω (( u k + + u k )/2) if one seeks a symmetric method.The DRGs (2.8) and (2.9) become the standard AVF and midpoint discrete gradients in thiscase. For the Itoh–Abe DRG, the practical choice for the orthogonal basis would be the set ofunit vectors, { e , ..., e n } , so that α i = v i − u i , and we get (2.10) with a j = H ( w j ) − H ( w j − ) v j − u j if u j (cid:54)= v j , ∂ H ∂ u j ( w j − ) if u j = v j , w j = (cid:88) ji = v i e i + (cid:88) ni = j + u i e i , which is a reformulation of the Itoh–Abe discrete gradient as it is given in [16], [2] and theliterature otherwise. In the Euclidean setting, a strategy to obtain energy preserving methods of higher order waspresented in [17] and later in [11], see also [18]. This technique is generalized to a Lie groupsetting in [10]. We will here formulate these methods in the context of Riemannian manifolds.
Let c , ..., c s be distinct real numbers. Consider the Lagrange basis polynomials, l i ( ξ ) = s (cid:89) j = j (cid:54)= i ξ − c j c i − c j , and let b i : = (cid:90) l i ( ξ ) d ξ . (3.1)We assume that c , . . . , c s are such that b i (cid:54)= for all i . A step of the energy-preserving collocation-like method, starting at u ∈ M , is defined via a polynomial σ : (cid:82) → T c M of degree s satisfying σ (0) = φ − c ( u ), (3.2) dd ξ σ ( ξ h ) (cid:175)(cid:175)(cid:175) ξ = c j = T U j φ − c (cid:179) Ω j grad j H (cid:180) , U j : = φ c (cid:161) σ ( c j h ) (cid:162) (3.3) u : = φ c ( σ ( h )) , (3.4)6here grad j H : = (cid:90) l j ( ξ ) b j (cid:161) T U j φ − c (cid:162) T ( T σ ( ξ h ) φ c ) T grad H (cid:161) φ c ( σ ( ξ h )) (cid:162) d ξ , and Ω j : = Ω ( U j ). Notice that with s = and independently on the choice of c , we reproduce the DRG method(2.6)-(2.7) with the AVF DRG (2.8).Using Lagrange interpolation and (3.3), the derivative of σ ( ξ h ) at every point ξ h is dd ξ σ ( ξ h ) = s (cid:88) j = l j ( ξ ) T U j φ − c (cid:179) Ω j grad j H (cid:180) , (3.5)from which by integrating we get σ ( τ h ) = φ − c ( u ) + h s (cid:88) j = (cid:90) τ l j ( ξ ) d ξ T U j φ − c (cid:179) Ω j grad j H (cid:180) . The defined method is energy preserving, which we see by using dd ξ (cid:161) φ c ( σ ( ξ h )) (cid:162) = T σ ( ξ h ) φ c (cid:181) dd ξ σ ( ξ h ) (cid:182) , and (3.5) to get H ( u ) − H ( u ) = (cid:90) g (cid:181) grad H (cid:161) φ c ( σ ( ξ h )) (cid:162) , dd ξ φ c ( σ ( ξ h )) (cid:182) d ξ = (cid:90) g (cid:195) grad H (cid:161) φ c ( σ ( ξ h )) (cid:162) , T σ ( ξ h ) φ c (cid:195) s (cid:88) j = l j ( ξ ) T U j φ − c (cid:179) Ω j grad j H (cid:180)(cid:33)(cid:33) d ξ = (cid:90) g (cid:195)(cid:161) T σ ( ξ h ) φ c (cid:162) T grad H (cid:161) φ c ( σ ( ξ h )) (cid:162) , s (cid:88) j = l j ( ξ ) T U j φ − c (cid:179) Ω j grad j H (cid:180)(cid:33) d ξ = s (cid:88) j = b j g (cid:181)(cid:90) l j ( ξ ) b j (cid:161) T U j φ − c (cid:162) T (cid:161) T σ ( ξ h ) φ c (cid:162) T grad H (cid:161) φ c ( σ ( ξ h )) (cid:162) d ξ , Ω j grad j H (cid:182) = s (cid:88) j = b j g (cid:179) grad j H , Ω j grad j H (cid:180) = and hence repeated use of (3.2)-(3.4) ensures H ( u k ) = H ( u ) for all k ∈ (cid:78) . From the Itoh–Abe DRG one can get a new DRG, also satisfying (2.4), by grad
SIA H ( u , v ) = (cid:179) grad IA H ( u , v ) + grad IA H ( v , u ) (cid:180) . (3.6)We call this the symmetrized Itoh–Abe DRG . Note that we need the base point c to be the samein the evaluation of grad IA H ( u , v ) and grad IA H ( v , u ) . When c ( u , v ) = c ( v , u ) and Ω ( u , v ) = Ω ( v , u ) ,we get a symmetric DRG method (2.6)-(2.7), which is therefore of second order.7lternatively, one can get a symmetric -stage method by a composition of the Itoh–AbeDRG method and its adjoint. Furthermore, one can get energy preserving methods of any orderusing a composition strategy. To ensure symmetry of an s -stage composition method, oneneeds c i ( u , v ) = c s + − i ( v , u ) for different center points c i belonging to each stage and, similarly, Ω i ( u , v ) = Ω s + − i ( v , u ) . In this section, ϕ t ( u ) is the t -flow of the ODE vector field F . The most standard discrete gradientmethods have a low or moderate order of convergence, and that is also the case for the DRGmethods unless special care is taken in designing Ω and grad H . We shall not pursue this approachhere, but refer to the collocation-like methods if high order of accuracy is required. We shall see,however, that the methods designed here are consistent and can be made symmetric. Analysisof the local error can be done in local coordinates, assuming that the step size is always chosensufficiently small, so that within a fixed step, u k , u k + , c ( u k , u k + ) and the exact local solution u ( t k + ) all belong to the same given coordinate chart. From the definition (2.6)-(2.7) it followsimmediately that the representation of u k + ( h ) satisfies u k + (0) = u k and ddh u k + (0) = F ( u k ) .Then by equivalence of local coordinate norms and the Riemannian distance, we may concludethat the local error in DRG methods satisfies d ( u k + , ϕ h ( u k )) ≤ C h . Similar to what was also observed in [10], the DRG methods (2.6)-(2.7) are symmetric whenever grad H ( u , v ) = grad H ( v , u ) , Ω ( u , v ) = Ω ( v , u ) , and c ( u , v ) = c ( v , u ) for all u , v ∈ M . In that casewe obtain an error bound for the local error of the form d ( u k + , ϕ h ( u k )) ≤ C h .The collocation-like methods of section 3 have associated nodes { c i } si = and weights { b i } si = defined by (3.1). The order of the local error depends on the accuracy of the underlying quadratureformula given by these nodes and weights. The following result is a simple consequence ofTheorem 4.3 in [18]. Theorem 1.
Let ψ h be the method defined by (3.2) - (3.4) . The order of the local error is at least p = min( r , 2 r − s + where r is the largest integer such that (cid:80) si = b i c q − i = q for all ≤ q ≤ r . This means that thereare positive constants C and h such that d ( ψ h ( u ), ϕ h ( u )) ≤ C h p + for h < h , u ∈ M . Proof.
Choose h small enough such that the solution can be represented in the form u ( h ξ ) = φ c ( γ ( ξ h )), ξ ∈ [0, 1], and consider the corresponding differential equation for γ in T c M : dd t γ ( t ) = (cid:161) φ ∗ c F (cid:162) ( γ ( t )) = (cid:161) T γ ( t ) φ c (cid:162) − Ω grad H (cid:161) φ c ( γ ( t )) (cid:162) . (4.1)8otice that (cid:161) T γ φ c (cid:162) − = T U φ − c where U = φ c ◦ γ and T U ( t ) φ − c : T U ( t ) M → T c M for every t . Weobtain dd t γ ( t ) = T U ( t ) φ − c Ω (cid:161) T U ( t ) φ − c (cid:162) T (cid:161) T γ ( t ) φ c (cid:162) T grad H (cid:161) φ c ( γ ( t )) (cid:162) . (4.2)Considering the Hamiltonian (cid:101) H : T c M → (cid:82) , (cid:101) H ( γ ) : = φ ∗ c H ( γ ) = H ◦ φ c ( γ ) , we can then rewrite(4.1) in the form dd t γ ( t ) = (cid:101) Ω ( γ ) grad (cid:101) H ( γ ), (cid:101) Ω ( γ ) : = T U ( t ) φ − c Ω (cid:161) T U ( t ) φ − c (cid:162) T , (4.3)where we have used that grad (cid:101) H = T γ ( t ) φ T c grad H ( φ c ( γ ( t ))) , which is now a gradient on the linearspace T c M with respect to the metric inherited from M , g | c . Locally in a neighborhood of c ,(3.2)-(3.4) applied to (4.3) coincides with the methods of Cohen and Hairer, and therefore theorder result [18, Thm 4.3] can be applied. Since the Riemannian distance d ( · , · ) and any norm inlocal coordinates are equivalent, the result follows. We prove the following result for the global error in DRG methods.
Theorem 2.
Let u ( t ) be the exact solution to (2.1) where F is a complete vector field on aconnected Riemannian manifold ( M , g ) with flow u ( t ) = ϕ t ( u ) . Let ψ h represent a numericalmethod u k + = ψ h ( u k ) whose local error can be bounded for some p ∈ (cid:78) as d ( ψ h ( u ), ϕ h ( u )) ≤ C h p + for all u ∈ M . Suppose there is a constant L such that (cid:107)∇ F (cid:107) g ≤ L , where ∇ is the Levi-Civita connection and (cid:107) · (cid:107) g is the operator norm with respect to the metric g . Then the global error is bounded as d (cid:179) u ( kh ), u k (cid:180) ≤ CL (e khL − h p for all k > Proof.
Denoting the global error as e k : = d ( u ( kh ), u k ), the triangle inequality yields e k + ≤ d (cid:161) ϕ h ( u ( kh ) (cid:162) , ϕ h ( u k )) + d (cid:179) ϕ h ( u k ), ψ h ( u k ) (cid:180) . The first term is the error at nh propagated over one step, the second term is the local error. Forthe first term, we find via a Grönwall type inequality of [19], d (cid:179) ϕ h ( u ( kh )), ϕ h ( u k ) (cid:180) ≤ e hL d (cid:179) u ( kh ), u k (cid:180) = e hL e k . Using the local error estimate for the second term, we get the recursion e k + ≤ e hL e k + C h p + , e k ≤ C e khL − hL − h p + ≤ CL (e khL − h p . Remark:
Following Theorem 1.4 in [19], the condition that F is complete can be relaxed if ϕ t ( u ) and { u k } k ∈ (cid:78) lie in a relatively compact submanifold N of M containing all the geodesicsfrom u k to ϕ kh ( u ) . This is the case if, for instance, H has compact, geodesically convex sublevelsets, since both ϕ t ( u ) and { u k } k ∈ (cid:78) are restricted to the level set M H ( u ) = { p ∈ M | H ( p ) = H ( u )} and hence lie in the sublevel set N H ( u ) = { p ∈ M | H ( p ) ≤ H ( u )} . We test our methods on two different variants of the classical spin system, whose solution evolveson the d -fold product of two-spheres, ( S ) d , d s i d t = s i × ∂ H ∂ s i , s i ∈ S , i =
1, . . . , d , H ∈ F (cid:179)(cid:161) S (cid:162) d (cid:180) . (5.1)The Riemannian metric g on ( S ) d restricts to the so-called round metric on each copy of thesphere. This metric coincides with the Euclidean inner product on the tangent planes of each ofthe spheres.Geometric integrators for such systems are discussed widely in the literature, see e.g. [20, 21,22, 23] and references therein. We study one or more bodies whose orientation is represented bya vector s i of unit length in (cid:82) , so that s i lies on the manifold M = S = (cid:169) s ∈ (cid:82) : (cid:107) s (cid:107) = (cid:170) . Hereand in what follows, (cid:107)·(cid:107) denotes the -norm. Starting with d = , our choice of retraction φ isgiven by its restriction to p , φ p ( x ) = p + x (cid:176)(cid:176) p + x (cid:176)(cid:176) , (5.2)with the inverse φ − p ( u ) = up T u − p defined when p T u > . We note that p T x = for all x ∈ T p S . The tangent map of the retractionand its inverse are given by T x φ p = (cid:176)(cid:176) p + x (cid:176)(cid:176) (cid:195) I − ( p + x ) ⊗ ( p + x ) (cid:176)(cid:176) p + x (cid:176)(cid:176) (cid:33) , T u φ − p = p T u (cid:181) I − u ⊗ pp T u (cid:182) , (5.3)where ⊗ denotes the outer product of the vectors. For d > , we use the retraction defined by Φ p ( x ) = ( φ p ( x i ), . . . , φ p d ( x d )) , where each φ p i ( x i ) is given by (5.2). If x and y are in (cid:82) , x ⊗ y is the matrix-matrix product of x taken as a × matrix and y taken as a × matrix. .1 Example 1: Perturbed spinning top We consider first a nonlinear perturbation of a spinning top, see [22]. This is a spin system withone spin s . Given the inertia tensor (cid:73) = diag ( (cid:73) , (cid:73) , (cid:73) ) , and denoting by s the component-wisesquare of s , we can define the Hamiltonian as H ( s ) =
12 ( (cid:73) − s ) T ( s + s ). The ODE system can be written in the form d s d t = Ω ( s ) grad H ( s ), Ω ( s ) = ˆ s , using the hat operator defined by ˆ s y = s × y . We approximate this system numerically, testingthe scheme (2.6)-(2.7) with different discrete Riemannian gradients: the AVF (2.8), the midpoint(2.9), the Itoh–Abe (2.10) and its symmetrized version (3.6). For the three symmetric methods,we have chosen c ( s , (cid:101) s ) = s + (cid:101) s (cid:107) s + (cid:101) s (cid:107) , so that φ − c ( (cid:101) s ) = − φ − c ( s ) . Using that grad H ( s ) = (cid:73) − ( s + s ) andconsidering the transpose of T γ ξ φ c from (5.3), the AVF DRG becomes grad AVF H ( s , (cid:101) s ) = (cid:90) (cid:107) l ξ (cid:107) (cid:181) I − l ξ ⊗ l ξ (cid:107) l ξ (cid:107) (cid:182) (cid:73) − ( φ c ( γ ξ ) + φ c ( γ ξ ) ) d ξ = (cid:90) (cid:107) l ξ (cid:107) (cid:161) (cid:73) − (cid:161) φ c ( γ ξ ) + φ c ( γ ξ ) (cid:162) − φ c ( γ ξ ) T (cid:73) − (cid:161) φ c ( γ ξ ) + φ c ( γ ξ ) (cid:162) φ c ( γ ξ ) (cid:162) d ξ , with γ ξ = (1 − ξ ) φ − c ( s ) + ξφ − c ( (cid:101) s ) = (1 − ξ ) φ − c ( s ) and l ξ = c + γ ξ . Similarly, the midpoint DRGbecomes grad MP H ( s , (cid:101) s ) = (cid:107) s + (cid:101) s (cid:107) (cid:195) (cid:73) − (cid:181) s + (cid:101) s + (cid:161) s + s (cid:101) s + (cid:101) s (cid:162)(cid:182) + (cid:107) s + (cid:101) s (cid:107) − (cid:107) (cid:101) s − s (cid:107) ( H ( (cid:101) s ) − H ( s )) ( (cid:101) s − s ) (cid:33) , where we have used that g ( s , s ) = s T s = for all s ∈ S . To obtain the basis of T c M for thedefinition of the Itoh–Abe DRG, we have used the singular-value decomposition. For the firstorder scheme, noting that φ − s ( s ) = , we choose c ( s , (cid:101) s ) = s , and get α j = φ s ( (cid:101) s ) T E j , for j =
1, 2 .Then the DRG (2.10) can be written as grad IA H ( s , (cid:101) s ) = H (cid:161) φ s (cid:161) φ − s ( (cid:101) s ) T E E (cid:162)(cid:162) − H ( s ) φ − s ( (cid:101) s ) T E E + H ( (cid:101) s ) − H (cid:161) φ s (cid:161) φ − s ( (cid:101) s ) T E E (cid:162)(cid:162) φ − s ( (cid:101) s ) T E E . (5.4)We solve the same problem using the 4th, 6th and 8th order variants of the collocation-likescheme (3.2)-(3.4). Choosing in the 4th order case the Gaussian nodes c = ∓ (cid:112) as collocationpoints and setting c ( s , (cid:101) s ) = s , we get the nonlinear system S = h φ s (cid:195) T S φ − s (cid:161) Ω grad H (cid:162) + (cid:195) − (cid:112) (cid:33) T S φ − s (cid:161) Ω grad H (cid:162)(cid:33) , S = h φ s (cid:195)(cid:195) + (cid:112) (cid:33) T S φ − s (cid:161) Ω grad H (cid:162) + T S φ − s (cid:161) Ω grad H (cid:162)(cid:33) , s = h φ s (cid:161) T S φ − s (cid:161) Ω grad H (cid:162) + T S φ − s (cid:161) Ω grad H (cid:162)(cid:162) , -3 -2 -1 Step size -10 -8 -6 -4 -2 N o r m o f e rr o r AVFMPIASIAComp-2 -1 Step size -15 -10 -5 N o r m o f e rr o r Coll-4Coll-6Coll-8Comp-SIAComp-4
Figure 1: Error norm at t = for the perturbed spinning top problem solved with differentschemes, plotted with black, dashed reference lines of order 1, 2, 4, 6 and 8. Initial condition s = ( − −
1, 1)/ (cid:112) and (cid:73) = diag (1, 2, 4) . Left:
The AVF, midpoint (MP), Itoh–Abe (IA) and sym-metrized Itoh–Abe (SIA) DRGs and a 3-stage composition of the IA DRG scheme (Comp-2).
Right:
Collocation-type schemes of order 4, 6 and 8, a 3-stage composition of the SIA DRGscheme (Comp-SIA), and a 6-stage composition of the IA DRG scheme (Comp-4).where σ ( ξ h ) = (cid:179)(cid:179) + (cid:112) (cid:180) φ − s ( S ) + (cid:179) − (cid:112) (cid:180) φ − s ( S ) (cid:180) ξ + (cid:179) (cid:179) (cid:112) − (cid:180) φ − s ( S ) − (cid:179) + (cid:112) (cid:180) φ − s ( S ) (cid:180) ξ and we use the transposes of (5.3) and grad H ( s ) = (cid:73) − ( s + s ) in the evaluation of grad H and grad H . The 6th and 8th order schemes are derived in a similar manner, using the standardGaussian nodes.A second order scheme is derived by composing the Itoh–Abe DRG method with its adjoint,and a 4th order scheme is obtained by composing this method again with itself, as well as oneby composition of the symmetrized Itoh–Abe DRG method with itself. In all stages of thesecomposition methods, a symmetric c ( u , v ) is used.Plots confirming the order of all methods can be seen in Figure 1, where solutions using thedifferent schemes are compared to a reference solution obtained using a very small step size. Seethe left hand panel of Figure 2 for numerical confirmation that our methods do indeed preservethe energy to machine precision, while the implicit midpoint method does not. In the right handpanel of Figure 2, the solution obtained by the Itoh–Abe DRG scheme with a step size h = isplotted together with a solution obtained using the symmetrized Itoh–Abe DRG method with amuch smaller time step. We observe, as expected for a method that conserves both the energy andthe angular momentum, that the solution stays on the trajectories of the exact solution, althoughnot necessarily at the right place on the trajectory at any given time.12
200 400 600 800 1000
Time -15 -10 -5 E ne r g y e rr o r AVFMPIAIMP -1-1 -1-0.500.5 001 11
Figure 2:
Left:
Energy error with increasing time for the AVF, midpoint (MP) and Itoh–Abe(IA) DRG methods, as well as the implicit midpoint (IMP) method, with step size h = , initialcondition s = ( − −
1, 1)/ (cid:112) and (cid:73) = diag (1, 2, 4) . Right:
Curves of constant energy on the sphere,found by our method with different starting values. The black solid line is the solution using thesymmetrized Itoh–Abe DRG method with step size h = , while the red dots are the solutionsobtained by the Itoh–Abe DRG method with step size h = . We now consider the Heisenberg spin chain of micromagnetics. This problem is considered in[20, 23], where different geometric integrators are tested. Here, s ∈ (cid:161) S (cid:162) d , and the Hamiltonian is H ( s ) = d (cid:88) i = s i T s i − , (5.5)with s = s d and s d + = s . The system (5.1) becomes, for this Hamiltonian, d s i d t = ˆ s i ( s i − + s i + ) , i =
1, . . . , d , and can be written in the block form d s d t = Ω ( s ) grad H ( s ), where Ω ( s ) = diag( ˆ s , . . . , ˆ s d ). (5.6)For such a d -particle system, we may write the DRGs as grad H ( s , (cid:101) s ) = (cid:179) grad H ( s , (cid:101) s ), . . . , grad d H ( s , (cid:101) s ) (cid:180) , where we note that grad i H ( s , (cid:101) s ) is a discrete approximation to ∂ H ∂ s i . We thus get the AVF DRG13efined by grad i AVF H ( s , (cid:101) s ) = (cid:90) (cid:107) l i , ξ (cid:107) (cid:181) I − l i , ξ ⊗ l i , ξ (cid:107) l i , ξ (cid:107) (cid:182) (cid:161) φ c i − ( γ i − ξ ) + φ c i + ( γ i + ξ ) (cid:162) d ξ = (cid:90) (cid:107) l i , ξ (cid:107) (cid:179) φ c i − ( γ i − ξ ) + φ c i + ( γ i + ξ ) − l T i , ξ (cid:161) φ c i − ( γ i − ξ ) + φ c i + ( γ i + ξ ) (cid:162) l i , ξ (cid:180) d ξ , with γ i , ξ = (1 − ξ ) φ − c i ( s i ) and l i , ξ = c i + γ i , ξ . For the midpoint DRG we get grad i MP H ( s , (cid:101) s ) = c i − + c i + + H ( (cid:101) s ) − H ( s ) − (grad H ( c ( s , (cid:101) s )) T ηη T η η i , where η = ( η , . . . , η d ) and η i = − φ − c i ( s i ). In the numerical experiments, however, we have useda small modification of this, grad i MMP H ( s , (cid:101) s ) = c i − + c i + + (cid:101) s T i (cid:101) s i − − s i T s i − − ( c i − + c i + ) T η i η i T η i η i . This DRG, which does indeed satisfy (2.4)-(2.5), leads to a more computationally efficientscheme than the original midpoint DRG. Each grad i IA H ( s , (cid:101) s ) in the Itoh–Abe DRG is found asin the previous example, by (5.4). Higher order schemes are also derived in the same manner asbefore.We test our schemes by comparing the numerical solutions with the exact solution s j ( t ) = ( a cos θ j + (cid:101) a sin θ j ) cos φ + ¯ a sin φ , θ j = j p − − cos p ) sin φ , for a choice of constants φ , p ∈ (cid:82) and orthogonal unit vectors a , (cid:101) a , ¯ a ∈ (cid:82) , see [20]. Order plotsfor the methods are provided in Figure 3, using d = , φ = π /3 , p = π / d , a = (1, 2, − (cid:112) , (cid:101) a = (2, 1, 4)/ (cid:112) and ¯ a = a × (cid:101) a . All schemes are shown to have the expected order.
We have presented a general framework for constructing energy preserving numerical integratorson Riemannian manifolds. The main tool is to generalize the notion of discrete gradients asknown from the literature. The new methods make use of an approximation to the Riemanniangradient coined the discrete Riemannian gradient, as well as a retraction map and a coordinatecenter function. An appealing feature of the new methods is that they do not depend on aparticular choice of local coordinates or on an embedding of the manifold into a (larger) Euclideanspace, but are of an intrinsic nature. Particular examples of discrete Riemannian gradient methodsare given as generalizations of well-known schemes, such as the average vector field method,the midpoint discrete gradient method and the Itoh–Abe method. Extensions to higher order areproposed via a collocation-like method. We have analysed the local and global error behaviourof the methods, and they have been implemented and tested for certain spin systems where thephase space is (cid:161) S (cid:162) d . 14 -3 -2 -1 Step size -6 -4 -2 N o r m o f e rr o r AVFMMPIASIAComp-2 -1 Step size -15 -10 -5 N o r m o f e rr o r Coll-4Coll-6Coll-8Comp-SIAComp-4
Figure 3: Error norm at t = for the Heisenberg spin chain problem solved with differentschemes, plotted with black, dashed reference lines of order 1, 2, 4, 6 and 8. Left:
The AVF,modified midpoint (MMP), Itoh–Abe (IA) and symmetrized Itoh–Abe (SIA) DRGs and a 3-stagecomposition of the IA DRG scheme (Comp-2).
Right:
Collocation-type schemes of order 4, 6and 8, a 3-stage composition of the SIA DRG scheme (Comp-SIA), and a 6-stage compositionof the IA DRG scheme (Comp-4).Possible directions for future research include a more detailed study of the stability andpropagation of errors, taking into account particular features of the Riemannian manifold; forinstance, it may be expected that the sectional curvature will play an important role. Moreexamples should also be tried out, and we belive, inspired by [13], that there is a potential formaking our implementations more efficient by tailoring them for the particular manifold, as wellas the ODE problem considered.
References [1] O. Gonzalez, “Time integration and discrete Hamiltonian systems,”
J. Nonlinear Sci. , vol. 6,no. 5, pp. 449–467, 1996.[2] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux, “Geometric integration using discretegradients,”
R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. , vol. 357, no. 1754,pp. 1021–1045, 1999.[3] P. E. Crouch and R. Grossman, “Numerical integration of ordinary differential equationson manifolds,”
Journal of Nonlinear Science , vol. 3, no. 1, pp. 1–33, 1993.[4] H. Munthe-Kaas, “Lie–Butcher theory for Runge–Kutta methods,”
BIT Numerical Mathe-matics , vol. 35, no. 4, pp. 572–587, 1995.155] H. Munthe-Kaas, “Runge–Kutta methods on Lie groups,”
BIT Numerical Mathematics ,vol. 38, no. 1, pp. 92–111, 1998.[6] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna, “Lie-group methods,”
ActaNumerica , vol. 9, pp. 215–365, 2000.[7] E. Celledoni, H. Marthinsen, and B. Owren, “An introduction to Lie group integrators –basics, new developments and applications,”
Journal of Computational Physics , vol. 257,pp. 1040–1061, 2014.[8] H. Z. Munthe-Kaas and A. Lundervold, “On post-Lie algebras, Lie–Butcher series andmoving frames,”
Found. Comput. Math. , vol. 13, no. 4, pp. 583–613, 2013.[9] B. Leimkuhler and G. W. Patrick, “A symplectic integrator for Riemannian manifolds,”
Journal of Nonlinear Science , vol. 6, no. 4, pp. 367–384, 1996.[10] E. Celledoni and B. Owren, “Preserving first integrals with symmetric Lie group methods,”
Discrete Contin. Dyn. Syst. , vol. 34, no. 3, pp. 977–990, 2014.[11] E. Hairer, “Energy-preserving variant of collocation methods,”
JNAIAM. J. Numer. Anal.Ind. Appl. Math. , vol. 5, no. 1-2, pp. 73–84, 2010.[12] R. L. Adler, J.-P. Dedieu, J. Y. Margulies, M. Martens, and M. Shub, “Newton’s methodon Riemannian manifolds and a geometric model for the human spine,”
IMA Journal ofNumerical Analysis , vol. 22, no. 3, pp. 359–390, 2002.[13] E. Celledoni, S. Eidnes, B. Owren, and T. Ringholm, “Dissipative schemes on Riemannianmanifolds,” arXiv preprint, arXiv:1804.08104 , 2018.[14] A. Harten, P. D. Lax, and B. van Leer, “On upstream differencing and Godunov-typeschemes for hyperbolic conservation laws,”
SIAM Rev. , vol. 25, no. 1, pp. 35–61, 1983.[15] G. Quispel and D. McLaren, “A new class of energy-preserving numerical integrationmethods,”
J. of Phys. A: Math. and Theor. , vol. 41, no. 4, pp. 045206, 7, 2008.[16] T. Itoh and K. Abe, “Hamiltonian-conserving discrete canonical equations based on varia-tional difference quotients,”
Journal of Computational Physics , vol. 76, no. 1, pp. 85–102,1988.[17] L. Brugnano, F. Iavernaro, and D. Trigiante, “Hamiltonian boundary value methods (energypreserving discrete line integral methods),”
J. Numer. Anal. Ind. Appl. Math , vol. 5, no. 1,pp. 17–37, 2010.[18] D. Cohen and E. Hairer, “Linear energy-preserving integrators for Poisson systems,”
BITNumerical Mathematics , vol. 51, no. 1, pp. 91–101, 2011.[19] M. Kunzinger, H. Schichl, R. Steinbauer, and J. A. Vickers, “Global Gronwall estimates forintegral curves on Riemannian manifolds,”
Rev. Mat. Complut. , vol. 19, no. 1, pp. 133–137,2006. 1620] J. Frank, W. Huang, and B. Leimkuhler, “Geometric integrators for classical spin systems,”
Journal of Computational Physics , vol. 133, no. 1, pp. 160–172, 1997.[21] D. Lewis and N. Nigam, “Geometric integration on spheres and some interesting applica-tions,”
Journal of Computational and Applied Mathematics , vol. 151, no. 1, pp. 141–170,2003.[22] R. I. McLachlan, K. Modin, and O. Verdier, “Symplectic integrators for spin systems,”
Physical Review E , vol. 89, no. 6, p. 061301, 2014.[23] R. McLachlan, K. Modin, and O. Verdier, “A minimal-variable symplectic integrator onspheres,”