# Properties of Hamiltonian variational integrators

Properties of Hamiltonian variational integrators Abstract Discrete Hamiltonian variational integrators are derived from type II and type III generating functions for symplectic maps, and in this article, we establish a variational error analysis result that relates the order of accuracy of the associated numerical methods with the extent to which these generating functions approximate the exact discrete Hamiltonians. We also introduce the notion of an adjoint discrete Hamiltonian and relate it to the adjoint of the associated symplectic integrator. We show that when constructing discrete Lagrangians and discrete Hamiltonians using the same approximation method, this does not necessarily lead to the same symplectic integrator. Numerical experiments also indicate that the resonance behavior of variational integrators based on averaging methods depends on the type of generating functions used, and we relate this resonance behavior to the ill-posedness of the boundary-value problems used to define the exact discrete Lagrangian and exact discrete Hamiltonian. 1. Introduction Geometric numerical integration is a field of numerical analysis that develops numerical methods with the goal of preserving geometric properties of dynamical systems (see Hairer et al., 2006). Variational integrators are geometric numerical integrators derived from discretizing Hamilton’s principle from classical mechanics (see Marsden and West, 2001). They have many desirable properties such as symplecticity, momentum preservation and near-energy preservation, which results in excellent long-term stability. While the Lagrangian formulation of variational integrators has been thoroughly investigated (see Marsden et al., 1998; Marsden and West, 2001; Lew et al., 2003; de León et al., 2007; Leyendecker et al., 2008; Leok and Shingel, 2012), only recently has the Hamiltonian formulation of variational integrators been established (see Lall and West, 2006; Leok and Zhang, 2011). In this article, we will continue the investigation of Hamiltonian variational integrators and establish theorems on error analysis and symmetry of the method, and provide numerical experiments to elucidate the relative numerical advantages and disadvantages of the Lagrangian and Hamiltonian formulations. In particular, evidence is presented to show that for oscillatory problems the discrete Lagrangian and discrete Hamiltonian variational integrators have differing resonance and conditioning properties. In addition, it is shown that some approximation methods will yield a symmetric method only when derived from a specific type of generating function. The upshot is that the numerical properties of a variational integrator are determined both by the approximation scheme used to construct it and by the type of the generating function being approximated. 1.1 Discrete mechanics Lagrangian variational integrators are based on a discrete analog of Hamilton’s principle, whereas Hamiltonian variational integrators are based on a discrete analog of Hamilton’s phase space variational principle. The fundamental objects in the discretization are generating functions of symplectic maps, and in the Hamiltonian case, they are obtained by approximating the exact type II generating function associated with a Hamiltonian flow, which we refer to as the exact discrete right Hamiltonian,   Hd+,E(q0,p1)=ext(q,p)∈C2([0,T],T∗Q)q(0)=q0,p(T)=p1(p1q1−∫0T[pq˙−H(q,p)]dt). (1.1) This can be viewed as the solution at time $$T$$ of the type II Hamilton–Jacobi equation   ∂S2(q0,p,t)∂t=H(∂S2∂p,p), (1.2) which more generally describes the type II generating function that generates the time-$$t$$ Hamiltonian flow map,   S2(q0,p,t)=ext(q,p)∈C2([0,t],T∗Q)q(0)=q0,p(t)=p(p(t)q(t)−∫0t[p(s)q˙(s)−H(q(s),p(s))]ds). (1.3) Similarly, the exact discrete left Hamiltonian is given by   Hd−,E(p0,q1;h)=ext(q,p)∈C2([0,T],T∗Q)q(0)=q0,p(T)=p1(−p0q0−∫0T[pq˙−H(q,p)]dt) (1.4) and it can be viewed as a solution at time $$T$$ of the type III Hamilton–Jacobi equation   ∂S3(p0,q,t)∂t=H(q,−∂S3∂q). (1.5) Given discrete Hamiltonians $$H_d^+(q_k,p_{k+1};h)$$ and $$H_d^-(p_k,q_{k+1};h)$$, the discrete Hamilton’s equations are given by   qk+1 =D2Hd+(qk,pk+1;h), (1.6)  pk =D1Hd+(qk,pk+1;h) (1.7) and   qk =−D1Hd−(pk,qk+1;h), (1.8)  pk+1 =−D2Hd−(pk,qk+1;h). (1.9) These can also be expressed in terms of the discrete Legendre transformations $$\mathbb{F}^\pm H_d^+:(q_k,p_{k+1})\rightarrow T^*Q$$,   F+Hd+(qk,pk+1;h) =(D2Hd+(qk,pk+1;h),pk+1), (1.10)  F−Hd+(qk,pk+1;h) =(qk,D1Hd+(qk,pk+1;h)) (1.11) and $$\mathbb{F}^\pm H_d^-:(p_k,q_{k+1})\rightarrow T^*Q$$,   F+Hd−(pk,qk+1;h) =(qk+1,−D2Hd−(pk,qk+1;h)), (1.12)  F−Hd−(pk,qk+1;h) =(−D1Hd−(pk,qk+1;h),pk). (1.13) We observe that the Hamiltonian maps $$\tilde{F}_{H_d^\pm}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$ can be expressed as   F~Hd±=F+Hd±∘(F−Hd±)−1. (1.14) 2. Error analysis and symmetric methods 2.1 Error analysis Variational integrators are able to benefit from and adopt many traditional techniques and methods of numerical analysis (see Leok and Shingel, 2012). This can be largely attributed to the following theorem from Marsden and West (2001). Theorem 2.1 (Marsden and West, 2001, Theorem 2.3.1). If a discrete Lagrangian, $$L_{\rm d}:Q \times Q \rightarrow \mathbb{R}$$ approximates the exact discrete Lagrangian $$L_{\rm d}^{\rm E}:Q \times Q \rightarrow \mathbb{R}$$ to order $$r$$, i.e.,   Ld(q0,q1;h)=LdE(q0,q1;h)+O(hr+1), then the discrete Hamiltonian map $$\tilde{F}_{L_d}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$, viewed as a one-step method, is order-$$r$$ accurate. Thus, in order to generate a variational integrator of a particular order, one can leverage techniques from numerical analysis with the goal of approximating the exact discrete Lagrangian, then the associated discrete Hamiltonian map yields the variational integrator. We first present the corresponding theorem for discrete Hamiltonian variational integrators, which draws much of its inspiration from the theorem and proof of the above result as detailed in Marsden and West (2001). Theorem 2.2 If a discrete right Hamiltonian $$H^+_d:T^*Q \rightarrow \mathbb{R}$$ approximates the exact discrete right Hamiltonian $$H_{\rm d}^{+,{\rm E}}:T^*Q \rightarrow \mathbb{R}$$ to order $$r$$, i.e.,   Hd+(q0,p1;h)=Hd+,E(q0,p1;h)+O(hr+1), and the Hamiltonian is continuously differentiable then the discrete map $$\tilde{F}^h_{H_d^+}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$, viewed as a one-step method, is order-$$r$$ accurate. Note that the following proof can be easily adjusted to prove an equivalent theorem for the discrete left Hamiltonian case. First, we will need the following lemma. Lemma 2.3 Let $$f_1,g_1,e_1,f_2,g_2,e_2 \in C^r$$ be such that   f1(x,h) =g1(x,h)+hr+1e1(x,h),f2(x,h) =g2(x,h)+hr+1e2(x,h). Then, there exist functions $$e_{12}$$ and $$\bar{e}_1$$ bounded on compact sets such that   f2(f1(x,h),h) =g2(g1(x,h),h)+hr+1e12(g1(x,h),h),f1−1(y) =g1−1(y)+hr+1e¯1(y). Proof. We have   f2(f1(x,h),h) =f2(g1(x,h)+hr+1e1(x,h),h) =g2(g1(x,h)+hr+1e1(x,h),h)+hr+1e2(g1(x,h)+hr+1e1(x,h),h) =g2(g1(x,h),h)+hr+1e~1(g1(x,h),h)+hr+1e2(g1(x,h)+hr+1e1(x,h),h), where $$\tilde{e}_1$$ is bounded on a compact set. This last line comes from combining the compactness of the set with the smoothness of the functions to obtain a Lipschitz property of the form   ‖g2(g1(x,h)+hr+1e1(x,h),h)−g2(g1(x,h),h)‖≤Chr+1. For each choice of $$(x,h)$$, equality holds for a particular choice of constant, which defines $$\tilde{e}_1$$ and establishes its smoothness as a function. Adding $$e_2$$ to $$\tilde{e}_1$$ we obtain a function $$e_{12}$$, which is also bounded on compact sets such that   f2(f1(x,h),h)=g2(g1(x,h),h)+hr+1e12(g1(x,h),h). Let $$y=f_1(x,h)$$ and note that by definition,   f1−1(f1(x,h))=g1−1(g1(x,h)). Since $$g_1^{-1}(y)=g_1^{-1}(g_1(x,h)+h^{r+1}e_1(x,h))$$, then   ‖g1−1(y)−f1−1(y)‖=‖g1−1(y)−g1−1(g1(x,h))‖≤C¯hr+1. From this, it follows that there exists a function $$\bar{e}_1$$ bounded on compact sets such that   f1−1(y)=g1−1(y)+hr+1e¯1(y). □ Now we are ready for the proof of the theorem. Proof. By assumption there is some bounded continuously differentiable function $$e$$ such that   Hd+(q(0),p(h);h)=Hd+,E(q(0),p(h),h)+hr+1e(q(0),p(h),h). Differentiating yields   D1Hd+(q(0),p(h);h)=D1Hd+,E(q(0),p(h);h)+hr+1D1e(q(0),p(h),h), where $$\|D_1e(q(0),p(h),h)\| \leq \tilde{C}$$. This implies   ‖F−Hd+(q(0),p(h);h)−F−Hd+,E(q(0),p(h);h)‖≤C~hr+1. Now combining this with the fact that $$\tilde{F}_{H^+_d}=\mathbb{F}^+H^+_d \circ (\mathbb{F}^-H^+_d)^{-1}$$ and applying Lemma 2.3, we have   F~Hd+h=F~Hd+,Eh+O(hr+1). □ Determining the order of a variational integrator is greatly simplified via the above theorems, which relate the order of the integrator to the order to which the associated discrete Lagrangian or discrete right Hamiltonian approximates the corresponding exact generating function. Similarly, it was shown in Marsden and West (2001) that one can determine whether the variational integrator is a symmetric method by examining the corresponding discrete Lagrangian. We would like to extend this result to the case of discrete Hamiltonians. 2.2 Symmetric methods Definition 2.4 (See Hairer et al., 2006, Chapters II.3 and V). A numerical one-step method $${\it{\Phi}}_h$$ is called symmetric or time reversible if it satisfies   Φh∘Φ−h=id or equivalently   Φh=Φ−h−1. The adjoint of a numerical one-step method, denoted $${\it{\Phi}}_h^*$$, is defined as   Φh∗=Φ−h−1. A numerical one-step method is a symmetric method if it is self-adjoint, i.e., $${\it{\Phi}}_h={\it{\Phi}}_h^*$$. The adjoint of a discrete Lagrangian $$L_d^*$$ is defined as   Ld∗(q0,q1,h)=−Ld(q1,q0,−h). The discrete Lagrangian is called self-adjoint if $$L_d^*(q_0,q_1,h)=L_d(q_0,q_1,h)$$. The following theorem from Marsden and West (2001) relates the self-adjointness of the discrete Lagrangian with the self-adjointness of the corresponding variational integrator. Theorem 2.5 (Marsden and West, 2001, Theorem 2.4.1). The discrete Lagrangian (or an equivalent discrete Lagrangian) $$L_d$$ is self-adjoint if and only if the method associated with the corresponding discrete Hamiltonian map is self-adjoint, i.e., symmetric. In many cases, it is easier to check whether the discrete Lagrangian is self-adjoint rather than checking the variational integrator itself. We seek a definition for the adjoint of a discrete right Hamiltonian. The adjoint of a one-step method $$(q_1,p_1)={\it{\Phi}}_h(q_0,p_0)$$ can be obtained by reversing the direction of time and reversing the roles of the initial data and terminal solution, i.e., $$(q_0,p_0)={\it{\Phi}}^*_{-h}(q_1,p_1)$$. This corresponds to swapping out $$(q_0,p_0,q_1,p_1,h)$$ for $$(q_1,p_1,q_0,p_0,-h)$$. This motivates the definition of the adjoint of a type II generating function as a type III generating function and vice versa. In particular, given a type II discrete Hamiltonian $$H_d^+$$, we seek a definition for the type III adjoint $$(H_d^+)^*$$ that will satisfy $$\tilde{F}^h_{(H_d^+)^*} = (\tilde{F}^h_{H_d^+})^*$$. Let $$\tilde{F}^h_{(H_d^+)^*}(q_0,p_0)=(q_1,p_1)$$. Then, we want   (q1,p1) =F~(Hd+)∗h(q0,p0) =(F~Hd+h)∗(q0,p0) =(F~Hd+−h)−1(q0,p0). This implies $$\tilde{F}^{-h}_{H_d^+}(q_1,p_1) = (q_0,p_0)$$, which, together with $$\tilde{F}^h_{(H_d^+)^*}(q_0,p_0)=(q_1,p_1),$$ yields the respective sets of equations   p1 =D1Hd+(q1,p0;−h),q0 =D2Hd+(q1,p0;−h) and   p1 =−D2(Hd+)∗(p0,q1;h),q0 =−D1(Hd+)∗(p0,q1;h). Comparing these equations we see that setting $$(H_d^+)^*(p_0,q_1;h) = -H_d^+(q_1,p_0;-h)$$ satisfies $$\tilde{F}^h_{(H_d^+)^*} = (\tilde{F}^h_{H_d^+})^*$$. A similar calculation yields an analogous expression for the adjoint of a type III generating function $$H_d^-$$. Definition 2.6 Given a type II/III generating function, $$H_d^\pm$$, define the adjoint as the type III/II generating function, $$(H_d^\pm)^*$$, where $$\tilde{F}_{(H_d^\pm)^*}^h(q_0,p_0) = (q_1,p_1)$$, as   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h), (2.1)  (Hd−)∗(q0,p1;h) =−Hd−(p1,q0;−h). (2.2) Example 2.7 The symplectic Euler-A method for a Lagrangian of the form $$L(q,\dot{q}) = \frac{1}{2}\dot{q}M\dot{q} - V(q)$$ is given by   p1 =p0−h∇V(q0),q1 =q0+hM−1p1. The corresponding discrete right Hamiltonian is given by   Hd+(q0,p1,h) =p1(q0+hM−1p1)−h[p1M−1p1−H(q0,p1)] =p1q0+hH(q0,p1). The adjoint of this method is given by symplectic Euler-B   q1=q0+hM−1p0,p1=p0−h∇V(q1). We now derive the corresponding adjoint of the discrete right Hamiltonian for symplectic Euler-A:   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h) =−p0(q1−hM−1p0)−h[p0M−1p0−H(q1,p0)] =−p0q1+hH(q1,p0). We can verify that this generates symplectic Euler-B by applying the discrete left Hamilton equations   q0 =−D1(Hd+)∗(p0,q1;h) =D2Hd+(q1,p0;−h) =q1−hM−1p0,p1 =−D2(Hd+)∗(p0,q1;h) =D1Hd+(q1,p0;−h) =p0−h∇V(q1). Solving the first equation for $$q_1$$ gives symplectic Euler-B, as expected. Theorem 2.8 $$(H_d^\pm)^{**}=H_d^\pm$$. Proof. We consider the case of the type II generating function $$H_d^+$$. Let $$\tilde{F}^h_{(H_d^+)^{**}}(q_0,p_0)=(q_1,p_1)$$. Since $$(H_d^+)^*$$ is a type III generating function, applying the definition of the adjoint twice gives   (Hd+)∗∗(q0,p1;h) =−(Hd+)∗(p1,q0;−h) =Hd+(q0,p1;h) and a similar calculation shows that this holds for the type III generating function $$H_d^-$$ as well. □ Since the notion of the adjoint that we introduced converts a type II to a type III generating function, for a discrete Hamiltonian to be self-adjoint, we need to compare the adjoint to the Legendre transformation of the discrete Hamiltonian, which is given by   Hd−(pk,qk+1;h)=−pkqk−pk+1qk+1+Hd+(qk,pk+1;h), where we view $$p_{k+1}$$ and $$q_k$$ as functions of $$p_k$$ and $$q_{k+1}$$. Then, the following calculation shows that these two generating functions generate the same symplectic map, i.e., $$\tilde{F}_{H_d^-}=\tilde{F}_{H_d^+}$$,   −D1Hd−(pk,qk+1;h) =qk+pk∂qk∂pk+∂pk+1∂pkqk+1−D1Hd+(qk,pk+1;h)∂qk∂pk−D2Hd+(qk,pk+1;h)∂pk+1∂pk =qk+(pk−D1Hd+(qk,pk+1;h))∂qk∂pk+(qk+1−D2Hd+(qk,pk+1;h))∂pk+1∂pk,−D2Hd−(pk,qk+1;h) =pk∂qk∂qk+1+∂pk+1∂qk+1qk+1+pk+1−D1Hd+(qk,pk+1;h)∂qk∂qk+1 −D2Hd+(qk,pk+1;h)∂pk+1∂qk+1 =pk+1+(pk−D1Hd+(qk,pk+1;h))∂qk∂qk+1+(qk+1−D2Hd+(qk,pk+1;h))∂pk+1∂qk+1. Definition 2.9 A type II/III generating function is self-adjoint if it is equal (up to equivalency) to the Legendre transform of its adjoint. Note that this definition implies that a discrete right Hamiltonian is self-adjoint if its adjoint is equal (up to equivalency) to the associated discrete left Hamiltonian, i.e., $$(H_d^+)^* = H_d^-$$. Theorem 2.10 Given a self-adjoint discrete right Hamiltonian, i.e., $$H_d^-=(H_d^+)^*$$, the method associated with the discrete right Hamiltonian map is self-adjoint. Likewise, if a method coming from a discrete right Hamiltonian map is self-adjoint, then the associated discrete right Hamiltonian is self-adjoint. Proof. Assume $$H_d^-=(H_d^+)^*$$. Then,   (F~Hd+)∗=F~(Hd+)∗=F~Hd−=F~Hd+, and so, by definition, the map is self-adjoint. Now assume $$\tilde{F}_{H_d^+}=(\tilde{F}_{H_d^+})^*$$. Then,   F~Hd−=F~Hd+=(F~Hd+)∗=F~(Hd+)∗, which implies $$(H_d^+)^*=H_d^-$$ (up to equivalency) and, by definition, the discrete right Hamiltonian is self-adjoint. □ A similar proof can be used to prove an identical theorem for the discrete left Hamiltonian. The previous theorem provides an easy way to check whether a variational integrator is self-adjoint. Assuming the Hamiltonian flow is time reversible, it follows that the exact discrete right Hamiltonian is self-adjoint. This can also be shown using the definition of a self-adjoint exact discrete right Hamiltonian, and the same result can be shown for the discrete left Hamiltonian with only minor adjustments. Corollary 2.11 The exact discrete right Hamiltonian $$H_{\rm d}^{+,{\rm E}}$$ is self-adjoint. Proof. A direct calculation shows that   (Hd+,E)∗(p0,q1;h) =−Hd+,E(q1,p0;−h) =−(p~(−h)q~(−h)−∫0−h[p~(τ)q~(τ)−H(q~(τ),p~(τ))]dτ) =−p(−h+h)q(−h+h)−∫−h0[p(τ+h)q(τ+h)−H(q(τ+h),p(τ+h))]dτ =−p(0)q(0)−∫0h[p(t)q(t)−H(q(t),p(t))]dt =Hd−,E(p0,q1;h), where we used the fact that the time-reversed solution $$(\tilde{q}(\tau),\tilde{p}(\tau))$$ over the time domain $$[-h,0]$$ with $$(q_1,p_0)$$ boundary data is related to the solution curve $$(q(t),p(t))$$ over the time domain $$[0,h]$$ with $$(q_0,p_1)$$ boundary data by $$(\tilde{q}(\tau),\tilde{p}(\tau))=(q(\tau+h),p(\tau+h))$$. □ The definition of the adjoint also provides a simple way to construct symmetric methods. Given any method defined by $$H_d$$, we can construct a symmetric method using composition, for example, $$\tilde{F}^{{h}/{2}}_{H_d} \circ \tilde{F}^{{h}/{2}}_{H_d^*}$$, which is nothing more than composing a half-step of the adjoint method with a half-step of the method. It is well known that this leads to a symmetric method, as the following calculation demonstrates,   (F~Hdh/2∘F~Hd∗h/2)∗ =(F~Hd∗h/2)∗∘(F~Hdh/2)∗ =F~Hd∗∗h/2∘F~Hd∗h/2 =F~Hdh/2∘F~Hd∗h/2, where the last line used Theorem 2.8. More generally, a composition method of the form,   F~Hdαsh∘F~Hd∗βsh∘⋯∘F~Hd∗β2h∘F~Hdα1h∘F~Hd∗β1h, where $$\alpha_{s+1-i}=\beta_i$$ for $$i=1,\ldots,s$$, will be symmetric. For a more in-depth discussion of symmetric composition methods, see Hairer et al. (2006, Chapter V.3). 3. Discrete Lagrangians vs. discrete Hamiltonians 3.1 Composition of discretization and the Legendre transform A transformation of one type of generating function into another type of generating function, which preserves the associated sympletic map, is given by the following Legendre transforms. Definition 3.1 (i) The Legendre transform of a type I generating function into a type II generating function is given by the equation   Hd+(qk,pk+1;h)=pk+1qk+1−Ld(qk,qk+1;h), where $$q_{k+1}$$ is implicitly defined by $$p_{k+1}=D_2L_d(q_k,q_{k+1};h)$$. The Legendre transform of a type II generating function into type I is given by the same equation, where $$p_{k+1}$$ is implicitly defined by $$q_{k+1}=D_2H^+_d(q_k,p_{k+1};h)$$. (ii) The Legendre transform of a type I generating function into a type III generating function is given by the equation   Hd−(pk,qk+1;h)=−pkqk−Ld(qk,qk+1;h), where $$q_k$$ is implicitly defined by $$p_k=-D_1L_d(q_k,q_{k+1};h)$$. The Legendre transform of a type III generating function into type I is given by the same equation, where $$p_k$$ is implicitly defined by $$q_k=-D_1H^-_d(p_k,q_{k+1};h)$$. (iii) The Legendre transform of a type II generating function into a type III generating function is given by the equation   Hd−(pk,qk+1;h)=−pkqk−pk+1qk+1+Hd+(qk,pk+1;h), where $$q_k$$ and $$p_{k+1}$$ are implicitly defined by the set of equations $$p_k=D_1H^+_d(q_k,p_{k+1};h)$$ and $$q_{k+1}=D_2H^+_d(q_k,p_{k+1};h)$$. The Legendre transform of a type III generating function into type II is given by the same equation, where $$q_{k+1}$$ and $$p_k$$ are implicitly defined by the set of equations $$q_k=-H^-_d(p_k,q_{k+1};h)$$ and $$p_{k+1}=-H^-_d(p_k,q_{k+1})$$. Variational integrators are derived by first discretizing the exact generating function, which results in a new generating function of the same type. Using the Legendre transforms defined above, the discretized generating function can be transformed into an equivalent generating function of a different type. This can be viewed as composing the Legendre transform and the discretization. Likewise, one could first take the Legendre transform of the given exact generating function and then apply the discretization. The question we address next is whether changing the order of this composition results in the same symplectic map, and for one particular type of discretization this question has already been answered. It was shown in Leok and Zhang (2011) that the Galerkin variational integrator construction leads to equivalent discrete Lagrangian and discrete Hamiltonian methods for the same choice of quadrature rule and finite-dimensional function space, and the result is given in the following theorem. Theorem 3.2 (Leok and Zhang, 2011, Proposition 4.1). If the continuous Hamiltonian $$H(q,p)$$ is hyperregular and we construct a Lagrangian $$L(q,\dot{q})$$ by the Legendre transformation then the generalized Galerkin Hamiltonian variational integrator (see Leok and Zhang, 2011) and the generalized Galerkin Lagrangian variational integrator, associated with the same choice of basis functions and numerical quadrature formula, are equivalent. Will this hold for other types of variational integrators? To begin to address this question, we propose the following examples. Example 3.3 Consider a Lagrangian of the form $$L(q,\dot{q})=\frac{1}{2}\dot{q}M\dot{q} - V(q)$$, where $$M$$ is symmetric positive definite and $$V$$ is sufficiently smooth. The exact discrete Lagrangian, which is defined as   LdE(q0,q1;h) =∫0hL(q01(t),q˙01(t))dt, where $$q_{01}(0)=q_0,$$$$q_{01}(h)=q_1,$$ and $$q_{01}$$ satisfies the Euler–Lagrange equation in the time interval $$(0,h)$$. Letting $$q_0$$ and $$q_1$$ be fixed, then a first-order finite difference approximation of the velocities yields   q˙0 ≈q1−q0h,q˙1 ≈q1−q0h. Using the rectangular quadrature rule about the initial point results in the following discrete Lagrangian (i.e., type I generating function):   Ld(q0,q1;h)=12(q1−q0h)M(q1−q0h)−V(q0). (3.1) Applying the implicit discrete Euler–Lagrange equations yields   p0=Mq1−q0h+h∇V(q0),p1=Mq1−q0h. Finally rearranging these equations results in the variational integrator, also known as symplectic Euler-A,   q1=q0+hM−1p1,p1=p0−h∇V(q0). Now we apply this same approximation scheme to the associated exact type II generating function. The boundary-value formulation of the exact discrete right Hamiltonian is given by   Hd+,E(q0,p1) =(p1q1−∫0h[pq˙−H(q,p)]dt), where $$(q(t),p(t))$$ satisfy Hamilton’s equations with boundary conditions $$q(0)=q_0$$, $$p(h)=p_1$$. The first-order finite difference approximations of the velocities yield the following equations for the momentum:   p0 ≈Mq1−q0h,p1 ≈Mq1−q0h. Solving for $$q_1$$ in terms of $$q_0$$ and $$p_1$$ yields the approximation $$q_1\approx q_0 + hM^{-1}p_1$$, and this simplifies the approximation to $$p_0$$ as $$p_0 \approx p_1$$. Now using these approximations in combination with the rectangular rule about the initial point yields the discrete right Hamiltonian   Hd+(q0,p1;h)=p1(q0+hM−1p1)−h(12p1M−1p1−V(q0)). After applying the implicit discrete Hamilton equations and rearranging terms, the resulting method is again symplectic Euler-A:   q1=q0+hM−1p1p1=p0−h∇V(q0). The composition of the Legendre transform and the discretization of a generating function will be called commutative if, regardless of the order of this composition, either of the resulting generating functions leads to the same symplectic map. In this context, the previous example shows that this particular discretization scheme, which includes both the trajectory approximation and the quadrature rule, commutes with the Legendre transform between type I and type II generating functions. Now we look at the same velocity approximation but using the rectangular rule about the end point rather than the initial point. Example 3.4 As before, we first build the discrete Lagrangian with $$q_0$$, $$q_1$$ fixed and velocity approximations   q˙0 ≈q1−q0h,q˙1 ≈q1−q0h. Applying the rectangular rule about the end point yields the discrete Lagrangian   Ld(q0,q1;h)=12(q1−q0h)M(q1−q0h)−V(q1) (3.2) and the associated variational integrator is symplectic Euler-B:   q1=q0+hM−1p0,p1=p0−h∇V(q1). Now applying the same approximation scheme to construct the discrete right Hamiltonian results in the following type II generating function:   Hd+(q0,p1;h)=p1(q0+hM−1p1)−h(12p1M−1p1−V(q0+hM−1p1)). Applying the implicit discrete Hamilton equations and rearranging the terms results in the method   q1=q0+hM−1p0,p1=p0−h∇V(q0+hM−1p1), which is not symplectic Euler-B. We have just proven the following theorem. Theorem 3.5 In general, the composition of the Legendre transform of a generating function and the discretization of a generating function do not commute. Therefore, the answer to our original question is that in general, a fixed approximation scheme used to construct a discrete Lagrangian will not generate the same method when it is used to construct a discrete Hamiltonian. In general, how might the two resulting methods differ? A complete characterization of this issue is subtle, and beyond the scope of this article, but it will be a topic of future work. For now, we will consider how the two approaches differ when combined with the method of averaging, which will also serve to illustrate how the type of generating function and the associated boundary data can affect the numerical properties of the method. 3.2 Averaged Lagrangians and Hamiltonians Averaging methods have played a role in solving differential equations since at least as far back as the time of Lagrange (see Verhulst, 2000), and they continue to play a key role particularly in the field of numerical differential equations applied to nearly integrable systems or problems with multiple time scales. We consider perturbed Hamiltonian systems with Hamiltonians of the form   H=H(A)+ϵH(B), (3.3) where $$\epsilon \ll 1$$ and the dynamics of the Hamiltonian system corresponding to $$H^{(A)}$$ is exactly solvable or at the very least cheap to approximate. We call this an almost-integrable system, the motivation being that the dynamics of the system are largely influenced by an integrable Hamiltonian with simpler dynamics, but smaller influences also play a role in the overall dynamics. An example is the classic $$n$$-body problem of the solar system, where a particular planet’s trajectory is largely influenced by the sun, but other planets and nearby objects also play a role. Averaging methods can be constructed to exploit the larger influence of $$H^{(A)}$$ on the dynamics of the system by averaging out the smaller influences. Ideally, averaging techniques will allow for larger time steps to be used while still yielding a reasonable approximation to the solution. A variational integrator for such a system was proposed in Farr (2009) using a discrete Lagrangian formulation, which drew inspiration from the kick-drift-kick leapfrog method (see Wisdom and Holman, 1991). We will discuss the original Lagrangian formulation (hereafter referred to as the averaged Lagrangian) and in addition construct an analogous method in terms of a discrete right Hamiltonian (referred to as the averaged Hamiltonian). The Lagrangian corresponding to (3.3) is given by   L=L(A)+ϵL(B) (3.4) and we will make the assumption that $$L^{(B)}(q(t),\dot{q}(t))=-V^{(B)}(q(t))$$. The averaging method of interest has a local truncation error of $$\mathcal{O}(\epsilon^2 h^3)$$ and is defined in terms of a discrete Lagrangian, $$L_d$$. This method, proposed in Farr (2009), uses a discrete Lagrangian of the form   Ld(q0,q1,h) =Ld(A),E(q0,q1;h)+ϵ∫0hL(B)(qA(q0,q1,t),q˙A(q0,q1,t))dt =Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt, where we denote the trajectory corresponding to $$L^{(A)}$$ with boundary conditions $$(q_0,q_1)$$ by $$(q_A(t),\dot{q}_A(t))$$. The idea is to use the dynamics of $$L^{(A)}$$, which is solved for either exactly, or efficiently and accurately approximated, to average the contribution of $$L^{(B)}$$ to the dynamics. The corresponding discrete Hamiltonian map is given implicitly by   −p0 =D1Ld(A),E(q0,q1;h)−ϵ∫0hD1V(B)(qA(q0,q1,t))dt, (3.5a)  p1 =D2Ld(A),E(q0,q1;h)−ϵ∫0hD2V(B)(qA(q0,q1,t))dt. (3.5b) As shown in Farr (2009), this method has local truncation error of size $$\mathcal{O}(\epsilon^2 h^3)$$. Using the notation $$p_0^A(q_0,q_1)=-D_1L_d^{(A),E}(q_0,q_1;h)$$ and $$p_1^A(q_0,q_1)=D_2L_d^{(A),E}(q_0,q_1;h)$$, we rearrange the above equations to get   p0−ϵ∫0hD1V(B)(qA(q0,q1,t))dt=p0A(q0,q1), (3.6a)  p1=p1A(q0,q1)−ϵ∫0hD2V(B)(qA(q0,q1,t))dt. (3.6b) In Farr (2009), it is noted that $$- \epsilon \int_0^h D_1V^{(B)}(q_A(q_0,q_1,t)) \,{\rm d}t$$ is an average along the trajectory generated by $$L^{(A)}$$ which, in general, gives more weight to the initial periods of the trajectory, while $$- \epsilon \int_0^h D_2V^{(B)}(q_A(q_0,q_1,t)) \,{\rm d}t$$ is an average along the trajectory generated by $$L^{(A)}$$ that, in general, favors the latter periods of the trajectory. Now let us consider the discrete right Hamiltonian given by the same form of approximation,   Hd+(q0,p1;h)=Hd(A),+,E(q0,p1;h)+ϵ∫0hV(B)(qA(q0,p1,t))dt. The discrete right Hamiltonian map is given implicitly by   p0 =D1Hd(A),+,E(q0,p1;h)+ϵ∫0hD1V(B)(qA(q0,p1,t))dt,q1 =D2Hd(A),+,E(q0,p1;h)+ϵ∫0hD2V(B)(qA(q0,p1,t))dt. Using the notation $$p_0^A(q_0,p_1) = D_1 H_d^{(A),+,E}(q_0,p_1;h)$$ and $$q_1^A(q_0,p_1) = D_2 H_d^{(A),+,E}(q_0,p_1;h)$$, we rearrange the equations to yield   p0−ϵ∫0hD1V(B)(qA(q0,p1,t))dt=p0A(q0,p1), (3.7a)  q1=q1A(q0,p1)+ϵ∫0hD2V(B)(qA(q0,p1,t))dt. (3.7b) Theorem 3.6 The method defined implicitly by (3.7) has local truncation error $$\mathcal{O}(\epsilon^2 h^3)$$. Proof. Using variational error analysis, we need to show   O(ϵ2h3) =HdE,+−Hd+ =ΔA+ϵΔB, where $${\it{\Delta}}_A$$ is given by   p(h)q(h)−∫0h[p(t)q˙(t)−H(A)(q(t),p(t))]dt−(pA(h)qA(h)−∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt), and $$\epsilon{\it{\Delta}}_B$$ is given by   ϵ∫0h[V(B)(q(t))−V(B)(qA(t))]dt. Using a functional Taylor expansion, $${\it{\Delta}}_A$$ becomes   ΔA =δδqA(∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt)δqA +δ2δqA2(∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt)δqA2+O(δqA3), where $$\delta q_A$$ is the difference between $$q$$ and $$q_A$$. Noting that $$q$$ and $$q_A$$ differ in forces of order $$\epsilon h^2$$ and $$p$$ differs from $$p_A$$ to first order in $$\epsilon h$$, implies that $$\delta q_A$$ is on the order of $$\mathcal{O}(\epsilon h)$$. This can be seen explicitly by comparing Taylor expansions about time zero. Since $$q_A$$ satisfies Hamilton’s equations for $$H^{(A)}$$, the first variation vanishes (see Leok and Zhang, 2011, Lemma 2.1) leaving a term on the order of $$h \delta q_A^2$$. Therefore, we have   ΔA=O(ϵ2h3). Likewise, a functional Taylor expansion for $${\it{\Delta}}_B$$ yields,   ΔB=δδqA[∫0hV(B)(qA(t))dt]δqA+O(δqA2). Noting that $$V^{(B)}$$ is a function of only $$q_A$$ and that $$q$$ differs from $$q_A$$ on the order of $$\epsilon h^2$$, implies $$\epsilon {\it{\Delta}}_B = \mathcal{O}(\epsilon^2 h^3)$$. □ Theorem 3.7 Assuming the flow associated with $$L^{(A)}$$ is time reversible, then both methods, defined respectively by (3.6) and (3.7), are symmetric methods. Proof. The discrete Lagrangian associated with (3.6) is given by   Ld(q0,q1;h)=Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt. The adjoint of the discrete Lagrangian is given by   (Ld(q0,q1;h))∗ =−Ld(q1,q0;−h) =−Ld(A),E(q1,q0;−h)+ϵ∫0−hV(B)(qA(q1,q0,t))dt =−Ld(A),E(q1,q0;−h)−ϵ∫0hV(B)(qA(q1,q0,t))dt =Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt =Ld(q0,q1;h). The third equality comes from the time reversibility of the flow associated with $$L^{(A)}$$, and the fourth equality uses this property together with the fact that the exact discrete Lagrangian is self-adjoint. The discrete right Hamiltonian associated with (3.7) is given by   Hd+(q0,p1;h)=Hd(A),+,E(q0,p1;h)+ϵ∫0hV(B)(qA(q0,p1,t))dt. The adjoint of the discrete right Hamiltonian is given by   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h) =−Hd(A),+,E(q1,p0;−h)−ϵ∫0−hV(B)(qA(q1,p0,t))dt =−Hd(A),+,E(q1,p0;−h)+ϵ∫0hV(B)(qA(q1,p0,t))dt =Hd(A),−,E(p0,q1;h)+ϵ∫0hV(B)(qA(p0,q1,t))dt =Hd−(p0,q1;h), where the third equality comes from the time reversibility of the flow associated with $$H^{(A)}$$, and the fourth equality uses this property together with the fact that the exact discrete Hamiltonian is self-adjoint. By Definition 2.9 and Theorem 2.10, the method is symmetric. □ It can be shown that methods (3.6) and (3.7) do not result in the same symplectic map. How do these respective maps differ? To gain insight into this question, we now turn to numerical experimentation. 4. Numerical results 4.1 Exact generating functions First, we consider the unperturbed harmonic oscillator boundary-value problem   q¨(t)+q(t)=0,q(0)=q0, q(h)=q1. (4.1) Analytically, the boundary-value problem is not well posed when $$h$$ is an integer multiple of $$\pi$$, and in particular there are infinitely many solutions. Recall that the exact discrete Lagrangian is given by   LdE(q0,q1;h)=∫0hL(q01(t),q˙01(t))dt, (4.2) where $$q_{01}(0)=q_0$$, $$q_{01}(h)=q_1$$, and $$q_{01}(t)$$ satisfies the Euler–Lagrange equation in the time interval $$(0,h)$$. Thus, the exact type I generating function is ultimately defined in terms of such a boundary-value problem. The integrator obtained from the exact discrete Lagrangian is given by   q1 =q0cos⁡(h)+p0sin⁡(h),p1 =q1cot⁡(h)−q0csc⁡(h). This integrator is analytically the true solution of the harmonic oscillator initial-value problem, where $$(q(0),p(0))=(q_0,p_0)$$, and consequently the local truncation error is zero. However, noting that $$\cot(h)$$ and $$\csc(h)$$ both involve dividing by $$\sin(h)$$, we expect increased round-off error for values of $$h$$ that are near integer multiples of $$\pi$$. Similarly, the exact discrete right Hamiltonian is given by   Hd+,E(q0,p1;h)=p1q1−∫0h[p01(t)q˙01(t)−H(q01(t),p01(t))]dt, (4.3) where $$q_{01}(0)=q_0$$, $$p_{01}(h)=p_1$$, and $$(q_{01}(t),p_{01}(t))$$ satisfy Hamilton’s equations in the time interval $$(0,h)$$. This is related to the unperturbed harmonic oscillator boundary-value problem given by,   q˙(t)=p(t),p˙(t)=−q(t),q(0)=q0, p(h)=p1. (4.4) This boundary-value problem is not well posed for values of $$h$$ that are odd multiples of $$\frac{\pi}{2}$$ and there are infinitely many solutions for such values of $$h$$. The integrator obtained from the exact discrete right Hamiltonian for the harmonic oscillator is given by   p1 =p0cos⁡(h)−q0sin⁡(h),q1 =p1tan⁡(h)+q0sec⁡(h). This integrator is analytically the true solution to the harmonic oscillator initial-value problem, where $$(q(0),p(0))=(q_0,p_0)$$ and the local truncation error will be zero. Noting that the method involves $$\tan(h)$$ and $$\sec(h)$$, we expect increased round-off error around odd multiples of $$\frac{\pi}{2}$$. Both of the integrators given by the exact discrete Lagrangian and the exact discrete right Hamiltonian have been implemented for the harmonic oscillator with initial conditions $$(q_0,p_0) = (1,0)$$ over the time interval $$[0,10000]$$, and the energy error is shown in Fig. 1. Note the jump in round-off error corresponding to values of $$h$$ that are odd multiples of $$\pi$$ (for the discrete Lagrangian) and odd multiples of $$\frac{\pi}{2}$$ (for the discrete right Hamiltonian). The bottom plot takes the minimum error of the two methods, and this indicates that a step size causing noticeable round-off error for one method will work just fine for the other method. Fig. 1. View largeDownload slide The first plot is the energy error vs. step size for the exact discrete right Hamiltonian applied to the harmonic oscillator. The second plot shows the energy error vs. step size for the exact discrete Lagrangian, while the third plot takes the minimum of the energy error from either method. Fig. 1. View largeDownload slide The first plot is the energy error vs. step size for the exact discrete right Hamiltonian applied to the harmonic oscillator. The second plot shows the energy error vs. step size for the exact discrete Lagrangian, while the third plot takes the minimum of the energy error from either method. In this particular case, we can conclude that the numerical difference between the symplectic maps generated by the respective exact discrete Lagrangian and exact discrete right Hamiltonian is a matter of numerical conditioning, which is inherited from the underlying ill-posedness of the associated boundary-value problem. Despite the fact that the methods are applied to an initial-value problem, numerical properties can be attributed to a boundary-value problem that is no longer visible in the methods themselves. Considering that many symplectic integrators are derived independently of the variational integrator formulation, perhaps some of their numerical properties can be better understood by reinterpreting them in the framework of variational integrators. 4.2 Averaged variational integrators for a nonlinearly perturbed harmonic oscillator Now we consider the previous averaging methods applied to a Hamiltonian of the form   H(q,p)=12(p2+q2)+ϵ3q3, (4.5) which is the Hamiltonian for a nonlinearly perturbed harmonic oscillator. The corresponding averaged Lagrangian is given by   Ld(q0,q1,h)=∫0h12(q˙A(t)2−qA(t)2)dt−∫0hϵ3qA(t)3dt, (4.6) where $$(q_A(t),\dot{q}_A(t))$$ is the solution corresponding to the Lagrangian $$L^{(A)}(q,\dot{q}) = \frac{1}{2}(\dot{q}^2 - q^2)$$ with boundary conditions $$(q_0,q_1)$$. Analogously, the averaged Hamiltonian is given by   Hd+(q0,p1,h)=p1qA(h)−∫0h12(pA(t)2−qA(t)2)dt+ϵ3∫0hqA(t)3dt, (4.7) where $$(q_A(t),p_A(t))$$ is the solution corresponding to the Hamiltonian $$H^{(A)}(q,p)$$ with boundary conditions $$(q_0,p_1)$$. Applying the discrete right and left Legendre transforms implicitly defines the discrete Hamiltonian map for $$L_d(q_0,q_1,h)$$ and the discrete right Hamiltonian map for $$H_d^+(q_0,p_1,h)$$, which yields the respective one-step methods. Numerical simulations were run over a time span from 0 to 10000 or the nearest integer value to 10000 for the respective time step. The initial conditions are given by $$(q_0,p_0) = (1,0)$$. Figures 2 and 3 show plots of the energy error vs. step size for two different values of $$\epsilon$$. The third plot in each of the figures hints that the discrete Lagrangian and discrete right Hamiltonian have numerical resonance that is nearly dual, in some sense, with respect to step size. The discrete Lagrangian exhibits excessive numerical resonance for step sizes near odd multiples of $$\pi$$, while the discrete right Hamiltonian exhibits excessive numerical resonance for step sizes near odd multiples of $$\frac{\pi}{2}$$. It should be noted that the arbitrary value of $$10^6$$ was substituted for output that was either near infinite or NaN. What is particularly striking is that the occurrence of the numerical resonance is intimately connected to the corresponding boundary values for each generating function. Fig. 2. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.1$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance around odd integer multiples of $$\frac{\pi}{2}$$ and exactly at odd multiples of $$\pi$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 2. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.1$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance around odd integer multiples of $$\frac{\pi}{2}$$ and exactly at odd multiples of $$\pi$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 3. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.001$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance at some odd integer multiples of $$\frac{\pi}{2}$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 3. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.001$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance at some odd integer multiples of $$\frac{\pi}{2}$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Now this by no means provides a rigorous analysis of the numerical resonances, nor does it fully explain all of the resonance effects, but it does provide motivation and insight into the numerical differences between the discrete Lagrangian and discrete right Hamiltonian. A more in-depth analysis might be provided by applying something similar to modulated Fourier expansions (see Hairer and Lubich, 2001, 2012; Hairer et al., 2006, Chapter XIII). Modulated Fourier expansions are particularly well suited for oscillatory problems when large step sizes are sought. The standard backward error analysis relies on $$h\omega \to 0$$, which is not the case for high oscillatory problems when seeking large step sizes. Modulated Fourier expansions can provide a tool for deriving many of the same results as backward error analysis, such as long-term energy preservation. Furthermore, it can be quite useful for examining the step sizes that lead to excessive numerical resonance. However, it should be noted that while modulated Fourier expansions have been used quite successfully to analyse explicit trigonometric integrators, it is not quite as clear how easily it can deal with implicit integrators such as those obtained from the discrete averaged Lagrangian and discrete averaged Hamiltonian. 5. Conclusion Error analysis and symmetry results have now been extended to cover discrete Hamiltonian variational integrators. Furthermore, examples have been presented indicating that the underlying well-posedness in terms of the boundary conditions of the exact generating function can be directly related to numerical resonance. In conclusion, it is clear that the numerical properties of variational integrators are dependent on both the approximation scheme used in constructing the generating function and the type of generating function being approximated. This article indicates that the class of variational integrators generated using the Hamiltonian formulation are not necessarily equivalent to the ones obtained from the Lagrangian formulation, and it would therefore be of interest to continue developing methods based on the discrete Hamiltonian variational integrator formulation. In particular, the results presented suggest that further work remains to be done to better understand the circumstances under which it is preferable to favor one approach over the other. Funding National Science Foundation grants DMS-1010687, CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, DMS-1345013. References De León, M., Martín de Diego, D. & Santamaría-Merino, A. ( 2007) Discrete variational integrators and optimal control theory. Adv. Comput. Math. , 26, 251– 268. Google Scholar CrossRef Search ADS   Farr, W. M. ( 2009) Variational integrators for almost-integrable systems. Celest. Mech. Dyn. Astron. , 102, 105– 118. Google Scholar CrossRef Search ADS   Hairer, E. & Lubich, C. ( 2001) Long-time energy conservation of numerical methods for oscillatory differential equations. SIAM J. Numer. Anal. , 38, 414– 441. Google Scholar CrossRef Search ADS   Hairer, E. & Lubich, C. ( 2012) Modulated Fourier expansions for continuous and discrete oscillatory systems. Foundations of Computational Mathematics, Budapest 2011  ( Cucker, F. Krick, T. Pinkus A. & Szanto A. eds). London Mathematical Society Lecture Note Series. Cambridge: Cambridge University Press, pp. 113– 128. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric Numerical Integration,  2nd edn. Springer Series in Computational Mathematics, vol. 31. Berlin: Springer. Lall, S. & West, M. ( 2006) Discrete variational Hamiltonian mechanics. J. Phys. A , 39, 5509– 5519. Google Scholar CrossRef Search ADS   Leok, M. & Shingel, T. ( 2012) General techniques for constructing variational integrators. Front. Math. China , 7, 273– 303. Google Scholar CrossRef Search ADS   Leok, M. & Zhang, J. ( 2011) Discrete Hamiltonian variational integrators. IMA J. Numer. Anal. , 31, 1497– 1532. Google Scholar CrossRef Search ADS   Lew, A., Marsden, J. E., Ortiz, M. & West, M. ( 2003) Asynchronous variational integrators. Arch. Ration. Mech. Anal. , 167, 85– 146. Google Scholar CrossRef Search ADS   Leyendecker, S., Marsden, J. E. & Ortiz, M. ( 2008) Variational integrators for constrained mechanical systems. Z. Angew. Math. Mech. , 88, 677– 708. Google Scholar CrossRef Search ADS   Marsden, J. E., Patrick, G. W. & Shkoller, S. ( 1998) Multisymplectic geometry, variational integrators, and nonlinear PDEs. Comm. Math. Phys. , 199, 351– 395. Google Scholar CrossRef Search ADS   Marsden, J. E. & West, M. ( 2001) Discrete mechanics and variational integrators. Acta Numer. , 10, 357– 514. Google Scholar CrossRef Search ADS   Verhulst, F. ( 2000) Non-linear Differential Equations and Dynamical Systems . Berlin: Springer. Wisdom, J. & Holman, M. ( 1991) Symplectic maps for the $$n$$-body problem. Astron. J. , 102, 1528– 1538. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Properties of Hamiltonian variational integrators

, Volume 38 (1) – Jan 1, 2018
22 pages

/lp/ou_press/properties-of-hamiltonian-variational-integrators-cuMmaVO1wx
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx010
Publisher site
See Article on Publisher Site

### Abstract

Abstract Discrete Hamiltonian variational integrators are derived from type II and type III generating functions for symplectic maps, and in this article, we establish a variational error analysis result that relates the order of accuracy of the associated numerical methods with the extent to which these generating functions approximate the exact discrete Hamiltonians. We also introduce the notion of an adjoint discrete Hamiltonian and relate it to the adjoint of the associated symplectic integrator. We show that when constructing discrete Lagrangians and discrete Hamiltonians using the same approximation method, this does not necessarily lead to the same symplectic integrator. Numerical experiments also indicate that the resonance behavior of variational integrators based on averaging methods depends on the type of generating functions used, and we relate this resonance behavior to the ill-posedness of the boundary-value problems used to define the exact discrete Lagrangian and exact discrete Hamiltonian. 1. Introduction Geometric numerical integration is a field of numerical analysis that develops numerical methods with the goal of preserving geometric properties of dynamical systems (see Hairer et al., 2006). Variational integrators are geometric numerical integrators derived from discretizing Hamilton’s principle from classical mechanics (see Marsden and West, 2001). They have many desirable properties such as symplecticity, momentum preservation and near-energy preservation, which results in excellent long-term stability. While the Lagrangian formulation of variational integrators has been thoroughly investigated (see Marsden et al., 1998; Marsden and West, 2001; Lew et al., 2003; de León et al., 2007; Leyendecker et al., 2008; Leok and Shingel, 2012), only recently has the Hamiltonian formulation of variational integrators been established (see Lall and West, 2006; Leok and Zhang, 2011). In this article, we will continue the investigation of Hamiltonian variational integrators and establish theorems on error analysis and symmetry of the method, and provide numerical experiments to elucidate the relative numerical advantages and disadvantages of the Lagrangian and Hamiltonian formulations. In particular, evidence is presented to show that for oscillatory problems the discrete Lagrangian and discrete Hamiltonian variational integrators have differing resonance and conditioning properties. In addition, it is shown that some approximation methods will yield a symmetric method only when derived from a specific type of generating function. The upshot is that the numerical properties of a variational integrator are determined both by the approximation scheme used to construct it and by the type of the generating function being approximated. 1.1 Discrete mechanics Lagrangian variational integrators are based on a discrete analog of Hamilton’s principle, whereas Hamiltonian variational integrators are based on a discrete analog of Hamilton’s phase space variational principle. The fundamental objects in the discretization are generating functions of symplectic maps, and in the Hamiltonian case, they are obtained by approximating the exact type II generating function associated with a Hamiltonian flow, which we refer to as the exact discrete right Hamiltonian,   Hd+,E(q0,p1)=ext(q,p)∈C2([0,T],T∗Q)q(0)=q0,p(T)=p1(p1q1−∫0T[pq˙−H(q,p)]dt). (1.1) This can be viewed as the solution at time $$T$$ of the type II Hamilton–Jacobi equation   ∂S2(q0,p,t)∂t=H(∂S2∂p,p), (1.2) which more generally describes the type II generating function that generates the time-$$t$$ Hamiltonian flow map,   S2(q0,p,t)=ext(q,p)∈C2([0,t],T∗Q)q(0)=q0,p(t)=p(p(t)q(t)−∫0t[p(s)q˙(s)−H(q(s),p(s))]ds). (1.3) Similarly, the exact discrete left Hamiltonian is given by   Hd−,E(p0,q1;h)=ext(q,p)∈C2([0,T],T∗Q)q(0)=q0,p(T)=p1(−p0q0−∫0T[pq˙−H(q,p)]dt) (1.4) and it can be viewed as a solution at time $$T$$ of the type III Hamilton–Jacobi equation   ∂S3(p0,q,t)∂t=H(q,−∂S3∂q). (1.5) Given discrete Hamiltonians $$H_d^+(q_k,p_{k+1};h)$$ and $$H_d^-(p_k,q_{k+1};h)$$, the discrete Hamilton’s equations are given by   qk+1 =D2Hd+(qk,pk+1;h), (1.6)  pk =D1Hd+(qk,pk+1;h) (1.7) and   qk =−D1Hd−(pk,qk+1;h), (1.8)  pk+1 =−D2Hd−(pk,qk+1;h). (1.9) These can also be expressed in terms of the discrete Legendre transformations $$\mathbb{F}^\pm H_d^+:(q_k,p_{k+1})\rightarrow T^*Q$$,   F+Hd+(qk,pk+1;h) =(D2Hd+(qk,pk+1;h),pk+1), (1.10)  F−Hd+(qk,pk+1;h) =(qk,D1Hd+(qk,pk+1;h)) (1.11) and $$\mathbb{F}^\pm H_d^-:(p_k,q_{k+1})\rightarrow T^*Q$$,   F+Hd−(pk,qk+1;h) =(qk+1,−D2Hd−(pk,qk+1;h)), (1.12)  F−Hd−(pk,qk+1;h) =(−D1Hd−(pk,qk+1;h),pk). (1.13) We observe that the Hamiltonian maps $$\tilde{F}_{H_d^\pm}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$ can be expressed as   F~Hd±=F+Hd±∘(F−Hd±)−1. (1.14) 2. Error analysis and symmetric methods 2.1 Error analysis Variational integrators are able to benefit from and adopt many traditional techniques and methods of numerical analysis (see Leok and Shingel, 2012). This can be largely attributed to the following theorem from Marsden and West (2001). Theorem 2.1 (Marsden and West, 2001, Theorem 2.3.1). If a discrete Lagrangian, $$L_{\rm d}:Q \times Q \rightarrow \mathbb{R}$$ approximates the exact discrete Lagrangian $$L_{\rm d}^{\rm E}:Q \times Q \rightarrow \mathbb{R}$$ to order $$r$$, i.e.,   Ld(q0,q1;h)=LdE(q0,q1;h)+O(hr+1), then the discrete Hamiltonian map $$\tilde{F}_{L_d}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$, viewed as a one-step method, is order-$$r$$ accurate. Thus, in order to generate a variational integrator of a particular order, one can leverage techniques from numerical analysis with the goal of approximating the exact discrete Lagrangian, then the associated discrete Hamiltonian map yields the variational integrator. We first present the corresponding theorem for discrete Hamiltonian variational integrators, which draws much of its inspiration from the theorem and proof of the above result as detailed in Marsden and West (2001). Theorem 2.2 If a discrete right Hamiltonian $$H^+_d:T^*Q \rightarrow \mathbb{R}$$ approximates the exact discrete right Hamiltonian $$H_{\rm d}^{+,{\rm E}}:T^*Q \rightarrow \mathbb{R}$$ to order $$r$$, i.e.,   Hd+(q0,p1;h)=Hd+,E(q0,p1;h)+O(hr+1), and the Hamiltonian is continuously differentiable then the discrete map $$\tilde{F}^h_{H_d^+}:(q_k,p_k) \mapsto (q_{k+1},p_{k+1})$$, viewed as a one-step method, is order-$$r$$ accurate. Note that the following proof can be easily adjusted to prove an equivalent theorem for the discrete left Hamiltonian case. First, we will need the following lemma. Lemma 2.3 Let $$f_1,g_1,e_1,f_2,g_2,e_2 \in C^r$$ be such that   f1(x,h) =g1(x,h)+hr+1e1(x,h),f2(x,h) =g2(x,h)+hr+1e2(x,h). Then, there exist functions $$e_{12}$$ and $$\bar{e}_1$$ bounded on compact sets such that   f2(f1(x,h),h) =g2(g1(x,h),h)+hr+1e12(g1(x,h),h),f1−1(y) =g1−1(y)+hr+1e¯1(y). Proof. We have   f2(f1(x,h),h) =f2(g1(x,h)+hr+1e1(x,h),h) =g2(g1(x,h)+hr+1e1(x,h),h)+hr+1e2(g1(x,h)+hr+1e1(x,h),h) =g2(g1(x,h),h)+hr+1e~1(g1(x,h),h)+hr+1e2(g1(x,h)+hr+1e1(x,h),h), where $$\tilde{e}_1$$ is bounded on a compact set. This last line comes from combining the compactness of the set with the smoothness of the functions to obtain a Lipschitz property of the form   ‖g2(g1(x,h)+hr+1e1(x,h),h)−g2(g1(x,h),h)‖≤Chr+1. For each choice of $$(x,h)$$, equality holds for a particular choice of constant, which defines $$\tilde{e}_1$$ and establishes its smoothness as a function. Adding $$e_2$$ to $$\tilde{e}_1$$ we obtain a function $$e_{12}$$, which is also bounded on compact sets such that   f2(f1(x,h),h)=g2(g1(x,h),h)+hr+1e12(g1(x,h),h). Let $$y=f_1(x,h)$$ and note that by definition,   f1−1(f1(x,h))=g1−1(g1(x,h)). Since $$g_1^{-1}(y)=g_1^{-1}(g_1(x,h)+h^{r+1}e_1(x,h))$$, then   ‖g1−1(y)−f1−1(y)‖=‖g1−1(y)−g1−1(g1(x,h))‖≤C¯hr+1. From this, it follows that there exists a function $$\bar{e}_1$$ bounded on compact sets such that   f1−1(y)=g1−1(y)+hr+1e¯1(y). □ Now we are ready for the proof of the theorem. Proof. By assumption there is some bounded continuously differentiable function $$e$$ such that   Hd+(q(0),p(h);h)=Hd+,E(q(0),p(h),h)+hr+1e(q(0),p(h),h). Differentiating yields   D1Hd+(q(0),p(h);h)=D1Hd+,E(q(0),p(h);h)+hr+1D1e(q(0),p(h),h), where $$\|D_1e(q(0),p(h),h)\| \leq \tilde{C}$$. This implies   ‖F−Hd+(q(0),p(h);h)−F−Hd+,E(q(0),p(h);h)‖≤C~hr+1. Now combining this with the fact that $$\tilde{F}_{H^+_d}=\mathbb{F}^+H^+_d \circ (\mathbb{F}^-H^+_d)^{-1}$$ and applying Lemma 2.3, we have   F~Hd+h=F~Hd+,Eh+O(hr+1). □ Determining the order of a variational integrator is greatly simplified via the above theorems, which relate the order of the integrator to the order to which the associated discrete Lagrangian or discrete right Hamiltonian approximates the corresponding exact generating function. Similarly, it was shown in Marsden and West (2001) that one can determine whether the variational integrator is a symmetric method by examining the corresponding discrete Lagrangian. We would like to extend this result to the case of discrete Hamiltonians. 2.2 Symmetric methods Definition 2.4 (See Hairer et al., 2006, Chapters II.3 and V). A numerical one-step method $${\it{\Phi}}_h$$ is called symmetric or time reversible if it satisfies   Φh∘Φ−h=id or equivalently   Φh=Φ−h−1. The adjoint of a numerical one-step method, denoted $${\it{\Phi}}_h^*$$, is defined as   Φh∗=Φ−h−1. A numerical one-step method is a symmetric method if it is self-adjoint, i.e., $${\it{\Phi}}_h={\it{\Phi}}_h^*$$. The adjoint of a discrete Lagrangian $$L_d^*$$ is defined as   Ld∗(q0,q1,h)=−Ld(q1,q0,−h). The discrete Lagrangian is called self-adjoint if $$L_d^*(q_0,q_1,h)=L_d(q_0,q_1,h)$$. The following theorem from Marsden and West (2001) relates the self-adjointness of the discrete Lagrangian with the self-adjointness of the corresponding variational integrator. Theorem 2.5 (Marsden and West, 2001, Theorem 2.4.1). The discrete Lagrangian (or an equivalent discrete Lagrangian) $$L_d$$ is self-adjoint if and only if the method associated with the corresponding discrete Hamiltonian map is self-adjoint, i.e., symmetric. In many cases, it is easier to check whether the discrete Lagrangian is self-adjoint rather than checking the variational integrator itself. We seek a definition for the adjoint of a discrete right Hamiltonian. The adjoint of a one-step method $$(q_1,p_1)={\it{\Phi}}_h(q_0,p_0)$$ can be obtained by reversing the direction of time and reversing the roles of the initial data and terminal solution, i.e., $$(q_0,p_0)={\it{\Phi}}^*_{-h}(q_1,p_1)$$. This corresponds to swapping out $$(q_0,p_0,q_1,p_1,h)$$ for $$(q_1,p_1,q_0,p_0,-h)$$. This motivates the definition of the adjoint of a type II generating function as a type III generating function and vice versa. In particular, given a type II discrete Hamiltonian $$H_d^+$$, we seek a definition for the type III adjoint $$(H_d^+)^*$$ that will satisfy $$\tilde{F}^h_{(H_d^+)^*} = (\tilde{F}^h_{H_d^+})^*$$. Let $$\tilde{F}^h_{(H_d^+)^*}(q_0,p_0)=(q_1,p_1)$$. Then, we want   (q1,p1) =F~(Hd+)∗h(q0,p0) =(F~Hd+h)∗(q0,p0) =(F~Hd+−h)−1(q0,p0). This implies $$\tilde{F}^{-h}_{H_d^+}(q_1,p_1) = (q_0,p_0)$$, which, together with $$\tilde{F}^h_{(H_d^+)^*}(q_0,p_0)=(q_1,p_1),$$ yields the respective sets of equations   p1 =D1Hd+(q1,p0;−h),q0 =D2Hd+(q1,p0;−h) and   p1 =−D2(Hd+)∗(p0,q1;h),q0 =−D1(Hd+)∗(p0,q1;h). Comparing these equations we see that setting $$(H_d^+)^*(p_0,q_1;h) = -H_d^+(q_1,p_0;-h)$$ satisfies $$\tilde{F}^h_{(H_d^+)^*} = (\tilde{F}^h_{H_d^+})^*$$. A similar calculation yields an analogous expression for the adjoint of a type III generating function $$H_d^-$$. Definition 2.6 Given a type II/III generating function, $$H_d^\pm$$, define the adjoint as the type III/II generating function, $$(H_d^\pm)^*$$, where $$\tilde{F}_{(H_d^\pm)^*}^h(q_0,p_0) = (q_1,p_1)$$, as   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h), (2.1)  (Hd−)∗(q0,p1;h) =−Hd−(p1,q0;−h). (2.2) Example 2.7 The symplectic Euler-A method for a Lagrangian of the form $$L(q,\dot{q}) = \frac{1}{2}\dot{q}M\dot{q} - V(q)$$ is given by   p1 =p0−h∇V(q0),q1 =q0+hM−1p1. The corresponding discrete right Hamiltonian is given by   Hd+(q0,p1,h) =p1(q0+hM−1p1)−h[p1M−1p1−H(q0,p1)] =p1q0+hH(q0,p1). The adjoint of this method is given by symplectic Euler-B   q1=q0+hM−1p0,p1=p0−h∇V(q1). We now derive the corresponding adjoint of the discrete right Hamiltonian for symplectic Euler-A:   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h) =−p0(q1−hM−1p0)−h[p0M−1p0−H(q1,p0)] =−p0q1+hH(q1,p0). We can verify that this generates symplectic Euler-B by applying the discrete left Hamilton equations   q0 =−D1(Hd+)∗(p0,q1;h) =D2Hd+(q1,p0;−h) =q1−hM−1p0,p1 =−D2(Hd+)∗(p0,q1;h) =D1Hd+(q1,p0;−h) =p0−h∇V(q1). Solving the first equation for $$q_1$$ gives symplectic Euler-B, as expected. Theorem 2.8 $$(H_d^\pm)^{**}=H_d^\pm$$. Proof. We consider the case of the type II generating function $$H_d^+$$. Let $$\tilde{F}^h_{(H_d^+)^{**}}(q_0,p_0)=(q_1,p_1)$$. Since $$(H_d^+)^*$$ is a type III generating function, applying the definition of the adjoint twice gives   (Hd+)∗∗(q0,p1;h) =−(Hd+)∗(p1,q0;−h) =Hd+(q0,p1;h) and a similar calculation shows that this holds for the type III generating function $$H_d^-$$ as well. □ Since the notion of the adjoint that we introduced converts a type II to a type III generating function, for a discrete Hamiltonian to be self-adjoint, we need to compare the adjoint to the Legendre transformation of the discrete Hamiltonian, which is given by   Hd−(pk,qk+1;h)=−pkqk−pk+1qk+1+Hd+(qk,pk+1;h), where we view $$p_{k+1}$$ and $$q_k$$ as functions of $$p_k$$ and $$q_{k+1}$$. Then, the following calculation shows that these two generating functions generate the same symplectic map, i.e., $$\tilde{F}_{H_d^-}=\tilde{F}_{H_d^+}$$,   −D1Hd−(pk,qk+1;h) =qk+pk∂qk∂pk+∂pk+1∂pkqk+1−D1Hd+(qk,pk+1;h)∂qk∂pk−D2Hd+(qk,pk+1;h)∂pk+1∂pk =qk+(pk−D1Hd+(qk,pk+1;h))∂qk∂pk+(qk+1−D2Hd+(qk,pk+1;h))∂pk+1∂pk,−D2Hd−(pk,qk+1;h) =pk∂qk∂qk+1+∂pk+1∂qk+1qk+1+pk+1−D1Hd+(qk,pk+1;h)∂qk∂qk+1 −D2Hd+(qk,pk+1;h)∂pk+1∂qk+1 =pk+1+(pk−D1Hd+(qk,pk+1;h))∂qk∂qk+1+(qk+1−D2Hd+(qk,pk+1;h))∂pk+1∂qk+1. Definition 2.9 A type II/III generating function is self-adjoint if it is equal (up to equivalency) to the Legendre transform of its adjoint. Note that this definition implies that a discrete right Hamiltonian is self-adjoint if its adjoint is equal (up to equivalency) to the associated discrete left Hamiltonian, i.e., $$(H_d^+)^* = H_d^-$$. Theorem 2.10 Given a self-adjoint discrete right Hamiltonian, i.e., $$H_d^-=(H_d^+)^*$$, the method associated with the discrete right Hamiltonian map is self-adjoint. Likewise, if a method coming from a discrete right Hamiltonian map is self-adjoint, then the associated discrete right Hamiltonian is self-adjoint. Proof. Assume $$H_d^-=(H_d^+)^*$$. Then,   (F~Hd+)∗=F~(Hd+)∗=F~Hd−=F~Hd+, and so, by definition, the map is self-adjoint. Now assume $$\tilde{F}_{H_d^+}=(\tilde{F}_{H_d^+})^*$$. Then,   F~Hd−=F~Hd+=(F~Hd+)∗=F~(Hd+)∗, which implies $$(H_d^+)^*=H_d^-$$ (up to equivalency) and, by definition, the discrete right Hamiltonian is self-adjoint. □ A similar proof can be used to prove an identical theorem for the discrete left Hamiltonian. The previous theorem provides an easy way to check whether a variational integrator is self-adjoint. Assuming the Hamiltonian flow is time reversible, it follows that the exact discrete right Hamiltonian is self-adjoint. This can also be shown using the definition of a self-adjoint exact discrete right Hamiltonian, and the same result can be shown for the discrete left Hamiltonian with only minor adjustments. Corollary 2.11 The exact discrete right Hamiltonian $$H_{\rm d}^{+,{\rm E}}$$ is self-adjoint. Proof. A direct calculation shows that   (Hd+,E)∗(p0,q1;h) =−Hd+,E(q1,p0;−h) =−(p~(−h)q~(−h)−∫0−h[p~(τ)q~(τ)−H(q~(τ),p~(τ))]dτ) =−p(−h+h)q(−h+h)−∫−h0[p(τ+h)q(τ+h)−H(q(τ+h),p(τ+h))]dτ =−p(0)q(0)−∫0h[p(t)q(t)−H(q(t),p(t))]dt =Hd−,E(p0,q1;h), where we used the fact that the time-reversed solution $$(\tilde{q}(\tau),\tilde{p}(\tau))$$ over the time domain $$[-h,0]$$ with $$(q_1,p_0)$$ boundary data is related to the solution curve $$(q(t),p(t))$$ over the time domain $$[0,h]$$ with $$(q_0,p_1)$$ boundary data by $$(\tilde{q}(\tau),\tilde{p}(\tau))=(q(\tau+h),p(\tau+h))$$. □ The definition of the adjoint also provides a simple way to construct symmetric methods. Given any method defined by $$H_d$$, we can construct a symmetric method using composition, for example, $$\tilde{F}^{{h}/{2}}_{H_d} \circ \tilde{F}^{{h}/{2}}_{H_d^*}$$, which is nothing more than composing a half-step of the adjoint method with a half-step of the method. It is well known that this leads to a symmetric method, as the following calculation demonstrates,   (F~Hdh/2∘F~Hd∗h/2)∗ =(F~Hd∗h/2)∗∘(F~Hdh/2)∗ =F~Hd∗∗h/2∘F~Hd∗h/2 =F~Hdh/2∘F~Hd∗h/2, where the last line used Theorem 2.8. More generally, a composition method of the form,   F~Hdαsh∘F~Hd∗βsh∘⋯∘F~Hd∗β2h∘F~Hdα1h∘F~Hd∗β1h, where $$\alpha_{s+1-i}=\beta_i$$ for $$i=1,\ldots,s$$, will be symmetric. For a more in-depth discussion of symmetric composition methods, see Hairer et al. (2006, Chapter V.3). 3. Discrete Lagrangians vs. discrete Hamiltonians 3.1 Composition of discretization and the Legendre transform A transformation of one type of generating function into another type of generating function, which preserves the associated sympletic map, is given by the following Legendre transforms. Definition 3.1 (i) The Legendre transform of a type I generating function into a type II generating function is given by the equation   Hd+(qk,pk+1;h)=pk+1qk+1−Ld(qk,qk+1;h), where $$q_{k+1}$$ is implicitly defined by $$p_{k+1}=D_2L_d(q_k,q_{k+1};h)$$. The Legendre transform of a type II generating function into type I is given by the same equation, where $$p_{k+1}$$ is implicitly defined by $$q_{k+1}=D_2H^+_d(q_k,p_{k+1};h)$$. (ii) The Legendre transform of a type I generating function into a type III generating function is given by the equation   Hd−(pk,qk+1;h)=−pkqk−Ld(qk,qk+1;h), where $$q_k$$ is implicitly defined by $$p_k=-D_1L_d(q_k,q_{k+1};h)$$. The Legendre transform of a type III generating function into type I is given by the same equation, where $$p_k$$ is implicitly defined by $$q_k=-D_1H^-_d(p_k,q_{k+1};h)$$. (iii) The Legendre transform of a type II generating function into a type III generating function is given by the equation   Hd−(pk,qk+1;h)=−pkqk−pk+1qk+1+Hd+(qk,pk+1;h), where $$q_k$$ and $$p_{k+1}$$ are implicitly defined by the set of equations $$p_k=D_1H^+_d(q_k,p_{k+1};h)$$ and $$q_{k+1}=D_2H^+_d(q_k,p_{k+1};h)$$. The Legendre transform of a type III generating function into type II is given by the same equation, where $$q_{k+1}$$ and $$p_k$$ are implicitly defined by the set of equations $$q_k=-H^-_d(p_k,q_{k+1};h)$$ and $$p_{k+1}=-H^-_d(p_k,q_{k+1})$$. Variational integrators are derived by first discretizing the exact generating function, which results in a new generating function of the same type. Using the Legendre transforms defined above, the discretized generating function can be transformed into an equivalent generating function of a different type. This can be viewed as composing the Legendre transform and the discretization. Likewise, one could first take the Legendre transform of the given exact generating function and then apply the discretization. The question we address next is whether changing the order of this composition results in the same symplectic map, and for one particular type of discretization this question has already been answered. It was shown in Leok and Zhang (2011) that the Galerkin variational integrator construction leads to equivalent discrete Lagrangian and discrete Hamiltonian methods for the same choice of quadrature rule and finite-dimensional function space, and the result is given in the following theorem. Theorem 3.2 (Leok and Zhang, 2011, Proposition 4.1). If the continuous Hamiltonian $$H(q,p)$$ is hyperregular and we construct a Lagrangian $$L(q,\dot{q})$$ by the Legendre transformation then the generalized Galerkin Hamiltonian variational integrator (see Leok and Zhang, 2011) and the generalized Galerkin Lagrangian variational integrator, associated with the same choice of basis functions and numerical quadrature formula, are equivalent. Will this hold for other types of variational integrators? To begin to address this question, we propose the following examples. Example 3.3 Consider a Lagrangian of the form $$L(q,\dot{q})=\frac{1}{2}\dot{q}M\dot{q} - V(q)$$, where $$M$$ is symmetric positive definite and $$V$$ is sufficiently smooth. The exact discrete Lagrangian, which is defined as   LdE(q0,q1;h) =∫0hL(q01(t),q˙01(t))dt, where $$q_{01}(0)=q_0,$$$$q_{01}(h)=q_1,$$ and $$q_{01}$$ satisfies the Euler–Lagrange equation in the time interval $$(0,h)$$. Letting $$q_0$$ and $$q_1$$ be fixed, then a first-order finite difference approximation of the velocities yields   q˙0 ≈q1−q0h,q˙1 ≈q1−q0h. Using the rectangular quadrature rule about the initial point results in the following discrete Lagrangian (i.e., type I generating function):   Ld(q0,q1;h)=12(q1−q0h)M(q1−q0h)−V(q0). (3.1) Applying the implicit discrete Euler–Lagrange equations yields   p0=Mq1−q0h+h∇V(q0),p1=Mq1−q0h. Finally rearranging these equations results in the variational integrator, also known as symplectic Euler-A,   q1=q0+hM−1p1,p1=p0−h∇V(q0). Now we apply this same approximation scheme to the associated exact type II generating function. The boundary-value formulation of the exact discrete right Hamiltonian is given by   Hd+,E(q0,p1) =(p1q1−∫0h[pq˙−H(q,p)]dt), where $$(q(t),p(t))$$ satisfy Hamilton’s equations with boundary conditions $$q(0)=q_0$$, $$p(h)=p_1$$. The first-order finite difference approximations of the velocities yield the following equations for the momentum:   p0 ≈Mq1−q0h,p1 ≈Mq1−q0h. Solving for $$q_1$$ in terms of $$q_0$$ and $$p_1$$ yields the approximation $$q_1\approx q_0 + hM^{-1}p_1$$, and this simplifies the approximation to $$p_0$$ as $$p_0 \approx p_1$$. Now using these approximations in combination with the rectangular rule about the initial point yields the discrete right Hamiltonian   Hd+(q0,p1;h)=p1(q0+hM−1p1)−h(12p1M−1p1−V(q0)). After applying the implicit discrete Hamilton equations and rearranging terms, the resulting method is again symplectic Euler-A:   q1=q0+hM−1p1p1=p0−h∇V(q0). The composition of the Legendre transform and the discretization of a generating function will be called commutative if, regardless of the order of this composition, either of the resulting generating functions leads to the same symplectic map. In this context, the previous example shows that this particular discretization scheme, which includes both the trajectory approximation and the quadrature rule, commutes with the Legendre transform between type I and type II generating functions. Now we look at the same velocity approximation but using the rectangular rule about the end point rather than the initial point. Example 3.4 As before, we first build the discrete Lagrangian with $$q_0$$, $$q_1$$ fixed and velocity approximations   q˙0 ≈q1−q0h,q˙1 ≈q1−q0h. Applying the rectangular rule about the end point yields the discrete Lagrangian   Ld(q0,q1;h)=12(q1−q0h)M(q1−q0h)−V(q1) (3.2) and the associated variational integrator is symplectic Euler-B:   q1=q0+hM−1p0,p1=p0−h∇V(q1). Now applying the same approximation scheme to construct the discrete right Hamiltonian results in the following type II generating function:   Hd+(q0,p1;h)=p1(q0+hM−1p1)−h(12p1M−1p1−V(q0+hM−1p1)). Applying the implicit discrete Hamilton equations and rearranging the terms results in the method   q1=q0+hM−1p0,p1=p0−h∇V(q0+hM−1p1), which is not symplectic Euler-B. We have just proven the following theorem. Theorem 3.5 In general, the composition of the Legendre transform of a generating function and the discretization of a generating function do not commute. Therefore, the answer to our original question is that in general, a fixed approximation scheme used to construct a discrete Lagrangian will not generate the same method when it is used to construct a discrete Hamiltonian. In general, how might the two resulting methods differ? A complete characterization of this issue is subtle, and beyond the scope of this article, but it will be a topic of future work. For now, we will consider how the two approaches differ when combined with the method of averaging, which will also serve to illustrate how the type of generating function and the associated boundary data can affect the numerical properties of the method. 3.2 Averaged Lagrangians and Hamiltonians Averaging methods have played a role in solving differential equations since at least as far back as the time of Lagrange (see Verhulst, 2000), and they continue to play a key role particularly in the field of numerical differential equations applied to nearly integrable systems or problems with multiple time scales. We consider perturbed Hamiltonian systems with Hamiltonians of the form   H=H(A)+ϵH(B), (3.3) where $$\epsilon \ll 1$$ and the dynamics of the Hamiltonian system corresponding to $$H^{(A)}$$ is exactly solvable or at the very least cheap to approximate. We call this an almost-integrable system, the motivation being that the dynamics of the system are largely influenced by an integrable Hamiltonian with simpler dynamics, but smaller influences also play a role in the overall dynamics. An example is the classic $$n$$-body problem of the solar system, where a particular planet’s trajectory is largely influenced by the sun, but other planets and nearby objects also play a role. Averaging methods can be constructed to exploit the larger influence of $$H^{(A)}$$ on the dynamics of the system by averaging out the smaller influences. Ideally, averaging techniques will allow for larger time steps to be used while still yielding a reasonable approximation to the solution. A variational integrator for such a system was proposed in Farr (2009) using a discrete Lagrangian formulation, which drew inspiration from the kick-drift-kick leapfrog method (see Wisdom and Holman, 1991). We will discuss the original Lagrangian formulation (hereafter referred to as the averaged Lagrangian) and in addition construct an analogous method in terms of a discrete right Hamiltonian (referred to as the averaged Hamiltonian). The Lagrangian corresponding to (3.3) is given by   L=L(A)+ϵL(B) (3.4) and we will make the assumption that $$L^{(B)}(q(t),\dot{q}(t))=-V^{(B)}(q(t))$$. The averaging method of interest has a local truncation error of $$\mathcal{O}(\epsilon^2 h^3)$$ and is defined in terms of a discrete Lagrangian, $$L_d$$. This method, proposed in Farr (2009), uses a discrete Lagrangian of the form   Ld(q0,q1,h) =Ld(A),E(q0,q1;h)+ϵ∫0hL(B)(qA(q0,q1,t),q˙A(q0,q1,t))dt =Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt, where we denote the trajectory corresponding to $$L^{(A)}$$ with boundary conditions $$(q_0,q_1)$$ by $$(q_A(t),\dot{q}_A(t))$$. The idea is to use the dynamics of $$L^{(A)}$$, which is solved for either exactly, or efficiently and accurately approximated, to average the contribution of $$L^{(B)}$$ to the dynamics. The corresponding discrete Hamiltonian map is given implicitly by   −p0 =D1Ld(A),E(q0,q1;h)−ϵ∫0hD1V(B)(qA(q0,q1,t))dt, (3.5a)  p1 =D2Ld(A),E(q0,q1;h)−ϵ∫0hD2V(B)(qA(q0,q1,t))dt. (3.5b) As shown in Farr (2009), this method has local truncation error of size $$\mathcal{O}(\epsilon^2 h^3)$$. Using the notation $$p_0^A(q_0,q_1)=-D_1L_d^{(A),E}(q_0,q_1;h)$$ and $$p_1^A(q_0,q_1)=D_2L_d^{(A),E}(q_0,q_1;h)$$, we rearrange the above equations to get   p0−ϵ∫0hD1V(B)(qA(q0,q1,t))dt=p0A(q0,q1), (3.6a)  p1=p1A(q0,q1)−ϵ∫0hD2V(B)(qA(q0,q1,t))dt. (3.6b) In Farr (2009), it is noted that $$- \epsilon \int_0^h D_1V^{(B)}(q_A(q_0,q_1,t)) \,{\rm d}t$$ is an average along the trajectory generated by $$L^{(A)}$$ which, in general, gives more weight to the initial periods of the trajectory, while $$- \epsilon \int_0^h D_2V^{(B)}(q_A(q_0,q_1,t)) \,{\rm d}t$$ is an average along the trajectory generated by $$L^{(A)}$$ that, in general, favors the latter periods of the trajectory. Now let us consider the discrete right Hamiltonian given by the same form of approximation,   Hd+(q0,p1;h)=Hd(A),+,E(q0,p1;h)+ϵ∫0hV(B)(qA(q0,p1,t))dt. The discrete right Hamiltonian map is given implicitly by   p0 =D1Hd(A),+,E(q0,p1;h)+ϵ∫0hD1V(B)(qA(q0,p1,t))dt,q1 =D2Hd(A),+,E(q0,p1;h)+ϵ∫0hD2V(B)(qA(q0,p1,t))dt. Using the notation $$p_0^A(q_0,p_1) = D_1 H_d^{(A),+,E}(q_0,p_1;h)$$ and $$q_1^A(q_0,p_1) = D_2 H_d^{(A),+,E}(q_0,p_1;h)$$, we rearrange the equations to yield   p0−ϵ∫0hD1V(B)(qA(q0,p1,t))dt=p0A(q0,p1), (3.7a)  q1=q1A(q0,p1)+ϵ∫0hD2V(B)(qA(q0,p1,t))dt. (3.7b) Theorem 3.6 The method defined implicitly by (3.7) has local truncation error $$\mathcal{O}(\epsilon^2 h^3)$$. Proof. Using variational error analysis, we need to show   O(ϵ2h3) =HdE,+−Hd+ =ΔA+ϵΔB, where $${\it{\Delta}}_A$$ is given by   p(h)q(h)−∫0h[p(t)q˙(t)−H(A)(q(t),p(t))]dt−(pA(h)qA(h)−∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt), and $$\epsilon{\it{\Delta}}_B$$ is given by   ϵ∫0h[V(B)(q(t))−V(B)(qA(t))]dt. Using a functional Taylor expansion, $${\it{\Delta}}_A$$ becomes   ΔA =δδqA(∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt)δqA +δ2δqA2(∫0h[pA(t)q˙A(t)−H(A)(qA(t),pA(t))]dt)δqA2+O(δqA3), where $$\delta q_A$$ is the difference between $$q$$ and $$q_A$$. Noting that $$q$$ and $$q_A$$ differ in forces of order $$\epsilon h^2$$ and $$p$$ differs from $$p_A$$ to first order in $$\epsilon h$$, implies that $$\delta q_A$$ is on the order of $$\mathcal{O}(\epsilon h)$$. This can be seen explicitly by comparing Taylor expansions about time zero. Since $$q_A$$ satisfies Hamilton’s equations for $$H^{(A)}$$, the first variation vanishes (see Leok and Zhang, 2011, Lemma 2.1) leaving a term on the order of $$h \delta q_A^2$$. Therefore, we have   ΔA=O(ϵ2h3). Likewise, a functional Taylor expansion for $${\it{\Delta}}_B$$ yields,   ΔB=δδqA[∫0hV(B)(qA(t))dt]δqA+O(δqA2). Noting that $$V^{(B)}$$ is a function of only $$q_A$$ and that $$q$$ differs from $$q_A$$ on the order of $$\epsilon h^2$$, implies $$\epsilon {\it{\Delta}}_B = \mathcal{O}(\epsilon^2 h^3)$$. □ Theorem 3.7 Assuming the flow associated with $$L^{(A)}$$ is time reversible, then both methods, defined respectively by (3.6) and (3.7), are symmetric methods. Proof. The discrete Lagrangian associated with (3.6) is given by   Ld(q0,q1;h)=Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt. The adjoint of the discrete Lagrangian is given by   (Ld(q0,q1;h))∗ =−Ld(q1,q0;−h) =−Ld(A),E(q1,q0;−h)+ϵ∫0−hV(B)(qA(q1,q0,t))dt =−Ld(A),E(q1,q0;−h)−ϵ∫0hV(B)(qA(q1,q0,t))dt =Ld(A),E(q0,q1;h)−ϵ∫0hV(B)(qA(q0,q1,t))dt =Ld(q0,q1;h). The third equality comes from the time reversibility of the flow associated with $$L^{(A)}$$, and the fourth equality uses this property together with the fact that the exact discrete Lagrangian is self-adjoint. The discrete right Hamiltonian associated with (3.7) is given by   Hd+(q0,p1;h)=Hd(A),+,E(q0,p1;h)+ϵ∫0hV(B)(qA(q0,p1,t))dt. The adjoint of the discrete right Hamiltonian is given by   (Hd+)∗(p0,q1;h) =−Hd+(q1,p0;−h) =−Hd(A),+,E(q1,p0;−h)−ϵ∫0−hV(B)(qA(q1,p0,t))dt =−Hd(A),+,E(q1,p0;−h)+ϵ∫0hV(B)(qA(q1,p0,t))dt =Hd(A),−,E(p0,q1;h)+ϵ∫0hV(B)(qA(p0,q1,t))dt =Hd−(p0,q1;h), where the third equality comes from the time reversibility of the flow associated with $$H^{(A)}$$, and the fourth equality uses this property together with the fact that the exact discrete Hamiltonian is self-adjoint. By Definition 2.9 and Theorem 2.10, the method is symmetric. □ It can be shown that methods (3.6) and (3.7) do not result in the same symplectic map. How do these respective maps differ? To gain insight into this question, we now turn to numerical experimentation. 4. Numerical results 4.1 Exact generating functions First, we consider the unperturbed harmonic oscillator boundary-value problem   q¨(t)+q(t)=0,q(0)=q0, q(h)=q1. (4.1) Analytically, the boundary-value problem is not well posed when $$h$$ is an integer multiple of $$\pi$$, and in particular there are infinitely many solutions. Recall that the exact discrete Lagrangian is given by   LdE(q0,q1;h)=∫0hL(q01(t),q˙01(t))dt, (4.2) where $$q_{01}(0)=q_0$$, $$q_{01}(h)=q_1$$, and $$q_{01}(t)$$ satisfies the Euler–Lagrange equation in the time interval $$(0,h)$$. Thus, the exact type I generating function is ultimately defined in terms of such a boundary-value problem. The integrator obtained from the exact discrete Lagrangian is given by   q1 =q0cos⁡(h)+p0sin⁡(h),p1 =q1cot⁡(h)−q0csc⁡(h). This integrator is analytically the true solution of the harmonic oscillator initial-value problem, where $$(q(0),p(0))=(q_0,p_0)$$, and consequently the local truncation error is zero. However, noting that $$\cot(h)$$ and $$\csc(h)$$ both involve dividing by $$\sin(h)$$, we expect increased round-off error for values of $$h$$ that are near integer multiples of $$\pi$$. Similarly, the exact discrete right Hamiltonian is given by   Hd+,E(q0,p1;h)=p1q1−∫0h[p01(t)q˙01(t)−H(q01(t),p01(t))]dt, (4.3) where $$q_{01}(0)=q_0$$, $$p_{01}(h)=p_1$$, and $$(q_{01}(t),p_{01}(t))$$ satisfy Hamilton’s equations in the time interval $$(0,h)$$. This is related to the unperturbed harmonic oscillator boundary-value problem given by,   q˙(t)=p(t),p˙(t)=−q(t),q(0)=q0, p(h)=p1. (4.4) This boundary-value problem is not well posed for values of $$h$$ that are odd multiples of $$\frac{\pi}{2}$$ and there are infinitely many solutions for such values of $$h$$. The integrator obtained from the exact discrete right Hamiltonian for the harmonic oscillator is given by   p1 =p0cos⁡(h)−q0sin⁡(h),q1 =p1tan⁡(h)+q0sec⁡(h). This integrator is analytically the true solution to the harmonic oscillator initial-value problem, where $$(q(0),p(0))=(q_0,p_0)$$ and the local truncation error will be zero. Noting that the method involves $$\tan(h)$$ and $$\sec(h)$$, we expect increased round-off error around odd multiples of $$\frac{\pi}{2}$$. Both of the integrators given by the exact discrete Lagrangian and the exact discrete right Hamiltonian have been implemented for the harmonic oscillator with initial conditions $$(q_0,p_0) = (1,0)$$ over the time interval $$[0,10000]$$, and the energy error is shown in Fig. 1. Note the jump in round-off error corresponding to values of $$h$$ that are odd multiples of $$\pi$$ (for the discrete Lagrangian) and odd multiples of $$\frac{\pi}{2}$$ (for the discrete right Hamiltonian). The bottom plot takes the minimum error of the two methods, and this indicates that a step size causing noticeable round-off error for one method will work just fine for the other method. Fig. 1. View largeDownload slide The first plot is the energy error vs. step size for the exact discrete right Hamiltonian applied to the harmonic oscillator. The second plot shows the energy error vs. step size for the exact discrete Lagrangian, while the third plot takes the minimum of the energy error from either method. Fig. 1. View largeDownload slide The first plot is the energy error vs. step size for the exact discrete right Hamiltonian applied to the harmonic oscillator. The second plot shows the energy error vs. step size for the exact discrete Lagrangian, while the third plot takes the minimum of the energy error from either method. In this particular case, we can conclude that the numerical difference between the symplectic maps generated by the respective exact discrete Lagrangian and exact discrete right Hamiltonian is a matter of numerical conditioning, which is inherited from the underlying ill-posedness of the associated boundary-value problem. Despite the fact that the methods are applied to an initial-value problem, numerical properties can be attributed to a boundary-value problem that is no longer visible in the methods themselves. Considering that many symplectic integrators are derived independently of the variational integrator formulation, perhaps some of their numerical properties can be better understood by reinterpreting them in the framework of variational integrators. 4.2 Averaged variational integrators for a nonlinearly perturbed harmonic oscillator Now we consider the previous averaging methods applied to a Hamiltonian of the form   H(q,p)=12(p2+q2)+ϵ3q3, (4.5) which is the Hamiltonian for a nonlinearly perturbed harmonic oscillator. The corresponding averaged Lagrangian is given by   Ld(q0,q1,h)=∫0h12(q˙A(t)2−qA(t)2)dt−∫0hϵ3qA(t)3dt, (4.6) where $$(q_A(t),\dot{q}_A(t))$$ is the solution corresponding to the Lagrangian $$L^{(A)}(q,\dot{q}) = \frac{1}{2}(\dot{q}^2 - q^2)$$ with boundary conditions $$(q_0,q_1)$$. Analogously, the averaged Hamiltonian is given by   Hd+(q0,p1,h)=p1qA(h)−∫0h12(pA(t)2−qA(t)2)dt+ϵ3∫0hqA(t)3dt, (4.7) where $$(q_A(t),p_A(t))$$ is the solution corresponding to the Hamiltonian $$H^{(A)}(q,p)$$ with boundary conditions $$(q_0,p_1)$$. Applying the discrete right and left Legendre transforms implicitly defines the discrete Hamiltonian map for $$L_d(q_0,q_1,h)$$ and the discrete right Hamiltonian map for $$H_d^+(q_0,p_1,h)$$, which yields the respective one-step methods. Numerical simulations were run over a time span from 0 to 10000 or the nearest integer value to 10000 for the respective time step. The initial conditions are given by $$(q_0,p_0) = (1,0)$$. Figures 2 and 3 show plots of the energy error vs. step size for two different values of $$\epsilon$$. The third plot in each of the figures hints that the discrete Lagrangian and discrete right Hamiltonian have numerical resonance that is nearly dual, in some sense, with respect to step size. The discrete Lagrangian exhibits excessive numerical resonance for step sizes near odd multiples of $$\pi$$, while the discrete right Hamiltonian exhibits excessive numerical resonance for step sizes near odd multiples of $$\frac{\pi}{2}$$. It should be noted that the arbitrary value of $$10^6$$ was substituted for output that was either near infinite or NaN. What is particularly striking is that the occurrence of the numerical resonance is intimately connected to the corresponding boundary values for each generating function. Fig. 2. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.1$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance around odd integer multiples of $$\frac{\pi}{2}$$ and exactly at odd multiples of $$\pi$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 2. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.1$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance around odd integer multiples of $$\frac{\pi}{2}$$ and exactly at odd multiples of $$\pi$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 3. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.001$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance at some odd integer multiples of $$\frac{\pi}{2}$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Fig. 3. View largeDownload slide Three plots of step size vs. energy error with fixed $$\epsilon = 0.001$$. The first plot corresponds to the averaged Hamiltonian, and it suffers from numerical resonance at some odd integer multiples of $$\frac{\pi}{2}$$. The second plot corresponds to the averaged Lagrangian which suffers from numerical resonance around odd multiples of $$\pi$$. The last plot takes the minimum error of the respective methods. Now this by no means provides a rigorous analysis of the numerical resonances, nor does it fully explain all of the resonance effects, but it does provide motivation and insight into the numerical differences between the discrete Lagrangian and discrete right Hamiltonian. A more in-depth analysis might be provided by applying something similar to modulated Fourier expansions (see Hairer and Lubich, 2001, 2012; Hairer et al., 2006, Chapter XIII). Modulated Fourier expansions are particularly well suited for oscillatory problems when large step sizes are sought. The standard backward error analysis relies on $$h\omega \to 0$$, which is not the case for high oscillatory problems when seeking large step sizes. Modulated Fourier expansions can provide a tool for deriving many of the same results as backward error analysis, such as long-term energy preservation. Furthermore, it can be quite useful for examining the step sizes that lead to excessive numerical resonance. However, it should be noted that while modulated Fourier expansions have been used quite successfully to analyse explicit trigonometric integrators, it is not quite as clear how easily it can deal with implicit integrators such as those obtained from the discrete averaged Lagrangian and discrete averaged Hamiltonian. 5. Conclusion Error analysis and symmetry results have now been extended to cover discrete Hamiltonian variational integrators. Furthermore, examples have been presented indicating that the underlying well-posedness in terms of the boundary conditions of the exact generating function can be directly related to numerical resonance. In conclusion, it is clear that the numerical properties of variational integrators are dependent on both the approximation scheme used in constructing the generating function and the type of generating function being approximated. This article indicates that the class of variational integrators generated using the Hamiltonian formulation are not necessarily equivalent to the ones obtained from the Lagrangian formulation, and it would therefore be of interest to continue developing methods based on the discrete Hamiltonian variational integrator formulation. In particular, the results presented suggest that further work remains to be done to better understand the circumstances under which it is preferable to favor one approach over the other. Funding National Science Foundation grants DMS-1010687, CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, DMS-1345013. References De León, M., Martín de Diego, D. & Santamaría-Merino, A. ( 2007) Discrete variational integrators and optimal control theory. Adv. Comput. Math. , 26, 251– 268. Google Scholar CrossRef Search ADS   Farr, W. M. ( 2009) Variational integrators for almost-integrable systems. Celest. Mech. Dyn. Astron. , 102, 105– 118. Google Scholar CrossRef Search ADS   Hairer, E. & Lubich, C. ( 2001) Long-time energy conservation of numerical methods for oscillatory differential equations. SIAM J. Numer. Anal. , 38, 414– 441. Google Scholar CrossRef Search ADS   Hairer, E. & Lubich, C. ( 2012) Modulated Fourier expansions for continuous and discrete oscillatory systems. Foundations of Computational Mathematics, Budapest 2011  ( Cucker, F. Krick, T. Pinkus A. & Szanto A. eds). London Mathematical Society Lecture Note Series. Cambridge: Cambridge University Press, pp. 113– 128. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric Numerical Integration,  2nd edn. Springer Series in Computational Mathematics, vol. 31. Berlin: Springer. Lall, S. & West, M. ( 2006) Discrete variational Hamiltonian mechanics. J. Phys. A , 39, 5509– 5519. Google Scholar CrossRef Search ADS   Leok, M. & Shingel, T. ( 2012) General techniques for constructing variational integrators. Front. Math. China , 7, 273– 303. Google Scholar CrossRef Search ADS   Leok, M. & Zhang, J. ( 2011) Discrete Hamiltonian variational integrators. IMA J. Numer. Anal. , 31, 1497– 1532. Google Scholar CrossRef Search ADS   Lew, A., Marsden, J. E., Ortiz, M. & West, M. ( 2003) Asynchronous variational integrators. Arch. Ration. Mech. Anal. , 167, 85– 146. Google Scholar CrossRef Search ADS   Leyendecker, S., Marsden, J. E. & Ortiz, M. ( 2008) Variational integrators for constrained mechanical systems. Z. Angew. Math. Mech. , 88, 677– 708. Google Scholar CrossRef Search ADS   Marsden, J. E., Patrick, G. W. & Shkoller, S. ( 1998) Multisymplectic geometry, variational integrators, and nonlinear PDEs. Comm. Math. Phys. , 199, 351– 395. Google Scholar CrossRef Search ADS   Marsden, J. E. & West, M. ( 2001) Discrete mechanics and variational integrators. Acta Numer. , 10, 357– 514. Google Scholar CrossRef Search ADS   Verhulst, F. ( 2000) Non-linear Differential Equations and Dynamical Systems . Berlin: Springer. Wisdom, J. & Holman, M. ( 1991) Symplectic maps for the $$n$$-body problem. Astron. J. , 102, 1528– 1538. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations