1 Introduction↩︎

The discretization of Hamiltonian systems has long been a central topic in geometric numerical integration [1]–[3]. Symplectic integrators, in particular, are well known for preserving the underlying geometric structure of Hamiltonian dynamics, providing consistent tools for long-term stability and energy analysis. Their use across different areas of physics, control, and applied mathematics has motivated a broad class of structure-preserving schemes that retain essential invariants of continuous dynamics, such as volume and momentum conservation.

Although the geometric theory of discrete Hamiltonian systems is well established for autonomous dynamics (see, e.g., [2], [4]), the case of non-autonomous systems remains less systematically treated. A possible way to study non-autonomous mechanical systems consists in considering the extended phase space \(T^*(Q \times \mathbb{R})\), with local coordinates \((q,t, p, p_t)\), where the time is treated as variable with its own conjugate momentum. In this setting, the dynamics is governed by the canonical extended symplectic form on \(T^*(Q\times \mathbb{R})\) \[\Omega_0 = dq\wedge dp + dt\wedge dp_t,\] while the time-dependent space \(T^*Q \times \mathbb{R}\) is endowed with a different geometry induced by the canonical cosymplectic structure \((\omega_0,\eta) = (dq\wedge dp,\,dt)\). Observe that now \(\omega_0\) is a presymplectic 2-form since \(\partial/\partial t\) is in \(\ker \omega_0\). A central question then arises: can one construct numerical integrators that preserve the canonical geometric structures associated to this cosymplectic framework?

Marthinsen and Owren in [5] address this problem by characterizing discrete canonical transformations as symplectomorphisms on the extended space \(T^*(Q \times \mathbb{R})\). Their analysis provides local matrix conditions that ensure that the extended symplectic form \(\Omega_0\) is preserved by the discrete flow. However, their construction remains entirely within the symplectic setting and does not explicitly describe how such transformations act on the projected, time-dependent space \(T^*Q \times \mathbb{R}\), where the relevant invariants are the cosymplectic volume and the Poisson bracket associated with \((\omega_0,\eta)\).

To better appreciate the role of these invariants, consider a time-dependent Hamiltonian \(H:T^*Q \times \mathbb{R}\to\mathbb{R}\). For sake of simplicity, take \(Q=\mathbb{R}\) and \(H(q,p)=\tfrac12e^{-t}(p^2+q^2)\), so that the continuous dynamics describes a harmonic oscillator with dissipation. The exact flow is given by \[\begin{pmatrix}q_{k+1}\\p_{k+1}\end{pmatrix} = \begin{pmatrix} \cos\theta_h & \sin\theta_h \\ -\sin\theta_h & \cos\theta_h \end{pmatrix} \begin{pmatrix} q_k \\ p_k \end{pmatrix},\] where \(\theta_h = e^{-t_k} - e^{-(t_k+h)},\) and it preserves the canonical Poisson bracket or the symplecticity on fibers on the projection of \(T^*Q \times \mathbb{R}\) to \({\mathbb{R}}\) since \[\begin{pmatrix} \cos\theta_h & \sin\theta_h \\ -\sin\theta_h & \cos\theta_h \end{pmatrix}^T\begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \cos\theta_h & \sin\theta_h \\ -\sin\theta_h & \cos\theta_h \end{pmatrix}= \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}\] However, if we use an explicit Euler discretization on \(T^*Q \times \mathbb{R}\): \[q_{k+1}=q_k+h\,\partial_p H(q_k, t_k, p_k),\quad t_{k+1}=t_k+h,\quad p_{k+1}=p_k-h\,\partial_q H(q_k, t_k, p_k),\] and we denote by \(\Phi_k:(q_k, t_k, p_k)\mapsto (q_{k+1}, t_{k+1}, p_{k+1})\) the one–step map defined by this explicit Euler update, then \[\begin{pmatrix}q_{k+1}\\p_{k+1}\end{pmatrix} = \begin{pmatrix} 1 & he^{-t_k} \\ - he^{-t_k} & 1 \end{pmatrix} \begin{pmatrix} q_k \\ p_k \end{pmatrix},\] but \[\det \begin{pmatrix} 1 & he^{-t_k} \\ - he^{-t_k} & 1 \end{pmatrix}=1+h^2e^{-2t_k}>1\] and then the integrator is not symplectic. This simple example highlights the importance of enforcing symplecticity on each time slice \(T^*Q \times \{t\}\) for the numerical methods: without “slice-simplecticity”, even elementary oscillators exhibit spurious energy drift, whereas structure-preserving schemes maintain bounded error and reproduce the correct geometric features of the continuous flow.

The approach developed in this paper directly addresses this limitation. We construct discrete canonical transformations on the extended phase space \(T^*(Q \times \mathbb{R})\) whose projections onto \(T^*Q \times \mathbb{R}\) define structure-preserving integrators for non-autonomous Hamiltonian systems[6]. This construction relies on the notion of discretization maps [7], which provide a geometric mechanism to approximate flows on tangent and cotangent bundles while preserving their intrinsic structures. By applying this formalism to the extended setting, we obtain a discrete flow that preserves both the volume \(\omega_0^{\,n}\wedge\eta\) and the Poisson structure associated with the cosymplectic structure \((\omega_0,\eta)\), thereby providing a faithful discrete analogue of the continuous non-autonomous Hamiltonian dynamics. Crucially, for each fixed \(t\), the induced map \(\phi_{t,h}:T^*Q\to T^*Q\) is a symplectomorphism, ensuring per–time-slice symplecticity and geometric consistency of the discrete evolution (see [8] for a discrete Lagrangian description).

Our formulation extends the local results of [5] to a global, coordinate-free framework that connects symplectic and cosymplectic discretizations through the projection \(\mu:T^*(Q \times \mathbb{R})\to T^*Q \times \mathbb{R}\). This provides the first geometric proof of structure preservation for discrete canonical transformations in the non-autonomous case, bridging the discretization-map formalism of [7] with the canonical transformation theory of [6].

Paper outline. Section 2 reviews the geometric background on extended phase spaces, cosymplectic structures, and Poisson morphisms. Section 3 introduces discrete canonical transformations generated by discretization maps and establishes their preservation properties. Section 4 presents numerical examples, including dissipative and mechanical systems with time-dependent potential, that illustrate the behaviour of the proposed integrators compared to standard Runge–Kutta schemes. Section 5 shows an example of how the system is well-suited for applications in gauge-invariant systems. Finally, Section 6 summarizes the main conclusions and outlines future research directions.

All the functions, vector fields, forms and objects treated in this paper are assumed to be smooth unless explicitly stated.

2 Geometry of non-autonomous Hamiltonian systems↩︎

2.1 Hamiltonian mechanics↩︎

Let \(Q\) be an \(n\)-dimensional configuration manifold of a mechanical system with local coordinates \((q^A)\), \(1\leq A\leq n\) [9]. If \(T_{q}Q\) denotes the tangent space of \(Q\) at the point \(q\), the tangent bundle \(TQ\) is defined by \(\displaystyle{TQ:=\sqcup_{q\in Q}T_{q}Q}\), with local coordinates \((q^A, \dot{q}^A)\). The bundle projection \(\tau_{Q}: TQ \rightarrow Q\), sends each vector \(v_{q}\) to the corresponding base point \(q\), in local coordinates \(\tau_{Q}(q^{A},\dot{q}^{A})=q^{A}\).

\(\displaystyle{T^{*}Q:=\sqcup_{q\in Q}T^{*}_{q}Q}\), endowed with local coordinates \((q^A,p_A)\), indicates the cotangent bundle. \(\pi_{Q}:T^{*}Q \rightarrow Q\) denotes the canonical projection and sends each momentum \(p_{q}\) to the corresponding base point \(q\). In local coordinates \(\pi_{Q}(q^{A},p_{A})=q^{A}\).

We associate a system of \(2n\) first order ordinary differential equations to every function \(H: T^{*}Q\to\mathbb{R}\), namely the associated Hamilton’s equations \[\label{hameq1} \dot{q}^{A}=\frac{\partial H}{\partial p_A},\quad\dot{p}_{A}=-\frac{\partial H}{\partial q^A},\quad 1\leq A\leq n.\tag{1}\] Their solutions are integral curves of the Hamiltonian vector field \(X_H\) associated with \(H\) that locally reads \(X_H(q,p)=\left(q, p, \frac{\partial H}{\partial p},-\frac{\partial H}{\partial q}\right)\).

A symplectic 2-form \(\omega\) on a manifold \(Q\) is a \((0,2)\)-type tensor field that is skew-symmetric and non-degenerate, i.e., \(\omega(X,Y)=-\omega(Y,X)\) for all vector fields \(X\) and \(Y\) and if \(\omega(X,Y)=0\) for all vector fields \(X\), then \(Y\equiv 0\). A pair \((Q,\omega)\) is called a symplectic manifold if \(Q\) is a manifold and \(\omega\) is a symplectic 2-form.

The set of vector fields and the set of 1-forms on \(Q\) are denoted by \(\mathfrak{X}(Q)\) and \(\Omega^{1}(Q)\), respectively.

As described in [10], the cotangent bundle \(T^*Q\) of a manifold \(Q\) is equipped with a canonical exact symplectic structure \(\omega_Q=-d\theta_Q\), where \(\theta_Q\) is the canonical 1-form on \(T^*Q\). In canonical bundle coordinates \((q^A, p_A)\) on \(T^*Q\), \(\theta_Q= p_A\, \mathrm{d}q^A\) and \(\omega_Q= \mathrm{d}q^A\wedge \mathrm{d}p_A.\)

The Hamiltonian dynamics is characterized by the following two essential properties [10]:

• the Hamiltonian function, hence the state of the system, is preserved along the dynamics: \(\omega_Q(X_H, X_H)=dH(X_H)=X_H(H)=0\);

• the symplectic form \(\omega\) is preserved along the dynamics, that is, if \(\Phi^t_{X_H}\) denotes the flow of \(X_H\), then the pull-back of the symplectic form \(\omega_Q\) by the flow is preserved, i.e., \((\Phi^t_{X_H})^*\omega_Q=\omega_Q\).

Let \(Q_1\) and \(Q_2\) be \(n\)-dimensional manifolds and \(F: Q_1\rightarrow Q_2\) be a map. The tangent lift \(TF: TQ_1\rightarrow TQ_2\) of \(F\) is defined by \(TF(v_q)=T_qF (v_q)\in T_{F(q)} Q_2\, \, where v_q\in T_qQ_1\), and \(T_qF\) is the tangent map of \(F\) whose matrix representation is the Jacobian matrix of \(F\) at \(q\in Q_1\).

As the tangent map \(T_qF\) is linear, the dual map \(T_{q}^*F\colon T^*_{F(q)}Q_2\rightarrow T^*_qQ_1\) is defined as follows: \[\langle(T^*_{q}F)(\alpha_2), v_{q}\rangle=\langle \alpha_2, T_{q}F(v_{q})\rangle for everyv_q\in T_qQ_1.\] Note that \((T^*_{q}F)(\alpha_2)\in T^*_qQ_1\).

2.2 Time-dependent Hamiltonian systems and cosymplectic geometry↩︎

To describe the dynamics of a non-autonomous Hamiltonian system, we consider the manifold \(T^*Q \times \mathbb{R}\) equipped with \((\omega_0, \eta)\), where \(\omega_0=pr_1^*\omega_Q\) and \(\eta=pr_2^* dt\). Here \(pr_1: T^*Q \times \mathbb{R} \rightarrow T^*Q\) and \(pr_2: T^*Q \times \mathbb{R} \rightarrow \mathbb{R}\) denote the projections.

Note that the 2-form \(\omega_0\) is presymplectic, since its kernel is generated by \(\partial_t\). This degeneracy reflects the fact that the time variable \(t\) plays the role of an external parameter rather than a canonical coordinate. Consequently, the pair \((\omega_0, \eta)\) provides the canonical cosymplectic structure on \(T^*Q \times \mathbb{R}\), which constitutes the natural geometric framework for time-dependent Hamiltonian systems (see [11] and references therein for details).

Since \(\omega_0^n\wedge dt\not= 0\), then \((T^*Q \times \mathbb{R}, \omega_0, \eta)\) is a cosymplectic manifold. In Darboux coordinates \((q^i, t, p_i)\), we have that \[\omega_0=dq^i\wedge dp_i, \qquad \eta=dt.\]

The cosymplectic structure defines a unique vector field \(R\), the Reeb vector field, that satisfies \[i_R\omega_0=0, \qquad i_R \eta=1.\] In the particular case studied in this work, the Reeb vector field \(R=\partial/\partial t=\partial_t\).

For any function \(f\in C^{\infty}(T^*Q \times \mathbb{R})\) we define the Hamiltonian vector field \(X_f\) by \[i_{X_f}\omega_0=df-R(f)\eta, \qquad i_{X_f}\eta=0,\] and the evolution vector field \(E_f\) by \(E_f=R+X_f\) or equivalently \[\label{def:evolution} i_{E_f}\omega_0=df-R(f)\eta, \qquad i_{E_f}\eta=1\;.\tag{2}\] The Lie derivative of \(\omega_0\) and \(\eta\) along \(E_f\) produces \[\mathcal{L}_{E_f}\omega_0=-d(R(f))\wedge \eta, \qquad \mathcal{L}_{E_f}\eta=0.\]

Locally, the evolution vector field of the Hamiltonian function \(H\) reads \[E_H=R+X_H=\frac{\partial}{\partial t}+\frac{\partial H}{\partial p_A}\frac{\partial }{\partial q^A}-\frac{\partial H}{\partial q^A}\frac{\partial }{\partial p_A}.\] Its integral curves coincide with the trajectories defined by Hamilton’s equations \[\dot{q}=\partial_p H, \quad \dot{t} = 1, \quad \dot{p}=-\partial_q H.\]

The cosymplectic structure \((\omega_0, \eta)\) induces a Poisson bracket \(\{\cdot,\cdot\}_0\) on \(T^*Q \times \mathbb{R}\) by \[\{f, g\}_0=\omega_0(X_f, X_g)=\omega_0(E_f, E_g)\] for all \(f, g\in C^{\infty}(T^*Q \times \mathbb{R})\), that in Darboux coordinates is given by the bivector \[\Lambda_0=\frac{\partial}{\partial q^A}\wedge \frac{\partial}{\partial p_A}\; .\] Observe that \(\mathcal{L}_{E_H}\Lambda_0=\mathcal{L}_{R}\Lambda_0+\mathcal{L}_{X_H}\Lambda_0=0\) since \(\mathcal{L}_R\Lambda_0=0\). The canonical volume form on \(T^*Q \times \mathbb{R}\) \(\mathrm{vol}_0=\omega_0^n\wedge \eta\) is preserved by the flow of \(E_H\) as well.

An interesting way to represent the dynamics of the evolution vector field \(E_H\) is to consider the modified cosymplectic structure by taking the cosymplectic pair \((\Omega_H, dt)\) where \[\Omega_H=\omega_0+dH\wedge \eta,\quad \eta=dt,\] with corresponding Reeb vector field precisely the evolution vector field: \[i_{E_H}\Omega_H=0,\quad i_{E_H}\eta=1.\] Therefore, the flow of \(E_H\) satisfies \[\mathcal{L}_{E_H}\Omega_H = 0, \quad \mathcal{L}_{E_H}dt = 0.\]

We stress the fact that this property is difficult to preserve using numerical methods, since \(\Omega_H\) depends on the Hamiltonian function \(H\). Therefore, we aim to obtain infinitesimal preservation properties in terms of canonical objects.

We also observe that the Evolution vector field \(E_H\) preserves:

1. the volume \(\mathrm{vol}_0=\omega_0^n\wedge \eta\), that is \(\mathcal{L}_{E_H}\mathrm{vol}_0=0\);

2. the Poisson structure \(\Lambda\), that is \(\mathcal{L}_{E_H}\Lambda_0=0\).

Denote now by \(\Psi^H_s\colon U\subset T^*Q \times \mathbb{R}\rightarrow T^*Q \times \mathbb{R}\) the flow of the evolution vector field \(E_H\), where \(U\) is an open subset of \(T^*Q \times \mathbb{R}\). Observe that if \(\alpha_q\in T^*_qQ\) \[\Psi^H_s(\alpha_q, t) = (\Psi^H_{t,s}(\alpha_q), t+s) \,,\] where \(\Psi^H_{t,s}(\alpha_q)={pr}_1(\Psi^H_s(\alpha_q, t))\). Therefore, we induce a map \[\Psi^H_{t,s}\colon U_t\subseteq T^*Q\to T^*Q\,,\] where \(U_t=\{ \alpha_q\in T^*Q \;|\; (t, \alpha_q)\in U \}\). In terms of the flow of \(E_H\) we have then that the flow of the evolution vector field \(E_H\) preserves:

1. volume: \((\Psi^H_s)^*\mathrm{vol}_0= \mathrm{vol}_0\);

2. Poisson bracket: \(\{f\circ \Psi^H_s, g\circ \Psi^H_s\}_0=\{f, g\}_0\circ \Psi^H_s\);

3. Symplectic form: \((\Psi^H_{t,s})^*\omega_Q=\omega_Q\).

These are the three related properties that are important to preserve when constructing appropriate geometric integrators.

2.3 Generalized canonical transformations↩︎

In order to ensure geometric preservations, we will use a symplectification method to work in the standard symplectic framework (see [6] and references therein for details). Given the Hamiltonian function \(H:T^*Q \times \mathbb{R}\rightarrow \mathbb{R}\) we can construct the extended Hamiltonian \(\tilde{H}:T^*(Q \times \mathbb{R})\rightarrow \mathbb{R}\) by taking \[\tilde{H}=H\circ \mu+p_t\] where \(\mu: T^* (Q \times \mathbb{R})\rightarrow T^*Q \times \mathbb{R}\) is the projection \(\mu( \alpha_q; t, p_t)=( \alpha_q, t)\), with \((t, p_t)\in T^*{\mathbb{R}}=\mathbb{R}\times {\mathbb{R}}^*\) and \(\alpha_q\in T^*_q Q\). In Darboux coordinates \((q^A,t,p_A,p_t)\) on \(T^*(Q \times \mathbb{R})\), the map \(\mu\) reads \[\mu(q^A,t,p_A,p_t) = (q^A,p_A,t) \in T^*Q \times \mathbb{R}.\]

\((T^*(Q \times \mathbb{R}), \omega_{Q \times \mathbb{R}})\) is a symplectic manifold with \(\omega_{Q \times \mathbb{R}}\) the canonical symplectic form, that in Darboux coordinates reads \[\omega_{Q \times \mathbb{R}} = dq^A\wedge dp_A + dt\wedge dp_t.\]

The proof is quite simple since the expression of \(X_{\tilde{H}}\) in coordinates is: \[\begin{align} X_{\tilde{H}}&=\frac{\partial \tilde{H}}{\partial p_A}\frac{\partial}{\partial q^A}-\frac{\partial \tilde{H}}{\partial q^A}\frac{\partial}{\partial p_i} +\frac{\partial \tilde{H}}{\partial p_t}\frac{\partial}{\partial t} -\frac{\partial \tilde{H}}{\partial t}\frac{\partial}{\partial p_t}\\ &=\frac{\partial}{\partial t}+\frac{\partial ({H}\circ \mu)}{\partial p_A}\frac{\partial}{\partial q^A}-\frac{\partial (H\circ \mu)}{\partial q^A}\frac{\partial}{\partial p_A} -\frac{\partial (H\circ \mu)}{\partial t}\frac{\partial}{\partial p_t}. \end{align}\] In \(T^*(Q \times \mathbb{R})\) the associated canonical bracket \(\{\, ,\,\}\) is given by \[\{\tilde{f}, \tilde{g}\}=\omega_{Q \times \mathbb{R}}(X_{\tilde{f}}, X_{\tilde{g}})\] for all \(\tilde{f}, \tilde{g}\in C^{\infty}(T^*(Q \times \mathbb{R}))\). In Darboux coordinates is represented by the bivector \[\Lambda=\frac{\partial}{\partial q^A}\wedge \frac{\partial}{\partial p_A}+ \frac{\partial}{\partial t}\wedge \frac{\partial}{\partial p_t} \; .\] Additionally, as an immediate consequence of the definition of the Poisson structures on \(T^*(Q \times \mathbb{R})\) and \(T^*Q \times \mathbb{R}\) it emerges that \(\mu: T^*(Q \times \mathbb{R})\rightarrow T^*Q \times \mathbb{R}\) is a Poisson map, that is, \[\{\mu^*f, \mu^*g\}=\mu^*\{f, g\}_0\] for all \(f, g\in C^{\infty}(T^*Q \times \mathbb{R})\).

The concept of canonical transformation does not apply directly to time-dependent mechanical systems in \(T^*Q \times \mathbb{R}\), but a similar notion can be formulated considering the extended phase space.

This definition requests that the diagram in Figure 1 commutes.

It is important to notice that this definition does not depend on the Hamiltonian function \(H\).

Now, consider \(\tau_s(q,p,t,p_t)=(q,p,t,p_t+s)\) the natural translation by \(s\in {\mathbb{R}}\) along the fibers. Then we have the following proposition.

Observe that the flow of the Hamiltonian \(\tilde{H}=H\circ \mu+p_t\) verifies that \((\tau_s)_* X_{\tilde{H}}=X_{\tilde{H}}\), therefore its flow commutes with the translation along the fibers. Now, our aim is to construct symplectic integrators that are canonical transformations and verify the commuting condition with respect to \(\tau_s\) obtaining preservation of Poisson bracket and volume form by applying Propositions 1 and 2.

3 Discrete canonical transformations and structure-preserving integrators↩︎

The first notion of retraction appearing in the literature can be found in [12], [13] from a topological viewpoint. Later on, the notion of retraction map as defined below is used to obtain Newton’s method on Riemannian manifolds [14], [15].

The condition ?? is known as local rigidity condition since, given \(\xi\in T_qQ\), the curve \(\gamma_\xi(t)=R_q(t\xi)\) has \(\xi\) as tangent vector at \(q\), i.e. \(\dot{\gamma}_\xi(t)= \langle DR_q(t\xi), \xi\rangle\; and, in consequence, \dot{\gamma}_\xi(0)= {\rm Id}_{T_qQ}(\xi)=\xi\).

A typical example of a retraction map is the exponential map, \(exp\), on Riemannian manifolds given for example in [16]. Therefore, the image of \(\xi\) through the exponential map is a point on the Riemannian manifold \((Q,g)\) obtained by moving along a geodesic a length equal to the norm of \(\xi\) starting with the velocity \(\xi/\|\xi\|\), that is, \(exp_q(\xi)=\sigma (\|\xi\|),\) where \(\sigma\) is the unit speed geodesic such that \(\sigma(0)=q\) and \(\dot{\sigma}(0)=\xi/\|\xi\|\).

Next, we define a generalization of the retraction map in Definition 4 that allows a discretization of the tangent bundle of the configuration manifold leading to the construction of numerical integrators as described in [7]. Given a point and a velocity, we obtain two nearby points that are not necessarily equal to the initial base point.

Thus, the discretization map \(R_d\) is a local diffeomorphism from some neighborhood of the zero section of \(TQ\).

If \(R^1(q,v)=q\), the two properties in Definition 5 guarantee that the both properties in Definition 4 are satisfied by \(R^2\). Thus, Definition 5 generalizes Definition 4.

3.1 Cotangent lift of discretization maps↩︎

As the Hamiltonian vector field takes value on \(TT^*Q\), the discretization map must be on \(T^*Q\), that is, \(R^{T^*}_d: TT^*Q \rightarrow T^*Q\times T^*Q\). Such a map is obtained by cotangently lifting a discretization map \(R_d\colon TQ\rightarrow Q\times Q\), so that the construction \(R^{T^*}_d\) is a symplectomorphism. In order to do that, we need the following three symplectomorphisms (see [7] for more details):

• The cotangent lift of the diffeomorphism \(R_d\colon TQ\rightarrow Q\times Q\) as described in Definition 2.

• The canonical symplectomorphism \(\alpha_Q\colon TT^*Q \longrightarrow T^*TQ\) such that
  \(\alpha_Q(q,p,v_q,v_p)= (q,v_q,v_p,p)\).

• The 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)\): \(\Phi:T^*Q\times T^*Q \longrightarrow T^*(Q\times Q)\), given by \(\Phi(q_0, p_0; q_1, p_1)=(q_0, q_1, -p_0, p_1).\)

Diagram in Figure 2 summarizes the construction process from \(R_d\) to \(R_d^{T^*}\).

3.2 Geometric integrators for non-autonomous systems↩︎

The aim of this section is to extend the construction of geometric integrators to the non-autonomous Hamiltonian setting. Our goal is to construct a discrete analogue of the continuous Hamiltonian flow that preserves the same geometric structures - namely, a canonical transformation on the extended phase space and its projection onto the cosymplectic manifold \(T^*Q \times \mathbb{R}\).

Given a smooth manifold \(M\), consider a discretization map \(R_d : TM \longrightarrow M\times M\), which associates to each tangent vector a discrete displacement between two nearby configurations. The differential of \(R_d\) allows us to define its cotangent lift \(\widehat{R}_d:T^*TM\longrightarrow T^*(M\times M)\), which provides the natural bridge between discrete Lagrangian and Hamiltonian formulations.

A key ingredient in the construction is the canonical anti-symplectomorphism \[{\mathcal{I}}_{TM}:T^*T^*M \longrightarrow T^*TM,\] which establishes a natural correspondence between the symplectic manifolds \((T^*T^*M,\omega_{T^*M})\) and \((T^*TM,\omega_{TM})\). Locally, this map is given by \[{\mathcal{I}}_{TM}(q,p,\mu_q,\mu_p) = (q,\mu_p,-\mu_q,p),\] and satisfies \({\mathcal{I}}_{TM}^*\omega_{TM} = -\omega_{T^*M}\), hence the name “anti-symplectomorphism”. This identification allows us to transfer the symplectic structure between tangent and cotangent bundles, which is essential for defining discrete flows that are consistent with the underlying continuous geometry.

The overall structure of the discrete construction is summarized in the commutative diagram of Figure 3. Each horizontal layer corresponds to a different geometric level: the bottom layer represents the discrete flow on the configuration manifold \(M\), the middle layer its tangent lift, and the top layer the corresponding lift to the cotangent bundle, where the discrete Hamiltonian evolution is defined as a Lagrangian submanifold of \(T^*T^*M\).

For non-autonomous systems, we take \(M = Q \times \mathbb{R}\) and construct the discretization map \(R^M_d\) as the product of two simpler maps: \[R^M_d = (R_d^{\mathbb{R}}, R_d^Q), \qquad R_d^{\mathbb{R}} : T^*\mathbb{R} \to \mathbb{R}\times\mathbb{R}, \quad R_d^Q : T^*Q \to Q\times Q.\] This decomposition allows us to encode the evolution of both the configuration variables and the time variable within a unified discrete geometric framework. This composition, together with the symplectic identifications described above, leads to a Lagrangian submanifold in \(T^*T^*(Q \times \mathbb{R})\) that implicitly defines the discrete evolution map on the extended phase space. The next result shows that the discrete flow constructed through this geometric framework defines a canonical transformation on \(T^*(Q \times \mathbb{R})\) whose projection onto \(T^*Q \times \mathbb{R}\) preserves the geometric properties (volume preservation, Poisson bracket...) associated to the continuous cosymplectic flow.

Proof: Let \(M = Q \times \mathbb{R}\) and consider the extended Hamiltonian \[\tilde{H} = H\circ\mu + p_t, \qquad \mu(q,t,p_q,p_t) = (q, p_q, t).\] We show that the submanifold \[\Lambda_h = (\Phi^{-1}\circ \widehat R_d^M \circ \mathcal{I}_{TM})(h\,d\tilde{H}) \subset T^*M \times T^*M\] is Lagrangian, and that it implicitly defines a symplectic map \(\Psi_h:T^*M\to T^*M\) that is \(\mu\)–projectable.

In local coordinates \((q,t,p_q,p_t)\) on \(T^*M\), the differential of \(\tilde{H}\) is \[d\tilde{H} = \left( \frac{\partial H}{\partial q}, \frac{\partial H}{\partial t}, \frac{\partial H}{\partial p_q}, 1 \right).\] Multiplying by \(h\) yields the discrete analogue of the Hamiltonian 1–form: \[h\,d\tilde{H} = \bigl( h\partial_q H,\; h\partial_t H,\; h\partial_{p_q}H,\; h \bigr).\] This is the cotangent vector that encodes a single time step of the discretized dynamics.

The canonical antisymplectomorphism \(\mathcal{I}_{TM}:T^*T^*M \longrightarrow T^*TM\) is locally given by \[\mathcal{I}_{TM}(x,p;\mu,\tilde{\mu}) = (x,\tilde{\mu};-\mu,p)\] and satisfies \(\mathcal{I}_{TM}^*\omega_{TM}=-\omega_{T^*M}\).

Applying \(\mathcal{I}_{TM}\) to \(h\,d\tilde{H}\) produces an element in \(T^*TM\): \[\mathcal{I}_{TM}(h\,d\tilde{H}) = \bigl( q,\;t,\; h\partial_{p_q}H,\; h;\; -h\partial_q H,\; -h\partial_t H,\; p_q,\; p_t \bigr).\]

Let \(R_d = R_d^{\mathbb{R}}\times R_d^Q : TM \to M\times M\). Its cotangent lift decomposes as \(\widehat R_d^M = \widehat R_d^{\mathbb{R}} \times \widehat R_d^Q\). In local coordinates \((q,t,\dot{q},\dot{t})\), \[R_d(q,t,\dot{q},\dot{t}) = \bigl( (R_d^{Q})^1(q,\dot{q}),\; (R_d^{\mathbb{R}})^1(t,\dot{t}),\; (R_d^{Q})^2(q,\dot{q}),\; (R_d^{\mathbb{R}})^2(t,\dot{t}) \bigr),\] and its cotangent lift transports momenta via \((T R_d)^{-1}\).

Applying \(\widehat R_d^M\) to \(\mathcal{I}_{TM}(h\,d\tilde{H})\) gives a point in \(T^*(M\times M)\). The map \(\Phi:T^*M\times T^*M\longrightarrow T^*(M\times M)\) is a symplectomorphism whose inverse reorganizes cotangent coordinates to produce pairs \((x,p_x;\; x',p_{x'})\).

Thus \[\Lambda_h = (\Phi^{-1}\circ\widehat R_d^M\circ\mathcal{I}_{TM}) (h\,d\tilde{H}) \subset T^*M\times T^*M.\]

Since the 1-form \(h\,d\tilde{H}\) is exact and both \(\mathcal{I}_{TM}\) and \(\Phi\) are symplectic or antisymplectic, the image \(\Lambda_h\) of \((\Phi^{-1}\circ\widehat R_d^M\circ\mathcal{I}_{TM})\) is by construction a Lagrangian submanifold of the product symplectic manifold \(T^*M\times T^*M\). It therefore encodes a symplectic map \(\Psi_h:T^*M\to T^*M\).

We must now show that: \[\mu\circ\Psi_h = \phi_h\circ\mu \quad\text{for some map }\phi_h:T^*Q \times \mathbb{R}\to T^*Q \times \mathbb{R}.\]

To show that \(\Psi_h\) is \(\mu\)–projectable, we first investigate the explicit coordinate expression of the Lagrangian submanifold \(\Lambda_h\). Using the coordinate formulas for \(\widehat R_d^M\) and \(\Phi^{-1}\) obtained earlier, we have \[\begin{align} (\Phi^{-1}\circ \widehat{R^M_d})(t, q, \dot{t}, \dot{q}; \, & p_t, p_q, p_{\dot{t}}, p_{\dot{q}})\\ &= \begin{pmatrix} (R^{\mathbb{R}}_d)^1(t, \dot{t}),\; -(\mathrm{pr}^{T^*\mathbb{R}}_1\circ T^*_{(t, \dot{t})}{(R^{\mathbb{R}}_d)^{-1}})(p_t, p_{\dot{t}})\\[0.25em] (R^{Q}_d)^1(q, \dot{q}),\; -(\mathrm{pr}^{T^*Q}_1\circ T^*_{(q, \dot{q})}{(R^{Q}_d)^{-1}})(p_q, p_{\dot{q}})\\[0.25em] (R^{\mathbb{R}}_d)^2(t, \dot{t}),\; (\mathrm{pr}^{T^*\mathbb{R}}_2\circ T^*_{(t, \dot{t})}{(R^{\mathbb{R}}_d)^{-1}})(p_t, p_{\dot{t}})\\[0.25em] (R^{Q}_d)^2(q, \dot{q}),\; (\mathrm{pr}^{T^*Q}_2\circ T^*_{(q, \dot{q})}{(R^{Q}_d)^{-1}})(p_q, p_{\dot{q}}) \end{pmatrix}. \end{align}\]

Substituting the value of \({\mathcal{I}}_{T(Q \times \mathbb{R})}(h\,d\tilde{H})\), namely \[(t,q,\; \dot{t},\dot{q};\; p_t,p_q,p_{\dot{t}},p_{\dot{q}}) = \bigl( t,q,\; h,\; h\partial_{p_q}H;\; -h\partial_tH,\; -h\partial_qH,\; p_t,\; p_q \bigr),\]

we obtain the explicit coordinate expression of \(\Lambda_h\): \[\Lambda_h= \begin{pmatrix} (R^{\mathbb{R}}_d)^1(t,h),\; -(\mathrm{pr}^{T^*\mathbb{R}}_1\circ T^*_{(t,h)}(R^{\mathbb{R}}_d)^{-1}) \bigl(-h\partial_t H,\; p_t \bigr) \\[0.4em] (R^{Q}_d)^1(q,h\partial_{p_q}H),\; -(\mathrm{pr}^{T^*Q}_1\circ T^*_{(q,h\partial_{p_q}H)}(R^{Q}_d)^{-1}) \bigl(-h\partial_q H,\; p_q \bigr) \\[0.4em] (R^{\mathbb{R}}_d)^2(t,h),\; (\mathrm{pr}^{T^*\mathbb{R}}_2\circ T^*_{(t,h)}(R^{\mathbb{R}}_d)^{-1}) \bigl(-h\partial_t H,\; p_t \bigr) \\[0.4em] (R^{Q}_d)^2(q,h\partial_{p_q}H),\; (\mathrm{pr}^{T^*Q}_2\circ T^*_{(q,h\partial_{p_q}H)}(R^{Q}_d)^{-1}) \bigl(-h\partial_q H,\; p_q \bigr) \end{pmatrix}.\]

From these expressions we see that the update of the configuration variables \(q\) and the momentum \(p_q\) depends on \(q,p_q,t\) and the Hamiltonian derivatives, but no term involving \(p_t\) appears in the \(Q\)-component updates. In particular, \[q_{k+1} = (R_d^Q)^2\bigl(q,\; h\partial_{p_q}H\bigr), \qquad p_{q,k+1} = (\mathrm{pr}^{T^*Q}_2\circ T^*R_d^Q)(-h\partial_q H,\; p_q),\] contains no \(p_t\). In addition, the variable \(p_t\) appears only inside the cotangent components of the \(\mathbb{R}\)-part, i.e., only in the expressions updating \(p_{t,k}\) and \(p_{t,k+1}\). Thus, the dependence of the discrete dynamics on \(p_t\) is vertical and does not couple into the \(Q\)–sector.

Hence, the full discrete map necessarily factorizes as \[\Psi_h(q,t,p_q,p_t) = \bigl( \phi_h(q,t,p_q),\; p_t + F(q,t,h,p_q) \bigr),\] for some smooth function \(F\) that does not depend on \(p_t\).

Applying \(\mu\) to both sides, \[\mu\circ\Psi_h(q,t,p_q,p_t) = \phi_h(q,t,p_q).\] Thus \(\Psi_h\) is \(\mu\)-projectable.\(\square\)

As a concrete example, consider \(M=\mathbb{R}^n\) and apply the midpoint discretization map \[R_d^M(q,t,\dot{q},\dot{t}) = \bigl(q-\tfrac{\dot{q}}{2}\,;\;t-\tfrac{\dot{t}}{2},\; q+\tfrac{\dot{q}}{2}\,;\; t+\tfrac{\dot{t}}{2} \bigr),\] which corresponds to evaluating the discrete step at the midpoint of the segment determined by \((q,t)\) and the tangent vector \((\dot{q},\dot{t})\). Substituting into the general construction of Theorem 1, the Lagrangian submanifold \(\Lambda_h\) yields the following implicit discrete update for the extended Hamiltonian system.

Indeed, since \[\dot{t} = h,\qquad \dot{q} = h\,\frac{\partial H}{\partial p}(q,t,p),\] and the cotangent components are transformed by the pullback of \((R_d^M)^{-1}\), the coordinate expression of \(\Lambda_h\) leads to the midpoint evaluation \[(q_{k+\frac{1}{2}}, t_{k+\frac{1}{2}}, p_{k+\frac{1}{2}}) = \Bigl( \tfrac{q_k+q_{k+1}}{2},\; \tfrac{t_k+t_{k+1}}{2},\; \tfrac{p_k+p_{k+1}}{2} \Bigr).\] Hence the discrete Hamilton equations produced by the geometric construction read: \[\begin{align} \frac{q_{k+1}-q_k}{t_{k+1}-t_k} &= \frac{\partial H}{\partial p}\bigl( q_{k+\frac{1}{2}},\; t_{k+\frac{1}{2}},\; p_{k+\frac{1}{2}} \bigr),\\[0.5em] \frac{p_{k+1}-p_k}{t_{k+1}-t_k} &= -\frac{\partial H}{\partial q}\bigl( q_{k+\frac{1}{2}},\; t_{k+\frac{1}{2}},\; p_{k+\frac{1}{2}} \bigr), \end{align}\] which is precisely the implicit midpoint rule for the non-autonomous Hamiltonian system.

Because this method arises from a Lagrangian submanifold defined through a discrete canonical transformation on the extended phase space \(T^*(Q \times \mathbb{R})\), it automatically inherits the geometric preservation properties established in the next corollary: volume preservation on \(T^*Q \times \mathbb{R}\), preservation of the Poisson bracket \(\{\cdot,\cdot\}_0\), and symplecticity on each time slice \(T^*Q \times \{t\}\). This is reflected with the following corollary.

4 Numerical Simulations↩︎

In this section we present two numerical examples to illustrate the performance and qualitative behavior of the proposed geometric integrators. These examples are chosen because they combine simplicity of interpretation with nontrivial dynamical features arising from dissipative effects and time-dependent potentials. Moreover, these are prototypical models in mechanical systems where energy dissipation and external potential play a central role, and the preservation of geometric structures (such as volume or Poisson structure) can significantly influence the accuracy and long-term stability of numerical simulations. The simulations are carried out using an implicit symplectic Euler-type geometric integrator, which can be schematized as follows \[\label{eq:32euler-type} \begin{align} q(i+1) &= q(i) + h \tilde{H}_p(q(i), t(i), p(i+1), p_{t}(i+1) ),\\ t(i+1) &= t(i) + h \tilde{H}_{p_t}(q(i), t(i), p(i+1), p_{t}(i+1) ),\\ p(i+1) &= p(i) - h \tilde{H}_q (q(i), t(i), p(i+1), p_{t}(i+1) ),\\ p_{t}(i+1) &= p_{t}(i) - h \tilde{H}_{t}(q(i), t(i), p(i+1), p_{t}(i+1) ). \end{align}\tag{3}\] for an arbitrary Hamiltonian \(\tilde{H}: T^*(Q\times \mathbb{R})\rightarrow \mathbb{R}\). These symplectic integrators are compared against a standard fourth-order Runge–Kutta (RK4) scheme. The RK4 method provides a reference in terms of local accuracy, but it does not preserve the underlying geometry. The comparison, therefore, highlights the qualitative differences in the long-term behavior.

4.1 Damped harmonic oscillator↩︎

We first consider a damped harmonic oscillator, a classical mechanical system used to describe oscillations under the influence of linear friction. This case provides a simple yet meaningful benchmark to test the consistency of the proposed discretization in the presence of an explicit time dependence in the Lagrangian. The formulation follows the approach in [19], where the damping is incorporated through a time-exponential factor that modifies both kinetic and potential terms. The Lagrangian of the system is given by \[\begin{align} L(t,q,\dot{q}) = e^{-\gamma t} \left( \frac{1}{2} m\dot{q}^2 - \frac{1}{2} m \omega_0 q^2 \right) \end{align}\] where \(m\) is the mass, \(\omega_0\) the undamped natural frequency, and \(\gamma\) the damping constant. The exponential factor \(e^{-\gamma t}\) models the continuous decay of energy due to frictional losses, making the system explicitly time-dependent and non-conservative. The corresponding extended Hamiltonian is \[\begin{align} H(q, t, p, p_t) = \frac{1}{2}e^{\gamma t} \frac{p^2}{m} + \frac{1}{2}e^{-\gamma t} m \omega_0 q^2 + p_t. \end{align}\] The dynamic equations for the system are \[\begin{align} \dot{q} = \frac{1}{m} e^{\gamma t} p, \quad \dot{t} = 1, \quad \dot{p} = -m \omega_0 e^{-\gamma t} q, \quad \dot{p}_t = -\frac{\gamma}{2m} e^{\gamma t} p^2 + \frac{\gamma m \omega_0}{2} e^{-\gamma t} q^2. \end{align}\] For these simulations, we set the parameters \(m = 1 \,\mathrm{kg}\), \(\omega_0 = 1 \,\mathrm{rad/s}\), \(\gamma = 0.2\); and the initial conditions \(q(0) = 1 \,\mathrm{rad}\), \(p(0) = 1 \,\mathrm{kg \, m^2 \, rad/s}\). Both integrators use a time step of \(T_s = 0.1 \,\mathrm{s}\) over a total integration horizon of \(T_{\text{fin}} = 50 \,\mathrm{s}\).

The time evolutions of the extended Hamiltonians are shown in Figure 4 (b). In Figure 4 (c) we show the evolution of the symplectic form \(\omega\) computed for the RK4 and the symplectic integrator. Two motions with close initial conditions are considered. The associated vector fields were approximated through finite differences between consecutive states of each trajectory \[X_{k+1}^i = \begin{bmatrix} q_{k+1}\\ p_{k+1}\\ \end{bmatrix}-\begin{bmatrix} q_{k}\\ p_{k}\\ \end{bmatrix}, \quad i = 1,2.\] The preservation of the discrete symplectic form was then verified by evaluating these vector fields on the symplectic form \(\omega (X_1, X_2)\). In Figure 4 (d) we present the evolution of the volume form \(\mathrm{vol}_0\). The comparison includes the proposed RK4 scheme and the symplectic integrator. For this analysis, three trajectories were integrated starting from slightly perturbed initial conditions. The associated vector fields were approximated through finite differences between consecutive states of each trajectory \[X_{k+1}^i = \begin{bmatrix} q_{k+1}\\ t_{k+1}\\ p_{k+1}\\ \end{bmatrix}-\begin{bmatrix} q_{k}\\ t_{k}\\ p_{k}\\ \end{bmatrix}, \quad i = 1,2,3.\] The discrete volume form is then evaluated as the determinant of the matrix composed of these approximate vector fields \[\mathrm{vol}_0 = \det \begin{bmatrix} X^1 \, X^2 \, X^3 \end{bmatrix}.\] For these simulations, we change the damping constant and final time horizon, choosing \(\gamma = 2 \cdot 10^{-4}\) and \(T_{\text{fin}} = 5 \cdot 10^5 \,\mathrm{s}\).

4.2 Harmonic oscillator with time-dependent potential↩︎

We consider a harmonic oscillator with a time-dependent potential energy as in [5]. In this example the time dependency affects only the potential energy of the system and does not impact the conjugate momenta. The Hamiltonian of the system has the form \[H(q, t, p) = \frac{1}{2} \left( (1 + \varepsilon \sin{(\alpha t)}) q^2 + p^2 \right)\] where \(\alpha, \varepsilon \in \mathbb{R}\) are constant parameters. The extended Hamiltonian becomes \[H(q, t, p, p_t) = \frac{1}{2} \left( (1 + \varepsilon \sin{(\alpha t)}) q^2 + p^2 \right) + p_t\] with dynamic equations \[\dot{q} = p,\quad \dot{t} = 1,\quad \dot{p} = -(1 + \varepsilon \sin{(\alpha t)}) q,\quad \dot{p}_t = -\frac{1}{2} \alpha\varepsilon \cos{(\alpha t)} q^2.\] We set as parameters \(\varepsilon = 0.0035\), \(\alpha = 0.0008\), and as initial conditions \(q(0) = 1 \,\mathrm{rad}\), \(p(0) = 1 \,\mathrm{kg \, m^2 \, rad/s}\). For the time parameters, we choose \(T_s = 0.1 \, \mathrm{s}\), and \(T_{fin} = 50 \, \mathrm{s}\).

Figure 5 (b) shows a comparison between the evolutions of the extended Hamiltonians computed using RK4 and the symplectic integrator. In Figure 5 (c) and Figure 5 (d), we show the time evolution of the symplectic and volume forms computed as in the previous example. For these simulations, we consider a longer time horizon of \(T_{fin} = 2\cdot 10^5 \, \mathrm{s}\).

In both the examples 4.1 and 4.2 the numerical simulations produced the expected results. The symplectic numerical integrator was designed to preserve the geometric structure, and this is verified by observing the behaviours of the symplectic and volume forms. The latter, indeed, present a dynamic that on average tend to maintain a constant value, with small oscillations due to numerical approximations. The Runge-Kutta method, instead, despite its high order of integration, shows dynamics that tend to diverge. Concerning the extended Hamiltonian, it is sufficient to consider a shorter time horizon to observe that, even in this case, the Runge-Kutta method produces a slight bias that can lead to incorrect values. It is known, instead, that the symplectic method does not preserve the Hamiltonian value exactly, but produces oscillations around an average value that coincides with the correct one.

5 Application to an electric particle in an electromagnetic field↩︎

In this section, we show how the theory developed in this work aligns with gauge-invariant systems. We present a particular application of the structure-preserving integrators to the electromagnetic field. In particular, we regard the case of different gauge choices for the system, and we evaluate which quantities are preserved under transformations from one gauge to another.

We consider a charged particle moving in space \(Q := \mathbb{R}^3\) under the action of a time-dependent electromagnetic field described in terms of a scalar \(\phi(q,t)\) and a vector potential \(A(q,t) = (A_1(q,t),A_2(q,t),A_3(q,t))\)², where \(q=(x,y,z)\) represents the three-dimensional coordinates. The Lagrangian of the system is given by \[\begin{align} L(q,\dot{q},t) &=\frac{1}{2} m\norm{\dot{q}}^2+eA(q,t)\cdot\dot{q}-e\phi(q,t). \end{align}\] We stress that the potentials are not uniquely determined, but can be chosen and modified according to the so-called gauge transformations. Let \(\chi=\chi(q,t)\) be a smooth function on \(Q \times \mathbb{R}\). A gauge transformation defined by \[\label{eq:gauge-transform-potentials} A'(q,t)=A(q,t)+\nabla\chi(q,t), \qquad \phi'(q,t)=\phi(q,t)-\partial_t\chi(q,t)\tag{4}\] leaves the electromagnetic fields unchanged and correspondingly the Lagrangian changes by a total derivative (see, e.g. [20]), therefore the Euler-Lagrange equations remain unchanged.

5.1 Gauge transformations as canonical transformations↩︎

The key observation for the extended Hamiltonian formalism is that a gauge transformation induces a canonical transformation on the extended phase space. Let \[\label{eq:extended-hamiltonian-general} \widetilde{H}(q,p,t,p_t)=\frac{1}{2m}\norm{p-eA(q,t)}^2+e\phi(q,t)+p_t\tag{5}\] denote the extended Hamiltonian, where \(p = (p_x, p_y, p_z) = \nabla_{\dot{q}} L\) indicate the canonical momenta. The canonical symplectic form on the extended phase space is \[\label{eq:32canonical32form32electric} \omega = \mathrm{d} x \wedge \mathrm{d} p_x + \mathrm{d} y \wedge \mathrm{d} p_y + \mathrm{d} z \wedge \mathrm{d} p_z + \mathrm{d} t \wedge \mathrm{d} p_t.\tag{6}\] Given a gauge function \(\chi(q,t)\), we define a gauge-dependent function on the phase state \[\label{eq:gauge-canonical-map} G_\chi(q,p,t,p_t)=\bigl(q,\,p+e\nabla\chi(q,t),\,t,\,p_t+e\partial_t\chi(q,t)\bigr).\tag{7}\] The following proposition explains the relation between gauge and canonical transformations.

Proposition 4 states that it’s possible to switch from a gauge choice to another while preserving the canonical structure. We state here the discrete counterpart.

Proposition 5 clarifies that, rather than discretizing independently in two gauge-equivalent coordinate descriptions, one may transport the discrete canonical map by conjugation. This operation preserves the volume form, the Poisson bracket and hence the symplectic structure.

5.2 Gauge analysis↩︎

In this section, we present several simulations that empirically validate the results obtained in the previous section. We start by considering the following vector and scalar potentials (we will refer to this choice as ‘original gauge’) \[A(q,t) = \beta(t)\begin{bmatrix} x-y \\ y+x \\ 0 \end{bmatrix}, \quad \phi(q,t) = 0,\] where \(\beta(t)\) is a time-dependent function. With these potentials, the electric and magnetic fields are given by \[E(q,t) = -\dot{\beta}(t) \begin{bmatrix} x-y \\ x+y \\ 0 \end{bmatrix}, \quad B(q,t) = \beta(t) \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix},\] and thus the model is non-autonomous and electromagnetic. The Lagrangian of the system reads \[\label{eq:lagrangian-initial} L(x,y,\dot{x},\dot{y},t) = \frac{1}{2} m(\dot{x}^2+\dot{y}^2) +e\beta(t)((x-y)\dot{x} + (x+y)\dot{y})\tag{8}\] with canonical momenta \[\label{eq:canonical-momenta-initial} \begin{align} p_x =& \frac{\partial L}{\partial \dot{x}} = m\dot{x} + e\beta(t)(x-y), \\ p_y =& \frac{\partial L}{\partial \dot{y}} = m\dot{y} + e\beta(t)(x+y), \end{align}\tag{9}\] and the mechanical momenta \[\begin{align} \Pi_x := p_x- eA_1 =& p_x - e\beta(t)(x-y), \\ \Pi_y := p_y- eA_2 =& p_y - e\beta(t)(x+y). \end{align}\] The non-autonomous Hamiltonian is then \[\label{eq:hamiltonian-initial} H(x,y,p_x,p_y,t) = \frac{1}{2m}\left(p_x - e\beta(t)(x-y)\right)^2 + \frac{1}{2m}\left(p_y - e\beta(t)(x+y)\right)^2,\tag{10}\] and the associated extended Hamiltonian is \[\label{eq:extended-hamiltonian-initial} \widetilde{H}(x,y,p_x,p_y,t,p_t) = H(x,y,p_x,p_y,t)+p_t.\tag{11}\]

We consider the following choices of gauges, with their respective vector and scalar potentials (we drop their coordinates and time dependencies).

1. Coulomb, characterized by the null gradient of the vector potential, i.e., \(\nabla \cdot A_{\text{c}} = 0\): \[\begin{align} \chi_{\text{c}} = -x^2 \beta(t), \quad A_{\text{c}} = \beta(t)\begin{bmatrix} -x-y \\ y+x \\ 0 \end{bmatrix}, \quad \phi_{\text{c}} = x^2 \dot{\beta}(t); \end{align}\]

2. symmetric, that requests \(\partial_y (A_{\text{s}})_1 + \partial_x (A_{\text{s}})_2 = 0\): \[\begin{align} \chi_{\text{s}} = -\frac{1}{2}(x^2+y^2) \beta(t), \quad A_{\text{s}} = \beta(t)\begin{bmatrix} -y \\ x \\ 0 \end{bmatrix}, \quad \phi_{\text{s}} = \frac{1}{2}(x^2+y^2) \dot{\beta}(t); \end{align}\]

3. Landau, with vector potential linear in \(x\) in the second component and zero in the others, i.e., \((A_{\text{La}})_2 = f(t)x\), \((A_{\text{La}})_1 = (A_{\text{La}})_3 = 0\): \[\begin{align} \chi_{\text{La}} = -\frac{1}{2}(x-y)^2\beta(t), \quad A_{\text{La}} = \beta(t)\begin{bmatrix} 0 \\ 2x \\ 0 \end{bmatrix}, \quad \phi_{\text{La}} = \frac{1}{2}(x-y)^2 \dot{\beta}(t); \end{align}\]

4. Lorenz, with scalar and vector potential invariant under Lorez transformations, i.e., \(\nabla \cdot A_{\text{Lo}} + \partial_t \phi_{\text{Lo}} = 0\): \[\begin{align} \chi_{\text{Lo}} = - 2 \int^t_0 \beta(\tau) d \tau, \quad A_{\text{Lo}} = \beta(t)\begin{bmatrix} x-y \\ y+x \\ 0 \end{bmatrix}, \quad \phi_{\text{Lo}} = 2 \beta(t). \end{align}\]

The corresponding non-autonomous extended Hamiltonian are \[\begin{align} \tilde{H}_{\mathrm{c}}(x,y,p_x,p_y,t,p_t) =& \frac{1}{2m}\left(p_x - e\beta(t)(-x-y)\right)^2\\ &+ \frac{1}{2m}\left(p_y - e\beta(t)(x+y)\right)^2 + ex^2 \dot{\beta}(t) + p_t, \\ \tilde{H}_{\mathrm{s}}(x,y,p_x,p_y,t,p_t) =& \frac{1}{2m}\left(p_x + e\beta(t)y\right)^2 + \frac{1}{2m}\left(p_y - e\beta(t)x\right)^2\\ &+ \frac{e}{2}(x^2+y^2) \dot{\beta}(t) + p_t, \\ \tilde{H}_{\mathrm{La}}(x,y,p_x,p_y,t,p_t) =& \frac{1}{2m}\left(p_x\right)^2+ \frac{1}{2m}\left(p_y - 2e\beta(t)x\right)^2\\ &+ \left(xy - \frac{1}{2}(x^2+y^2)\right) \dot{\beta}(t) + p_t,\\ \tilde{H}_{\mathrm{Lo}}(x,y,p_x,p_y,t,p_t) =& \frac{1}{2m}\left(p_x - e\beta(t)(x-y)\right)^2\\ &+ \frac{1}{2m}\left(p_y - e\beta(t)(x+y)\right)^2 + 2e \beta(t) + p_t. \end{align}\]

5.3 Numerical experiments↩︎

We compare the results obtained by applying the symplectic integrator to Hamiltonian systems with the different gauge choices. We consider the vector and scalar potentials defined in Section 5.2 with the function \(\beta\) defined as \[\beta(t) = B_0 (1 + \varepsilon \sin(\omega t)).\] We choose the parameters \(B_0 = 0.05\), \(\varepsilon = 10\), \(\omega = 2\), the mass and the charge are set to \(m = 1 \mathrm{kg}\), \(e = 1 \mathrm{C}\), while for the initial conditions, we put \(x(0) = 0.1 \mathrm{m}\), \(y(0) = 0.4 \mathrm{m}\), \(p_x(0) = 0.5 \mathrm{kg \, m/s}\), and \(p_y(0) = -0.2 \mathrm{kg \, m/s}\). We implement the Euler-type integrator as defined in 3 with an integration step \(T_s = 0.001 \mathrm{s}\), and for a time horizon of \(T_{fin} = 20 \mathrm{s}\). The trajectory simulations for the different gauge choices are shown in Figure 6 (a). Here, we compare the trajectories for \(x\) and \(y\), the canonical momenta \(p_x\) and \(p_y\), and the mechanical momenta \(\Pi_x\) and \(\Pi_y\). The experimental results show the behaviors that we expect. The trajectories on the \(xy\) plane are the same for all the choices of the gauge, as well as the mechanical momenta \(\Pi_x\) and \(\Pi_y\), with a small deviation depending on the integration order. Figure 6 (b) highlights these differences by showing how the various trajectories, obtained with different gauge choices, differ from the original one, as prescribed by the Remark 2. Moreover, we can observe clearly different behaviors regarding canonical momenta as highlighted by Remark 1.

6 Conclusions and Future Work↩︎

In this paper we have presented a geometric framework for discrete canonical transformations in the non-autonomous Hamiltonian setting. Starting from the (presymplectic) cosymplectic formulation on \(T^*Q \times \mathbb{R}\), we introduced a generalized canonical transformation on the extended phase space \(T^*(Q \times \mathbb{R})\) and showed that its projection acts as a Poisson morphism. This construction naturally induces a discrete flow that preserves the cosymplectic volume, the Poisson bracket \(\{\cdot,\cdot\}_0\), and the symplectic form on each time fiber \(T^*Q \times \{t\}\). The resulting integrator provides a geometric and structure-preserving alternative to traditional discretization schemes for non-autonomous systems. Numerical simulations confirm that the proposed method avoids artificial energy drift and accurately reproduces the geometric structure of the phase space over long integration times.

Moreover, by expressing the discrete map as the pullback of the canonical symplectic form on \(T^*(Q \times \mathbb{R})\), we established a global, coordinate-free criterion for structure preservation that generalizes the local, matrix-based approach of Marthinsen and Owren [5] and extends it to the cosymplectic setting.

Future work will first explore the application of the proposed integrator to the discretization of geometric filters and observers on Lie groups, particularly those based on the second-order minimum-energy framework developed in [21] and further explored in [22]–[25]. In that formulation, the underlying estimation problem is a time-dependent optimal control problem on the tangent or cotangent bundle of a Lie group, where the Hamiltonian depends explicitly on time through the control and measurement processes. This makes the filtering dynamics inherently non-autonomous, and therefore naturally suited to the cosymplectic framework introduced in this paper. By discretizing the extended Hamiltonian system associated with the filtering problem, we aim to obtain discrete update maps that remain symplectic (or cosymplectic) on each time slice, thereby preserving group invariance, geometric consistency, and energy properties of the estimation dynamics. This approach would yield a new class of structure-preserving recursive filters on Lie groups, combining the theoretical rigor of non-autonomous Hamiltonian discretization with the practical robustness required for real-time estimation.

Beyond this direct application, the discrete geometric framework developed here could also be extended toward variational formulations of estimation, leading to discrete variational filters derived from a discrete action principle. Such extensions would provide a unified variational–geometric approach to filtering and control on manifolds, potentially integrating estimation, prediction, and learning within the same symplectic–cosymplectic framework.

References↩︎