Construction and analysis of higher order Galerkin variational integrators
CConstruction and analysis of higher order Galerkin variationalintegrators
Sina Ober-Bl¨obaum ∗ and Nils Saake † Computational Dynamics and Optimal Control, Department of Mathematics , University of Paderborn, Warburger Str. 100, 33098 Paderborn, Germany
31 October 2018
Abstract
In this work we derive and analyze variational integrators of higher order for the structure-preservingsimulation of mechanical systems. The construction is based on a space of polynomials together withGauss and Lobatto quadrature rules to approximate the relevant integrals in the variational principle.The use of higher order schemes increases the accuracy of the discrete solution and thereby decreasethe computational cost while the preservation properties of the scheme are still guaranteed. The orderof convergence of the resulting variational integrators are investigated numerically and it is discussedwhich combination of space of polynomials and quadrature rules provide optimal convergence rates.For particular integrators the order can be increased compared to the Galerkin variational integratorspreviously introduced in [MW01]. Furthermore, linear stability properties, time reversibility, structure-preserving properties as well as efficiency for the constructed variational integrators are investigated anddemonstrated by numerical examples. discrete variational mechanics; numerical convergence analysis;symplectic methods; variational integrators
During the last years the development of geometric numerical integrators has been of high interest in numeri-cal integration theory. Geometric integrators are structure-peserving integrators with the goal to capture thedynamical system’s behavior in a most realistic way ([MW01, HLW02, Rei94]). Using structure-preservingmethods for the simulation of mechanical systems, specific properties of the underlying system are handeddown to the numerical solution, for example, the energy of a conservative system shows no numerical drift orthe first integrals induced by symmetries are preserved exactly. One particular class of structure-preservingintegrators is the class of variational integrators, introduced in [MW01, Sur90], and which has been furtherdeveloped and extended to different systems and applications during the last years. Variational integrators([MW01]) are based on a discrete variational formulation of the underlying system, e.g. based on a discreteversion of Hamilton’s principle for conservative mechanical systems. The resulting integrators are symplecticand momentum-preserving and have an excellent long-time energy behavior.By choosing different variational formulations (e.g. Hamilton, Lagrange-d’Alembert, Hamilton-Pontryagin,etc.), variational integrators have been developed for a large class of problems: These involve classical ∗ Corresponding author. Email: [email protected] † Email: [email protected] a r X i v : . [ m a t h . NA ] A p r onservative mechanical systems (for an overview see [LMOW04a, LMOW04b]), forced and controlled sys-tems ([OBJM11, KMOW00]) , constrained systems (holonomic ([LMO08, LOBMO10]) and nonholonomicsystems ([KMS10, CM01])), nonsmooth systems ([FMOW03]), stochastic systems ([BRO08]), multiscale sys-tems ([TOM10, LOB13, SG09]) and Lagrangian PDE systems ([LMOW03, MPS98]). The applicability ofvariational integrators is not restricted to mechanical systems. In [OBTC +
13] variational integrators havebeen developed for the structure-preserving simulation of electric circuits.Of special interest is the construction of higher order symplectic integrators: To ensure moderate computa-tional costs for long-time simulations, typically first or second order integrators are used. However, manyapplications, in particular problems in space mission design, demand more accurate discretization schemes.There are mainly three different ways of constructing symplectic integrators of higher order (cf. [LR04]):(i) By applying composition methods higher order symplectic schemes can be constructed in a system-atic way based on a splitting of the Hamiltonian into explicitly solvable subproblems (for an overview see[MQ02, Yos90]). (ii) For (partitioned) Runge-Kutta methods there is a well-developed order theory whichcan be used to identify higher order symplectic Runge-Kutta methods (order conditions on the coefficientshave first been introduced by [SS88, Las88, Sur89, Sun93]). However, the identification of the coefficientsfor a symplectic Runge-Kutta scheme of desired order is not trivial and quite involved. (iii) In contrast,generating functions (see e.g. [Arn99]) can be constructed that automatically guarantees the symplecticityof the associated numerical method. Since variational integrators rely on the approximation of the action, agenerating function of first kind (cf. [MW01]), we focus on the latter approach for the construction of higherorder methods.Of particular interest are
Galerkin variational integrators which have been already studied in e.g. [MW01,LS11, HL12]. They rely on the approximation of the action based on a choice of a finite-dimensional functionspace and a numerical quadrature formula. In [LS11] Galerkin and shooting-based constructions for discreteLagrangian are presented. Rather than choosing a infinite-dimensional function space, the shooting-basedconstruction depends on a choice of a numerical quadrature formula together with a one-step method. In[HL12] a convergence analysis for Galerkin type variational integrators is established showing that undersuitable assumptions the integrators inherit the order of convergence given by the finite-dimensional approx-imation space used in the construction. More detailed, it is shown that a Galerkin variational integratorbased on an s + 1-dimensional function space (e.g. the space of polynomials of degree s ) and a quadraturerule of order s has convergence oder s .In this contribution we numerically demonstrate that the convergence order of the variational integratorcan even be increased if higher order quadrature rules are used. We focus on a particular class of Galerkinvariational integrators: As finite-dimensional function space we choose the space of polynomials of degree s . Furthermore, as quadrature rules we focus on the Gauss and Lobatto quadrature formula. However, incontrast to [MW01] we do not restrict the number r of quadrature points being equal to the polynomialorder s . For two numerical examples we investigate which combination of space of polynomials, quadraturerules and number of quadrature points provide optimal convergence rates. In particular, the numericalresults indicate that the order of the higher order variational integrator constructed by a polynomial ofdegree s and a quadrature rule of order u is min (2 s, u ). Thus, the integrator order can be increased to2 s for sufficient accurate quadrature rules. While the focus on this work lies on numerical investigations,a formal proof of this superconvergence result is still subject of ongoing research. Based on the numericalresults, we perform a numerical analysis regarding efficiency versus accuracy (see also [Saa12]). Furthermore,we investigate analytically and numerically preservation properties, time reversibility and linear stability ofthe constructed integrators. Whereas preservation properties and time reversibly has also been subject ofprevious works for particular Galerkin variational integrates (see e.g. [LS11]), the stability analysis providesanother new contribution. The generalization r (cid:54) = s is also described in [HL12], however the influence of the relation of the number of quadraturepoints and the polynomial degree is not investigated. utline In Section 2 we recall the basic definitions and concepts of variational mechanics and variationalintegrators. In Section 3 the higher order integrators are constructed following the Galerkin approach intro-duced in [MW01]. Properties of the Galerkin variational integrators are presented in Section 4. In Section 4.1preservation properties, such as symplecticity and preservation of momentum maps are discussed. In Sec-tion 4.2 it is shown under which conditions the constructed integrators are time-reversible. Furthermore, inSection 4.3 a linear stability analysis is performed for specific examples showing in which region the con-structed higher order variational integrators are asymptotically stable (in the sense that the growth of thesolution is asymptotically bounded (cf. [LR04]). A-stability for a particular class of variational integrators isshown. In Section 5 the numerical convergence analysis by means of two numerical examples, the harmonicoscillator and the Kepler problem, is performed using different combinations of the polynomial degree s andthe number r of quadrature points. Furthermore, the relation of computational efficiency and accuracy fortwo different classes of variational integrators is investigated. Finally, we conclude with a summary of theresults and an outlook for future work in Section 6. Consider a mechanical system defined on the n -dimensional configuration manifold Q with correspondingtangent bundle T Q and cotangent bundle T ∗ Q . Let q ( t ) ∈ Q and ˙ q ( t ) ∈ T q ( t ) Q , t ∈ [0 , T ] be the time-dependent configuration and velocity of the system.The action S : C ([0 , T ] , Q ) → R of a mechanical system is defined as the time integral of the Lagrangian,i.e., S ( q ) = (cid:90) T L ( q ( t ) , ˙ q ( t )) dt (2.1)where the C -Lagrangian L : T Q → R consists of kinetic minus potential energy. Hamilton’s principle seekscurves q ∈ C ([0 , T ] , Q ) with fixed initial value q (0) and fixed final value q ( T ) satisfying δ S ( q ) = 0 (2.2)for all variations δq ∈ T q C ([0 , T ] , Q ). This leads to the Euler-Lagrange equations ddt ∂L∂ ˙ q − ∂L∂q = 0 (2.3)which are second-order differential equations describing the dynamics for conservative systems. The concept of variational integrators is based on a discretization of the variational principle (2.2). Considera time grid ∆ t = { t k = kh | k = 0 , . . . , N } , N h = T , where N is a positive integer and h the step size. Wereplace the configuration q ( t ) by a discrete curve q d = { q k } Nk =0 with q k = q d ( t k ) as approximations to q ( t k ).We define a discrete Lagrangian L d : Q × Q → R L d ( q k , q k +1 ) ≈ (cid:90) t k +1 t k L ( q ( t ) , ˙ q ( t )) dt (2.4)that approximates the action on [ t k , t k +1 ] based on two neighboring discrete configurations q k and q k +1 . The discrete action S d : Q N +1 → R is defined as S d ( q d ) = N − (cid:88) k =0 L d ( q k , q k +1 ) . discrete Hamilton principle is formulated by finding stationary points of the discrete action given by δ S d ( q d ) = 0 (2.5)with δq = δq N = 0. This gives the discrete Euler-Lagrange equations (DEL) D L d ( q k , q k +1 ) + D L d ( q k − , q k ) = 0 (2.6)for k = 1 , . . . , N − D i being the derivative w.r.t. the i -th argument. Equation (2.6) providesa discrete iteration scheme for (2.3) that determines q k +1 for given q k − and q k . It is also known as the discrete Lagrangian map F L d : Q × Q → Q × Q , given by F L d ( q k − , q k ) = ( q k , q k +1 ) and ( q k − , q k ) , ( q k , q k +1 )satisfy (2.6). The discrete iteration schemes derived by a discrete variational principle are called variationalintegrators and are well-known to be symplectic and momentum-preserving and exhibit excellent long-timeenergy behavior (cf. Section 4.1).The discrete Legendre transforms F ± L d : Q × Q → T ∗ Q provide discrete expressions for the conjugatemomenta by F − L d : ( q k , q k +1 ) → ( q k , p − k ) = ( q k , − D L ( q k , q k +1 )) and F + L d : ( q k − , q k ) → ( q k , p + k ) = ( q k , D L ( q k − , q k )) . Note that (2.6) can be equivalently written as p − k = p + k and the discrete Hamiltonian map ˜ F L d : T ∗ Q → T ∗ Q defined by ˜ F L d : ( q k , p k ) → ( q k +1 , p k +1 ) = F ± L d ◦ F L d ◦ ( F ± L d ) − ( q k , p k ) . is equivalent to the discrete Lagrangian map. The approximation of the action integral consists of two approximation steps: the approximation of the spaceof trajectories and the approximation of the integral of the Lagrangian by appropriate quadrature rules.We approximate the space of trajectories C ([0 , h ] , Q ) = { q : [0 , h ] → Q | q (0) = q a , q ( h ) = q b } by a finite-dimensional approximation C s ([0 , h ] , Q ) ⊂ C ([0 , h ] , Q ) of the trajectory space given by C s ([0 , h ] , Q ) = { q ∈ C ([0 , h ] , Q ) | q ∈ Π s } , with Π s being the space of polynomials of degree s . Given s + 1 control points 0 = d < d < · · · < d s − 1. To obtain a continuous approximation on [0 , T ], we set q sk = q k +1 for all k = 0 , . . . , N − q ( t ) and ˙ q ( t ) by the piecewisepolynomials q d,k ( t ; q k , h ) and ˙ q d,k ( t ; q k , h ), k = 0 , . . . , N , and approximate (cid:90) ( k +1) hkh L ( q d,k ( t ; q k , h ) , ˙ q d,k ( t ; q k , h )) dt (3.1)on each time interval [ kh, ( k + 1) h ] , k = 0 , . . . , N − 1, by choosing a numerical quadrature rule ( b i , c i ) ri =1 w.r.t. the time interval [0 , 1] with quadrature points c i ∈ [0 , 1] and weights b i , i = 1 , . . . , r . The choice ofquadrature rule should be adapted to the desired order of accuracy of the integrator since the order of thequadrature rule provides an upper bound for the order of the variational integrator (cf. e.g. [LS11]). Byapplying the quadrature rule ( b i , c i ) ri =1 to the integral (3.1) we define the discrete Lagrangian L d,k as L d,k = L d ( q k = ( q k , ..., q sk ) , h ) = h r (cid:88) i =1 b i L ( q d,k ( c i h + kh ; q k , h ) , ˙ q d,k ( c i h + kh ; q k , h ))= h r (cid:88) i =1 b i L ( q d ( c i h ; q k , h ) , ˙ q d ( c i h ; q k , h ))which provides an approximation of the action on the interval [ kh, ( k + 1) h ] as L d,k ≈ (cid:82) ( k +1) hkh L ( q ( t ) , ˙ q ( t )) dt .Note that the discrete Lagrangian depends on s + 1 configurations q k and the step size h . In the following,we write L d ( q k ) for L d ( q k , h ). Finally, we define the discrete action sum over the entire trajectory to be S d ( q , . . . , q N − ) = N − (cid:88) k =0 L d ( q k ) (3.2)which is an approximation of the action sum on [0 , T ] as S d ( q , . . . , q N − ) ≈ S ( q ). For the discrete action S d defined in (3.2) we can apply discrete Hamilton’s principle as described in Sec-tion 2.2. Since we want to determine discrete approximations of curves for which the discrete action is5tationary, the derivatives of the action w.r.t. q νk have to vanish for all k = 0 , . . . , N − ν = 0 , . . . , s .This leads for k = 0 , . . . , N − ν = 1 , . . . , s − ∂ S d ∂q νk ( q , ..., q N − ) = ∂L d ∂q νk ( q k )= h r (cid:88) i =1 b i (cid:18) ∂L∂q ( c i h ; q k ) ∂q d ∂q νk + ∂L∂ ˙ q ( c i h ; q k ) ∂ ˙ q d ∂q νk (cid:19) = h r (cid:88) i =1 b i (cid:18) ∂L∂q ( c i h ; q k ) l ν,s ( c i ) + ∂L∂ ˙ q ( c i h ; q k ) 1 h ˙ l ν,s ( c i ) (cid:19) . Note that we use the short notation ∂L∂q ( c i h ; q k ) = ∂L∂q ( q d ( c i h ; q k , h ) , ˙ q d ( c i h ; q k , h )) and that we have ∂q d ∂q νk = ∂q d ( c i h ; q k ,h ) ∂q νk = l ν,s ( c i ). The analog holds for the other two terms. With q sk − = q k for all k = 1 , . . . , N − ν = 0 and ν = s ∂ S d ∂q sk − ( q , ..., q N − ) = ∂ S d ∂q k ( q , ..., q N − ) = ∂L d ( q k − ) ∂q sk − + ∂L d ( q k ) ∂q k = h r (cid:88) i =1 b i (cid:18) ∂L∂q ( c i h ; q k − ) l s,s ( c i ) + ∂L∂ ˙ q ( c i h ; q k − ) 1 h ˙ l s,s ( c i ) (cid:19) + h r (cid:88) i =1 b i (cid:18) ∂L∂q ( c i h ; q k ) l ,s ( c i ) + ∂L∂ ˙ q ( c i h ; q k ) 1 h ˙ l ,s ( c i ) (cid:19) . With the notation D i L d ( q k , . . . , q sk ) := ∂L d ( q k ) ∂q i − k we obtain the discrete Euler-Lagrange equations D s +1 L d ( q k − , . . . , q sk − ) + D L d ( q k , . . . , q sk ) = 0 , (3.3) D i L d ( q k , . . . , q sk ) = 0 ∀ i = 2 , . . . , s, (3.4)for k = 1 , . . . , N − 1. A sequence { q k } N − k =0 = { ( q k , ..., q sk ) } N − k =0 that satisfies (3.3)-(3.4) and the transitioncondition is a solution of the discrete Euler-Lagrange equations . We denote the left hand side of (3.3) and(3.4) by D DEL L d ( q k − , q k ) such that we have( D DEL L d ( q k − , q k )) = D s +1 L d ( q k − , ..., q sk − ) + D L d ( q k , ..., q sk )( D DEL L d ( q k − , q k )) = D L d ( q k , ..., q sk )...( D DEL L d ( q k − , q k )) s = D s L d ( q k , ..., q sk ) . As in [MW01] we can introduce the standard discrete Lagrangian that depends only on two configurationsas L d ( q k , q sk ) = L d ( q k ) , where q k , . . . , q s − k are implicitly determined by satisfying the internal stage equations (3.4). Alternatively,one can characterize the discrete Lagrangian in the following way (see [MW01]), L d ( q k , q sk ) = ext q νk ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b i L ( q d ( c i h ; q k )) , ˙ q d ( c i h ; q k )) , meaning that ( s − 1) configurations q k , . . . , q s − k are determined by extremizing the discrete Lagrangian. TheLagrangian L d ( q k , q sk ) provides the same iteration scheme as the discrete Lagrangian L d ( q k ).6et q k = ( q k , . . . , q sk ) and let α ( q k , q k +1 ) = q k +1 be the translation operator and π ( q k , q k +1 ) = q k theprojection operator. The discrete Lagrangian evolution operator X L d : Q s +1 → Q s +1 × Q s +1 satisfies X L d ( q k − ) = ( q k − , q k ) with q sk − = q k ,π ◦ X L d ( q k − ) = q k − and D DEL L d ◦ X L d ( q k − ) = 0 . The discrete Lagrangian map F L d : Q s +1 → Q s +1 is defined by F L d ( q k − ) = α ◦ X L d ( q k − ) = q k and generates the sequence of configurations that is denoted as the solution of the Euler-Lagrange equations.The discrete Legendre transforms F ± L d : Q × Q → T ∗ Q are defined as F − L d ( q k , q sk ) = ( q k , p − k ) = ( q k , − D L d ( q k , q sk )) , F + L d ( q k − , q sk − ) = ( q k , p k ) = ( q k , D s +1 L d ( q k − , q sk − )) . From the discrete Euler-Lagrange equations it follows that along the solution of the discrete Euler-Lagrangeequations we have p k := p − k = p k . Note that the discrete Lagrangian flow is well-defined, if the discrete Lagrangian is regular, i.e., if the discreteLegendre transforms are local isomorphisms what is assumed in the following. As shown in Fig. 3.1, the discrete Hamiltonian map ˜ F L d : T ∗ Q → T ∗ Q is given by˜ F L d = F ± L d ◦ F L d ◦ ( F ± L d ) − . ( q , . . . , q s ) ( q , . . . , q s )( q , p ) ( q , p ) ( q , p ) F L d F − L d F + L d F − L d F + L d ˜ F L d ˜ F L d Figure 3.1: Correspondence between the discrete Lagrangian and the discrete Hamiltonian map.With the proposed method different variational integrators can be constructed. We use the following no-tation: ( P sN rQu ) is an integrator constructed as described above with a polyomial of degree s with s + 1control points ( d i ) si =0 and a quadrature formula of order u with r quadrature points. Note that u depends on r . If explicitly given, we denote by three letters the quadrature rule in use, i.e. Lob for Lobatto quadrature( u = 2 r − 2) and Gau for Gauss quadrature ( u = 2 r ). Example 3.1. (i) The integrator ( P N Q Gau ) is based on a polynomial with control points d = 0 , d = 1 and the quadrature approximation (cid:90) h L ( q ( t ) , ˙ q ( t )) dt ≈ L d (( q , q ) , h ) = hL ( q d ( h/ , ˙ q d ( h/ hL (cid:18) q + q , q − q h (cid:19) and thus results in the midpoint rule discrete Lagrangian. ii) If the trapezoidal rule is applied as quadrature formula, we obtain the variational integrator ( P N Q Lob ) with discrete Lagrangian L d (( q , q ) , h ) = h L ( q d (0) , ˙ q d (0)) + h L ( q d ( h ) , ˙ q d ( h ))= h L (cid:18) q , q − q h (cid:19) + h L (cid:18) q , q − q h (cid:19) which is the discrete Lagrangian of the St¨ormer-Verlet method with factor for a Lagrangian of theform L ( q, ˙ q ) = ˙ q T M ˙ q − V ( q ) .(iii) For a second order polynomial ( s = 2 ) with control points d = 0 , d = and d = 1 and byapplying Simpson’s rule, which is a Lobatto quadrature of order four, we obtain the variational integrator ( P N Q Lob ) with discrete Lagrangian given by L d (( q , q , q ) , h ) = h L ( q d (0) , ˙ q d (0)) + 2 h L ( q d ( h/ , ˙ q d ( h/ h L ( q d ( h ) , ˙ q d ( h ))= h L (cid:18) q , − q + 4 q − q h (cid:19) + 2 h L (cid:18) q , q − q h (cid:19) + h L (cid:18) q , q − q + 3 q h (cid:19) . Remark 3.2 (Implementation) . For given configurations q , . . . , q s − , q we compute the unknown configu-ration q by performing one step of the discrete Lagrangian evolution ( q , ..., q s − , q ) F Ld −→ ( q , ..., q s − , q ) . To this end, the system of s nonlinear equations (3.3)-(3.4) is solved for the s unknowns ( q , ..., q s − , q ) . Forthe numerical solution, a Newton method can be used. For given initial configuration q (0) and momentum p (0) , in the first step, (3.3) is replaced by the discrete Legendre transform p (0) = − D L d ( q , . . . , q s ) and the system of equations is solved for ( q , . . . , q ) . For the quadrature rules considered here, (3.3)-(3.4)typically provide an implicit scheme for nonlinear systems that has to be solved by an iterative solver all atonce. Remark 3.3 (Galerkin methods) . For the Galerkin variational integrators as introduced in [MW01], Sec-tion 2.2.6, the number of quadrature points of the quadrature formula ( b i , c i ) ri =1 is fixed to r = s , where s is the degree of the polynomial q d . In our notation that means that only the methods ( P sN sQu ) are inves-tigated, which are shown to be equivalent to partitioned symplectic Runge-Kutta methods. In particular, itis pointed out that the integrator ( P sN sQ sGau ) , which uses the Gauss quadrature formula, corresponds tothe collocation Gauss-Legendre rule , whereas ( P sN sQ s − Lob ) yields the standard Lobatto IIIA-IIIB par-titioned Runge-Kutta method. For both methods the order is determined by the order of quadrature rule, i.e., ( P sN sQ sGau ) is of order s and ( P sN sQ s − Lob ) is of order s − (cf. [HLW02]). Although the Gaussquadrature formula leads to higher order schemes, in particular for stiff systems the choice of a quadraturerule involving c s = 1 leads to better numerical performance (cf. [HLW02, MW01]). If, in addition, onewishes to use a symmetric quadrature rule, i.e., c = 0 , the Legendre-Lobatto quadrature rule provides thehighest possible order. In this section properties of the Galerkin variational integrators, such as symplecticity, momentum-preservation,time reversibility and linear stability are studied. 8 .1 Preservation properties of variational integrators As already mentioned in Section 2 variational integrators are structure-preserving, in particular they aresymplectic and momentum-preserving. In this section we briefly repeat the notion of symplecticity and firstintegrals and their discrete counterparts. In the case of conservative systems (as considered here), the flow on T ∗ Q of the Euler-Lagrange equationspreserves the canonical symplectic form Ω = dq i ∧ dp i = d θ of the Hamiltonian system, where θ = p i dq i isthe canonical one-form. It is well known (cf. e.g. [MW01]) that variational integrators are symplectic, thatis the same property holds for the discrete flow of the discrete Euler-Lagrange equations. As a consequence,the canonical discrete symplectic form Ω = dq i ∧ dp i is exactly preserved for the discrete solution which isin particular true for the variational integrators presented in this work.By using techniques from backward error analysis it is shown (cf. e.g. [HLW02]) that symplectic integratorshave excellent energy properties, meaning that for long-time integrations there is no artificial energy growthor decay due to numerical errors. This property is demonstrated numerically in Section 5 and illustratesthe advantage of variational integrators over e.g. nonsymplectic Runge-Kutta integrators in particular forlong-time simulations. The Noether theorem provides first integrals of the Euler-Lagrange equations which are also called momentummaps . The discrete Noether theorem states that these invariants are also preserved for the discrete solution.The following two theorems are taken from [HLW02]. Theorem 4.1 (Noether theorem) . Let L ( q, ˙ q ) be a regular Lagrangian. Suppose G = { g v : v ∈ R } is aone-parameter group of transformations ( g v ◦ g w = g v + w ) which leaves the Lagrangian invariant such that L ( g v ( q ) , g (cid:48) v ( q ) ˙ q ) = L ( q, ˙ q ) ∀ v ∈ R ∀ ( q, ˙ q ) ∈ T Q. Let a ( q ) = ddv g v ( q ) | v =0 be defined as the vector field with flow g v ( q ) . Then I ( q, p ) = p T a ( q ) (4.1) is a first integral of the Euler-Lagrange equations. The discrete analog of the Noether theorem is stated as follows. Theorem 4.2 (Discrete Noether theorem) . Suppose the one-parameter group of transformations G = { g v : v ∈ R } leaves the discrete Lagrangian L d ( q , q ) invariant, that means L d ( g v ( q ) , g v ( q )) = L d ( q , q ) ∀ v ∈ R ∀ ( q , q ) ∈ Q × Q. Then (4.1) is an invariant of the discrete Hamiltonian map ˜ F L d , i.e., I ◦ ˜ F L d ( q k , p k ) = I ( q k , p k ) . Proofs of Theorems 4.2 and 4.2 can be found in [HLW02].Note that the invariant of the discrete Hamiltonian map only equals the first integral of the Euler-Lagrangeequations if the discrete Lagrangian L d inherits the invariance of the Lagrangian L . In the following we showunder which condition this invariance is inherited. 9 efinition 4.3 (Equivariance) . Let { g v : v ∈ R } be a one-parameter group of transformation. The interpo-lation polynomial q d of ( P sN sQu ) is equivariant w.r.t. { g v : v ∈ R } if g v ( q d ( t ; ( q , . . . , q s ))) = q d ( t ; ( g v ( q ) , . . . , g v ( q s ))) , (4.2) g (cid:48) v ( q d ( t ; ( q , . . . , q s ))) ˙ q ( t ; ( q , . . . , q s )) = ˙ q d ( t ; ( g v ( q ) , . . . , g v ( q s ))) . (4.3)Note that (4.3) follows from (4.2) by applying the chain rule. Theorem 4.4 (Invariance of discrete Lagrangian) . Suppose that the interpolation polynomial q d of ( P sN rQu ) is equivariant and that the regular Lagrangian L is invariant w.r.t. the one-parameter group G = { g v : v ∈ R } .Then the discrete Lagrangian L d of ( P sN rQu ) is invariant w.r.t. G . Proof. Let q = ( q , . . . , q s ) with q s = q and g v · q = ( g v ( q ) , ..., g v ( q s − ) , g v ( q )). Let ( b i , c i ) ri =1 be thequadrature formula that corresponds to ( P sN rQu ). With the invariance of L and the equivariance of q d wehave that L d ( g v ( q ) , g v ( q )) = ext g v ( q ν ) ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b i L ( q d ( c i h ; g v · q ) , ˙ q d ( c i h ; g v · q ))= ext q ν ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b i L ( g v ( q d ( c i h ; q )) , g (cid:48) v ( q d ( c i h ; q )) ˙ q d ( c i h ; q ))= ext q ν ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b i L ( q d ( c i h ; q )) , ˙ q d ( c i h ; q ))= L d ( q , q ) . (cid:4) A general form of the statement of Theorem 4.4 for Galerkin Lie group variational integrators can also befound in [LS11]. Remark 4.5. From Theorem 4.4 it follows that for linear group transformations the discrete Lagrangian L d of ( P sN rQu ) inherits the invariance of the Lagrangian L (cf. [MW01]). The properties described in Section 4.1 are valid for all variational integrators. In the following we presentspecial properties of the variational integrators constructed in this work. A further geometric property of Hamiltonian systems is the time reversibility. It seems likely to use numericalmethods that produce a reversible numerical flow when they are applied to a reversible Hamiltonian system.Furthermore, the numerical solution has a long-time behavior similar to the exact solution (see [HLW02]),and thus time reversibility is a desirable property also for the variational integrators presented in this work.In the following, we repeat the definitions of time reversibility and the adjoint discrete Lagrangian (following[HLW02] and [MW01]) und show, under which conditions the variational integrators ( P sN rQu ) are time-reversible. Definition 4.6 (Time reversibility and adjoint ([HLW02])) . A numerical one-step method Φ d is called symmetric or time-reversible if it satisfies Φ hd ◦ Φ − hd = id or equivalently Φ hd = (Φ − hd ) − . he adjoint method denoted by Φ ∗ h is defined by (Φ hd ) ∗ = (Φ − hd ) − . A method is called self-adjoint if we have Φ ∗ d = Φ d . Thus, symmetric, time-reversible, and self-adjoint are equivalent notions. In the following, we show that theintegrators ( P sN rQ rGau ) and ( P sN rQ r − Lob ) are self-adjoint and thus time-reversible methods. Definition 4.7 (Adjoint of the discrete Lagrangian ([MW01])) . The adjoint discrete Lagrangian L ∗ d of thediscrete Lagrangian is given by L ∗ d ( q , q , h ) = − L d ( q , q , − h ) . The discrete Lagrangian is self-adjoint if L d ( q , q , h ) = L ∗ d ( q , q , h ) . The following well-known theorem (see e.g. [MW01]) connects the adjoint of the discrete Lagrangian withthe adjoint of the discrete Hamiltonian flow. Theorem 4.8. If the discrete Lagrangian L d has a discrete Hamiltonian map ˜ F L d , then the discrete Hamil-tonian map of the adjoint discrete Lagrangian L ∗ d equals the adjoint map, i.e. ˜ F L ∗ d = ˜ F ∗ L d . If the discreteLagrangian is self-adjoint, then the method is self-adjoint. Conversely, if the method is self-adjoint, then thediscrete Lagrangian is equivalent to a self-adjoint discrete Lagrangian. A proof of this theorem can be found in [MW01], Theorem 2.4.1. To show the time reversibility, we firstshow that the discrete Lagrangian is self-adjoint. Theorem 4.9. Let L d be the discrete Lagrangian of ( P sN rQu ) with symmetric quadrature formula ( b i , c i ) ri =1 and interpolation polynomial q d with symmetrically distributed control points ( d i ) si =0 , i.e., b i = b r +1 − i , c i + c r +1 − i = 1 , i = 1 , . . . , r , and d i = 1 − d s − i , i = 0 , . . . , s . Then L d is self-adjoint. Proof. We have q d ( − t ; q , − h ) = s (cid:88) ν =0 q ν l ν,s ( t/h ) = q d ( t ; q , h ) , (4.4)˙ q d ( − t ; q , − h ) = 1 − h s (cid:88) ν =0 q ν ˙ l ν,s ( t/h ) = − ˙ q d ( t ; q , h ) . (4.5)Due to the symmetry of the control points d i it follows for the Lagrange-polynomials l ν,s with τ ∈ [0 , l ν,s (1 − τ ) = (cid:89) ≤ i ≤ s,i (cid:54) = ν (1 − τ ) − d i d ν − d i = (cid:89) ≤ i ≤ s,i (cid:54) = ν ( − τ ) + (1 − d i )( − d ν ) + (1 − d i )= (cid:89) ≤ i ≤ s,i (cid:54) = ν ( − τ ) + d s − i − d s − ν + d s − i = (cid:89) ≤ i ≤ s,i (cid:54) = ν τ − d s − i d s − ν − d s − i = l s − ν,s ( τ ) . Two discrete Lagrangian are equivalent if their discrete Hamiltonian map are the same. q = ( q , ..., q s ), q s = q , and ˜ q = ( q s , ..., q ) we have for t ∈ [0 , h ] q d ( t, ˜ q , h ) = s (cid:88) ν =0 q s − ν l ν,s (cid:18) th (cid:19) = s (cid:88) k =0 q k l s − k,s (cid:18) th (cid:19) = s (cid:88) k =0 q k l k,s (cid:18) − th (cid:19) = s (cid:88) k =0 q k l k,s (cid:18) h − th (cid:19) = q d ( h − t ; q , h )and by taking the time derivative we have˙ q d ( h − t ; q , h ) = − ˙ q d ( t ; ˜ q , h ) . Thus, together with (4.4)-(4.5) we obtain q d ( − t ; ˜ q , − h ) = q d ( t ; ˜ q , h ) = q d ( h − t ; q , h ) , ˙ q d ( − t ; ˜ q , − h ) = − ˙ q d ( t ; ˜ q , h ) = ˙ q d ( h − t ; q , h ) . By substituting these expressions in the adjoint discrete Lagrangian we have with the symmetry of thequadrature formula ( b i , c i ) ri =0 − L d ( q , q , − h ) = ext q ν ∈ Qν ∈{ ,...,s − } − ( − h ) r (cid:88) i =1 b i L ( q d ( − c i h ; ˜ q , − h )) , ˙ q d ( − c i h ; ˜ q , − h ))= ext q ν ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b i L ( q d ( h − c i h ; q , h )) , ˙ q d ( h − c i h ; q , h ))= ext q ν ∈ Qν ∈{ ,...,s − } h r (cid:88) i =1 b r +1 − i L ( q d ( c r +1 − i h ; q , h )) , ˙ q d ( c r +1 − i h ; q , h ))= L d ( q , q ) . (cid:4) Theorem 4.10. The discrete Hamiltonian maps of ( P sN rQ rGau ) and ( P sN rQ r − Lob ) are self-adjointand thus time-reversible. Proof. Since the Gauss and the Lobatto quadrature are symmetric and since the control points of q d arechosen symmetrically, it follows with Theorem 4.9 that the discrete Lagrangian of ( P sN rQ rGau ) and( P sN rQ r − Lob ) are self-adjoint. The statement follows with Theorem 4.8. (cid:4) In the following, we investigate the stability properties of the constructed variational integrators. We restrictourselves to a linear stability analysis. Following [LR04], we consider the Lagrangian of a harmonic oscillator L ( q, ˙ q ) = 12 ˙ q − ω q (4.6)with ω, q ∈ R . The Hamiltonian equations read˙ p = − ω q, ˙ q = p (cid:18) p ( t ) ωq ( t ) (cid:19) = (cid:18) cos ( ωt ) − sin ( ωt )sin ( ωt ) cos ( ωt ) (cid:19) (cid:18) p (0) ωq (0) (cid:19) = M ω (cid:18) p (0) ωq (0) (cid:19) with det( M ω ) = 1. The eigenvalues of M ω are λ , = e ± iwt and thus, we have that | λ , | = 1 . By applying a variational integrator ( P sN rQu ) to the Lagrangian (4.6), the discrete Euler-Lagrange equa-tions form a linear system of equations such that the discrete Hamiltonian map ˜ F L d of ( P sN rQu ) can bewritten as linear map (cid:18) p wq (cid:19) = M h,w (cid:18) p ωq (cid:19) with matrix M h,ω .Following [LR04], we call a numerical solution asymptotically stable if the growth of the solution is asymp-totically bounded. A sufficient condition for asymptotic stability is that the eigenvalues of M h,ω are in theunit desk of the complex plane and are simple if they lie on the unit circle. We investigate this property forselected variational integrators.1. For the midpoint rule ( P N Q Gau ) we have M h,ω = (cid:32) − h ω − h ω +4 − hωh ω +44 hωh ω +4 − h ω − h ω +4 (cid:33) . Since M h,ω M Th,ω = Id, M h,ω is orthogonal with | λ ( M h,ω ) | = 1. Thus, the midpoint rule is asymptoti-cally stable for all h, ω ∈ R .2. The St¨ormer-Verlet method ( P N Q Lob ) has the iteration matrix M h,ω = (cid:32) − h ω h ω − hωhω − h ω (cid:33) . with eigenvalues λ , = 1 − h ω ± hω √ h ω − , thus, the method is stable for ( hω ) < P N Q Gau ) method we have M h,ω = h ω − h ω +144 h ω +12 h ω +144 − hω ( h ω − ) h ω +12 h ω +14412 hω ( h ω − ) h ω +12 h ω +144 h ω − h ω +144 h ω +12 h ω +144 , and the method is stable for all q, ω ∈ R since M h,ω is orthogonal.4. The ( P N Q Lob ) scheme results in the iteration matrix M h,ω = (cid:32) h ω − h ω +482 h ω +48 24 hω − h ω h ω +24 − hω ( h ω − h ω +288 ) h ω +288 h ω − h ω +482 h ω +48 (cid:33) with eigenvalues λ , = h ω − h ω + 48 ± hω √ h ω − h ω + 576 h ω − h ω + 48 . The stability region is shown in Fig. 4.1. The integrator is stable for | hω | < √ λ , | hw Figure 4.1: Modulus of the eigenvalues of M h,ω for the ( P N Q Lob ) integrator. The method is stable for | hω | < √ P N Q Gau ) reads M h,ω = − h ω − h ω +6480 h ω − h ω +24 h ω +720 h ω +14400 24 hω ( h ω − h ω +600 ) h ω +24 h ω +720 h ω +14400 − hω ( h ω − h ω +600 ) h ω +24 h ω +720 h ω +14400 − h ω − h ω +6480 h ω − h ω +24 h ω +720 h ω +14400 . The scheme is again stable for all q, ω ∈ R due to the orthogonality of M h,ω .6. The integrator ( P N Q Lob ) gives M h,ω = − h ω − h ω +840 h ω − h ω +60 h ω +1800 6 hω ( h ω − h ω +300 ) h ω +60 h ω +1800 hω ( h ω − h ω +5760 h ω − ) h ω +60 h ω +1800) − h ω − h ω +1680 h ω − h ω +120 h w +3600 . The eigenvalues are given as λ , = 3600 − x + 92 x − x ± x √ x − x + 11820 x − x + 3456000 x − x + 120 x + 3600with x := hω . In Figs 4.2 and 4.3 (zoom of Fig. 4.2) the stability region is shown. | λ , | hw Figure 4.2: Modulus of the eigenvalues of M h,ω for the ( P N Q Lob ) integrator.An interesting observation is that for the integrators ( P sN sQ sGau ), s = { , , } the iteration matrices M h,ω are orthogonal independent of the step size h and thus asymptotically stable. Indeed, we can showthat this unrestricted stability property holds for general integrators of this type.14.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.180.980.991.001.011.02 hw | λ , | Figure 4.3: Modulus of the eigenvalues of M h,ω for the ( P N Q Lob ) integrator (zoom). Lemma 4.11 (A-stability of ( P sN sQ sGau )) . The variational integrator ( P sN sQ sGau ) is A-stable. Proof. In [MW01] it is shown that the integrator ( P sN sQ sGau ) is equivalent to the collocation Gauss-Legendre rule (see also Remark 3.3) which is A-stable (see [HW10], Theorem 5.2). (cid:4) In this section we numerically analyze the convergence order of the constructed variational integrators( P sN rQu ) for s, r, u ∈ N and the theoretical results on the preservation properties are evaluated. To thisend, we consider two examples, the harmonic oscillator and the Kepler problem and we want to numericallydetermine the maximal order of the variational integrator ( P sN rQu ).In Section 3.2 we assume that the discrete Lagrangian is regular to obtain a well-defined discrete Lagrangianflow. In [HL12] it is shown that the well-posedness depends on the order of the quadrature rule thatis used to approximate the action integral. In particular, it is shown that for a Lagrangian of the form L = ˙ q T M ˙ q − V ( q ) with M symmetric positive-definite and ∇ V Lipschitz continuous a unique solution ofthe internal stage equations (3.4) exists if the used quadrature rule is of order at least 2 s − s thedegree of the interpolation polynomial), i.e. u ≥ s − 1. For the Gauss and Lobatto quadrature this meansthat we have to choose r ≥ s and r ≥ s + 1, respectively, since this yields quadrature rules of order 2 s . Notethat the variational integrator ( P sN sQ s − Lob ) with r = s which yields the standard Lobatto IIIA-IIIBpartitioned Runge-Kutta method does not satisfy the condition u ≥ s − 1, however, it is only a sufficientnot a necessary condition for uniqueness of solutions. Thus, to ensure a well-defined discrete Lagrangianflow, we restrict to variational integrators with quadrature rule of Gauss or Lobatto type and for whichthe polynomial degree is smaller or equal to r ( s ≤ r ). The implementation is performed as described inRemark 3.2. Consider the two-dimensional harmonic oscillator with mass being equal to one and the Lagrangian L ( q, ˙ q ) = ˙ q T ˙ q − q T q with q, ˙ q ∈ R . The Euler-Lagrange equations are¨ q ( t ) = − q ( t ) . P sN rQ rGau ) r \ s s = 1 s = 2 s = 3 s = 4 s = 5 r = 2 2 4 r = 3 2 4 6 r = 4 2 4 6 8 r = 5 2 4 6 8 10By the Legendre transform we obtain the Hamiltonian system˙ q ( t ) = p ( t ) , ˙ p ( t ) = − q ( t ) . (5.1)The total energy of the system is given by the Hamiltonian H ( q, p ) = p T p + q T q with q, p ∈ R . Let ( q e , p e )( t ) denote the exact solution consisting of configuration and momentum of the Hamiltoniansystem (5.1). For the error calculations we use the global error determined bymax k ∈{ ,...,N } i ∈{ , } | q k,i − q ei ( hk ) | , max k ∈{ ,...,N } i ∈{ , } | p k,i − p ei ( hk ) | (5.2)with step size h and the index i denotes the components of the configuration and momenta, respectively.( q k , p k ) Nk =0 is the discrete solution computed by a variational integrator with q N = q sN − and p N = p sN − . InFigs 5.1 and 5.3 the global error for the different variational integrators ( P sN rQu ) is shown in dependenceof the step size h . The numerically determined order is given in Figs 5.2 and 5.4 and summarized in theTables 1 and 2 . Note that the values for r = s in both tables correspond to the convergence orders of thecollocation Gauss-Legendre rule and the Lobatto IIIA-IIIB partitioned Runge-Kutta method, respectively.The first observation by considering the values in the Tables 1 and 2 is that the convergence order for allcombinations of polynomial degree s and number r of quadrature points can be determined as min (2 s, u )with u = 2 r and u = 2 r − s the order of convergence can not be improved for any quadrature rule based on morequadrature points than s . Thus, a reasonable combination of polynomial degree and number of quadraturepoints is r ≥ s . However, for the Galerkin variational integrators based on the Lobatto quadrature formula(Table 2), the order of convergence can be improved if the number of quadrature points is increased by atleast one, i.e. the best combinations satisfy r ≥ s +1. This demonstrates that the order of ( P sN sQ s − Lob )increases to 2 s if a Lobatto quadrature with s + 1 quadrature points is used, or, the order does not decreaseif a polynomial of degree s − s is used.Based on this numerical results we can conclude that for a variational integrator based on an approximationspace of degree s polynomials the convergence order 2 s is possible, i.e. the integrator is superconvergent. Thismaximal order is reduced if the quadrature rules used for the approximation of the action is not accurateenough. Thus, to guarantee the maximal convergence order we have to choose u ≥ s . The numerical results also indicate which combination of quadrature rule and polynomial degree leads tolowest computational effort for a given order. In Fig. 5.13 the run-time compared with the global error is Note that in Table 2 the numerical results for r = 2 are not shown in Fig. 5.3 and 5.4. −2 −1 −15 −10 −5 P1N2Q4Gauss quadrature of order 4h E rr o r P2N2Q4 10 −2 −1 −15 −10 −5 P1N3Q6Gauss quadrature of order 6h E rr o r P2N3Q6 P3N3Q6 10 −2 −1 −15 −10 −5 P1N4Q8Gauss quadrature of order 8h E rr o r P2N4Q8 P3N4Q8 P4N4Q8 10 −2 −1 −15 −10 −5 P1N5Q10Gauss quadrature of order 10h E rr o r P2N5Q10 P3N5Q10 P4N5Q10 P5N5Q10 Figure 5.1: Harmonic oscillator: Log-log plot of the error for the configurations q and the momenta p (superimposed) by step size h ; with the use of variational integrators ( P sN rQ rGau ); divided into foursubplots, separated by applied Gauss quadrature. P1 P20246810 O r de r Gauss quadrature of order 4 P1 P2 P30246810 O r de r Gauss quadrature of order 6 P1 P2 P3 P40246810 O r de r Gauss quadrature of order 8 P1 P2 P3 P4 P50246810 O r de r Gauss quadrature of order 10 Figure 5.2: Harmonic oscillator: Order of the variational integrators ( P sN rQ rGau ) with respect to theconfigurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Gaussquadrature which is 2 r . The x -axis indicates the degree of the used polynomial.given for the different integrators ( P sN rQ r − Lob ) with s = { , ..., } and r = { , ..., } and for differentstep sizes h ∈ { , . , . , . , . , . , . } . The step size is reduced until the error is below 10 − .The integrators ( P N Q Lob ), ( P N Q Lob ), ( P N Q Lob ), and ( P N Q Lob ), which are of order four,demonstrate that for an increasing number of nodes also the run-time increases. The same behavior is alsoobservable for the other integrators.Further, a higher polynomial degree leads to an increasing number of Euler-Lagrange equations that haveto be solved for and thus to higher computational effort. In Fig. 5.13 it is shown that this leads to anincreasing run-time (compare for example ( P N Q Lob ) and ( P N Q Lob ) which are both of order four; or( P N Q Lob ) and ( P N Q Lob ), which are of order ten). Since the order of the integrator is min (2 s, u ),a reasonable choice for the polynomial degree is half of the order of the quadrature formula, i.e., u = 2 s .This guarantees a minimal number of discrete Euler-Lagrange equations without an order reduction of thevariational integrator. For the Gauss and Lobatto quadrature rule the optimal combinations of polynomialdegree s and number r of quadrature points is r = s and r = s + 1, respectively. That is, the integrators( P sN sQ sGau ) and ( P sN s + 1 Q sLob ) are the most efficient integrators with view of run-time per order.For a clearer illustration we omit in Fig. 5.14 the less efficient integrators which are displayed in Fig. 5.13 andinclude the most efficient integrators constructed with the Gauss and Lobatto quadrature. It can be readoff which of these integrators provide the desired accuracy with shortest run-time. Larger s means longer17 −2 −1 −15 −10 −5 P1N3Q4Lobatto quadrature of order 4h E rr o r P2N3Q4 P3N3Q4 10 −2 −1 −15 −10 −5 P1N4Q6Lobatto quadrature of order 6h E rr o r P2N4Q6 P3N4Q6 P4N4Q6 10 −2 −1 −15 −10 −5 P1N5Q8Lobatto quadrature of order 8h E rr o r P2N5Q8 P3N5Q8 P4N5Q8 P5N5Q8 10 −2 −1 −15 −10 −5 P1N6Q10Lobatto quadrature of order 10h E rr o r P2N6Q10 P3N6Q10 P4N6Q10 P5N6Q10 P6N6Q10 Figure 5.3: Harmonic oscillator: Log-log plot of the error for the configurations q and the momenta p (superimposed) by step size h ; with the use of variational integrators ( P sN rQ r − Lob ); divided into foursubplots, separated by applied Lobatto quadrature. Note that ( P s − N sQ s − Lob ) and ( P sN sQ s − Lob )are of the same order. P1 P2 P30246810 O r de r Lobatto quadrature of order 4 P1 P2 P3 P40246810 O r de r Lobatto quadrature of order 6 P1 P2 P3 P4 P50246810 O r de r Lobatto quadrature of order 8 P1 P2 P3 P4 P5 P60246810 O r de r Lobatto quadrature of order 10 Figure 5.4: Harmonic oscillator: Order of the variational integrators ( P sN rQ r − Lob ) with respect to theconfigurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Lobattoquadrature which is 2 r − 2. The x -axis indicates the degree of the used polynomial.processing time but also higher order. Notice that ( P N Q Gau ) and ( P N Q Lob ) are the midpoint ruleand the St¨ormer-Verlet method, respectively. Positioning the system of the two-dimensional harmonic oscillator in the ( x, y )-plane of the three-dimensionalspace, the z -component of the angular momentum I ( p, q ) = − p q + p q is a conserved quantity. This follows from Noether’s theorem because the Lagrangian of the two-dimensionalharmonic oscillator L = ˙ q T ˙ q − q T q is invariant under the group of rotations SO (2) = { B ∈ R × | B T B =Id , det( B ) = 1 } . With the linearity of R v ∈ SO (2) and Remark 4.5 it follows that all variational integrators( P sN rQu ) conserve the z -component of the angular momentum. In Fig. 5.5 it is shown that the z -componentof the angular momentum, if the variational integrators ( P N Q Lob ), ( P N Q Lob ) and ( P N Q Lob )are used, is preserved up to an error less than 10 − . In Fig. 5.5 (right) the behavior of a Runge-Kutta Note that the accuracy is limited to machine precision and the accuracy of the applied Newton method to solve the discreteEuler-Lagrange equations. P sN rQ r − Lob ) r \ s s = 1 s = 2 s = 3 s = 4 s = 5 s = 6 r = 2 2 2 r = 3 2 4 4 r = 4 2 4 6 6 r = 5 2 4 6 8 8 r = 6 2 4 6 8 10 10method of order four is included which is not symplectic nor momentum-preserving. Thus, Noether’s theoremdoes not apply and the angular momentum is not preserved. In Section 4.1 we introduced the good long-time −15 t A ngu l a r M o m en t u m E rr o r −11 P2N3Q4; h=0.5t A ngu l a r M o m en t u m E rr o r P3N4Q6; h=0.5 P4N5Q8; h=0.5RK4; h=0.015625 Figure 5.5: Harmonic oscillator: Left: Error of the z -component of the angular momentum for the variationalintegrators ( P N Q Lob ) (dots), ( P N Q Lob ) (crosses) and ( P N Q Lob ) (squares) with step size h =0 . 5. Right: Same plot as on the left with non-variational integrator included (Runge-Kutta method of orderfour with step size h = 2 − (red solid)) .energy behavior of symplectic integrators, in particular of variational integrators. In Fig. 5.6 the error ofthe total energy of the harmonic oscillator simulated by different integrators is shown. While the use of thevariational integrators ( P N Q Lob ) and ( P N Q Lob ) leads to an oscillating but stable energy behavior,the use of a nonsymplectic Runge-Kutta method of order four clearly shows an energy drift. The two-body problem, also known as Kepler problem, is to determine the motion of two point particles withmasses m , m ∈ R . By assuming m = 1 with gravitational constant γ and k = γm m ∈ R we constructthe Lagrangian of the Kepler problem L ( q, ˙ q ) = 12 ˙ q T ˙ q + k (cid:112) q + q (5.3)and the Hamiltonian H ( q, p ) = 12 p T p − k (cid:112) q + q with q = ( q , q ) T , ˙ q, p ∈ R . The Hamiltonian equations provide the system of differential equationswhich we want to solve with respect to the initial condition ( q , p ). For the simulations we set k =19 50 100 150 200−4−3−2−101234 x 10 −9 t E ne r g y − E rr o r P3N4Lob P4N5LobRK4 Figure 5.6: Harmonic oscillator: Error of the energy for the integrators ( P N Q Lob ) with step size h = 0 . P N Q Lob ) with step size h = 0 . h = 2 − (solid line).1 . · , q = (5 , T and p = (0 , T since this results in a motion where the mass m describes an elliptic orbit around mass m with period T = 5 . 0. Therefore, after h steps of step size h thesimulated mass should return the fifth time to the initial point ( q , p ). As error we compute the maximal difference between the given initial value and the value after N = h stepsof integration given as max i ∈{ , } | q N,i − q ,i | and max i ∈{ , } | p N,i − p ,i | for configuration and momentum, respectively. The error is shown in Fig. 5.7 and Fig. 5.9 for the differ-ent variational integrators ( P sN rQu ) and for step sizes h ∈ { , . , . , . , . , . , . } . Thenumerically determined order is given in Fig. 5.8 and Fig. 5.10 and coincides nicely with the orders for theharmonic oscillator as given in Tables 1 and 2. It can be observed in the plots that the error of the integratorsdecreases not far below 10 − . Since the iteration is implicit and has to be solved by a Newton method, theaccuracy is limited by the machine precision and the used solver (we used fsolve implemented in Matlab).20 −2 −1 −12 −10 −8 −6 −4 −2 P1N2Q4Gauss quadrature of order 4h E rr o r P2N2Q4 10 −2 −1 −10 −5 P1N3Q6Gauss quadrature of order 6h E rr o r P2N3Q6 P3N3Q6 10 −2 −1 −10 −5 P1Gauss quadrature of order 8h E rr o r P2 P3 P4 10 −2 −1 −10 −5 P1Gauss quadrature of order 10h E rr o r P2 P3 P4 P5 Figure 5.7: Kepler problem: Log-log plot of the error for the configurations q (crosses) and momenta p (dots)with step size h for the integrators ( P sN rQ rGau ); divided into four subplots, separated by applied Gaussquadrature. In the third subplot P s denotes ( P sN Q Gau ) and in the fourth P s denotes ( P sN Q Gau ). P1 P20246810 O r de r Gauss quadrature of order 4 P1 P2 P30246810 O r de r Gauss quadrature of order 6 P1 P2 P3 P40246810 O r de r Gauss quadrature of order 8 P1 P2 P3 P4 P5 P5*0246810 O r de r Gauss quadrature of order 10 Figure 5.8: Kepler problem: Order of the variational integrators ( P sN rQ rGau ) with respect to the config-urations q (left bar) and the momenta p (right bar). The title gives the order of the applied Gauss quadraturewhich is 2 r . The x -axis indicates the degree of the used polynomial. P P ∗ gives the numerically determined order with use of thefirst five values. Since the Lagrangian (5.3) of the Kepler problem is invariant under the group of rotations SO (2) = { B ∈ R × | B T B = Id , det( B ) = 1 } , the angular momentum I ( p, q ) = − p q + p q is a conserved quantity in the system. The angular momentum is (up to numerical accuracy) also preservedin the discrete solution using the variational integrators ( P sN rQu ) (cf. Remark 4.5) as shown in Fig. 5.11for the integrators ( P N Q Gau ), ( P N Q Gau ) and ( P N Q Gau ). However, using a Runge-Kuttaintegrator of order four, angular momentum is not preserved anymore as shown in Fig. 5.11 (right). Theerror in the energy is given in Fig. 5.12. As for the harmonic oscillator, due to the symplecticity of thevariational integrators the error of the integrators ( P N Gau ) and ( P N Gau ) is oscillating but boundedwhereas the Runge-Kutta solution exhibits an energy drift.21 −2 −1 −12 −10 −8 −6 −4 −2 P1N3Q4Lobatto quadrature of order 4h E rr o r P2N3Q4 P3N3Q4 10 −2 −1 −10 −5 P1N4Q6Lobatto quadrature of order 6h E rr o r P2N4Q6 P3N4Q6 P4N4Q6 10 −2 −1 −10 −5 P1N5Q8Lobatto quadrature of order 8h E rr o r P2N5Q8 P3N5Q8 P4N5Q8 P5N5Q8 10 −2 −1 −10 −5 P1Lobatto quadrature of order 10h E rr o r P2 P3 P4 P5 P6 Figure 5.9: Kepler problem: Log-log plot of the error for the configurations q (crosses) and momenta p (dots)with step size h for the integrators ( P sN rQ r − Lob ); divided into four subplots, separated by appliedLobatto quadrature. The values of ( P s − N sQ s − Lob ) and ( P sN sQ s − Lob ) are superimposed, theyare of the same order. In the last subplot P s denotes ( P sN Q Lob ). P1 P2 P30246810 O r de r Lobatto quadrature of order 4 P1 P2 P3 P40246810 O r de r Lobatto quadrature of order 6 P1 P2 P3 P4 P50246810 O r de r Lobatto quadratur of order 8 P1 P2 P3 P4 P5 P60246810 O r de r Lobatto quadrature of order 10 Figure 5.10: Kepler problem: Order of the variational integrators( P sN rQ r − Lob ) with respect to theconfigurations q (left bar) and the momenta p (right bar). The title gives the order of the applied Lobattoquadrature which is 2 r − 2. The x -axis indicates the degree of the used polynomial. In this work variational integrators of higher order have been constructed by following the approach ofGalerkin variational integrators introduced in [MW01]. Thereby, the solution of the Euler-Lagrange equa-tions is approximated by a polynomial of degree s and the action by a quadrature formula based on r quadrature points. The restriction to r = s quadrature points (as assumed in [MW01]), which leads to sym-plectic partitioned Runge-Kutta methods, is dropped. For the resulting methods the order of convergencehas been determined numerically for two numerical examples. It has been numerically demonstrated thatthe order of convergence can be increased by adapting the number of quadrature points to the polynomialdegree. In particular, if the Lobatto quadrature is used, the choice of r = s + 1 leads to an integrator oforder 2 s which is of two orders higher compared to an integrator with r = s quadrature points (order 2 s − r and the polynomial degree s could bedetermined for different quadrature rules and the numerical results predict that 2 s is the maximal possibleorder of the constructed variational integrators. The structure-preserving properties such as symplecticity,momentum-preservation and good long-time energy behavior have been demonstrated by numerical exam-ples. In addition, for symmetrically distributed control points of the polynomial, the variational integrators( P sN rQ rGau ) and ( P sN rQ r − Lob ) have been shown to be time-reversible. Furthermore, a linear sta-22 10 20 30 40 50−2−1.5−1−0.500.511.522.53 x 10 −12 t A ngu l a r M o m en t u m E rr o r −6 t A ngu l a r M o m en t u m E rr o r P3N3Q6 h=0.25 P4N4Q8 h=0.25 P2N2Q4 h=0.25 RK4 Figure 5.11: Kepler problem: Left: Error of the angular momentum for the variational integrators( P N Gau ) (dots), ( P N Gau ) (stars) and ( P N Gau ) (squares) with step size h = 0 , 25. Right: Sameplot as on the left with non-variational integrator added (Runge-Kutta method of order four with step size h = 2 − (red solid)). −5 P4N4Q8; h=0.25t E ne r g y − E rr o r P3N3Q6; h=0.125RK4; h=0.015625 Figure 5.12: Kepler problem: Error of the energy for ( P N Q Gau ) with step size h = 0 . 125 (squares),( P N Q Gau ) with step size h = 0 . 25 (diamonds) and for a Runge-Kutta method of order four with stepsize h = 2 − (solid line).bility analysis has been performed for selected integrators and stability regions have been determined. It hasbeen shown that the integrators ( P sN sQ sGau ) are A-stable, i.e. there are no stability restrictions on thestep size h .In the future, a formal proof of the numerically determined convergence order of min (2 s, u ) has to be per-23ormed. To this end, techniques of variational error analysis [MW01, HL12] can be applied which are basedon the following main idea: Rather than considering how closely the trajectory of the discrete Hamiltonianmap matches the exact trajectory, one considers how closely the discrete Lagrangian matches the actionintegral. It is shown in [MW01] and [PC09] that both order concepts are equivalent. The analytical compu-tation of the variational error defined in this way for the general Galerkin variational integrators introducedin this work is still subject of ongoing research. The application of the constructed higher order variationalintegrators to holonomic and nonholonomic integrators is straightforward, however, a careful analysis hasto be carried out to see if the predicted orders also hold for these systems. Currently, the approach isused for the optimal control of mechanical systems ([CJOB12]) and numerical results confirm that also theadjoint resulting from the necessary optimality condition inherits the same order of the variational schemeas it also does for symplectic partitioned Runge-Kutta methods (cf. [OBJM11]). Another topic of interestis the construction of time-adaptive variational integrators since naive time-adapting strategies destroy thestructure-preserving properties (see e.g. [LR04]). The higher order integrators could be applied to adaptthe order rather than to adapt the step size, e.g. a higher order scheme can be deployed if higher accu-racy requirements have to be matched. Furthermore, for systems involving slow and fast time scales (seee. g. [LOB13]) variational integrators of different orders for the different subsystems can be used to increasethe efficiency of the simulations. Thereby, an investigation regarding preserved quantities and long-timebehavior is essential. References [Arn99] V. I. Arnold. Mathematical Methods of Classical Mechanics . Springer New York, 1999.[BRO08] N. Bou-Rabee and H. Owhadi. Stochastic variational integrators. IMA Journal on Numerical Analysis ,29:421–443, 2008.[CJOB12] C. M. Campos, O. Junge, and S. Ober-Bl¨obaum. Higher order variational time discretization of optimalcontrol problems. In ,Melbourne, Australia, 9-13 July 2012.[CM01] J. Cort´es and S. Mart´ınez. Non-holonomic integrators. Nonlinearity , 14(5):1365–1392, 2001.[FMOW03] R. C. Fetecau, J. E. Marsden, M. Ortiz, and M. West. Nonsmooth Lagrangian Mechanics and Varia-tional Collision Integrators. SIAM Journal on Applied Dynamical Systems , 2(3):381–416, 2003.[HL12] J. Hall and M. Leok. Spectral Variational integrators. (preprint, arXiv:1211.4534 ), 2012.[HLW02] E. Hairer, C. Lubich, and G. Wanner. Geometric numerical integration , volume 31 of Springer Seriesin Computational Mathematics . Springer, 2002.[HW10] E. Hairer and G. Wanner. Solving ordinary differential equations. II. Stiff and differential-algebraicproblems . Springer series in computational mathematics. Springer, Heidelberg, New York, 2010.[KMOW00] C. Kane, J. E. Marsden, M. Ortiz, and M. West. Variational integrators and the Newmark algorithmfor conservative and dissipative mechanical systems. International Journal for Numerical Methods inEngineering , 49(10):1295–1325, 2000.[KMS10] M. Kobilarov, J. E. Marsden, and G. S. Sukhatme. Geometric discretization of nonholonomic systemswith symmetries. Discrete and Continuous Dynamical Systems - Series S , 1(1):61–84, 2010.[Las88] F. M. Lasagni. Canonical Runge-Kutta methods. Zeitschrift f¨ur Angewandte Mathematik und PhysikZAMP , 39:952–953, 1988.[LMO08] S. Leyendecker, J. E. Marsden, and M. Ortiz. Variational integrators for constrained dynamical systems. Journal of Applied Mathematics and Mechanics , 88(9):677–708, 2008.[LMOW03] A. Lew, J. E. Marsden, M. Ortiz, and M. West. Asynchronous variational integrators. Archive forRational Mechanics and Analysis , 167:85–146, 2003.[LMOW04a] A. Lew, J. E. Marsden, M. Ortiz, and M. West. An overview of variational integrators. In L. P. Franca,T. E. Tezduyar, and A. Masud, editors, Finite Element Methods: 1970’s and Beyond , pages 98–115.CIMNE, 2004.[LMOW04b] A. Lew, J. E. Marsden, M. Ortiz, and M. West. Variational time integrators. International Journal forNumerical Methods in Engineering , 60(1):153–212, 2004. LOB13] S. Leyendecker and S. Ober-Bl¨obaum. A variational approach to multirate integration for constrainedsystems. In Jean-Claude Samin and Paul Fisette, editors, Multibody Dynamics , volume 28 of Compu-tational Methods in Applied Sciences , pages 97–121. Springer Netherlands, 2013.[LOBMO10] S. Leyendecker, S. Ober-Bl¨obaum, J. E. Marsden, and M. Ortiz. Discrete mechanics and optimal controlfor constrained systems. Optimal Control, Applications and Methods , 31(6):505–528, 2010.[LR04] B. Leimkuhler and S. Reich. Simulating Hamiltonian dynamics . Cambridge University Press, 2004.[LS11] M. Leok and T. Shingel. General Techniques for Constructing Variational Integrators. Frontiers ofMathematics in China , 2011.[MPS98] J. E. Marsden, G. W. Patrick, and S. Shkoller. Multisymplectic geometry, variational integrators, andnonlinear PDEs. Communication in Mathematical Physics , 199:351–395, 1998.[MQ02] R. I. McLachlan and G. R. W. Quispel. Splitting methods. Acta Numerica , 11:341–434, 2002.[MW01] J. E. Marsden and M. West. Discrete mechanics and variational integrators. Acta Numerica , 10:357–514,2001.[OBJM11] S. Ober-Bl¨obaum, O. Junge, and J. E. Marsden. Discrete mechanics and optimal control: an analysis. Control, Optimisation and Calculus of Variations , 17(2):322–352, 2011.[OBTC + 13] S. Ober-Bl¨obaum, M. Tao, M. Cheng, H. Owhadi, and J. E. Marsden. Variational integrators for electriccircuits. Journal of Computational Physics , 242(C):498–530, 2013.[PC09] G. W. Patrick and C. Cuell. Error analysis of variational integrators of unconstrained Lagrangiansystems. Numerische Mathematik , 113:243–264, 2009.[Rei94] S. Reich. Momentum conserving symplectic integrations. Physica D , 76(4):375–383, 1994.[Saa12] N. Saake. Konstruktion und Analyse variationeller Integratoren h¨oherer Ordnung. Diploma thesis,Paderborn, 2012.[SG09] A. Stern and E. Grinspun. Implicit-explicit variational integration of highly oscillatory problems. SIAMMultiscale Modeling and Simulation , 7:1779–1794, 2009.[SS88] J. M. Sanz-Serna. Runge-Kutta schemes for Hamiltonian systems. BIT Numerical Mathematics , 28:877–883, 1988.[Sun93] G. Sun. Symplectic partitioned Runge-Kutta methods. Journal of Computational Mathematics ,11(4):365–372, 1993.[Sur89] Y. B. Suris. The canonicity of mappings generated by Runge-Kutta type methods when integrating thesystems ¨ x = − ∂u/∂x . U.S.S.R. Computational Mathematics and Mathematical Physics , 29(1):138–144,1989.[Sur90] Y. B. Suris. Hamiltonian methods of Runge-Kutta type and their variational interpretation. Mathe-matical Modelling , 2:78–87, 1990.[TOM10] M. Tao, H. Owhadi, and J. E. Marsden. Nonintrusive and structure preserving multiscale integration ofstiff ODEs, SDEs, and Hamiltonian systems with hidden slow dynamics via flow averaging. MultiscaleModeling and Simulation , 8(4):1269–1324, 2010.[Yos90] H. Yoshida. Construction of higher order symplectic integrators. Physics Letters A , 150:262 – 268,1990. − − − − − − − P N Q R un − T i m e Error P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P N Q P s N Q P s N Q P s N Q P s N Q Figure 5.13: Harmonic oscillator: Global error against run-time for the integrators ( P sN rQ r − Lob ) anddifferent step sizes h . 26 − − − − − − − P N Q R un − T i m e Error P N Q P N Q P N Q P N Q P N Q G au P N Q G au P N Q G au P N Q G au P N Q G au Figure 5.14: Harmonic oscillator: Global error against run-time for the integrators ( P sN sQ sGau ) (dottedline) and ( P sN s + 1 Q sLob ) (solid line) for s = { , ..., } and different step sizes hh