1 Introduction↩︎

In Numerical Analysis, Geometric Numerical Integration refers to the construction of algorithms that approximate the solution of ordinary differential equations while preserving the geometric characteristics of the given problem [1]. When the differential equations arise as equations of motion of mechanical systems there is a well known way of constructing geometric integrators using what are known as Discrete Mechanical Systems. The solution of the equation of motion of a (continuous) mechanical system, a trajectory, can be seen as a critical point of a variational problem in a space of paths. Similarly, trajectories of a discrete mechanical system are critical points of a certain functional —defined in a space of discrete paths that, crucially, is finite dimensional— and are characterized by equations of motion that are algebraic. Solving these equations gives rise to a numerical integrator of the original differential problem (see [2]), known as a variational integrator, provided that

1. the discrete mechanical system is “close enough” to the continuous one, and

2. that [it:goal-approximation] suffices to conclude that the trajectories of the systems are “close enough”.

Also of importance,

1. how well do variational integrators preserve the geometric characteristics of the system?

The description of trajectories of mechanical systems defined on a configuration space \(Q\) as critical points of a functional is characteristic of the Lagrangian formulation of Mechanics —variational formulation would be a better name—, where the functional, called the action, is computed using a Lagrangian function \(L\) over the tangent bundle \(TQ\). Alternatively, it is possible to give a characterization of trajectories as critical points of a functional on curves in the cotangent bundle \(T^*Q\), computed using a Hamiltonian function \(H\) on \(T^*Q\) (see [3] or [4]). In most cases, the two descriptions are equivalent. Both approaches have discrete versions: by far, the most common is the Lagrangian approach, where the discrete action is defined using a discrete Lagrangian function \(L_d:Q\times Q\rightarrow \mathbb{R}\) (see [5],  [2] and [6]). The discrete Hamiltonian approach considers discrete paths in \(T^*Q\) and the corresponding action functional is constructed using a discrete Hamiltonian function \(H_d:T^*Q\rightarrow \mathbb{R}\) (see [7] and [8]). These variational integrators are known to satisfy points [it:goal-approximation] and [it:goal-close95implies95flow95close] and, regarding [it:goal-properties], preserve some natural symplectic structures whereas, if symmetry is present, the associated momenta are conserved (see [2], [9] and [10]).

Unfortunately, in many real world applications, it is necessary to consider mechanical systems that are subjected to external forces. These systems rarely have the nice conservation properties mentioned above. Still, the forced discrete Lagrangian systems have been successfully used and studied for some time (see [2], [11], [9], and [12]).

On the other hand, the study of forced discrete Hamiltonian systems lags behind. The purpose of this paper is to study such systems, with focus on the case where the configuration manifold \(Q\) is a finite dimensional real vector space. We define trajectories of a forced discrete Hamiltonian system (FDHS) as the solutions of a discrete variational problem, similar to what is done in Section 3.2 of [13], and, also, extending the idea used in Section 3 of [8] to the forced case. Actually, for a technical reason, we introduce a slightly more general notion that we call extended trajectory of the system. Trajectories, extended or plain, are, also, characterized as solutions of a system of algebraic equations. As in the Lagrangian case, we introduce a notion of forced discrete Legendre transformation —in fact, two of them, a \(+\) and a \(-\) version— and call regular the FDHSs where they are local diffeomorphisms. We prove that, if the space of extended trajectories of length \(N\) is non-empty, it is an open set in \(Q^{N+1}\).

We also study some of the structural properties of FDHSs. As expected, we see that the canonical symplectic form \(\omega_Q\) on \(T^*Q\) is not, in general, preserved by the flow; an exception is the case when the forces are closed. Similarly, if a Lie group is a symmetry group of the system, we obtain a formula describing the evolution of the corresponding canonical momentum map and give a condition that ensures that it is conserved. These properties are the current, partial, answer to the point [it:goal-properties] above.

Even though the theory that we develop in this paper is purely Hamiltonian, we see that given a “good” Lagrangian system (meaning hyperregular and satisfying a certain “modified hyperregularity condition”) it is possible to construct an FDHS such that there is a bijective correspondence between the trajectories of the former system and the extended trajectories of the latter.

A central part of the paper is the error analysis of the variational integrators constructed using FDHSs. We first prove that, given a (continuous) forced Hamiltonian system \(\mathcal{M}\), there is an \(h\)-dependent family of FDHSs \(\mathcal{M}_{d,h}^E\) (\(h\) is a scalar parameter defined in a neighborhood of \(0\)) such that the trajectories of \(\mathcal{M}_{d,h}^E\) are the trajectories of \(\mathcal{M}\) evaluated at times that are multiples of \(h\), and conversely. We call \(\mathcal{M}_{d,h}^E\) the exact FDHS associated to \(\mathcal{M}\). Unfortunately, this system \(\mathcal{M}_{d,h}^E\) has no direct practical use because it cannot be computed effectively in most real cases. The true interest in \(\mathcal{M}_{d,h}^E\) comes from using approximations \(\mathcal{M}_{d,h}\) of \(\mathcal{M}_{d,h}^E\), usually thought of as discretizations of \(\mathcal{M}\). The main result we prove in this area is that if \(\mathcal{M}_{d,h} = \mathcal{M}_{d,h}^E + \mathcal{O}(h^{r+1})\) for some \(r\in\mathbb{N}\) (notation to be explained), then the corresponding flows satisfy \(\mathbf{F}^{{\mathcal{M}}}_{{h}} = \mathbf{F}^{{\mathcal{M}_{d,h}^E}} = \mathbf{F}^{{\mathcal{M}_{d,h}}} + \mathcal{O}(h^{r+1})\). Thus, if we choose a discretization of \(\mathcal{M}\) that is accurate to order \(r\), the corresponding variational integrator has, at least, the same order of accuracy (local error of the same order). This is the answer to the point [it:goal-close95implies95flow95close].

We also touch on the practical matter of constructing discretizations for a given (continuous) forced Hamiltonian system \(\mathcal{M}\). We propose a method based on expanding \(\mathcal{M}_{d,h}^E\) as a Taylor series around \(h=0\); we provide explicit formulas for the expansion up to orders \(1\) and \(2\). An alternative method using Gaussian quadrature and the shooting method for integrating boundary value problems for ODEs is discussed. Anyone of these methods provides a concrete way to satisfy [it:goal-approximation].

The paper is structured as follows: Section 2 reviews some notions and results on (continuous) forced Hamiltonian systems from the variational point of view. Forced discrete Hamiltonian systems and their dynamics are introduced in Section 3. Section 4 analyzes the evolution of the symplectic form and momentum by the flow of an FDHS; it also discusses the construction of an FDHS from a given forced Lagrangian system, such that the trajectories of the two systems are in bijective correspondence. Sections 5 and 6 are dedicated to the error analysis of the corresponding variational integrators: the former focuses on the exact FDHS while the latter proves that a discretization of order \(r\) leads to an integrator of order, at least, \(r\). Last, Section 7 discusses two methods that can be used to construct FDHSs in practice.

2 Forced Hamiltonian systems: variational approach↩︎

Since we are looking for a variational version of forced discrete Hamiltonian system, we will work, following [8], on a real vector space \(Q\). Hence, its cotangent bundle can be trivialized as \(T^*Q \simeq Q \times Q^*\) and expressions such as \(p \dot{q}\) or \((q,p)\) for a curve in it are adequate.

We begin with a brief review of the variational formulation of Hamiltonian mechanics.

The equations ?? are called forced Hamilton equations (see, for example, [2, Sec. 3.1.2]).

The previous formulation works for trajectories whose initial and final positions are known. There are other formulations in which the known data are \(q(0)\) and \(p(T)\) or \(p(0)\) and \(q(T)\). In this direction, and inspired by [8], we have the following definition.

We close this section with a result that proves that, locally, solutions to the forced Hamilton equations ?? with boundary conditions do exist. We assume that \(Q=\mathbb{R}^n\) or, more generally, \(Q\subset \mathbb{R}^n\) is an open subset. Of course, this means no loss of generality as this identification corresponds to the choice of a basis in the \(\mathbb{R}\)-vector space \(Q\). Also, for \(x\in\mathbb{R}^n\), we use the norm \(\lVert x\rVert_\infty\mathrel{\vcenter{:}}=\max\{\left\lvert{x_1}\right\rvert,\ldots, \left\lvert{x_n}\right\rvert\}\) and denote the corresponding balls by \(B^\infty\).

3 Forced discrete Hamiltonian systems↩︎

Given a forced Hamiltonian system \(\mathcal{M}\mathrel{\vcenter{:}}=(Q,H,\phi)\) and \(h > 0\), we introduce a notion of discrete Hamiltonian as an approximation \[H_d(q_0,p_1) \approx p(h) q(h) - \int_0^h (p(t) \dot{q}(t) - H(q(t),p(t))) \;dt,\] where \((q(t),p(t)) \mathrel{\vcenter{:}}=\mathbf{F}^{{\mathcal{M}}}_{{h},{t}}(q_0,p_1)\) is the trajectory of \(\mathcal{M}\) that satisfies the boundary conditions \(q(0) = q_0\) and \(p(h) = p_1\). Similarly, the discrete force is an approximation \[\phi_d(q_0,p_1)(\delta q_0, \delta p_1) \approx \int_0^h \phi(q(t),p(t)) T_{(q_0,p_1)} \mathbf{F}^{{\mathcal{M}}}_{{h},{t}}(\delta q_0, \delta p_1) \;dt.\]

In what follows, unless explicitly stated otherwise, we will consider systems with type II trajectories. Thus, we will drop the “type II” in the name.

A similar approach to the variational principle used in Definition 4 has, also, been considered in [13]. A difference between the two approaches is that we consider \(q_\cdot\) and \(p_\cdot\) to be independent, whereas in the cited work, they are functionally related. This is reflected in the difference between our equations ?? and their (4a).

Notice that the equations ?? that characterize the trajectories of an FDHS do not involve neither \(p_0\) nor \(q_N\). We could use a notion of trajectory in which those points are absent, considering curves such as \(((q_0,p_1),\ldots,(q_{N-1},p_N))\), but we opt for keeping them and asking the points to satisfy an additional condition, as it is done in [8, Sec. 3.1].

For practical purposes, it is important to notice that the equations of motion for trajectories and extended trajectories of FDHSs can be decoupled (in \(k\)) and, so, can be solved iteratively instead of having to solve them concurrently. We have the following result, whose proof is a straightforward verification.

Algebraic equations do not always have solutions and the equations of motion of an FDHS ?? (or 4 ) are no exception, even for regular systems. The next result establishes that, given a (length-\(1\)) extended trajectory of an FDHS, a flow for the system can be defined near that trajectory.

We can depict the extended discrete trajectories as follows. \[\xymatrixcolsep{1.5pc}\xymatrix{ {} & {(q_0,p_1)} \ar@{|->}[dl]_{\mathbb{F}_{\phi_d}^-H_d} \ar@{|->}[dr]^{\mathbb{F}_{\phi_d}^+H_d} & {} & {(q_1,p_2)} \ar@{|->}[dl]_{\mathbb{F}_{\phi_d}^-H_d} & {\cdots} & {(q_{N-1},p_N)} \ar@{|->}[dr]^{\mathbb{F}_{\phi_d}^+H_d} & {}\\ {(q_0,p_0)} \ar@{|->}[rr]_{F^{\mathcal{M}_d}} & {} & {(q_1,p_1)} \ar@{|-}[r] & {} & {\cdots} & {} \ar@{->}[r] & {(p_N,q_N)} }\]

4 Structural properties↩︎

4.1 Relation with Lagrangian systems↩︎

It is well known that, under certain regularity conditions, there exists a relation between the trajectories of (continuous) Lagrangian and Hamiltonian systems. Specifically, given a Lagrangian system, the Legendre transform can be used to construct a Hamiltonian system whose trajectories are in correspondence with those of the original Lagrangian system. We begin this section observing how this is reflected in the discrete setting.

In this context, let us recall a few definitions. See Part 3 of [2] and [12] for additional information.

Notice that in 9 we are taking advantage of the Cartesian product structure of \(Q\times Q\) to decompose \(f_d = f_d^- + f_d^+\).

Notice that the forced discrete Euler–Lagrange equations 9 may be written as \[\label{eq:forced95EL95equations-ff} \mathbb{F}^+_{f_d}L_{d}(q_{k-1},q_{k}) = \mathbb{F}^-_{f_d}L_{d}(q_{k},q_{k+1}) \quad\text{ for }\quad k=1,...,N-1.\tag{10}\]

In what follows, we will use the second components of both transforms, so it will be useful to define the functions \[\begin{align} \overline{\mathbb{F}^+_{f_d}L_{d}}(q_0,q_1) &\mathrel{\vcenter{:}}=(\operatorname{pr}_2 \circ \mathbb{F}^+_{f_d}L_{d})(q_0,q_1) \in Q^* \\ \overline{\mathbb{F}^-_{f_d}L_{d}}(q_0,q_1) &\mathrel{\vcenter{:}}=(\operatorname{pr}_2 \circ \mathbb{F}^-_{f_d}L_{d})(q_0,q_1) \in Q^*, \end{align}\] where \(\operatorname{pr}_2 : Q \times Q^* \longrightarrow Q^*\) is the projection onto the second factor.

Let us consider the function \(\widetilde{\mathbb{F}}_{f_d}L_d : Q \times Q \longrightarrow Q \times Q^*\) given by \[\widetilde{\mathbb{F}}_{f_d} L_d(q,q') \mathrel{\vcenter{:}}= \left(q,\overline{\mathbb{F}_{f_d}^+ L_d}(q,q')\right).\]

It is easy to find examples that show that the regularity of a FDLS does not guarantee the satisfaction of the MRC.

Let \((Q,L_d,f_d)\) be a FDLS that satisfies the MHC. We define \(H_d : T^*Q \longrightarrow\mathbb{R}\) by \[\label{eq:H95d95from95L95d} H_d(q,p) \mathrel{\vcenter{:}}= p q^+ - L_d(q,q^+), \quad\text{where q^+ \in Q satisfies}\quad p = \overline{\mathbb{F}_{f_d}^+ L_d}(q,q^+).\tag{11}\]

Since \(\widetilde{\mathbb{F}}_{f_d} L_d\) is a diffeomorphism, \(q^+ \mathrel{\vcenter{:}}=(\operatorname{pr}_2 \circ (\widetilde{\mathbb{F}}_{f_d} L_d)^{-1})(q,p)\), so that \(H_d\) is well defined. In other words, we have a function \(q^+ : Q \times Q^* \longrightarrow Q\) that satisfies \((\widetilde{\mathbb{F}}_{f_d} L_d)^{-1} (q,p) = (q,q^+(q,p))\).

The other ingredient required to construct an FDHS associated to \((Q,L_d,f_d)\) is a \(1\)-form on \(T^*Q\). With this in mind, we consider \[\label{eq:phi95d95from95f95d} \phi_d \mathrel{\vcenter{:}}=((\widetilde{\mathbb{F}}_{f_d} L_d)^{-1})^* f_d.\tag{12}\]

Unraveling the definitions, we obtain \[\label{eq:phi95d95from95f95d-explicit} \begin{align} \phi_d^q(q,p)(\delta q) &= \left( f_d^-(q,q^+(q,p)) + f_d^+(q,q^+(q,p)) D_1 q^+ (q,p) \right) (\delta q) \\ \phi_d^p(q,p)(\delta p) &= f_d^+(q,q^+(q,p)) D_2q^+(q,p)(\delta p). \end{align}\tag{13}\]

In summary, we have constructed an FDHS \((Q,H_d,\phi_d)\) from a FDLS \((Q,L_d,f_d)\) that satisfies the MHC.

Next, we study the relation between the trajectories of both systems. Given a discrete curve \(q_\cdot \mathrel{\vcenter{:}}=(q_0,\ldots,q_N)\) in \(Q\), we define a discrete curve in \(T^*Q\) by \[\label{eq:curve95in95T9442Q95from95curve95in95Q} (q_k,p_k) \mathrel{\vcenter{:}}= \begin{cases} \mathbb{F}_{f_d}^+ L_d (q_{k-1},q_k), \quad\text{ if }\quad k = 1,\ldots,N,\\ \mathbb{F}_{\phi_d}^-H_d(q_0,p_1), \quad\text{ if }\quad k=0. \end{cases}\tag{14}\]

4.2 The canonical symplectic structure↩︎

We now study the evolution of the canonical symplectic structure of \(T^*Q\) by the flow of an FDHS. In [8], the authors study this for unforced systems, where one expects this structure to be conserved.

Given an FDHS \((Q,H_d,\phi_d)\), let us define the \(2\)-forms \(\omega_{\phi_d}^\pm \mathrel{\vcenter{:}}=(\mathbb{F}_{\phi_d}^\pm H_d)^* \omega_Q\), where \(\omega_Q\) is the canonical \(2\)-form on \(T^*Q\). In coordinates, using the decomposition ?? , \[\phi_d = \phi_d^q + \phi_d^p = \phi_{d,j}^q \;dq^j + \phi_d^{p,j} \;dp_j\] we have \[(\mathbb{F}_{\phi_d}^+ H_d)^*(q^i \;dp_i) = \left( \frac{\partial H_d}{\partial p_i} - \phi_d^{p,i} \right) \;dp_i\] and, therefore, \[\begin{align} (\mathbb{F}_{\phi_d}^+ H_d)^*(\omega_Q) &= \frac{\partial^2 H_d}{\partial q^j \partial p_i} \;dq^j \wedge dp_i + \frac{\partial^2 H_d}{\partial p_j \partial p_i} \;dp_j \wedge dp_i - d\phi_d^{p,i} \wedge dp_i \\ &= \frac{\partial^2 H_d}{\partial q^j \partial p_i} \;dq^j \wedge dp_i - d\phi_d^{p,i} \wedge dp_i. \end{align}\]

A similar computation shows that \[(\mathbb{F}_{\phi_d}^- H_d)^*(\omega_Q) = \frac{\partial^2 H_d}{\partial q^j \partial p_i} \;dq^j \wedge dp_i + d\phi_{d,j}^q \wedge dq^j.\]

Defining \[\omega_{H_d} \mathrel{\vcenter{:}}= \frac{\partial^2 H_d}{\partial q^j \partial p_i} \;dq^j \wedge dp_i,\] we have proved the following result:

Let \(\mathcal{M}_d\mathrel{\vcenter{:}}=(Q,H_d,\phi_d)\) be a regular FDHS. We are interested in the evolution of \(\omega_Q\) by the discrete flow \((q_k,p_k)\mapsto (q_{k+1},p_{k+1})\) of the system. If we denote this flow by \(\mathbf{F}^{{\mathcal{M}_d}}\), recalling 5 and 6 , we can express \(\mathbf{F}^{{\mathcal{M}_d}} = \left( \mathbb{F}_{\phi_d}^+ H_d \right) \circ \left( \mathbb{F}_{\phi_d}^- H_d \right)^{-1}\) and, using the previous computations, derive the following result.

4.3 The canonical momentum map↩︎

Let \(\mathcal{M}_d\mathrel{\vcenter{:}}=(Q,H_d,\phi_d)\) be an FDHS and let \(G\) be a Lie group acting on the left on \(Q\), by a free and proper action \(l^Q\). Let \(l^{Q^*}\) be the induced action on \(Q^*\). Following the ideas for the unforced case found in [8, Sec. 5], \(G\) is said to be a symmetry group for the system \((Q,H_d,\phi_d)\) if the function \(R_d : Q \times Q \times Q^* \longrightarrow\mathbb{R}\) given by \[R_d(q_0,q_1,p_1) \mathrel{\vcenter{:}}= p_1 q_1 - H_d(q_0,p_1)\] is \(G\)-invariant for the action \(l_g^{Q \times Q \times Q^*}(q_0,q_1,p_1) \mathrel{\vcenter{:}}= (l_g^Q(q_0),l_g^Q(q_1),l_g^{Q^*}(p_1))\).

Let \(J : T^*Q \longrightarrow\mathfrak g^*\) be the canonical momentum map, defined by \(J(q,p)\xi \mathrel{\vcenter{:}}= p(\xi_Q(q))\) for any \(\xi \in \mathfrak g\). Given an extended trajectory \(((q_0,p_0),(q_1,p_1))\) of \(\mathcal{M}_d\) and \(\varepsilon > 0\), let \(q_i^\varepsilon \mathrel{\vcenter{:}}= l^Q_{\exp(\varepsilon \xi)}(q_i)\) and \(p_i^\varepsilon \mathrel{\vcenter{:}}= l^{Q^*}_{\exp(\varepsilon \xi)}(p_i)\), with \(i = 0,1\). Since \(R_d\) is \(G\)-invariant, a computation similar to that of [8] yields \[\begin{align} 0 & = \frac{d}{d\varepsilon} \left( p_1^\varepsilon q_1^\varepsilon - H_d(q_0^\varepsilon,p_1^\varepsilon) \right) \bigg|_{\varepsilon=0} \\ & = p_1 \left. \frac{d}{d\varepsilon} q_1^\varepsilon \right|_{\varepsilon=0} + q_1 \left. \frac{d}{d\varepsilon} p_1^\varepsilon \right|_{\varepsilon=0} - D_1H_d(q_0,p_1) \left. \frac{d}{d\varepsilon} q_0^\varepsilon \right|_{\varepsilon=0} - D_2H_d(q_0,p_1) \left. \frac{d}{d\varepsilon} p_1^\varepsilon \right|_{\varepsilon=0} \\ & = p_1 \left. \frac{d}{d\varepsilon} q_1^\varepsilon \right|_{\varepsilon=0} + q_1 \left. \frac{d}{d\varepsilon} p_1^\varepsilon \right|_{\varepsilon=0} - \left(p_0 + \phi^q_d(q_0,p_1) \right) \left. \frac{d}{d\varepsilon} q_0^\varepsilon \right|_{\varepsilon=0} - \left(q_1 + \phi^p_d(q_0,p_1) \right) \left. \frac{d}{d\varepsilon} p_1^\varepsilon \right|_{\varepsilon=0} \\ & = p_1 \left. \frac{d}{d\varepsilon} q_1^\varepsilon \right|_{\varepsilon=0} - p_0 \left. \frac{d}{d\varepsilon} q_0^\varepsilon \right|_{\varepsilon=0} - \phi^q_d(q_0,p_1) \left. \frac{d}{d\varepsilon} q_0^\varepsilon \right|_{\varepsilon=0} - \phi^p_d(q_0,p_1) \left. \frac{d}{d\varepsilon} p_1^\varepsilon \right|_{\varepsilon=0} \\ & = p_1 \xi_Q(q_1) - p_0 \xi_Q(q_0) - \phi_d(q_0,p_1) \xi_{Q \times Q^*}(q_0,p_1). \end{align}\]

Rewriting the previous identity in terms of the canonical momentum map \(J\), we have the following result.

5 Error analysis: discrete exact systems↩︎

So far, we have considered forced discrete Hamiltonian systems as discrete-time dynamical systems, mostly independent of their, more common, continuous-time counterparts. In this section we want to see how, given a (continuous) forced Hamiltonian system \(\mathcal{M}\), there is a one parameter family of FDHSs \(\mathcal{M}_{d,h}^E\) —known as the discrete exact systems—, whose trajectories (for a given \(h\)) interpolate the trajectories of \(\mathcal{M}\) at discrete time steps (usually \(h\)).

It is not clear that the family of systems \(\mathcal{M}_{d,h}^E\) described in Example 9 is well defined because it relies on the existence of the flow associated to boundary value problems. Proposition 19 shows that, under certain conditions, \(\mathcal{M}_{d,h}^E\) is well defined.

Let \(\mathcal{M}\mathrel{\vcenter{:}}=(Q,H,\phi)\) be a forced Hamiltonian system and \((q_0,p_1)\in T^*Q = Q\times Q^*\). Then, by Proposition 5, there exist constants \(r,r'>0\) and \(T_-<T_+\) such that \(0\in (T_-,T_+)\) and, for any \(T\in [T_-,T_+]\) and \((q_0',p_1') \in \overline{B^\infty_{r'}(q_0,p_1)} \subset Q\times Q^*\), the boundary value problem ?? has a unique solution \((q(t),p(t))\) that satisfies \((q(0),p(0))\in B^\infty_r(q_0,p_1)\). Fix \(T\in [T_-,T_+]\); given \((q_0',p_1')\in B^\infty_{r'}(q_0,p_1)\), let \((q_T(t),p_T(t))\) be the unique solution of ?? such that \((q_T(0),p_T(0))\in B^\infty_r(q_0,p_1)\). With this information we define \[\label{eq:exact95discrete95hamiltonian95system-II-hamiltonian-local} H_{d,T}^E(q_0',p_1') \mathrel{\vcenter{:}}= p_T(T) q_T(T) - \int_0^T (p_T(t) \dot{q_T}(t) -H(q_T(t),p_T(t))) dt,\tag{20}\] and \[\label{eq:exact95discrete95hamiltonian95system-II-force-local} \phi_{d,T}^E(q_0',p_1')(\delta q_0',\delta p_1') \mathrel{\vcenter{:}}= \int_0^T \phi(q_T(t),p_T(t))T_{(q_0',p_1')} F^{\mathcal{M}}_t(\delta q_0',\delta p_1') dt,\tag{21}\] over \(B^\infty_{r'}(q_0,p_1)\).

Next, we explain the “exact” part of the name of \(\mathcal{M}_{d,h}^E\) introduced in Example 9. We first have a technical result involving forced discrete Legendre transforms and, then, the real goal, relating discrete trajectories of \(\mathcal{M}_{d,h}^E\) to the continuous ones of \(\mathcal{M}\).

6 Error analysis: approximations↩︎

Given a (continuous) forced Hamiltonian system \(\mathcal{M}\), we saw in Section 5 that the discrete exact system \(\mathcal{M}_{d,h}^E\) has the very attractive property of having its trajectories be the trajectories of \(\mathcal{M}\), evaluated at discrete times that are multiples of \(h\).

Still, the family \(\mathcal{M}_{d,h}^E\) is almost never constructible in practice. Thus, we consider approximations of \(\mathcal{M}_{d,h}^E\) by families of FDHSs \(\mathcal{M}_{d,h}\) —known as discretizations of \(\mathcal{M}\)— that, ideally, also approximate the trajectories of \(\mathcal{M}_{d,h}^E\) which, in turn, interpolate the trajectories of \(\mathcal{M}\).

6.1 Contact order of maps↩︎

Here we introduce the notion of contact order for families of maps between (open subsets of) finite-dimensional real vector spaces. In what follows, \(Q\) and \(Q'\) are two such spaces. For any map \(f:Q\times\mathbb{R}\rightarrow Q'\), we define \(\overline{f}:Q\times \mathbb{R}\rightarrow Q'\times \mathbb{R}\) by \(\overline{f}(q,h)\mathrel{\vcenter{:}}=(f(q,h),h)\). In what follows, we will consider smooth maps, with respect to the natural structures of smooth manifold on \(Q\) and \(Q'\). Also, \(\mathcal{W},\mathcal{W}_a\subset Q\times \mathbb{R}\) are open neighborhoods of \(Q\times \{0\}\).

The following result provides a convenient characterization for the contact order.

Proposition 23 is an application of Taylor’s Formula (with a parameter). See Lemma 2.24 in [9].

In the next two results, we use open neighborhoods \(\mathcal{W}_a\subset Q\times\mathbb{R}\) and \(\mathcal{W}'_a\subset Q'\times\mathbb{R}\) of \(Q\times\{0\}\) and \(Q'\times \{0\}\), respectively.

When \(\mathcal{W}\subset Q\times \mathbb{R}\) is an open subset, \(T\mathcal{W}= (TQ\oplus T\mathbb{R})|_\mathcal{W}\). Then, if \(f:\mathcal{W}\rightarrow Q'\) is smooth, we consider \(D_1f\) to be the restriction of \(Tf\) to \((TQ\oplus\{0\})|_\mathcal{W}\). As \((TQ\oplus\{0\})|_\mathcal{W}\simeq \mathcal{W}\times Q\) is an open subset in a vector space, we can apply our order of contact notions to \(D_1f\).

6.2 Contact order of systems↩︎

Here we extend the notion of contact order between maps to contact order between FDHSs.

The next result is the key to the error analysis of forced Hamiltonian systems.

6.3 Discretizations of FDHSs↩︎

Now we come back to the relationship between a family of forced discrete Hamiltonian systems and a given, continuous, one.

We can apply Theorem 27 to discretizations as follows.

7 Construction of forced discrete Hamiltonian systems↩︎

In this section, we consider the practical construction of forced discrete Hamiltonian systems \(\mathcal{M}_{d,h}\mathrel{\vcenter{:}}=(Q,H_{d,h},\phi_{d,h})\) that approximate a given forced Hamiltonian system \(\mathcal{M}\mathrel{\vcenter{:}}=(Q,H,\phi)\). There are different possible approaches. One of them relies on the approximations of orders \(1\) and \(2\) of the exact discrete system \(\mathcal{M}_{d,h}^E\) provided by Example 10. More precisely, defining \[\label{eq:H95d94E95and95phi95d94E-order951} H_{d,h}(q_0,p_1) = p_1q_0 + h H(q_0,p_1) \quad\text{ and }\quad \phi_{d,h}(q_0,p_1) = h \check{\phi}(q_0,p_1)\tag{28}\] we have a discretization of \(\mathcal{M}\) that is accurate to order \(1\), whereas defining \(H_{d,h}\) and \(\phi_{d,h}\) using 26 and 27 we have a discretization of \(\mathcal{M}\) that is accurate to order \(2\).

Next we describe an alternative, using the shooting method, following [19] and [20]. Both \(H_{d,h}^E\) and \(\phi_{d,h}^E\) are defined in 17 and 18 as integrals of certain functions that involve the exact trajectories of \(\mathcal{M}\). Thus, one can use a quadrature formula to approximate the integral and a numerical integrator for the (initial value problem) that approximates the exact trajectory of \(\mathcal{M}\). More precisely, fix a Gaussian quadrature formula such that, for \(f:[0,h]\rightarrow \mathbb{R}\), it is accurate to order \(a\), that is, \[\label{eq:shooting-quadrature95precission} \sum_{j=0}^n b_j f(c_jh) = \int_0^h f(x) dx +\mathcal{O}(h^{a+1}).\tag{29}\] We also fix a numerical integrator \(\Phi_h:T^*Q\rightarrow T^*Q\) for Hamilton’s equations ?? that is accurate to order \(b\), that is, \[\label{eq:shooting-integrator95precission} \Phi_h(q,p) = \mathbf{F}^{{\mathcal{M}}}_{{h}}(q,p) + \mathcal{O}(h^{b+1}).\tag{30}\]

Ideally, we would proceed as follows: in order to define \(H_{d,h}(q_0,p_1)\) we would find \(p_0\in Q^*\) such that \(\operatorname{pr}_2(\mathbf{F}^{{\mathcal{M}}}_{{h}}(q_0,p_0))=p_1\). Then, we would use \(\Phi_{c_j h}(q_0,p_0)\) to define a discrete trajectory (that approximates \(\mathbf{F}^{{\mathcal{M}}}_{{c_j h}}(q_0,p_0)\)) and, last, use the quadrature formula and 17 to define \(H_{d,h}(q_0,p_1)\). Unfortunately, we don’t know \(\mathbf{F}^{{\mathcal{M}}}_{{h}}(q_0,p_0)\), so that we replace \(p_0\) by \(\tilde{{p_0}}\in Q^*\) such that \(\operatorname{pr}_2(\Phi_h(q_0,\tilde{{p_0}})) = p_1\). It can be seen that \(\tilde{{p_0}} = p_0 +\mathcal{O}(h^{b+1})\) (this is analogous to Lemma 3.1 in [20]). Then, we define \[(q^j,p^j)\mathrel{\vcenter{:}}=\Phi_{c_j h}(q_0,\tilde{{p_0}}) \quad\text{ and }\quad v^j \mathrel{\vcenter{:}}= D_2H(q^j,p^j) \quad\text{ for }\quad j=0,\ldots,n.\] Notice that, for \((q(t),p(t)) \mathrel{\vcenter{:}}=\mathbf{F}^{{\mathcal{M}}}_{{t}}(q_0,p_0)\), \[\begin{align} (q(c_j h),p(c_j h)) =& \mathbf{F}^{{\mathcal{M}}}_{{t}}(q_0,p_0) = \mathbf{F}^{{\mathcal{M}}}_{{t}}(q_0,\tilde{{p_0}} + \mathcal{O}(h^{b+1})) = \mathbf{F}^{{\mathcal{M}}}_{{t}}(q_0,\tilde{{p_0}}) + \mathcal{O}(h^{b+1}) \\=& \Phi_{c_j h}(q_0,\tilde{{p_0}}) +\mathcal{O}(h^{b+1}) = (q^j,p^j) +\mathcal{O}(h^{b+1}), \end{align}\] where the third equality follows by a first order Taylor expansion of the flow and the fourth is due to the accuracy of the numerical integrator. On the other hand, as \((q(t),p(t))\) solves ?? , we have \[\begin{align} \dot{q}(c_j h) =& D_2H(q(c_j h),p(c_j h)) = D_2H((q^j,p^j)+\mathcal{O}(h^{b+1})) \\=& D_2H(q^j,p^j) + \mathcal{O}(h^{b+1}) = v^j +\mathcal{O}(h^{b+1}). \end{align}\] Last, we define \[\label{eq:shooting95H95phi-def}\begin{gather}[t] H_{d,h}(q_0,p_1) \mathrel{\vcenter{:}}= p^nq^n-\sum_{j=0}^n b_j(p^jv^j-H(q^j,p^j)),\\ \phi_{d,h}(q_0,p_1)(\delta q_0,p_1) \mathrel{\vcenter{:}}=\sum_{j=1}^n b_j\check{\phi}(q^j,p^j)(T_{(q_0,p_1)}q^j(\delta q_0,\delta p_1)). \end{gather} \tag{31}\] Notice that, by construction, all \((q^j,p^j)\) are functions of \(q_0\) and \(p_1\) and, as long as \(\Phi_h\) is \(C^1\), the derivative \(T_{(q_0,p_1)}q^j\) that appears in \(\phi_{d,h}\) is well defined. In addition, it follows from \(q^j = q(c_j h) +\mathcal{O}(h^{b+1})\) that \(T_{(q_0,p_1)}q^j = T_{(q_0,p_1)}q(c_j h) + \mathcal{O}(h^{b+1})\).

The following result, whose proof follows by unraveling the definitions, shows that the order \(1\) discretization considered at the beginning of the section has some interesting properties in connection with the commutativity of the discretization and the passage from the Hamiltonian to the Lagrangian formalisms.

Acknowledgments↩︎

This document is the result of research partially supported by grants from the Universidad Nacional de Cuyo [code 06/80020240100069UN], Universidad Nacional de La Plata [code SX007], and CONICET.

References↩︎