Retraction Maps: A Seed of Geometric Integrators

Abstract

The classical notion of retraction map used to approximate geodesics is extended and rigorously defined to become a powerful tool to construct geometric integrators and it is called discretization map. Using the geometry of the tangent and cotangent bundles, we are able to tangently and cotangent lift such a map so that these lifts inherit the same properties as the original one and they continue to be discretization maps. In particular, the cotangent lift of a discretization map is a natural symplectomorphism, what plays a key role for constructing geometric integrators and symplectic methods. As a result, a wide range of (higher-order) numerical methods are recovered and canonically constructed by using different discretization maps, as well as some operations with Lagrangian submanifolds.

1 Introduction

The notion of retraction map is an essential tool in different research areas like optimization theory, numerical analysis and interpolation (see [3] and references therein).

In optimization theory, the goal is to find a value x in a differentiable manifold M such that f(x) is the minimum of a real-valued function $f: M\rightarrow {{\mathbb {R}}}$. In the case that M is a linear space, as ${{\mathbb {R}}}^n$ equipped with the standard inner product, the notions of gradient or Hessian of the function f are properly defined and give us useful local information to localize the possible candidates to minimize f. Moreover, gradient descent or Newton’s method can also be used to search for a solution.

Riemannian geometry allows us to introduce similar concepts to gradient and Hessian in a differentiable manifold paving the way for optimization. But we need another important ingredient: how to move on a manifold. In Riemannian geometry this notion is given by the exponential map. On a Riemannian manifold (M, g) (or more generally a semi-Riemannian manifold), we can define the Riemannian exponential map ${\mathrm{exp}}_x: T_x M \rightarrow M$ at the point x. As mentioned, for instance, in [22], ${\mathrm{exp}}_x(\xi ) = \sigma (1),$ for $\xi \in T_x M$, where $\sigma : [0, 1] \rightarrow M$ is the unique geodesic in M with initial velocity $\xi$, that is, $\sigma (0) = x$ and ${\dot{\sigma }}(0) =\xi$. Moreover, there exist open subsets ${{\mathcal {U}}} \subseteq T_x M$ and $U \subseteq M$, with ${{\mathcal {U}}}$ star-shaped about $0_x \in {{\mathcal {U}}}$ and $x \in U$, such that ${\mathrm{exp}}_x: {{\mathcal {U}}} \rightarrow U$ is a diffeomorphism and ${\mathrm{exp}}_x(t\xi ) = \sigma (t)$ and $T_{0_x}\mathrm{exp}_x = Id_{T_x M}$. However, only in simple examples it is possible to explicitly compute the exponential map of a Riemannian manifold. Therefore, efficient approximations of geodesics are crucial for designing algorithms on manifolds. Here is where retraction maps play an important role (see [2, 3] and references therein).

Roughly speaking, retraction maps provide a way to select a smooth curve on a differentiable manifold given an initial position and velocity. Such a curve is an approximation of the Riemannian exponential map. More specifically, a retraction is typically defined as a local $C^1$-map $R_x: U_x\subset T_xM\rightarrow M$ such that $R_x(0_x) = x$ and $\left. \frac{d}{dt}\right| _{t=0}R_x(t\xi )=T_{0_x}R_x(\xi )= \xi$ for all $\xi \in T_x M$, where we use the identification $T_{0_x}T_x M \equiv T_xM$ and $T_{0_x}R_x: T_{0_x}(T_xM)\rightarrow T_xM$ denotes the tangent map of $R_x$ at $0_x$ (see [1]). Observe that, since we are using first-order approximations, this definition is independent of the initial Riemannian metric . However, for second- or higher-order retractions the particular Riemannian metric does play a role. The property $\left. \frac{d}{dt}\right| _{t=0}R_x(t\xi )=\xi$ implies that $d_g(R_x(t\xi ),\sigma (t))=O(t^2)$, where $d_g$ denotes the Riemannian distance (see [48]).

For our purposes we will need a more general definition of a retraction map. We construct a discretization map $R_d: U\subset TM\rightarrow M\times M$ in Definition 2.2, where the image of $\xi \in U$ is now two “nearby" points of M. We understand such a map as a discretization of the tangent bundle because TM is locally diffeomorphic to two copies of the manifold M. As they can be related to retraction maps, they are denoted by $R_d$, where the subscript d stands for discretization.

As an example, if we have a Riemannian manifold (M, g), with associated exponential map exp, then a discretization map is

$$\begin{aligned} R_d(\xi )=\left( exp_{\tau _M(\xi )} \left( -\frac{1}{2}\xi \right) , exp_{\tau _M(\xi )}\left( \frac{1}{2}\xi \right) \right) , \end{aligned}$$

where $\tau _M: TM\rightarrow M$ is the canonical projection of the tangent bundle. This particular map $R_d$ applied in Eq. (1) will lead to the implicit midpoint method on Euclidean spaces. We precisely discuss the properties of these discretization maps in Sect. 2.

In numerical analysis, if we have a vector field X on M, that is, a section $X: M\rightarrow TM$ (that is, $\tau _M\circ X=Id_{M}$), and we want to find a numerical approximation of the integral curves, an idea is to use a discretization map and consider the following first-order discrete equation:

$$\begin{aligned} h X\left( \tau _M\left( R_d^{-1}(x_k, x_{k+1})\right) \right) =R_d^{-1}\left( x_k, x_{k+1}\right) \,. \end{aligned}$$

(1)

We prove in Proposition 2.1 that $R_d$ is a local diffeomorphism and the inverse map can be computed. Given an initial condition $x_0$, we might be able to solve the implicit system (1) to find a sequence $\{x_k\}$ which is an approximation of $\{x(kh)\}$, where x(t) is the integral curve of X with initial condition $x_0$ and h is the time step. For instance, if M is the vector space ${{\mathbb {R}}}^n$ and $R_d(x, v)=\left( x-\frac{v}{2}, x+\frac{v}{2}\right)$, then Eq. (1) becomes

$$\begin{aligned} \frac{x_{k+1}-x_k}{h}=X \left( \frac{x_{k}+x_{k+1}}{2}\right) . \end{aligned}$$

Our main interest in this article consists of designing numerical methods for second-order differential equations (SODEs) and mainly for Hamilton’s equations. For instance, a second-order differential equation $\ddot{x}=f(x, {\dot{x}})$ is geometrically represented by a special vector field

$$\begin{aligned} \Gamma (x,{\dot{x}})={\dot{x}}\frac{\partial }{\partial x}+ f(x, {\dot{x}})\frac{\partial }{\partial {\dot{x}}}\; , \end{aligned}$$

which is now defined on the tangent bundle TM of M [1]. These vector fields are called SODEs.

On the other hand, it is well-known that the classical Hamilton’s equations are defined on the cotangent bundle $T^*M$ of the manifold M. Therefore, we face the problem of how, given a discretization map on M, we can lift it to the tangent and cotangent bundles. Besides, we define in Proposition 2.5 adjoint discretization maps by inversion with the objective to construct symplectic symmetric numerical methods of higher order.

In Sect. 3.1, we lift a discretization map on a manifold to the tangent bundle using the canonical involution. This tangent lift makes possible to define geometric discretizations of SODEs in Sect. 4.

In Sect. 3.2 we lift a discretization map on a manifold to the cotangent bundle using well-known constructions from symplectic geometry. We show in Sect. 3.3 that this cotangent lift is nothing else than the dual construction of the above-mentioned tangent lift. Moreover, it is essential to prove that the cotangent lift of a discretization map is always a symplectomorphism because that makes possible to construct symplectic integrators for Hamilton’s equations and Euler–Lagrange equations in Sect. 5.

In Sect. 3.4 we carefully work out a few examples of the discretization maps on different manifolds.

Section 5 describes how to obtain numerical methods for Euler–Lagrange equations and Hamilton’s equations using the tools described in the previous sections. In particular, in Sect. 5.3 we compare the geometric integrators from Sect. 5.1 with the theory of discrete variational calculus [41]. When the symplectic numerical methods in Sect. 5.1 are understood as Lagrangian submanifolds, we can start to compose Lagrangian submanifolds coming from different discretization maps to construct general symplectic methods for Hamilton’s equation in Sect. 6. We study how to define higher-order geometric methods by composing symmetric symplectic methods in Sect. 6.2.

Along the paper we show how well-known geometric methods (Newmark, Störmer-Verlet, etc.) are obtained using the new tools described here. Hence, the work developed in this paper opens the path to define, even higher-order, geometric integrators for more complex mechanical systems that may include forced systems, system with constraints, optimal control problems, Dirac systems, etc. We describe specific future research lines in Sect. 7.

2 Retraction Maps

A retraction map plays the role of generalizing the linear-search methods in Euclidean spaces to general manifolds. On a manifold with nonzero curvature to move along the tangent line does not guarantee that the motion stays on the manifold. The retraction map provides the tool to define the notion of moving in a direction of a tangent vector while staying on the manifold. That is why retraction maps have been widely used to construct numerical integrators of ordinary differential equations, since it allows us to move from a point and a velocity to one nearby point so that the differential equation can be discretized.

The first notion of retraction that appears in the literature can be found in [9] from a topological viewpoint. Later on, the notion of retraction map as defined below is used to define Newton’s method on Riemannian manifolds [4, 48].

Definition 2.1

A retraction map on a manifold M is a smooth mapping R from the tangent bundle TM onto M. Let $R_x$ denote the restriction of R to $T_xM$, the following properties are satisfied:

• 1. $R_x(0_x)=x$ where $0_x$ denotes the zero element of the vector space $T_xM$.

• 2. With the canonical identification $T_{0_x}T_xM\simeq T_xM$, $R_x$ satisfies

  $$\begin{aligned} {\mathrm{D}}R_x(0_x)=T_{0_x}R_x={\mathrm{Id}}_{T_xM}, \end{aligned}$$

  (2)

  where ${\mathrm{Id}}_{T_xM}$ denotes the identity mapping on $T_xM$.

The condition (2) is known as local rigidity condition since, given $\xi \in T_xM$, the curve $\gamma _\xi (t)=R_x(t\xi )$ has $\xi$ as tangent vector at x, i.e.,

$$\begin{aligned} {\dot{\gamma }}_\xi (t)= \langle DR_x(t\xi ), \xi \rangle \; \hbox { and, in consequence, } {\dot{\gamma }}_\xi (0)= {\mathrm{Id}}_{T_xM}(\xi )=\xi \; . \end{aligned}$$

This notion connects with the geometric interpretation of the exponential map exp on Riemannian manifolds given in [22, Chapter 3.2]. Therefore, the image of $\xi$ through the exponential map is a point on the Riemannian manifold obtained by moving along a geodesic a length equal to the norm of $\xi$ starting with the velocity $\xi /\Vert \xi \Vert$, that is,

$$\begin{aligned} exp_x(\xi )=\sigma (\Vert \xi \Vert )\; , \end{aligned}$$

where $\sigma$ is the unit speed geodesic such that $\sigma (0)=x$ and ${\dot{\sigma }}(0)=\xi /\Vert \xi \Vert$.

Remember that the exponential map is a typical example of a retraction map. With all that in mind we are able to generalize the property of local rigidity in Definition 2.1 that allows a discretization of the tangent bundle of the configuration manifold opening a new path to construct numerical integrators.

After studying the contribution given in [18, 39] we define a generalization of the retraction map in Definition 2.1. Given a point and a velocity, we obtain two nearby points that are not necessarily equal to the initial base point. As discussed in the sequel, numerical methods will be recovered from this new map.

Definition 2.2

A map $R_d:U\subset TM\rightarrow M\times M$ given by

$$\begin{aligned} R_d(x,v)=(R^1(x,v),R^2(x,v)), \end{aligned}$$

where U is an open neighborhood of the zero section $0_x$ of TM, defines a discretization map on M if it satisfies

• 1. $R_d(x,0)=(x,x)$,

• 2. $T_{0_x}R^2_x-T_{0_x}R^1_x:T_{0_x}T_xM\simeq T_xM\rightarrow T_xM$ is equal to the identity map on $T_xM$ for any x in M, where $R^a_x$ denotes the restrictions of $R^a$, $a=1,2$, to $T_xM$.

If $R^1(x,v)=x$, the two properties in Definition 2.2 guarantee that the both properties in Definition 2.1 are satisfied by $R^2$. Thus, as mentioned, Definition 2.2 generalizes Definition 2.1.

Proposition 2.1

Let $R_d$ be an discretization map on M, $R_d$ is a local diffeomorphism from some neighborhood of the zero section of TM.

Proof

Let $(x^i)$ be local coordinates for M centered at x, and $(x^i, v^i)$ be the corresponding induced coordinates on TM centered at $0_x$. By the definition of the discretization map, the Jacobian matrix of $R_d$ at $(x^i,0)$ is locally written as

$$\begin{aligned} \left( \begin{array}{cc} {\mathrm{Id}}&{} \frac{\partial R^1}{\partial v}(x, 0)\\ {\mathrm{Id}} &{} \frac{\partial R^2}{\partial v}( x, 0) \end{array} \right) \, , \end{aligned}$$

where ${\mathrm{Id}}$ denotes the identity matrix. Note that the regularity of that Jacobian matrix is equivalent to the invertibility of the following matrix

$$\begin{aligned} \left( \begin{array}{cc} {\mathrm{Id}}&{} \frac{\partial R^1}{\partial v}(x, 0)\\ 0 &{}\frac{\partial R^2}{\partial v}(x, 0)-\frac{\partial R^1}{\partial v}(x, 0) \end{array} \right) =\left( \begin{array}{cc} {\mathrm{Id}}&{} \frac{\partial R^1}{\partial v}(x, 0)\\ 0&{} {\mathrm{Id}}\end{array} \right) \, , \end{aligned}$$

due to the property 2 in Definition 2.2. Therefore, the inverse function theorem guarantees that $R_d$ is a local diffeomorphism from some neighborhood of the identity section to its image. $\square$

There is a general and interesting way to obtain discretization maps from the usual retraction maps. The following result is very useful for Examples 2.2 and 2.3.

Proposition 2.2

Let $R:TM \rightarrow M$ be a retraction map as in Definition 2.1. For any $\theta \in [0,1]$ the map $R_d:TM \rightarrow M\times M$ given by

$$\begin{aligned} R_d(x,v)=\left( R(x,-\theta v),R(x,(1-\theta )v) \right) \end{aligned}$$

is a discretization map on M.

Proof

From the definition of retraction map, it is immediate that $R_d(x,0)=(R(x,0),R(x,0))=(x,x)$ and $T_{0_x}R^2_x-T_{0_x}R^1_x=(1-\theta )\, T_{0_x}R_x+ \theta \, T_{0_x}R_x\equiv {\mathrm{Id}}_{T_xM}$. $\square$

Starting from a retraction map, we may define different discretization maps, as shown in the above proposition. In the sequel, we will see that these different maps will lead to known numerical methods. For step size h and retraction map $R^h(x,v)=x+h\,v$ on the Euclidean space, one possible discretization map is $R_d^h(x,v)=(x,x+h\, v)$ that corresponds with a first-order integrator method as described, for instance, in [42]. However, other discretization maps may be defined from the same retraction map to construct different integrators. For example, for step size h and the above retraction map $R^h$ we define:

$$\begin{aligned} R_d^h(x,v)=(R^{-h/2}(x,v), R^{h/2}(x,v))\, , \end{aligned}$$

that corresponds with a second-order method as described in [42].

Let us describe another method to generate more discretization maps from a given one that will be useful in Sect. 6.2 to obtain some higher-order numerical methods. Define the inversion map $I_M: M\times M\rightarrow M\times M$ by $I_M (x,y)=(y, x)$ for all $x, y\in M$.

Proposition 2.3

If $R_d:U\subset TM\rightarrow M\times M$ is a discretization map, then $R^*_d: {\overline{U}}\subset TM\rightarrow M\times M$ with ${\overline{U}}=\{ (x, v)\in TM\; \mid \, (x,-v)\in U\}$ and defined by

$$\begin{aligned} R^*_d(x, v)=\left( I_M \circ R_d\right) (x, -v) \end{aligned}$$

is also a discretization map. The map $R^*_d$ is called the adjoint discretization map of $R_d$.

Proof

Using the notation

$$\begin{aligned} R_d^*(q, v)=((R^*)^1(x,v), (R^*)^2(x,v))=((R^*)^1_x(v), (R^*)^2_x(v))\, , \end{aligned}$$

the two properties in Definition 2.2 are satisfied by $R^*_d$:

• 1. $R^*_d(x, 0)=\left( I_M\circ R_d\right) (x, 0)=I_M(x, x)=(x,x).$

• 2. $T_{0_x}(R^*)_x^2-T_{0_x}(R^*)_x^1$ is the identity map on $T_xM$ because $T_{0_x}(R^*)_x^2=-T_{0_x}R^1_x$ and $T_{0_x}(R^*)_x^1=-T_{0_x}R^2_x$.$\square$

Definition 2.3

A discretization map is symmetric if $R^*_d=R_d$.

Example 2.1

Let us provide some examples of retraction maps typically used in the literature for the construction of numerical methods, see [30], that can be used to define discretization maps satisfying the properties in Definition 2.2.

• 1. The explicit Euler method: $R_d(x,v)=(x,x+v)$ being its adjoint $R^*_d(x,v)=(x-v,x)$ .

• 2. The implicit midpoint rule is a symmetric discretization map: $R_d(x,v)=\left( x-\dfrac{v}{2},x+\dfrac{v}{2}\right)$.

• 3. The $\theta$-method: $R_d(x,v)=\left( x-\theta \, v,x+(1-\theta )\, v\right)$ where $\theta \in [0,1]$.

As known, for $\theta \in \{0, 1/2\}$, we recover the first two maps from the third one. All these methods are defined on the Euclidean vector space ${{\mathbb {R}}}^n$. $\square$

Example 2.2

Given a Riemannian manifold (M, g) and the associated exponential map $exp_x: T_xM\rightarrow M$ we can define the following map

$$\begin{aligned} R_d(x,\xi )=\left( exp_x(-\xi /2), exp_x(\xi /2)\right) \, , \end{aligned}$$

(3)

that satisfies the properties in Definition 2.2 and it is a symmetric discretization map. Let us give some specific examples of discretization maps that can be associated with the exponential map.

For instance, on the sphere $S^2$ with the Riemannian metric induced by the restriction of the standard metric on ${{\mathbb {R}} }^3$ we have that

$$\begin{aligned} exp_x(\xi )=\cos ( \Vert \xi \Vert ) \, x+\sin (\Vert \xi \Vert ) \, \frac{\xi }{\Vert \xi \Vert }, \qquad \xi \in T_x S^2. \end{aligned}$$

Thus, we move along the greatest circle that are the geodesics on the sphere. Remember that $exp_x(0_x)=x$ and the exponential map is a continuous map. Hence, we can define the following discretization map on M:

$$\begin{aligned} R_d (x, \xi )=\left( \cos \left( \frac{\Vert \xi \Vert }{2}\right) x-\sin \left( \frac{\Vert \xi \Vert }{2}\right) \frac{ \xi }{\Vert \xi \Vert }, \cos \left( \frac{\Vert \xi \Vert }{2}\right) x+\sin \left( \frac{\Vert \xi \Vert }{2}\right) \frac{ \xi }{\Vert \xi \Vert }\right) . \nonumber \\ \end{aligned}$$

(4)

Another option is to use as a retraction map on the sphere the projection $R_x(\xi )= \frac{x+\xi }{\Vert x+\xi \Vert }$ that leads to the following discretization map:

$$\begin{aligned} R_d (x, \xi )=\left( \frac{x-\xi /2}{\Vert x-\xi /2\Vert }, \frac{x+\xi /2}{\Vert x+\xi /2\Vert }\right) . \end{aligned}$$

Proposition 2.2 for $\theta =1/2$ guarantees that both maps are discretization maps. $\square$

Example 2.3

Consider a Lie group G and denote by ${{\mathfrak {g}}}$ its Lie algebra. It is a fact that any element $\xi$ in the Lie algebra is in one-to-one correspondence with a left-invariant vector field on G, $X_{\xi }=T_eL_g(\xi )$, where e is the identity element of G and $L_g: G\rightarrow G$ denotes the left-translation map. If $\gamma _{\xi }: {{\mathbb {R}}}\rightarrow G$ is an integral curve of $X_{\xi }$ with initial condition $\gamma _{\xi }(0)=e$, then we can generate a map between the Lie algebra and the Lie group called the exponential map: $\text {exp}(\xi )=\gamma _{\xi }(1)$. It is possible to check that the map $R: TG\rightarrow G$ given by

$$\begin{aligned} (g, X)\longrightarrow g\, \text {exp}(T_g L_{g^{-1}}(X)) \end{aligned}$$

is a retraction map where $X(g)\in T_gG$. Then, we define a symmetric discretization map on the Lie group G, $R_d: TG\rightarrow G\times G$, as follows

$$\begin{aligned} R_d (g, X)=\left( g\, \text {exp}\left( - \frac{1}{2}T_g L_{g^{-1}}(X)\right) , \, g\,\text {exp}\left( \frac{1}{2}T_g L_{g^{-1}}(X)\right) \right) \,. \end{aligned}$$

The properties in Definitions 2.1 and 2.2 are satisfied because the tangent map $T_e\exp$ is the identity map.

In the case of $SO(3)=\{ A \in GL(3, {{\mathbb {R}}})\; \mid \, AA^T=A^TA={\mathrm{Id}}_3,\ \det A=1\}$ we have that an element $(A, X)\in TSO(3)$ is given by a pair of matrices such that $A\in SO(3)$ and $XA^T+AX^T=0$. Therefore, the Lie algebra ${\mathfrak so}(3)$ is the set of skew-symmetric matrices: $\xi =A^TX\in {\mathfrak so}(3)$. The above retraction map for SO(3) becomes:

$$\begin{aligned} R(A,X)=A\,\text {exp} (A^TX) \,. \end{aligned}$$

The exponential map could be replaced by the Cayley transformation:

$$\begin{aligned} \begin{array}{rcl} \text {cay}: {{\mathfrak {s}}}{{\mathfrak {o}}}(3) &{}\longrightarrow &{} SO(3)\\ \xi &{}\longmapsto &{} \text {cay} (\xi )=({\mathrm{Id}}_3-\xi /2)^{-1}({\mathrm{Id}}_3+\xi /2)\, , \end{array} \end{aligned}$$

where ${\mathrm{Id}}_3$ stands for the identity matrix. Then we define the following retraction map $R_{{\mathrm{cay}}}:TSO(3)\rightarrow SO(3)$:

$$\begin{aligned} R_{{\mathrm{cay}}}(A,X)= A\,\text {cay} (A^TX)=A({\mathrm{Id}}_3-A^TX/2)^{-1}({\mathrm{Id}}_3+A^TX/2)\,. \end{aligned}$$

(5)

Using Proposition 2.2, we obtain the following discretization map

$R_{d,{\mathrm{cay}}}:TSO(3)\rightarrow SO(3)\times SO(3)$:

$$\begin{aligned} R_{d,{\mathrm{cay}}}(A,X)&= \left( R_{{\mathrm{cay}}}(A,-X/2), R_{{\mathrm{cay}}}(A,X/2)\right) \\&= \left( A({\mathrm{Id}}_3+A^TX/4)^{-1}({\mathrm{Id}}_3-A^TX/4), A({\mathrm{Id}}_3-A^TX/4)^{-1} \right. \\&\left. \qquad ({\mathrm{Id}}_3+A^TX/4)\right) \,. \end{aligned}$$

$\square$

3 Lift of Discretization Maps

We can construct discretization maps, as described in Definition 2.2, on any manifold. When studying mechanical systems, it may be useful to define discretization maps on the tangent bundle for the Lagrangian framework or on the cotangent bundle for the Hamiltonian framework. As discretization maps can be defined on different manifolds, we introduce the notation $R_d^{TM}$ so that the superscript tells us the domain of such a map. Thus, the map $R_d^{TM}:TM \rightarrow M\times M$ is called a discretization map on M. Note that “on M" emphasizes where the image takes values. The manifold M could be equal to the tangent bundle TQ or to the cotangent bundle $T^*Q$ depending on the dynamics under study.

Here, we are interested in constructing specific discretization maps on the tangent and cotangent bundles obtained from discretization maps on the base manifold. The objective is to generate geometric integrators for mechanical systems by using a suitable notion of lifted discretization maps to the tangent and cotangent bundles to encompass both the Lagrangian and the Hamiltonian framework.

We first review the notion of tangent and cotangent lift of a map between manifolds, see [40].

Let $M_1$ and $M_2$ be n-dimensional manifolds and $F: M_1\rightarrow M_2$ be a smooth map. The tangent lift $TF: TM_1\rightarrow TM_2$ of F is defined by

$$\begin{aligned} TF(v_x)=T_xF (v_x)\in T_{F(x)} M_2\, , \qquad \text{ where } v_x\in T_xM_1\; , \end{aligned}$$

and $T_xF$ is the tangent map of F whose matrix is the Jacobian matrix of F at $x\in M_1$ in a local chart.

As the tangent map $T_xF$ is linear, the dual map $T_{x}^*F:T^*_{F(x)}M_2\rightarrow T^*_xM_1$ is defined as follows:

$$\begin{aligned} \langle (T^*_{x}F)(\alpha _2), v_{x}\rangle =\langle \alpha _2, T_{x}F(v_{x})\rangle \text{ for } \text{ every } v_x\in T_xM_1. \end{aligned}$$

Note that $(T^*_{x}F)(\alpha _2)\in T^*_xM_1$.

To define the cotangent lift in Sect. 3.2, we need the cotangent lift of the inverse of the discretization map. Thus, we fix the notation for such a cotangent lift.

Definition 3.1

Let $F: M_1\rightarrow M_2$ be a diffeomorphism. The vector bundle morphism $\widehat{F}: T^*M_1\rightarrow T^*M_2$ defined by

$$\begin{aligned} \widehat{F}=T^*F^{-1} \end{aligned}$$

is called the cotangent lift of $F^{-1}$.

In other words, $\widehat{F}(\alpha _x)= T^*_{F(x)}F^{-1} (\alpha _x)$ where $\alpha _x\in T^*_x M_1$. Obviously, $(T^*F^{-1})\circ (T^*F)={\mathrm{Id}}_{T^*M_2}$.

We quickly review here some notions from symplectic geometry, see [37]. Denote by $\pi _M: T^*M\rightarrow M$ the canonical projection of the cotangent bundle and define the Liouville 1-form $\theta _M$ on $T^*M$ by $\langle \theta _M(\alpha _x), X_{\alpha _x}\rangle = \langle \alpha _x, T_{\alpha _x}\pi _M( X_{\alpha _x})\rangle$ where $X_{\alpha _x}\in T_{\alpha _x} T^*M$ and denote by $\omega _{M}=-d\Theta _M$ the canonical symplectic 2-form on $T^*M$. Thus, $(T^*M,\omega _M)$ is a symplectic manifold. For a diffeomorphism $F: M_1\rightarrow M_2$, we recall the well-known proposition for symplectic manifolds in [37].

Proposition 3.1

Let $F: M_1\rightarrow M_2$ be a diffeomorphism. The cotangent lift $\widehat{F}: T^*M_1\rightarrow T^*M_2$ of $F^{-1}$ is a symplectomorphism for the symplectic manifolds $(T^*M_1, \omega _{M_1})$ and $(T^*M_2, \omega _{M_2})$. In other words, the symplectic 2-form is preserved by the pull-back of ${\hat{F}}$:

$$\begin{aligned} \widehat{F}^*(\omega _{M_2})=\omega _{M_1}\, \text{ where } \widehat{F}^*:\Omega ^2(T^*M_2) \rightarrow \Omega ^2(T^*M_1). \end{aligned}$$

Equivalently, the inverse of the cotangent lift $\widehat{F}^{-1}:T^*M_2\rightarrow T^*M_1$ is also a symplectomorphism.

Some expressions in coordinates will be useful in the sequel. Take local coordinates $q=(q^1,\ldots , q^n)$ on $M_1$ and $x=(x^1,\ldots , x^m)$ on $M_2$ and induced coordinates (q, v) on $TM_1$ and (x, u) on $TM_2$, respectively. If $F: M_1\rightarrow M_2$ is written in local coordinates as $(q^1, \ldots , q^n)\rightarrow (F^1(q), \ldots , F^m( q))$ Then

$$\begin{aligned} TF({q}, {v})= & {} \left( F^i( q)\; ;\; \frac{\partial F^i}{\partial q^j}(q)v^j\right) \; . \end{aligned}$$

Taking now induced coordinates (q, p) on $T^*M_1$ and (x, r) on $T^*M_2$, we have

$$\begin{aligned} \widehat{F}(q, p)= & {} \left( F^i( q)\; ;\; p_j\frac{\partial (F^{-1})^j}{\partial q^i} (F( q))\right) \; . \end{aligned}$$

We could also use the matrix notation:

$$\begin{aligned} \begin{aligned} D_{q}F=&{} \left( \frac{\partial F^i}{\partial q^j}(q)\right) _{1\le i, j\le \dim M_1} \quad \text { and } \quad \\ D_{F(q)}F^{-1}=&{} \left( \frac{\partial (F^{-1})^i}{\partial x^j}(F( q))\right) _{1\le i, j\le \dim M_2}\,. \end{aligned} \end{aligned}$$

Note that

$$\begin{aligned} D_{F( q)}F^{-1}=\left[ D_{ q}F\right] ^{-1}\; . \end{aligned}$$

When we restrict the previous maps TF and $\widehat{F}$ to a fiber, we induce the maps

$$\begin{aligned} \begin{array}{rrcl} T_{q}F:&{} T_{ q} M_1&{}\rightarrow &{} T_{F( q)}M_2\\ &{} { v}&{}\longmapsto &{} D_{ q}F\; {v}^T \end{array} \end{aligned}$$

and

$$\begin{aligned} \begin{array}{rrcl} \widehat{F}_{ q}:&{} T^*_{ q} M_1&{}\rightarrow &{} T^*_{F(q)}M_2\\ &{} { p}&{}\longmapsto &{} ((D_{ q} F)^{-1})^T { p}^T= \left( { p} (D_{ q} F)^{-1}\right) ^T. \end{array} \end{aligned}$$

Consequently,

$$\begin{aligned} \begin{array}{rrcl} \widehat{F}^{-1}_{F(q)}:&{} T^*_{F( q)} M_2&{}\rightarrow &{} T^*_{ q}M_1\\ &{} {r}&{}\longmapsto &{} { r} D_{ q} F. \end{array} \end{aligned}$$

(6)

3.1 Tangent Lift of Discretization Maps

We prove that if we suitably lift the discretization map $R_d:TQ \rightarrow Q\times Q$ on Q in Definition 2.2, we obtain a new discretization map on the tangent bundle TQ. These constructions are able to provide a geometric framework to obtain numerical integrators for second-order differential equations (SODEs), see Sect. 4, and for the dynamics of mechanical systems as shown in Sects. 5 and 6.

Remember that the notation $R_d^{TTQ}$ for a discretization map on TQ makes clear the manifold to be discretized, that is, $R_d^{TTQ}:TTQ \rightarrow TQ\times TQ$. To define it from a discretization map $R_d:{TQ}\rightarrow Q\times Q$ on Q is necessary to use the canonical involution map $\kappa _Q$ that shows the double vector bundle structure of the vector bundle TTQ and defines a vector bundle isomorphism, as described, for instance, in [50, 52].

Let us recall here the definition of the canonical involution. Let Q be a smooth manifold of dimension n, $\tau _{Q}:TQ\rightarrow Q$ be the canonical tangent bundle projection and TTQ the double tangent bundle of Q. The manifold TTQ naturally admits two vector bundle structures. The first vector bundle structure is the canonical one with vector bundle projection $\tau _{TQ}: TTQ\rightarrow TQ$. For the second vector bundle structure of TTQ, the vector bundle projection is given by the tangent map $T\tau _{Q}:TTQ\rightarrow TQ$. The canonical involution $\kappa _{Q}: TTQ\rightarrow TTQ$ is a vector bundle isomorphism (over the identity of TQ) between the two previous vector bundles. In fact, $\kappa _{Q}$ is characterized by the following condition: let $\Phi :U\subseteq {{\mathbb {R}}}^{2}\rightarrow Q$ be a smooth map on an open subset U of ${{\mathbb {R}}}^{2}$ defined by

$$\begin{aligned} (t,s)\mapsto \Phi (t,s)\in Q, \end{aligned}$$

then

$$\begin{aligned} \kappa _{Q}\left( \frac{\partial }{\partial t}\frac{\partial }{\partial s} \Phi (t,s) \right) = \frac{\partial }{\partial s}\frac{\partial }{\partial t} \Phi (t,s). \end{aligned}$$

Note that $\kappa _{Q}$ is an involution of TTQ, that is, $\kappa _{Q}^{2}={\mathrm{Id}}_{TTQ}$. If (q, v) are canonical fibered coordinates of TQ and $(q,v,{\dot{q}},{\dot{v}})$ are the corresponding local fibered coordinates of TTQ, then

$$\begin{aligned} \kappa _{Q}(q,v, {\dot{q}},{\dot{v}})=(q,{\dot{q}}, v, {\dot{v}}). \end{aligned}$$

Having all this in mind, remember that the tangent lift of a vector field X on Q does not define a vector field on TQ. It is necessary to consider the composition $\kappa _Q\circ TX$ to obtain a vector field on TQ that is called complete lift $X^c$ of the vector field X. A similar trick must be used to lift a discretization map from TQ to TTQ as shown in the following diagram.

Note that $T(Q\times Q)$ and $TQ\times TQ$ are trivially identified since any vector on $T(Q\times Q)$ is given as a tangent vector at 0 of a curve $\sigma : {{\mathbb {R}}}\rightarrow Q\times Q$, that is, ${\dot{\sigma }}(0)\in T_{\sigma (0)}(Q\times Q)$. As $\sigma$ has two components $\sigma (t)=(\sigma _1 (t), \sigma _2(t))$ where $\sigma _i: {{\mathbb {R}}}\rightarrow Q$, $i=1,2$, the identification ${\dot{\sigma }}(0)\equiv ({\dot{\sigma }}_1(0),{\dot{\sigma }}_2(0))\in T_{\sigma _1(0)}Q\times T_{\sigma _2(0)}Q$ is made. The following proposition shows that ${\mathrm{T}}R_d\circ \kappa _Q$ is a discretization map on TQ. From now on, such a map is denoted by $R_d^T$ to emphasize it is obtained by tangently lifting $R_d$.

Proposition 3.2

If $R_d$ is a discretization map on Q, then $R_d^{T}={\mathrm{T}}R_d\circ \kappa _Q$ is a discretization map on TQ.

Proof

In local coordinates $(q,v,{\dot{q}},{\dot{v}})$ of TTQ we have that ${\mathrm{T}}R_d(q,v,{\dot{q}},{\dot{v}})=(R_d(q,v),{\mathrm{D}}_{(q,v)} R_d(q,v) \, ({\dot{q}},{\dot{v}})^T)$ and

$$\begin{aligned} R_d^{T}(q,{\dot{q}},v,{\dot{v}})=(R_d(q,v),{\mathrm{D}}_{(q,v)} R_d \; ({\dot{q}},{\dot{v}})^T)\; . \end{aligned}$$

Remember the abuse of notation because $T(Q\times Q)$ and $TQ \times TQ$ are trivially identified.

Let us prove that the properties in Definition 2.2 are satisfied by $R_d$ knowing that $R_d(q,v)=(R^1(q,v),R^2(q,v))$.

• 1. We know that $R_d(q,0)=(q,q)$ for all $q\in Q$. Consequently,

  $$\begin{aligned} R_d^{T}(q,{\dot{q}},0,0)=(R_d(q,0),{\mathrm{D}}_{(q,0)} R_d \, ({\dot{q}},0)^T)=(q,q,{\dot{q}},{\dot{q}})\equiv (q,{\dot{q}}; q,{\dot{q}})\, , \end{aligned}$$

  where we use the natural identification between $T(Q\times Q)$ and $TQ\times TQ$.

• 2. For the second property, we know that

  $$\begin{aligned} R_d^{T}(q,{\dot{q}},v,{\dot{v}})=((TR^1)(q,v; {\dot{q}},{\dot{v}}), (TR^2)(q, v; {\dot{q}},{\dot{v}})). \end{aligned}$$

  We need to compute

  $$\begin{aligned} {\mathrm{T}}_{(0,0)_{(q,{\dot{q}})}}({\mathrm{T}}R^a)_{(q, {\dot{q}})}: T_{(0,0)_{(q,{\dot{q}})}}T_{(q, {\dot{q}})}TQ\equiv T_{(q, {\dot{q}})}TQ\rightarrow T_{(q, {\dot{q}})}TQ\; , \end{aligned}$$

  for $a=1,2$, to prove that the map ${\mathrm{T}}_{(0,0)_{(q,{\dot{q}})}}({\mathrm{T}}R^2)_{(q, {\dot{q}})}-{\mathrm{T}}_{(0,0)_{(q,{\dot{q}})}}({\mathrm{T}}R^1)_{(q,{\dot{q}})}$ is the identity map understood as an application from $T_{(q,{\dot{q}})}TQ$ to itself.

  At $(q,{\dot{q}},0,0)$, the linear map ${\mathrm{T}}_{(0,0)_{(q,{\dot{q}})}}({\mathrm{T}}R^a)_{(q, {\dot{q}})}$ is given by the following matrix

  $$\begin{aligned} \begin{pmatrix} \partial _{v^j}(R^a)^i(q, 0) &{} 0\\ \partial _{q^k} \partial _{v^j}(R^a)^i(q,0){\dot{q}}^k &{} \partial _{v^j}(R^a)^i(q,0) \end{pmatrix}\, , \end{aligned}$$

  after calculating

  $$\begin{aligned} \left. \frac{d}{dt}\right| _{t=0} \left( R_d^a(q, tv), \partial _{q^j} R_d^a (q, tv){\dot{q}}^j+\partial _{v^j} R_d^a (q, tv)t{\dot{v}}^j\right) \, . \end{aligned}$$

  Using again the properties of the discretization map $R_d$, the Jacobian matrix of $({\mathrm{T}}R^2)_{(q, {\dot{q}})}-({\mathrm{T}}R^1)_{(q,{\dot{q}})}$ at $(0,0)_{(q,{\dot{q}})}$ is:

  $$\begin{aligned} \begin{pmatrix} \partial _vR^2(q, 0)- \partial _vR^1(q, 0)&{} 0\\ \partial _q(\partial _vR^2(q,0)-\partial _vR^1(q,0)){\dot{q}} &{} \partial _vR^2(q, 0)- \partial _vR^1(q, 0) \end{pmatrix}= {\mathrm{Id}}_{2n\times 2n}\, , \end{aligned}$$

  as needed. Note that $\partial _q\partial _v(R^2-R^1)(q,0)=0$ because $\partial _v(R^2-R^1)(q,0)={\mathrm{Id}}_{n\times n}$.$\square$

Remark 3.1

If we use the discretization map obtained from the exponential map of a Riemannian metric g as in Eq. (3), then the tangent lift of this specific discretization map is associated with the complete lift of g, denoted by $g^C$, which is a semi-Riemannian metric on TQ (see details in [5, 55]). $\square$

Proposition 3.3

Let $R_d$ be a discretization map and $R_d^*$ be the adjoint discretization map. Then the tangent lift of a symmetric discretization map is also symmetric, that is, $(R_d^{*})^T= (R_d^T)^{*}$.

Proof

It is simple to check that

$$\begin{aligned} (R_d^*)^T(q, {\dot{q}}, v, {\dot{v}})= & {} TR_d^*(q, v, {\dot{q}}, {\dot{v}}) =(T(R^*)^1(q, v, {\dot{q}}, {\dot{v}}), T(R^*)^2(q, v, {\dot{q}}, {\dot{v}}))\\= & {} (TR^2(q, v, -{\dot{q}}, -{\dot{v}}), TR^1(q, v, -{\dot{q}}, -{\dot{v}}))\\= & {} I_{TQ}(TR^1(q, v, -{\dot{q}}, -{\dot{v}}), TR^2(q, v, -{\dot{q}}, -{\dot{v}}))\\= & {} (R_d^T)^*(q, {\dot{q}}, v, {\dot{v}}) \end{aligned}$$

where $I_{TQ}(u_q, v_q)=(v_q, u_q)$ for all $u_q, v_q\in T_qQ$. $\square$

3.2 Cotangent Lift of Discretization Maps

To encompass the Lagrangian and Hamiltonian dynamics together to build numerical integrators, we are interested in defining a very particular notion of discretization map on the cotangent bundle.

Given a discretization map $R_d: TQ\rightarrow Q\times Q$ we know that the cotangent lift $\widehat{R_d}: T^*TQ\rightarrow T^*(Q\times Q)$ is a symplectomorphism between the symplectic manifolds $(T^*TQ, \omega _{TQ})$ and $(T^*(Q\times Q), \omega _{Q\times Q})$ as mentioned in Proposition 3.1.

According to Definition 3.1, in local coordinates $(q,v,p_q,p_v)$ for $T^*TQ$ the cotangent lift of $R_d$ is given by:

$$\begin{aligned} \widehat{R_d}:T^*TQ\longrightarrow & {} T^*(Q\times Q) \\ (q,v,p_q,p_v)\longmapsto & {} \left( R_d(q,v), \left( p_q, \; p_v\right) \, (D_{(q,v)}R_d)^{-1}\right) \end{aligned}$$

where $(D_{(q,v)}R_d)^{-1}$ is the inverse of the Jacobian matrix of $R_d$.

We use the cotangent lift $\widehat{R_d}$ of the discretization map on Q to define a discretization map on $T^*Q$ that must be a map from $TT^*Q$ to $T^*Q\times T^*Q$.

For this purpose it is necessary to use the canonical symplectomorphism $\alpha _Q: TT^*Q\rightarrow T^*TQ$ between double vector bundles (see [51, 52]). Locally,

$$\begin{aligned} \alpha _Q:TT^*Q\longrightarrow & {} T^*TQ\\ (q,p,{\dot{q}},{\dot{p}})\longrightarrow & {} (q,{\dot{q}},{\dot{p}},p). \end{aligned}$$

As described in [51], the symplectomorphism $\alpha _Q$ is between the symplectic manifold $(TT^*Q, {\mathrm{d}}_T\omega _Q)$ and the natural symplectic manifold $(T^*TQ, \omega _{TQ})$. Recall that in local coordinates $(q,p,{\dot{q}},{\dot{p}})$ for $TT^*Q$, the symplectic form $d_T\omega _Q$ has the following expression: ${\mathrm{d}}_T \omega _Q= {\mathrm{d}}q\wedge {\mathrm{d}}{\dot{p}}+{\mathrm{d}}{\dot{q}}\wedge {\mathrm{d}}p$. Moreover, we need the diffeomorphism

$$\begin{aligned} \Phi :T^*Q\times T^*Q\longrightarrow & {} T^*(Q\times Q)\\ (q_0, p_0; q_1, p_1)\longmapsto & {} (q_0, q_1, -p_0, p_1) \end{aligned}$$

which is also a symplectomorphism between $(T^*(Q\times Q), \omega _{Q\times Q})$ and $(T^*Q\times T^*Q, \Omega _{12}=pr_2^*\omega _Q-pr^*_1\omega _Q)$, where ${\mathrm{pr}}_i:T^*(Q\times Q) \rightarrow T^*Q\times T^*Q$ denotes the projection into the i–th factor of the Cartesian product in the image.

The following diagram shows how to define the discretization map on $T^*Q$ from the one on Q.

Now we prove that $R_d^{T^*}$ is a discretization map on $T^*Q$ according to Definition 2.2. From now on, it will be called the cotangent lift of $R_d$.

Proposition 3.4

Let $R_d:TQ\rightarrow Q\times Q$ be a discretization map on Q as in Definition 2.2. Then $R_d^{T^*}=\Phi ^{-1}\circ \widehat{R_d}\circ \alpha _Q:TT^*Q\rightarrow T^*Q\times T^*Q$ is a discretization map on $T^*Q$.

Proof

Let us compute the cotangent lift of the tangent map of $R_d$ for local coordinates $(q,v,p_q,p_v)$ of $T^*TQ$:

$$\begin{aligned} \widehat{R_d}(q,v,p_q,p_v)= \left( R_d(q,v), (p_q,\; p_v) ({\mathrm{D}}_{(q,v)}R_d)^{-1}\right) \,. \end{aligned}$$

Expressing the inverse of $({\mathrm{D}}R_d)^{-1}$ as a matrix with two blocks $({\mathrm{D}}R_d)_i^{-1}$ of size $2n\times n$, $i=1,2$, that is

$$\begin{aligned} ({\mathrm{D}}_{(q,v)}R_d)^{-1}= \begin{pmatrix} ({\mathrm{D}}_{(q,v)}R_d)_1^{-1}&({\mathrm{D}}_{(q,v)}R_d)_2^{-1} \end{pmatrix} \end{aligned}$$

and

$$\begin{aligned} ({\mathrm{D}}_{(q,v)}R_d)^{-1}=\begin{pmatrix} \partial _q R^1(q,v)&{}\partial _v R^1(q,v) \\ \partial _q R^2(q,v)&{}\partial _v R^2(q,v) \end{pmatrix}^{-1}. \end{aligned}$$

We can write

$$\begin{aligned}&R_d^{T^*}(q,p,{\dot{q}},{\dot{p}})\\&\quad = \left( R^1(q,{\dot{q}}),- ({\dot{p}}, \; p) ({\mathrm{D}}_{(q,{\dot{q}})}R_d)_1^{-1}; R^2(q,{\dot{q}}),({\dot{p}}, \; p) ({\mathrm{D}}_{(q,{\dot{q}})}R_d)_2^{-1}\right) \,. \end{aligned}$$

Let us check if it satisfies the properties in Definition 2.2:

• 1. Note that the Jacobian matrix of $R_d$ at (q, 0) is

  $$\begin{aligned} {\mathrm{D}}_{(q,0)} R_d=\begin{pmatrix} {\mathrm{Id}} &{} \partial _vR^1(q,0)\\ {\mathrm{Id}} &{} \partial _vR^2(q,0)\end{pmatrix}. \end{aligned}$$

  As $\partial _vR^2(q,0)-\partial _vR^1(q,0)={\mathrm{Id}}$, the inverse is

  $$\begin{aligned} ({\mathrm{D}}_{(q,0)} R_d)^{-1}=\begin{pmatrix} {\mathrm{Id}}+\partial _vR^1(q,0)&{} -\partial _vR^1(q,0)\\ -{\mathrm{Id}}&{}{\mathrm{Id}}\end{pmatrix}. \end{aligned}$$

  Thus,

  $$\begin{aligned}&R_d^{T^*}(q,p,0,0)\nonumber \\&\quad = (R^1(q,0),-(0, \; p) ({\mathrm{D}}_{(q,0)}R_d)_1^{-1}; R^2(q,0),(0, \; p) ({\mathrm{D}}_{(q,0)}R_d)_2^{-1})\,, \end{aligned}$$

  and it is straightforward that $R_d^{T^*}(q,p,0,0)=(q,p; q,p)$.

• 2. We must prove that $T_{(q,p,0,0)}\left( R^{T^*}_d\right) ^2_{(q,p)}-T_{(q,p,0,0)}\left( R^{T^*}_d\right) _{(q,p)}^1$ is the identity map from $T_{(q,p,0,0)}TT^*Q\simeq T_{(q,p)}T^*Q$ to itself.

  Let us compute the following derivatives for $a=1,2$:

  $$\begin{aligned} \left. \frac{d}{dt}\right| _{t=0}\left( R^{T^*}_d\right) ^a(q,p, t{\dot{q}}, t{\dot{p}}). \end{aligned}$$

  For instance, for $i=1$ we have

  $$\begin{aligned} \left. \frac{d}{dt}\right| _{t=0}\left( R^{T^*}_d\right) ^1(q,p, t{\dot{q}}, t{\dot{p}})= \left. \frac{d}{dt}\right| _{t=0}\left[ R^1(q,t{\dot{q}}),- (t{\dot{p}}, \; p) ({\mathrm{D}}_{(q,t{\dot{q}})}R_d)_1^{-1}\right] \; . \end{aligned}$$

  Using the expression for the derivative of an inverse matrix, we have that $\left. \frac{d}{dt}\right| _{t=0}({\mathrm{D}}_{(q, t{\dot{q}})}R_d)^{-1}$ is equal to

  $$\begin{aligned}&-\begin{pmatrix} {\mathrm{Id}}+A&{} -A\\ -{\mathrm{Id}}&{}{\mathrm{Id}}\end{pmatrix} \begin{pmatrix} \partial _{v}\partial _qR^1(q,0)&{}\partial _{v}\partial _vR^1(q,0)\\ \partial _{v}\partial _qR^2(q,0)&{}\partial _{v}\partial _vR^2(q,0) \end{pmatrix} \begin{pmatrix} {\mathrm{Id}}+A&{} -A\\ -{\mathrm{Id}}&{}{\mathrm{Id}}\end{pmatrix}\\&\quad =\begin{pmatrix} {\mathbf *}&{}{\mathbf *}\\ \partial _{v}\partial _vR^1(q,0)-\partial _{v}\partial _vR^2(q,0) &{}\partial _{v}\partial _vR^2(q,0)-\partial _{v}\partial _vR^1(q,0)\end{pmatrix} \end{aligned}$$

  where $A=\partial _vR^1(q,0)$ and $(*)$ denotes terms that are not explicitly needed in the computations. We have used that $\partial _q\partial _v(R^2-R^1)(q,0)=0$ since $\partial _v(R^2-R^1)(q,0)={\mathrm{Id}}_{n\times n}$. Thus,

  $$\begin{aligned}&\left. \frac{d}{dt}\right| _{t=0}\left[ R^1(q,t{\dot{q}}),- (t{\dot{p}}, \; p) ({ D}_{(q,t{\dot{q}})}R_d)_1^{-1}\right] \\&\quad = \begin{pmatrix} \partial _vR^1(q,0)&{}0\\ p(\partial _{v}\partial _vR^1(q,0)-\partial _{v}\partial _vR^2(q,0))&{}-{\mathrm{Id}}-\left( \partial _v R^1(q,0)\right) ^T \end{pmatrix} \begin{pmatrix} {\dot{q}}\\ {\dot{p}} \end{pmatrix} \end{aligned}$$

  Analogously,

  $$\begin{aligned}&\left. \frac{d}{dt}\right| _{t=0}\left[ R^2(q,t{\dot{q}}), (t{\dot{p}}, \; p) ({ D}_{(q,t{\dot{q}})}R_d)_2^{-1}\right] \\&\quad = \begin{pmatrix} \partial _vR^2(q,0)&{}0\\ p(\partial _{v}\partial _vR^1(q,0)-\partial _{v}\partial _vR^2(q,0))&{}-\left( \partial _v R^1(q,0)\right) ^T \end{pmatrix} \begin{pmatrix} {\dot{q}}\\ {\dot{p}} \end{pmatrix} \end{aligned}$$

  As a result,

  $$\begin{aligned} \begin{pmatrix} \partial _vR^2(q,0)&{}0\\ C&{}-\partial _v R^1(q,0) \end{pmatrix} -\begin{pmatrix} \partial _vR^1(q,0)&{}0\\ C&{}-{\mathrm{Id}}-\partial _v R^1(q,0) \end{pmatrix} =\begin{pmatrix} {\mathrm{Id}}&{}0\\ 0&{}{\mathrm{Id}} \end{pmatrix} \end{aligned}$$

  where $C=p(\partial _{v}\partial _vR^1(q,0)-\partial _{v}\partial _vR^2(q,0))$.$\square$

As the composition of symplectomorphisms is a symplectomorphism [37], the following result is straightforward.

Proposition 3.5

Let $R_d:TQ \rightarrow Q\times Q$ be a retraction map on Q, then $R_d^{T^*}=\Phi ^{-1}\circ \widehat{R_d} \circ \alpha _Q:T(T^*Q)\rightarrow T^*Q\times T^*Q$ is a symplectomorphism between $(T(T^*Q), {\mathrm{d}}_T \omega _Q)$ and $(T^*Q\times T^*Q, \Omega _{12})$.

As a consequence,

$$\begin{aligned} \left( R_d^{T^*}\right) ^*(\Omega _{12})=d_T\omega _Q. \end{aligned}$$

The above result is essential to obtain symplectic methods in the following sections.

When constructing numerical integrators in Sect. 5 for Hamiltonian systems, the inverse map of $R_d^{T^*}: TT^*Q\rightarrow T^*Q\times T^*Q$ is useful. Using Proposition 3.4, we specifically write the inverse map

$$\begin{aligned} \left( R_d^{T^*}\right) ^{-1}=\alpha ^{-1}_Q\circ \widehat{R_d}^{-1} \circ \Phi : T^*Q\times T^*Q\rightarrow TT^*Q . \end{aligned}$$

In local coordinates $(q_0,p_0;q_1,p_1)$ for $T^*Q\times T^*Q$ and using (6), it is quite simple to compute the inverse map

$$\begin{aligned} \left( R_d^{T^*}\right) ^{-1}(q_0,p_0; q_1,p_1)= \alpha _Q^{-1}\left( R_d^{-1} (q_0,q_1), (-p_0, p_1)\, { D}_{R_d^{-1} (q_0,q_1)}R_d\right) . \end{aligned}$$

(7)

Remember that $\alpha _Q^{-1}(q,v, p_q, p_v)=(q, p_v, v, p_q)$.

3.3 Duality Between the Cotangent and the Tangent Lift of Discretization Maps

After introducing both the tangent and cotangent lift of discretization maps, we show here the existing duality between the two maps.

For a discretization map on Q, we consider the tangent lift $R_d^{T}: TTQ\rightarrow TQ\times TQ$ defined by $R_d^{T}=TR_d\circ \kappa _Q$ and the corresponding cotangent lift $R_d^{T^*}=\Phi ^{-1}\circ \widehat{R_d}\circ \alpha _Q:TT^*Q\rightarrow T^*Q\times T^*Q$. As mentioned earlier, $TT^*Q$ is a symplectic manifold with the 2-form ${\mathrm{d}}_T \omega _Q$ that induces a natural pairing as follows. Let $v\in TT^*Q$ and let $w\in TTQ$ such that $\tau _{TQ}(w)=T\pi _Q(v)$, the pairing $\langle \cdot , \cdot \rangle ^T$ induced by the symplectic structure of $TT^*Q$ is given by

$$\begin{aligned} \langle v, \kappa _Q(w)\rangle ^T= \frac{d}{dt}\langle \sigma _v(t), \gamma _{{\tilde{w}}}(t)\rangle (0)=\langle \alpha _Q(v), w\rangle \, , \end{aligned}$$

where $\alpha _Q: TT^*Q\rightarrow T^*TQ$, $\sigma _v: I\rightarrow T^*Q$ and $\gamma _{{\tilde{w}}}: I\rightarrow TQ$ satisfy ${\dot{\sigma }}_v(0)=v$ and ${\dot{\gamma }}_{{\tilde{w}}}(0)={\tilde{w}}$ with ${\tilde{w}}=\kappa _Q(w)$ and $\pi _Q\circ \sigma _v=\tau _Q \circ \gamma _{{\tilde{w}}}$.

Proposition 3.6

The tangent lift and the cotangent lift of a discretization map on Q satisfy the following equality:

$$\begin{aligned} \left\langle \Phi (\alpha _{q_0}, \alpha _{q_1}), R_d^{T} (w)\right\rangle =\left\langle \left( R_d^{T^*}\right) ^{-1}(\alpha _{q_0}, \alpha _{q_1}), w\right\rangle ^T \, , \end{aligned}$$

where $w\in TTQ$, $(R_d)^{-1}(q_0, q_1)=T\tau _Q(w)$ and the pairing $\langle \cdot , \cdot \rangle ^T$ is induced by the symplectic structure of $TT^*Q$.

Proof

Observe that

$$\begin{aligned} \left\langle \Phi (\alpha _{q_0}, \alpha _{q_1}), R_d^{T} (w)\right\rangle= & {} \left\langle (\widehat{R_d}^{-1}\circ \Phi )(\alpha _{q_0}, \alpha _{q_1}), \kappa _Q(w)\right\rangle \\= & {} \left\langle (\alpha _Q)^{-1}\left( \left( \widehat{R_d}^{-1}\circ \Phi \right) (\alpha _{q_0}, \alpha _{q_1})\right) , w \right\rangle ^T\\= & {} \left\langle \left( R_d^{T^*}\right) ^{-1}(\alpha _{q_0}, \alpha _{q_1}), w\right\rangle ^T. \end{aligned}$$

$\square$

Using Propositions 3.3 and 3.6, it is easy to prove the following relation between the cotangent lift of the adjoint discretization and the adjoint of the cotangent lift discretization map.

Proposition 3.7

Let $R_d$ be a discretization map and $R_d^*$ be the adjoint discretization map. Then the cotangent lift of a symmetric discretization map is also symmetric, that is, $(R_d^*)^{T^*}= (R_d^{T^*})^*$.

3.4 Examples

We resume Examples 2.1, 2.2, 2.3 to construct the lifts of discretization maps described in the previous sections. In other words, we define discretization maps on TQ and $T^*Q$ starting from a discretization map on Q.

Example 3.1

We focus now on the mid-point rule described in Example 2.1 to define the tangent and cotangent lift of that symmetric discretization map. Assume that Q is a vector space and let $R_d:TQ \rightarrow Q\times Q$ be the discretization map induced by the mid-point rule as follows $R_d(q, v)=\left( q-\frac{1}{2}v, q+\frac{1}{2}v\right)$. If we compute the inverse map $R_d^{-1}(q_0,q_1)=\left( \dfrac{q_0+q_1}{2},q_1-q_0 \right)$, we construct the sequence of points that will be used either for optimization or numerical integration as the discrete flow

$$\begin{aligned} \phi _d :Q\times Q\rightarrow & {} Q\\ (q_0,q_1)\mapsto & {} \left( \tau _Q\circ R_d^{-1} \right) (q_0,q_1)=\tau _Q \left( \dfrac{q_0+q_1}{2},q_1-q_0\right) =\dfrac{q_0+q_1}{2}\,. \end{aligned}$$

Thus, the mid-point rule is recovered.

To define the tangent lift of the discretization map $R_d^{T}:TTQ \rightarrow TQ\times TQ$ on TQ, we first need to compute the tangent map whose matrix is

$$\begin{aligned} {\mathrm{D}} R_d= \begin{pmatrix} {\mathrm{Id}} &{} -\dfrac{1}{2}\, { \mathrm Id}\\ { \mathrm Id} &{} \dfrac{1}{2}\, {\mathrm{Id}} \end{pmatrix}. \end{aligned}$$

The tangent lift of $R_d$ is given by:

$$\begin{aligned} R_d^{T}(q,{\dot{q}},v,{\dot{v}})&=\left( TR_d \circ \kappa _Q\right) (q,{\dot{q}},v,{\dot{v}})=TR_d(q,v;{\dot{q}},{\dot{v}})\\&=\left( q-\dfrac{1}{2}\, v, q+\dfrac{1}{2}\, v; \;{\dot{q}}-\dfrac{1}{2}\, {\dot{v}}, {\dot{q}}+\dfrac{1}{2}\, {\dot{v}},\right) \\&\equiv \left( q-\dfrac{1}{2}\, v, {\dot{q}}-\dfrac{1}{2}\, {\dot{v}}; \; q+\dfrac{1}{2}\, v, {\dot{q}}+\dfrac{1}{2}\, {\dot{v}},\right) \,, \end{aligned}$$

where we naturally identify elements of $T(Q\times Q)$ with elements of $TQ\times TQ$.

We can also compute the inverse map:

$$\begin{aligned} \left( R_d^{T}\right) ^{-1}(q_0,v_0;q_1,v_1) =\left( \dfrac{q_0+q_1}{2},\dfrac{v_0+v_1}{2}; q_1-q_0,v_1-v_0 \right) . \end{aligned}$$

To compute the cotangent lift of $R_d$, that it, $R_d^{T^*}:TT^*Q\rightarrow T^*Q\times T^*Q$, we first need the tangent map of the inverse map $R^{-1}_d(q_0,q_1)=\left( \dfrac{q_0+q_1}{2},q_1-q_0\right)$:

$$\begin{aligned} D R^{-1}_d=\begin{pmatrix} \dfrac{1}{2}\, { \mathrm Id} &{} \dfrac{1}{2}\, { \mathrm Id} \\ - {\mathrm{Id}} &{} {\mathrm{Id}} \end{pmatrix}. \end{aligned}$$

Thus, the cotangent lift of $R_d^{-1}$ is given by:

$$\begin{aligned} \widehat{R_d}(q,v,p_q,p_v)&=\left( R_d(q,v), (p_q,\, p_v) (D_{(q,v)}R_d^{-1})\right) \\&= \left( q-\dfrac{1}{2}\, v, q+\dfrac{1}{2}\, v; \; \dfrac{p_q}{2}-p_v, \dfrac{p_q}{2}+p_v\right) . \end{aligned}$$

Finally, the cotangent lift of $R_d$ is the following discretization map on $T^*Q$:

$$\begin{aligned} R_d^{T^*}(q,p,{\dot{q}},{\dot{p}})&=\left( \Phi ^{-1}\circ \widehat{R_d} \circ \alpha _Q \right) (q,p,{\dot{q}},{\dot{p}})=\left( \Phi ^{-1}\circ \widehat{R_d}\right) (q,{\dot{q}},{\dot{p}},p)\nonumber \\&= \Phi ^{-1} \left( q-\dfrac{1}{2}\,{\dot{q}}, q+\dfrac{1}{2}\, {\dot{q}}; \; \dfrac{{\dot{p}}}{2}-p, \dfrac{{\dot{p}}}{2}+p\right) \nonumber \\&=\left( q-\dfrac{1}{2}\,{\dot{q}}, p-\dfrac{{\dot{p}}}{2}; \; q+\dfrac{1}{2}\, {\dot{q}}, p+\dfrac{{\dot{p}}}{2}\right) . \end{aligned}$$

(8)

The inverse map $\left( R_d^{T^*}\right) ^{-1}:T^*Q\times T^*Q \rightarrow TT^*Q$ is given by:

$$\begin{aligned} \left( R_d^{TT^*Q}\right) ^{-1}(q_0,p_0; q_1,p_1)=\left( \dfrac{q_0+q_1}{2},\frac{p_0+p_1}{2},q_1-q_0, p_1-p_0\right) . \end{aligned}$$

(9)

$\square$

Example 3.2

Let us lift the discretization map in Example 2.2. To simplify the computations we consider the discretization map that fixes the first point (compared with (4)) as follows

$$\begin{aligned} R_d (x, \xi )=\left( x, \cos \left( \Vert \xi \Vert \right) x+\sin \left( \Vert \xi \Vert \right) \frac{\xi }{\Vert \xi \Vert }\right) . \end{aligned}$$

Remember that $x\,\cdot x^T=1$ and $x\, \cdot \xi ^T=0$ because the manifold Q is the sphere $S^2$.

The tangent map of $R_d$ is given by the matrix

where

$$\begin{aligned} N_{ij}=-\sin \left( \Vert \xi \Vert \right) \, \dfrac{\xi _j x_i}{\Vert \xi \Vert }+ \cos \left( \Vert \xi \Vert \right) \, \dfrac{\xi _j \xi _i}{\Vert \xi \Vert }+\sin \left( \Vert \xi \Vert \right) \, \cdot {\left\{ \begin{array}{ll} \dfrac{\Vert \xi \Vert ^2-\xi _i\xi _i}{\Vert \xi \Vert ^3}, \quad \text {for} \; i=j\, ,\\ \\ \dfrac{-\xi _i\xi _j}{\Vert \xi \Vert ^3}, \quad \text {for} \; i\ne j\, ,\\ \end{array}\right. } \end{aligned}$$

are the entries of the invertible matrix $N(x,\xi )$. Thus, the tangent lift of $R_d$ is the following discretization map on TQ:

$$\begin{aligned} R_d^{T}(x,{\dot{x}},\xi ,{\dot{\xi }})&= \left( TR_d \circ \kappa _Q\right) (x,{\dot{x}},\xi ,{\dot{\xi }})=TR_d(x,\xi ; {\dot{x}},{\dot{\xi }})\\&=\left( x, \cos \left( \Vert \xi \Vert \right) x+\sin \left( \Vert \xi \Vert \right) \frac{\xi }{\Vert \xi \Vert }, {\dot{x}}, \cos (\Vert \xi \Vert ) \, {\dot{x}}+ N(x,\xi ) \, {\dot{\xi }}\right) \\&\equiv \left( x, {\dot{x}} ; \; \cos \left( \Vert \xi \Vert \right) x+\sin \left( \Vert \xi \Vert \right) \frac{\xi }{\Vert \xi \Vert }, \cos (\Vert \xi \Vert ) \, {\dot{x}}+ N(x,\xi ) \, {\dot{\xi }}\right) . \end{aligned}$$

To compute the cotangent lift of $R_d$, $R_d^{T^*}:TT^*Q\rightarrow T^*Q\times T^*Q$, we first need the tangent map of the inverse map $R^{-1}_d$ or, equivalently, the inverse of the tangent map:

$$\begin{aligned} D R^{-1}_d=\left( D R_d\right) ^{-1}=\begin{pmatrix} {\mathrm{Id}} &{} 0 \\ -\cos (\Vert \xi \Vert ) \, N^{-1} &{} N^{-1} \end{pmatrix}. \end{aligned}$$

Thus, the cotangent lift of $R_d^{-1}$ is given by:

$$\begin{aligned} \widehat{R_d}(x,\xi ,p_x,p_{\xi })&=\left( R_d(x,\xi ), (p_x,p_{\xi }) \left( D_{(x,\xi )} R_d\right) ^{-1}\right) . \end{aligned}$$

Finally, the discretization map on $T^*Q$ is obtained as follows:

$$\begin{aligned}&R_d^{T^*}(x,p,{\dot{x}},{\dot{p}})\\&\quad =\left( \Phi ^{-1}\circ \widehat{R_d} \circ \alpha _Q \right) (x,p,{\dot{x}},{\dot{p}})=\left( \Phi ^{-1}\circ \widehat{R_d}\right) (x,{\dot{x}},{\dot{p}},p)\\&\quad = \Phi ^{-1} \left( x, \cos \left( \Vert {\dot{x}}\Vert \right) x+\sin \left( \Vert {\dot{x}}\Vert \right) \frac{{\dot{x}}}{\Vert {\dot{x}}\Vert }\, ; \, {\dot{p}}-\cos \left( \Vert {\dot{x}}\Vert \right) \, p N^{-1}, p\, N^{-1}\right) \, \\&\quad \equiv \left( x,-{\dot{p}}+\cos \left( \Vert {\dot{x}}\Vert \right) \, p N^{-1}\, ; \, \cos \left( \Vert {\dot{x}}\Vert \right) x+\sin \left( \Vert {\dot{x}}\Vert \right) \frac{{\dot{x}}}{\Vert {\dot{x}}\Vert }\, , p\, N^{-1} \right) \,. \end{aligned}$$

$\square$

Example 3.3

Let us lift the discretization map in Example 2.3. As in the previous example, to simplify the computations we consider the discretization map $R_{d,{\mathrm{cay}}}:TSO(3)\rightarrow SO(3)\times SO(3)$ that fixes the first point (compared with (5)) as follows:

$$\begin{aligned} R_{d,{\mathrm{cay}}}(A,X)= \left( \!A, A\, {\mathrm{cay}} (A^T\,X)\!\right) =\left( \!A,A \left( \!{\mathrm{Id}}_3-A^TX/2\!\right) ^{-1}\left( {\mathrm{Id}}_3+A^TX/2\right) \!\right) . \end{aligned}$$

The tangent map of $R_{d,{\mathrm{cay}}}$ is given by the matrix

Thus,

$$\begin{aligned} R_{d,{\mathrm{cay}}}^{T}(A,{\dot{A}},X,{\dot{X}})&= \left( TR_{d,{\mathrm{cay}}} \circ \kappa _Q\right) (A,{\dot{A}},X,{\dot{X}})=TR_d(A,X,{\dot{A}},{\dot{X}})\\&=\left( A, A\, {\mathrm{cay}} (A^TX), {\dot{A}}, M {\dot{A}}+ N \,{\dot{X}}\right) \\&\equiv \left( A, {\dot{A}}; A \,{\mathrm{cay}} (A^TX), M {\dot{A}}+ N \,{\dot{X}}\right) \, , \end{aligned}$$

where N is an invertible matrix. To compute the discretization map $R_{d,{\mathrm{cay}}}^{T^*}:TT^*Q\rightarrow T^*Q\times T^*Q$ on $T^*Q$ as the cotangent lift of $R_d$, we first need the tangent map of the inverse map $R^{-1}_{d,{\mathrm{cay}}}$ or, equivalently, the inverse of the tangent map:

$$\begin{aligned} D R^{-1}_{d,{\mathrm{cay}}}=\left( D R_{d,{\mathrm{cay}}}\right) ^{-1}=\begin{pmatrix} {\mathrm{Id}} &{} 0 \\ -N^{-1}M &{} N^{-1} \end{pmatrix}. \end{aligned}$$

Thus, the cotangent lift of $R^{-1}_{d,{\mathrm{cay}}}$ is given by:

$$\begin{aligned} \widehat{R_{d,{\mathrm{cay}}}}(A,X,p_A,p_X)&=\left( R_{d,{\mathrm{cay}}}(A,X), (p_A,p_X) \left( D R_{d,{\mathrm{cay}}}\right) ^{-1}(A,X)\right) . \end{aligned}$$

Finally, the discretization map on $T^*Q$ is obtained as follows:

$$\begin{aligned} R_{d,{\mathrm{cay}}}^{T^*}(A,p,{\dot{A}},{\dot{p}})&=\left( \Phi ^{-1}\circ \widehat{R_{d,{\mathrm{cay}}}} \circ \alpha _Q \right) (A,p,{\dot{A}},{\dot{p}})\\&=\left( \Phi ^{-1}\circ \widehat{R_{d,{\mathrm{cay}}}}\right) (A,{\dot{A}},{\dot{p}},p)\\&= \Phi ^{-1} \left( A, A\, {\mathrm{cay}} (A^T{\dot{A}}), {\dot{p}}-pN^{-1}M,pN^{-1}\right) \\&= \left( A, pN^{-1}M-{\dot{p}};A \, {\mathrm{cay}} (A^T{\dot{A}}),pN^{-1}\right) \,. \end{aligned}$$

$\square$

4 Discretization Maps Associated to SODEs

The tangent lift of discretization maps defined in Sect. 3 appears naturally when geometrically designing discretizations of second-order differential equations (SODEs). Remember that a second-order differential equation is a vector field $\Gamma$ such that $\tau _{TQ}(\Gamma )=T\tau _{Q}(\Gamma )$. This implies that the vector field $\Gamma$ on TQ is a section of the second-order tangent bundle $T^{(2)}Q$, as described in [19]. Locally, if we take coordinates $(q^i)$ on Q and induced coordinates $(q^i, {\dot{q}}^i)$ on TQ, then

$$\begin{aligned} \Gamma ={\dot{q}}^i\frac{\partial }{\partial q^i}+\Gamma ^i(q, {\dot{q}})\frac{\partial }{\partial {\dot{q}}^i}\; . \end{aligned}$$

To find the integral curves of $\Gamma$ is equivalent to solve the following system of second-order differential equations:

$$\begin{aligned} \frac{d^2 q^i}{dt^2}=\Gamma ^i\left( q, \frac{dq}{dt}\right) \; . \end{aligned}$$

Now, we want to discretize these equations using the notion of discretization map defined on TQ as in Definition 2.2. Here we have two options: we could directly define a discretization map on TQ denoted by $R_d^{TTQ}:TTQ \rightarrow TQ\times TQ$ or we could tangently lift a discretization map on Q to obtain $R_d^T:TTQ \rightarrow TQ\times TQ$ as defined in Proposition 3.2.

Let us consider in general that we have a discretization map on TQ,

$$\begin{aligned} R_d^{TTQ}: TTQ\rightarrow TQ\times TQ\, , \end{aligned}$$

given by $R_d^{TTQ}(q, v, {\dot{q}}, {\dot{v}})= \left( \left( R^{TTQ}\right) ^1 (q, v, {\dot{q}}, {\dot{v}}), \left( R^{TTQ}\right) ^2 (q, v, {\dot{q}}, {\dot{v}})\right)$. Note that $\left( R^{TTQ}\right) ^i (q, v, {\dot{q}}, {\dot{v}})\in TQ$ for $i=1,2$.

A first option for discretizing a SODE $\Gamma$ consists of the following implicit discrete equation:

$$\begin{aligned} \left( \left( R^{TTQ}\right) ^2\circ h\Gamma \right) (q_k, v_k)=\left( \left( R^{TTQ}\right) ^1\circ h\Gamma \right) (q_{k+1}, v_{k+1})\, , \end{aligned}$$

(10)

where h is a positive small real number that determines the step size. The numerical method starts from the initial data $v_k\in T_{q_k} Q$, then the Eq. (10) is solved implicitly to obtain $v_{k+1}\in T_{q_{k+1}} Q$. Sect. 4.1 shows that a discretization map on TQ, not coming from a tangent lift, recovers Newmark method using the discretization method in Eq. (10). Geometrically, these methods given in Eq. (10) are based on the structure of groupoid of an implicit difference equation, in this case $TQ\times TQ\rightrightarrows TQ$ (see [29] for more details).

A second option for discretizing a SODE consists of the following numerical scheme:

$$\begin{aligned}&h\, \Gamma \left( \left( \tau _{TQ}\circ \left( R^{TTQ}_d\right) ^{-1}\right) (q_k, v_k; q_{k+1}, v_{k+1})\right) \nonumber \\&\quad =\left( R^{TTQ}_d\right) ^{-1}(q_k, v_k; q_{k+1}, v_{k+1})\,. \end{aligned}$$

(11)

As in Eq. (10), the numerical method is usually implicit. We will focus on this discretization process in Sects. 4.2 and 5 when constructing geometric integrators for mechanical systems.

Let us do a simple example to show that the numerical schemes in Eqs. (10) and (11) are usually different.

Example 4.1

Consider the discretization map on TQ obtained from a tangent lift in Example 3.1 and the inverse map:

$$\begin{aligned} R_d^T(q,{\dot{q}},v,{\dot{v}})&=\left( q-\dfrac{1}{2}\, v, {\dot{q}}-\dfrac{1}{2}\, {\dot{v}}\, ; q+\dfrac{1}{2}\, v\, , {\dot{q}}+\dfrac{1}{2}\, {\dot{v}}\right) \,,\\ \left( R_d^T\right) ^{-1}(q_0,v_0;q_1,v_1)&=\left( \dfrac{q_0+q_1}{2},\dfrac{v_0+v_1}{2}; q_1-q_0,v_1-v_0 \right) . \end{aligned}$$

The method in Eq. (10) becomes:

$$\begin{aligned} \frac{q_{k+1}-q_{k}}{h}&= \frac{v_k+v_{k+1}}{2}\, ,\\ \frac{v_{k+1}-v_{k}}{h}&= \frac{1}{2}\left( \Gamma (q_k, v_k) +\Gamma (q_{k+1}, v_{k+1})\right) . \end{aligned}$$

However, for the same discretization map on TQ the method in Eq. (11) is given by the following equations:

$$\begin{aligned} \frac{q_{k+1}-q_{k}}{h}&= \frac{v_k+v_{k+1}}{2}\, ,\\ \frac{v_{k+1}-v_{k}}{h}&= \Gamma \left( \frac{q_k+q_{k+1}}{2}, \frac{v_k+v_{k+1}}{2}\right) . \end{aligned}$$

$\square$

4.1 Newmark Method from a Discretization Map

An example of discretization using Eq. (10) is the Newmark method [46], a classical time-stepping method very common in structural mechanical codes. For simplicity, we consider a typical mechanical Lagrangian $L: T{{\mathbb {R}}}^n\longrightarrow {{\mathbb {R}}}$:

$$\begin{aligned} L(q, {\dot{q}})=\frac{1}{2}{\dot{q}} M{\dot{q}}^T-V(q)\, , \end{aligned}$$

where $(q, {\dot{q}})\in T{{\mathbb {R}}}^n$, M is a positive definite constant matrix and V is a potential function. The corresponding Euler–Lagrange equations are:

$$\begin{aligned} \ddot{q}=-M^{-1}\nabla V(q)\,, \end{aligned}$$

(12)

where $\nabla$ denotes the gradient of the potential function.

The Newmark methods are widely used in simulations of such mechanical systems, including even external forces [32]. To construct the method, two real parameters $\alpha$ and $\beta$ are selected so that the algorithm determines $(q_{k+1}, {\dot{q}}_{k+1})$ in terms of $(q_{k}, {\dot{q}}_{k})$ as follows:

$$\begin{aligned} q_{k+1}&=q_k+h{\dot{q}}_k+\frac{h^2}{2}\left( (1-2\beta ) a_{k}+2\beta a_{k+1}\right) \\ {\dot{q}}_{k+1}&={\dot{q}}_k+h\left( (1-\gamma ) a_k+\gamma a_{k+1}\right) \, , \nonumber \end{aligned}$$

(13)

where $a_k=-M^{-1}\nabla V(q_k)$ and $a_{k+1}=-M^{-1}\nabla V(q_{k+1})$.

We show here that the family of Newmark methods can be obtained from a discretization map on the tangent bundle TQ. Let us define $\left( R^{TTQ}_d\right) : TT{{\mathbb {R}}}^n\equiv {{\mathbb {R}}}^{4n}\rightarrow T{{\mathbb {R}}}^n\times T {{\mathbb {R}}}^n\equiv {{\mathbb {R}}}^{2n}\times {{\mathbb {R}}}^{2n}$ by

$$\begin{aligned} \left( R^{TTQ}\right) ^1(q, v, {\dot{q}}, {\dot{v}})&=\left( q-\frac{1}{2}{\dot{q}}+\frac{h}{2}(\gamma -2\beta ){\dot{v}}, v-\gamma {\dot{v}}\right) \, , \\ \left( R^{TTQ}\right) ^2(q, v, {\dot{q}}, {\dot{v}})&=\left( q+\frac{1}{2}{\dot{q}}+\frac{h}{2}(\gamma -2\beta ){\dot{v}}, v+(1-\gamma ) {\dot{v}}\right) . \end{aligned}$$

The Jacobian matrix of $R_d^{TTQ}$ is

$$\begin{aligned} \left( \begin{array}{rrrr} {\mathrm{Id}} &{}0&{}-\frac{1}{2}\,{\mathrm{Id}} &{}\frac{h}{2}(\gamma -2\beta ) \,{\mathrm{Id}}\\ 0&{} {\mathrm{Id}}&{}0&{}-\gamma \,{\mathrm{Id}}\\ {\mathrm{Id}} &{}0&{}\frac{1}{2} \,{\mathrm{Id}}&{}\frac{h}{2}(\gamma -2\beta ) \,{\mathrm{Id}}\\ 0&{}{\mathrm{Id}}&{}0&{}(1-\gamma ) \,{\mathrm{Id}} \end{array} \right) . \end{aligned}$$

It is straightforward that $R_d^{TTQ}$ satisfies both properties in Definition 2.2. Hence, $R_d^{TTQ}$ is a discretization map on TQ.

The Euler–Lagrange equations (12) can be rewritten as the submanifold S of $T^{(2)}Q\subset TTQ$,

$$\begin{aligned} S=\{(q, {\dot{q}}, a ) \; \mid \; a=-M^{-1}\nabla V(q) \}\, , \end{aligned}$$

with the natural inclusion $i: T^{(2)}Q\hookrightarrow TTQ$, $i(q, {\dot{q}}, a )=(q, {\dot{q}}, {\dot{q}}, a )$.

Hence, the dynamics induced by the Newmark method is equivalent to the following algorithm:

• 1. Take an initial position and velocity $(q_k, {\dot{q}}_k)$.

• 2. Evaluate $a_k=-M^{-1}\nabla V(q_k)$.

• 3. Solve the system obtained from Eq. (10):

  $$\begin{aligned} \left( R^{TTQ}\right) ^2(q_k, {\dot{q}}_k; h{\dot{q}}_k, ha_k)= \left( R^{TTQ}\right) ^1(q_{k+1}, {\dot{q}}_{k+1}; h{\dot{q}}_{k+1}, h a_{k+1})\, , \end{aligned}$$

  (14)

  where $a_{k+1}=-M^{-1}\nabla V(q_{k+1})$.

Observe that Eq. (14) is equal to

$$\begin{aligned} q_k+\frac{h}{2}{\dot{q}}_k+\frac{h^2}{2}(\gamma -2\beta )a_k= & {} q_{k+1}-\frac{h}{2}{\dot{q}}_{k+1}-\frac{h^2}{2}(\gamma -2\beta )a_{k+1}\, ,\\ {\dot{q}}_k+h(1-\gamma )a_k= & {} {\dot{q}}_{k+1}-h\gamma a_{k+1}. \end{aligned}$$

After algebraic manipulations, the above equations are equivalent to the well-known Newmark method in Eq. (13).

Note that if $\gamma =1/2$ and $\beta =1/4$, then $R_d^{TTQ}$ is precisely the tangent lift of the discretization map on Q coming from the mid-point rule as described in Example 3.1.

4.2 Discretization Maps Associated with Discrete Second-Order Equations

In this section we briefly discuss the possibility to find a discrete version of a second-order differential equation (SODE) using a second-order discrete equation (SOdE).

According to [41], a SOdE is given by a map $\Gamma _d: Q\times Q\rightarrow Q\times Q\times Q\times Q$ such that

$$\begin{aligned} \Gamma _d(q_{k-1}, q_k)=(q_{k-1}, q_k, q_k, {\tilde{\Gamma }}_d(q_{k-1}, q_k))\,, \end{aligned}$$

in other words, $q_{k+1}={\tilde{\Gamma }}_d(q_{k-1}, q_k)$. From two initial conditions $q_0$, $q_1$ this equation defines the discrete evolution as the sequence $\{q_0, q_1, q_2, \ldots \}$.

Given a discretization map on Q and a second-order vector field $\Gamma$, we wonder if, under any assumption, the tangent lift of the discretization map, $R_d^T$, could define a discrete second-order equation $\Gamma _d$. The specific question is: When does a discretization map $R_d$ make Diagram (15) commutative?

(15)

Proposition 4.1

Let Q be a vector space and $\Gamma$ be a SODE. If $R_d:TQ \rightarrow Q\times Q$ is the discretization map on Q defined from the $\theta$-method:

$$\begin{aligned} R_d(q,v)=(q-\theta \, v, q+ (1-\theta )\, v), \end{aligned}$$

(16)

then Diagram (15) is commutative, that is,

$$\begin{aligned} (R_d,R_d)\circ R_d^{T}\circ \Gamma \circ R_ d^{-1}:Q\times Q \rightarrow Q\times Q \times Q \times Q \end{aligned}$$

defines a second-order discrete equation (SOdE).

Proof

Let (q, v) local coordinates for TQ. Let $R_d(q,v)=\left( R^1_d(q,v),R^2_d(q,v)\right)$, we compute

$$\begin{aligned} \left( R_d^{T}\circ \Gamma \right) (q,v)=\left( T_{(q,v)}R^1_d (\Gamma (q,v)), T_{(q,v)}R^2_d (\Gamma (q,v))\right) . \end{aligned}$$

If we apply now $(R_d,R_d):TQ\times TQ \rightarrow Q\times Q \times Q\times Q$, the resulting expression defines a SOdE if and only if the second and third component are equal, that is,

$$\begin{aligned} R^2_d(T_{(q,v)}R^1_d (\Gamma (q,v)))=R^1_d(T_{(q,v)}R^2_d (\Gamma (q,v))). \end{aligned}$$

(17)

It is a straightforward computation to verify that the discretization maps defined from the $\theta$-method in (16) satisfy Eq. (17). $\square$

In fact, the above proposition could be stated more generally. Any discretization map that satisfies Eq. (17) defines a SOdE by the tangent lift of that map.

Equation (17) is equivalent to the commutativity of the following diagram:

(18)

In particular, if the discretization map $R_d: TQ\rightarrow Q\times Q$ is defined from a standard retraction map as in Definition 2.1, that is, $R_d^1=\tau _Q$, then Diagram (18) is always commutative

since $(\tau _Q\circ TR_d^2)(\Gamma )= (R_d^2\circ T\tau _Q)(\Gamma )$. Therefore, any standard retraction map defines a SOdE $\Gamma _d$.

5 Construction of Geometric Integrators from Discretization Maps

In this section we describe how geometric integrators are obtained for both Hamiltonian and Euler–Lagrange equations by discretizing their equations using discretization maps. In Sect. 5.3, we establish the relation with discrete variational calculus where the variational principles are discretized to obtain the discrete flow (see [41]).

In Sect. 5.1 we look at the Hamiltonian framework [1]. Hamiltonian systems have the property that the associated flow is a symplectic transformation. As described in [7, 27, 47], it is important to define numerical methods that also preserve that property. Remember that a numerical one-step method is called symplectic if the one-step map, in other words, the discrete flow, is symplectic whenever the method is applied to a smooth Hamiltonian system.

Second, we describe geometric integrators obtained from the Lagrangian viewpoint in Sect. 5.2.

In order to describe Hamiltonian and Lagrangian mechanics, we consider the symplectic manifold $(T^*Q,\omega _Q)$ that has the musical isomorphisms $\omega _Q^\flat : {{\mathfrak {X}}}(T^*Q)\rightarrow \Omega ^1 (T^*Q)$ defined by $\omega _Q^\flat (X)=\alpha$ where $i_X\omega _Q=\alpha$ (see, for instance, [37]). The inverse of $\omega _Q^\flat$ is denoted by $\omega _Q^\sharp$, that is, $\omega _Q^\sharp =(\omega _Q^\flat )^{-1}$.

5.1 Geometric Integrators in Hamiltonian Framework

Let $H: T^*Q\rightarrow {{\mathbb {R}}}$ be a Hamiltonian function with corresponding Hamiltonian vector field $X_H$ derived from Hamilton’s equations:

$$\begin{aligned} i_{X_H}\omega _Q=dH. \end{aligned}$$

The triple $(T^*Q,\omega _Q,H)$ defines a Hamiltonian system. Equivalently, an integral curve of $X_H$ is solution to Hamilton’s equations:

$$\begin{aligned} \frac{dq^i}{dt}=\frac{\partial H}{\partial p_i}\; ,\qquad \frac{dp_i}{dt}=-\frac{\partial H}{\partial q^i}\, , \end{aligned}$$

where $(q^i, p_i)$ are canonical coordinates on $T^*Q$ (see [1]). In other words, a solution $\gamma :I \rightarrow T^*Q$ of Hamilton’s equations must satisfy

$$\begin{aligned} \omega _Q^\flat \left( {\dot{\gamma }}(t)\right) ={\mathrm{d}} H(\gamma (t))\,, \text{ equivalently } {\dot{\gamma }}(t) =\omega ^\sharp \left( {\mathrm{d}} H(\gamma (t))\right) . \end{aligned}$$

A discretization map on $T^*Q$, that is, $R^{TT^*Q}_d:TT^*Q\rightarrow T^*Q\times T^*Q$ defines the following numerical integrator for step size h:

$$\begin{aligned}&\left( R^{TT^*Q}_d\right) ^{-1}(q_0, p_0; q_1, p_1)\nonumber \\&\quad =\omega ^\sharp \left( h\, {\mathrm{d}}H\left( \left( \tau _{T^*Q}\circ \left( R^{TT^*Q}_d\right) ^{-1}\right) (q_0, p_0; q_1, p_1)\right) \right) . \end{aligned}$$

(19)

Equivalently, similar to Eq. (11), we have

$$\begin{aligned} h\, X_H \left( \left( \tau _{T^*Q}\circ \left( R^{TT^*Q}_d\right) ^{-1}\right) (q_0, p_0; q_1, p_1)\right) = \left( R^{TT^*Q}_d\right) ^{-1}(q_0, p_0; q_1, p_1) .\nonumber \\ \end{aligned}$$

(20)

This numerical integrator may be defined for any discretization map on $T^*Q$. However, if such a map is the cotangent lift of a discretization map on Q (see Sect. 3.2), then the numerical integrator is symplectic as stated in the following proposition.

For the proof we need to recall the notion of a Lagrangian submanifold of a symplectic manifold $(M, \omega )$.

An immersed submanifold N of M, or immersion, $f : N \rightarrow M$ is Lagrangian if so is the space $Tf (T_x N)$ as a subspace of $T_{f(x)} M$ for each point $x \in N$, that is, $Tf (T_x N)=(Tf (T_x N))^\perp$ where $^\perp$ denotes the orthogonal complement of the subspace with respect to the symplectic form. Note that an immersion $f:N\rightarrow M$ is Lagrangian if and only if $f^*\omega =0$ and the dimension of N is half the dimension of M. The most common way to define a Lagrangian submanifold of a symplectic manifold is as graph of a closed one-form.

Proposition 5.1

Let $R_d:TQ\rightarrow Q\times Q$ be a discretization map on Q and $H: T^*Q\rightarrow {{\mathbb {R}}}$ be a Hamiltonian function. Equation (19) written for the cotangent lift of $R_d$, that is, $R_d^{T^*}$, defines a symplectic integrator of the Hamiltonian system $(T^*Q,\omega _Q, H)$.

Proof

As the submanifold ${\mathrm{d}}H\left( \left( \tau _{T^*Q}\circ \left( R_d^{T^*}\right) ^{-1}\right) (q_0, p_0; q_1, p_1)\right)$ in Eq. (19) is Lagrangian on $(T^*T^*Q, \omega _{T^*Q})$, Proposition 3.5 guarantees that Eq. (19) determines a Lagrangian submanifold of $(T^*Q\times T^*Q, \Omega _{12})$. Locally, such a manifold can be expressed as the graph of a local symplectomorphism $\varphi _h:T^*Q\rightarrow T^*Q$, see [37] for more details. Consequently, the geometric method obtained from Eq. (19) is symplectic. $\square$

Let us use the above result to obtain some of the symplectic numerical methods known in the literature.

Example 5.1

Let $H:T^*Q\rightarrow {\mathbb {R}}$ be a Hamiltonian function, the cotangent lift (8), (9) of the discretization map associated to the mid-point rule in Example 3.1 used in Eq. (19) leads to the following equations:

$$\begin{aligned}&\left( \frac{q_0+q_1}{2}, \right. \left. \frac{p_0+p_1}{2}, q_1-q_0, p_1-p_0\right) \\&\quad =\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}, \right. \\&\qquad \left. h\frac{\partial H}{\partial p}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) , -h\frac{\partial H}{\partial q}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) \right) \; . \end{aligned}$$

Equivalently, the equations describe the following symplectic integrator:

$$\begin{aligned} \frac{q_1-q_0}{h}= & {} \frac{\partial H}{\partial p}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) \, ,\\ \frac{p_1-p_0}{h}= & {} -\frac{\partial H}{\partial q}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) . \end{aligned}$$

The above integrator corresponds with an implicit second-order symplectic method with initial condition $(q_0,p_0)$. $\square$

Example 5.2

The discretization map on Q, $R_d(q,v)=(q-v,q)\,,$ is lifted to the cotangent bundle as follows

$$\begin{aligned} \begin{array}{cccc} R_d^{T^*}:&{} TT^*Q &{} \longrightarrow &{} T^*Q \times T^*Q \\ &{} (q,p,{\dot{q}},{\dot{p}}) &{}\longmapsto &{} (q-{\dot{q}}, p, q,p+{\dot{p}}).\end{array} \end{aligned}$$

As $\left( R_d^{T^*}\right) ^{-1}(q_0,p_0,q_1,p_1)=(q_1,p_0,q_1-q_0,p_1-p_0)$, Eq. (19) leads to the following symplectic method:

$$\begin{aligned} \frac{q_1-q_0}{h}= & {} \frac{\partial H}{\partial p}\left( q_1,p_0\right) \, ,\\ \frac{p_1-p_0}{h}= & {} -\frac{\partial H}{\partial q}\left( q_1,p_0\right) . \end{aligned}$$

For a Hamiltonian function $H(p,q)=\dfrac{1}{2}p M p^T+V(q)$, with a constant positive definite matrix M, the integrator is an explicit symplectic method. $\square$

Example 5.3

Now consider a Hamiltonian function $H: T^*S^2\rightarrow {{\mathbb {R}}}$ on $T^*S^2$ that we identify with the tangent bundle $T S^2$:

$$\begin{aligned} T^*S^2\equiv \{(x, p)\in {{\mathbb {R}}}^3\times {{\mathbb {R}}}^3\; \mid \; \left\Vert x\right\Vert =1, \; x\cdot p=0\}. \end{aligned}$$

For the discretization of the corresponding Hamiltonian equations, we will use the discretization map $R_d: TS^2\rightarrow S^2\times S^2$ given by

$$\begin{aligned} R_d(x, \xi )=\left( x, \frac{x+\xi }{\left\Vert x+\xi \right\Vert }\right) , \end{aligned}$$

whose inverse is precisely:

$$\begin{aligned} R_d^{-1}(x_0, x_1)=\left( x_0, \frac{x_1}{x_0\cdot x_1} -x_0 \right) \, , \end{aligned}$$

whenever it is well defined. Now, we will compute the inverse of the cotangent lift of the discretization map given in Eq. (7), that is,

$$\begin{aligned} \left( R_d^{T^*}\right) ^{-1}(x_0,p_0; x_1,p_1)= \alpha _Q^{-1}\left( R_d^{-1} (x_0,x_1), (-p_0, p_1)\, { D}_{R_d^{-1} (x_0,x_1)}R_d\right) . \end{aligned}$$

Having in mind the definition of $T^*S^2$, it can be computed that the matrix ${D}_{R_d^{-1} (x,y)}R_d$ is equal to:

$$\begin{aligned} {D}_{R_d^{-1} (x,y)}R_d = \left( \begin{array}{rr} {\mathrm{Id}}_{3\times 3}&{}0\\ (x\cdot y)\, {\mathrm{Id}}_{3\times 3}&{} C \end{array} \right) \end{aligned}$$

where C is the matrix with entries

$$\begin{aligned} c_{ij}=\left\{ \begin{array}{ll} (x\cdot y)\left[ 1+(x\cdot y)y_ix_i-y_i^2\right] &{}\quad \text{ if } i=j, \\ (x\cdot y)\left[ (x\cdot y)y_ix_j-y_iy_j\right] &{}\quad \text{ if } i\not =j. \end{array}\right. \end{aligned}$$

Therefore,

$$\begin{aligned} \left( R_d^{T^*}\right) ^{-1}(x_0,p_0; x_1,p_1)= \left( x_0, p_1 C\, ; \, \frac{1}{x_0\cdot x_1} x_1-x_0, -p_0+(x_0\cdot x_1)p_1\right) . \end{aligned}$$

As a result, we obtain the following symplectic integrator for Hamilton’s equations:

$$\begin{aligned} \frac{1}{x_k\cdot x_{k+1}} x_{k+1}-x_k&=h\frac{\partial H}{\partial p}(x_k, p_{k+1}C)\, ,\\ -p_k+(x_k\cdot x_{k+1})p_{k+1}&=-h\frac{\partial H}{\partial q}(x_k, p_{k+1}C) . \end{aligned}$$

$\square$

Remark 5.1

Another option to construct geometric integrators is to use an expression similar to Eq. (10) but now adapted for Hamiltonian vector fields, that is,

$$\begin{aligned} \left( \left( R_d^{TT^*Q}\right) ^2\circ h X_H\right) (q_k, p_k)=\left( \left( R_d^{TT^*Q}\right) ^1\circ h X_H\right) (q_{k+1}, p_{k+1})\; . \end{aligned}$$

Note that here the discretization map on $T^*Q$ does not have to be the cotangent lift of one on Q. However, even if the cotangent lift is considered, the method is not necessarily symplectic. For instance, for the discretization map coming from the mid-point rule in Examples 2.1 and 3.1, we obtain the symmetric second-order method:

$$\begin{aligned} \frac{q_{k+1}-q_k}{h}&= \frac{1}{2}\left( \frac{\partial H}{\partial p}(q_k, p_k)+\frac{\partial H}{\partial p}(q_{k+1}, p_{k+1}) \right) \; ,\\ \frac{p_{k+1}-p_k}{h}&= -\frac{1}{2}\left( \frac{\partial H}{\partial q}(q_k, p_k)+\frac{\partial H}{\partial q}(q_{k+1}, p_{k+1}) \right) \; .\\ \end{aligned}$$

However, this method is not symplectic because, in general, $dq_{k+1}\wedge d p_{k+1}-dq_{k}\wedge d p_{k}\not = 0$ when restricted to the numerical scheme. $\square$

Remark 5.2

Observe that our method gives us a constructive way to derive symplectic integrators for Hamiltonian systems. It will be interesting to compare our methods with other previous approaches ([36]), specially when the configuration space is a Lie group and we can use well-known retraction maps such as exponential maps and other approximations. See [8, 10, 12, 31]. $\square$

5.2 Geometric Integrators in Lagrangian Framework

Let us consider a regular Lagrangian function $L: TQ\rightarrow {{\mathbb {R}}}$ so that there exists a second-order vector field $\Gamma _L$ on TQ and Euler–Lagrange equations are given by

$$\begin{aligned} {\mathrm{i }}_{\Gamma _L}\Omega _L={\mathrm{d}} E_L, \end{aligned}$$

where $E_L$ is the energy function and $\Omega _L$ is the symplectic Lagrange 2-form obtained by the pull-back of the Legendre map ${{\mathcal {F}}}L:TQ \rightarrow T^*Q$ of the natural symplectic form on $T^*Q$, that is, $\Omega _L=({{\mathcal {F}}}L)^*\omega _Q$ (see [1] for more details).

As in Eq. (11), a discretization map on TQ, that is, $R^{TTQ}_d:TTQ \rightarrow TQ\times TQ$, defines the following numerical integrator:

$$\begin{aligned} R^{TTQ}_d\left( h\, \Gamma _L\left( \left( \tau _{TQ}\circ \left( R_d^{TTQ}\right) ^{-1}\right) (q_0,v_0;q_1,v_1)\right) \right) =(q_0,v_0;q_1,v_1). \end{aligned}$$

(21)

Equivalently,

$$\begin{aligned} h\, \Gamma _L\left( \left( \tau _{TQ}\circ \left( R_d^{TTQ}\right) ^{-1}\right) (q_0,v_0;q_1,v_1) \right) =\left( R_d^{TTQ}\right) ^{-1}(q_0,v_0;q_1,v_1). \end{aligned}$$

As the Lagrangian function is regular, we could move to the Hamiltonian framework and construct a symplectic numerical integrator using the cotangent lift of a discretization map on Q as in Proposition 5.1. It remains to prove if the obtained numerical integrator is a discretization map on TQ as described in Definition 2.2.

Remember that the manifold $(TQ\times TQ, \Omega _L^1-\Omega _L^0)$ is symplectic. Locally, the symplectic 2-form is given by $\Omega _L^1-\Omega _L^0=({{\mathcal {F}}}L,{{\mathcal {F}}}L)^*( {\mathrm{d}}q_i^1\wedge {\mathrm{d}} p_i^1\,- {\mathrm{d}}q_i^0\wedge {\mathrm{d}} p_i^0)$.

Proposition 5.2

Let $R_d:TQ\rightarrow Q\times Q$ be a discretization map on Q and $L: TQ\rightarrow {{\mathbb {R}}}$ be a regular Lagrangian function. The two following facts are satisfied:

• (a) the map $R_d^{L}=({{\mathcal {F}}}L^{-1},{{\mathcal {F}}}L^{-1})\circ R_d^{T^*}\circ T{{\mathcal {F}}}L:TTQ \rightarrow TQ\times TQ$ defines a symplectic numerical integrator of the Euler–Lagrange equations for L;

• (b) the above-mentioned map $R_d^{L}$ is a discretization map on TQ.

Proof

First, we prove property (a). As the Lagrangian function is regular, the Legendre map is a local diffeomorphism. Propositions 3.1 and 5.1 guarantee that $R_d^L$ is a symplectomorphism because it is a composition of symplectomorphisms. Hence, $R_d^{L}$ defines a symplectic numerical integrator in Eq. (21).

The diagram below shows the constructive process for $R_d^{L}$:

In other words, the Lagrangian submanifold ${\mathrm{Im}}\, \Gamma _L$ of $\left( TTQ,\left( T{{\mathcal {F}}}L\right) ^*({\mathrm{d}}_T \omega _Q)\right)$ is preserved by $R_d^{L}$ and the numerical method in Eq. (21) is symplectic.

For (b), we must prove the properties in Definition 2.2 for the map $R^L_d :TTQ \rightarrow TQ\times TQ$.

• 1. Note that $R^L_d(v_q,0_{v_q})=(v_q,v_q)$ because

  $$\begin{aligned} R^L_d(v_q,0_{v_q})&=\left( ({{\mathcal {F}}}L^{-1},{{\mathcal {F}}}L^{-1})\circ R_d^{T^*}\right) \left( {{\mathcal {F}}}L(v_q); 0_{{{\mathcal {F}}}L(v_q)}\right) \\&=({{\mathcal {F}}}L^{-1},{{\mathcal {F}}}L^{-1})\left( {{\mathcal {F}}}L(v_q); {{\mathcal {F}}}L(v_q)\right) \\&=(v_q,v_q). \end{aligned}$$

  The second equality is true because $R_d^{T^*}$ is a discretization map on $T^*Q$ as shown in Proposition 3.4.

• 2. We must prove that $T_{(q,v,0,0)}\left( R^L_d\right) _{(q,v)}^2-T_{(q,v,0,0)}\left( R^L_d\right) _{(q,v)}^1$ is the identity map from $T_{(q,v,0,0)}TTQ\simeq T_{(q,v)}TQ$ to itself.

  Let us first compute it for $i=1,2$:

  $$\begin{aligned}&\left. \dfrac{\mathrm{d}}{{\mathrm{d}}\, t}\right| _{t=0}\, {{\mathcal {F}}}L^{-1}\left( \left( R_d^{T^*} \right) ^i \left( T {{\mathcal {F}}} L\right) (q,v,t\, {\dot{q}},t\, {\dot{v}})\right) \\&\quad = D {{\mathcal {F}}}L^{-1}_{(q,\frac{\partial L}{\partial v})} T_{(q,\frac{\partial L}{\partial v},0,0)} \left( R_d^{T^*}\right) ^i D_{3,4}\left( T {{\mathcal {F}}} L\right) _{(q,v,0,0)} \end{aligned}$$

  Note that $D_{3,4}\left( T {{\mathcal {F}}} L\right) _{(q,v,0,0)}$ is the fiber derivative of the tangent map $T{{\mathcal {F}}} L$. Knowing that the tangent map is linear on the fiber, together with the fact that $R_d^{T^*}$ is a discretization map on $T^*Q$ and it satisfies the second property in Definition 2.2, we can conclude that the map $R^L_d$ satisfies the second property for being a discretization map on TQ.$\square$

It can be proved that only for a very specific discretization map $R_d$ on Q and Lagrangian function, the discretization map $R^L_d$ of TTQ is the tangent lift of $R_d$.

Corollary 5.3

Let $L(q,v)=\dfrac{1}{2}\, v^T \, M v-V(q)$ be the Lagrangian function, being M a positive-definite symmetric constant mass matrix and V the potential function. If $R_d:TQ\rightarrow Q \times Q$ is the discretization map on Q given by the mid-point rule in Example 2.1, then $R^L_d$ is the tangent lift of $R_d$.

Proof

The Legendre transformation for L in the corollary is ${{\mathcal {F}}}L(q,v)=(q,Mv)$ and the inverse map is ${{\mathcal {F}}}L^{-1}(q,p)=(q,M^{-1}p)$. Thus,

$$\begin{aligned} R^L_d(q,v,{\dot{q}},{\dot{v}})=\left( R^1_d(q,{\dot{q}}), -({\dot{v}},v){\mathrm{D}}R_d^{-1}(q,{\dot{q}})_{*,1};R^2_d(q,{\dot{q}}), ({\dot{v}},v){\mathrm{D}}R_d^{-1}(q,{\dot{q}})_{*,2} \right) \end{aligned}$$

where $A_{*,i}$ denotes the ith column of the matrix A. This expression is equal to the tangent lift $R_d^T$ if

$$\begin{aligned} -({\dot{v}},v){\mathrm{D}}R_d^{-1}(q,{\dot{q}})_{*,1}&= D_{(q,v)}R_d^1(q,v) ({\dot{q}},\;{\dot{v}})^T\\ ({\dot{v}},v){\mathrm{D}}R_d^{-1}(q,{\dot{q}})_{*,2}&=D_{(q,v)}R_d^2(q,v) ({\dot{q}},\;{\dot{v}})^T. \end{aligned}$$

Both equalities are satisfied if the discretization map $R_d$ is given by the mid-point rule, see Example 2.1. $\square$

Example 5.4

Let us consider a Lagrangian second-order vector field given by $\ddot{q}=-M^{-1} \, \nabla V(q)$. The numerical method in Eq. (21) for the mid-point rule described in Example 3.1 becomes:

$$\begin{aligned}&R^{T}_d\left( h\, \Gamma _L\left( \dfrac{q_0+q_1}{2}, \dfrac{v_0+v_1}{2} \right) \right) =(q_0,v_0,q_1,v_1) \\&R^{T}_d\left( \dfrac{q_0+q_1}{2}, \dfrac{v_0+v_1}{2} , h(q_1-q_0), -h\, M^{-1}\, \nabla V\left( \dfrac{q_0+q_1}{2} \right) \right) =(q_0,v_0,q_1,v_1). \end{aligned}$$

Given $(q_0,v_0)$, the numerical integrator is defined implicitly by

$$\begin{aligned} \dfrac{v_0+v_1}{2}= & {} \dfrac{q_1-q_0}{h}\; ,\\ \dfrac{v_1-v_0}{h}= & {} -M^{-1}\, \nabla V\left( \dfrac{q_0+q_1}{2}\right) . \end{aligned}$$

After some straightforward computations, we obtain that this discrete method is rewritten as an implicit second-order discrete equation given by:

$$\begin{aligned} \dfrac{q_2-2q_1+q_0}{h^2}=-\dfrac{1}{2}M^{-1}\left( \nabla V\left( \dfrac{q_0+q_1}{2}\right) + V\left( \dfrac{q_1+q_2}{2}\right) \right) \end{aligned}$$

In the next subsection we will explore the relation of these methods with discrete variational calculus. $\square$

5.3 Discrete Variational Calculus

As mentioned in [41], a usual way to design symplectic integrators from a Lagrangian system consists of discretizing the variational principle using a discrete Lagrangian map. Many of these discrete maps are obtained from a continuous Lagrangian map and a discretization map on Q by discretizing the continuous action as follows

$$\begin{aligned} {{\mathcal {S}}}(q_0, q_1)=\int _0^h L(q(t), {\dot{q}}(t))\; dt\approx h L\left( \frac{1}{h}R_d^{-1}(q_0, q_1)\right) =L^h_d(q_0,q_1)\, , \end{aligned}$$

where q(t) is the unique solution of the Euler–Lagrange equations such that $q(0)=q_0$ and $q(h)=q_1$ with h enough small. Observe that if $R_d^{-1}(q_0, q_1)=v_q\in T_qQ$ then $\frac{1}{h}R_d^{-1}(q_0, q_1)=\frac{1}{h}v_q\in T_qQ$. Therefore, the discrete Lagrangian $L^h_d: Q\times Q\rightarrow {{\mathbb {R}}}$ is defined by $L^h_d=h\left( L\circ \frac{1}{h}R_d^{-1}\right)$.

If we consider the Hamiltonian function $H(p,q)=\langle p, {\dot{q}} \rangle -L(q,{\dot{q}})$, then we can simultaneously consider the discretization of both the Lagrangian and Hamiltonian framework. The following diagram is commutative by construction (see for the left-hand side of the diagram):

In this diagram we understand that the multiplication by 1/h in $\frac{1}{h}R^{-1}_d$ is with respect to the vector bundle structure given by $\tau _Q: TQ\rightarrow Q$, $\frac{1}{h}\widehat{R_d}^{-1}$ with respect to the vector bundle structure given by $\pi _{TQ}: T^*TQ\rightarrow TQ$ and $\frac{1}{h}(R_d^{T^*})^{-1}$ with respect to $\tau _{T^*Q}: TT^*Q\rightarrow T^*Q$.

Therefore, we have that

$$\begin{aligned} R_d^{T^*}\left( hX_H\left( {{\mathcal {F}}} L \left( \frac{1}{h}R_d^{-1}(q_k, q_{k+1})\right) \right) \right)&=\Phi ^{-1}\widehat{R_d}\left( h {\mathrm{d}} L\left( \frac{1}{h}( R_d^{-1})(q_k,q_{k+1})\right) \right) \nonumber \\&=\Phi ^{-1} \left( {\mathrm{d}}\, L^h_d (q_k, q_{k+1})\right) . \end{aligned}$$

(22)

Using the previous equation and (20), we obtain

$$\begin{aligned} \Phi ^{-1} \left( {\mathrm{d}}\, L^h_d (q_k, q_{k+1})\right) = (q_k, p_{k}; q_{k+1}, p_{k+1}) \end{aligned}$$

(23)

and the discrete variational equations in [41] are recovered:

$$\begin{aligned} p_k= & {} -D_1L^h_d(q_k, q_{k+1})\, ,\\ p_{k+1}= & {} D_2L^h_d(q_k, q_{k+1}). \end{aligned}$$

These equations lead to the well-known discrete Euler–Lagrange equations:

$$\begin{aligned} D_1L^h_d(q_k,q_{k+1})+D_2L^h_d(q_{k-1},q_k)=0. \end{aligned}$$

Example 5.5

Given a regular Lagrangian $L: TQ\rightarrow {{\mathbb {R}}}$, where Q is a vector space, consider the discrete Lagrangian $L^h_d(q_0, q_1)=h L\left( \frac{q_0+q_1}{2}, \frac{q_1-q_0}{h}\right)$. Using Example 3.1, the left-hand side of Eq. (22) and the right-hand side of (23), we obtain the following integrator:

$$\begin{aligned} \frac{p_{k+1}-p_k}{h}= & {} \frac{\partial L}{\partial q}\left( \frac{q_k+q_{k+1}}{2}, \frac{q_{k+1}-q_k}{h}\right) ,\\ \frac{p_k+p_{k+1}}{2}= & {} \frac{\partial L}{\partial {\dot{q}}}\left( \frac{q_k+q_{k+1}}{2}, \frac{q_{k+1}-q_k}{h}\right) . \\ \end{aligned}$$

$\square$

6 Composition of Geometric Integrators

The construction of symplectic integrators based on discretization maps is closely related to the notion of Lagrangian submanifolds, as already appears in Sect. 5. For instance, Equation (20) defines the following Lagrangian submanifold of the symplectic manifold $(T^*Q\times T^*Q, \Omega _{12})$

$$\begin{aligned} {{\mathcal {L}}}^h=\left\{ (\alpha _q, \beta _{q'})\in T^*Q\times T^*Q\; \mid \; (\alpha _q, \beta _{q'})=R_d^{T^*}( h\, X_H \left( \gamma _{q''}\right) ) \right\} \end{aligned}$$

(24)

where $\gamma _{q''}=(\tau _{T^*Q}\circ \left( R^{T^*}_d\right) ^{-1})(\alpha _q, \beta _{q'})\in T^*Q$.

Now, we will use some well-known properties of Lagrangian submanifolds as the composition of Lagrangian submanifolds (see [25] for more details) to describe a particularly elegant method to construct high-order methods from a given low-order integrator (see [27]). To be more precise, we are going to geometrically describe the composition of two (or more) geometric integrators defined by different discretization maps. As a particular example, we will recover the well-known Störmer–Verlet method, a second-order symplectic method.

Let $R_{d,1}$ and $R_{d,2}:TQ \rightarrow Q\times Q$ be two discretization maps on Q and $H: T^*Q\rightarrow {{\mathbb {R}}}$ be a Hamiltonian function, using Equation (24) we define two Lagrangian submanifolds of $(T^*Q\times T^*Q, \Omega _{12})$ as follows:

$$\begin{aligned} {{\mathcal {L}}}^{h/2}_1&=\bigg \{ (q_k, p_k;q_{k+1/}, p_{k+1/2})\in T^*Q\times T^*Q\; \mid \; \exists \; \gamma _{k, k+{1/2}}\in T^*Q \text{ s. } \text{ t. } \\&\qquad (q_k, p_k ; q_{k+1/2}, p_{k+1/2})=R_{d,1}^{T^*}\left( \frac{h}{2}\, X_H \left( \gamma _{k, k+{1/2}}\right) \right) \bigg \}\, ,\\ {{\mathcal {L}}}^{h/2}_2&= \bigg \{ (q_{k+1/2}, p_{k+1/2}; q_{k+1}, p_{k+1})\in T^*Q\times T^*Q\; \mid \; \exists \; \gamma _{k+{1/2}, k+1}\in T^*Q \text{ s. } \text{ t. } \\&\qquad (q_{k+1/2}, p_{k+1/2}; q_{k+1}, p_{k+1})=R_{d,2}^{T^*}\left( \frac{h}{2}\, X_H \left( \gamma _{ k+{1/2}, k+1}\right) \right) \bigg \} \, , \end{aligned}$$

where

$$\begin{aligned} \gamma _{k, k+{1/2}}&=\left( \tau _{T^*Q}\circ \left( R^{T^*}_{d,1}\right) ^{-1}\right) (q_k, p_k; q_{k+1/2}, p_{k+1/2})\in T^*Q\, , \\ \gamma _{ k+{1/2}, k+1}&=\left( \tau _{T^*Q}\circ \left( R^{T^*}_{d,2}\right) ^{-1}\right) (q_k, p_k; q_{k+1/2}, p_{k+1/2})\in T^*Q. \end{aligned}$$

Under the assumption of clean intersection (see [26]), we compose the above two Lagrangian submanifolds as follows

$$\begin{aligned} {{\mathcal {L}}}^{h/2}_2\circ {{\mathcal {L}}}^{h/2}_1&= \left\{ (\alpha _q,\beta _{q''})\in T^*Q\times T^*Q\; \mid \, \exists \; \gamma _{q'}\in T^*Q \hbox { with } (\alpha _q,\gamma _{q'})\in {{\mathcal {L}}}_1^{h/2}, \right. \\&\quad \left. (\gamma _{q'},\beta _{q''})\in {{\mathcal {L}}}_2^{h/2}\right\} \, , \end{aligned}$$

obtaining an immersed Lagrangian submanifold. Thus, it generates a new symplectic integrator. Moreover, it is possible to compose more than two Lagrangian submanifolds to generate more involved methods where the intermediate points $\gamma _{q'}$ play the role of micro-nodes (see [11, 35, 41]).

Example 6.1

Let $Q={{\mathbb {R}}}^n$, we consider the two discretization maps $R_{d,1}(q, v)=(q, q+v)$ and $R_{d,2}(q, v)=(q-v, q)$. Then, we compute their corresponding cotangent lifts, $R^{T^*}_{d,1}$ and $R^{T^*}_{d,2}$ as described in Sect. 3.2, and obtain

$$\begin{aligned} {{\mathcal {L}}}^{h/2}_1= & {} \left\{ (q_k, p_k; q_{k+1/2}, p_{k+1/2})\;\left| \; \begin{array}{c} p_{k+1/2}=p_k-\frac{h}{2} \nabla _q H(q_k, p_{k+1/2})\\ q_{k+1/2}=q_k+\frac{h}{2}\nabla _p H(q_k, p_{k+1/2}) \end{array} \right. \right\} \\ {{\mathcal {L}}}^{h/2}_2= & {} \left\{ (q_{k+1/2}, p_{k+1/2}; q_{k+1}, p_{k+1})\; \left| \; \begin{array}{c} p_{k+1}=p_{k+1/2}-\frac{h}{2} \nabla _q H(q_{k+1}, p_{k+1/2})\\ q_{k+1}=q_{k+1/2}+\frac{h}{2} \nabla _p H(q_{k+1}, p_{k+1/2}) \end{array} \right. \right\} \end{aligned}$$

The composition ${{\mathcal {L}}}^{h/2}_2\circ {{\mathcal {L}}}^{h/2}_1$ gives a new symplectic integrator that corresponds with the Störmer–Verlet method [27]:

$$\begin{aligned}&p_{k+1/2} = p_k-\frac{h}{2} \nabla _q H(q_k, p_{k+1/2})\, ,\\&q_{k+1}-\frac{h}{2} \nabla _p H(q_{k+1}, p_{k+1/2}) = q_k+\frac{h}{2}\nabla _p H(q_k, p_{k+1/2})\, ,\\&p_{k+1}= p_{k+1/2}-\frac{h}{2} \nabla _q H(q_{k+1}, p_{k+1/2}). \end{aligned}$$

$\square$

When discrete Lagrangian functions are given as in Sect. 5.3, Eqs. (20) and (23) can be expressed as Lagrangian submanifolds of $(T^*Q\times T^*Q, \Omega _{12})$ and many of the methods described in [35, 41] are recovered.

For instance, for a small positive step size h, we consider the following two discretization maps on Q, $R_{d,i}:TQ \rightarrow Q\times Q$:

$$\begin{aligned} R_{d,1}\left( q, \frac{h}{2}v\right) =\left( q, q+\frac{h}{2}v\right) ,&\text{ with } \text{ inverse }&R_{d,1}^{-1}(q_0,q_1)=\left( q_0, \dfrac{q_1-q_0}{h/2}\right) \, , \\ R_{d,2}\left( q, \frac{h}{2}v\right) =\left( q-\frac{h}{2}v, q\right) ,&\text{ with } \text{ inverse }&R_{d,2}^{-1}(q_0,q_1)=\left( q_1, \dfrac{q_1-q_0}{h/2}\right) . \end{aligned}$$

We define the discrete Lagrangian functions $L_{d,i}=\left( h\, L\circ \left( R_{d,i}\right) ^{-1}\right) :Q\times Q \rightarrow {\mathbb {R}}$ such that the image of $(\Phi ^{-1}\circ \, {\mathrm{d}} L_{d,i})$ define the Lagrangian submanifolds ${{\mathcal {L}}}_i$ of the symplectic manifold $(T^*Q \times T^*Q, \Omega _{12})$. The submanifolds ${{\mathcal {L}}}_i$ define a discrete dynamical system whose equations are locally described by

$$\begin{aligned} {{\mathcal {L}}}_i\!=\!\{(q_0,p_0;q_1,p_1) \in T^*Q \!\times \! T^*Q \; \mid \; p_0\!=\!-{\mathrm{D}}_1 {\mathrm{L}}^i_d(q_0,q_1), \quad p_1\!=\!{\mathrm{D}}_2 {\mathrm{L}}^i_d(q_0,q_1)\}. \end{aligned}$$

The composition

$$\begin{aligned} {{\mathcal {L}}}_{2}\circ {{\mathcal {L}}}_1=\{(\alpha _1,\alpha _2)\; \mid \; \exists \; \alpha _{1/2}\in T^*Q \text{ s. } \text{ t. } (\alpha _1,\alpha _{1/2})\in {{\mathcal {L}}}_1 \, , \, (\alpha _{1/2},\alpha _2)\in {{\mathcal {L}}}_2\}. \end{aligned}$$

has associated the dynamics given by the discrete Lagrangian ${\mathrm{L}}^3_d(q_0,q_2)={\mathrm{L}}^1_d(q_0,q_1)+{\mathrm{L}}^2_d(q_1,q_2)$ that plays the role of generating function (see also [21]). The discrete equations are

$$\begin{aligned} p_0= & {} -{\mathrm{D}}_1 {\mathrm{L}}_d^1(q_0,q_1)\, ,\\ 0= & {} {\mathrm{D}}_2 {\mathrm{L}}^1_d(q_0,q_1)+{\mathrm{D}}_1 {\mathrm{L}}_d^2(q_1,q_2)\, ,\\ p_2= & {} {\mathrm{D}}_2 {\mathrm{L}}^2_d(q_1,q_2). \end{aligned}$$

6.1 Symplectic Symmetric Methods

If we have a Lagrangian submanifold ${{\mathcal {L}}}$ of $(T^*Q\times T^*Q, \Omega _{12})$, then the transpose ${{\mathcal {L}}}^{\dagger }$ defined by

$$\begin{aligned} {{\mathcal {L}}}^{\dagger }=\{(\alpha _q, \beta _{q'})\in T^*Q\times T^*Q\; \mid \; (\beta _{q'}, \alpha _q)\in {{\mathcal {L}}}\} \end{aligned}$$

is also a Lagrangian submanifold of $(T^*Q\times T^*Q, \Omega _{12})$.

For a Hamiltonian function and a discretization map on $T^*Q$, we consider the following Lagrangian submanifold used in the previous section:

$$\begin{aligned} {{\mathcal {L}}}^h=\left\{ (\alpha _q, \beta _{q'})\in T^*Q\times T^*Q\; \mid \; \exists \; \gamma _{q''}\in T^*Q \text{ s. } \text{ t. } (\alpha _q; \beta _{q'})=R_d^{T^*}( h\, X_H \left( \gamma _{q''}\right) ) \right\} . \end{aligned}$$

As described in [27, 41], the composition of symplectic methods (seen here as Lagrangian submanifolds) gives rise to new symplectic methods. For instance, the Lagrangian submanifold

$$\begin{aligned} \left( {{\mathcal {L}}}^{h/2}\circ {{\mathcal {L}}}^{h/2}\right) ^{\dagger }\, \end{aligned}$$

is another way to interpret the Störmer–Verlet method considered in the previous section.

Definition 6.1

A symplectic method defined by ${{\mathcal {L}}}^h$ is symmetric if

$$\begin{aligned} \left( {{\mathcal {L}}}^h\right) ^{\dagger }={{\mathcal {L}}}^{-h}. \end{aligned}$$

Proposition 6.1

Let $\iota : Q\times Q\rightarrow Q\times Q$ be the inversion map defined by $\iota (q, q')=(q', q)$. If $R_d(v_q)=\iota (R_d(-v_q))$ for all $v_q\in T_qQ$, then ${{\mathcal {L}}}^h$ is symmetric.

Proof

Observe that

$$\begin{aligned} ({{\mathcal {L}}}^h)^{\dagger }=\left\{ (\alpha _q,\beta _{q'})\in T^*Q\times T^*Q\; \mid \; \exists \; \gamma _{q''}\in T^*Q \text{ s. } \text{ t. } (\beta _{q'},\alpha _q)=R_d^{T^*}( h\, X_H \left( \gamma _{q''}\right) )\right\} . \end{aligned}$$

By Proposition 3.7, $R_d^{T^*}$ is symmetric and if we apply the inversion ${\iota }_{T^*Q}: T^*Q\times T^*Q\rightarrow T^*Q\times T^*Q$ on $T^*Q$ on both sides of the equation defining $({{\mathcal {L}}}^h)^{\dagger }$ we get:

$$\begin{aligned} \iota _{T*Q}(\beta _{q'},\alpha _q)=\iota _{T^*Q}\left( R_d^{T^*}( h\, X_H \left( \gamma _{q''}\right) )\right) =R_d^{T^*}( -h\, X_H \left( \gamma _{q''}\right) ). \end{aligned}$$

Thus, we immediately deduce that $({{\mathcal {L}}}^h)^{\dagger }={{\mathcal {L}}}^{-h}$. $\square$

It is well-known that the order of a symmetric method is necessarily even, then using symmetric discretization maps as in Definition 2.3 we always obtain a second-order method.

6.2 Construction of Higher-Order Symplectic Methods

In the previous sections we have introduced first- and second-order symplectic methods starting with different discretization maps. Now, we will show that the composition of Lagrangian submanifolds is a geometric tool to produce higher-order symplectic methods equivalent to the composition of numerical methods (see [7, 27, 33, 56])

From an initial discretization map $R_d: TQ\rightarrow Q\times Q$ and a Hamiltonian system $X_H$, we construct the Lagrangian submanifold ${{\mathcal {L}}}^h$ as in Eq. (24). For real numbers $\gamma _1, \ldots , \gamma _s$, we define the Lagrangian submanifold

$$\begin{aligned} {{\mathcal {L}}}^{\gamma _sh}\circ \ldots \circ {{\mathcal {L}}}^{\gamma _1h} \end{aligned}$$

(25)

that generates a symplectic composition method.

If ${{\mathcal {L}}}^h$ generates a method of order two and the coefficients $\gamma _1, \ldots , \gamma _s$ verify

$$\begin{aligned}&\gamma _1+ \ldots + \gamma _s = 1\, ,\\&\gamma ^3_1+ \ldots + \gamma ^3_s = 0, \end{aligned}$$

then the symplectic composition method (25) is at least of order 3.

As described in [27] for $s=3$, if we start with a method of order two

$$\begin{aligned} \gamma _1 =\gamma _3= \frac{1}{2-2^{1/3}}, \quad \gamma _2= -\frac{2^{1/3}}{2-2^{1/3}}\, , \end{aligned}$$

we obtain a method of order 4 due to the symmetry of the coefficients. By repeating this procedure, we obtain methods of order 6, 8, etc. Of course, other choices of the parameters produce different higher-order numerical methods (see, for instance, [45, 49]). Additionally, other techniques like splitting methods perfectly fit in our framework as composition of Lagrangian submanifolds.

Another interesting family of methods are the symplectic Runge–Kutta methods. For instance, the diagonally implicit Runge–Kutta methods (that is, the coefficients verify $a_{ij}=0$ if $i<j$) and $b_i\not =0$ are derived from the symmetric discretization map $R_d(q, v)=(q-v/2, q+v/2)$. For a Hamiltonian function $H: T^*Q \rightarrow {{\mathbb {R}}}$, the map $R_d$ produces the Lagrangian submanifold of $T^*Q\times T^*Q$:

$$\begin{aligned} {{\mathcal {L}}}^h=\left\{ (q_k, p_k, q_{k+1}, p_{k+1})\in T^*Q \times T^*Q \; \left| \; \ \begin{array}{r} \frac{q_1-q_0}{h}=\frac{\partial H}{\partial p}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) \, \\ \\ \frac{p_1-p_0}{h}= -\frac{\partial H}{\partial q}\left( \frac{q_0+q_1}{2}, \frac{p_0+p_1}{2}\right) \, \end{array}\right. \right\} . \end{aligned}$$

This Lagrangian submanifold corresponds with the implicit midpoint rule. It is well-known that the symplectic diagonally implicit Runge–Kutta methods are equivalent to the following composition [27]

$$\begin{aligned} {{\mathcal {L}}}^{b_1h}\circ {{\mathcal {L}}}^{b_2h}\circ \ldots {{\mathcal {L}}}^{b_sh}. \end{aligned}$$

Similar argument also holds for partitioned Runge–Kutta method based on two diagonally implicit methods

To sum up, our constructions allow to reinterpret other well-known techniques for designing higher-order methods. Specifically, from second-order methods obtained by discretization maps in Sect. 6.1 we can derive higher-order methods using composition of Lagrangian submanifolds. These techniques are not limited to standard Hamiltonian systems but can also be used for Hamiltonian systems defined on cotangent bundles of manifolds (as Lie groups [8], for instance) or more general situations as we will describe in Sect. 7.

7 Conclusions and Future Work

In this paper we have introduced the lift of a discretization map to tangent and cotangent bundles which are the phase spaces of mechanical systems. These lifts allow us to derive geometric integrators for systems defined by a Lagrangian or Hamiltonian function. Standard constructions in symplectic geometry, as well as properties of Lagrangian submanifolds, create a geometric framework to obtain several well-known symplectic integrators. Our geometric point of view opens the door for new types of applications of discretization maps, as well as for the construction of geometric (symplectic) integrators following simple rules (lifting of retractions, composition and generating functions for Lagrangian submanifolds, etc.). Now, we will mention some promising future research lines.

7.1 Reduced Systems and Systems with Holonomic Constraints

The notion of lift of discretization maps can be easily extended to the Lie algebroid setting using the Lie algebroid prolongation [20]. This theory covers all the examples of reduced systems by symmetry groups. Therefore, combining both constructions we can directly apply our results to the construction of geometric integrators for Lagrangian or Hamiltonian functions invariant under the action of a symmetry Lie group [38, 54]. It is also well-known how to produce geometric integrators for systems subject to holonomic constraints. Thus, it will be interesting to produce constrained geometric integrators using our approach and compare them with [33, 44].

7.2 Discrete Gradient Methods

In general for ordinary differential equations in ${{\mathbb {R}}}^n$ in skew-gradient form, i.e., ${\dot{x}}=\Pi (x)\nabla H(x)$ where $x\in {\mathbb {R}}^n$ and $\Pi (x)$ is a skew-symmetric matrix, it is clear that H is a first integral. Using discretizations of the gradient $\nabla H(x)$, it is possible to define a class of integrators that preserve exactly the first integral H (see [24, 43]). They are defined as follows: let $H:{\mathbb {R}}^n\rightarrow {\mathbb {R}}$ be a differentiable function, then ${\overline{\nabla }}H:{\mathbb {R}}^{2N}\longrightarrow {\mathbb {R}}^N$ is a discrete gradient of H if it is continuous and satisfies

$$\begin{aligned} {\overline{\nabla }}H(x,x')^T(x'-x)&=H(x')-H(x)\, , \quad \, \text{ for } \text{ all } x,x' \in {\mathbb {R}}^n \, , \nonumber \\ {\overline{\nabla }}H(x,x)&=\nabla H(x)\, , \quad \quad \quad \quad \text{ for } \text{ all } x \in {\mathbb {R}}^n . \end{aligned}$$

(26)

For a Hamiltonian system $H: T^*Q\rightarrow {{\mathbb {R}}}$, we can generalize the previous construction by using a discretization map $R_d^{TT^*Q}: TT^*Q \rightarrow T^*Q\times T^*Q$ on a general differentiable manifold $T^*Q$. We define a discrete gradient as a map $\overline{{d}{H}}: T^*Q\times T^*Q \longrightarrow T^{*}T^*Q$ that makes the following diagram commutative

Similar to (26), ${\overline{d}}{H}$ must verify the following properties:

$$\begin{aligned} \begin{array}{cl} \langle \overline{{d}{H}}(x,x'), \left( R_d^{TT^*Q}\right) ^{-1}(x, x')\rangle =H(x')-H(x)\, , &{}\quad \text{ for } \text{ all } x,x' \in T^*Q \, , \\ \overline{{d}}H(x,x) =d H(x)\, , &{}\quad \text{ for } \text{ all } x \in T^*Q . \end{array} \end{aligned}$$

In this case, an energy preserving integrator would be

$$\begin{aligned} \left( R_d^{TT^*Q}\right) ^{-1}(x, x')=\omega ^{\sharp }(x'')(\overline{{d}{H}}(x, x')) \end{aligned}$$

where $x''=\tau _{T^*Q}^*\left( \left( R_d^{TT^*Q}\right) ^{-1}(x, x')\right)$. We will explore this possibility in a future paper (see also [12, 13], for the case of extension to manifolds, in special, Riemannian manifolds).

7.3 Higher-Order Lagrangian Systems

Another topic of interest is related with the discretization of higher-order Lagrangian systems $L: T^{(k)}Q\rightarrow {{\mathbb {R}}}$ using appropriate higher-order lifts of discretization maps. These constructions will be useful for interpolation problems on manifolds and for optimal control problems (see [14, 17, 23]). For instance, some interesting optimal control problems for mechanical systems are describe using a second-order Lagrangian system $L: T^{(2)}Q\rightarrow {{\mathbb {R}}}$ (position, velocities and accelerations). The corresponding Hamiltonian formulation is given in the cotangent bundle $T^*T Q$. To obtain a symplectic integrator, we could use first the tangent lift of a discretization map $R_d: TQ\rightarrow Q\times Q$ and finally its cotangent lift.

7.4 Higher-Order Geometric Methods

Section 6.2 is the starting point of a research line that seeks to obtain higher-order geometric methods by using discretization maps. In this paper, we have only briefly discussed composition and splitting methods, but in the future, it would be interesting to explore other possibilities, as, for instance, the well-known family of implicit symplectic Runge–Kutta methods and symplectic partitioned Runge–Kutta method. Moreover, the use of discretization maps on the cotangent bundle that are not lifted from the base manifold may be useful for obtaining higher-order methods.

7.5 Geometric Integration of Dirac Systems

Dirac structures were introduced in [15, 16] as a way to unify presymplectic and Poisson geometries giving a way to collect in the same geometric framework many situations of interest in mechanics and mathematical physics. As an example we can think of the Dirac structure $D\subset TT^*Q\oplus T^*T^*Q$ induced by the canonical symplectic structure $\omega _Q$, that is

$$\begin{aligned} D=\{(X_{\alpha _q}, \lambda _{\alpha _q})\in TT^*Q\oplus T^*T^*Q\; \mid \; i_{X_{\alpha _q}}\omega _Q=\lambda _{\alpha _q}\}. \end{aligned}$$

For a Hamiltonian system $H:T^*Q\rightarrow {{\mathbb {R}}}$, we can write Hamilton’s equations as

$$\begin{aligned} X_{\alpha _q}\oplus {\mathrm{d}} H(\alpha _q)\in D_{\alpha _q}\, , \end{aligned}$$

(27)

with $\alpha _q\in T_q^*Q$. As studied in the literature, Dirac structures are more general than the above example. They can be given by a Poisson tensor, for instance, or could not satisfy the integrability condition admitting new generalizations as in the case of nonholonomic constraints (see [6]). Moreover, it is also interesting to study the case of Dirac systems where the dynamics is not induced by a function on the cotangent bundle (as in the case of standard Hamiltonian dynamics), but for a general Lagrangian submanifold ${{\mathcal {S}}}$ of $(T^*T^*Q, \omega _{T^*Q})$ (as in the case of singular Lagrangians, optimal control theory, etc.). Now, Eq. (27) must be replaced by

$$\begin{aligned} X_{\alpha _q}\oplus S_{\alpha _q}\in D_{\alpha _q}. \end{aligned}$$

discretization maps could also be used here to deduce geometric integrators. Let $R_d: TQ\rightarrow Q\times Q$ be a discretization map on Q, the cotangent lift of $R_d$ defines the following geometric integrator

$$\begin{aligned} \left( R_d^{T^*}\right) ^{-1}(q_k, p_k; q_{k+1}, p_{k+1}) \oplus S_{\gamma _{k, k+1}}\in D_{\gamma _{k, k+1}} \end{aligned}$$

where $\tau _{T^*Q}\left( (R_d^{T^*})^{-1}(q_k, p_k; q_{k+1}, p_{k+1})\right) =\gamma _{k, k+1}$. We will study in a forthcoming paper the design of Dirac integrators using the lift of the discretization map, their geometrical properties (preservation of the associated presymplectic foliation, etc.) and compare them with other approaches on this topic ([34, 35]).

7.6 Morse Families for Lagrangian Submanifolds and Symplectic Integration

In this paper it is clear the close relationship between the design of different symplectic methods and the construction of Lagrangian submanifolds. The notion of a Morse family or phase function was introduced in [28] (see also [53]) and it is possible to prove that locally any Lagrangian submanifold is the image of a Lagrangian immersion generated by a Morse family. In a recent paper [6], we have combined Dirac structures and Morse families to obtain a geometric formalism that unifies most of the scenarios in mechanics (constrained calculus, nonholonomic systems, optimal control theory, higher-order mechanics, etc.), as the examples in the paper show. Employing the techniques introduced here, we aim to study the construction of geometric integrators for all the above-mentioned cases, as well as the notion of Morse family to construct new geometric integrators.