Variational order for forced Lagrangian systems
D. Martín de Diego R. Sato Martín de Almagro
math-ph math.DG math.MP math.NA math.SG 17B66, 22A22, 70G45, 70Hxx





Variational order for forced Lagrangian systems
     D. Martín de Diego, R. Sato Martín de Almagro


Instituto de Ciencias Matemáticas (CSIC-UAM-UC3M-UCM) 

C/Nicolás Cabrera 13-15, 28049 Madrid, Spain
    December 30, 2023
=========================================================================================================================================================


We are able to derive the equations of motion for forced mechanical systems in a purely variational setting, both in the context of Lagrangian or Hamiltonian mechanics, by duplicating the variables of the system as introduced by <cit.>. We show that this construction is useful to design high-order integrators for forced Lagrangian systems and, more importantly, we give a characterization of the order of a method applied to a forced system using the corresponding variational order of the duplicated one.






§ INTRODUCTION

In recent years, a new area of numerical analysis has been developed with the idea of designing numerical methods to integrate the evolution of ordinary and partial differential equations preserving, as much as possible, the qualitative features of the original dynamics (see <cit.>). This area is known as geometric integration. 

An important property to preserve for variational Lagrangian and Hamiltonian systems is symplecticity. A particularly elegant way to derive symplectic integrators is by discretizing Hamilton's principle, i.e. discretizing the variational principle first instead of the equations of motion. Numerical methods derived using this procedure are called variational integrators and automatically preserve symplecticity and momentum and exhibit quasi-energy conservation for exponentially long times (see <cit.> and references therein).

One of the difficulties is obtaining high-order methods since usually many of the typical examples derived in the literature are second order schemes, derived from the midpoint or trapezoidal rules essentially, but it is important to study the construction of high-order symplectic schemes (see <cit.>). This is the objective of some recent papers in the last years (see for instance <cit.>).

A key ingredient that considerably lowers the difficulty of proving the order of an integrator is the variational error theorem stated in <cit.>. The proof given there was correctly formulated and proved by <cit.>. In consequence, the equivalent proof for forced systems given in <cit.> was an open question. 

In this paper, apart from laying a geometric framework for Galley's theory both in Lagrangian and Hamiltonian sides, we observe that by using the duplication technique introduced by <cit.> we can give an answer to the question of the order of forced Lagrangian systems using the proof for free Lagrangian systems given by <cit.> (see also <cit.>).  

In particular, we can derive numerical methods for forced systems using the same methods that were previously developed for standard variational integrators. Preliminary numerical tests of this principle under the name of slimplectic integrators can be found in <cit.> and our approach can be used to rigorously derive the order of the methods previously developed for forced systems (see, for instance, <cit.> and references therein). 



§ FORCED LAGRANGIAN AND HAMILTONIAN SYSTEMS

A Lagrangian system is defined by a smooth manifold Q, the configuration manifold, and a smooth function L: TQ→ℝ, the Lagrangian function, which determines the dynamics. Here TQ denotes the tangent bundle of Q with canonical projection τ_Q: TQ→ Q. 
Let us consider local coordinates (q^i) on Q, 1, …, n= Q. The corresponding fibered coordinates on TQ will be denoted by (q^i, v^i) where τ_Q(q^i, v^i)=(q^i). 

Consider curves c:[a, b]⊆ℝ→ Q of class C^2 connecting two fixed points q_0, q_1∈ Q and the collection of all these curves

    C^2(q_0, q_1, [a,b])={ c: [a, b] → Q  |   c∈ C^2([a, b]), c(a)=q_0, c(b)=q_1}  .

whose tangent space is

    T_c C^2(q_0, q_1, [a,b])={ X: [a, b]→ TQ  |  X∈ C^1([a, b]), τ_Q∘ X=c    and    X(a)=X(b)=0}  .

Given the Lagrangian function L: TQ→ℝ we define the functional: 

    [                     𝒥:   C^2(q_0, q_1, [a,b])                      ⟶                      ℝ;                                             c                      ⟼ ∫_a^b L(c(t), ċ(t)) dt ]

(Hamilton's principle)  A curve c∈ C^2(q_0, q_1, [a, b]) is a solution of the Lagrangian system defined by L: TQ→ℝ if and only if c is a critical point of the functional 𝒥, i.e. d𝒥(c)=0
 
 
Using standard techniques from variational calculus, it is easy to show that the curves c(t)=(q^i(t)) solutions of the Lagrangian system defined by L are the solutions (c(t), ċ(t)), with = d/dt, of the following system of second-order implicit  differential equations:

    d/dt( ∂ L/∂ v^i)-∂ L/∂q^i=0  .

which are the well-known Euler-Lagrange equations (see <cit.>).

Let us denote by
S = dq^i ⊗∂/∂ v^i
and
Δ = v^i  ∂/∂ v^i
the vertical endomorphism and the Liouville vector field on TQ (see <cit.> for intrinsic definitions).

The Poincaré-Cartan 2-form is defined by
ω_L = - dθ_L   ,  θ_L = S^*(dL)
and the energy function
E_L = Δ(L)-L
which in local coordinates read as

    θ_L    =   ∂ L/∂ v^i d q^i 
    ω_L    =   dq^i ∧d( ∂ L/∂ v^i)
    
    E_L    =    v^i ∂ L/∂ v^i - L (q, v)  .

 Here S^* denotes the adjoint operator of S.

The Lagrangian L is said to be regular if the Hessian matrix whose entries are

    W_ij = ∂^2 L/∂ v^i ∂ v^j

is regular, and in this case, ω_L is a symplectic form on TQ. We will assume in the sequel that L is regular. 

The Euler-Lagrange equations are geometrically encoded as the equations for the flow of the vector field X_E_L: 

    _X_E_Lω_L = dE_L  .


An external force is expressed as a map F: TQ→ T^*Q such that π_Q ∘ F=τ_Q where τ_Q: TQ → Q is the canonical projection. In coordinates F(q^i, v^i)=(q^i, F_i(q, v)). 

    TQ [rr]^F[rd]_τ_Q       T^*Q[ld]^π_Q
       Q

where T^*Q is the cotangent bundle with canonical projection π_Q: T^*Q→ Q.

Given a force we can construct a semibasic 1-form μ_F∈Λ^1(TQ) by ⟨μ_F(u_q) , W_u_q⟩=⟨ F(u_q), Tτ_Q (W_u_q)⟩, for all W_u_q∈ T_u_qTQ. In coordinates

    μ_F=F_i(q^i, v^i) dq^i  .

Construct the vertical vector field Z^v_F given by _ Z^v_Fω_L=-μ_F. In coordinates

    Z^v_F=W^ij F_j∂/∂ v^i  ,

where W^ij denote the components of the inverse of the Hessian matrix W.
The dynamics of the forced Lagrangian system is determined by the vector field X_E_L+Z^v_F: 

    _X_E_L+Z^v_Fω_L = dE_L -μ_F  .

The integral curves of X_E_L+Z^v_F are the solutions of the forced Euler-Lagrange equations: 

    d/dt( ∂ L/∂ v^i)-∂ L/∂q^i=F_i  .


Now, we move on to the Hamiltonian description of systems subjected to external forces. The cotangent bundle T^*Q of a differentiable manifold Q is equipped with a canonical exact symplectic structure ω_Q=dθ_Q, where θ_Q is the canonical 1-form on T^*Q defined point-wise by

    (θ_Q)_α_q(X_α_q)=⟨α_q, T_α_qπ_Q(X_α_q)⟩

where X_α_q∈ T_α_qT^*Q, α_q∈ T_q^*Q.

In canonical bundle coordinates (q^i, p_i) on T^*Q the projection reads as π_Q(q^i, p_i)=(q^i), and 

    θ_Q= p_i dq^i  , ω_Q= dq^i∧dp_i  .


Given a Hamiltonian function H: T^*Q→ℝ we define the Hamiltonian vector field

    _X_Hω_Q=dH

which can be written locally as

    X_H=∂ H/∂ p_i∂/∂ q^i- ∂ H/∂ q^i∂/∂ p_i  .

 Its integral curves are determined by Hamilton's equations: 

    dq^i/dt   =   ∂ H/∂ p_i  ,
    dp_i/dt   =   -∂ H/∂ q^i  .


Let TQ be the tangent bundle of Q with canonical projection τ_Q: TQ→ Q and canonical coordinates (q^i, v^i) where τ_Q(q^i, v^i)=(q^i). 

Given a function H: T^*Q→ℝ, a Hamiltonian function, we may construct the transformation 𝔽H: T^*Q→ TQ where ⟨𝔽H(α_q), β_q⟩=d/dt |_ t=0H(α_q+tβ_q).
In coordinates, 𝔽H(q^i, p_i)=(q^i, ∂ H/∂ p_i(q,p)). 
We say that the Hamiltonian is regular if 𝔽H is a local diffeomorphism, which in local coordinates is equivalent to the regularity of the Hessian matrix whose entries are:

    M^i j = ∂^2 H/∂ p_i∂ p_j.

Consider now the external forcing previously defined in the Lagrangian description and denote F^H = 𝔽H^* F: T^*Q→ T^*Q. 

    T^*Q [rr]^F^H[rd]_π_Q       T^*Q[ld]^π_Q
       Q

It is possible to modify Hamiltonian vector field X_H to obtain the forced Hamilton's equations as the integral curves of the vector field X_H+Y^v_F where the vector field
Y^v_F∈𝔛(T^*Q) is defined by

    Y^v_F(α_q)=d/dt|_t=0(α_q+t F^H(α_q))  .

We will say the the forced Hamiltonian system is determined by the pair (H, F^H).


Locally, 

    Y^v_F=F_i(q^j, ∂ H/∂ p_j(q, p))∂/∂ p_i=F^H_i(q, p)∂/∂ p_i

modifying Hamilton's equations as follows:
 
    dq^i/dt   =   ∂ H/∂ p_i(q,p)  ,
    dp_i/dt   =   -∂ H/∂ q^i(q,p)+F^H_i(q,p)  .




§ FORCED HAMILTONIAN DYNAMICS AS FREE DYNAMICS OBTAINED BY DUPLICATION
  
In this section we define a new unforced Hamiltonian system whose dynamical equations are related with the forced system (H,F^H). In a final appendix we will expand on the interesting geometry (that of a cotangent groupoid) that arises with the duplication of variables. 

The idea is to duplicate variables in such a way that system (<ref>) transforms into a different free Hamiltonian system. This opens up the opportunity to analyze the original system using standard Hamiltonian techniques. 

Consider the Cartesian product Q × Q and its cotangent bundle T^*(Q× Q) ≡ T^*Q × T^*Q with canonical projections pr_1: T^*Q × T^*Q→ T^*Q and pr_2: T^*Q× T^*Q→ T^*Q given respectively by pr_1(α_q, β_q')=α_q and pr_2(α_q, β_q')=β_q' where α_q∈ T^*_qQ and β_q'∈ T^*Q_q'. Also in this space we have the inversion map ι: T^*Q × T^*Q → T^*Q × T^*Q (cf. section <ref>), which in local coordinates acts as ι(α_q, β_q') = (β_q', α_q).

Let us endow T^*Q × T^*Q with the following symplectic structure:

    Ω_Q× Q =pr_2^*ω_Q-pr_1^*ω_Q  .

Observe that the mapping 

    [           Ψ:   T^*Q× T^*Q            ⟶    T^*(Q× Q);               (α_q, β_q')            ⟼ (-α_q, β_q') ]

is a symplectomorphism, i.e. Ψ^*(Ω_Q× Q)=ω_Q× Q where ω_Q× Q is the standard symplectic form on T^*(Q× Q).

Consider the Hamiltonian 𝐇: T^*Q× T^*Q→ℝ defined by

    𝐇(α_q, β_q') = (H ∘pr_2 - H ∘pr_1)(α_q, β_q') = H(β_q')- H(α_q)  ,
 
satisfying 𝐇∘ι = - 𝐇, and the identity map (cf. section <ref>)

    [         ϵ:       T^*Q          → T^*Q× T^*Q;                   α_q          ⟼ (α_q, α_q) ]

Observe that the spaces T^*Q and ϵ(T^*Q) are canonically diffeomorphic. We have the following 
 

The Hamiltonian vector field X_𝐇 given by

    _X_𝐇Ω_Q× Q=d𝐇

 verifies that:


  *  X_𝐇 is tangent to ϵ(T^*Q);

  *  .X_𝐇|_ϵ(T^*Q)=ϵ_* (X_H).



See Proposition <ref>.

 
Observe that the proof of Lemma <ref> is quite straightforward using coordinates. In fact, if we take adapted coordinates (q^i, p_i; Q^i, P_i) on T^*Q × T^* Q then 
     
    𝐇(q^i,  p_i; Q^i, P_i) = H(Q^i, P_i)- H(q^i, p_i)  ,
                  
and so                         
  
    X_𝐇= ∂ H/∂ p_i(q,p)∂/∂ q^i- ∂ H/∂ q^i(q,p)∂/∂ p_i+  ∂ H/∂ P_i(Q,P)∂/∂ Q^i- ∂ H/∂ Q^i(Q,P)∂/∂ P_i  .

which is obviously tangent to ϵ(T^*Q), this last space is locally given by the vanishing of the 2n-constraints Q^i-q^i=0 and P_i-p_i=0 and moreover <ref><ref> follows immediately since ϵ (q^i, p_i)=(q^i, p_i; q^i, p_i). 



 Define  𝔽𝐇^×(α_q, β_q')=(-𝔽H(α_q), 𝔽H(β_q')) with Hessian matrix: 
  
    (
    [ -∂^2 H/∂ p_i∂ p_j (q,p)                  0_n× n;                  0_n× n  ∂^2 H/∂ P_i∂ P_j (Q,P) ])  .


The following lemma is trivial but it will be interesting for us when going to the Lagrangian side

The transformation 𝔽𝐇^×: T^*Q× T^*Q→ TQ× TQ is a local diffeomorphism if and only if 𝔽H: T^*Q→ TQ is a local diffeomorphism 

Given a Hamiltonian 𝐇, we may want to add a function K: T^*Q× T^*Q →ℝ. In light of the result of lemma <ref>, if this function has the property K(β_q', α_q) = -K(α_q, β_q'), then 𝐇_K = 𝐇 + K will still preserve that property and the trajectories of the resulting dynamics at the identities will remain bound to the identities. This allows a richer behaviour of the original system in T^*Q, and for this, functions K satisfying said property are called generalized potentials <cit.>. 

The previous results are only a preparation for our real objective, to find a purely Hamiltonian representation of systems with forces using this duplication of variables. For that, we need to introduce the definition of retraction (see <cit.>): 

 
A retraction on a manifold Q is a smooth mapping R: TQ→ Q satisfying the following properties. If we denote by R_q=.R|_T_qQ the restriction of R to T_qQ then: 


  * R_q(0_q)=q, 

  * Identifying T_0_q T_q Q ≡  T_q Q then 
T_0_qR_q =id_T_qQ.



In what follows, it will be interesting to introduce σ: TQ→ Q× Q defined by
σ(v_q)=(q, R_q(v_q)). From the previous definition it is easy to show that σ is invertible in some neighborhood of 0_q∈ T_qQ for any q∈ Q. Denote this local inverse by τ: Q× Q→ TQ which in coordinates will take the form

    τ(q^i, Q^i)=(q^i, τ^i (q,Q))  .



Consider a chart (U, φ) around a point q∈ Q then 
the map
σ: TU→ Q× Q defined by σ(v_q)=(q, R_q(v_q)) where q∈ U verifies that the map
T_0_qσ in coordinates (TU, Tφ) and (U× U, φ×φ) is 

    (
    [ I 0; I I ])

and, in consequence, the map T_(q,q)τ is represented by the matrix

    (
    [  I  0; -I  I ])



We have that in the chosen coordinates σ(q^i, v^i)=(F^i(q, v), G^i(q,v))=(q, R_q(v)). Therefore, in coordinates, 

    ∂ F^i/∂ q^j(q, 0)=δ^i_j,   ∂ F^i/∂ v^j(q, 0)=0,

where δ^i_j is the Kronecker delta, and observe that

    ∂ G^i/∂ q^j(q, 0)=δ^i_j

since R_q'(0)=q' for q'∈ U and 

    ∂ G^i/∂ v^j(q, 0)=δ^i_j

since T_0_qR_q =id_T_qQ.


Typically, we can induce this kind of mappings using an auxiliary Riemannian metric g on Q with associated geodesic spray Γ_g (see <cit.>). The associated exponential for a small enough neighborhood U⊂ TQ of 0_q.

    [      exp^Γ_g:         U⊂ TQ             ⟶          Q× Q;                         v_q             ⟼ (q, γ_v_q(1)) ]

where t→γ_v_q(t) is the unique geodesic such that γ'_v_q(0)=v_q. 
For instance when Q=ℝ^n and we take the Euclidean metric, we induce the map

    τ(q, q')=(q, q'-q)  and  σ(q, v)=(q, q + v)  .


Given a forced Hamiltonian system (H, F^H) we can construct the function K_F: T^*Q× T^*Q→ℝ as: 

    K_F(α_q, β_q')= 1/2⟨ F^H(β_q'), τ(q', q)⟩ - 1/2⟨ F^H(α_q), τ(q, q')⟩

Observe that this function satisfies the important property K_F(β_q', α_q)=-K_F(α_q, β_q'), and thus is a generalized potential.
Consider the Hamiltonian 𝐇_K_F: T^*Q× T^*Q→ℝ defined by
 
    𝐇_K_F(α_q, β_q') = 𝐇(α_q, β_q') + K_F(α_q, β_q')  ,
 
which also satisfies 𝐇_K_F(β_q', α_q) = -𝐇_K_F(α_q, β_q') by construction.

 The Hamiltonian vector field X_𝐇_K_F given by
 
    _X_𝐇_K_FΩ_Q× Q=d𝐇_K_F

 verifies that:


  *  X_𝐇_K_F is tangent to ϵ(T^*Q);

  *  .X_𝐇_K_F|_ϵ(T^*Q)=ϵ_* (X_H+Y^v_F).

 
 
Part <ref> is again a direct consequence of proposition <ref>.

To deduce part <ref> observe that 
 
    X_𝐇_K_F   =   (∂ H/∂ p_i(q,p)+1/2∂ F_j^H/∂ p_i(q,p)τ^j(q, Q) )∂/∂ q^i
          -( ∂ H/∂ q^i(q,p)+1/2∂ F_j^H/∂ q_i(q,p)τ^j(q, Q).
          .+1/2 F_j^H(q,p)∂τ^j/∂ q^i(q, Q)-1/2 F_j^H(Q,P)∂τ^j/∂ q^i (Q, q) )
      ∂/∂ p_i
          + (∂ H/∂ P_i(Q,P)+1/2∂ F_j^H/∂ P_i(Q,P)τ^j(Q, q) )∂/∂ Q^i
          -( ∂ H/∂ Q^i(Q,P)-1/2 F_j^H(q,p)∂τ^j/∂ Q^i(q, Q).
          .+1/2∂ F_j^H/∂ Q_i(Q,P)τ^j(Q, q)+1/2 F_j^H(Q,P)∂τ^j/∂ Q^i (Q, q) )
      ∂/∂ P_i


Now using Lemma <ref> we have that along the identities ϵ(T^*Q)
 
    X_𝐇_K_F   =   ∂ H/∂ p_i(q,p) ∂/∂ q^i
      -( ∂ H/∂ q^i(q,p)
       -F_j^H(q,p) )
      ∂/∂ p_i
          + ∂ H/∂ p_i(q,p) ∂/∂ Q^i
      -( ∂ H/∂ q^i(q,p)
      - F_j^H(q,p) )
      ∂/∂ P_i

and thus .X_𝐇_K_F|_ϵ(T^*Q)=ϵ_* (X_H+Y^v_F) as we wanted to prove. 
 

 Define the mapping 𝔽𝐇^×_K_F: T^*Q× T^*Q→ TQ× TQ given in local coordinates as

    𝔽𝐇^×_K_F(q^i, p_i, Q^i, P_i)=(q^i, -∂𝐇_K_F/∂ p_i,  Q^i, ∂𝐇_K_F/∂ P_i).



If H is regular then the transformation 𝔽𝐇^×_K_F: T^*Q× T^*Q→ TQ× TQ is a local diffeomorphism in a neighborhood of ϵ(T^*Q). 
 
  

Locally if we take coordinates (q^i, p_i, Q^i, P_i) then from the definition of K_F we observe that 

    .( 
    [ ∂^2 K_F/∂ p_i∂ p_j ])|_ϵ (T^*Q)=.( 
    [ ∂^2 K_F/∂ P_i∂ P_j ])|_ϵ (T^*Q)=.( 
    [ ∂^2 K_F/∂ p_i∂ P_j ])|_ϵ (T^*Q)
    =0

therefore from the regularity of H it is trivial to derive the regularity of 𝐇_K on a tubular neighborhood of ϵ(T^*Q).









§ FORCED LAGRANGIAN DYNAMICS AS FREE DYNAMICS OBTAINED BY DUPLICATION

Now we will define a new free Lagrangian system whose dynamical equations are related with the forced system (L,F). This is precisely the path chosen by <cit.> although we will not enter into details as to how and why he arrived at his formulation.

Consider the Cartesian product Q × Q and its tangent bundle T(Q× Q) ≡ TQ × TQ with canonical projections pr_1: TQ× TQ→ TQ and pr_2: TQ× TQ→ TQ. In local coordinates we have pr_1(v_q, V_Q) = v_q and pr_2(u_q, V_Q) = V_Q. Consider also the maps ι̃: TQ × TQ → TQ × TQ and ϵ̃: TQ → TQ × TQ locally defined as ι̃(v_q,V_Q)=(V_Q, v_q) and ϵ̃(v_q)=(v_q, v_q), respectively. In this new bundle TQ × TQ ≡ T(Q × Q) we also have a vertical endomorphism and a Liouville field, whose local presentation in adapted coordinates (q^i,v^i; Q^i,V^i) would be S = dq^i ⊗∂/∂ v^i + dQ^i ⊗∂/∂ V^i, = v^i ∂_v^i + V^i ∂_V^i.

Define a new Lagrangian 𝐋: TQ× TQ→ℝ as:

    𝐋(v_q, V_Q) = (L ∘pr_2 - L ∘pr_1)(v_q, V_Q) = L(V_Q) - L(v_q)  ,

Much like 𝐇, this new Lagrangian satisfies that 𝐋∘ι̃ = - 𝐋.

From this new Lagrangian we can define corresponding Poincaré-Cartan forms on TQ × TQ, θ_𝐋 and ω_𝐋, as in Section 2. We can also define a new fibre derivative 𝔽𝐋^×(v_q, V_Q)=(-𝔽L(v_q), 𝔽L(V_Q)). This leads us to state an analogue of lemma <ref>:

The transformation 𝔽𝐋^×: TQ × TQ → T^*Q × T^*Q is a local diffeomorphism if and only if 𝔽L: TQ→ T^*Q is a local diffeomorphism.


It is easy to check that with this definition of fibre derivative the following diagram commutes:

    TQ × TQ [rr]^𝔽𝐋^×       T^*Q × T^*Q
    TQ [u]^ϵ̃[rr]^𝔽L       T^*Q [u]^ϵ

We can also define the energy of the system as in the usual case, with E_𝐋 = (𝐋) - 𝐋, but in order to relate this with the Hamiltonian formulation, it will be useful to rewrite it as:

    E_𝐋(v_q, V_Q) = ⟨𝔽𝐋^×(v_q, V_Q), (v_q, V_Q)⟩_× - 𝐋(v_q, V_Q).

where ⟨·,·⟩_×: (T^*Q × T^*Q) × (TQ × TQ) →ℝ is the inner product defined as:

    ⟨α , v ⟩_× = ⟨pr_2(α), pr_2(v)⟩ - ⟨pr_1(α), pr_1(v)⟩.


The following results will help us prove the analogue of lemma <ref>.


Let 𝐋: TQ × TQ →ℝ be such that 𝐋∘ι̃ = - 𝐋 then
the following diagram commutes:

    TQ × TQ [d]^ι̃[r]^𝔽𝐋^×    T^*Q × T^*Q [d]^ι
    TQ × TQ [r]^𝔽𝐋^×    T^*Q × T^*Q



 We need to show that:

    𝔽𝐋^×∘ι̃ = ι∘𝔽𝐋^×.

On the left-hand side we have

    (𝔽𝐋^×∘ι̃)(q,v,Q,V)    = 𝔽𝐋^× (Q,V,q,v)
       = (Q, -∂𝐋/∂ V(Q,V,q,v),q,∂𝐋/∂ v(Q,V,q,v)),

while on the right-hand side we have

    (ι∘𝔽𝐋^×)(q,v,Q,V)    = ι(q,-∂𝐋/∂ v(q,v,Q,V),Q,∂𝐋/∂ V(q,v,Q,V))
       = (Q,∂𝐋/∂ V(q,v,Q,V),q,-∂𝐋/∂ v(q,v,Q,V))


Now, using that 𝐋∘ι̃ = - 𝐋 we find:

    ∂𝐋/∂ V(Q,V,q,v)    = - ∂𝐋/∂ V(q,v,Q,V),
    ∂𝐋/∂ v(Q,V,q,v)    = - ∂𝐋/∂ v(q,v,Q,V).


Applying this we immediately arrive at the desired result.




The inner product ⟨·,·⟩_× satisfies that

    ⟨ι(α),ι̃(v)⟩_× = - ⟨α,v⟩_×.




First note that pr_1 ∘ι = pr_2 and pr_1 ∘ι̃ = pr_2, and that the same holds under the exchange 1 ↔ 2.

Clearly:

    ⟨ι(α),ι̃(v)⟩_×   = ⟨ (pr_2 ∘ι)(α), (pr_2 ∘ι̃)(v)⟩ - ⟨ (pr_1 ∘ι)(α), (pr_1 ∘ι̃)(v)⟩
       = ⟨pr_1(α), pr_1(v)⟩ - ⟨pr_2(α), pr_2(v)⟩
        = - ⟨α , v ⟩_×





Let 𝐋 be a Lagrangian satisfying  such that 𝐋∘ι̃ = - 𝐋 then also E_𝐋∘ι̃ = - E_𝐋.



Applying the inversion to the definition of the energy given in terms of 𝔽𝐋^× we have:

    (E_𝐋∘ι̃)(·)    = ⟨ (𝔽𝐋^×∘ι̃) (·), ι̃(·)⟩_× - (𝐋∘ι̃)(·)
       = ⟨ (ι∘𝔽𝐋^×) (·), ι̃(·)⟩_× - (𝐋∘ι̃)(·)

Applying lemma <ref> and the inversion property of 𝐋 we get:

    (E_𝐋∘ι̃)(·)    = - ⟨𝔽𝐋^×(·), ·⟩_× + 𝐋(·)
       = - E_𝐋(·)





Let 𝐋 be a regular Lagrangian satisfying the hypothesis of proposition <ref>, and define its associated Lagrangian by the expression 𝐇∘𝔽𝐋^× = E_𝐋. Then 𝐇∘ι = - 𝐇.



Using proposition <ref> we have that 𝐇∘ι∘𝔽𝐋^× = 𝐇∘𝔽𝐋^×∘ι̃ = E_𝐋∘ι̃. Applying proposition <ref> the result follows immediately.


Finally we can state the following result:

Assume L is a regular Lagrangian, then the Hamiltonian vector field X_𝐋 associated to 𝐋 given by

    _X_𝐋ω_𝐋=dE_𝐋

verifies that:


  *  X_𝐋 is tangent to ϵ̃(TQ);

  *  .X_𝐋|_ϵ̃(TQ)=ϵ̃_* (X_L).




It is easy to check that (𝔽𝐋^×)^* Ω_Q × Q = ω_𝐋. Defining H = E_L∘(𝔽L)^-1 we then get that 𝐇 = E_𝐋∘(𝔽𝐋^×)^-1. Thus the results of lemma <ref> also apply to X_𝐋 and these results can be brought back to TQ × TQ, proving our claim.


As with the Hamiltonian formulation we may also include potentials in our description. Again, let K : TQ × TQ →ℝ be a function such that K∘ι̃ = -K, then K is a generalized potential and 𝐋_K = 𝐋 - K satisfies 𝐋_K∘ι̃ = -𝐋_K.

Given a forced Lagrangian system (L, F) we may define the generalized potential K_F: TQ× TQ→ℝ explicitly written as:

    K_F(v_q, V_Q) = 1/2⟨ F(V_Q), τ(Q, q)⟩ - 1/2⟨ F(v_q), τ(q, Q)⟩.


Note that if 𝐋_K_F is regular we may obtain a Hamiltonian from its energy as 𝐇_K_F = E_𝐋_K_F∘(𝔽𝐋^×_K_F)^-1 but in general it will not be the same Hamiltonian as we defined in the previous section, i.e. 

    𝐇_K_F≠ H ∘pr_2 - H ∘pr_1 + 1/2⟨ F^H ∘pr_2, τ∘π_Q × Q∘ι⟩ - 1/2⟨ F^H ∘pr_1, τ∘π_Q × Q⟩

in general. This will only be equal if F does not depend on velocities, which means 𝔽𝐋^×_K_F = 𝔽𝐋^×. This means we cannot directly invoke the result from theorem <ref> to prove that the resulting Euler-Lagrange field coincides with the forced dynamics at the identities and instead we must work a bit more to get the same result. Still at the end of the section we will show that we can actually relate both dynamics obtained from 𝐇_K_F and 𝐇_K_F.

Let us begin with this regularity result:

If L is regular then 𝐋_K_F is regular in a neighborhood of ϵ̃(TQ).



The proof is essentially the same as that of proposition <ref>.


Given 𝐋_K = 𝐋 - K, with 𝐋 = L ∘pr_2 - L ∘pr_1, let us reserve X_𝐋 for the free Euler-Lagrange field, i.e., the vector field that satisfies 
    _X_𝐋ω_𝐋 = dE_𝐋,
 and define the vector field Y_K as the one resulting from the decomposition X_𝐋_K = X_𝐋 - Y_K. Similarly, let us define θ_K by θ_𝐋_K = θ_𝐋 - θ_K and E_K by E_𝐋_K = E_𝐋 - E_K. Clearly if ω_K = - dθ_K then ω_𝐋_K = ω_𝐋 - ω_K, and it is not difficult to show that:

    _Y_Kω_𝐋_K   = d E_K - _X_𝐋ω_K
       = ℒ_X_𝐋θ_K - dK  .



We can now prove the following theorem:


Given 𝐋_K = 𝐋 - K regular in a neighbourhood of ϵ(TQ), with 𝐋 = L ∘pr_2 - L ∘pr_1, whose Hamiltonian vector field X_𝐋_K = X_𝐋 - Y_K satisfies:

    _X_𝐋_Kω_𝐋_K = d E_𝐋_K

then:


  *  X_𝐋_K is tangent to ϵ̃(TQ);

  *  Y_K is vertical and such that ι̃_* Y_K = Y_K.

Furthermore, if K = K_F, then we have that:


  *  ._Y_K_Fω_K_F|_ϵ̃(TQ) = 0;

  *  ._Y_K_Fω_𝐋|_ϵ̃(TQ) = ( pr_1^* F-pr_2^* F )|_ϵ̃(TQ).

  *  .X_E_𝐋_K|_ϵ̃(TQ) = ϵ̃_* (X_E_L + Z^v_F)






  * 
Geometrically, given the regularity in a neighborhood of the identities ω_𝐋_K is non-degenerate there. As E_𝐋_K∘ι̃ = - E_𝐋_K, we are in a position to apply proposition <ref>, rendering this equivalent to theorem <ref>.<ref>. Variationally the proof is a consequence of <ref>.

  *  We know that X_𝐋_K and X_𝐋 are second order vector fields, as they solve their respective Euler-Lagrange equations. Thus S(X_𝐋_K) = S(X_𝐋) =. Then necessarily S(Y_K) = 0, which means it is vertical. As by proposition <ref> both fields satisfy the symmetry property with respect to ι̃, Y_K must necessarily satisfy it too.

  *  For a general K, using Cartan's magic formula we get that _Y_Kω_K = d( _Y_Kθ_K) - ℒ_Y_Kθ_K, but as we have just shown in <ref>, Y_K is vertical, and so _Y_Kθ_K vanishes identically. Thus, _Y_Kω_K = - ℒ_Y_Kθ_K. Now, if Y_K = Y_v^i ∂_v^i + Y_V^i ∂_V^i and θ_K = θ^q_idq^i + θ^Q_idQ^i, then:

    ℒ_Y_Kθ_K = (Y_v^j ∂_v^j + Y_V^j ∂_V^j) θ^q_idq^i + (Y_v^j ∂_v^j + Y_V^j ∂_V^j) θ^Q_idQ^i


In the special case of K_F we have that ∂_V^jθ^q_i = ∂_v^jθ^Q_i= 0, which reduces the former expression to:

    ℒ_Y_K_Fθ_K_F = Y_v^j ∂_v^jθ^q_idq^i + Y_V^j ∂_V^jθ^Q_idQ^i


Furthermore, ∂_v^jθ^q_i = -1/2⟨∂^2 F/∂ v^j ∂ v^i, τ(q,Q)⟩, and similarly with ∂_V^jθ^Q_i. As τ∘ϵ̃ = 0, all these terms vanish at the identities.

  *  We know that _Y_K_Fω_𝐋 = ℒ_X_𝐋_K_Fθ_K_F - dK_F + _Y_K_Fω_K_F. We have just proven in <ref> that at the identities the last term on the right-hand side vanishes. Thus we only need to worry about the first and second terms. Proceeding as before, we can expand the Lie derivative with X_𝐋_K = v^i ∂_q^i + V^i ∂_Q^i + X_v^i ∂_v^i + X_V^i ∂_V^i:

    ℒ_X_𝐋_Kθ_K   = (v^j ∂_q^j + V^j ∂_Q^j + X_v^j ∂_v^j + X_V^j ∂_V^j) θ^q_idq^i
       + (v^j ∂_q^j + V^j ∂_Q^j + X_v^j ∂_v^j + X_V^j ∂_V^j) θ^Q_idQ^i 
       + θ^q_j∂_v^i v^j dv^i + θ^Q_j∂_V^i V^j dV^i

Again, for K = K_F we have that ∂_V^jθ^q_i = ∂_v^jθ^Q_i= 0, so simplifying this expression we get:

    ℒ_X_𝐋_K_Fθ_K_F   = (v^j ∂_q^j + V^j ∂_Q^j + X_v^j ∂_v^j) θ^q_idq^i 
       + (v^j ∂_q^j + V^j ∂_Q^j + X_V^j ∂_V^j) θ^Q_idQ^i
       + θ^q_idv^i + θ^Q_idV^i

Under the same argument as in <ref>, the terms with derivatives in the v and V variables vanish at the identities. We can also see that v^j ∂_q^jθ^q_i + V^j ∂_Q^jθ^q_i and its θ^Q counterpart must also vanish at the identities because:

    ∂_Qτ(q,q) = - ∂_qτ(q,q)

Thus .ℒ_X_𝐋_K_Fθ_K_F|_ϵ̃(TQ) = .(∂K_F/∂ v^idv^i + ∂K_F/∂ V^idV^i)|_ϵ̃(TQ) = 0. These also coincide with the v and V components of dK, so we only need to check what happens with the q and Q components.

    dK_F    = 1/2(⟨ F(Q,V), ∂τ/∂ q^i(Q,q)⟩ - ⟨∂ F/∂ q^i(q,v), τ(q,Q)⟩.
       - .⟨ F(q,v), ∂τ/∂ q^i(q,Q)⟩) dq^i
       + 1/2(⟨∂ F/∂ Q^i(Q,V), τ(Q,q)⟩ + ⟨ F(Q,V), ∂τ/∂ Q^i(Q,q)⟩.
       - .⟨ F(q,v), ∂τ/∂ Q^i(q,Q)⟩) dQ^i + ...


At the identities all terms with a bare τ vanish, and using the properties of its derivatives, the remaining terms add up together forming F_i(q,v) dq^i - F_i(q,v) dQ^i.

  *  Let us develop the left-hand side of <ref>:

    _Y_K_Fω_𝐋 = -∂^2 L/∂ v^j ∂ v^i Y_v^j d q^i + ∂^2 L/∂ V^j ∂ V^i Y_V^j d Q^i

Restricting to ϵ̃(TQ) and using <ref> we get that:

    .Y_K_F|_ϵ̃(TQ) = -W^ij F_j ∂_v^j - W^i j F_j ∂_V^j = ϵ̃_*(- Z^v_F)

where W^i j are the entries of the inverse of the Hessian matrix of L, as defined in section <ref>.



Finding the integral curves of X_E_L+Z^v_F is equivalent to solving the forced Euler-Lagrange equations

    d/dt(∂ L/∂q̇^i)- ∂ L/∂q^i=F_i(q, q̇^i)  .

Then from theorem <ref> we have that this is also equivalent to solving the unforced Lagrangian system derived by duplication given by 𝐋_K_F: TQ× TQ→ℝ and restricting the dynamics to ϵ̃(TQ); in other words, it is equivalent to solving this system's Euler-Lagrange equations

    d/dt(∂𝐋_K_F/∂q̇^i)- ∂𝐋_K_F/∂q^i   =   0
    d/dt(∂𝐋_K_F/∂Q̇^i)- ∂𝐋_K_F/∂Q^i   =   0

when restricted to ϵ̃(TQ)={(v_q, v_q)∈ TQ× TQ}. 

After theorem <ref> the following result does not add much more, but gives us a better picture of the difference between the Hamiltonian and Lagrangian side and why the dynamics at the identities coincide:


Let (L, F) and (H, F^H) be a regular forced Lagrangian system and its associated forced Hamiltonian system, and denote by 𝐇_K_F = E_𝐋_K_F∘(𝔽𝐋^×_K_F)^-1 and 𝐇_K_F the corresponding generalized Hamiltonians. Then their respective Hamiltonian vector fields X_𝐇_K_F and X_𝐇_K_F satisfy that .X_𝐇_K_F|_ϵ(T^*Q) = .X_𝐇_K_F|_ϵ(T^*Q).



Working on the Lagrangian side we know that 𝐇_K_F∘𝔽𝐋^× = E_𝐋 + K_F, where 𝔽𝐋^× is the Legendre transformed induced by the free Lagrangian. This means that we have the following two concurrent dynamics:

    _X_𝐋_K_Fω_𝐋_K_F   = d E_𝐋_K_F
    _X_𝐋_K_Fω_𝐋   = d( E_𝐋 + K_F)


The respective transformed versions of the vector fields X_𝐋_K_F and X_𝐋_K_F are X_𝐇_K_F and X_𝐇_K_F. Clearly both X_𝐋_K_F and X_𝐋_K_F can be decomposed into X_𝐋 - Y_K_F and X_𝐋 - Y_K_F respectively. We are left with:

    _Y_K_Fω_𝐋_K_F   = ℒ_X_𝐋θ_K_F - dK_F
    _Y_K_Fω_𝐋   = - dK_F


As we saw in proposition <ref>.<ref>, ._Y_K_Fω_K_F|_ϵ̃(TQ) = 0, so we are left with:

    _Y_K_Fω_𝐋   = ℒ_X_𝐋θ_K_F - dK_F

and we saw in <ref>.<ref> .ℒ_X_𝐋θ_K_F|_ϵ̃(TQ) = 0, thus .Y_K_F|_ϵ̃(TQ) = .Y_K_F|_ϵ̃(TQ). As both 𝔽𝐋^×_K_F and 𝔽𝐋^× coincide at the identities, then so will X_𝐇_K_F and X_𝐇_K_F, proving our claim.


The two terms that must vanish at the identities for both dynamics to coincide, _Y_K_Fω_K_F and ℒ_X_𝐋θ_K_F, amount to the condition that ℒ_X_𝐋_K_Fθ_K_F vanishes at the identities. This is still true for any K such that K ∘𝔽𝐋^× = K.

















































§ VARIATIONAL ORDER FOR FORCED DISCRETE LAGRANGIAN SYSTEMS



 §.§ Introduction to discrete mechanics


We will consider Q× Q as a discrete version of TQ and therefore Q× Q× Q × Q as a discrete analogue of TQ × TQ, see <cit.>. Instead of curves on Q, the solutions are replaced by sequences of points on Q. If we fix some N∈ℕ then we use the notation 

    𝒞_d(Q)={ q_d:{ k }_k=0^N ⟶ Q }

for the set of possible solutions, which can be identified with the manifold Q×(N+1)⋯× Q.
Define a functional, the discrete action map, on the space of sequences 𝒞_d(Q) by

    S_d (q_d)=∑_k=0^N-1L_d(q_k,q_k+1),    q_d∈𝒞_d(Q) .

If we consider variations of q_d with fixed end points q_0 and q_N and extremize S_d over q_1,…,q_N-1, we obtain the discrete Euler-Lagrange equations (DEL)

    D_1 L_d (q_k,q_k+1)+D_2L_d(q_k-1,q_k)=0      k=1,…,N-1   ,

where D_1L_d(q_k-1,q_k)∈ T^*_q_k-1Q and D_2L_d(q_k-1,q_k)∈ T^*_q_kQ correspond to 
dL_d(q_k-1,q_k) under the identification T^*_(q_k-1,q_k)(Q× Q)≅ T^*_q_k-1Q× T^*_q_kQ.

If L_d is regular, that is, D_12L_d is regular, then we obtain a well defined discrete Lagrangian map

    [                     F_L_d:                       Q× Q                          ⟶                      Q × Q;                                           (q_k-1,q_k)                          ⟼ (q_k,q_k+1(q_k-1,q_k))   , ]

where q_k+1 is the unique solution of (<ref>) for the given pair (q_k-1,q_k). We can further assure that the discrete Lagrangian map is invertible so that it is possible to write q_k-1=q_k-1(q_k,q_k+1), see <cit.>.

In this setting we can define two discrete Legendre transformations

    𝔽^+L_d,𝔽^-L_d:Q× Q⟶ T^*Q  ,

since each projection is equally eligible for the base point. They can be defined as

    𝔽^+L_d(q_k-1,q_k)= (q_k,D_2L_d(q_k-1,q_k))   ,
          𝔽^-L_d(q_k-1,q_k) = (q_k-1,-D_1L_d(q_k-1,q_k))   .

It holds that (𝔽^+L_d)^*ω_Q=(𝔽^-L_d)^*ω_Q=:Ω_L_d, with local expression

    Ω_L_d(q_k-1,q_k)=∂^2 L_d/∂ q_k-1^i ∂ q_k^jdq_k-1^i ∧dq_k^j .


We can also define the evolution of the discrete system on the Hamiltonian side, F_L_d:T^*Q ⟶ T^*Q, by any of the formulas

    F_L_d=𝔽^+L_d∘ (𝔽^-L_d)^-1=𝔽^+L_d∘ F_L_d∘ (𝔽^+L_d)^-1=𝔽^-L_d∘ F_L_d∘ (𝔽^-L_d)^-1  ,

because of the commutativity of the following diagram:

    (q_k-1,q_k) [rr]^F_L_d[dr]_𝔽^+L_d       (q_k,q_k+1) [ld]_𝔽^-L_d[rd]_𝔽^+L_d[rr]^F_L_d        (q_k+1,q_k+2) [ld]_𝔽^-L_d
        (q_k,p_k) [rr]_F_L_d       (q_k+1,p_k+1)

The discrete Hamiltonian map F_L_d:(T^*Q,ω_Q) ⟶ (T^*Q,ω_Q) is symplectic. Therefore the submanifold 

    ( q_k, p_k, q_k+1, p_k+1)    =   ( q_k, 𝔽^-L_d(q_k,q_k+1), q_k+1, 𝔽^-L_d(q_k+1,q_k+2) ) 
       =   ( q_k, 𝔽^+L_d(q_k-1,q_k), q_k+1, 𝔽^+L_d(q_k,q_k+1) )

is Lagrangian in (T^*Q× T^*Q,Ω_Q), where Ω_Q:=β_T^*Q^*ω_Q-α_T^*Q^*ω_Q is a symplectic form and α_T^*Q,β_T^*Q:T^*Q× T^*Q ⟶ T^*Q denote the projections onto the first and second factor respectively.


So far we have taken as the starting point a discrete Lagrangian L_d : Q× Q ⟶ℝ. However, if we start with a continuous Lagrangian and take an appropriate discrete Lagrangian then the DEL equations become a geometric integrator for the continuous Euler-Lagrange system, known as a variational integrator. 
Hence, given a regular Lagrangian function L:TQ ⟶ℝ, we define a discrete Lagrangian L_d: Q× Q ×ℝ⟶ℝ as an approximation to the action of the continuous Lagrangian. More precisely, for a regular Lagrangian L, and appropriate h,q_0,q_1, we can define the exact discrete Lagrangian as

    L_d^e(q_0,q_1,h) = ∫_0^h L(q_0,1(t),q̇_0,1(t))dt,

where q_0,1(t) is the unique solution of the Euler-Lagrange equations for L satisfying q_0,1(0)=q_0 and q_0,1(h)=q_1, see <cit.>. 
Then for a sufficiently small h, the solutions of the DEL equations for L_d^e lie on the solutions of the Euler-Lagrange equations for L, see <cit.>.

In practice, L_d^e(q_0,q_1,h) will not be explicitly given. Therefore we will take 

    L_d(q_0,q_1,h) ≈ L_d^e(q_0,q_1,h)  ,

using some quadrature rule. We obtain symplectic integrators in this way, see <cit.>.


Now we recall the result of <cit.> and 
<cit.> for a discrete Lagrangian
L_d Q× Q→ℝ.


Let L_d Q× Q→ℝ be a discrete Lagrangian. We say that
L_d is a discretization of order r if there exist an
open subset U_1⊂ TQ with compact closure and
constants C_1>0, h_1>0 so that


    ‖ L_d(q(0),q(h))-L_d^e(q(0),q(h))‖≤ C_1h^r+1
 for all solutions q(t) of the second-order Euler–Lagrange equations with initial conditions (q_0,q̇_0)∈ U_1 and for all h≤ h_1.


Following <cit.>, we have the next result about the order of a variational integrator.


If F_L_d is the evolution map of an order r
discretization L_d Q× Q→ℝ of the exact discrete
Lagrangian L_d^e Q× Q→ℝ, then

    F_L_d=F_L_d^e+𝒪(h^r+1).

In other words, F_L_d gives an integrator of order
r for F_L_d^e=F_H^h.


Note that given a discrete Lagrangian L_d Q× Q→ its
order can be calculated by expanding the expressions for
L_d(q(0),q(h)) in a Taylor series in h
and comparing this to the same expansions for the exact Lagrangian.
If the series agree up to r terms, then the discrete Lagrangian is
of order r.




§ DISCRETE LAGRANGIAN DYNAMICS OBTAINED BY DUPLICATION

We have a regular system defined by 𝐋_K: TQ× TQ→ℝ, now we consider a discretization of this Lagrangian 

    𝐋^d_K: Q× Q× Q× Q→ℝ

such that 𝐋^d_K=-𝐋^d_K∘ι̃_d where ι̃_d:  Q× Q× Q× Q → Q× Q× Q× Q is the inversion defined by

    ι̃_d(q_k, q_k+1, Q_k, Q_k+1)=(Q_k, Q_k+1, q_k, q_k+1)  .


Additionally define the identity map ϵ̃_d: Q× Q→ Q× Q× Q× Q by

    ϵ̃_d (q_k, q_k+1)=(q_k, q_k+1, q_k, q_k+1)  .



The flow F_𝐋^d_K:Q× Q× Q× Q→ Q× Q× Q× Q defined by a discrete Lagrangian 𝐋^d_K: Q× Q× Q× Q→ℝ verifying that 𝐋^d_K=-𝐋^d_K∘ι̃_d restricts to ϵ̃_d(Q× Q), that is, 

    F_𝐋^d_K∘ϵ̃_d(Q× Q)∈ϵ̃_d(Q× Q)  .



The proof is a consequence of Proposition <ref>.






§ MAIN RESULT. VARIATIONAL ORDER FOR FORCED LAGRANGIAN SYSTEMS


Now, we are in a position to state the main result of this paper. 

Let (L,F) be a forced Lagrangian system. Derive from it the extended regular Lagrangian 𝐋_K_F: TQ× TQ→ℝ and consider an order r discretization 𝐋^d_K_F: TQ× TQ→ℝ of the exact Lagrangian 

    𝐋^e_K_F (q_0, Q_0, q_1, Q_1)=∫_0^h𝐋_K_F (q_0,1(t), q̇_0,1(t), Q_0,1(t), Q̇_0,1(t))  dt

where t→ (q_0,1(t), q̇_0,1(t), Q_0,1(t)) the unique solution of the Euler-Lagrange equations for 𝐋_K_F satisfying 
q_0,1(0)=q_0, Q_0,1(0)=Q_0, q_0,1(h)=q_1, Q_0,1(h)=Q_1 and satisfying additionally that 𝐋_K_F^d∘ĩ_d=-𝐋_K_F^d. 
Then, the discrete Euler-Lagrange equations of 𝐋^d_K_F restricted to ϵ̃_d (Q× Q) give us a numerical integrator of order r for the flow of the forced Lagrangian system (L,F). 







As an example consider a Lagrangian L: ℝ^2n→ℝ: 

    L(q, v)=1/2 v^T M v-1/2q^T K q
    ; ,


and a dissipation force F(q^i, v^i)=(q^i, -D_i j v^j)

The forced Euler-Lagrange equations are: 

    Mq̈(t)+Dq̇(t)+Kq(t)=0  .


We will derive the extended Lagrangian 𝐋_K_F: ℝ^4n→ℝ. For that, we consider the function: 

    K_F(q, v, Q, V)= - 1/2 D V· (q-Q) + 1/2 D v· (Q-q) = [D/2(V+v)]^T(Q-q)

Then

    𝐋_K_F (q, v, Q, V)    = L(Q, V) - L(q, v)  - K_F(q, v, Q, V)
       =1/2 V^T M V-1/2Q^T K Q-1/2 v^T M v+1/2q^T K q - [D/2(V+v)]^T(Q-q)

By construction 𝐋_K_F(q, v, Q, V) = -𝐋_K_F(Q, V, q, v). 

Let us discretize by using the so-called “midpoint rule'"

    q≈q_0+q_1/2  , q̇≈q_1-q_0/h

which leads to:

    𝐋_K_F^d    =   h/2(Q_k+1-Q_k/h)^T M (Q_k+1-Q_k/h)
    -h/2(Q_k+Q_k+1/2)^T K (Q_k+Q_k+1/2)
          -h/2(q_k+1-q_k/h)^T M (q_k+1-q_k/h)
    +h/2(q_k+q_k+1/2)^T K (q_k+q_k+1/2)
           -h/2[D (Q_k+1-Q_k/h+q_k+1 - q_k/h)]^T (Q_k+Q_k+1/2-q_k+q_k+1/2)

Observe that 

    𝐋_K_F^d(q_k, Q_k, q_k+1, Q_k+1)=-𝐋_K_F^d(Q_k, q_k, Q_k+1, q_k+1)

Therefore, from theorem <ref> this leads to a second order method restricting the discrete Euler-Lagrange equations to 
ϵ̃_d(Q× Q). The resulting equations are not very surprising:

    M (Q_k+2-2Q_k+1+Q_k/h^2) + D(Q_k+2-Q_k/2h) + K (Q_k+2+2Q_k+1+Q_k/4) =0  .




The following results provide a purely variational base for the exact discrete forcing offered by Marsden and West, and show that usual Runge-Kutta type discretization schemes provide the same results as in their article.



The exact discrete Lagrangian defined by 𝐋^e_K_F(u_q, v_q') at the identities is equivalent to two copies of the one defined in <cit.>.



The corresponding parts for L need not be checked as they correspond trivially to those of eq.(3.2.7a) with the adequate change of notation. It remains to show that K_F generates the exact discrete forces f_d^e+, f_d^e-.

The contribution of some K to the exact discrete Lagrangian is:

    K_d^e(q_0,Q_0,q_1,Q_1) = ∫_0^h K(q(t),v(t),Q(t),V(t))dt

where t→ (q(t), Q(t))∈ Q× Q is the unique solution for 𝐋_K_F with boundary conditions q(0)=q_0, Q(0)=Q_0, q(h)=q_1, Q(h)=Q_1. 


In the case where K = K_F, differentiating K_d with respect to q_0 we get:

    D_1 K_F,d^e(q_0,Q_0,q_1,Q_1)
       = ∫_0^h [D_1 K_F·∂ q(t)/∂ q_0 + D_2 K_F·∂ v(t)/∂ q_0 + D_3 K_F·∂ Q(t)/∂ q_0 + D_4 K_F·∂ V(t)/∂ q_0]dt

where:

    D_1 K_F   = 1/2[⟨ F(Q,V), D_2 τ(Q,q)⟩ - ⟨ D_1 F(q,v), τ(q,Q)⟩ - ⟨ F(q,v), D_1 τ(q,Q)⟩]
    
    D_2 K_F   = - 1/2⟨ D_2 F(q,v), τ(q,Q)⟩
    
    D_3 K_F   = 1/2[⟨ D_1 F(Q,V), τ(Q,q)⟩ + ⟨ F(Q,V), D_1 τ(Q,q)⟩ - ⟨ F(q,v), D_2 τ(q,Q)⟩]
    
    D_4 K_F   = 1/2⟨ D_2 F(Q,V), τ(Q,q)⟩


Similar expressions are found after differentiation with respect to q_0', q_1 and q_1'. Now, when restricted to the identities, we find that:

    ϵ̃^* D_1 K_F   = ⟨ F(q, v), id_T_q Q⟩
    ϵ̃^* D_2 K_F   = 0
    ϵ̃^* D_3 K_F   = -⟨ F(q, v), id_T_q Q⟩
    ϵ̃^* D_4 K_F   = 0

where we used the fact that τ(q,q) = 0_q and ϵ̃^* D_2 τ = - ϵ̃^* D_1 τ = - id_T_q Q. This leads to:

    ϵ̃_d^ * D_1 K_F,d^e    = - ϵ̃_d^ * D_2 K_F,d^e = ∫_0^h F(q(t),u(t)) ·∂ q(t)/∂ q_0dt = f_d^e-
    ϵ̃_d^ * D_3 K_F,d^e    = - ϵ̃_d^ * D_4 K_F,d^e = ∫_0^h F(q(t),u(t)) ·∂ q(t)/∂ q_1dt = f_d^e+


Putting everything together we find two copies of the forced discrete equations with opposite sign, which is what we set to prove.






Let γ^i : ℝ× Q × Q → TQ, i = 1, ..., s, be differentiable discretisation functions, and let us use for convenience the notation τ_Q ∘γ^i(t,q_0, q_1) = (f^i(t,q_0,q_1), g^i(t,q_0,q_1)). Let also (b_i, c_i) be some quadrature coefficients such that the conservative discrete Lagrangian is approximated as:

    L_d(q_0,q_1) = h ∑_i=1^s b_i L ∘γ^i(c_i,q_0,q_1).

Then the contribution of K_F, as defined in Proposition <ref>, to the discrete Lagrangian 𝐋_K_F at the identities becomes:

    ϵ̃_d^ *∂K_F^d/∂ q_0   = -h ∑_i=1^s b_i ⟨ F ∘γ^i(c_i,q_0,q_1), ∂ f^i(c_i,q_0,q_1)/∂ q_0⟩
    ϵ̃_d^ *∂K_F^d/∂ Q_0   = h ∑_i=1^s b_i ⟨ F ∘γ^i(c_i,q_0,q_1), ∂ f^i(c_i,q_0,q_1)/∂ q_0⟩
    ϵ̃_d^ *∂K_F^d/∂ q_1   = -h ∑_i=1^s b_i ⟨ F ∘γ^i(c_i,q_0,q_1), ∂ f^i(c_i,q_0,q_1)/∂ q_1⟩
    ϵ̃_d^ *∂K_F^d/∂ Q_1   = h ∑_i=1^s b_i ⟨ F ∘γ^i(c_i,q_0,q_1), ∂ f^i(c_i,q_0,q_1)/∂ q_1⟩




For the contribution of K_F to the discrete Lagrangian we have:

    K_F^d(q_0,Q_0,q_1,Q_1) = h ∑_i=1^s b_i K(γ^i(c_i,q_0,q_1),γ^i(c_i,Q_0,Q_1))


Differentiating with respect to q_0, q_0', q_1 and q_1' we have:

    ∂K_F^d/∂ q_0   = h ∑_i=1^s b_i [ D_1 K ·∂ f^i(c_i,q_0,q_1)/∂ q_0 + D_2 K ·∂ g^i(c_i,q_0,q_1)/∂ q_0]
    ∂K_F^d/∂ Q_0   = h ∑_i=1^s b_i [ D_3 K ·∂ f^i(c_i,Q_0,Q_1)/∂ Q_0 + D_4 K ·∂ g^i(c_i,Q_0,Q_1)/∂ Q_0]
    ∂K_F^d/∂ q_1   = h ∑_i=1^s b_i [ D_1 K ·∂ f^i(c_i,q_0,q_1)/∂ q_1 + D_2 K ·∂ g^i(c_i,q_0,q_1)/∂ q_1]
    ∂K_F^d/∂ Q_1   = h ∑_i=1^s b_i [ D_3 K ·∂ f^i(c_i,Q_0,Q_1)/∂ Q_1 + D_4 K ·∂ g^i(c_i,Q_0,Q_1)/∂ Q_1]

where D_i K are the same as those of proposition <ref> with q^i(t) = f^i(t,q_0,q_1), v^i(t) = g^i(t,q_0,q_1),  (Q)^i(t) = f^i(t,Q_0,Q_1), V^i(t) = g^i(t,Q_0,Q_1). Restriction to the identities proves our claim.


 Let us choose our discretisation to be:

    L_d^α(q_0,q_1) = h L((1-α)q_0 + α q_1, q_1 - q_0/h)

as in <cit.>. This results in:

    ϵ̃_d^ *∂K_F,d^α/∂ q_0   = h (1 - α) F((1-α)q_0 + α q_1, q_1 - q_0/h)
    ϵ̃_d^ *∂K_F,d^α/∂ Q_0   = - h (1 - α) F((1-α)q_0 + α q_1, q_1 - q_0/h)
    ϵ̃_d^ *∂K_F,d^α/∂ q_1   = h α F((1-α)q_0 + α q_1, q_1 - q_0/h)
    ϵ̃_d^ *∂K_F,d^α/∂ Q_1   = - h α F((1-α)q_0 + α q_1, q_1 - q_0/h)

which coincides with their result.




 §.§ Numerical tests

For our numerical tests we have chosen a well-known system composed of two coupled van der Pol oscillators (cf. <cit.>). Remember that a single dimensionless van der Pol oscillator is described by the differential equation:

    q̈ - (ϵ - q^2) q̇ + q = 0

where ϵ is a parameter related to the damping of the system.

The dimensionless system we are going to study can be thought to be composed of two coupled harmonic oscillators with slightly differing natural frequencies under the action of non-linear forcing. Its configuration manifold is 𝕋×𝕋 = 𝕋^2, with velocity phase space T𝕋^2 where we will use local coordinates (q_1, q_2, v_1, v_2), and the Lagrangian describing the non-forced part L: T𝕋^2 →ℝ is:

    L = 1/2(v_1^2 + v_2^2) - 1/2[q_1^2 + (1 + ρ)q_2^2] - λ(q_1 - q_2)^2

where ρ accounts for the deviation of q_2 from the natural frequency of q_1, and λ measures the intensity of the coupling between both oscillators. The van der Pol force acting on this system is F = (ϵ - q_1^2) v_1 d q_1 + (ϵ - q_2^2) v_2 d q_2. As our configuration space is flat, τ(q,Q) = Q - q, and the generalized potential K is:

    K_F = 1/2∑_i = 1^2 [(ϵ - q_i^2) v_i + (ϵ - Q_i^2) V_i] (q_i-Q_i)

Note that for such an L and K_F, at the identities we have that v_i = p_i, i = 1,2, so they are interchangeable.







We chose to discretize the corresponding generalized Lagrangian, 𝐋_K_F, using Lobatto schemes of 2, 3, 4 and 5 stages. The order of an s-stage Lobatto method is p = 2 s - 2 so the resulting numerical methods are of order 2, 4, 6 and 8 respectively. The parameters used for the numerical simulations shown here are (ϵ,ρ,λ) = (0.5, 0.02, 0.8), for no particular reason. The other choices of parameters that were tested showed essentially the same behaviour. We run each simulation for a total of 1 unit of simulation time with several different choices of step-size h ranging between 5 · 10^-5 and 1 and measure numerical error as the difference between the final value of the magnitude in study found for a reference simulation and the corresponding one for that we want to study. In this case our reference is taken as the simulation with the finest step-size. The initial values chosen for the results on diagrams <ref> and <ref> are (q_1, q_2, v_1, v_2) = (-1/2, -1/4, 0, 4). The results shown on diagram <ref> were found as the average from an ensemble of 25 random initial values in the square [-4,4]×[-4,4] and the pointed lines represent the maximum and minimum values found in said ensemble.

For the resolution of the resulting non-linear system of equations derived for each method, we used MATLAB's  with  and , which explains the flat tails.

Diagram <ref> is a composite plot showing the error in the energy for the different Lobatto methods tested. The results are in agreement with the result of theorem <ref>. We have chosen to show only the energy to avoid clutter, but the same holds for each of the components of the system.



§ GEOMETRY OF THE METHOD OF DUPLICATION OF VARIABLES

In this section, we analyze the interesting geometry related with the proposed method of duplication of variables that is related with some results about symplectic groupoids. Additionally, this section will allow us in the future extend our results for reduced systems, and in general for Hamiltonian systems defined on Lie algebroids. 
We will review the definition of Lie groupoid and its associated Lie algebroid and then we will introduce the notion of symplectic Lie algebroid.









 §.§ Lie groupoids and algebroids


First of all, we will recall some definitions related with Lie groupoid and Lie algebroids.
 (for more details, see <cit.>).



A groupoid over a set Q is a set G together with the
following structural maps:


  * A pair of maps α: G → Q, the source, and
β: G → Q, the target. Thus, we can think an element g ∈ G an arrow from x= α(g) to y = β(g) in Q

    *=0∙@/^1pc/@<1ex>[rrr]_g         *=0∙

The source and target mappings define the set of composable pairs

    G_2={(g,h) ∈ G × G / β(g)=α(h)}.


  * A multiplication on composable elements μ: G_2→ G, denoted simply by μ(g,h)=gh, such that


  * α(gh)=α(g) and β(gh)=β(h).

  * g(hk)=(gh)k.

If g is an arrow from x = α(g) to y = β(g) =
α(h) and h is an arrow from y to z = β(h) then
gh is the composite arrow from x to z

    *=0∙@/^2pc/@<2ex>[rrrrrr]_gh@/^1pc/@<2ex>[rrr]_g         *=0∙@/^1pc/@<2ex>[rrr]_h         *=0∙


  * An identity section ϵ: Q → G of α and β, such that


  * ϵ(α(g))g=g and gϵ(β(g))=g.


  * An inversion map ι: G → G, to be denoted simply by ι(g)=g^-1, such that


  * g^-1g=ϵ(β(g)) and gg^-1=ϵ(α(g)).


    *=0∙@/^1pc/@<2ex>[rrr]_g         *=0∙@/^1pc/@<2ex>[lll]_g^-1



A groupoid G over a set Q will be denoted simply by the symbol
G ⇉ Q.

The groupoid G ⇉ Q is said to be a Lie
groupoid if G and Q are differentiable manifolds and all the structural maps
are differentiable with α and β differentiable
submersions. If G ⇉ M is a Lie groupoid then μ
is a submersion, ϵ is an immersion and ι is a
diffeomorphism. Moreover, if x ∈ M, α^-1(x) (resp.,
β^-1(x)) will be said the α-fiber (resp.,
the β-fiber) of x.



Typical examples of Lie groupoids are: the pair or banal groupoid Q × Q over Q (the example that we have used along all this paper), a Lie group G (as a Lie groupoid over a single point), the Atiyah groupoid (Q × Q)/G (over Q/G) associated with a free and proper action of a Lie group G on Q , etc.



 If G ⇉ Q is a Lie groupoid and g ∈ G then the left-translation by
g ∈ G and the right-translation by g are the
diffeomorphisms

    [ l_g: α^-1(β(g)) ⟶α^-1(α(g))                           ;            h ⟶
    l_g(h) = gh,; r_g: β^-1(α(g)) ⟶β^-1(β(g))                           ;            h ⟶ r_g(h) = hg. ]


Note that l_g^-1 = l_g^-1 and r_g^-1 = r_g^-1.


A vector field ξ∈𝔛(G) is said to be
left-invariant (resp., right-invariant) if it is
tangent to the fibers of α (resp., β) and
ξ(gh) = (T_hl_g)(ξ_h) (resp.,
ξ(gh) = (T_gr_h)(ξ(g))), for (g,h) ∈
G_2.


The infinitesimal version of a Lie groupoid is a Lie algebroid which is defined as follows.


A Lie algebroid is a real vector bundle A→ Q equipped with a Lie bracket 
· ,· on its sections Γ(A) and a bundle map ρ: A→ TQ called the anchor map such that 
the homomorphism of C^∞(Q)-modules induced by the anchor map, that we also denote by ρ: Γ(A)→𝔛(Q), verifies 

    X, f Y =f X, Y +ρ(X)(f) Y,

for X, Y∈Γ (A) and f∈ C^∞(Q).


With this definition the anchor map ρ: Γ(A)→𝔛(Q)
is a Lie algebra homomorphism, where 𝔛(Q) is endowed with the usual Lie bracket of vector field [·, ·]. 


Given a Lie groupoid G ⇉ Q, the associated Lie algebroid AG→ Q is given by its fibers 
 A_qG = V_ϵ(q)α = Ker
(T_ϵ(q)α).
There is a bijection between the space Γ (AG) and the set of
left-invariant vector fields on G. If
X is a section of τ: AG → Q, the corresponding
left-invariant vector field on G will
be denoted X (resp., X), where

    X(g) = (T_ϵ(β(g))l_g)(X(β(g))),

for g ∈ G. Using the above facts, one may introduce a 
bracket · , · on the space of sections Γ(AG) and a bundle map ρ: AG → TQ, which
are defined by

    X, Y = [X, Y], [.3cm]ρ(X)(q) = (T_ϵ(x)β)(X(q)),

for X, Y ∈Γ(AG) and q ∈ Q. 

Using that [·, ·] induces a Lie algebra structure on the space of vector fields on G, it is easy to prove that ·, · also defines a Lie algebra structure on Γ(AG). In addition, it follows that

    X, f Y = f  X, Y + ρ(X)(f) Y,

for X, Y ∈Γ(AG) and f ∈ C^∞(Q).

One can also stablish a bijection between sections X∈Γ(AG) and right invariant vector fields X∈𝔛(G) defined by 

    X(g) = -(T_ϵ(α(g))r_g)((T_ϵ
    (α(g))ι)( X(α(g)))),

which yields the Lie bracket relation

    X, Y =- [X, Y]  .



The following proposition will be useful for the results in this paper. 


Let G ⇉ Q be a Lie groupoid and Z∈𝔛(G) a vector field invariant by the inversion, that is, 

    T_g ι(Z(g))=Z(g^-1),  for all  g∈ G  .
 
Then, for all q∈ Q, 

    Z(ϵ(q))∈ T_ϵ(q)ϵ (Q)  .



For all v_q∈ A_qG consider an α-vertical curve g: I→ G such that v=dg/dt(0). Then

    T_(ϵ(q), ϵ(q))μ(0_q, v_q)=d/dtμ(ϵ(q), g(t))|_t=0=dg/dt(0)=v

Also, for the β-vertical curve g^-1: I→ G we have 

    T_(ϵ(q), ϵ(q))μ(T_ϵ(q)ι(v), 0_q)=d/dtμ(ι(g(t)), ϵ(q))|_t=0=T_ϵ(q)ι(v)  .

Therefore, 

    T_(ϵ(q), ϵ(q))μ(T_ϵ(q)ι(v), v)=v+T_ϵ(q)ι(v)  .

Since μ(g^-1(t), g(t))=ϵ(β(g(t))), then

    (T_ϵ(q)ι)(v)=-v+T_ϵ(q)(ϵ∘β)(v)  .

Using that

    Z(ϵ(q))-T_ϵ(q)(ϵ∘α)(Z(ϵ(q)))∈ A_qG  ,

then from expression (<ref>):

    T_ϵ(q)ι(Z(ϵ(q)))-T_ϵ(q)(ϵ∘α)(Z(ϵ(q)))+Z(ϵ(q))∈ T_ϵ(q)ϵ(Q)  .

but from the hypothesis about Z, we have that

    T_ϵ(q)ι(Z(ϵ(q)))=Z(ϵ(q))  .

Therefore

    Z(ϵ(q))∈ T_ϵ(q)ϵ(Q)  .




 §.§ Symplectic groupoid



A symplectic groupoid is a Lie groupoid G ⇉ Q, such that


  * (G, ω) is a symplectic manifold,

  * the graph of μ: G_2→ G is a Lagrangian submanifold of G^-× G^-× G, where G^-=(G, -ω) has the negative symplectic structure.



If G ⇉ Q is a symplectic groupoid with symplectic form ω on G then one may prove that ( T_gα)^ω= T_gβ, for g∈ G, where 

    ( T_gα)^ω={v∈ T_gG  | ω(v, u)=0, for all   u∈ T_gα}  ,

that is, the symplectic orthogonal of T_gα. Moreover there exists a unique Poisson structure on Q such that α: G→ Q (respectively, 
β: G→ Q) is a Poisson (respectively, anti-Poisson) morphism. Moreover, the inversion map is an antisymplectomorphism, that is, 
ι^*ω=-ω  .


Let G ⇉ Q be a Lie groupoid, an let A^*G→ Q be the dual vector bundle of the associated Lie algebroid AG. Then, the cotangent groupoid T^*G ⇉ A^*G is a symplectic groupoid with the canonical symplectic form ω_G. Given μ∈ T_g^*G, the source and target mappings are defined 

    ⟨α̃(μ), X(α(g))⟩=⟨μ, X(g)⟩ ,   ⟨β̃(μ), X(β(g))⟩=⟨μ, X(g)⟩  .

for all X∈Γ(AG).
(See <cit.> for more details and the definition of the remaining structure maps of this Lie groupoid).





Let G ⇉ Q be a symplectic groupoid with symplectic form ω_G symplectic groupoid and E: G→ℝ a function such that E∘ι=-E. Then, the corresponding Hamiltonian vector field X_E

    _X_Eω=dE  ,

verifies that X_E(ϵ(q))∈ T_ϵ(q)ϵ (Q) for all q∈ Q.


Since ι^*ω=-ω then for all Y∈𝔛(G)

    ⟨dE, Y⟩=ω(X_E, Y)=-ι^*ω(X_E, Y)=-ω( ι_*X_E, ι_*(Y))  .

but form the hypothesis we have that

    ⟨dE, Y⟩=-⟨d(E∘ι), Y⟩=-⟨dE, ι_*(Y)⟩=-ω(X_E, ι_*(Y))  .

Therefore, from Proposition <ref> we deduce that X_E (ϵ (q))∈ T_ϵ(q)ϵ(Q).




§ CONCLUSIONS AND FUTURE WORK


The main contributions of this paper are:


  * Using the duplication of variables we have rigorously deduced the error analysis of forced lagrangian systems in terms of variational error.

  * With this technique it is possible to design efficient numerical methods for forced lagrangian systems  using previous results for variational integrators including high-order methods.
 
  * We have completely elucidated the geometry of  the procedure of duplication of variables connecting with the concept of symplectic groupoid.
 
  * Moreover, we have separately study the hamiltonian and lagrangian formalism and stablished the relation between both. 
 
 
 Future work includes:


  * Let G ⇉ M be a Lie groupoid with structural maps

    α, β: G → M,   ϵ: M → G,   ι: G → G,
      μ: G_2→ G.

Suppose that τ: AG → M is the Lie algebroid of G and that
𝒫^τG is the prolongation of G over the
fibration τ: AG → M, that is,

    𝒫^τG = AG _τ-3pt×_α
    G _β-3pt×_τ AG.

It is clear that 𝒫^τG is equipped with a  Lie groupoid structure over AG but also the vector bundle π^τ:𝒫^τG → G admits an integrable Lie algebroid structure (see <cit.>). This is the corresponding version for reduced systems of the space with “duplicated variables". We will check in a future paper how apply this methodology to analyze the order 
of geometric integrators for forced systems using  this method. 

  * An example interesting will be the case of Euler-Poincaré equations and double bracket dissipation <cit.>. We will study the possibility of constructing geometric integrators preserving some of the geometric structure. For instance, it can be checked that in particular examples,  the energy is dissipated but the angular momentum is not. 




§ ACKNOWLEDGEMENTS

The authors have been partially supported by Ministerio de Economía, Industria y Competitividad (MINEICO, Spain) under grants MTM 2013-42870-P, MTM 2015-64166-C2-2P, MTM2016-76702-P and “Severo Ochoa Programme for Centres of Excellence” in R&D (SEV-2015-0554). R. Sato has been financially supported by a FPI scholarship from MINEICO, Spain. 






§ INVARIANT SETS DEFINED FROM DISCRETE SYMMETRIES OF A LAGRANGIAN FUNCTION


The Lagrangian 𝐋_K: TQ× TQ→ℝ verifies from construction that 𝐋_K∘ι̃=-𝐋_K 
where 
 ι̃: TQ× TQ→ TQ× TQ is the inversion mapping given by ι̃(u_q, v_q')=(v_q', u_q).
Then the identity set ϵ̃(TQ) is invariant set for the flow of X_E_𝐋_K as a consequence of the following Proposition. Observe that ι̃ is the tangent lift of the map 
Q× Q→ Q× Q given by (q_k, q_k+1)↦ (q_k+1, q_k).

Let L: TQ→ℝ be a regular Lagrangian and φ: Q→ Q a diffeomorphism verifying that 
L∘φ_*=± L. Denote by M_φ={v_q∈ Q  | φ_*(v_q)=v_q}. Then M_φ is an invariant set for any solution of the Euler-Lagrange equations. 


Consider the action sum 

    [                   𝒥_L:   C^2(q_0, q_1, [a,b])                      ⟶                      ℝ;                                             c                      ⟼ ∫_a^b L(c(t), ċ(t)) dt ]

where 
C^2(q_0, q_1, [a,b])={c: [a, b]→ Q  |  c  is C^2, c(a)=q_0, c(b)=q_1}.
Its tangent space is

    T_c C^2(q_0, q_1, [a,b])={ X: [0, T]→ TQ  |   X  is C^1, τ_Q(X(t))=c(t), X(a)=0, X(b)=0}

Typically, we express X∈ T_c C^2(q_0, q_1, [a,b]) as the tangent vector to a curve at s=0 in C^2(q_0, q_1, [a,b]),

    s∈ (-ϵ, ϵ)⊂ℝ⟼ c_s∈  C^2(q_0, q_1, [a,b])

with c_0=c. That is, 

    X=d c_s/ds|_s=0


We then have

    d𝒥_L(c)(X)   =   d/ds|_s=0 (𝒥_L(c_s))
       =   ∫_a^b d/ds|_s=0 L(c_s(t), ċ_s(t))  dt

Using that 
L∘φ=± L

    d𝒥_L(c)(X)   =   ±∫_0^T d/ds|_s=0 L(φ(c_s(t)), d/dt(φ(c_s(t))))  dt
       =   ±d/ds|_s=0 (𝒥_L(φ∘ c_s))
       =   ±d𝒥_L(φ∘ c)(φ_*∘ X)

Observe that φ_*∘ X∈ T_φ∘ c C^2(φ(q_0), φ(q_1), [a,b]). Since φ is a diffeomorphism then 
c is a critical point of 𝒥_L iff φ∘ c is a critical point of 𝒥_L.

Now, if c: [a, b]→ Q is a solution of the Euler-Lagrange equations (d𝒥_L(c)=0) with ċ(a)∈ M_φ then also φ∘ c is a solution of the Euler-Lagrange equations. Observe that 

    d (φ∘ c)/dt(a)=φ_*(ċ(a))=ċ(a)

Then, φ∘ c and c are solutions of the Euler-Lagrange equations with the same initial conditions. Since L is regular, it implies that c=φ∘ c. 
and ċ(t)∈ M_φ, for all t.





§ INVARIANT SETS DEFINED FROM DISCRETE SYMMETRIES OF A DISCRETE LAGRANGIAN FUNCTION


The discrete Lagrangian 𝐋^d_K: TQ× TQ→ℝ verifies that 𝐋^d_K∘ĩ_d=-𝐋^d_K . The following proposition gives the required result as a particular case. 




Let L_d: Q× Q→ℝ be a regular discrete Lagrangian and φ_d: Q→ Q a diffeomorphism verifying that 
L_d∘ (φ_d×φ_d)=± L_d. Denote by M_φ_d={(q, q')∈ Q  |  φ_d(q)=q, φ_d(q')=q'}. Then F_φ_d is an invariant set for any solution of the discrete Euler-Lagrange equations. 


The proof is similar to that of Proposition <ref>. 
Consider the space

    𝒞_d(q_0, q_N)={q_d: k_k=0^N⟶ Q  |  q_0,q_N  fixed}

and the discrete action sum 

    [                     S_d:            𝒞_d(q_0, q_N)                        ⟶                        ℝ;                                               q_d                        ⟼ ∑_k=0^N-1L_d(q_k, q_k+1) ]

The extremals are characterized as the solutions of  the discrete Euler-Lagrange equations:

    D_1L_d(q_k, q_k+1)+D_2 L_d(q_k-1, q_k)=0,  k=1, …, N-1  .

Then, it is clear that if {q_k}_k=0, …,  N is a solution of the discrete Euler-Lagrange equations, then from the invariance of L_d we easily derive that 
{φ_d(q_k)}_k=0, …,  N is also a solution with boundary conditions φ_d(q_0) and φ_d(q_N). 

Therefore, if L_d is regular we have defined its discrete flow  or discrete Lagrangian map: 

    [          F_L_d:            Q× Q               ⟶           Q × Q;                     (q_k-1,q_k)               ⟼ (q_k,q_k+1)   , ]

observe that also F_L_d(φ_d(q_k-1), φ_d(q_k))=(φ_d(q_k), φ_d(q_k+1).
Now starting from initial conditions (q_0, q_1)∈ M_φ_d, that is, φ_d(q_0)=q_1, φ_d(q_1)=q_1 from the unicity of solutions of the discrete Euler-Lagrange equations we obtain that  (q_k-1, q_k)∈ M_φ_d, k=1, …, N and, as a consequence, M_φ_d is an invariant set of the discrete Euler-Lagrange equations. 







 






