Lagrangian and Hamiltonian Taylor Variational Integrators

TL;DR

This paper introduces a Taylor-based variational integrator that leverages the structure of Taylor's method for higher-order accuracy and efficiency in solving Euler-Lagrange boundary-value problems, with demonstrated numerical benefits.

Contribution

It presents a novel Taylor variational integrator that achieves higher order accuracy and efficiency by exploiting Taylor's method structure, including a symmetric version.

Findings

Higher order accuracy compared to existing shooting methods

Enhanced efficiency through polynomial evaluations

Numerical experiments confirm effectiveness

Abstract

In this paper, we present a variational integrator that is based on an approximation of the Euler--Lagrange boundary-value problem via Taylor's method. This can viewed as a special case of the shooting-based variational integrator. The Taylor variational integrator exploits the structure of the Taylor method, which results in a shooting method that is one order higher compared to other shooting methods based on a one-step method of the same order. In addition, this method can generate quadrature nodal evaluations at the cost of a polynomial evaluation, which may increase its efficiency relative to other shooting-based variational integrators. A symmetric version of the method is proposed, and numerical experiments are conducted to exhibit the efficacy and efficiency of the method.

Equations

$$\delta {\mathbb{S}}_d=0,$$

$${\mathbb{S}}_d(q_0,q_1,\dots ,q_n)={\sum }_{i=0}^{n-1}{L}_d(q_i,q_{i+1}).$$

$$L_d^E(q_0,q_1;h)={\int }_0^{h}L(q_{01}(t),{\dot{q}}_{01}(t))dt,$$

$$D_2{L}_d(q_{k-1},q_k)+D_1{L}_d(q_k,q_{k+1})=0,$$

$$p_k=-D_1{L}_d(q_k,q_{k+1}),p_{k+1}=D_2{L}_d(q_k,q_{k+1}),$$

$$L_d(q_0,q_1;h)=L_d^E(q_0,q_1;h)+\mathcal{O}(h^{p+1}),$$

$$L_d^E(q_0,q_1;h)={\int }_0^{h}L({\Phi }_t(q_0,v_0))dt,$$

$$q_1={\pi }_Q\circ {\Psi }_h^{(r+1)}(q_0,{\tilde{v}}_0),$$

$$(q_{c_i},v_{c_i})$$

$$L_d(q_0,q_1;h)=h\sum _{i=1}^m{b}_{i}L(q_{c_i},v_{c_i}).$$

$$p_0$$

$$p_1$$

$$q_{01}(c_{i}h)$$

$$v_{01}(c_{i}h)$$

$$L_d^E(q_0,q_1;h)$$

$$=\left[h{\sum }_{i=1}^m{b}_{i}L(q_{01}(c_{i}h),v_{01}(c_{i}h))\right]+\mathcal{O}(h^{s+1})$$

$$=\left[h{\sum }_{i=1}^m{b}_{i}L(q_d(c_{i}h)+\mathcal{O}(h^{r+1}),v_d(c_{i}h)+\mathcal{O}(h^{r+1}))\right]+\mathcal{O}(h^{s+1})$$

$$=\left[h{\sum }_{i=1}^m{b}_{i}L(q_d(c_{i}h),v_d(c_{i}h))\right]+\mathcal{O}(h^{r+2})+\mathcal{O}(h^{s+1})$$

$$=L_d(q_0,q_1;h)+\mathcal{O}(h^{r+2})+\mathcal{O}(h^{s+1})$$

$$=L_d(q_0,q_1;h)+\mathcal{O}(h^{\min (r+1,s)+1}),$$

$$q_1=q_0+h{\tilde{v}}_0.$$

$$L_d(q_0,q_1;h)$$

$$\displaystyle=h\left[\frac{1}{2}\Big{(}\frac{q_{1}-q_{0}}{h}\Big{)}^{T}M\Big{(}\frac{q_{1}-q_{0}}{h}\Big{)}-V(q_{0})\right].$$

$$p_0$$

$$p_1$$

$$q_1$$

$$p_1$$

$${\tilde{v}}_0=\frac{q_1-q_0}{h}-\frac{h}{2}{M}^{-1}\nabla V(q_0).$$

$$q_1$$

$$p_1$$

$$H_d^{+,E}(q_0,p_1;h)$$

$$p_1={\pi }_{T^{\ast }Q}\circ {\Psi }_h^{(r)}(q_0,{\tilde{p}}_0),$$

$$(q_{c_i},p_{c_i})={\Psi }_{c_{i}h}^{(r)}(q_0,{\tilde{p}}_0),$$

$${\tilde{q}}_1={\pi }_Q\circ {\Psi }_h^{(r+1)}(q_0,{\tilde{p}}_0).$$

$$H_{d}^{+}(q_{0},p_{1};h)=p_{1}^{T}\tilde{q}_{1}-h\sum^{m}_{i=1}\Big{[}p_{c_{i}}^{T}\dot{q}_{c_{i}}-H\Big{(}\Psi^{(r)}_{c_{i}h}(q_{0},\tilde{p}_{0})\Big{)}\Big{]},$$

$$q_1=D_2{H}_d^{+}(q_0,p_1),p_0=D_1{H}_d^{+}(q_0,p_1).$$

Full text

1. Introduction

This paper is concerned with the systematic construction and analysis of Lagrangian and Hamiltonian variational integrators of arbitrarily high-order derived from an underlying Taylor integrator. This can be viewed, on the Lagrangian side, as a special case of the shooting-based variational integrators introduced in [11], which provided a general framework for constructing a Lagrangian variational integrator from a given one-step method.

The main limitation of the shooting-based variational integrator approach is that in order to achieve higher-order accuracy, one requires multiple steps of the underlying one-step method in order to obtain approximations of the solution of the Euler–Lagrange boundary-value problem at the quadrature points. This is of course the best one can hope to achieve given a generic one-step method, but for one-step methods such as collocation methods or Taylor methods, one obtains a continuous approximation that can be evaluated at multiple points. As such, these methods only require a single step of the one-step method in order to obtain a continuous approximation of the Euler–Lagrange boundary-value problem that can be used to construct discrete Lagrangians and discrete Hamiltonians that generate symplectic integrators.

We focus on the use of Taylor integrators as the underlying one-step method, since they can be efficiently implemented to arbitrarily high-order for a broad range of problems by leveraging automatic differentiation techniques, and the resulting solution can be evaluated at additional quadrature points at the cost of a polynomial evaluation.

2. Discrete Mechanics

Discrete Lagrangian mechanics [13] is based on a discrete analogue of Hamilton’s principle, referred to as the discrete Hamilton’s principle,

[TABLE]

where the discrete action sum, $\mathbb{S}_{d}:Q^{n+1}\rightarrow\mathbb{R}$, is given by

[TABLE]

The discrete Lagrangian, $L_{d}:Q\times Q\rightarrow\mathbb{R}$, is a generating function of the symplectic flow, and is an approximation to the exact discrete Lagrangian,

[TABLE]

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)$.

The discrete variational principle yields the discrete Euler–Lagrange (DEL) equation,

[TABLE]

which implicitly defines the discrete Lagrangian map $F_{L_{d}}:(q_{k-1},q_{k})\mapsto(q_{k},q_{k+1})$ for initial conditions $(q_{0},q_{1})$ that are sufficiently close to the diagonal of $Q\times Q$. This is equivalent to the implicit discrete Euler–Lagrange (IDEL) equations,

[TABLE]

which implicitly defines the discrete Hamiltonian map $\tilde{F}_{L_{d}}:(q_{k},p_{k})\mapsto(q_{k+1},p_{k+1})$, where the discrete Lagrangian is the Type I generating function of the symplectic transformation. Furthermore, the discrete Hamiltonian map associated with the exact discrete Lagrangian $\tilde{F}_{L_{d}^{E}}$ is the time-$h$ flow map of the Hamiltonian vector field. These observations serve as the basis by which the variational error analysis result of §2.1 is proven in [13]. In particular, variational error analysis relates the order to which a computable discrete Lagrangian approximates the exact discrete Lagrangian with the order of accuracy of the discrete Hamiltonian map when viewed as a one-step method for approximating the flow of Hamilton’s equations.

2.1. Variational error analysis

The natural setting for analyzing the order of accuracy of a variational integrator is the variational error analysis framework introduced in [13]. In particular, Theorem 2.3.1 of [13] states that if a discrete Lagrangian, $L_{d}:Q\times Q\rightarrow\mathbb{R}$, approximates the exact discrete Lagrangian, $L_{d}^{E}:Q\times Q\rightarrow\mathbb{R}$, given in (1) to order $p$, i.e.,

[TABLE]

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 $p$ accurate.

3. Lagrangian Taylor Variational Integrator

The exact discrete Lagrangian (1) is given by the action integral evaluated along the solution of the Euler–Lagrange boundary-value problem. In turn, the boundary-value problem with boundary data $(q_{0},q_{1})$ can be related to an initial-value problem with initial data $(q_{0},v_{0})$, which satisfies the condition $q_{1}=\pi_{Q}\Phi_{h}(q_{0},v_{0})$, where $\pi_{Q}:TQ\rightarrow Q$ is the canonical projection onto $Q$ and $\Phi_{h}:TQ\rightarrow TQ$ is the exact time-$h$ flow map. This yields the following characterization of the exact discrete Lagrangian,

[TABLE]

where $q_{1}=\pi_{Q}\Phi_{h}(q_{0},v_{0})$. The Taylor variational integrator is generated by a computable discrete Lagrangian obtained when the integral is approximated by a quadrature rule, and the Taylor method is used to approximate the flow map that relates the boundary data $(q_{0},q_{1})$ with the initial-value data $(q_{0},v_{0})$, and the trajectory associated with the initial data. The following summarizes the construction of the Taylor variational integrator.

1. (i)

The approximation to $\dot{q}(0)=v_{0}$, denoted as $\tilde{v}_{0}$, is defined via the inverse problem,

[TABLE]

where $\pi_{Q}:TQ\rightarrow Q$ is the canonical projection onto $Q$ and $\Psi^{(r+1)}_{h}:TQ\rightarrow TQ$ denotes a $(r+1)$-order Taylor method. 2. (ii)

Generate approximations to the quadrature nodal values, $q_{c_{i}}\approx q(c_{i}h)$ (excluding $q_{1}$ if needed, which is assumed to be given) and $v_{c_{i}}\approx\dot{q}(c_{i}h)$, via Taylor’s method using $\tilde{v}_{0}$,

[TABLE] 3. (iii)

Apply the quadrature rule to construct the associated discrete Lagrangian,

[TABLE] 4. (iv)

Applying the discrete Legendre transforms implicitly defines the method,

[TABLE]

Remark**.**

It may seem like a waste to solve for $\tilde{v}_{0}$ using a $(r+1)$-order Taylor method, and then to use only a $r$-order method to solve the Euler–Lagrange boundary-value problem, but from an implementation perspective, no additional derivative evaluations are needed to solve (5), other than those already required in implementing the $r$-order Taylor method on $TQ$. In fact, it is an efficient use of the higher-derivative information we already needed to compute in order to construct the $r$-order Taylor method on $TQ$.

This apparent discrepancy can be resolved by thinking of equation (5) as being a $(r+1)$-order Taylor method for the second-order differential equation on $Q$, and (6) as a $r$-order Taylor method on the first-order differential equation on $TQ$. In particular, notice that because of the canonical projection $\pi_{Q}$ in equation (5), we only need to compute up to $q^{(r+1)}(0)$ in order to solve for $\tilde{v}_{0}$, instead of the up to $q^{(r+2)}(0)$ that is necessary to define $\Psi^{(r+1)}_{h}$. But, we needed to compute up to $v^{(r)}(0)=q^{(r+1)}(0)$ in order to construct $\Psi^{(r)}_{h}$, the $r$-order Taylor method on $TQ$.

The following lemmas are needed for a theorem on the accuracy of the method. These lemmas can be proved using Lipschitz continuity and triangle inequalities (see Appendices for their proofs).

Lemma 1**.**

$\tilde{v}_{0}$* as defined by, (5), approximates $v_{0}$ to at least $\mathcal{O}(h^{r+1})$.*

Lemma 2**.**

A $r$-order Taylor method with initial conditions $(q_{0},\tilde{v}_{0})$, where $\tilde{v}_{0}$ is defined by (5), is accurate to at least $\mathcal{O}(h^{r+1})$ for the Euler–Lagrange boundary-value problem with boundary conditions $(q_{0},q_{1})$.

Theorem 1**.**

Assuming a Lagrangian $L$ that is Lipschitz continuous in both variables, then for a $r$-order accurate Taylor method, $\Psi^{(r)}_{h}$, and a $s$-order accurate quadrature formula, the associated Taylor discrete Lagrangian (7) has order of accuracy at least $\min(r+1,s)$.

Proof.

$(q_{d}(t),v_{d}(t))$, associated with the Taylor method $\Psi_{h}$ of order $r$ and initial data $(q_{0},\tilde{v}_{0})$, approximates the exact solution $(q_{01}(t),v_{01}(t))$ of the Euler–Lagrange boundary-value problem with the following error,

[TABLE]

If the numerical quadrature formula is order $s$ accurate, then

[TABLE]

where we used the quadrature approximation error, the error estimates on the shooting solution, and the assumption that $L$ is Lipschitz continuous in both variables. ∎

The choice of the Taylor method as the underlying one-step method has the advantage that it only requires one to precompute the prolongation of the Euler–Lagrange vector field once at the initial time, and the computational cost is not increased appreciably by having to compute the numerical solution at multiple quadrature nodes, since that only requires a polynomial evaluation. This efficiency in evaluation improves upon the methods outlined in [10] and [11], which utilized collocation and the shooting-method, respectively.

Example 1**.**

Consider a first-order Taylor variational integrator that uses the rectangular quadrature rule about the initial point. We assume a Lagrangian of the form $L(q,\dot{q})=\frac{1}{2}\dot{q}^{T}M\dot{q}-V(q)$. Then the integrator is constructed as follows:

1. (i)

The inverse problem is,

[TABLE]

This implies $\tilde{v}_{0}=\frac{q_{1}-q_{0}}{h}$, where $q_{0}$, $q_{1}$ are the given boundary conditions. 2. (ii)

The quadrature nodal values are $q_{c_{1}}=q_{0}$ and $v_{c_{1}}=\tilde{v}_{0}=\frac{q_{1}-q_{0}}{h}$. 3. (iii)

The corresponding discrete Lagrangian is given by,

[TABLE] 4. (iv)

The discrete Legendre transforms are given by,

[TABLE]

With some rearranging and substitution we see that this is symplectic Euler-A,

[TABLE]

If we use the rectangular quadrature rule about the end point, then the resulting method would be symplectic Euler-B. If instead we choose the trapezoid quadrature rule, then the resulting method will be Störmer–Verlet. All three of these classic symplectic integrators can be derived as Taylor variational integrators. However, there are also novel methods that come from the Taylor variational integrator framework, as the next example illustrates.

Example 2**.**

Consider a second-order Taylor variational integrator, which utilizes a first-order Taylor method combined with the trapezoid rule to approximate the discrete Lagrangian. The approximate initial velocity is given by,

[TABLE]

The resulting method is an explicit second-order method given by,

[TABLE]

As demonstrated above, Lagrangian Taylor variational integrators provide a very general family of symplectic integrators that include not only classic symplectic integrators, but also novel symplectic integrators. The Taylor variational integrator is amenable to the construction of higher-order symplectic integrators that can benefit from many of the numerical techniques that have enhanced the classical Taylor method (see [8], [19]). In particular, automatic differentiation allows for accurate and relatively cheap derivative evaluations (see [5], [16], [14]). In general, higher-order Taylor variational integrators will require solving a system of nonlinear equations, which can be dealt with using standard methods (see [6]). While it is clear that Taylor variational integrators will have a higher computational cost than the Taylor method, in many cases the Taylor variational integrator can preserve accuracy and structure for larger step sizes, which may justify the higher cost per step. We will further examine these topics in section 5. Next, we consider discrete Hamiltonian formulations and symmetric formulations of the Taylor variational integrator.

4. Hamiltonian and Symmetric Taylor Variational Integrators

4.1. Hamiltonian Taylor Variational Integrators

Thus far, we have derived the Taylor variational integrator by approximating the discrete Lagrangian, which is a type I generating function of the symplectic map/integrator. However, we will also consider the discrete right and discrete left Hamiltonians (see [9], [12]), which are type II and type III generating functions, respectively. The motivation being that for a degenerate Hamiltonian there may be no corresponding Lagrangian formulation, in which case the discrete Hamiltonian formulation may be the only way to construct a variational integrator. Also, it has recently been shown in [17] that even when the Legendre transform is a diffeomorphism, the discrete Lagrangian and discrete Hamiltonian formulation generated by a fixed approximation scheme can lead to different variational integrators.

The boundary-value formulation of the exact discrete right Hamiltonian is given by,

[TABLE]

where $(q(t),p(t))$ satisfy Hamilton’s equations with boundary conditions $q(0)=q_{0}$, $p(T)=p_{1}$. Now let us consider the construction of a Taylor discrete right Hamiltonian.

1. (i)

Construct a $r$-order Taylor expansion on the cotangent bundle, $T^{*}Q$, and solve for $\tilde{p}_{0}$,

[TABLE]

where $\pi_{T^{*}Q}:(q,p)\mapsto p$. 2. (ii)

Pick a quadrature rule of order $s$ with quadrature weights and nodes given by $(b_{i},c_{i})$ for $i=1,\ldots,m$. 3. (iii)

Use a $r$-order Taylor method to generate approximations of $(q(t),p(t))$ at the quadrature nodes,

[TABLE]

and use a $(r+1)$-order Taylor method on the configuration manifold to generate the approximation to the boundary term $q_{1}$,

[TABLE] 4. (iv)

Use the quadrature rule and approximate boundary term, $\tilde{q}_{1}$, to construct the discrete right Hamiltonian of order $\min(r+1,s)$,

[TABLE]

where $\dot{q}_{c_{i}}$ is obtained by inverting the continuous Legendre transform, $(q_{c_{i}},p_{c_{i}})=\mathbb{F}L(q_{c_{i}},\dot{q}_{c_{i}})$. 5. (v)

The method is implicitly defined by the implicit discrete right Hamilton’s equations,

[TABLE]

The boundary-value formulation of the exact discrete left Hamiltonian is given by,

[TABLE]

where $(q(t),p(t))$ satisfy Hamilton’s equations with boundary conditions $q(T)=q_{1}$, $p(0)=p_{0}$. Now let us consider the construction of a Taylor discrete left Hamiltonian.

1. (i)

Construct a $(r+1)$-order Taylor expansion on the cotangent bundle, $T^{*}Q$, and solve for $\tilde{q}_{0}$,

[TABLE] 2. (ii)

Pick a quadrature rule of order $s$ with quadrature weights and nodes given by $(b_{i},c_{i})$ for $i=1,\ldots,m$. 3. (iii)

Use a $r$-order Taylor method to generate approximations of $(q(t),p(t))$ at the quadrature nodes,

[TABLE] 4. (iv)

Use the quadrature rule and approximate boundary term, $\tilde{q}_{0}$, to construct the discrete left Hamiltonian of order $\min(r+1,s)$,

[TABLE]

where $\dot{q}_{c_{i}}$ is obtained by inverting the continuous Legendre transform, $(q_{c_{i}},p_{c_{i}})=\mathbb{F}L(q_{c_{i}},\dot{q}_{c_{i}})$. 5. (v)

The method is implicitly defined by the implicit discrete left Hamilton’s equations,

[TABLE]

The Störmer–Verlet method can be derived as a Lagrangian Taylor variational integrator by choosing $r=0$ for the respective Taylor methods and using the trapezoid rule for the quadrature rule. This yields a discrete Lagrangian corresponding to the Störmer–Verlet method,

[TABLE]

Choosing $r=0$ and the trapezoid rule to construct a Hamiltonian Taylor variational integrator results in a discrete right Hamiltonian given by,

[TABLE]

and a discrete left Hamiltonian given by,

[TABLE]

The corresponding methods are not Störmer–Verlet, in fact they are neither symmetric nor explicit. However, a simple calculation shows that these discrete Hamiltonians are adjoint to each other (see [17] for info on adjoint discrete Hamiltonians), i.e. $-H_{d}^{+}(q_{1},p_{0};-h)=H_{d}^{-}(q_{1},p_{0};h)$. Therefore, a symmetric method can be constructed by composing the two methods. We will denote the resulting symmetric method by SVHd, and we compare it to Störmer–Verlet in section 5 (see Figure 7 and Figure 9).

It should be noted that some approximations schemes do yield the same method when applied to a discrete Lagrangian and a discrete right/left Hamiltonian. For instance, choosing $r=0$ and the rectangular rule about the end point will yield symplectic Euler-B for both the discrete Lagrangian and discrete right Hamiltonian approximation. When can we expect a fixed approximation scheme applied to a discrete Lagrangian and a discrete right Hamiltonian to yield the same method? The following theorem answers this question.

Theorem 2**.**

Assuming a regular Lagrangian, we consider a fixed approximation scheme used to construct a discrete Lagrangian, $L_{d}$, and a discrete right Hamiltonian, $H_{d}^{+}$. This results in two integrators, $\tilde{F}_{L_{d}}:(q_{0},p_{0})\mapsto(q_{1,L_{d}},p_{1,L_{d}})$ and $\tilde{F}_{H_{d}^{+}}:(q_{0},p_{0})\mapsto(q_{1,H_{d}^{+}},p_{1,H_{d}^{+}})$. If the discrete right Hamiltonian approximation satisfies $p_{1,H_{d}^{+}}=D_{2}L_{d}(q_{0},\hat{q}_{1})$, where $\hat{q}_{1}$ is the approximated value of $q_{1}$, then the integrators represent the same map, i.e., $(q_{1,L_{d}},p_{1,L_{d}})=(q_{1,H_{d}^{+}},p_{1,H_{d}^{+}})$.

We have placed the proof of the above theorem in the appendix. It is important to note that even though the theorem guarantees the analytical equivalence of the integrators, this does not guarantee numerical equivalence (see [17]).

4.2. Symmetric Lagrangian Taylor Variational Integrators

Consider the following variational derivation of the Störmer–Verlet method. Construct the discrete Lagrangian by using the trapezoid rule and approximating $\tilde{v}_{0}$ and $\tilde{v}_{1}$ with the inverse problems given by,

[TABLE]

where $\Psi^{r}_{h}$ denotes the $r$-th order Taylor method with step size $h$. Then, the velocity approximations are given by,

[TABLE]

and the resulting discrete Lagrangian yields the Störmer–Verlet method,

[TABLE]

It is well-known that Störmer–Verlet is a symmetric method, and that symmetric methods preserve important structure of time-reversible equations and are desirable for highly-oscillatory problems (see chapters V and XI of [4]). We can generalize the above approximation to yield a class of symmetric Taylor variational integrators. The approximation scheme uses a symmetric quadrature rule with weights and nodes $\{b_{i},c_{i}\}_{i=1}^{m}$, and the Taylor method, and it is outlined as follows:

1. (i)

Solve the inverse problems for $\tilde{v}_{0}$ and $\tilde{v}_{1}$,

[TABLE] 2. (ii)

Generate approximations to the quadrature nodes $(q_{c_{i}},v_{c_{i}})$ via,

[TABLE]

Note $q_{0}$, $q_{1}$, $\tilde{v}_{0}$, and $\tilde{v}_{1}$ are used as the approximations for their respective quadrature nodal values. Also, since the quadrature rule is assumed to be symmetric, $c_{i}=1-c_{m-i+1}$ and $b_{i}=b_{m-i+1}$. 3. (iii)

Construct the discrete Lagrangian,

[TABLE] 4. (iv)

Apply the discrete Legendre transforms to implicitly define the variational integrator,

[TABLE]

Theorem 3**.**

The symmetric Taylor variational integrator is a symmetric method.

Proof.

By theorem 2.4.1 of [13], it is sufficient and necessary to show that the discrete Lagrangian of the symmetric Taylor variational integrator is self-adjoint, i.e., $L_{d}(q_{0},q_{1};h)=-L_{d}(q_{1},q_{0};-h)$. We will use (*) to denote the approximated values resulting from exchanging $(q_{0},q_{1},h)$ for $(q_{1},q_{0},-h)$. Exchanging $(q_{0},q_{1},h)$ for $(q_{1},q_{0},-h)$ transforms (10) into,

[TABLE]

so that $\tilde{v}^{*}_{0}=\tilde{v}_{1}$ and $\tilde{v}^{*}_{1}=\tilde{v}_{0}$. Therefore,

[TABLE]

The second to last line follows from the fact that the quadrature rule is symmetric and therefore satisfies $1-c_{i}=c_{m-i+1}$. The same steps show that $v^{*}_{c_{i}}=v_{c_{m-i+1}}$. The symmetric quadrature rule also implies that $b_{i}=b_{m-i+1}$, so that we have the following,

[TABLE]

∎

Theorem 4**.**

Given a regular Lagrangian, an odd $r$-order Taylor method, and a symmetric quadrature rule of order $r+1$, then the resulting symmetric Taylor variational integrator is of order $r+1$.

Proof.

Lemmas 1 and 2 imply that the nodal value approximations, $q_{c_{i}}$ and $v_{c_{i}}$, are of order $r$ and $r-1$ respectively. Therefore,

[TABLE]

We have used the order of the nodal approximations, the error order of the quadrature rule, and the Lipschitz continuity of a regular Lagrangian. By theorem 2.3.1 of [13], the resulting variational integrator, denoted by $\tilde{\Psi}_{h}$, is at least of order $r$, i.e.,

[TABLE]

where $\Phi_{h}$ is the true flow of the Euler–Lagrange equations, and the last equality is a consequence of the implicit function theorem.

Finally, since the variational integrator is symmetric and $r+1$ is even, the method will be of order $r+1$ as the following implies.

[TABLE]

which implies $C(q_{0},v_{0})h^{r+1}=0$, and the method is of order $r+1$ as claimed. ∎

The symmetric Taylor variational integrator is of order $r+1$, but only requires the derivatives of a $r$-order Taylor method, which makes it more efficient than the non-symmetric Taylor variational integrator, in addition to the qualitative benefits associated with its symmetry. However, applying this approximation scheme to generate a discrete Hamiltonian will not directly lead to a symmetric method. Recall that the symmetric Taylor variational integrator was inspired by Störmer–Verlet, so it is likely that using this approximation scheme to generate a discrete right and left Hamiltonian will result in the discrete left and right Hamiltonian methods that are adjoint to each other. In that case, the composition of these methods should yield a symmetric method from the discrete Hamiltonian formulation. We conjecture that if an approximation scheme yields a symmetric discrete Lagrangian, then the corresponding discrete right and left Hamiltonians will be adjoint. We will explore this further in future work.

5. Numerical Implementation and Experiments

We now discuss the numerical implementation of the methods introduced in this paper. Below, we present the algorithm for the Lagrangian Taylor variational integrator, and we discuss some of our observations about the implementation details. Additionally, we compare the methods to other kinds of variational integrators, and discuss their relative merits.

Algorithm Given $(q_{0},p_{0})$, $h$, $L(q(t),\dot{q}(t))$, the Euler–Lagrange vector field, quadrature weights and nodes $\{(b_{i},c_{i})\}_{i=1:m}$, and the desired order of the method $r+1$, then the Taylor variational integrator will output $(q_{1},p_{1})$ and is implemented as follows:

1. (1)

Prolongate the Euler–Lagrange vector field to obtain derivatives $q^{(j)}(q(t),v(t))$ for $j=1,\dots,r+1$. 2. (2)

Compute the partial derivatives $\frac{\partial q^{(j)}(q,v)}{\partial q}$ and $\frac{\partial q^{(j)}(q,v)}{\partial v}$. 3. (3)

Solve the following nonlinear system for $q_{1}$ and $\tilde{v}_{0}$:

[TABLE] 4. (4)

Finally, $p_{1}$ is given explicitly by,

[TABLE]

When solving the nonlinear system that arises above, the following points should be noted:

1. (1)

In general, the nonlinear system is not amenable to a fixed-point iteration, so a form of Newton’s method is preferable. 2. (2)

Each iteration will require evaluation of

[TABLE] 3. (3)

The following requires computing $\frac{\partial\tilde{v}_{0}}{\partial q_{0}}$,

[TABLE]

Fortunately, this can be found explicitly and need only be computed once at the beginning of the iteration,

[TABLE] 4. (4)

Likewise, when solving $p_{1}=\frac{\partial L_{d}(q_{0},q_{1})}{\partial q_{1}}$, it will be necessary to compute

[TABLE]

which is explicit and is composed of terms that have already been computed.

Observe that good initial guesses for the nonlinear system are provided with little computational cost, by using a $(r+1)$-order Taylor method for $q_{1}$ and the Legendre transform of $p_{0}$ for $\tilde{v}_{0}$. Since this yields an approximate solution that is comparable in accuracy to the one obtained by the corresponding Taylor variational integrator, this yields a predictor-corrector implementation, where the Taylor variational integrator applies a symplectic correction that converges very rapidly. In general, when solving a nonlinear system as part of a symplectic method, the method becomes an almost symplectic method (see [22]) unless it is solved to within machine precision. This implies that the error tolerance of the nonlinear solver will dictate to what order the symplectic structure is preserved and consequently, how well near-energy conservation is preserved (see Figure 1).

In practice, setting the nonlinear solver tolerance one or two orders above the order of the integrator is sufficient to maintain symplecticity. For most Taylor variational integrators, the nonlinear solver with moderate tolerance converges in a few iterations, and often in one or no iterations. The symmetric Taylor variational integrator showed excellent nonlinear convergence, and only required one iteration of the nonlinear solver for the various experiments we ran.

5.1. Automatic Differentiation

As with the Taylor method, an efficient general purpose implementation will require an efficient means of computing derivatives, such as automatic differentiation. For the following simulations, we used the AdiGator automatic differentiation package for MATLAB (see [15]). Implementation of a high-order Taylor variational integrator requires both the evaluation of higher time derivatives, $q^{(p+1)}(q_{0},\tilde{v}_{0})$, and the evaluation of the Jacobians of the time derivatives w.r.t. $q_{0}$ and $\tilde{v}_{0}$. The Jacobian evaluations are the most expensive part of the method (see Figure 2), especially for higher-dimensional systems, and for efficient high-order methods, the cost of Jacobian evaluations will need to be reduced to a level comparable to the time derivative. There appears to be some relationships between the Jacobians and the time derivatives that could potentially be exploited to decrease the evaluation costs. For instance,

[TABLE]

which allows us to replace expensive Jacobian evaluations with cheaper time derivative evaluations. Additionally, Jacobians of higher-order time derivatives appear to have some relations to Jacobians of lower-order time derivatives, such as,

[TABLE]

Hopefully, a good implementation of automatic differentiation will already take advantage of such relationships.

Automatic differentiation greatly benefits from the way it is compiled, which means the more efficient implementations will be in languages such as Fortran or C++. Another aspect to consider is parallel implementation. Combining automatic differentiation and parallel computing techniques has been shown to significantly reduce computational time (see [1]).

One possible implementation for the algorithm would be to construct the Taylor discrete Lagrangian, then apply automatic differentiation to the discrete Lagrangian in combination with a nonlinear solver to recover the discrete Legendre transforms and consequently $(q_{1},p_{1})$. In fact, this could provide a more general framework for the derivation of all implicit variational integrators.

5.2. Comparison of Methods

The simulations compare the discrete Lagrangian form of the Taylor variational integrator (TVI), the discrete right Hamiltonian form of the Taylor variational integrator (HTVI), the symmetric Taylor variational integrator of 4th order (SV4), Taylor’s method, and the Runge–Kutta shooting variational integrators (ShVI) (see [11]). Overall, high-order Taylor methods perform quite well in terms of computational time versus global error. However, as the length of integration time becomes very large, the variational integrators begin to show their strength. Of the three variational integrators, the symmetric Taylor variational integrator is the most efficient.

Comparison of the Lagrangian or Hamiltonian Taylor variational integrator to the Runge–Kutta shooting variational integrator does not result in a clear winner in terms of computational efficiency. It is well known that beyond 4th-order, Runge–Kutta (RK) methods require a higher number of stages/function evaluations, and the number of stages grows faster for vector differential equations as compared to scalar differential equations (see [2]). The number of order conditions grows quite quickly. For instance a 4th-order RK method has 8 order conditions, a 7th-order RK method has 85 order conditions, and a 25th-order method has 3,231,706,871 order conditions (see [20]). However, a 25th-order RK method only has 313 stages, so the function evaluations grow at a much slower rate. The Taylor method must contend with the increasing cost of evaluating higher-order derivatives, which for our implementation grows at a rate of $2^{n}$, where $n$ is the order of the derivative. For methods less than order 10 the difference in computational cost of the Taylor variational integrator and the Runge–Kutta based shooting variational integrator did not seem significant. However, the symmetric Taylor variational integrator did exhibit lower evaluation costs than the other methods. It should be noted that the most efficient implementations of the Taylor method involve variable stepsizes, and symplectic integrators are not predisposed to variable stepsizes.

The following simulations were implemented in MATLAB.

5.3. Simple Pendulum

Consider the simple pendulum with unit mass and length in a gravitational field with $g=-9.8\mathrm{m/s}^{2}$, where $q$ is parametrized by the angle between the y-axis and the pendulum. The corresponding Lagrangian is,

[TABLE]

The Euler–Lagrange equation yields,

[TABLE]

In Figure 3, the level sets of the corresponding Hamiltonian are compared to the trajectories generated by a 2nd-order Taylor variational integrator (TVI2) (see Example 2). The numerical solutions appear nearly identical to the level sets of the Hamiltonian, which indicates that the variational integrator exhibited good energy behavior for a variety of initial conditions.

The simulation in Figure 4 used initial conditions $(q_{0},p_{0})=(\frac{\pi}{2},0)$. The 6th-order Taylor variational integrator performed well at a stepsize of $h=0.5$, while the 6th-order Taylor method failed to generate a reasonable approximation for this stepsize. The ability of the Taylor variational integrator to perform well at larger stepsizes may gives it an advantage over traditional Taylor methods.

In Figure 5, we compare various types of Taylor variational integrators against the shooting-based variational integrator (ShVI). The plots compare the energy error versus computational time for methods of various order. It is clear the the symmetric Taylor variational integrator (SV4) is the most efficient in this respect, but it is not so clear whether the non-symmetric Taylor variational integrators (TVI and HTVI) are more efficient than ShVI.

5.4. Kepler’s Planar $2$-Body Problem

Consider two bodies interacting under mutual gravity and set one body as the center of the coordinate system (see [4]). Thus, constraining them to lie in a plane, we have Kepler’s planar 2-body problem with corresponding Lagrangian,

[TABLE]

Note here we are using $q_{1}$ and $q_{2}$ as the first and second components of $\mathbf{q}$. This in turn gives us the Euler–Lagrange equations,

[TABLE]

Our simulations used initial conditions $\mathbf{q}_{0}=\begin{bmatrix}1\\ 0\end{bmatrix}$ and $\mathbf{p}_{0}=\begin{bmatrix}0\\ 0.8\end{bmatrix}$. Figure 6 compares various Taylor variational integrators to Taylor methods of the same order using a stepsize of $h=0.25$. The trajectories of the Taylor methods for this stepsize behave poorly, while variational integrators show good qualitative performance.

Figure 7 compares the Störmer–Verlet method (SV) to the discrete Hamiltonian composition method (SVHd) discussed in section 4.1. Given that the Störmer–Verlet method is explicit, while SVHd is implicit, it is no surprise that the Störmer–Verlet method has lower computational cost. However, SVHd does exhibit lower energy error and performs slightly better qualitatively, so when the problem is non-separable (and SV is implicit), SVHd may be a better alternative.

5.5. Henon-Heiles Model

The Henon–Heiles model attempts to capture the dynamics of a galaxy with cylindrical symmetry (see [4] for more info). The Hamiltonian is given by, $H(p,q)=\frac{1}{2}(p_{1}^{2}+p_{2}^{2})+U(q)$, where $U(q)=\frac{1}{2}(q_{1}^{2}+q_{2}^{2})+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3}$. The corresponding Euler–Lagrange equation is,

[TABLE]

It is known that the dynamics become chaotic at higher energy levels. The following simulations were conducted with an initial energy level of $H_{0}=\frac{1}{12}$ (see Figure 8) and $H_{0}=\frac{1}{8}$ (see Figure 5). The second energy value corresponds to a chaotic system.

In Figure 8, we compare the 6th-order Taylor variational integrator (TVI6), the 6th-order Runge–Kutta shooting-based variational integrator (ShVI6), and the 4th-order symmetric Taylor variational integrator (SV4) applied to the Henon-Heiles model with $H_{0}=\frac{1}{12}$. For global errors between $10^{-1}$ and $10^{-5}$, SV4 is the more efficient method. Amongst the higher-order methods, TVI6 and ShVI6 appear to be the more efficient methods. A 6th-order symmetric Taylor variational integrator would be even more efficient for higher-order accuracy.

5.6. Fermi-Pasta-Ulam Model

The Fermi-Pasta-Ulam (FPU) model has a particularly distinguished place in the history of numerical simulations and nonlinear dynamics (see [3]). We apply the modified model as outlined in [4], consisting of a sequence of 6 mass points, fixed at both ends connected on opposite sides by a series of soft nonlinear springs and stiff linear springs. Letting $\{q_{i},p_{i}\}_{i=1}^{6}$ denote the displacements and velocities of the mass points, the corresponding Hamiltonian is given by,

[TABLE]

where $\omega=50$. By using the change of variables,

[TABLE]

the resulting Hamiltonian system has a nearly conserved quantity $I=I_{1}+\dots+I_{3}$, where

[TABLE]

is the energy of the $j$th stiff spring. Despite the significant energy exchange between individual springs, the total oscillatory energy, $I$, remains near constant. Our simulations used initial values of,

[TABLE]

Figure 9 compares the Störmer–Verlet method to SVHd. The first couple of plots use a stepsize of $h=0.03$, which is on the boundary of the linear stability of Störmer–Verlet (i.e. $h\omega=1.5$). SVHd does appear to be qualitatively more accurate, but neither method does well at this stepsize. For $h=0.01$, both methods give a much better qualitative representation of the system, but their global errors are still too large to be considered accurate. None of the methods in this paper are appropriate for a highly-oscillatory model such as the FPU model. For an accurate solution, one should consider either the IMEX method (see [21]) or Filon-type methods (see [7]). The combination of exponential type integrators with symplectic and energy-preserving integrators was also recently considered in [18].

5.7. Outer Solar System

Consider the motion of the five outer planets (including Pluto) relative to the sun. The corresponding Hamiltonian for this N-body problem is given by,

[TABLE]

where $G=2.95912208286\cdot 10^{-4}$. The initial data and masses is taken from Section 1.2.4 of [4], and corresponds to September 5, 1994 at 0h00. In Figure 10, we compare the 4th and 6th-order Taylor variational integrators to the 4th and 6th-order Taylor methods. The simulations was over the time period $[0,200000]$, and the stepsize was $h=400$ (days). The 4th-order methods did not produce a useful simulation at this stepsize, but both 6th-order integrators give a good representation of the system.

6. Conclusions and Future Directions

The Taylor variational integrators provide a way to build high-order symplectic integrators and include many of the classic symplectic integrators as special cases, i.e., symplectic Euler and Störmer–Verlet. This provides a framework for importing the large body of literature on the efficient construction of high-order Taylor integrators in order to construct similarly high-order symplectic integrators.

In particular, these methods can be viewed as a symplectic correction to higher-order Taylor methods that typically converges in a small number of iterations. By viewing these as predictor-corrector methods, one can interpolate between Taylor methods and Taylor variational integrators, and it would be interesting to see the extent to which a fixed number of iterations of the symplectic corrector can improve upon the performance of Taylor integrators for realistic problems.

The numerical simulations demonstrate that the geometric structure-preserving properties of symplectic integrators can be important for achieving numerical stability of long time simulations, so it should be of great interest to the computational astrophysics community to combine the high-order accuracy of high-order Taylor integrators with the geometric structure-preserving properties of variational integrators.

The most efficient implementations of the Taylor method utilize a variable stepsize, and extending variable stepsizes to the variational integrator framework is an area that deserves continued research. We are currently considering an approach based on the combination of Hamiltonian variational integrators and the Poincaré transformation that is quite promising. In particular, we note that the use of Hamiltonian as opposed to Lagrangian variational integrators is critical, as the Poincaré transformed Hamiltonian is degenerate, and there is no corresponding Lagrangian formulation.

Acknowledgements

The authors were supported in part by NSF grants DMS-1001521, CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, and NSF CAREER award DMS-1010687.

Appendix A Detailed Proofs

Given an Euler–Lagrange equation of the form,

[TABLE]

we denote the exact solution of the Euler–Lagrange boundary-value problem with boundary conditions $(q_{0},q_{1})$ by $(q(t),v(t))$. We seek an estimate of the true initial velocity, ${v}_{0}$, for the corresponding Euler–Lagrange initial-value problem, with order of accuracy $r$. Let us denote this estimate by $\tilde{v}_{0}$. Given a one-step method, $\hat{\Psi}_{h}:TQ\rightarrow TQ$, with order of accuracy $r+1$, we solve for the initial velocity $\tilde{v}_{0}$, such that,

[TABLE]

where $\pi_{Q}:TQ\rightarrow Q$ is the canonical projection. Let $\Phi_{h}:TQ\rightarrow TQ$ be the exact time-$h$ flow map of the Euler–Lagrange initial-value problem. By definition, the exact Euler–Lagrange flow applied to the initial condition $(q_{0},v_{0})$ is a solution of the Euler–Lagrange boundary-value problem with boundary conditions $(q_{0},q_{1})$, where

[TABLE]

Consider a Taylor method with order of accuracy $r$ and $r+1$,

[TABLE]

and

[TABLE]

where $q^{(k)}(0)$ is calculated by considering the prolongations of the Euler–Lagrange vector field, and evaluating it at $(q_{0},\tilde{v}_{0})$. An analogous approach, involving the prolongation of the Euler–Lagrange vector field at both the initial and final time, which can be viewed as a two-point Taylor method, was used to develop a prolongation-collocation variational integrator in [10].

Lemma 3**.**

$\tilde{v}_{0}$* as defined by, (11) and (14), approximates $v_{0}$ to at least $\mathcal{O}(h^{r+1})$.*

Proof.

Solving $\pi_{Q}\circ\hat{\Psi}_{h}(q_{0},\tilde{v}_{0})=q_{1}$ for $\tilde{v}_{0}$ yields,

[TABLE]

Since the exact solution $q(t)\in C^{r+2}([0,h])$, using Taylor’s Theorem, we have,

[TABLE]

Solving for $v_{0}$ yields,

[TABLE]

Now evaluating the norm of the difference we have,

[TABLE]

Since $q(t)\in C^{r+2}([0,h])$ each of $f^{(i-1)}$ is Lipschitz continuous in its arguments for $i=1,2,\ldots,r$. Let $M_{i}$ be the Lipschitz constant for $f^{(i-1)}$ over the compact interval $[0,C]$ with respect to velocity, and $C>0$ can be chosen so that $M_{i}$, $i=1,2,\ldots,p,$ is bounded. Using the triangle inequality, we have,

[TABLE]

Rearranging, we have,

[TABLE]

By continuity, there exists $\tilde{C}$ satisfying $0<\tilde{C}<C$, such that for all $h$ satisfying $0<h<\tilde{C}$, the term inside the parenthesis on the leftmost expression is positive and bounded away from zero. That concludes the proof.∎

Remark**.**

It is worth noting that a similar proof may be given for any $(r+1)$-order one-step method. This is due to the fact that any $(r+1)$-order one-step method agrees with the $(r+1)$-order Taylor’s method up to a local truncation error of order $\mathcal{O}(h^{r+2})$. Thus, the only change in the error term in the proof would be to replace $R_{r+1}(h)$ by the sum of the local truncation error of the one-step method and $R_{r+1}(h)$, which are both $\mathcal{O}(h^{r+2})$. Thus, this result can be generalized to any one-step method of the desired order.

Using this result, we can show that starting our $r$-order Taylor method at $\tilde{v_{0}}$, rather than at $v_{0}$, will not affect the order of accuracy of the method.

Lemma 4**.**

A $r$-order Taylor method, defined by (13), with initial conditions $(q_{0},\tilde{v}_{0})$, where $\tilde{v}_{0}$ is defined by (11), is accurate to at least $\mathcal{O}(h^{r+1})$ for the Euler–Lagrange boundary-value problem with boundary conditions $(q_{0},q_{1})$.

Proof.

As before, we denote the solution to the Euler–Lagrange boundary-value problem with boundary condition $(q_{0},q_{1})$ by $(q(t),v(t))$ for $t\in[0,h]$. This solution also satisfies the Euler–Lagrange initial-value problem with initial conditions $(q_{0},v_{0})$, where $v_{0}$ satisfies (12). We denote the solution of the Euler–Lagrange initial-value problem with initial conditions $(q_{0},\tilde{v}_{0})$ by $(\tilde{q}(t),\tilde{v}(t))$. Let $(q_{d}(t),v_{d}(t))$ denote the values generated by $r$-order Taylor method with initial conditions $(q_{0},\tilde{v}_{0})$. Noting that the Euler–Lagrange initial-value problem is well-posed, we denote the Lipschitz constant with respect to initial velocity by $M$.

[TABLE]

Combining this inequality with our $r$-order method yields,

[TABLE]

where we used the triangle inequality, and the fact that the local truncation error of a $r$-order Taylor method is $\mathcal{O}(h^{r+1})$ to bound the second term in line two, since $(\tilde{q}(t),\tilde{v}(t))$ and $(q_{d}(t),v_{d}(t))$ correspond to the exact solution and $r$-th order Taylor approximation, respectively, of the Euler–Lagrange initial-value problem with initial data $(q_{0},\tilde{v}_{0})$. ∎

Theorem 5**.**

Assuming a regular Lagrangian, we consider a fixed approximation scheme used to construct a corresponding discrete Lagrangian, $L_{d}$, and a discrete right Hamiltonian, $H_{d}^{+}$. This results in two integrators, $\tilde{F}_{L_{d}}:(q_{0},p_{0})\mapsto(q_{1,L_{d}},p_{1,L_{d}})$ and $\tilde{F}_{H_{d}^{+}}:(q_{0},p_{0})\mapsto(q_{1,H_{d}^{+}},p_{1,H_{d}^{+}})$. If the discrete right Hamiltonian approximation satisfies $p_{1,H_{d}^{+}}=D_{2}L_{d}(q_{0},\hat{q}_{1})$, where $\hat{q}_{1}$ is the approximated value of $q_{1}$, then the integrators represent the same map, i.e., $(q_{1,L_{d}},p_{1,L_{d}})=(q_{1,H_{d}^{+}},p_{1,H_{d}^{+}})$.

Proof.

Let $\hat{p}_{0}$ be defined by $-\hat{p}_{0}=D_{1}L_{d}(q_{0},\hat{q}_{1})$, where we consider $\hat{q}_{1}$ as an independent variable. The discrete right Hamiltonian is given by,

[TABLE]

Note that here $\hat{q}_{1}$ is being considered as a function of $q_{0}$ and $p_{1,H_{d}^{+}}$, as defined implicitly by the assumption $p_{1,H_{d}^{+}}=D_{2}L_{d}(q_{0},\hat{q}_{1})$. Then

[TABLE]

where the last line follows by the assumption $p_{1,H_{d}^{+}}=D_{2}L_{d}(q_{0},\hat{q}_{1})$. Therefore, $p_{0}=\hat{p}_{0}$, which then implies $-p_{0}=D_{1}L_{d}(q_{0},\hat{q}_{1})$ and consequently $\hat{q}_{1}=q_{1,L_{d}}$. Applying the next discrete Legendre transform yields,

[TABLE]

Therefore, $\hat{q}_{1}=q_{1,H_{d}^{+}}$, which implies $q_{1,L_{d}}=q_{1,H_{d}^{+}}$. Now we have,

[TABLE]

∎