Projected Variational Integrators for Degenerate Lagrangian Systems
Michael Kraus


=1 

di-men-sio-nal phase-space La-gran-gian La-grange Ha-mil-to-nian Ha-mil-ton sub-ma-ni-fold ma-ni-fold


















arrows
matrix



left=2.5cm,right=2.5cm,top=2.5cm,bottom=2.5cm






ifundefinedchapter
   theoremTheorem[section]
   theoremTheorem[chapter]


axiom[theorem]Axiom
claim[theorem]Claim
comment[theorem]Comment
conjecture[theorem]Conjecture
corollary[theorem]Corollary
definition[theorem]Definition
exerciseExercise
*exercise*Exercise
example[theorem]Example
lemma[theorem]Lemma
problem[theorem]Problem
proposition[theorem]Proposition
remark[theorem]Remark
summary[theorem]Summary








    
    
    
    
    



    δ÷÷OT1pzcmitProjected Variational Integrators for
    Degenerate Lagrangian Systems
        Michael Kraus
    (michael.kraus@ipp.mpg.de)
    
    Max-Planck-Institut für Plasmaphysik
    Boltzmannstraße 2, 85748 Garching, Deutschland
    
    
    and
    
    Technische Universität München, Zentrum Mathematik
    Boltzmannstraße 3, 85748 Garching, Deutschland
    
    
    
        August 24, 2017
    =============================================================================================================================================================================================================================================
    
    
    
    We propose and compare several projection methods applied to variational integrators for degenerate Lagrangian systems, whose Lagrangian is of the form L = ϑ() · - H() and thus linear in velocities.
    While previous methods for such systems only work reliably in the case of ϑ being a linear function of , our methods are long-time stable also for systems where ϑ is a nonlinear function of .
    We analyse the properties of the resulting algorithms, in particular with respect to the conservation of energy, momentum maps and symplecticity.
    In numerical experiments, we verify the favourable properties of the projected integrators and demonstrate their excellent long-time fidelity.
    In particular, we consider a two-dimensional Lotka–Volterra system, planar point vortices with position-dependent circulation and guiding centre dynamics.
    
    
    
    § INTRODUCTION
    
    
    In various areas of physics, we are confronted with degenerate Lagrangian systems, whose Lagrangian is of the form L = ϑ() · - H(), that is L is linear in the velocities .
    Examples for such systems are planar point vortices, certain Lotka–Volterra models and guiding centre dynamics.
    
    In order to derive structure-preserving integrators <cit.> for such systems, it seems natural to apply the variational integrator method <cit.>.
    This, however, does not immediately lead to stable integrators, as the resulting numerical methods will in general be multi-step methods which are subject to parasitic modes <cit.>. Moreover, we face the initialization problem, that is how to determine initial data for all previous time steps used by the method without introducing a large error into the solution.
    
    A potential solution to the first problem is provided by the discrete fibre derivative while a solution to the second problem is provided by the continuous fibre derivative. Using the discrete fibre derivative, we can rewrite the discrete Euler–Lagrange equations in position-momentum form, which constitutes a one-step method for numerically computing the phasespace trajectory in terms of the generalized coordinates  together with their conjugate momenta . The resulting system can be solved, as in general the discrete Lagrangian will not be degenerate even though the continuous Lagrangian is.
    The continuous fibre derivative ( (t),  (t) = ∂ L / ∂ ( (t),  (t))) can then be used in order to obtain an initial value for the conjugate momenta as functions of the coordinates as (t) = ϑ ( (t)).
    Tyranowski et al. <cit.> show that this is a viable strategy when ϑ is a linear function. Unfortunately, for the general case of ϑ being a nonlinear function, this idea does in general not lead to stable integrators as the numerical solution will drift away from the constraint submanifold defined by the continuous fibre derivative, ϕ((t), (t)) =  (t) - ϑ ( (t)) = 0. 
    
    A standard solution for such problems is to project the solution back to the constraint submanifold after each time step <cit.>. This, however, renders the integrator non-symmetric (assuming the variational integrator itself is symmetric), which leads to growing errors in the solution and consequently a drift in the total energy of the system.
    Improved long-time stability is achieved by employing a symmetric projection <cit.>, where the initial data is perturbed away from the constraint submanifold before the variational integrator is applied and then projected back to the manifold.
    While these projection methods are standard procedures for holonomic and non-holonomic constraints, there are only few traces in the literature on their application to Dirac constraints ϕ(,)=0.
    Some authors consider general differential algebraic systems of index two <cit.>, but do not go into the details of Lagrangian or Hamiltonian systems. 
    Or they consider symplectic integrators for Hamiltonian systems subject to holonomic constraints <cit.>, which simplifies the situation dramatically compared to the case of Dirac constraints.
    Most importantly, we are not aware of any discussion of the influence of such a projection on the symplecticity of the algorithm assuming that the underlying numerical integrator is symplectic. As symplecticity is a crucial property of Lagrangian and Hamiltonian systems, which is often important to preserve in a numerical simulation, we will analyse all of the proposed methods regarding its preservation.
    We will find that the well-known projection methods, both standard projection and symmetric projection, are not symplectic. However, we can introduce small modifications to the symmetric projection method which make it symplectic.
    
    
    
    The outline of the paper is as follows.
    In Section <ref> we provide an overview of degenerate Lagrangian systems and Dirac constraints and their various formulations and discuss symplecticity and momentum maps.
    In Section <ref> we review the discrete action principle leading to variational integrators and problems that arise when this method is applied to degenerate Lagrangians.
    This is followed by a discussion of the proposed projection methods in Section <ref> and numerical experiments in Section <ref>.
    
    
    
    
    
    § DEGENERATE LAGRANGIAN SYSTEMS
    
    
    Degenerate Lagrangian systems have attracted quite some interest in the geometric mechanics literature <cit.> due to their interesting properties.
    They are also relevant for practical applications like the study of population models, point vortex dynamics or reduced charged particle models like the guiding centre system.
    In the following, we will consider degenerate Lagrangian systems characterized by a Lagrangian that is linear or singular in the velocities. In particular, we consider the class of systems whose Lagrangian is of the form
    
        L (, ) = ϑ() · - H() .
    
    The Lagrangian L is a function on the tangent bundle M,
    
        L : M→ ,
    
    where M denotes the configuration manifold of the system which is assumed to be of dimension d. The cotangent bundle of the configuration manifold M is denoted by M.
    Further, we denote the coordinates of a point m ∈M by (m) = (^1 (m), , ^d (m)) and similarly coordinates of points in M by (^i, ^i) and coordinates of points in M by (^i, ^i).
    In the following, we will always assume the existence of a global coordinate chart, so that M can be identified with the Euclidean space ^d. For simplicity, we often use short-hand notation where we write  to refer to both a point in M as well as its coordinates. Similarly, we often denote points in the tangent bundle M by (, ).
    In local coordinates, the Lagrangian (<ref>) is thus written as a map (, ) ↦ L(, ).
    
    In Equation (<ref>), ϑ = ϑ_i ()  ^i is a differential one-form ϑ : M→M, whose components ϑ_i : M→ are general, possibly nonlinear functions of , some of which (but not all) could be identically zero.
    For details on differential forms, tangent and cotangent bundles the interested reader may consult any modern book in mathematical physics or differential geometry. We recommend <cit.> for more physics oriented accounts and <cit.> for more mathematics oriented accounts. In the following we assume a basic understanding of these concepts. To see their usefulness for classical mechanics we refer to <cit.>.
    
    
    
    
     §.§ Hamilton's Action Principle
    
    
    The evolution of Lagrangian systems is described by curves  on M.
    To make this precise, let us fix two points _1, _2∈M and an interval [t_1, t_2] ⊂ and define the path space connecting _1 and _2 as
    
        Q ( t_1, t_2, _1, _2 ) = { : [t_1, t_2] →M |  is a C^2 curve (t_1) = _1,  (t_2) = _2} .
    
    Elements  of Q ( t_1, t_2, _1, _2 ) map the time interval [t_1, t_2] to curves on M, whereby the first and last points, (t_1) and (t_2), take fixed values, _1 and _2, respectively.
    Such a curve : [t_1, t_2] →M with (t) = (q^1 (t), , ^d (t)) can be lifted to a curve : [t_1, t_2] →M, which in coordinates is given by
    
        (t) = ( ^1 (t), , ^d (t), d^1dt (t), , d^ddt (t) ) .
    
    In the following the derivative of the curve with respect to the parameter t is denoted by = d/dt.
    This constitutes slight abuse of notation as we also denote the tangent bundle coordinates that way (note that not all curves in the tangent bundle are lifts of curves in the configuration manifold), but it should be clear at any time if we refer to the derivative of the curve or to the coordinates.
    
    
    
    In order to determine the equations of motion of a Lagrangian system, we consider infinitesimal variations of the action integral
    
        A [] = ∫_0^T L ( (t),  (t))   dt ,
    
    where, without loss of generality and in order to simplify the discrete treatment, we choose an interval [0,T].
    Infinitesimal variations (c.f., Figure <ref>) are defined in terms of C^2 maps : [0,T] →^d which vanish at the boundaries of the interval [0,T], that is (0) = 0 and (T) = 0, and are such that c(t) ∈[q(t)]M for 0 ≤ t ≤ T with [q(t)]M the tangent space to M at the point q(t).
    In coordinates, we can write
    
        (t) = ( ^1 (t), , ^d (t) ) .
    
    We can now define a one-parameter family of trajectories ^∈Q ( 0, T, _1, _2 ), for which
    
        = dd^|_ = 0 ,
    
    with ∈ (-r,+r) and r ∈_+. That is ^ is a differentiable mapping (-r,+r) × [0,T] →M, such that ^ (0) = q (0) and ^ (T) = q (T) for all values of  and q^0 (t) = q(t) for all values of t.
    If Q is a vector space, the simplest example of such a family of trajectories is given by
    
        ^ (t) = ( ^1 (t) + ^1 (t), , ^d (t) + ^d (t) ) .
    
    On a general manifold, the corresponding expressions are usually more complicated.
    However, if the manifold is embedded in ^d, each member of the family of trajectories can be expanded in a power series, whose leading-order terms are those given in (<ref>).
    
    Hamilton's principle states that in order to determine the equations of motion, we need to find a curve = ^0, such that the action (<ref>) takes a stationary point with respect to all curves ^.
    A necessary condition for  making the action stationary is that 
    
        ddA [^] |_ = 0   = dd∫_0^T L (^ (t), ^ (t))   dt |_ = 0
           = ∫_0^T[ ∂ L∂(  (t),  (t) ) · (t) + ∂ L∂(  (t),  (t) ) · (t) ]   dt = 0 .
    
    Let us note that the derivatives of L with respect to  and  are sometimes also written as D_1 and D_2, respectively, where D_i denotes the slot-derivative with respect to the ith argument of L.
    Assuming that the operations of computing variations of  and computing the time derivative of  commute (which is a fair assumption, see e.g. <cit.>), we can integrate the second term of (<ref>) by parts to obtain
    
        ddA [q^] |_ = 0
        = ∫_0^T[ ∂ L∂(  (t),  (t) ) - ddt( ∂ L∂(  (t),  (t) ) ) ] · (t)   dt 
        
         + [ ∂ L∂(  (t),  (t) ) · (t) ]_0^T = 0 .
    
    As the infinitesimal variations  are required to vanish at the boundaries, the boundary term vanishes, and as the variations  are otherwise arbitrary, vanishing of the variations of A implies vanishing of the term in square brackets in the integrand.
    This leads us to the Euler–Lagrange equations, that is the equations of motion,
    
        ∂ L∂(  (t),  (t) ) - ddt( ∂ L∂(  (t),  (t) ) ) = 0 .
    
    For the Lagrangian (<ref>) the Euler–Lagrange equations yield
    
        ∇ϑ ( (t)) · (t) - ∇ H( (t)) - ϑ̇ ( (t)) = 0 ,
    
    which after computing the time-derivative of ϑ can be written as
    
        Ω^T ( (t))   (t)    = ∇ H( (t))       with   Ω_ij   = ∂ϑ_j∂^i - ∂ϑ_i∂^j .
    
    The skew-symmetric matrix Ω plays an important role as it holds the components of the symplectic form ω on M. 
    Let us note that in principle Ω can be of odd dimension, in which case the corresponding two-form ω is degenerate and therefore does not classify as a symplectic form. However, most degenerate Lagrangians of the form (<ref>), including all of the examples discussed in Section <ref>, originate from some kind of coordinate transformation of a canonical system, possibly followed by some reduction procedure, which always results in a system of even degree. We therefore assume in the following that the system under consideration is of even dimension and hence symplectic.
    Details of the symplectic structure will be discussed in Section <ref>.
    
    Equation (<ref>) has the structure of a noncanonical Hamiltonian system on M, characterized by the skew-symmetric matrix Ω and the Hamiltonian H.
    For such noncanonical Hamiltonian systems no general geometric integrators are known.
    In principal it is possible to use the Darboux theorem to find canonical coordinates and then apply some canonical symplectic integrator.
    In practice, however, the construction of such Darboux coordinates tends to be a non-trivial task and is often possible only locally but not globally.
    Our strategy will thus be to reformulate the system as a canonical Hamiltonian system by adding canonical conjugate momenta, thus doubling the size of the solution space, and restricting the numerical solution to the physical subspace of this extended solution space.
    The geometrical foundation of this procedure is the theory of Dirac constraints.
    
    
    
    
     §.§ Dirac Constraints
    
    
    Degenerate systems of the form (<ref>) can also be formulated in terms of the phasespace trajectory (, ) in the cotangent bundle M, subject to a primary constraint in the sense of Dirac, determined by the function ϕ : M→^d, given by
    
        ϕ (, ) =  - ϑ() = 0 ,
    
    and originating from the fibre derivative F L : M→M,
    
        F L (_) ·_ = dd|_=0 L(_ + _) ,
    
    where _ = (, ) and _ = (, ) denote two points in M which share the same base point  and are thus elements of the same fibre of M.
    By acting point-wise for each t, the fibre derivative maps the curve (, ) in the tangent bundle M into the curve (, ) in the cotangent bundle M,
    
        ( (t),  (t)) = (  (t), ∂ L∂ ( (t),  (t)) ) = ( (t), ϑ ( (t))) ,
    
    where the last equality follows for Lagrangians of the form (<ref>).
    The Dirac constraint arising from the degenerate Lagrangian restricts the dynamics to the submanifold
    
        Δ = { (, ) ∈M | ϕ (, ) = 0  }⊂M .
    
    In the preceding and the following, we assume that the Lagrangian is degenerate in all velocity components, that is, the Lagrangian is either linear or singular in each component of , so that
    
        ∂^2 L∂^i ∂^j = 0
        for all 1 ≤ i,j ≤ d.
    
    For instructive reasons, however, assume for a moment that the Lagrangian is degenerate in only m<d components of  and, e.g., quadratic in the other d-m components. That is to say we can write
    
        (t)
        = ( β_1 ( (t),  (t)), , β_d-m ( (t),  (t)),  ϑ_d-m+1 ( (t)), , ϑ_d ( (t)) )^T ,
    
    where
    
        ∂ L∂^i ( (t),  (t)) = β_i ( (t),  (t))     1 ≤ i ≤ d - m , 
        ϑ_i ( (t))     d-m+1 ≤ i ≤ d .
    
    We can then denote coordinates in Δ by (^i, π^j) with 1 ≤ i ≤ d and 1 ≤ j ≤ d-m, where the π^i denote those momenta which are “free”, i.e., not determined by the Dirac constraint. The inclusion map i : Δ→M can then be written as
    
        i : (, π) ↦ (, π, ϑ ()) .
    
    In the fully degenerate case, however, we have m=d, so that the configuration manifold M and the constraint submanifold Δ are isomorphic and we can label points in Δ by the same  we use to label points in M.
    The inclusion map i : Δ→M simplifies accordingly and reads
    
        i : ↦ (, ϑ ()) ,
    
    where it is important to keep in mind that  denotes a point in Δ.
    The inverse operation is given by the projection π_Δ : M→Δ, defined such that π_Δ∘ i = 𝕀.
    
    As we are lacking a general framework for constructing structure-preserving numerical algorithms for noncanonical Hamiltonian systems on M, we will construct such algorithms on i(Δ). This can be achieved by using canonically symplectic integrators on M and assuring that their solution stays on i(Δ). To this end we will employ various projection methods, as discussed in Section <ref>.
    
    
    
    
     §.§ Augmented Hamiltonian Approach
    
    
    Hamilton's form of the equations of motion for a degenerate Lagrangian system (<ref>) can be derived from the phasespace action
    
        A [q,p,λ] = ∫_0^T[  (t) · (t) - H ( (t),  (t),  (t)) ] dt ,
    
    with the augmented Hamiltonian H : M×^d→ defined as
    
        H (, , λ) = H() + ϕ (,) · .
    
    Applying Hamilton's principle of stationary action to (<ref>) results in the following index two differential-algebraic system of equations (see e.g. <cit.> for a definition of the notion of index),
    
        (t)    = ϕ_^T ( (t),  (t))   (t) , 
         (t)    = - H_ ( (t)) - ϕ_^T ( (t),  (t))   (t) , 
        
        0    = ϕ ( (t),  (t)) .
    
    Here, subscripts  and  denote partial derivatives with respect to the coordinates of M. For the constraint function ϕ these derivatives are explicitly written as
    
        ϕ_   = [ ∂ϕ_1∂^1 ∂ϕ_1∂^2         ∂ϕ_1∂^d;       ⋮       ⋮       ⋱       ⋮; ∂ϕ_d∂^1 ∂ϕ_d∂^2         ∂ϕ_d∂^d;         ]      and   ϕ_   = [ ∂ϕ_1∂^1 ∂ϕ_1∂^2         ∂ϕ_1∂^d;       ⋮       ⋮       ⋱       ⋮; ∂ϕ_d∂^1 ∂ϕ_d∂^2         ∂ϕ_d∂^d;         ] ,
    
    so that in ϕ_^Tλ and ϕ_^Tλ, the components of λ are contracted with the components of ϕ, not with the derivatives.
    As the Hamiltonian H does not depend on  and ϕ =  - ϑ(), we find that the first equation reduces to (t) =  (t), that is the Lagrange multiplier takes the role of the velocity.
    Denoting the trajectory in the cotangent bundle M by =(, ), the equations of motion (<ref>) can be rewritten more compactly as
    
        (t)    = Ω^-1 ∇ H ( (t)) + Ω^-1 ∇ϕ^T ( (t))   (t) , 
        
        0    = ϕ ( (t)) ,
    
    where ∇ denotes the derivatives with respect to = (, ).
    
    
    
    
     §.§ Hamilton–Pontryagin Principle
    
    
    The phasespace action principle of the previous section is equivalent to the Hamilton–Pontryagin principle <cit.> on M⊕M, given by
    
        δ∫_0^T[ L ( (t),  (t)) +  (t) · (  (t) -  (t) ) ]   dt = 0 .
    
    Here, the dynamics of the system are described by the evolution of (, , ), which constitute a trajectory in the Pontryagin bundle M⊕M.
    If the Lagrange multiplier  is replaced by the velocity  it is easy to verify that the Lagrangian is related to the augmented Hamiltonian (<ref>) by L (, ) = · - H (, , ), and that the Hamilton–Pontryagin principle (<ref>) is equivalent to the phasespace action principle δA (, , ) = 0 with the augmented action A given in (<ref>).
    Computing variations of (<ref>), where ,  and  are all varied independently and only restricted in that the variations of  have to vanish at the endpoints, we obtain the implicit Euler–Lagrange equations,
    
        (t)    =  (t) ,     (t)    = ∂ L∂v ( (t),  (t)) ,     (t)    = ∂ L∂ ( (t),  (t)) ,
    
    which are easily seen to be equivalent to (<ref>) and (<ref>).
    Here, the Dirac constraint ϕ( (t),  (t)) = 0 appears quite naturally as one of the equations of motion, which suggests that the Hamilton–Pontryagin principle might be the natural starting point for the discretization of degenerate Lagrangian systems.
    That this is not necessarily the case will be discussed in Section <ref>.
    
    
    
    
     §.§ Symplecticity
    
    
    Our aim is to construct methods which retain the symplecticity of the integrator as well as its momentum maps.
    Care has to be taken, when stating that the variational integrator and the projection are symplectic.
    The continuous system preserves two symplectic forms, the canonical two-form ω on Δ, but also the noncanonical two-form ω on M defined by
    
        ω   = ϑ = 12Ω_ij ^i∧^j      with   Ω_ij   = ∂ϑ_j∂^i - ∂ϑ_i∂^j .
    
    The matrix Ω is the noncanonical symplectic matrix which we already encountered in the equations of motion (<ref>).
    The function ϑ is interpreted as a one-form on M, in coordinates given by ϑ = ϑ_i ()  ^i.
    In principle it is possible that ω is degenerate, namely in the case of a system of odd dimension d. Then ω is not a symplectic form but a presymplectic form. Most of the following discussion also holds in this case. However, in almost all examples of practical relevance the configuration space is even-dimensional. For this reason, we will always refer to ω as symplectic form.
    
    Note that ω is not the symplectic form ω_L on M originating from the boundary terms in the action principle (<ref>).
    Besides leading to the equations of motion, the variational principle provides a direct and natural way to derive the fundamental geometric structures of classical mechanics.
    For this derivation, the boundary conditions (0) =  (T) = 0 are relaxed, while the time interval [0,T] is kept fixed.
    Thus the variational principle reads
    
        ddA [^] |_ = 0   = ∫_0^T[ ∂ L∂(  (t),  (t) ) · (t) + ∂ L∂(  (t),  (t) ) · (t) ]   dt  + [ ∂ L∂·]_t_1^t_2 .
    
    where the variations (t) do not vanish at the boundary point, so that the last term on the right hand side does not vanish.
    This last term corresponds to a linear pairing of the function ∂ L / ∂, which in general is a function of (, ), with the tangent vector .
    The boundary term in (<ref>) can be written as Θ_Lδ q|_0^T, where Θ_L is the so-called Lagrangian one-form or Cartan one-form, in coordinates given by
    
        Θ_L = ∂ L∂^i ^i ,
    
    One could be tempted to regard ∂ L / ∂ as a one-form on M as it only has a component in . The same way  could be regarded as a tangent vector on M.
    However, in general ∂ L / ∂ is a function of (, ) and therefore clearly a function on M.
    The exterior derivative of the Lagrangian one-form gives the Lagrangian two-form, also referred to as the symplectic two-form,
    
        ω_L = Θ_L ,
    
    given in coordinates by
    
        ω_L   = ∂^2 L∂^i ∂^j  d ^i∧ d ^j
         + ∂^2 L∂^i ∂^j  d ^i∧ d ^j .
    
    For more details on this derivation see e.g. <cit.>.
    As the Lagrangian (<ref>) is degenerate, so is the corresponding symplectic matrix Ω_L, which can be written in block form as
    
        Ω_L =
        [ Ω 0; 0 0;   ] .
    
    We recognize the upper left block, which corresponds to the noncanonical symplectic matrix Ω on M. When we discuss symplecticity in the following, we are always referring to the noncanonical symplectic form ω or its matrix representation Ω.
    
    Preserving the noncanonical symplectic form ω on M is equivalent to preserving the canonical symplectic form ω = Θ on the embedding of Δ in M.
    Denoting coordinates on M by = (, ), the canonical one-form Θ and the symplectic two-form ω can be written in coordinates as
    
        Θ = p_i  q^i ,     ω = 12Ω_ij ^i∧^j = _i∧^i ,
    
    with Ω the canonical symplectic matrix, given by
    
        Ω = [    0 -  1;    1    0 ] .
    
    On the constraint submanifold we have that _i = ϑ_i (), and therefore _i = ϑ_i (), so that ω restricted to Δ reads
    
        ω|_Δ
         = ∂ϑ_i∂^j ^j∧^i   = 12∂ϑ_i∂^j ^j∧^i
         + 12∂ϑ_j∂^i ^i∧^j
           = 12( ∂ϑ_j∂^i - ∂ϑ_i∂^j) ^i∧^j
         = ω .
    
    Using the inclusion (<ref>), we can write ω = i^*ω.
    The preceding arguments thus suggest that in order to construct a numerical algorithm that preserves the noncanonical symplectic form ω on M, a viable strategy could be to construct a canonically symplectic algorithm on M whose solution stays on the constraint submanifold Δ.
    
    
    
    
     §.§ Noether Theorem and Conservation Laws
    
    
    On of the most influential results of classical mechanics in the 20th century is the correspondence of point-symmetries of the Lagrangian and conservation laws of the Euler–Lagrange equations established by Emmy Noether <cit.>.
    In the following we will summarize her famous theorem.
    
    Consider a Lagrangian system L : M→ and a one-parameter group of transformations {σ^ : ∈ B_0^r ,  σ_0 = 𝕀}, where B_0^r denotes the open ball with radius r > 0 centred at 0.
    We denote the transformed trajectory by ^ = σ^∘ and its time derivative by ^ = d (σ^∘) / dt such that ^0 = q and ^0 =.
    We have a symmetry if the transformation σ^ leaves the Lagrangian L invariant, that is
    
        L ( ^ (t), ^ (t) ) = L (  (t),  (t) )
        for all  and all .
    
    Taking the  derivative of (<ref>), we obtain the infinitesimal invariance condition,
    
        dd L ( ^ (t), ^ (t) ) |_=0 = 0 ,
    
    which is equivalent to (<ref>).
    Explicitly computing this  derivative, we obtain
    
        dd L ( ^ (t), ^ (t) ) |_=0   = ∂ L∂(  (t),  (t) ) ·d σ^d|_=0 + ∂ L∂(  (t),  (t) ) ·d σ̇^d|_=0 = 0 .
    
    Denoting by V the vector field with flow σ^, defined as follows,
    
        V    = d σ^d|_=0 ,
    
    and if  solves the Euler–Lagrange equations (<ref>), we can rewrite (<ref>) as
    
        [ ddt∂ L∂(  (t),  (t) ) ] · V (  (t) ) + ∂ L∂(  (t),  (t) ) ·dVdt(  (t),  (t) )    = 0 .
    
    The time derivative of the vector field V is simply computed by the chain rule, so that assuming that the transformation σ^ does not explicitly depend on time it is given by
    
        dVdt   = ^j ∂ V^i∂^j∂∂^i .
    
    The expression in (<ref>) amounts to a total time derivative of the so-called Noether current, which constitutes the preserved quantity,
    
        ddt[ P (  (t),  (t) ) ]    = 0 ,    
        P  (  (t),  (t) )    = ∂ L∂(  (t),  (t) ) · V ( (t)) .
    
    Thus, the momentum ∂ L / ∂ in the direction V is conserved along solutions  of the Euler–Lagrange equations (<ref>) obtained from L for all times t.
    In Section <ref>, we will apply the Noether theorem several times in order to determine the conservation laws for the various examples we will consider.
    
    
    
    
    § VARIATIONAL INTEGRATORS
    
    
    Variational integrators can be seen as the Lagrangian equivalent of symplectic integrators for Hamiltonian systems. Instead of discretizing the equations of motion, the action integral is discretized, followed by the application of a discrete version of Hamilton's principle of stationary action. This leads to discrete Euler–Lagrange equations (the discrete equations of motion) at once. 
    The evolution map that corresponds to the discrete Euler–Lagrange equations is what is called a variational integrator. Such a numerical scheme preserves a discrete symplectic form which originates from the boundary terms in the variation of the discrete action.
    
    The seminal work in the development of a discrete equivalent of classical mechanics was presented by <cit.>.
    His method, based on a discrete variational principle, leads to symplectic integration schemes that automatically preserve constants of motion <cit.>.
    Comprehensive reviews of variational integrators and discrete mechanics can be found in <cit.> and <cit.>, including thorough accounts on the historical development preceding and following the work of .
    
    In the following we collect some material on variational integrators, specifically on the discrete action principle, the position-momentum form, and on variational Runge–Kutta methods, before discussing the problems that arise when trying to apply the method to degenerate Lagrangians.
    
    
    
    
     §.§ Discrete Action Principle
    
    
    Time will be discretized uniformly, i.e., the time step h ∈_+ is constant. We thus split the interval [0,T] into a finite sequence of times { t_n = nh  |  n = 0, , N }, where h = T/N, so that t_N = T.
    Let us denote the configuration of the discrete system at time t_n by _n, so that _n≈ (t_n), where (t_n) is the configuration of the continuous system at time t_n.
    Then a discrete trajectory can be written as _d = {_n}_n=0^N.
    
    The discrete Lagrangian is defined as an approximation of the time integral of the continuous Lagrangian over the interval I_n = (t_n, t_n+1), i.e.,
    
        L_d (_n, _n+1) ≈∫_I_n L ( _n,n+1 (t) ,  _n,n+1 (t) )   dt ,
    
    where _n,n+1 denotes the solution of the Euler–Lagrange equations (<ref>) in I_n.
    The specific expression of the discrete Lagrangian is determined by the polynomial approximation of the trajectory and the quadrature rule used to approximate the integral.
    The discrete action then becomes merely a sum over the time index of discrete Lagrangians,
    
        A_d [_d] = ∑_n=0^N-1 L_d (_n, _n+1) .
    
    Using linear interpolation between _n and _n+1 to describe the discrete trajectory, thus approximating _n,n+1 by
    
        _n,n+1 (t)    ≈_n t_n+1 - tt_n+1 - t_n + _n+1 t - t_nt_n+1 - t_n ,
    
    the velocity _n,n+1 will be approximated by a simple finite-difference expression, namely
    
        _n,n+1 (t)    ≈_n+1 - _nt_n+1 - t_n .
    
    As we assume the time step to be constant, in the following we will just write h instead of t_n+1 - t_n.
    The quadrature approximating the integral in (<ref>) is most often realized by either the trapezoidal rule, leading to the discrete Lagrangian
    
        L_d^tr (_n, _n+1) = h2  L ( _n, _n+1 - _nh) + h2  L ( _n+1, _n+1 - _nh) ,
    
    or the midpoint rule, leading to the discrete Lagrangian
    
        L_d^mp (_n, _n+1) = h   L ( _n + _n+12, _n+1 - _nh) .
    
    The configuration manifold of the discrete system is still M, but the discrete state space is M×M instead of M,
    such that the discrete Lagrangian L_d is a function
    
        L_d : M×M→ .
    
    
    The discrete equations of motion are determined in the same way as the continuous equations of motion (<ref>), that is by applying Hamilton's principle of stationary action.
    Infinitesimal variations of the discrete trajectories _d are given in terms of maps _d : { t_n}_n=0^N→^d, which vanish at t_0 and t_N, that is _d(t_0) = 0 and _d(t_N) = 0, and are such that _d(t_n) ∈[q_d (t_n)]M, where _d(t_n) = _n and q_d (t_n) = q_n and we denote such maps by _d = {_n}_n=0^N.
    Discrete one-parameter families of trajectories _d^ = {_n^}_n=0^N, are defined by
    
        _d = dd_d^|_ = 0 ,
    
    with the simplest example (c.f., Figure <ref>) given by
    
        _d^ = {_n + _n }_n=0^N .
    
    Such trajectories are elements of the discrete path space, defined as
    
        Q_d (t_0, t_N, _0, _N) = {_d : { t_n}_n=0^N→M | _d (t_0) = _0,  _d (t_N) = _N} .
    
    A necessary condition for the discrete trajectory _d = _d^0 making the discrete action (<ref>) stationary with respect to all curves ^, is that 
    
        ddA_d [_d^] |_ = 0
        = dd∑_n=0^N-1 L_d (_n^, _n+1^) |_ = 0
        = 0 .
    
    Computing the -derivative of the discrete action explicitly, we obtain
    
        ddA_d [_d^] |_ = 0   = ∑_n=0^N-1[ D_1 L_d (_n, _n+1) ·_n + D_2 L_d (_n, _n+1) ·_n+1] ,
    
    where D_i denotes the slot-derivative with respect to the ith argument of L_d.
    What follows corresponds to a discrete integration by parts, i.e., a reordering of the summation.
    The n=0 term is removed from the first part of the sum and the n=N-1 term is removed from the second part,
    
        ddA_d [_d^] |_ = 0   = D_1 L_d (_0, _1) · c_0 + ∑_n=1^N-1 D_1 L_d (_n, _n+1) ·_n
           + ∑_n=0^N-2 D_2 L_d (_n, _n+1) ·_n+1 + D_2 L_d (_N-1, _N) ·_N .
    
    As the variations at the endpoints, _0 and _N, vanish, the corresponding terms in the above sum also vanish.
    At last, the summation range of the second sum is shifted upwards by one with the arguments of the discrete Lagrangian adapted correspondingly, so that
    
        ddA_d [_d^] |_ = 0   =  ∑_n=1^N-1[ D_1 L_d (_n, _n+1) + D_2 L_d (_n-1, _n) ] ·_n .
    
    Hamilton's principle of least action requires the variation of the discrete action A_d to vanish for any choice of _n. Consequently, the expression in the square brackets of (<ref>) has to vanish. This defines the discrete Euler–Lagrange equations
    
        D_1 L_d (_n, _n+1) + D_2 L_d (_n-1, _n) = 0 .
    
    Denoting coordinates on M×M by (_1^1, , _1^d, _2^1, , _2^d), this can also be written as
    
        ∂ L_d∂_1 (_n, _n+1) + ∂ L_d∂_2 (_n-1, _n) = 0 .
    
    The discrete Euler–Lagrange equations (<ref>) define an evolution map
    
        F_L_d  :  M×M→M×M  :   ( _n-1, _n ) ↦ ( _n, _n+1 ) .
    
    Starting from two configurations, _0≈ (t_0) and _1≈ (t_1 = t_0 + h), the successive solution of the discrete Euler–Lagrange equations (<ref>) for _2, _3, etc., up to _N, determines the discrete trajectory _d.
    
    For the class of degenerate Lagrangians (<ref>) under consideration, the prescription of two sets of initial conditions, _0 and _1, is not natural. The continuous Euler–Lagrange equations are ordinary differential equations of first order and therefore need only one set of initial conditions in order to solve the equations. In practice, we face the problem that there is no unique way of determining a second set of initial conditions. All methods will introduce some error that will propagate to the solution and eventually most often lead to a break down of the solution.
    
    
    
    
     §.§ Position-Momentum Form
    
    
    A viable way around this problem appears to be to use the discrete fibre derivative to reformulate the discrete Euler–Lagrange equations (<ref>) in position-momentum form. For regular Lagrangians L this is equivalent to rewriting the continuous Euler–Lagrange equations in the form of Hamilton's equations by using the continuous fibre derivative and Legendre transform.
    
    Given a point (_n, _n+1) in M×M and a discrete Lagrangian L_d : M×M→, we define two discrete fibre derivatives, F^- L_d and F^+ L_d, in analogy to the continuous case (<ref>) by
    
        F^- L_d : (_n, _n+1)    ↦ (_n,   _n)   = ( _n  , - D_1 L_d (_n, _n+1) ) , 
        F^+ L_d : (_n, _n+1)    ↦ (_n+1, _n+1) = ( _n+1,   D_2 L_d (_n, _n+1) ) .
    
    With this, the discrete Euler–Lagrange equations (<ref>) can be written as
    
        F^+ L_d (_n-1, _n) - F^- L_d (_n, _n+1) = 0 ,
    
    which motivates the introduction of the position-momentum form F_L_d : M→M of the variational integrator (<ref>) by
    
        _n     =           -  D_1 L_d (_n, _n+1) , 
        _n+1   =  D_2 L_d (_n, _n+1) .
    
    Given (_n, _n), Equation (<ref>) can be solved for _n+1. This is generally a nonlinearly implicit equation that has to be solved by some iterative technique like Newton's method.
    Equation (<ref>) is an explicit function, so to obtain _n+1 we merely have to plug in _n and _n+1.
    The corresponding Hamiltonian evolution map is
    
        F_L_d  :  M→M  :   ( _n, _n ) ↦ ( _n+1, _n+1 ) .
    
    In terms of the discrete fibre derivatives it can be equivalently expressed as
    
        F_L_d   = ( F^- L_d) ∘ F_L_d∘( F^- L_d)^-1 , 
        F_L_d   = ( F^+ L_d) ∘ F_L_d∘( F^+ L_d)^-1 , 
        F_L_d   = ( F^+ L_d) ∘( F^- L_d)^-1 .
    
    In position-momentum form, the variational integrator can be initialized by prescribing an initial position _0 in conjunction with the corresponding momentum _0 = ϑ (_0). We thus have an exact initialization mechanism, as _0 constitutes a well-defined second set of initial conditions which can be exactly determined.
    Starting with an initial position _0 and an initial momentum _0, the repeated solution of (<ref>) gives the same discrete trajectory _d = {_n}_n=0^N as (<ref>).
    
    
    
    
     §.§ Discrete Symplectic Structure
    
    
    In the following we will show why variational integrators can be considered symplectic integrators and shed some light on the relation with symplectic integrators.
    As in the continuous case, we can obtain a discrete Lagrangian one-form by computing the variation of the action for varying endpoints,
    
        ddA_d [_d^] |_ = 0   = ∑_n=0^N-1[ D_1  L_d (_n, _n+1) ·_n + D_2  L_d (_n, _n+1) ·_n+1] 
           = ∑_n=1^N-1[ D_1  L_d (_n, _n+1) + D_2  L_d (_n-1, _n) ] ·_n
           
        + D_1  L_d (_0, _1) ·_0 + D_2  L_d (_N-1, _N) ·_N .
    
    The two latter terms originate from the variation at the boundaries. They form the discrete counterpart of the Lagrangian one-form.
    However, there are two boundary terms that define two distinct one-forms on M×M,
    
    
    
        [    Θ_L_d^- ( q_0  , q_1 ) · ( _0   , _1 ) ≡           -  D_1 L_d (q_0,   q_1) ·_0 ,;    Θ_L_d^+ ( q_N-1, q_N ) · ( _N-1 , _N )              ≡ D_2 L_d (q_N-1, q_N) ·_N . ]
    
    In general, these one-forms are defined as
    
        Θ_L_d^- ( _n , _n+1 )    ≡           -  D_1 L (_n, _n+1)  _n , 
        Θ_L_d^+ ( _n , _n+1 )    ≡ D_2 L (_n, _n+1)  _n+1 .
    
    As L_d = Θ_L_d^+ - Θ_L_d^- and ^2 L_d = 0 one observes that
    
        Θ_L_d^+ = Θ_L_d^-
    
    such that the exterior derivative of both discrete one-forms defines the same discrete Lagrangian two-form or discrete symplectic form
    
        ω_L_d   = Θ_L_d^+
         = Θ_L_d^-
         = ∂^2 L_d (_n, _n+1)∂_n  ∂_n+1 _n∧_n+1      (no summation over n) .
    
    Now consider the exterior derivative of the discrete action (<ref>).
    Upon insertion of the discrete Euler–Lagrange equations (<ref>), it becomes
    
        A_d   = D_1 L_d (_0, _1) ·_0 + D_2 L_d (_N-1, _N) ·_N
        = Θ_L_d^+ (_N-1, _N) - Θ_L_d^- (_0, _1) .
    
    On the right hand side we find the just defined Lagrangian one-forms (<ref>).
    Taking the exterior derivative of (<ref>) gives
    
        ω_L_d (_0, _1) = ω_L_d (_N-1, _N) ,
    
    where _N-1 and _N are connected with _0 and _1 through the discrete Euler–Lagrange equations (<ref>).
    Therefore, (<ref>) implies that the discrete symplectic structure ω_L_d is preserved while the system advances from t=0 to t=Nh according to the discrete equations of motion (<ref>).
    As the number of time steps N is arbitrary, the discrete symplectic form ω_L_d is preserved at all times of the simulation.
    Note that this does not automatically imply that the continuous symplectic structure ω_L is preserved by the variational integrator.
    However, as can be seen by comparing (<ref>) and (<ref>), the discrete one-forms (<ref>) correspond to the canonical one-form _i ^i under pullback by the discrete fibre derivatives (<ref>). Thus conservation of the discrete symplectic form ω_L_d by the discrete Euler–Lagrange equations (<ref>) on M×M is equivalent to conservation of the canonical symplectic form Ω by the position momentum form (<ref>) on M.
    
    
    
    
    
     §.§ Variational Runge–Kutta Methods
    
    
    The derivation of higher-order variational integrators in either standard or position-momentum form is rather cumbersome.
    A convenient framework for the derivation of integrators of arbitrary order is provided by variational Runge–Kutta methods, which can be seen as a generalization of the position-momentum form.
    These methods constitute a special family of symplectic-partitioned Runge–Kutta methods for Lagrangian systems, which are of the form
    
        _n,i   = ∂ L∂ (_n,i, _n,i) ,    _n,i   = ∂ L∂ (_n,i, _n,i) , 
        _n,i   = _n + h ∑_j=1^s a_ij _n,j ,    _n,i   = _n + h ∑_j=1^sa̅_ij _n,j , 
        _n+1   = _n + h ∑_i=1^s b_i _n,i ,    _n+1   = _n + h ∑_i=1^sb̅_i Ṗ_n,i ,
    
    with coefficients satisfying the symplecticity conditions,
    
        b_ia̅_ij + b̅_j a_ji   = b_ib̅_j      and   b̅_i   = b_i .
    
    Here, s denotes the number of internal stages, a_ij and a̅_ij are the coefficients of the Runge–Kutta method and b_i and b̅_i the corresponding weights.
    Note that while _n,i and _n,i represent velocities and forces at the internal stages, they strictly speaking do not correspond to the time derivatives of _n,i and _n,i, respectively.
    As (_n,i, _n,i) is nothing else than a point in M, the concept of a time derivative of these quantities does not make any sense.
    Instead, _n,i and _n,i as well as _n,i and _n,i, respectively, denote independent degrees of freedom, which are however related by (<ref>).
    
    <cit.> show that variational Runge–Kutta methods (<ref>) correspond to the position-momentum form (<ref>) of the discrete Lagrangian
    
        L_d (_n, _n+1) = h ∑_i=1^s b_i  L ( _n,i, _n,i) .
    
    They can also be obtained from a discrete action principle similar to the Hamilton–Pontryagin principle presented in Section (<ref>).
    For discretizations of Gauss–Legendre type, like the midpoint Lagrangian (<ref>), this is achieved by extremizing a discrete action of the following form <cit.> (see also <cit.>),
    
        𝒜_d =
        ∑_n=0^N-1
        h ∑_i=1^s b_i [ L ( _n,i, _n,i)
        + _n,i·( _n,i - _n - h ∑_j=1^s a_ij _n,j) ] 
        
        - _n+1·( _n+1 - _n - h ∑_i=1^s b_i _n,i)
         .
    
    Here, the definition of the generalized coordinates at the internal stages _n,i and the update rule determining _n+1 are added as constraints with the corresponding momenta _n,i and _n+1 taking the role of Lagrange multipliers.
    Requiring stationarity of the discrete action (<ref>) for arbitrary variations of _n, _n+1, _n,i, _n,i and _n,i, we recover (<ref>) with the conditions (<ref>) automatically satisfied.
    For discretizations of Lobatto–IIIA type like the trapezoidal Lagrangian (<ref>), where the first internal stage coincides with the solution at the previous time step, the velocities _n,i are not linearly independent and the discrete action (<ref>) needs to be augmented by an additional constraint to take this dependence into account (for details see <cit.>),
    
        𝒜_d =
        ∑_n=0^N-1
        h ∑_i=1^s b_i [ L ( _n,i, _n,i)
        + _n,i·( _n,i - _n - h ∑_j=1^s a_ij _n,j) ] 
        
        - _n+1·( _n+1 - _n - h ∑_i=1^s b_i _n,i)
        + μ_n·( ∑_i=1^s d_i_n,i)
         .
    
    Requiring stationarity of (<ref>), we obtain a modified system of equations,
    
        _n,i   = ∂ L∂ (_n,i, _n,i) ,    _n,i   = ∂ L∂ (_n,i, _n,i) , 
        _n,i   = q_n + h ∑_j=1^s a_ij _n,j ,    _n,i   = p_n + h ∑_j=1^sa̅_ij _n,j - μ_nd_ib_i , 
        _n+1   = q_n + h ∑_i=1^s b_i _n,i ,    _n+1   = p_n + h ∑_i=1^sb̅_i _n,i , 
        
        0    = ∑_i=1^s d_i_n,i ,
    
    accounting for the linear dependence of the Q̇_n,i and consequently also of the P_n,i.
    The particular values of d_i depend on the number of stages s and the definition of the _n,i <cit.>. For two stages, we have d_1 = - d_2, so that we can choose, for example, d_1 = 1 and d_2 = -1, and (<ref>) becomes equivalent to the variational integrator of the trapezoidal Lagrangian (<ref>). For three stages, we can choose d = (12, -1, 12), and for four stages we can use d = (+1, -√(5), +√(5), -1).
    
    
    
    
     §.§ Variational Integrators and Degenerate Lagrangians
    
    
    In this section, we want to give an overview of some approaches for discretizing degenerate Lagrangian systems which fail and discuss why they fail.
    
    The most obvious option for obtaining a geometric integrator for any Lagrangian systems is to directly discretize the Lagrangian (<ref>) and compute the corresponding discrete Euler–Lagrange equations (<ref>), followed by a discrete fibre derivative (<ref>) in order to obtain the position-momentum form (<ref>) of a variational integrator.
    Indeed, it has been shown by  <cit.> that this is a viable strategy for those cases where ϑ is a linear function.
    For cases where ϑ is a nonlinear function, however, we observe in simulations with such integrators that the numerical solution will in general not satisfy the constraint (<ref>).
    Thus the discrete trajectory { (q_n, p_n) }_n=0^N will drift away from the constraint submanifold (<ref>), i.e., even though (q_0, p_0) ∈ i(Δ) we usually find that (q_n, p_n) ∉ i(Δ) for n ≥ 1. Whence the solution becomes unphysical, c.f., Figure <ref>.
    In the standard form of the variational integrators, which are multi-step methods, this behaviour can be explained in terms of parasitic modes <cit.>.
    
    In some cases, the solution stays close to the constraint submanifold, which is to say that ϕ (q_n, p_n) although not zero at least stays bounded for very long times. In such cases variational integrators might still be a viable solution method. However, in general it is not clear to which extend this behaviour depends on the initial conditions. It is easily perceivable that the deviation from the constraint submanifold is bounded for some initial conditions but not for others. And indeed, we observe such behaviour for the example of guiding centre dynamics, described in Section <ref>, where certain particles are found to stay close to the constraint submanifold for very long times while other particles diverge further and further from the constraint submanifold as the simulation proceeds until it eventually crashes.
    
    
    
    
    Let us note that for variational Runge–Kutta methods (<ref>) and (<ref>) the constraint is automatically satisfied for the internal stages of the method, so that ϕ(Q_n,i, P_n,i) = 0 for all 1 ≤ i ≤ s, but not for the solution at the next time step (q_n+1, p_n+1). This is due to the fact that the internal stages obey discrete versions of the equations of motion, whereas the final step merely amounts to a numerical quadrature.
    It appears, though, that the drift-off problem could be avoided by choosing particular coefficient matrices a_ij and a̅_ij and weights b_i and b̅_i in (<ref>) or (<ref>) such that the last internal stage corresponds to the solution at the next time step, that is Q_n,s = q_n+1 and P_n,s = p_n+1.
    It turns out, however, that such a choice is incompatible with the symplecticity conditions (<ref>), i.e., variational Runge–Kutta methods for which the coefficients and weights are such that the Dirac constraint is automatically satisfied do not exist <cit.>.
    
    It seems then natural to augment the discrete action (<ref>) or (<ref>) with the constraint evaluated at the solution at the next time step via a Lagrange multiplier similar to (<ref>).
    This approach and why it fails will be discussed in some detail in the next section.
    
    Yet another seemingly natural approach is to apply a discrete version of the Hamilton–Pontryagin principle (<ref>) as proposed by <cit.>.
    We will see in Section <ref> that the resulting integrator is exactly equivalent to the position-momentum form (<ref>) and therefore shares the same problems.
    
    After discussing in more detail why these approaches fail, we will present several projection methods for enforcing the constraint ϕ (q_n+1, p_n+1) = 0 in Section <ref>.
    In these, we take the solution of a variational integrator and project it to the constraint submanifold (<ref>).
    Although these methods are not strictly-speaking variational or geometric, they lead to useful long-time stable integration algorithms.
    
    
    
    
     §.§ Augmented Variational Runge–Kutta Methods
    
    
    For degenerate Lagrangian systems, we see that the constraint p = ϑ(q) is automatically enforced at the internal stages. But as the constraint is not enforced at the time steps n, the solution tends to drift away from the constraint submanifold (<ref>).
    It seems to be natural to add the constraint to the discrete action (<ref>) or (<ref>), e.g.,
    
        𝒜_d =
        ∑_n=0^N-1
        h ∑_i=1^s b_i [ L ( Q_n,i, Q̇_n,i)
        - Ṗ_n,i·( Q_n,i - q_n - h ∑_j=1^s a_ij Q̇_n,j) ]
        
        
        + p_n+1·( q_n+1 - q_n - h ∑_i=1^s b_i Q̇_n,i) + ϕ (q_n+1, p_n+1) ·λ_n+1 .
    
    The momenta are denoted by p_n+1 instead of p_n+1 in order to highlight the problem with this approach as discussed below.
    Variation of (<ref>) leads to a modified integrator, which upon defining
    
        p_n   = p_n - ϕ_^T (q_n, p_n)  λ_n      and   q_n   = q_n - ϕ_^T (q_n, p_n)  λ_n ,
    
    can be written as
    
        P_n,i   = ∂ L∂ (Q_n,i, Q̇_n,i) ,    Ṗ_n,i   = ∂ L∂ (Q_n,i, Q̇_n,i) , 
        
        Q_n,i   = q_n + h ∑_j=1^s a_ij Q̇_n,j ,    
        P_n,i   = p_n + h ∑_j=1^sa̅_ij Ṗ_n,j , 
        q_n+1   = q_n + h ∑_i=1^s b_i Q̇_n,i ,    p_n+1   = p_n + h ∑_i=1^sb̅_i Ṗ_n,i ,
    
    with projection
    
        q_n+1   = q_n+1 + ϕ_^T (q_n+1, p_n+1)  λ_n+1 , 
        
        p_n+1   = p_n+1 - ϕ_^T (q_n+1, p_n+1)  λ_n+1 , 
        
        0    = ϕ (q_n+1, p_n+1) ,
    
    where we assume that p_0 = ϑ(q_0) so that p_0 = p_0 and λ_0 = 0.
    We observe that the constraint is enforced at (q_n+1, p_n+1), that is the projected coordinate but the unprojected momentum. Practically, we fix the momentum p_n+1 and change the coordinate _n+1 until it matches the constraint ϕ (q_n+1, p_n+1).  Then, the momentum is shifted using λ_n+1 as determined from the projection of _n+1. The result is that while (q_n+1, p_n+1) is guaranteed to lie on the constraint submanifold, this is not the case for (q_n+1, p_n+1).
    
    
    
    
     §.§ Discrete Hamilton–Pontryagin Principle
    
    
    As the Hamilton–Pontryagin principle of Section <ref> provides a very natural setting for degenerate Lagrangian systems, where the Dirac constraint appears as one of the equations of motion, it appears as an appropriate starting point for discretization.
    To that end, consider the (+) and (-) discrete Lagrange–Pontryagin principles proposed by <cit.> and given by
    
        δ∑_n=0^N-1[ L_d (q_n, q_n^+) + p_n+1· (q_n+1 - q_n^+) ] = 0 ,
    
    and
    
        δ∑_n=0^N-1[ L_d (q_n+1^-, q_n+1) - p_n· (q_n - q_n+1^-) ] = 0 ,
    
    respectively.
    Computing the variations and assuming that the variations of  are fixed at the endpoints, δ q_0 = δ q_N = 0, we obtain
    
        q_n+1   = q_n^+ ,    
        p_n   = - D_1 L_d (q_n, q_n^+) ,    
        p_n+1   = D_2 L_d (q_n, q_n^+) ,
    
    as well as
    
        q_n   = q_n+1^- ,    
        p_n   = - D_1 L_d (q_n+1^-, q_n+1) ,    
        p_n+1   = D_2 L_d (q_n+1^-, q_n+1) ,
    
    which in both cases are immediately recognized as being equivalent to the position-momentum form (<ref>) and therefore subject to the same instabilities.
    
    
    
    
    
    § PROJECTION METHODS
    
    
    Projection methods are a standard technique for the integration of ordinary differential equations on manifolds <cit.>.
    The problem of constructing numerical integrators on manifolds with complicated structure is often difficult and thus avoided by embedding the manifold into a larger space with simple, usually Euclidean structure, where standard integrators can be applied.
    Projection methods are then used to ensure that the solution stays on the correct subspace of the extended solution space, as that is usually not guaranteed by the numerical integrator itself.
    
    In the standard projection method, a projection is applied after each step of the numerical algorithm.
    Assuming that the initial condition lies in the manifold, the solution of the projected integrator will stay in the manifold.
    The problem with this approach is that even though assuming that the numerical integrator is symmetric, the whole algorithm comprised of the integrator and the projection will not be symmetric.
    This often leads to growing errors in the solution and consequently a drift in the total energy of the system.
    This can be remedied by symmetrizing the projection <cit.>, where the initial data is first perturbed out of the constraint submanifold, before the numerical integrator is applied, and then projected back to the manifold.
    This leads to very good long-time stability and improved energy behaviour.
    
    While such projection methods, both standard and symmetric ones, are standard procedures for conserving energy, as well as holonomic and non-holonomic constraints, not much is known about their application to Dirac constraints. 
    Some authors consider general differential algebraic systems of index two <cit.>, the class to which the systems considered here belong, but a discussion of symplecticity seems to be mostly lacking from the literature, aside from some remarks on the conservation of quadratic invariants by the post-projection method of <cit.>.
    
    In the following, we apply several projection methods (standard, symmetric, symplectic, midpoint) to variational integrators in position-momentum form.
    As it turns out, both the standard projection and the symmetric projection are not symplectic. The symmetric projection nevertheless shows very good long-time stability, as it can be shown to be pseudo-symplectic. The symplectic projection method, as the name suggests, is indeed symplectic, although in a generalized sense. The midpoint projection method is symplectic in the usual sense but only for particular integrators.
    
    The general procedure is as follows.
    We start with initial conditions _n on Δ (recall that for the particular Lagrangian (<ref>) considered here, the configuration manifold M and the constraint submanifold Δ are isomorphic, so that we can use the same coordinates on Δ as we use on M).
    We compute the corresponding momentum p_n by the continuous fibre derivative (<ref>), which yields initial conditions (q_n, p_n = ϑ(q_n)) on M satisfying the constraint ϕ(q_n, p_n) = 0. This corresponds to the inclusion map (<ref>).
    Then, we may or may not perturb these initial conditions off the constraint submanifold by applying a map (q_n, p_n) ↦ (q_n, p_n) which is either the inverse P^-1 of a projection P : M→ i(Δ) or, in the case of the standard projection of Section <ref>, just the identity.
    The perturbation is followed by the application of some canonically symplectic algorithm Ψ_h on M, namely a variational integrator in position-momentum form (<ref>) or a variational Runge–Kutta method (<ref>) or (<ref>), in which cases we have that Ψ_h = ( F^+ L_d) ∘( F^- L_d)^-1.
    In general, the result of this algorithm, (q_n+1, p_n+1) = Ψ_h (q_n, p_n), will not lie on the constraint submanifold (<ref>). Therefore we apply a projection (q_n+1, p_n+1) ↦ (q_n+1, p_n+1) which enforces ϕ (q_n+1, p_n+1) = p_n+1 - ϑ(q_n+1) = 0. As this final result is a point in i(Δ) it is completely characterized by the value _n+1.
    
    
        < g r a p h i c s >
    Gradient of the constraint function ϕ orthogonal and Ω-orthogonal to constant surfaces of ϕ(, ) =  - √(_0^2 - ^2) for _0∈{ 1, 2, 3 }.
    
    Let us emphasize that in contrast to standard projection methods, where the solution is projected orthogonal to the constrained submanifold, along the gradient of ϕ, here the projection has to be Ω-orthogonal (c.f., Figure <ref>), where Ω is the canonical symplectic matrix (<ref>). That is, denoting by λ the Lagrange multiplier, the projection step is given by Ω^-1∇ϕ^T instead of an orthogonal projection ∇ϕ^T. This appears quite natural when comparing with (<ref>).
    
    Let us also note that, practically speaking, the momenta _n and _n+1 are merely treated as intermediate variables much like the internal stages of a Runge–Kutta method.
    The Lagrange multiplier , on the other hand, is determined in different ways for the different methods and can be the same or different in the perturbation and the projection. It thus takes the role of an internal variable only for the standard, symmetric projection and midpoint projection, but not for the symplectic projection.
    
    
    
    
     §.§ Projected Fibre Derivatives
    
    
    In the following, we will try to underpin the construction of the various projection methods with some geometric ideas.
    We already mentioned several times that the position-momentum form of the variational integrator (<ref>) suffers from the problem that it does not preserve the constraint submanifold Δ defined in (<ref>).
    That is, even though it is applied to a point in i(Δ), it usually returns a point in M, but outside of i(Δ).
    In order to understand the reason for this, let us define Δ_M^- and Δ_M^+ as the subsets of M×M which are mapped into the constraint submanifold i(Δ) by the discrete fibre derivatives F^- L_d and F^+ L_d, respectively, i.e.,
    
        Δ_M^-   = { (_n, _n+1) ∈M×M | F^- L_d (_n, _n+1) = (_n, _n) ∈ i(Δ) } , 
        Δ_M^+   = { (_n, _n+1) ∈M×M | F^+ L_d (_n, _n+1) = (_n+1, _n+1) ∈ i(Δ) } ,
    
    or more explicitly,
    
        Δ_M^-   = { (_n, _n+1) ∈M×M |  - D_1 L_d (_n, _n+1) = ϑ(_n) } , 
        Δ_M^+   = { (_n, _n+1) ∈M×M |  D_2 L_d (_n, _n+1) = ϑ(_n+1) } .
    
    A sufficient condition for the position-momentum form of the variational integrator (<ref>) to preserve the constraint submanifold (<ref>) would be that Δ_M^- and Δ_M^+ are identical. 
    Slightly weaker necessary conditions can be formulated depending on the formulation of the position-momentum form in terms of the discrete Euler–Lagrange equations (<ref>) and the discrete fibre derivative (<ref>), c.f., Equation (<ref>).
    For example, considering (<ref>), a necessary condition for the position-momentum form to preserve Δ is that the image of the inverse of F^- L_d, namely Δ_M^-, is in Δ_M^+,
    
        ( F^- L_d)^-1 (i(Δ)) = Δ_M^-⊂Δ_M^+ .
    
    Further, from (<ref>) and (<ref>) it follows that the image of the variational integrator F_L_d applied to Δ_M^- must be in Δ_M^- and the image of F_L_d applied to Δ_M^+ must be in Δ_M^+,
    
        F_L_d( Δ_M^-)    ⊂Δ_M^- ,    
        F_L_d( Δ_M^+)    ⊂Δ_M^+ .
    
    None of these conditions can be guaranteed and they are in general not satisfied. Although Δ_M^- and Δ_M^+ might have some overlap, they are usually not identical, and the variational integrator, applied to a point in Δ_M^- or Δ_M^+, does not necessarily result in a point in Δ_M^- or Δ_M^+, respectively.
    
    In order to construct a modified algorithm which does preserve the constraint submanifold, we compose the discrete fibre derivatives F^± with appropriate projections P^±,
    
        (_n, _n  )
           = ( P^-∘F^- L_d) (_n, _n+1)
         = P__n^-^-( _n, -  D_1 L_d (_n, _n+1) ) , 
        
        (_n+1, _n+1)
           = ( P^+∘F^+ L_d) (_n, _n+1)
         = P__n+1^+^+( _n+1, D_2 L_d (_n, _n+1) ) ,
    
    so that they take any point in M×M to the constraint submanifold Δ.
    The Lagrange multiplier λ is indicated as subscript and implicitly determined by requiring that the constraint ϕ is satisfied by the projected values of  and .
    These projected fibre derivatives will not be a fibre-preserving map anymore, but they will change both  and , mimicking the continuous equations (<ref>).
    Noting that the nullspace of P_ is the span of Ω^-1∇ϕ, a natural candidate for the projection P_ is given by
    
        P_λ^± (, ) : (, )    = (, ) ± h  Ω^-1∇ϕ^T (, ) λ ,    
        0    = ϕ(, ) ,
    
    so that ( P^-∘F^- L_d ) (_n, _n+1) explicitly reads
    
        _n   = _n - h  ϕ_^T (_n, _n) _n^- , 
        _n   = - D_1 L_d (_n, _n+1) + h  ϕ_^T (_n, _n) _n^- , 
        
        0    = ϕ(_n, _n) ,
    
    and ( P^+∘F^+ L_d ) (_n, _n+1) explicitly reads
    
        _n+1   = _n+1 + h  ϕ_^T (_n+1, _n+1) _n+1^+ , 
        _n+1   = D_2 L_d (_n, _n+1) - h  ϕ_^T (_n+1, _n+1) _n+1^+ , 
        
        0    = ϕ(_n+1, _n+1) .
    
    The signs in front of the projections have been chosen in correspondence with the signs of the discrete forces in <cit.>.
    With these projections we obtain all of the algorithms introduced in the following sections, except for the midpoint projection, in a similar fashion to the definition of the position-momentum form of the variational integrator (<ref>), as a map Δ→Δ which can formally be written as
    
        Φ_h = ( π_Δ∘P^+∘F^+ L_d) ∘( π_Δ∘P^-∘F^- L_d)^-1 .
    
    In total, we obtain algorithms which map _n into _n+1 via the steps
    
        Δ
        i(Δ)
        MM×MM
        i(Δ)
        Δ ,
    
    where π_Δ^-1 is identical to the inclusion (<ref>).
    The difference of the various algorithms lies in the choice of _n^- and _n+1^+ as follows
    
    
    Projection    _n^-   _n+1^+
    
    Standard      0   _n+1
    
    Symplectic    _n   R (∞)  _n+1
    
    Symmetric     _n+1/2   R (∞)  _n+1/2
    
    Midpoint      _n+1/2   R (∞)  _n+1/2
    
    
    For the symmetric, symplectic and midpoint projections, it is important to adapt the sign in the projection according to the stability function R(∞) of the basic integrator (for details see e.g. <cit.>).
    For the methods we are interested in, namely Runge–Kutta methods, the stability function is given by R(z) = 1 + z b^T ( - zA)^-1 e with e = (1, 1, ..., 1)^T∈^s, and we have R(∞)=1 or, more specifically, for Gauss–Legendre methods R(∞) = (-1)^s and for partitioned Gauss–Lobatto IIIA–IIIB and IIIB–IIIA methods we have R(∞) = (-1)^s-1.
    
    Let us remark that for the standard projection, the basic integrator and the projection step can be applied independently.
    Similarly, for the symplectic projection, the three steps, namely perturbation, numerical integrator, and projection, decouple and can be solved consecutively, as we use different Lagrange multipliers _n in the perturbation and _n+1 in the projection.
    For the symmetric projection and the midpoint projection, however, this is not the case. 
    There, we used the same Lagrange multiplier _n+1/2 in both the perturbation and the projection, so that the whole system has to be solved at once, which is more costly.
    This also implies that for the projection methods where _n^- and _n+1^+ are the same (possibly up to a sign due to R(∞)), strictly speaking we cannot write the projected algorithm in terms of a composition of two steps as we did in (<ref>). Instead the whole algorithm has to be treated as one nonlinear map.
    The idea of the construction of the methods is still the same, though. Only the midpoint projection of Section <ref> needs special treatment. There, the operator P_ is defined in a slightly more complicated way than in (<ref>), using different arguments in the projection step, which does not quite fit the general framework outlined here.
    
    
    [][tbh]
    	
        < g r a p h i c s >
    Illustration of the standard projection method. The solution is projected to the constraint submanifold Δ after each step of the numerical integrator Ψ_h.
    
     §.§ Standard Projection
    
    
    The standard projection method <cit.> is the simplest projection method.
    Starting from _n, we use the continuous fibre derivative (<ref>) to compute _n = ϑ (_n). Then we apply some symplectic one-step method Ψ_h to _n = (_n, _n) to obtain an intermediate solution _n+1,
    
        _n+1   = Ψ_h (_n) ,
    
    which is projected onto the constraint submanifold (<ref>) by
    
        _n+1   = _n+1 + h  Ω^-1∇ϕ^T (_n+1) _n+1 ,
    
    enforcing the constraint
    
        0    = ϕ (_n+1) .
    
    This projection method, combined with the variational integrators (<ref>), is not symmetric, and therefore not reversible. Moreover, it exhibits a drift of the energy, as has been observed before, e.g., for holonomic constraints <cit.>.
    
    
    
    
      §.§.§ Symplecticity
    
    
    In order to verify the symplecticity condition, we write the projection (<ref>) in terms of (, ), that is
    
        _n+1^i   = _n+1^i - h  ϕ_^i^k (_n+1, _n+1)  _n+1^k , 
        _n+1^i   = _n+1^i + h  ϕ_^i^k (_n+1, _n+1)  _n+1^k ,
    
    and assume that Ψ_h is a symplectic integrator so that
    
        _n^i∧_n^i = _n+1^i∧_n+1^i .
    
    We start by taking the exterior derivative of _n+1 and _n+1,
    
        _n+1^i   = _n+1^i + ( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) , 
        _n+1^i   = _n+1^i - ( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) .
    
    Take the wedge product of the two equations,
    
        _n+1^i∧_n+1^i
         = _n+1^i   ∧_n+1^i
         - _n+1^i∧( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) 
           + ( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) ∧_n+1^i
           - ( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) ∧( h  ϕ_^i^l (_n+1, _n+1)  _n+1^l) .
    
    The second and third term on the right-hand side become
    
        ( h  _n+1^kϕ_^i^k (_n+1   , _n+1) ) ∧_n+1^i
         + ( h  _n+1^kϕ_^i^k (_n+1  , _n+1) ) ∧_n+1^i
           = ( h  _n+1^k[ ϕ_^i^k (_n+1, _n+1)  _n+1^i + ϕ_^i^k (_n+1, _n+1)  _n+1^i] ) 
           = ( h  ^k [ ϕ^k (_n+1, _n+1) ] )
         = 0 .
    
    The term in square brackets vanishes as ϕ (_n+1, _n+1) = 0 and therefore ϕ (_n+1, _n+1) = 0.
    Further we have
    
        ( h  _n+1^kϕ_^i^k (_n+1, _n+1) ) 
           = h  _n+1^kϕ_^i^j^k (_n+1, _n+1)  _n+1^j
           + h  _n+1^kϕ_^i^j^k (_n+1, _n+1)  _n+1^j
           + h  ϕ_^i^k (_n+1, _n+1)  _n+1^k
           = - h  _n+1^kϑ_k,ij (_n+1)  _n+1^j
           - h  ϑ_k,i (_n+1)  _n+1^k ,
    
    and
    
        ( h  _n+1^lϕ_^i^l (_n+1, _n+1) )
           = h  _n+1^lϕ_^i^j^l (_n+1, _n+1)  _n+1^j
           + h  _n+1^lϕ_^i^j^l (_n+1, _n+1)  _n+1^j
           + h  ϕ_^i^l (_n+1, _n+1)  _n+1^l
           = h  _n+1^i ,
    
    where the terms involving ϕ_^i^j^k or ϕ_^i^j^k vanish as ϕ(p,q) is separable and the terms involving ϕ_^i^j^k vanish as ϕ is linear in .
    The wedge product of the two expressions becomes
    
        - ( h  ϕ_^i^k (_n+1, _n+1)  _n+1^k) ∧( h  ϕ_^i^l (_n+1, _n+1)  _n+1^l) = 
        
        = ( h  _n+1^kϑ_k,ij (_n+1)  _n+1^j
           + h  ϑ_k,i (_n+1)  _n+1^k) ∧( h  _n+1^i) = 
        
        = h^2 ϑ_j,i (_n+1)  _n+1^j∧_n+1^i .
    
    The result can be anti-symmetrized so that by using (<ref>) as well as (<ref>), we obtain
    
        _n^i∧_n^i
         = _n+1^i∧_n+1^i
         + h^22Ω_ij (_n+1)  _n+1^i∧_n+1^j
        
         + h^2 _n+1^kϑ_k,ij (_n+1)  _n+1^i∧_n+1^j .
    
    Using that the constraint ϕ(, ) = p - ϑ(q) = 0 holds for both (_n, _n) and (_n+1, _n+1), this can be rewritten as
    
        12Ω_ij (_n)  _n^i∧_n^j
         = 12Ω_ij (_n+1)  _n+1^i∧_n+1^j
         - h^22Ω_ij (_n+1)  _n+1^i∧_n+1^j
        
         - h^2 _n+1^kϑ_k,ij (_n+1)  _n+1^i∧_n+1^j ,
    
    and we see that the noncanonical symplectic form (<ref>) is not preserved, but in each step accumulates an error h^2 _n+1^kϑ_k,ij (_n+1)  _n+1^i∧_n+1^j + 12 h^2 Ω_ij (_n+1)  _n+1^i∧_n+1^j.
    In numerical simulations, this error accumulation usually manifests itself in form of a drift of the solution and the energy.
    
    
    
    
     §.§ Symmetric Projection
    
    
    To overcome the shortcomings of the standard projection, we consider a symmetric projection of the variational Runge–Kutta integrators following <cit.>, c.f., Figure <ref> (see also <cit.>). Here, one starts again by computing the momentum _n as a function of the coordinates _n according to the continuous fibre derivative, which can be expressed with the constraint function as
    
        0    = ϕ (_n) .
    
    Then the initial value _n is first perturbed,
    
        _n   = _n + h  Ω^-1∇ϕ^T (_n)  _n+1/2 ,
    
    followed by the application of some one-step method Ψ_h,
    
        _n+1   = Ψ_h (_n) ,
    
    and a projection of the result onto the constraint submanifold,
    
        _n+1   = _n+1 + h   R(∞)  Ω^-1∇ϕ^T (_n+1) _n+1/2 ,
    
    which enforces the constraint
    
        0    = ϕ (_n+1) .
    
    Here, it is important to note that Lagrange multiplier _n+1/2 is the same in both the perturbation and the projection step, and to account for the stability function R(∞) of the basic integrator, as mentioned before.
    The algorithm composed of the symmetric projection and some symmetric variational integrator in position-momentum form, constitutes a symmetric map
    
        Φ_h : _n↦_n+1 ,
    
    where, from a practical point of view, _n, _n+1 and _n+1/2 are treated as intermediate variables.
    
    
    
    
      §.§.§ Symplecticity
    
    
    In the following, we assume that R(∞) = 1. Then, the considerations of symplecticity for the symmetric projection follow along the very same lines as for the standard projection. 
    In addition to the projection, we also have to consider the perturbation. Assuming the integrator Ψ_h is such that
    
        _n^i∧_n^i = _n+1^i∧_n+1^i ,
    
    we obtain
    
        _n^i∧_n^i
         - h^22Ω_ij (_n)  ^i_n+1/2∧^j_n+1/2
         - h^2_n+1/2^kϑ_k,ij (_n)  _n^i∧_n+1/2^j = 
        
        = _n+1^i∧_n+1^i
         - h^22Ω_ij (_n+1)  ^i_n+1/2∧^j_n+1/2
         - h^2_n+1/2^kϑ_k,ij (_n+1)  _n+1^i∧_n+1/2^j .
    
    The symmetrically projected integrator admits a certain symmetry in the error terms and can be shown to be pseudo-symplectic <cit.>.
    It is worth to go one step back, and reconsider the derivation that leads to (<ref>). By the same considerations as for the standard projection, we obtain
    
        _n^i∧_n^i
         - ( h  ϕ_^i^k (_n, _n)  _n+1/2^k) ∧( h  ϕ_^i^l (_n, _n)  _n+1/2^l) = 
        
         = _n+1^i∧_n+1^i
         - ( h  ϕ_^i^k (_n+1, _n+1)  _n+1/2^k) ∧( h  ϕ_^i^l (_n+1, _n+1)  _n+1/2^l) .
    
    We see, that in general the symmetric projection is not symplectic unless
    
        ( h  ϕ_^i^k (_n, _n)  _n+1/2^k)    = ( h  ϕ_^i^k (_n+1, _n+1)  _n+1/2^k)
        for all i,
    
    as well as
    
        ( h  ϕ_^i^k (_n, _n)  _n+1/2^k)    = ( h  ϕ_^i^k (_n+1, _n+1)  _n+1/2^k)
        for all i,
    
    that is the initial perturbation is exactly the same as the final projection.
    While the first condition (<ref>) is not obvious, the second condition is immediately seen to be satisfied for ϕ (, ) = p - ϑ (q), as ϕ_ =, so that (<ref>) reduces to _n+1/2 = _n+1/2. 
    If the first condition is not satisfied, though, the method is not symplectic. However, as the error terms to the symplecticity condition appear on both sides of (<ref>), the accumulated error is much smaller than with the standard projection.
    
    Again, using that the constraint ϕ(, ) = p - ϑ(q) = 0 holds for both (_n, _n) and (_n+1, _n+1), the symplecticity condition (<ref>) can be rewritten as
    
        12Ω_ij (_n)  ( _n^i∧_n^j
         - h^2 _n+1/2^i∧_n+1/2^j)
         - h^2_n+1/2^kϑ_k,ij (_n)  _n^i∧_n+1/2^j = 
        
         = 12Ω_ij (_n+1)  ( _n+1^i∧_n+1^j
         - h^2 _n+1/2^i∧_n+1/2^j)
         - h^2_n+1/2^kϑ_k,ij (_n+1)  _n+1^i∧_n+1/2^j .
    
    This formulation suggests the following construction.
    
    [][bth]
    	
        < g r a p h i c s >
    Illustration of the post projection method. Starting on the constraint submanifold Δ, the numerical integrator Ψ_h moves the solution away from Δ in the first step. After each step, the solution is projected back onto Δ, but the perturbation at the beginning of each consecutive step is exactly the inverse of the previous projection, so that, practically speaking, the solution is projected back onto Δ only for output purposes.
    
     §.§ Symplectic Projection
    
    
    If we modify the perturbation (<ref>) to use the Lagrange multiplier at the previous time step, _n, instead of _n+1, that is we replace (<ref>) by
    
        0    = ϕ (_n) , 
        _n   = _n + h  Ω^-1∇ϕ^T (_n)  _n , 
        _n+1   = Ψ_h (_n) , 
        _n+1   = _n+1 + h   R(∞)  Ω^-1∇ϕ^T (_n+1) _n+1 , 
        
        0    = ϕ (_n+1) ,
    
    the symplecticity condition (<ref>) is modified as follows,
    
        12Ω_ij (_n)  ( _n^i∧_n^j
         - h^2 _n^i∧_n^j)
         - h^2_n^kϑ_k,ij (_n)  _n^i∧_n^j = 
        
         = 12Ω_ij (_n+1)  ( _n+1^i∧_n+1^j
         - h^2 _n+1^i∧_n+1^j)
         - h^2_n+1/2^kϑ_k,ij (_n+1)  _n+1^i∧_n+1^j ,
    
    implying the conservation of a modified symplectic form ω_ defined on an extended phasespace M×^d with coordinates (, ) by
    
        ω_
         = 12Ω_ij ()  ^i∧^j
         - h^22Ω_ij ()  ^i∧^j
         - h^2^kϑ_k,ij ()  ^i∧^j ,
    
    with matrix representation
    
        Ω_ =
        [        Ω - h^2·ϑ_;   h^2·ϑ_   - h^2Ω;          ] .
    
    To this corresponds a modified one-form ϑ_, such that ω_ = ϑ_, given by
    
        ϑ_ = ( ϑ_i (q) - h  λ^kϑ_k,i (q) )   (  q^i - h  λ^i )
    
    As noted by <cit.>, the modified perturbation (<ref>)-(<ref>) can be viewed as a change of variables from (, ) on M×^d to (, ) on M, and the projection (<ref>)-(<ref>) as a change of variables back from (, ) to (, ).
    The symplectic form ω_ on M×^d thus corresponds to the pullback of the canonical symplectic form ω on M by this variable transformation.
    
    Let us note that the sign in in front of the projection in (<ref>), given by the stability function of the basic integrator, has very important implications on the nature of the algorithm.
    If it is the same as in (<ref>), the character of the method is very similar to the symmetric projection method described before.
    If the sign is the opposite of the one in (<ref>), like for Gauss–Legendre Runge–Kutta methods with an odd number of stages, the perturbation reverses the projection of the previous step, so that we effectively apply the post-projection method of <cit.>. That is, the projected integrator Φ_h is conjugate to the unprojected integrator Ψ_h by
    
        Φ_h = P^-1∘Ψ_h∘P ,
    
    so that the following diagram commutes
    
     (m) [matrix of math nodes,row sep=3em,column sep=4em,minimum width=2em] _n   _n+1
    
    (_n, _n)     (_n+1, _n+1) 
    ;
    
    [-stealth, line width=.4mm]
    (m-1-1) edge node [above] Ψ_h       (m-1-2)
    (m-2-1) edge node [left]  P^-1   (m-1-1)
    (m-1-2) edge node [right] P        (m-2-2)
    (m-2-1) edge node [below] Φ_h       (m-2-2);
    
    
    and the projection is effectively only applied for the output of the solution, but the actual advancement of the solution in time happens outside of the constraint submanifold (c.f., Figure <ref>).
    In other words, applying n times the algorithm Φ_h to a point (_0, 0) is equivalent to applying the perturbation P^-1, then applying n times the algorithm Ψ_h and projecting the result with P.
    
    Potentially, this might degrade the performance of the algorithm. If the accumulated global error drives the solution too far away from the constraint submanifold, the projection step might not have a solution anymore. Interestingly, however, post-projected Gauss–Legendre Runge–Kutta methods retain their optimal order of 2s <cit.>. Moreover, for methods with an odd number of stages, the global error of the unprojected solution is O(h^s+1), compared to O(h^s) for methods with an even number of stages. In practice this seems to be at least part of the reason of the good long-time stability of these methods.
    
    
    
    
      §.§.§ Symplecticity
    
    
    While in the continuous case, the symplectic form on M is always degenerate, thus not symplectic but presymplectic, in the discretization this is changed. The discrete Lagrangian on M×M is in general not degenerate, thus the symplectic form on M×M is non-degenerate as well. Composing the usual position-momentum form with the projection to Δ⊂M, thus enforcing ϕ(, )=0 in the way outlined before, we effectively obtain an algorithm mapping M×^d into M×^d instead of the original variational integrator, which mapped M×M into M×M. However, the new algorithm preserves a true symplectic form on M×^d, which is not the same as the presymplectic form of the continuous dynamics, and also not the same as the discrete symplectic form on M×M.
    This change of the presymplectic form to a symplectic form appears to be due to the initial “non-conservation” of degeneracy when discretizing the Lagrangian in conjunction with the projection.
    
    
    
    
     §.§ Midpoint Projection
    
    
    For certain variational Runge–Kutta methods, we can also modify the symmetric projection in a different way in order to obtain a symplectic projection, namely by evaluating the projection at the midpoint
    
        _n+1/2   = (_n+1/2, _n+1/2) ,    _n+1/2   = 12( _n + _n+1) ,    _n+1/2   = 12( _n + _n+1) ,
    
    so that the projection algorithm becomes
    
        0    = ϕ (_n) , 
        _n   = _n + h  Ω^-1∇ϕ^T (_n+1/2)  _n+1/2 , 
        _n+1   = Ψ_h (_n) , 
        _n+1   = _n+1 + h  Ω^-1∇ϕ^T (_n+1/2) _n+1/2 , 
        
        0    = ϕ (_n+1) .
    
    This method can be shown to be symplectic with respect to the original noncanonical symplectic form on M if the integrator Ψ_h is a symmetric, symplectic Runge–Kutta method with an odd number of stages s, for which the central stage with index (s+1)/2 corresponds to _n+1/2.
    This is obviously the case for the implicit midpoint rule, that is the Gauss–Legendre Runge–Kutta method with s=1, but unfortunately not for higher-order Gauss–Legendre or for Gauss-Lobatto methods.
    However, following <cit.> and <cit.>, higher-order methods similar to Gauss–Legendre methods but satisfying the requested property can be obtained. See for example the method with three stages given in Table <ref>.
    
    
    
      §.§.§ Symplecticity
    
    
    In order to show symplecticity, we follow a similar path as before for the standard projection method.
    We start by computing the exterior derivative of the perturbation and projection steps,
    
        _n^i   = _n^i + ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) , 
        _n^i   = _n^i - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ,
    
    and
    
        _n+1^i   = _n+1^i - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) , 
        _n+1^i   = _n+1^i + ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) .
    
    Then we compute the wedge products _n^i∧_n^i,
    
        _n^i∧_n^i
         = _n^i∧_n^i   + ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧_n^i
           - _n^i∧( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) 
           - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ,
    
    and _n+1^i∧_n+1^i,
    
        _n+1^i∧_n+1^i
         = _n+1^i∧_n+1^i   - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧_n+1^i
           + _n+1^i∧( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) 
           - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) .
    
    Now assume that the integrator Ψ_h is symplectic and thus satisfies _n+1^i∧_n+1^i = _n^i∧_n^i, which allows us to insert the second equation into the first to obtain
    
        _n+1^i∧_n+1^i
         = _n^i∧_n^i   - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧ ( _n^i + _n+1^i ) 
           - ( h  ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^k) ∧ ( _n^i + _n+1^i ) .
    
    Noting that _n^i + _n+1^i = 2 _n+1/2 and _n^i + _n+1^i = 2 _n+1/2, we can rewrite the previous expression as
    
        _n+1^i∧_n+1^i
         = _n^i∧_n^i   - ( h  _n+1/2^k ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^i) 
           - ( h  _n+1/2^k ϕ_^i^k (_n+1/2, _n+1/2)  _n+1/2^i) .
    
    The last two terms can be combined, so that the symplecticity condition reads
    
        _n+1^i∧_n+1^i = _n^i∧_n^i
        - ( h  _n+1/2^k ϕ^k (_n+1/2, _n+1/2) ) .
    
    The additional terms vanish under the assumption that _n+1/2 = (_n+1/2, _n+1/2) is equivalent to one of the internal stages of the variational Runge–Kutta method.
    We pointed out before that for the internal stages the Dirac constraint ϕ(, )=0 is automatically satisfied by the first equation in (<ref>).
    Therefore, if _n+1/2 corresponds to one of the internal stages, we have that ϕ(_n+1/2, _n+1/2) = 0 and thus also ϕ (_n+1/2, _n+1/2) = 0 so that
    
        _n+1^i∧_n+1^i = _n^i∧_n^i .
    
    It is worth pointing out that this holds for arbitrary constraints ϕ(, ) = 0 and that we did not use the particular structure of (<ref>) like separability or ϕ_ =.
    Therefore, the midpoint projection method is applicable to arbitrary Hamiltonian systems with Dirac constraints, not just the degenerate Lagrangian systems discussed in this paper.
    
    
    
    
    
    § NUMERICAL EXPERIMENTS
    
    
    The projection methods described in the previous section have all been implemented in the  package, which is a library of geometric integrators for ordinary differential equations and differential algebraic equations in the Julia programming language <cit.> freely available on GitHub <cit.>.
    We use Newton's method with quadratic line search for solving the nonlinear systems and LU decomposition for solving the linear systems. The Jacobian is computed via automatic differentiation via the  package <cit.> and updated in every time step but only every five nonlinear iterations.
    If possible, the numerical integration step and the projection step are solved separately (that is for the standard and symplectic projection, but not for the symmetric and midpoint projection).
    The updates of the solution are computed using compensated summation (Kahan's algorithm) in order to reduce the propagation of round-off errors.
    
    The examples we will consider are a two-dimensional Lotka–Volterra model, planar point vortices with varying circulation and guiding centre dynamics.
    The first two examples are implemented in the  package. The latter is implemented in the  package. Both packages are also available on GitHub <cit.>.
    Except for the first example, all systems possess Noether symmetries and some related conservation law, whose preservation will be monitored in the simulations.
    
    We perform simulation with Gauss–Legendre Runge–Kutta methods with one to six stages as well as Gauss–Lobatto–IIIA, IIIB, IIIC, IIID and IIIE methods <cit.> with two, three and four stages. 
    Here, the referenced method always provides the coefficients a and the coefficients a̅ are chosen, such that the symplecticity condition (<ref>) is satisfied. That is Gauss–Lobatto–IIIA denotes the IIIA–IIIB pair, Gauss–Lobatto–IIIB denotes the IIIB–IIIA pair, and Gauss–Lobatto–IIIC denotes the IIIC–IIIC* pair. For the Gauss–Legendre as well as Gauss–Lobatto–IIID and IIIE methods, we have a̅ = a.
    The Gauss–Lobatto–IIIC* method is sometimes also referred to as Gauss–Lobatto–III. Similar inconsistent naming is found for the IIID and IIIE methods. Here, we denote by the Gauss–Lobatto–IIID method the special case of the Gauss–Lobatto–IIIS method with σ=1.0 and by the IIIE method the special case of the IIIS method with σ=0.5. The Gauss–Lobatto–IIIS methods are interpolations of the IIIA, IIIB, IIIC and IIIC* methods with coefficients given by
    
        a_ij^S (σ) = (1-σ) ( a_ij^A + a_ij^B) + (σ - 12)  ( a_ij^C + a_ij^C*) ,
    
    so that
    
        a_ij^D = a_ij^S (1) = 12 (a_ij^A + a_ij^B ) ,
    
    and
    
        a_ij^E = a_ij^S (12) = 12 (a_ij^C + a_ij^C* ) .
    
    We compare the results of the variational Runge–Kutta methods with simulations of Radau–IIA methods, which have the advantage that they automatically preserve the Dirac constraint but also have the disadvantage of dissipating energy.
    
    For all methods, we perform simulations both without projection and with standard, symmetric, symplectic and midpoint projection. Due to the limited space we will only show some selected examples. The collection of all simulation results can be found in the documentation of the  package <cit.>.
    
    For most examples the simulations with the Gauss–Lobatto–IIIA, IIIB and IIIC methods break down after very few time steps. Even when reducing the time step by an order of magnitude, the IIIA, IIIB and IIIC methods perform rather poorly in almost all of the experiments. For the unprojected integrator this was already shown as a motivating example in Figure <ref> for the Lotka–Volterra model.
    The origin of this behaviour is most likely related to the fact that for the IIIA, IIIB and IIIC methods, different Runge–Kutta coefficient a_ij and a_ij are used for the integration of the trajectory  and the conjugate momenta p. Even though the nodes of the stages c_i are the same for both  and p, the definition of the values at the nodes _n,i and _n,i in terms of the corresponding vector fields _n,i and _n,i is different. While this is usually fine for regular problems, especially with separable Hamiltonians, it does not seem appropriate for degenerate problems where there is a functional relationship between the momenta and the position along the trajectory given at the internal stages by _n,i = ϑ (_n,i) for 1 ≤ i ≤ s. This particular property of degenerate systems suggests that the same coefficient matrices should be used for the definition of the internal stages of both  and p.
    
    For the two- and four-stage Gauss–Lobatto–IIIA, IIIB and IIIC methods, the symplectic projection amounts to a post-projection. Therefore, if the simulation without projection breaks down, so does the simulation with symplectic projection.
    The midpoint projection is only symplectic for the Gauss–Legendre method with one stage and the SRK3 method whose tableau was given in Table <ref>. Nevertheless, we run experiments with this projection and all integrators to study the long time behaviour.
    
    
    
    
     §.§ Lotka–Volterra Model
    
    
    Lotka–Volterra models <cit.> are used in mathematical biology for modelling population dynamics of animal species, as well as many other fields where predator-prey and similar models appear.
    The dynamics of the growth of two interacting species can be modelled by the following Lagrangian system <cit.>,
    
        L (, ) = ( log_2_1 + _2) _1 + _1_2 - H () ,
    
    with the Hamiltonian H given by
    
        H () = a_1_1 + a_2_2 - b_1log_1 - b_2log_2 .
    
    The noncanonical symplectic form (<ref>) is computed as
    
        ω = - 1_1_2 _1∧_2 .
    
    In the position-momentum form, which is the basis for the variational Runge-Kutta methods we employ in the numerical experiments, we obtain the following functions for the momenta and forces,
    
        ϑ_1 ()    = log_2_1 + _2 ,    
        f_1 (, )    = _2 - log_2_1^2 _1 - a_1 + b_1_1 , 
        ϑ_2 ()    = _1 ,    
        f_2 (, )    = ( 1 + 1_1_2) _1 - a_2 + b_2_2 .
    
    In the simulations, we use a time step of h = 0.1 and consider initial conditions (q_1,0, q_2,0) = (1, 1) with parameters (a_1, a_2, b_1, b_2) = (1, 1, 1, 2), which give a periodic solution.
    
    
    
    
    
    We make the following observations:
    
      *  The Gauss–Legendre Runge–Kutta methods with an odd number of stages  (Figure <ref>) as well as the SRK3 method (Figure <ref>) are stable even without projection. Even though they do not preserve the Dirac constraint exactly, the error in the constraint oscillates about zero and the amplitude of that oscillation appears to be bounded or at least grows only slowly. A similar behaviour is observed for the Gauss–Lobatto–IIID and IIIE methods (not shown).
    
    
      *  The Gauss–Legendre Runge–Kutta methods with an even number of stages (Figure <ref>) show an increasing error in the Dirac constraint and also in the energy, which eventually renders the simulation unstable (after about 250 000 time steps for the two-stage method and after about 1 000 000 time steps for the four-stage method).
    
    
      *  The Gauss–Lobatto–IIIA, IIIB and IIIC methods (not shown) are unstable without projection. For the integrator with two stages, the simulation crashes after about 25 time steps. For the integrator with three stages, it crashes immediately on the first time step. Decreasing the time step to h = 0.01 both integrator run for a short period. The integrator with two stages crashes after about 350 time steps and the integrator with three stages after about 1.000 time steps.
    
    
      *  The standard projection leads to very good results with all Gauss–Legendre methods (Figures <ref>, <ref>), the SRK3 method (Figure <ref>), as well as the Gauss–Lobatto–IIID and IIIE methods (not shown), but not with the Gauss–Lobatto–IIIA, IIIB and IIIC methods (not shown), whose solution deteriorates quickly. We observe small drifts in the energy error, but over 10 000 000 time steps this drift is of the order of 10^-12.
    
    
      *  For the Gauss–Legendre methods (Figures <ref>, <ref>), the SRK3 method (Figure <ref>), and the Gauss–Lobatto–IIID and IIIE methods (not shown), the symmetric projection leads to similar results as the standard projection. In some cases the drift in the energy seems to be slightly larger than with the standard projection. This, however, is due to round-off errors (c.f., Section <ref>). The errors of the Gauss–Lobatto–IIIA, IIIB and IIIC methods (not shown) are smaller than with the standard projection, but there still is a substantial drift in the energy.
    
    
      *  The symplectic projection (Figures <ref>, <ref>, <ref>) leads to very good results with all Gauss–Legendre methods, the SRK3 method, as well as the Gauss–Lobatto–IIID and IIIE methods, comparable to the results obtained with the symmetric projection. For the Gauss–Lobatto–IIIA, IIIB and IIIC methods with an even number of stages, the symplectic projection corresponds to a post-projection method. For this reason and as these methods are unstable without projection, the symplectic projection is also unstable. Although for an odd number of stages, the projection does not correspond to a post-projection, simulations still tend to crash as quickly as without projection.
    
    
      *  The midpoint projection (Figures <ref>, <ref>, <ref>) leads to good results with all Gauss–Legendre methods, the SRK3 method as well as the Gauss–Lobatto–IIID and IIIE methods, again comparable to the results obtained with the symmetric projection. However, it is only symplectic for the Gauss–Legendre method with one stage and the SRK3 method.
    
    
      *  With the Radau methods, we observe exact conservation of the Dirac constraint, as expected, but dissipation of energy, which is related to large errors in the solution.
    
    
    In summary, the numerical experiments for the Lotka–Volterra problem suggest that the Gauss–Legendre, the SRK and the Gauss–Lobatto–IIID and IIIE methods lead to good results with all projection methods, whereas (for the time step used) the the results with the Gauss–Lobatto–IIIA, IIIB and IIIC methods are never satisfactory, even with projection. For these methods, the time step needs to be reduced by at least a factor 10 in order to obtain stable simulations. With such small time steps, however, the Gauss–Lobatto–IIIA, IIIB and IIIC methods are not competitive anymore and one should rather use the Gauss–Legendre or Gauss–Lobatto–IIID and IIIE methods.
    For the Gauss–Legendre methods with an odd number of stages, the simulation appears to be stable even without projection, at least for very long times (10 million time steps), although the order of the integrators is decreased in this case (c.f., Section <ref>).
    It was already reported by <cit.> that for index two differential-algebraic equations Gauss–Legendre Runge–Kutta methods with an odd number of stages behave much better than those with an even number of stages, which is related to the stability function R(∞) being +1 for the former and -1 for the latter.
    
    
    
    
     §.§ Planar Point Vortices with Varying Circulation
    
    
    Systems of planar point vortices <cit.> provide a challenging problem for numerical integrators. Such systems are integrable for up to three vortices but produce chaotic behaviour for a minimum number of four vortices.
    An interesting phenomenon is that of leapfrogging, which is usually observed only for two pairs of point vortices.
    However, also one pair of point vortices can leapfrog by itself (see Figure <ref>) if the circulation is position dependent <cit.>. In this case, the function ϑ in the Lagrange is nonlinear, hence this provides an interesting test case for our integrators.
    
    We denote coordinates on M by (,) = (^1, , ^d, ^1, , ^d) and correspondingly coordinates on M by (,,,) = (^1, , ^d, ^1, , ^d, ^1, , ^d, ^1, , ^d). The coordinates on M are sometimes also collectively referred to by q.
    The general form of the Lagrangian for point vortices is 
    
        L
        = 12∑_i,j=1^dΓ_ij (^i^j - ^i^j )
        - 14π∑_i ≠ k^d∑_j ≠ l^dΓ_ijΓ_kl log( (^i - ^k)^2 + (^j - ^l)^2) ,
    
    with d the number of vortices and Γ the matrix of vortex strengths, which is assumed to be of the form Γ_ij = γ_iδ_ij, where γ_i is the circulation of the ith vortex.
    Here, we consider the special case of Γ being position-dependent, specifically Γ_ij (^i, ^i, ^j, ^j) = γ_i  S(^i, ^i)  δ_ij, where we assume later on that S is such that it has rotational symmetry.
    For d=2 the Lagrangian thus becomes
    
        L    = 12( γ^1  S(^1, ^1)   (^1^1 - ^1^1 ) + γ^2  S(^2, ^2)   (^2^2 - ^2^2 ) ) 
        - H (, ) , 
        
        H    = 12πγ^1γ^2  S(^1, ^1)   S(^2, ^2)  log( (^1 - ^2)^2 + (^1 - ^2)^2) .
    
    The noncanonical symplectic form (<ref>) of this system reads
    
        ω
         = ∑_i=1^2γ_i  S(^i, ^i)  ^i∧^i
         + 12∑_i=1^2γ_i( ^i ∂ S∂ (^i, ^i) + ^i ∂ S∂ (^i, ^i) ) ^i∧^i .
    
    Assuming that the function S is of the form S(a,b) = s(a^2 + b^2) with some function s : →, the Lagrangian is invariant under rotations of all coordinates by a constant angle ϑ, that is, the following transformation of the coordinates,
    
        σ^ : 
        [ ^i; ^i;    ]↦[ ^icos (ϑ) - ^isin (ϑ); ^icos (ϑ) + ^isin (ϑ);                       ] ,
    
    together with the corresponding transformation of the velocities, leaves the Lagrangian L invariant.
    The generating vector field is computed as
    
        V = d σ^d|_=0
          = - ^i∂∂^i + ^i∂∂^i ,
    
    and the corresponding conserved quantity (<ref>) is obtained as
    
        P = ∂ L∂· V
           = 12∑_i=1^2γ_i ( ^i^2 + ^i^2)   S(^i, ^i) 
           = q^1ϑ_2 (, ) - q^2ϑ_1 (, )
          + q^3ϑ_4 (, ) - q^4ϑ_3 (, )  .
    
    We are particularly interested in the behaviour of this angular momentum under the various projection methods.
    
    We consider the simple case of s(r) = 1 + r, so that S(a,b) = 1 + a^2 + b^2 and the functions for the momenta are computed as
    
        ϑ_1 (,)    = - 12γ^1 ^1  S(^1, ^1) ,    ϑ_2 (,)    = 12γ^1 ^1  S(^1, ^1) , 
        ϑ_3 (,)    = - 12γ^2 ^2  S(^2, ^2) ,    ϑ_4 (,)    = 12γ^2 ^2  S(^2, ^2) ,
    
    and those for the forces as
    
        f^1 (,,,)
           = 12γ^1  (^1^1 - ^1^1 )   S^() (^1, ^1)
         + 12γ^1^1  S(^1, ^1)
         - ∇_1 H (,) , 
        
        f^2 (,,,)
           = 12γ^2  (^2^2 - ^2^2 )   S^() (^2, ^2)
         + 12γ^2^2  S(^2, ^2)
         - ∇_2 H (,) , 
        
        f_3 (,,,)
           = 12γ^1  (^1^1 - ^1^1 )   S^() (^1, ^1)
         - 12γ^1^1  S(^1, ^1)
         - ∇_3 H (,) , 
        
        f_4 (,,,)
           = 12γ^2  (^2^2 - ^2^2 )   S^() (^2, ^2)
         - 12γ^2^2  S(^2, ^2)
         - ∇_4 H (,) ,
    
    with the gradient of the Hamiltonian being
    
        ∇_1 H (,)
           = γ^1γ^22π  S^() (^1, ^1)   S(^2, ^2)  log( (^1 - ^2)^2 + (^1 - ^2)^2) 
           + γ^1γ^2π  S(^1, ^1)   S(^2, ^2)  ^1 - ^2(^1 - ^2)^2 + (^1 - ^2)^2 , 
        ∇_2 H (,)
           = γ^1γ^22π  S^() (^2, ^2)   S(^1, ^1)  log( (^1 - ^2)^2 + (^1 - ^2)^2) 
           - γ^1γ^2π  S(^1, ^1)   S(^2, ^2)  ^1 - ^2(^1 - ^2)^2 + (^1 - ^2)^2 , 
        ∇_3 H (,)
           = γ^1γ^22π  S^() (^1, ^1)   S(^2, ^2)  log( (^1 - ^2)^2 + (^1 - ^2)^2) 
           + γ^1γ^2π  S(^1, ^1)   S(^2, ^2)  ^1 - ^2(^1 - ^2)^2 + (^1 - ^2)^2 , 
        ∇_4 H (,)
           = γ^1γ^22π  S^() (^2, ^2)   S(^1, ^1)  log( (^1 - ^2)^2 + (^1 - ^2)^2) 
           - γ^1γ^2π  S(^1, ^1)   S(^2, ^2)  ^1 - ^2(^1 - ^2)^2 + (^1 - ^2)^2 ,
    
    where S^() and S^() denote the  and  derivative of S, respectively.
    
    We use the time step h = 0.1, circulations γ_1 = γ_2 = 0.1 and initial conditions _0 = (1.0, 0.1, 1.0, -0.1). This setup leads to a circular leapfrogging of the two point vortices. 
    
    
    We make the following observations:
    
      *  All methods except the two-stage Gauss–Lobatto–IIIA and IIIB method and all of the Gauss–Lobatto–IIIC methods are stable even without projection, although with reduced order (c.f., Section <ref>).
    
    
      *  For all methods except the two-stage Gauss–Lobatto–IIIA and IIIB and the Gauss–Lobatto–IIIC methods method we observe that the angular momentum oscillates about its initial value where the amplitude of the oscillation seems bounded, that is the angular momentum seems to be preserved in a nearby sense, similar to the energy with symplectic integrators.
    
    
      *  The standard projection worsens the result for all methods except the two-stage Gauss–Lobatto–IIIA and IIIB and the Gauss–Lobatto–IIIC methods method, which do not crash when applying the projection, but the projected methods still show large errors and do not provide satisfactory results.
    
    
      *  The symplectic, symmetric and midpoint projections lead to very good results with almost all methods, restoring the original order of the methods and showing good long-time behaviour of both the energy and the angular momentum. There are some exceptions, however:
    
    
      *  The symplectic projection applied to the two-stage Gauss–Lobatto–IIIA and IIIB and the Gauss–Lobatto–IIIC methods is just as unstable as the corresponding unprojected methods.
    
      *  Both, the symmetric and midpoint projection applied to all of the Gauss–Lobatto–IIIC methods lead to an improved behaviour compared to the unprojected case, but exhibit a strong drift in the energy.
    
      *  The Gauss–Lobatto–IIID and IIIE methods with an even number of stages together with the midpoint projection exhibit a rather erratic behaviour in the energy error.
    
      *  For the symmetric projection and the higher-order methods (e.g. Gauss–Legendre with four or more stages or Gauss–Lobatto-IIID with four stages), we observe a small drift in the angular momentum, but over 1 000 000 time steps this drift is of the order of 10^-12. This drift is most likely caused by round-off errors (see Section <ref> for more details).
    
    
    In summary, the numerical experiments suggest that the combination of almost all integration methods and all projection methods excluding the standard projection provide suitable integration algorithms for the point vortex example. Exceptions are the Gauss–Lobatto–IIIC methods with any of the projection methods and the combination of the midpoint projection with the two-stage Gauss–Lobatto–IIIA and IIIB method and Gauss–Lobatto–IIID and IIIE methods with an even number of stages.
    
    
    
    
    
     §.§ Guiding Centre Dynamics
    
    
    In plasma physics, the search for geometric integrators for guiding centre dynamics and gyrokinetics is currently of great interest.
    As the Hamiltonian structure of the guiding centre system is noncanonical, there are practically no standard methods which can be easily applied.
    As the guiding centre equations can also be obtained from a Lagrangian, the application of variational integrators seems natural and has recently been tried by various researchers <cit.>.
    However, the guiding centre Lagrangian is degenerate, leading to all the problems discussed so far. We will see in the following if our projection methods can overcome these deficits.
    
    Guiding centre dynamics <cit.> is a reduced version of charged particle dynamics, where the motion of the particle in a strong magnetic field B is reduced to the motion of the guiding centre, that is the centre of the gyro motion of the particle about a magnetic field line.
    The dynamics of the guiding centre can be described in terms of only four coordinates (as compared to six for the full motion of the charged particle), the position of the guiding centre  and the parallel velocity , where parallel refers to the direction of the magnetic field.
    Denoting coordinates on M by (, ) = (^1, ^2, ^3, ) and correspondingly coordinates on M by (, , , ) = (^1, ^2, ^3, , ^1, ^2, ^3, ), the guiding centre Lagrangian <cit.> can be written as
    
        L    = (A () +   b ()) · - H(,) ,    
        H    = 12^2 + μB () ,
    
    where b = B / B is the unit vector of the magnetic field B = ∇× A with A the magnetic vector potential and μ is the magnetic moment. The first term in H denotes the parallel part of the kinetic energy and the second term the perpendicular part (parallel and perpendicular to the direction of the magnetic field). Here, we consider the case of only a magnetic field with vanishing electrostatic potential.
    
    Denoting a curve in M by t ↦ (x(t), u(t)), the Euler–Lagrange equations (<ref>) are computed as follows,
    
        ∇ϑ^T( x(t), u(t) ) ·ẋ (t) - ϑ̇( x(t), u(t) )    = ∇ H ( x(t), u(t) ) , 
        
        b ( x(t) ) ·ẋ (t)    = u (t) ,
    
    with ϑ(,) = A () +   b () and the gradient denoting the derivative with respect to .
    This can be rewritten in an explicit form as
    
        ẋ (t)    = u (t)  β (x (t))b (x (t)) ·β (x (t)) + B (x (t))B (x (t)) ·β (x (t))×∇ H (x (t),u (t)) , 
        u̇ (t)    = - β (x (t))b (x (t)) ·β (x (t))·∇ H (x (t),u (t)) ,
    
    where β = ∇×ϑ.
    The noncanonical symplectic form (<ref>) is given by
    
        ω = 12 ( ϑ_j,i (, ) - ϑ_i,j (, ) )  ^i∧^j - b_i ()  ^i∧ .
    
    Let us assume that the magnetic field B is not uniform, but that both A and B do not depend on one of the coordinates, say _3.
    Than we have a symmetry for the transformation
    
        σ^ : ^3↦^3 +  ,
    
    with generating vector field
    
        V = d σ^d|_=0
          = ∂∂_3 .
    
    The corresponding conserved momentum map (<ref>),
    
        P = ∂ L∂· V
          = ϑ_3 (,) ,
    
    which, depending on the actual form of ϑ, can be quite complicated and is therefore a good test for our algorithms.
    Although the basic integrator will preserve this toroidal momentum if the discrete Lagrangian preserves the corresponding symmetry, the projection could potentially modify its value.
    The projection guarantees preservation of the constraint p_n+1 = ϑ(q_n+1) but it does not guarantee that p_n+1 = p_n.
    
    In the numerical experiments, we use toroidal coordinates = (R, Z, ), where R, Z and  denote the radial, vertical and toroidal direction, respectively.
    For the magnetic field B and the vector potential A we will use analytic expressions following <cit.>.
    The vector potential is given as
    
        A_R   = B_0 R_0 Z2 R ,    
        A_Z   = - ln( RR_0)  B_0 R_02 ,    
        A_   = - B_0 r^22 q_0 R .
    
    The magnetic field B = ∇× A is computed as
    
        B_R      = - B_0 Zq_0 R ,    
        B_Z      =   B_0  (R - R_0)q_0 R ,    
        B_   = - B_0 R_0R ,    B   = B_0 Sq_0 R ,
    
    and the normalized magnetic field as
    
        b_R      = - ZS ,    
        b_Z      = R - R_0S ,    
        b_   = - q_0 R_0S .
    
    Here, R_0 is the radial position of the magnetic axis, B_0 is the magnetic field at R_0, and _0 is the safety factor, regarded as constant. In all of the examples, these constants are set to R_0 = 2, B_0 = 5 and q = 2, respectively.
    The functions r and S are given by
    
        r    = √( (R - R_0)^2 + Z^2) ,    
        S    = √( r^2 + q_0^2 R_0^2) .
    
    In toroidal coordinates, the functions for the momenta are
    
        ϑ_1 (,)    = A_R () +   b_R () ,    ϑ_2 (,)    = A_Z () +   b_Z () , 
        ϑ_3 (,)    = R ( A_ () +   b_ () ) ,    ϑ_4 (,)    = 0 ,
    
    and those for the forces are computed as
    
        f_1 (,,,)
           = ϑ_1,1 (,)  ^1 + ϑ_2,1 (,)  ^2 + ϑ_3,1 (,)  ^3
         - ∇_1 H (,) , 
        
        f_2 (,,,)
           = ϑ_1,2 (,)  ^1 + ϑ_2,2 (,)  ^2 + ϑ_3,2 (,)  ^3
         - ∇_2 H (,) , 
        
        f_3 (,,,)
           = ϑ_1,3 (,)  ^1 + ϑ_2,3 (,)  ^2 + ϑ_3,3 (,)  ^3
         - ∇_3 H (,) , 
        
        f_4 (,,,)
           = b(x) ·
         - ∇_4 H (,) ,
    
    with the gradient of the Hamiltonian being
    
        ∇_1 H (,)    = ∇_1B() ,    ∇_2 H (,)    = ∇_2B() , 
        ∇_3 H (,)    = ∇_3B() ,    ∇_4 H (,)    =  .
    
    
    We consider four different initial conditions, that of a deeply trapped particle, a barely trapped particle, a barely passing particle and a deeply passing particle (c.f., Figure <ref>). The second and third case are particularly challenging. In all three cases we choose (R,Z,) = (2.5, 0, 0) and set the magnetic moment μ = 0.01. The parallel velocity u and the time step h are chosen as follows:
    
        deeply trapped     barely trapped     barely passing     deeply passing 
    u   0.1   0.3375   0.3425   0.5
    h   5.0   3.0   2.5   2.5
    
    
    
    
    We make the following observations:
    
      *  The passing particles (Figure <ref>) seem to be the more challenging ones of the considered examples. Without projection almost no method is long-time stable and we observe drifts in the energy as well as in the toroidal momentum. The Gauss–Lobatto–IIIA, IIIB and IIIC methods are particularly unstable, crashing after only a few or at best a few hundred time steps. This can be remedied by a smaller time step. Then, however, the computational effort is much larger than with the other methods.
    For the Gauss–Legendre Runge–Kutta methods we observe that for the integrators with an odd number of stages, the solution has two branches, while for the integrators with an even number of stages, the solution has only one branch.
    
    
      *  For the trapped particles (Figure <ref>), we obtain good results without projection for the Gauss–Legendre methods and the SRK3 method except for the six-stage Gauss–Legendre method. For the barely trapped particle, this is also true for the two-stage Gauss–Legendre method. The Gauss–Lobatto–IIIA, IIIB and IIIC methods crash quickly, that is after only about one or a few thousand time steps. The Gauss–Lobatto–IIID and IIIE methods work mostly well for the deeply trapped particle but not for the barely trapped particle.
    
    
      *  The standard projection seems to improve the situation in case of the passing particles (left column of Figure <ref>) and worsen the situation in case of the trapped particles (left column of Figure <ref>). In all cases, though, the results are not satisfactory as the solution always exhibits a drift in the energy error.
    
    
      *  The symplectic projection (left column of Figures <ref> and <ref>) leads to good results for the trapped particles but not for the passing particles. As before, we observe that in those cases where the unprojected solution is stable, the symplectic projection is also stable, but in those cases where the unprojected solution is unstable, the symplectic projection is also unstable.
    
    
      *  The symmetric (right column of Figures <ref> and <ref>) and midpoint (right column of Figures <ref> and <ref>) projections lead to good results with almost all methods except for the Gauss–Lobatto–IIIC methods. For the barely trapped particle both projections are unstable for almost all Gauss–Lobatto methods.
    For the symmetric projection and the higher-order Gauss–Legendre method we observe a drift in the energy error, which again is only of the order of 10^-12 for 1 250 000 time steps. This drift is not due to the non-symplecticity of the symmetric projection method, as one could expect, but it is caused by round-off errors and disappears in simulations in quadruple precision (c.f., Section <ref>). 
    
    
      *  With the Radau methods we observe exact conservation of the toroidal momentum, which constitutes one of the components of the Dirac constraint and is thus expected to be preserved, but we also observe dissipation of energy. For the two-stage Radau method these are related to large errors in the solution (Figures <ref> and <ref>). For the three-stage Radau method (Figures <ref> and <ref>) the errors in the solution are less pronounced.
    
    
    
    The numerical experiments suggest that the symmetric and midpoint projection lead to good results with all integrators, while the symplectic projection is only stable for those integrators and those examples which are stable even without projection.
    Similar to the point vortex example, the solution with the standard projection exhibits a drift in the energy and correspondingly a degradation of the numerical solution, as is expected.
    
    
    
    
     §.§ Energy Drift
    
    
    For several examples we observed a drift of the energy error with the higher-order variational integrators, e.g., for Gauss–Legendre Runge–Kutta methods with four or more stages.
    This behaviour was particularly prominent with the symmetric and midpoint projection when the energy error approached the machine accuracy.
    A natural question to ask is if the origin for this phenomenon lies in the non-symplecticity of the two projection methods or if it is just due to round-off errors.
    To this end, we repeated some simulations, namely those for a barely passing guiding centre particle, in quadruple precision. The resulting energy error is plotted in Figure <ref>. We also show the drift of the energy error, which is obtained by splitting the time interval of the simulation in 10 sub-intervals and computing the maximum of the absolute value of the energy error in each interval. As can be seen, there is no trend like linear growth, etc., visible, indicating that all errors due to either the non-symplecticity of the projection methods or due to round-off are much smaller than the energy error in these simulations.
    
    <cit.> attribute this drift behaviour to inaccuracies in the Runge-Kutta coefficients and weights, leading to an only approximate satisfaction of the symplecticity conditions (<ref>). This leads to a linear growth of the energy error with time, even though one would expect a growth with the square root of the time as follows from random walk arguments (Brouwer's law <cit.>).
    In order to reduce the influence of round-off errors for Runge–Kutta methods of high precision, <cit.> suggest to apply a method inspired by compensated summation, where the Runge-Kutta coefficients and weights are split into two parts,
    
        b_i   = b_i^* + b_i ,    
        a_ij   = a_ij^* + a_ij .
    
    Here, b_i^* and a_ij^* are approximations of b_i and a_ij, which are exact to machine precision, and b_i and a_ij are corrections e.g. to the values of b_i and a_ij in quadruple precision.
    With these coefficients and weights, the definition of the internal stages and the update rule e.g. in (<ref>) are modified to
    
        Q_n,i   = q_n + h ∑_j=1^s a_ij^* Q̇_n,j + h ∑_j=1^sa_ij Q̇_n,j , 
        
        q_n+1   = q_n + h ∑_i=1^s b_i^*Q̇_n,i + h ∑_i=1^sb_i Q̇_n,i ,
    
    which comes at practically no cost but allows to recover the missing accuracy in the last few digits of the solution.
    
    Another technique, recently proposed by <cit.>, consists in choosing the coefficients a̅_ij such that the symplecticity condition (<ref>) holds exactly in floating point precision. This means that even for methods like Gauss–Legendre, where a and a̅ are usually the same, slightly different coefficients are used.
    If this technique is applied, e.g., with Gauss–Legendre Runge–Kutta methods with five and six stages and symmetric projection, the acquired solutions exhibit no drift in the energy error (not shown).
    
    
    
    
     §.§ Convergence
    
    
    We analyse the convergence behaviour of the various integrators together with the various projection methods for the planar point vortex example (c.f., Figure <ref>).
    The orders for the solution error and the energy error are obtained as listed in Table <ref>, while the orders for the angular momentum error are obtained as listed in Table <ref>.
    
    We see that without projection the orders of the integrators are reduced to half the order for ODE problems or half the order plus one in the case of half the order being odd (as all considered methods are symmetric and therefore have even order).
    The standard, symmetric and symplectic projection restore the usual order of 2s for the Gauss–Legendre methods, 2s-2 for the Gauss–Lobatto methods and order 4 for the SRK3 method.
    These are mostly well known results (see e.g. <cit.>).
    The midpoint projection restores the usual order for the Gauss–Legendre method with one stage and the SRK3 method, that is for the integrators for which the midpoint projection is symplectic, as well as the Gauss–Lobatto methods. For the other Gauss–Legendre methods the order is either s+1 or s+2 depending on the method being of odd or even s, respectively.
    
    Interestingly, for the standard, symmetric and symplectic projection methods, the convergence order of the angular momentum error is increased compared to the convergence order for the solution and the energy error, by one for the standard projection and by two for the symmetric and symplectic projection.
    Without projection and with the midpoint projection the convergence of the angular momentum error is the same as the convergence order of the solution and the energy error.
    Only the SRK3 method forms an exception, as here the convergence order of the angular momentum error is decreased by one for the standard projection and by two for the midpoint projection with respect to the order of the other errors.
    
    
    
     §.§ Poincaré Integral Invariants
    
    
    It has recently been established that Poincaré integral invariants provide useful diagnostics for analysing the long-time accuracy of numerical integrators for Hamiltonian dynamics and for distinguishing between symplectic and non-symplectic integrators <cit.>.
    In the following we will apply this diagnostics to the guiding centre example with a simple, symmetric magnetic field.
    
    In particular, we consider the first and second Poincaré integral invariant ϑ and ω, which are given by the Lagrangian one- and two-form, respectively.
    The one-form ϑ is a relative integral invariant, which means that the integral
    
        I_1   = ∫_γϑ_i (q)   dq^i
    
    stays constant in time when γ is a closed loop in the configuration space M (a compact one-dimensional parametrized submanifold of M without boundary), that is advected along the solution of the dynamics.
    Figure <ref> shows examples of single trajectories of some samples of such a loop, as well as the temporal evolution of the whole loop following the dynamics of the guiding centre system.
    The two-form ω is an absolute integral invariant, which means that the integral
    
        I_2   = ∫_Sω_ij (q)   dq^i  dq^j
    
    stays constant in time when S is any compact two-dimensional parametrized submanifold of M, advected along the solution of the dynamics.
    Figure <ref> shows how an initially rectangular area in phasespace is advected by the dynamics of the guiding centre system.
    
    
    
    The loop γ is parametrized by τ∈ [0,1), so that
    
        I_1 (t)    = ∫_0^1ϑ_i (q_(τ)(t))  d q_(τ)^idτ  dτ .
    
    In order to compute this integral, we use N equidistant points in [0,1), so that the derivatives d q_(τ) / dτ can be efficiently computed via discrete Fourier transforms. The integral is approximated with the trapezoidal quadrature rule, which has spectral convergence on periodic domains <cit.>.
    The area S is parametrized by (σ, τ) ∈ [0,1]^2, so that
    
        I_2 (t)    = ∫_0^1∫_0^1ω_ij (q_(σ, τ) (t)  dq_(σ, τ)^idσ dq_(σ, τ)^jdτ  dσ  dτ .
    
    Here, we represent the surface in terms of Chebyshev polynomials and thus use Chebyshev points for the discretization of the domain [0,1]^2. The Chebyshev polynomials allow for an extremely accurate approximation of the surface, even if the latter becomes severely deformed. Moreover, they allow for the use of the  package <cit.> for the easy and accurate computation of the derivatives and the integral.
    As initial conditions we use
    
        q_(τ) (0) = (r_xcos (2πτ),   r_ysin (2πτ),   r_zsin (2π s),   u_0 + u_1cos (2π s))
    
    for the loop γ with r_x = 0.5, r_y = 0.3, r_z = 0.1, u_0 = 0.5, u_1 = 0.05 and
    
        q_(σ, τ) (0) = (r_0 (σ-0.5),   r_0 (τ-0.5),   r_zcos (2πσ) cos (2πτ),   u_0 + u_1sin (2πσ) sin (2πτ) )
    
    for the surface S with r_0 = 0.1, r_z = 0.1, u_0 = 0.5 and u_1 = 0.01. In both cases, the magnetic field is given by B = (0,   0,   B_0  (1 + x^2 + y^2)) with B_0 = 1 and we set μ = 0.01. We use a simpler configuration as in Section <ref> as for the magnetic field used there not all integrators are stable.
    We use 2000 points to discretize the loop γ and 100 × 100 points to discretize the surface S. The time step is h = 10 in both cases. 
    
    
    We make the following observations:
    
      *  The unprojected variational integrators preserve the Poincaré integral invariants with respect to the canonical one-form p_i  q^i and two-form p_i∧ q^i (Figures <ref> and <ref>) but not the invariants with respect to the noncanonical one-form ϑ_i  q^i (Figures <ref>-<ref>) and two-form ω_ij (q)   q^i∧ q^j (Figures <ref>-<ref>).
    
    This behaviour is expected. The variational integrators preserve the canonical forms by construction but not the noncanonical forms as the solution does not satisfy the Dirac constraint under which the two are equivalent.
    However, in cases like this one, where the solution stays close to the constraint submanifold, conservation of the canonical forms implies also approximate conservation of the noncanonical forms. This manifests itself in the error of the integral invariant appearing to be bounded.
    
    In addition, conservation of the canonical forms implies conservation of discrete noncanonical forms, which are obtained by pulling back the canonical forms with the discrete fibre derivatives (<ref>).
    
    
      *  For the standard projection, we observe a clear drift in both the first and second Poincaré integral invariant for the one- and two-stage Gauss–Legendre Runge–Kutta integrators (Figures <ref>-<ref> and <ref>-<ref>) as well as the SRK3 method (Figures <ref> and <ref>). 
    For the higher-order Gauss–Legendre Runge–Kutta integrators (not shown) we see no drift over the runtime of the simulation as the error in the Poincaré integral invariants introduced by the standard projection is smaller than machine accuracy and therefore not measurable.
    
    
      *  For the symmetric projection, both the first and second integral invariant is exactly preserved by all three of the considered methods (Figures <ref>-<ref> and <ref>-<ref>). 
    We only observe some “diffusion” of the error. This is expected due to the deformation of the loop and the surface, over which the integrals are evaluated, leading to a degradation of accuracy in the integrals.
    
    
      *  The symplectic projection preserves neither the first (Figures <ref>-<ref>) nor the second (Figures <ref>-<ref>) integral invariant exactly. Instead it shows a behaviour of the error of the invariant similar to the unprojected methods albeit at greatly reduced amplitude (for the two-stage Gauss–Legendre Runge–Kutta integrator the error is already reduced below the machine accuracy, so the invariants appear to be preserved).
    This behaviour is expected, though, as the symplectic projection does not preserve the standard noncanonical one- and two-form but the modified forms (<ref>) and (<ref>). 
    Taking the corresponding “correction terms” to the standard one- and two-form into account, we find exact conservation of both the first and the second invariant (Figures <ref> and <ref>).
    
    
      *  For the midpoint projection, we find exact conservation of both the first and second invariant for all three integrators (Figures <ref>-<ref> and <ref>-<ref>), even for the two-stage Gauss–Legendre Runge–Kutta method, for which the midpoint projection is not symplectic.
    
    
    We have seen that the Poincaré invariant diagnostics provides a useful tool in experimentally judging the symplecticity of a numerical algorithm.
    We see that the unprojected variational integrators only preserve the invariants with respect to the canonical one- and two-form, but not with respect to the noncanonical forms.
    With the standard projection all integrators exhibit a clear drift in both invariants.
    With the symmetric and the midpoint projection, this drift is essentially gone, even though the former is not symplectic.
    Similarly, with the symplectic projection and all the integrators, we find exact conservation of the corrected invariant.
    Yet, these results should be taken with a grain of salt. The test case considered here is relatively simple. It is conceivable that in more complicated situations larger errors in the Poincaré invariants will be observed with the various projection methods. It is still expected, though, that the errors are bounded with all three of the symmetric, midpoint and symplectic projection methods.
    If it is conceived necessary, it should always be possible to reduce the errors in the invariants to the level of the machine accuracy by reducing the time step, where in most cases we expect only a mild reduction to be necessary.
    When the errors in the invariants are smaller than the machine accuracy, the projected methods essentially behave like any other nonlinearly implicit symplectic integrator, which needs to be solved in finite precision and is thus never exactly symplectic (see <cit.> for details).
    
    
    
    
     §.§ Symplectic Projection
    
    
    For many integrators we observed that the symplectic projection is unstable, even when it does not amount to the post projection method.
    The reason for this is that in general there are no bounds on the amplitude of Lagrange multiplier λ that prevent the unprojected solution from moving further and further away from the constraint submanifold Δ.
    This behaviour is illustrated in Figure <ref> and prominently exemplified by the passing guiding centre particles as shown in Figure <ref>. Here, the Lagrange multiplier oscillates between positive and negative values while its amplitude grows without bound.
    Eventually, the projection step becomes larger than the integration step and the simulation becomes unstable.
    A deeper understanding of this behaviour is very much desirable. One could, for example, apply backward error analysis and try to understand under which conditions the modified equation for the evolution of λ is stable. There are also indications of a connection with index reduction of the continuous problem. Such investigations, however, are left for future work.
    
    
    
    
    § SUMMARY
    
    
    We have devised several projection methods and analyzed their influence on the symplecticity and long-time stability of variational integrators applied to degenerate Lagrangian systems.
    The corresponding system of equations constitutes a system of differential-algebraic equations of index two, for which standard symplectic integrators like Gauss–Legendre or Gauss–Lobatto Runge–Kutta methods are well known to deliver poor performance. In particular their order of accuracy is reduced severely compared to their order when applied to ordinary differential equations.
    
    This underperformance can be remedied by the application of appropriate projection methods. In the context of symplectic or variational integrators, this approach raises the question what influence such a projection method has on the symplecticity and long-time stability of the resulting integrator.
    While for simple problems all projection methods, and even some of the unprojected integrators, lead to long-time stable simulations, the only universally stable methods have been found to be the symmetric and midpoint projection.
    In most examples, the standard projection admits a drift in the solution and thus the energy error, rendering the simulation unstable in finite time. When exactly that happens depends strongly on the order of the underlying numerical integrator and thus the error in the algebraic constraint.
    Among the two stable projection methods, the symmetric projection is usually preferable, as the midpoint projection leads to integration methods of reduced order, whereas the symmetric projection restores the order of the underlying numerical integrator.
    On the other hand, the midpoint projection is exactly symplectic for the midpoint and SRK3 methods.
    Even though the symmetric projection method is not exactly symplectic, we found the error in the symplecticity condition to be so small, that it has no practical influence on the long-time stability of the projected variational integrator.
    
    We also analyzed a modification of the symmetric projection method that is symplectic in a generalized sense, namely in that it is preserving a generalized symplectic structure on a larger space. While for some problems this method leads to similarly good results as the symmetric method, it fails for others. The reason for this failure is that the Lagrange multiplier can grow without bounds. If that happens, the simulation becomes unstable.
    However, in those cases where the symplectic projection can be applied, it is preferable over the symmetric projection as it has lower computational cost. For the symplectic projection, the underlying numerical integrator and the projection step can be performed independently, whereas in the symmetric case the whole system has to be solved at once.
    
    
    
    
    
    § ACKNOWLEDGMENTS
     Helpful discussions with Joshua Burby, Leland Ellison, Melvin Leok and Tomasz Tyranowski are gratefully acknowledged.
    Moreover, the author is indebted to Omar Maj and Hiroaki Yoshimura for inspiring conversations as well as reading a draft of the paper. 
    The author has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska–Curie grant agreement No 708124. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
    
    
    
    tocsectionReferencesplainnat
