A fixed time-step variational integrator cannot preserve momentum, energy, and symplectic form simultaneously for nonintegrable systems. This barrier can be overcome by treating time as a discrete dynamic variable and deriving adaptive time-step variational integrators that conserve the energy in addition to being symplectic and momentum-preserving. Their utility, however, is still an open question due to the numerical difficulties associated with solving the discrete governing equations. In this work, we investigate the numerical performance of energy-preserving, adaptive time-step variational integrators. First, we compare the time adaptation and energy performance of the energy-preserving adaptive algorithm with the adaptive variational integrator for Kepler's two-body problem. Second, we apply tools from Lagrangian backward error analysis to investigate numerical stability of the energy-preserving adaptive algorithm. Finally, we consider a simple mechanical system example to illustrate the backward stability of this energy-preserving, adaptive time-step variational integrator.

Keywords: Energy-preserving integrators; Variational integrators; Adaptive time-step integrators; Backward stability.

§ INTRODUCTION

Variational integrators are a class of structure-preserving integrators, that are derived by discretizing the action principle rather than the governing differential equation.  The basic idea behind these methods is to obtain an approximation of the action integral called the discrete action. Stationary points of the discrete action give discrete time trajectories of the mechanical system. These integrators have good long-time energy behavior and they conserve the invariants of the dynamics in the presence of symmetries.  The basics of variational integrators can be reviewed in  [1, 2]  and more detailed theory can be found in [3, 4]. A more general framework encompassing variational integrators, asynchronous variational integrators, and symplectic-energy-momentum integrators is discussed in [5]. Due to their symplectic nature,  variational integrators are ideal for long-time simulation of conservative or weakly dissipative systems [6] found in astrophysics and molecular dynamics. The discrete trajectories obtained using variational integrators display excellent energy behavior for exponentially long times.

The fixed time-step variational integrators derived from the discrete variational principle cannot preserve the energy of the system exactly.  To conserve the energy in addition to preserving the symplectic structure and conserving the momentum, time adaptation needs to be used.  These energy-preserving, adaptive time-step variational integrators were first developed for conservative systems in [7]  by imposing an additional energy preservation equation to compute the time-step. The same integrators were derived through a variational approach for a more general case of time-dependent Lagrangian systems in [3] by preserving the discrete energy obtained from the discrete variational principle in the extended phase space. Guibout and Bloch [8] have discussed the subtle differences between the underlying discrete variational principles. Recently, Sharma et al. [9, 10, 11] extended these methods to systems with external nonconservative forcing and rigid body motion.

Despite their excellent conservation properties, the utility of these methods is still an open question due to the difficulties associated with their numerical implementation. These adaptive algorithms are implicit and require solving a coupled nonlinear system of equations at every time-step to update both the configuration variables and the time variable. Unlike traditional adaptive algorithms, the adaptive time-step computation in these methods is inherently coupled to the discrete dynamics. In fact, existence of solutions for these discrete trajectories is not always guaranteed. Shibberu [12, 13] has discussed the well-posedness of these adaptive algorithms and identified points in the extended state space where the adaptive algorithm has no solutions. Even for points where solutions exist, the governing update equations are ill-conditioned. This introduces a discrete energy error at every iteration and can accumulate

in long-time simulations. Therefore, the backward stability of energy-preserving, adaptive time-step variational integrator  algorithms needs to be investigated.

The early development of backward error analysis in 1960s was motivated by numerical linear algebra applications and it was mainly used as a tool for evaluating the propagation of rounding errors in matrix algorithms. Although Wilkinson [14] first used it in the context of eigenvalue problems, the underlying idea of backward error analysis has also been applied to numerical integration of dynamical systems [15, 16]. Truncation errors introduced by the numerical integrator play the same role as rounding errors in the numerical linear algebra setting. The focus in the backward error analysis is to answer the question: what dynamical system does the numerical integrator solve? Thus, instead of studying the difference between the trajectories obtained using the numerical integrator and the exact solution, we look for a nearby dynamical system which would be solved exactly by the numerical integrator. For dynamical systems with underlying geometric structure, we would like the modified system to also possess similar features. The existence of such a nearby dynamical system has direct implications for the accurate long-time behavior of the simulation.

Sanz-Serna [17] first applied tools from backward error analysis to symplectic integrators and showed that numerical trajectories from symplectic integrators can be interpreted as the exact solution of a perturbed Hamiltonian system. Reich [18] used a recursive definition of modified vector fields to provide a unifying framework for the backward error analysis of geometric numerical integration (GNI) methods. Hairer [19] developed variable time-step symplectic integrators and proved backward stability of the adaptive algorithm by using time transformations. Recently, Vermeeran [20] studied the backward stability of fixed time-step variational integrators from the Lagrangian perspective.

The main goal of this work is to understand the numerical performance of energy-preserving, adaptive time-step variational integrators. First, we provide a comparison between different adaptive time-stepping approaches through a numerical example of Kepler's two-body problem. We also obtain a conservative bound on the discrete energy error that is accumulated from inexact solutions to the residual of the discrete energy equation. This shows that the error can be controlled through the residual tolerance and the numerical precision that is used. Variable precision arithmetic in the energy-preserving adaptive algorithm is used to show how the energy performance can be improved when using more significant digits in the computation. Motivated by the improvement in energy performance, we theoretically justify these results by investigating the stability of the adaptive algorithm. We first use a time-transformation to represent the adaptive algorithm as a fixed time-step algorithm and then apply backward error analysis tools in the transformed setting.

The remainder of the paper is organized as follows. In Section <ref>, we review the basics of extended Lagrangian mechanics and derive the energy-preserving, adaptive time-step variational integrators. In Section <ref>, we compare the numerical performance of energy-preserving, adaptive time-step variational integrators and adaptive variational integrators for the Kepler's two-body problem, with special focus on their energy behavior and time adaptivity. We also study the energy behavior of the adaptive algorithm and the effect of higher precision arithmetic on the conservation properties.  In Section <ref>, we use tools from backward error analysis to study algorithmic stability of energy-preserving, adaptive time-step variational integrators. We first introduce a time transformation to represent the energy-preserving, adaptive time-step variational integrator as a fixed time-step variational integrator. We then use these backward stability concepts to derive modified equations and the corresponding modified Lagrangian for a simple mechanical system. Finally, in Section <ref> we provide concluding remarks and suggest future research directions.

§ BACKGROUND

In this section, we review the basics of extended Lagrangian mechanics and the derivation of energy-preserving, adaptive time-step variational integrators. Drawing on the work of Marsden et al. [7, 21, 3], we first derive equations of motion in continuous-time from the variational principle and then derive the corresponding variational integrators by considering the discretized variational principle in the discrete-time domain.  To achieve this, we first derive extended Euler-Lagrange equations from Hamilton's principle. After deriving equations of motion, we use concepts of extended discrete mechanics developed in [3] to derive extended discrete Euler-Lagrange equations and write these in a time-marching form to obtain energy-preserving, adaptive time-step variational integrators for Lagrangian systems.

§.§ Extended Lagrangian Mechanics

Hamilton's principle of stationary action is one of the most fundamental results of classical mechanics and is commonly used to derive equations of motion for a variety of systems. The forces and interactions that govern the dynamical evolution of the system are easily determined through Hamilton's principle in a formulaic and elegant manner.  Hamilton's principle [22] states that: The motion of the system between two fixed points from $t_0$ to $t_f$ is such that the action integral has a stationary value for the actual path of the motion. To derive the Euler-Lagrange equations via Hamilton's principle, we start by defining the configuration space, tangent space and path space.
Consider a time-dependent Lagrangian system with configuration manifold $\Q$ and time space $\Real$. In the extended Lagrangian mechanics framework [3], we treat time as a dynamic variable and define the extended configuration manifold $\Qbar=\Real \times \Q$; the corresponding state space $\T\Qbar$ is $\Real \times \T\Q$, where the tangent bundle $\T\Q$ is the union of all tangent spaces to $\Q$. The extended Lagrangian is $L: \Real \times \T\Q \to \Real$. In the extended Lagrangian mechanics framework, $t$ and $\q$ are both parametrized by an independent variable $a$. The two components of a trajectory $c$ are $c(a)=(c_t(a),c_{\q}(a))$. The extended path space is

\begin{equation}
\bar{\mathcal{C}}=\left\{c:[a_0,a_f] \to \Qbar \ |\  \mbox{$c$ is a $C^2$ curve and $c'_t(a) >0$} \right\}.
\end{equation}

For a given path $c(a)$, the initial time is $t_0=c_t(a_0)$ and the final time is $t_f=c_t(a_f)$. The extended action $\bar{\mathfrak{B}} : \bar{\mathcal{C}} \to \Real$ is

\begin{equation}
\bar{\mathfrak{B}} = \int ^{t_f} _{t_0} L(t, \q(t), \qdot(t)) dt.
\label{eq:action}
\end{equation}

Since time is a dynamic variable in this framework, we substitute $(t,\q(t),\dot{\q}(t))=\left(c_t(a),c_{\q}(a),\frac{c'_{\q}(a)}{c'_t(a)} \right)$ in the above equation to get

\begin{equation}
\bar{\mathfrak{B}} = \int ^{a_f} _{a_0} L\left(c_t(a),c_{\q}(a),\frac{c'_{\q}(a)}{c'_t(a)} \right)c'_t(a) \  da.
\end{equation}

We compute the first variation of the action as

\begin{equation}
\delta \bar{\mathfrak{B}} = \int^{a_f}_{a_0} \left[ \frac{\partial L}{\partial t} \delta c_t + \frac{\partial L}{\partial \q}\cdot \delta c_{\q}  +  \frac{\partial L}{\partial \qdot}\cdot \left( \frac{\delta c'_{\q}(a)}{c'_t}  -\frac{c'_{\q} \delta c'_t(a) }{(c'_t)^2} \right) \right] c'_t(a)da + \int ^{a_f} _{a_0} L \delta c'_t(a) \  da .
\end{equation}

Using integration-by-parts and setting the variations at the endpoints to zero gives

\begin{equation}
\delta \bar{\mathfrak{B}} = \int^{a_f}_{a_0} \left[ \frac{\partial L}{\partial \q} c'_t - \frac{d}{da} \frac{\partial L}{\partial \qdot} \right] \cdot \delta c_{\q}(a)da + \int ^{a_f} _{a_0}  \left[ \frac{\partial L}{\partial t}c'_t + \frac{d}{da}\left( \frac{\partial L}{ \partial \qdot}\cdot \frac{c'_{\q}}{c'_t} - L \right) \right] \delta c_t(a) \  da .
\end{equation}

Using $dt=c'_t(a) da$ in the above expression gives two equations of motion. The first is the Euler-Lagrange equation of motion

\begin{equation}
\frac{\partial L}{\partial \q}   - \frac{d}{dt} \left( \frac{\partial L}{\partial \qdot}  \right) =\bzero,
\label{eq:contp}
\end{equation}

which is the same as the equation obtained using the classical Lagrangian mechanics framework. The second equation is

\begin{equation}
\frac{\partial L}{ \partial t} + \frac{d}{dt} \left( \frac{\partial L}{ \partial \qdot} \cdot \dot{\q} - L \right) =0,
\label{eq:conte}
\end{equation}

which describes how the energy of the system evolves with time.

It should be noted that both (<ref>) and (<ref>) depend only on the associated curve $\q(t)$. The time component $c_t(a)$ of the extended path cannot be determined from the governing equations. This comes from the fact that the energy evolution equation (<ref>) is redundant as it is a consequence of the Euler-Lagrange equation (<ref>). Thus, the "velocity" of time, i.e. $c'_t(a)$, is indeterminate in the continuous-time formulation.

§.§ Energy-preserving, Adaptive Time-step Variational Integrators

For the extended Lagrangian mechanics, we define the extended discrete state space $\Qbar \times \Qbar$. The extended discrete path space is

\begin{equation}
\bar{\mathcal{C}}_d = \{\ c: \{\ 0,...,N \} \to \Qbar\ | \ \mbox{$c_t(k+1) > c_t(k)$  for all $k$}  \}\ .
\end{equation}

The extended discrete action map $\bar{\mathcal{B}}_d : \bar{\mathcal{C}}_d \to \Real$ is

\begin{equation}
\bar{\mathcal{B}}_d = \sum _{k=0}^{N-1} L_d(t_k,\q_k,t_{k+1},\q_{k+1}),
\end{equation}

where $L_d : \Qbar \times \Qbar \to \Real$ is the extended discrete Lagrangian function that approximates the action integral between two successive configurations. Taking variations of the extended discrete action map leads to

\begin{multline}
\delta \bar{B}_d = \sum_{k=1}^{N-1} \left[ \D_4L_d(t_{k-1},\q_{k-1},t_{k},\q_{k}) + \D_2L_d(t_k,\q_k,t_{k+1},\q_{k+1}) \right] \cdot \delta \q_k  \\  +\sum_{k=1}^{N-1} \left[ \D_3L_d(t_{k-1},\q_{k-1},t_{k},\q_{k}) + \D_1L_d(t_k,\q_k,t_{k+1},\q_{k+1}) \right]  \delta t_k  = 0,
\end{multline}

where $\D_i$ denotes differentiation with respect to the $i^{th}$ argument of the discrete Lagrangian $L_d$. Applying Hamilton's principle of least action and setting variations at endpoints to zero results in the extended discrete Euler-Lagrange equations

\begin{equation}
\D_4L_d(t_{k-1},\q_{k-1},t_{k},\q_{k}) + \D_2L_d(t_k,\q_k,t_{k+1},\q_{k+1}) = \bzero,
\label{eq:vi_pvv}
\end{equation}

\begin{equation}
\D_3L_d(t_{k-1},\q_{k-1},t_{k},\q_{k}) + \D_1L_d(t_k,\q_k,t_{k+1},\q_{k+1}) = 0.
\label{eq:vi_ev}
\end{equation}

Given $(t_{k-1},\q_{k-1},t_{k},\q_{k})$, the extended discrete Euler-Lagrange equations can be solved to obtain $\q_{k+1}$ and $t_{k+1}$. This  extended discrete Lagrangian system can be seen as a numerical integrator of the continuous-time nonautonomous Lagrangrian system with adaptive time-steps.

In the extended discrete mechanics framework, we define the discrete momentum $\p_k$ by

\begin{equation}
\p_{k}= \D_4L_d(t_{k-1},\q_{k-1},t_{k},\q_{k}).
\end{equation}

We also introduce the discrete energy

\begin{equation}
\label{eq:dE}
\end{equation}

Using the discrete momentum and discrete energy definitions, we can re-write the extended discrete Euler-Lagrange equations (<ref>) and (<ref>) in the following form

\begin{align}
\label{eq:pimp}
-\D_2L_d(t_k,\q_k,t_{k+1},\q_{k+1}) &= \p_k,\\
\label{eq:eimp}
\D_1L_d(t_k,\q_k,t_{k+1},\q_{k+1}) &= E_k,\\
\label{eq:pexp}
\p_{k+1}&= \D_4L_d(t_k,\q_k,t_{k+1},\q_{k+1}),\\
\label{eq:eexp}
E_{k+1}&= -\D_3L_d(t_k,\q_k,t_{k+1},\q_{k+1}).
\end{align}

Given $(t_k,\q_k,\p_k,E_k)$, the coupled nonlinear equations (<ref>) and (<ref>) are solved implicitly to obtain $\q_{k+1}$ and $t_{k+1}$. The configuration $\q_{k+1}$ and time $t_{k+1}$ are then used in (<ref>) and (<ref>) to obtain $(\p_{k+1},E_{k+1})$ explicitly. The extended discrete Euler-Lagrange equations were first written in the  time-marching form in [3] and are also known as symplectic-energy-momentum integrators.

§ TIME TRANSFORMATION AND ADAPTIVITY

In this section, we compare the numerical performance of the energy-preserving, adaptive time-step variational integrators discussed in Section <ref> to the adaptive variational integrators from [23]. Due to their variational derivations, both classes of integrators are symplectic and momentum-preserving but motivation for the time adaptation is very different. Adaptive variational integrators are designed for computational efficiency whereas the energy-preserving, adaptive time-step variational integrators are focused on preserving the discrete energy in addition to being symplectic and momentum-preserving.

§.§ Adaptive Variational Integrators

For a Hamiltonian system $H(\q, \p)$, Schmitt and Leok  [23] prescribed an a priori time transformation $t \to a$, given by $\frac{dt}{da}=g(\q, \p)$, with auxilary states $(q^t,p^t)$ to form a new Hamiltonian

\begin{equation}
\bar{H}(\qbar,\bar{\p})=g(\q,\p) (H(\q,\p) + p^t),
\label{eq: newH}
\end{equation}

where $(\qbar,\pbar)=\left( \begin{bmatrix}
\q \\ q^t
\end{bmatrix},\begin{bmatrix}
\p \\ p^t
\end{bmatrix}\right)$. Subsequently, the physical time $t$ is treated as a configuration variable by choosing $q^t=t$ with the corresponding momentum variable $p^t=-H(\q_0,\p_0)$.  Hamilton's equations for the new Hamiltonian (<ref>) are given by

\begin{equation}
\frac{d\qbar}{da} =
\begin{bmatrix}
g(\q,\p) \frac{\partial H}{\partial \p} \\ g(\q,\p)
\end{bmatrix}, \quad \quad \quad  \frac{d\pbar}{da} =
\begin{bmatrix}
-g(\q,\p) \frac{\partial H}{\partial \q} \\ 0
\end{bmatrix}.
\end{equation}

Applying Hamiltonian variational integrators with fixed time-step $\Delta a$ to the transformed Hamiltonian system result in adaptive variational integrators for the original Hamiltonian system $H(\q,\p)$. We apply the implicit midpoint rule for time-integration of the transformed Hamiltonian system resulting in

\begin{equation*}
\begin{bmatrix}
\frac{\q_{k+1}-\q_k}{\da} \\
\frac{q^t_{k+1}-q^t_k}{\da}
\end{bmatrix} = \begin{bmatrix}
g(\q_{av})H_{\p}(\q_{av},\p_{av}) \\
\end{bmatrix}, \quad \begin{bmatrix}
\frac{\p_{k+1}-\p_k}{\da} \\
\frac{p^t_{k+1}-p^t_k}{\da}
\end{bmatrix} = \begin{bmatrix}
- g(\q_{av})H_{\q}(\q_{av},\p_{av}) \\
\end{bmatrix},
\end{equation*}

where $\q_{av}=\frac{\q_{k+1} + \q_k}{2}$, $\p_{av}=\frac{\p_{k+1} + \p_k}{2}$, and  we have assumed that the monitor functions only depend on the configuration, i.e. $g(\q,\p)=g(\q)$. Given $(\q_k, q^t_k, \p_k)$, we need to solve the following coupled nonlinear equations to  obtain $(\q_{k+1}, q^t_{k+1}, \p_{k+1})$

\begin{equation}
\label{eq:adap_qk}
\frac{\q_{k+1}-\q_k}{\da} =  g\left( \frac{\q_{k+1} + \q_k}{2}\right)\frac{\partial H}{\partial \p}\left(\frac{\q_{k+1} + \q_k}{2},\frac{\p_{k+1} + \p_k}{2}\right),
\end{equation}

\begin{equation}
\label{eq:adap_pk}
\frac{\p_{k+1}-\p_k}{\da} =  -g\left( \frac{\q_{k+1} + \q_k}{2}\right)\frac{\partial H}{\partial \q}\left(\frac{\q_{k+1} + \q_k}{2},\frac{\p_{k+1} + \p_k}{2}\right),
\end{equation}

\begin{equation}
\label{eq:adap_tk}
\frac{q^t_{k+1}-q^t_k}{\da} =  g\left( \frac{\q_{k+1} + \q_k}{2}\right).
\end{equation}

§.§ Numerical Comparison for Kepler's Two-body Problem

For the comparison of the above adaptive variational integrator to the energy-preserving, adaptive time-step variational integrator presented in Section <ref>, we consider Kepler's two-body problem with two different initial conditions. The Hamiltonian for this conservative system is

\begin{equation}
H(q^1,p^1,q^2,p^2)=\frac{(p^1)^2 +(p^2)^2}{2} - \frac{1}{\sqrt{(q^1)^2 +(q^2)^2}}.
\end{equation}

For initial conditions, we  consider $q^1(0)=1-e, q^2(0)=0, p^1(0)=0, p^2(0)=\sqrt{\frac{1+e}{1-e}}$ for $e=0.1$ and $0.7$, where $e$ is the eccentricity of the elliptical trajectory.

§.§.§ Implementation Details

* For energy-preserving, adaptive time-step variational integrator (EpAVI), we use the midpoint rule to obtain the extended discrete Lagrangian

\begin{equation}
L_d=(t_{k+1}-t_k)L\left( \frac{\q_k + \q_{k+1}}{2},\frac{\q_{k+1}-\q_k}{t_{k+1}-t_k}\right).
\label{eq:Ld_comp}
\end{equation}

The time-marching equations for the EpAVI algorithm can be derived by substituting the discrete Lagrangian from (<ref>) in (<ref>)-(<ref>).

* For adaptive variational integrator (AVI), we consider two monitor functions used in [23]. The first is the arclength parametrization

\begin{equation}
g_1(\q)= (2(H_0-V(\q)) + \nabla V(\q)^T M^{-1} \nabla V(\q) )^{-\frac{1}{2}} ,
\end{equation}

and the second is the problem-specific monitor function

\begin{equation}
\end{equation}

which is motivated by Kepler's second law. We denote the adaptive algorithms based on the first monitor function $g_1$ and second monitor function $g_2$ by AVI1 and AVI2, respectively.

* All numerical studies are done with MATLAB 2020b. To compute the trajectory error plots, we compute the reference solution using the in-built ode45 solver with $reltol=1e-12$ and $abstol=1e-14$.

* For a fair comparison between the two adaptive approaches, $\Delta a$ for the adaptive algorithm (<ref>)-(<ref>) is chosen such that the resulting initial time step $h_0$ is same as the one used for the energy-preserving adaptive algorithm.

§.§.§ Results

We compare the energy error, trajectory error, and adaptive time-step behavior for $e=0.1$ and $e=0.7$ in Figure <ref> and Figure <ref>, respectively. For $e=0.1$, we observe that EpAVI algorithm and AVI1 algorithm give nearly same adaptive time-step behavior, while AVI2 algorithm achieves marginally higher adaptive time-step values. The energy error plots in Figure  <ref> demonstrate that AVI1 and AVI2 exhibit bounded energy error magnitude around $10^{-7}$ whereas EpAVI algorithm has energy error magnitude around $10^{-14}$. Trajectory error plots for $q_1$ and $q_2$ demonstrate the improved accuracy achieved by EpAVI algorithm. All three algorithms presented exhibit nearly  the same accuracy, with the EpAVI algorithm performing marginally better.
For the higher eccentricity value $e=0.7$, the adaptive variational integrators (AVI1 and AVI2) exhibit significantly higher adaptive time-step values compared with the EpAVI algorithm. Unlike the $e=0.1$ case where the adaptive time step $h_k$ does not increase substantially from the initial time step $h_0=0.001$, all three adaptive algorithms achieve significant speedup with the adaptive time-step value $h_k$ increasing by almost 12-15 times the initial value $h_0$. Both adaptive variational integrators have energy error around $10^{-4}$ whereas the EpAVI algorithm has energy error around $10^{-13}$. Similar to the $e=0.1$ case, the AVI1 and AVI2 algorithms exhibit similar trajectory error, and the EpAVI algorithm performs marginally better.

[subfigure]oneside,margin=2.0cm,0 cm

Adaptive time step

Energy error

$q^1$ trajectory error

$q^2$ trajectory error

Comparison between adaptive variational integrators for $e=0.1$

[subfigure]oneside,margin=2.0cm,0 cm

Adaptive time step

Energy error

$q^1$ trajectory error

$q^2$ trajectory error

Comparison between adaptive variational integrators for $e=0.7$

[subfigure]oneside,margin=1.5cm,0 cm

Time transformation $t(a)$ for $h_0=0.001$

Based on the numerically computed discrete time values $\{t_k\}_{k=0}^{N}$,  we have compared the time transformation $t(a)$ in Figure <ref>  for both adaptive approaches. These plots reveal the numerical relation between the physical time $t$ and the transformed time $a$ over one cycle for both eccentricity cases. In Figure <ref>, both energy-preserving adaption and monitor function approaches have a similar $t(a)$ relation for $e=0.1$. For an initial time-step of $h_0=0.001$, the AVI2 algorithm achieves an average time-step of $1.23h_0$ whereas the average time-step for the EpAVI algorithm is $1.17h_0$. By contrast, Figure <ref> shows how both monitor function and energy-preserving adaptation approaches have significantly higher average time-step for $e=0.7$. The AVI2 algorithm achieves an average time-step of $7.93h_0$ whereas the average time-step for the EpAVI algorithm is $6.22h_0$. Furthermore, $t(a)$ plots in Figure <ref> numerically validate the bounded nature of $t'(a)$, i.e. $ 0<t'(a)\leq const$. The fact that $t'(a)$ remains bounded will be used as an assumption in the backward error analysis in Section <ref>.

[subfigure]oneside,margin=1.5cm,0 cm

Adaptive time step ratio

Discrete energy error

Effect of initial time-step $h_0$ on adaptive variational integrators for $e = 0.7$

We have also studied the effect of initial time-step $h_0$ on the numerical performance of the adaptive algorithms. For the AVI2 algorithm, increasing the initial time-step from $h_0=0.001$ to $h_0=0.01$ leads to an increase in the energy error. The EpAVI algorithm, on the other hand, exhibits small decrease in the energy error for the smaller initial time-step $h_0=0.001$. This unexpected decrease in energy error for larger time-steps is due to residual error in not solving the discrete energy equation exactly. We discuss the dependence of the residual error on energy behavior of the EpAVI algorithm in detail in Section <ref>. Both algorithms demonstrate similar numerical behavior of $\frac{h_k}{h_0}$ for both  $h_0$ values which shows that the speedup achieved up the adaptive algorithms is independent of initial time-step and is mainly governed by the discrete dynamics.

We observe slow linear error growth for all three adaptive algorithms in the trajectory error plots for both eccentricity values in Figure <ref> and Figure <ref>. This error growth is due to the integrable nature of Kepler's two-body problem. Completely integrable systems form an interesting class of Hamiltonian systems, and symplectic integrators are known to exhibit linear error growth for such systems. More details about this result can be found in [24].

§.§ Energy Behavior

As shown in Section <ref>, the extended discrete equations for the adaptive time-step variational integrators conserve the discrete energy $E_k$ exactly. Unfortunately, numerical algorithms do not solve the discrete energy equation exactly and the residual from approximating solutions to this equation affects the errors in subsequent iterations. It is possible that the accumulation of the discrete energy error leads to a significant energy error that overcomes the energy-preservation advantage of the adaptive time-stepping algorithm. Thus, we need to investigate the long-term behavior of energy error.
For a given initial configuration $(t_0,\q_0,\p_0,E_0)$, the discrete energy $E_k$ after $k$ adaptive time-steps can be written as the telescoping series

\begin{equation}
E_k= E_0 + \sum_{i=1}^k  (E_i -E_{i-1}).
\end{equation}

We can use this expression to obtain a conservative bound on the discrete energy error

\begin{equation}
| E_k- E_0| \leq  \sum_{i=1}^k | E_k -E_{k-1}| \leq k \max_i |E_i -E_{i-1} |,
\label{eq: dEbound}
\end{equation}

where $\max_i|E_i-E_{i-1}|$ is the maximum discrete energy error introduced in one iteration during the entire numerical simulation. Using this conservative bound, we can compare the energy behavior of fixed time-step variational integrators with the worst case scenario for the EpAVI algorithm. For example, if we consider an example where the fixed time-step algorithm exhibits discrete energy error around $10^{-6}$ then the EpAVI algorithm, with $10^{-15}$ accuracy in solving discrete energy equation, will exhibit better energy performance for $10^{9}$ steps

\begin{equation}
k=\frac{10^{-6}}{\max_i |E_i -E_{i-1} |}= \frac{10^{-6}}{10^{-15}}=10^{9}  \ steps .
\end{equation}

[subfigure]oneside,margin=1.5cm,0 cm

Energy error

Adaptive time step

EpAVI algorithm for $e=0.7$ using VPA for different tolerance values.

It is important to note that the discrete energy error bound in (<ref>) is very conservative because it assumes the worst case scenario where  the energy error is maximal and of the same sign in every step. We have studied Kepler's two-body problem using Variable Precision Arithmetic (VPA) to demonstrate that the EpAVI algorithm can resolve the issue of accumulation of discrete energy error by solving the discrete equations with higher precision. We have compared the numerical performance of the EpAVI algorithm using VPA with 18 significant digits for different tolerance values in Figure <ref>. The higher precision with VPA allows us to solve the discrete update equations with lower tolerance values. The energy error plots in Figure <ref> show that the accumulation of discrete energy error decreases with decreasing tolerance, and VPA with $tolerance=10^{-17}$ achieves exact energy conservation up to precision. The adaptive time-step behavior in Figure <ref> also matches our previous results from Section <ref>.
Based on these results, we have shown that increasing the precision can lead to nearly exact energy behavior without requiring time-step sizes to vanish. It is important to note that using VPA makes the EpAVI algorithm computationally prohibitive. The key takeaway from these results obtained using VPA is that, if solved with enough precision, the adaptive algorithm can conserve the discrete energy exactly. In fact, the energy error results obtained by using VPA strongly indicate that the EpAVI algorithm is a backward stable algorithm.

§ BACKWARD ERROR ANALYSIS

There are two challenges in showing backward stability for energy-preserving, adaptive time-step variational integrators. First, these algorithms are derived from the Lagrangian perspective and hence, the governing Euler-Lagrange equations are second-order.  The second challenge is the unique nature of the adaptive algorithm where the adaptive time-step is computed from the discrete Euler-Lagrange equations (<ref>)-(<ref>) and thus the time transformation map $t(a)$ in (<ref>) is not explicitly given in the continuous-time setting. This is related to the fact that the Euler-Lagrange equation (<ref>) corresponding to the time variable $t$ is a redundant equation.
Vermeeren [20] developed the theoretical framework for backward error analysis from the Lagrangian side and derived modified equations for fixed time-step variational integrators. Our goal is to apply backward stability results from [20] to energy-preserving, adaptive time-step variational integrators. We first introduce a time transformation and show that the energy-preserving adaptive algorithm is equivalent to a fixed time-step variational integrator in the transformed setting. We then review some concepts required for the backward error analysis of variational integrators from a Lagrangian perspective and also modify some definitions to incorporate the time transformation.

[node distance = 1cm, auto,font=,

every node/.style=node distance=3cm,

comment/.style=rectangle, inner sep= 5pt, text width=5cm, node distance=0.25cm, font=,

force/.style=rectangle, draw, fill=black!10, inner sep=5pt, text width=5cm, text badly centered, minimum height=1.2cm, font=]

[force] (VI) Fixed time-step variational integrators for the transformed Lagrangian $\bar{L}(\qbar(a),\qbar'(a))$;
[force, left=3 cm of VI] (EpVI) Energy-preserving, adaptive time-step variational integrators for Lagrangian $L(\q(t), \dot{\q}(t))$;

[force, below of=VI] (Emod) Consistent discretization $\bar{L}_{\rm disc} (\qbar_k,\qbar_{k+1},\Delta a)$;
[force, below of=Emod] (Lmesh) Meshed modified Lagrangian $\mathfrak L_{\rm mesh}([\qbar(a)],\Delta a)$;
[force, left=3 cm of Lmesh] (Lmod) Truncated modified Lagrangian $\mathfrak L_{mod,r}(\q(a),\q'(a),[t(a)],,\Delta a)$;

(EpVI) edge (VI)
(VI) edge (Emod)
(Emod) edge (Lmesh)
(Lmesh) edge (Lmod);
(EpVI) edge (Lmod);

Schematic overview of the backward error analysis of energy-preserving, adaptive time-step variational integrators.

§.§ Time Transformation

In this section, we introduce a time transformation from the physical time $t$ to a fictitious time $a$ and represent the energy-preserving, adaptive time-step variational integrator (EpAVI) as a fixed time-step variational integrator in this transformed setting.

Using ideas from extended Lagrangian mechanics, we interpret the parameterization of $t=t(a)$ in terms of the independent variable $a$ as a time transformation. For such a time transformation, i.e. $t \to a$, the configuration in the transformed setting is given by $\qbar(a)=(\q(a),t(a))$. For this new set of states, we introduce

\begin{equation}
\bar{L}((\qbar(a),\qbar'(a)))=t'(a)L(\q(t),\qdot(t))=t'(a)L\left( \q(a),\frac{\q'(a)}{t'(a)}\right),
\label{eq:t(a)}
\end{equation}

where we have used $\qdot(t)=\frac{\q'(a)}{t'(a)}$. The Euler-Lagrange equations in the transformed setting are

\begin{align}
\frac{\partial \bar{L}}{\partial \q} - \frac{d}{da} \left(  \frac{\partial \bar{L}}{\partial \q'} \right)&= t'(a)  \frac{\partial L}{\partial \q} - \frac{d}{da} \left( t'(a)  \frac{\partial L}{\partial \qdot} \frac{1}{t'(a)} \right)   = t'(a) \left( \frac{\partial L}{\partial \qdot} - \frac{d}{dt} \left(  \frac{\partial L}{\partial \qdot} \right)   \right)  =\bzero,
\label{eq:CEL}\\
\frac{\partial \bar{L}}{\partial t} - \frac{d}{da} \left(  \frac{\partial \bar{L}}{\partial t'} \right)  &= -\frac{d}{da} \left( L + t'(a) \frac{\partial L}{\partial \qdot} \left(-\frac{\qdot}{t'(a)}\right) \right) = t'(a) \frac{d}{dt} \left(  \frac{\partial L}{\partial \qdot} \cdot \dot{\q} - L \right) = 0.
\label{eq:CEE}
\end{align}

Thus, the Euler-Lagrange equations in the transformed setting lead to the original Euler-Lagrange equation (<ref>) along with the redundant energy equation (<ref>). Although, these equations are the same as the extended Euler-Lagrange equations, it is important to interpret the above equations as Euler-Lagrange equations in the transformed setting. This interpretation will play a key role in the following Lemma about the equivalence of trajectories of the original Lagrangian system and the transformed system.

Lemma 1: Let $\q(t)$ and $\q(a)$ be the trajectory solutions to the Euler-Lagrange equations corresponding to the physical Lagrangian $L$ and the transformed Lagrangian $\bar{L}$, respectively. Then,

\begin{equation*}
\qbar(T)=(\q(\alpha(T)),\alpha(T)),
\end{equation*}

where $\alpha(T)=\int_0^T t'(a) da$.

Proof: This follows from the fact that $\qbar(T)$ and $\q(\alpha(T))$ satisfy the original Euler-Lagrange equation with same initial conditions.

Consider a discrete trajectory in the transformed setting with $\Delta a=a_{k+1}-a_k$. The discrete action $\bar{\mathfrak B}_d$ approximates the action integral

\begin{equation}
\bar{\mathfrak B}_d =  \sum_{i=1}^N    \bar{L}_d(\qbar_k,\qbar_{k+1}) \approx \int_{a_0}^{a_f} \bar{L}(\qbar(a),\qbar'(a)) da,
\end{equation}

where $\qbar_k=(\q_k,t_k)$ and the discrete Lagrangian is given by

\begin{equation}
\bar{L}_d(\qbar_k,\qbar_{k+1})=\Delta a  \left[ \bar{L}\left(\frac{\qbar_k+ \qbar_{k+1}}{2},\frac{\qbar_{k+1}-\qbar_k}{\Delta a}\right)  \right] \approx \int_{a_k}^{a_{k+1}} \bar{L}(\qbar(a),\qbar'(a)) \ da .
\end{equation}

The resulting discrete Euler-Lagrange equations are given by

\begin{equation}
\frac{\partial \bar{L}_d(\qbar_{k-1},\qbar_{k})}{\partial \q_k} + \frac{\partial \bar{L}_d(\qbar_k,\qbar_{k+1})}{\partial \q_k} =  \bzero,
\label{eq:DEL}
\end{equation}

\begin{equation}
\frac{\partial \bar{L}_d(\qbar_{k-1},\qbar_{k})}{\partial t_k} + \frac{\partial \bar{L}_d(\qbar_k,\qbar_{k+1})}{\partial t_k} =0.
\label{eq:DEE}
\end{equation}

These variational integrators with a fixed time-step $\Delta a$ in the transformed setting are equivalent to the energy-preserving, adaptive time-step variational integrators. This equivalence will allow us to apply established results from backward error analysis of fixed time-step variational integrators.

§.§ Consistent Discretization and Modified Equations

Since variational integrators are written in terms of the discrete  Lagrangian $\bar L_d$, the first step is to define a consistent discretization of the continuous Lagrangian $\bar{L}(\qbar,\qbar')$. We define $\bar{L}_{\rm disc}$ as follows

\begin{equation}
\bar{ \mathfrak{B}}_d=\sum_{k=0}^{N-1}\Delta a \bar L_{\rm disc}(\qbar_k,\qbar_{k+1},\Delta a),
\end{equation}

where $\bar{\mathfrak{B}}_d$ is the discrete action in the transformed setting. It is important to note that $\bar{L}_{disc}$ is different from the discrete Lagrangian $\bar  L_d$ defined earlier. In fact, $\bar L_d$ and $\bar L_{disc}$ are related by a nonzero scaling factor of $\Delta a$, i.e. $\bar L_d=\Delta a \bar L_{disc}$. A consistent discretization $\bar L_{disc}$ of $\bar L$ ensures that the discrete Euler-Lagrange equations (<ref>)-(<ref>)  are a consistent discretization of continuous Euler-Lagrange equations (<ref>)-(<ref>).  For example, the midpoint rule used in this work is an example of consistent discretization with

\begin{equation*}
\bar{L}_{\rm disc}= \left[ \bar{L}\left(\frac{\qbar_k+ \qbar_{k+1}}{2},\frac{\qbar_{k+1}-\qbar_k}{\Delta a}\right)  \right] .
\end{equation*}

As an example, consider the Kepler's two-body problem from Section <ref>. The governing discrete equations based on the midpoint rule

\begin{equation}
\left(\frac{q^1_{k+1}-q^1_k}{t_{k+1}-t_k} \right) - (t_{k+1}-t_k) \left( \frac{q^1_k + q^1_{k+1}}{4\left(\left(q^1_k + q^1_{k+1}\right)^2 +\left(q^2_k + q^2_{k+1}\right)^2\right)^{\frac{3}{2}}}\right)=p^1_k,
\end{equation}

\begin{equation}
\left(\frac{q^2_{k+1}-q^2_k}{t_{k+1}-t_k} \right) - (t_{k+1}-t_k) \left( \frac{q^2_k + q^2_{k+1}}{4\left(\left(q^1_k + q^1_{k+1}\right)^2 +\left(q^2_k + q^2_{k+1}\right)^2\right)^{\frac{3}{2}}}\right)=p^2_k,
\end{equation}

\begin{equation}
\frac{1}{2} \left( \frac{q^1_{k+1}-q^1_k}{t_{k+1}-t_k} \right)^2 + \frac{1}{2} \left( \frac{q^2_{k+1}-q^2_k}{t_{k+1}-t_k} \right)^2 - \frac{2}{\sqrt{\left(q^1_k + q^1_{k+1}\right)^2 +\left(q^2_k + q^2_{k+1}\right)^2}}=E_k.
\end{equation}

These are a consistent discretization of the continuous Euler-Lagrange equations.
Since the governing equations from the Lagrangian perspective, i.e. Euler-Lagrange equations, are second-order equations, the modified equation for the discrete Euler-Lagrange equations should also be second-order. As discussed in Section <ref>, the energy-preserving, adaptive time-step variational integrators are based on the extended Lagrangian mechanics and its discrete analogue.  To use tools from standard backward analysis of numerical integrators, we first need to define what constitutes a modified equation in the transformed setting.  We define $ \pmb{\psi}=[\pmb{\psi}_{EL},\psi_E]$  where $\pmb \psi_{EL}$ is the discrete Euler-Lagrange equation

\begin{equation}
\pmb \psi_{EL}(\qbar_{k-1},\qbar_{k},\qbar_{k+1},\Delta a) =  \frac{\partial \bar{L}_d(\qbar_{k-1},\qbar_{k})}{\partial {\q}_k} + \frac{\partial \bar{L}_d(\qbar_k,\qbar_{k+1})}{\partial {\q}_k} =  \bzero,
\end{equation}

and $\psi_E$ is the discrete energy equation

\begin{equation}
\psi_E(\qbar_{k-1},\qbar_{k},\qbar_{k+1},\Delta a) =  \frac{\partial \bar{L}_d(\qbar_{k-1},\qbar_{k})}{\partial t_k} + \frac{\partial \bar{L}_d(\qbar_k,\qbar_{k+1})}{\partial t_k} =  0.
\end{equation}

Definition: The differential equation $\qbar''=\mathbf f(\qbar(a),\qbar'(a),\Delta a)$ is a modified equation for the discrete equation  $\psi(\qbar_{k-1},\qbar_{k},\qbar_{k+1},\Delta a)=\bzero$ if, for every $r$, every admissible family $(\qbar_{\Delta a})$ of solutions of

\begin{equation}
\qbar''_{\Delta a}=\tau_r  (\mathbf f(\qbar_{\Delta a},\qbar'_{\Delta a},\Delta a)) ,
\end{equation}

where $\tau_r$ denotes truncation after the $\Delta a^r$- term, satisfies

\begin{equation*}
\pmb  \psi(\qbar_{\Delta a}(a-\Delta a),\qbar_{\Delta a}(a),\qbar_{\Delta a}(a+\Delta a),\Delta a ) = O(\Delta a^{r+1}),
\end{equation*}

as $\Delta a\to 0$ for all $a$.

A family of curves $(\qbar_{\Delta a})_{\Delta a \in (0,\infty)}$ of smooth curves $\qbar_{\Delta a}:[a_h,b_h] \to \Real^N$ is called admissible if there exists $\Delta a_{max}>0$ such that for each $k\geq 0,||\qbar_{\Delta a}^{(k)}||_{\infty} $ is bounded as a function of $\Delta a \in (0,\Delta a_{max}]$, where $|| \cdot ||_{\infty}$ denotes the supremum norm. These admissible families are families of curves parametrized by $\Delta a$ whose derivatives do not blow up when $\Delta a \to 0$. Admissibility of a family of smooth curves guarantees that in power series expansions like $\qbar_{\Delta a}(a+\Delta a)=\qbar_{\Delta a}(a) + \Delta a \qbar'_{\Delta a}(a) + \cdots$ the asymptotic behavior of each term is determined by the exponent of $\Delta a$ in that term.

§.§ Modified Lagrangian

For a variational integrator in the transformed setting with a modified equation $\pmb \psi$, ideally we would like to construct a modified Lagrangian $\bar L_{\rm mod}$ such that modified equation $\pmb \psi$ is found through the Euler-Lagrange equations. However, the modified equations generally have the form of non-converging power series so we will focus on the $r-$truncation of the modified equation. Thus for the backward error analysis, our goal is to obtain a modified Lagrangian $\bar{L}_{\rm mod}$ such that the $r-$truncation of the Euler-Lagrange equations for the modified Lagrangian are the $r-$truncation of the modified equations corresponding to the energy-preserving, adaptive algorithm.
If $\qbar$ is a curve such that $(\qbar(a_0 + k\Delta a))_{k\in [0,\cdots, n]}$ is critical for the discrete action $\bar{\mathfrak{B}}_d$, then $\qbar$ solves the meshed variational problem for $\bar{\mathfrak L}_{\rm mesh}$

\begin{equation}
\bar{  \mathfrak{B}}_d((\qbar(a_0 + k\Delta a)))= \sum_{k}\bar{L}_d(\qbar_k,\qbar_{k+1},\Delta a) \simeq \int^{a_0 + n\Delta a}_{a_0}\bar{\mathfrak L}_{\rm mesh} ([\qbar(a)],\Delta a) \ da ,
\end{equation}

where $[\qbar(a)]$ contains $\qbar(a)$ and its derivatives w.r.t fictitious time $a$. The meshed modified Lagrangian $\bar{\mathfrak L}_{\rm mesh}$ is defined by the following formal power series

\begin{equation}
\bar{\mathfrak L}_{\rm mesh}([\qbar],\Delta a):=  \bar{\mathfrak L}_{\rm disc}([\qbar],\Delta a) -\frac{\Delta a^2}{24}\frac{d^2}{dt^2}\bar{\mathfrak L}_{\rm disc}([\qbar],\Delta a) + \frac{7\Delta a^4}{5760}\frac{d^4}{dt^4}\bar{\mathfrak L}_{\rm disc}([\qbar],\Delta a) + \cdots  \ .
\label{eq:Lmeshmod}
\end{equation}

Since $\bar{\mathfrak L}_{\rm disc}([\qbar],\Delta a)$ is a consistent discretization of $\bar L$, it is clear that zeroth order term of the meshed modified is the original continuous Lagrangian, i.e. $\bar {\mathfrak L}_{\rm mesh}([\qbar],\Delta a)=\bar L(\qbar, \qbar') + O(\Delta a)$.

It is evident from (<ref>) that the meshed modified Lagrangian $\bar{\mathfrak L}_{mesh}$ depends on higher derivatives of $\qbar$ instead of just $\qbar$ and $\qbar'$. This is related to the unconventional nature of the meshed variational problem that defines $\bar{\mathfrak L}_{mesh}$. Unlike the classical variational problem that seeks critical curves of the action integral (<ref>) within the set of smooth curves $\mathcal C^{\infty}$, the meshed variational problem with mesh size $\Delta a$ seeks smooth curves that for every $a_0 \in \Real$ are critical for the action integral within the set of piecewise smooth curves $\mathcal M^{a_0,\Delta a}$ whose nonsmooth points lie in the mesh $a_0 + \Delta a \mathbb{Z}$. For a thorough exposition, the interested reader may consult
Section 4 in [20].

Based on Lemma $8$ from [20], the modified equation for variational integrators can be calculated from the Euler-Lagrange equations of the meshed modified Lagrangian. Due to the presence of higher order derivatives in $\bar {\mathfrak L}_{\rm mesh}$, these modified equations are written as a power series in the fixed time-step $\Delta a$ and are usually non-convergent. In such cases, we are often only interested in truncated power series expressions so we also need to relax our notion of critical curves for both continuous and discrete variational problems. Consider the following definitions:

* Classical $r$–critical: An admissible family $(\qbar_{\Delta a})_{\Delta a \in (0, \infty)}$ of curves $\qbar_{\Delta a}:[0,1] \to \Real$ is $r$–critical for a family of actions $\bar{\mathfrak{B}}$ if every family of smooth variations $\delta \qbar_{\Delta a} $ satisfies $\delta \bar{\mathfrak{B}}_{\Delta a}=O(\Delta a^{r+1})|| \delta \qbar_{\Delta a} ||$. An admissible family $(\qbar_{\Delta a})_{\Delta a \in(0,\infty)}$ of curves $\qbar_{\Delta a}:[0,1] \to \Real$ is $r-$critical for a family of actions $\bar{\mathfrak{B}}_{\Delta a}=\int_0^1\bar{L}_{\Delta a} dt$ if and only if

\begin{equation}
\bigg| \bigg| \frac{\delta \bar{L}_{\Delta a}}{\delta \qbar_{\Delta a}} \bigg|\bigg|_{\infty}=O(\Delta a^{r+1}) .
\end{equation}

* Meshed $r$–critical: An admissible family $(\qbar_{\Delta a})_{\Delta a \in (0, \infty)}$ of curves $\qbar_{\Delta a}:[0,1] \to \mathbb{R}$ is  meshed $r$–critical for a family of actions $\bar{\mathfrak{B}}$ if for every family of smooth variations $\delta \qbar_{\Delta a} \in \mathcal{M}^{a_{\Delta a},\Delta a} $ (i.e. piecewise smooth $\delta \qbar$ with nonsmooth points in a mesh of size $\Delta a$) we have $\delta \bar{\mathfrak{B}}_{\Delta a}=O(\Delta a^{r+1})|| \delta \qbar_{\Delta a} ||$. Similarly, an admissible family  $(\qbar_{\Delta a})_{\Delta a \in(0,\infty)}$ of curves $\qbar_{\Delta a}:[0,1] \to \mathbb{R}$ is meshed  $r$–critical for a family of actions $\bar{\mathfrak{B}}_{\Delta a}=\int_0^1\bar{L}_{\Delta a} dt$  if and only if

\begin{equation}
\bigg| \bigg|  \frac{\delta \bar{L}_{\Delta a}}{\delta \qbar_{\Delta a}} \bigg|\bigg|_{\infty}=O(\Delta a^{r+1}) \quad \mbox{and} \quad  \bigg| \bigg|  \frac{\delta \bar{L}_{\Delta a}}{\delta \qbar_{\Delta a}^{(l)}} \bigg|\bigg|_{\infty}=O(\Delta a^{r+l+1}) \ \mbox{for each} \ l \geq 2 .
\end{equation}

* Discrete $r$–critical: A family $(\qbar_{\Delta a})_{\Delta a \in (0, \infty)}$ of discrete curves $\qbar_{\Delta a}:[0,1] \to \mathbb{R}$ is $r-$critical for a family of discrete actions $\bar{\mathfrak{B}}_d$ if for every family of discrete variations $\delta \qbar_{\Delta a} $ we have $\delta \bar{\mathfrak{B}}_{d, \Delta a}=O(\Delta a^{r+1})|| \delta \qbar_{\Delta a} ||$. A family of discrete curves $\qbar_{\Delta a}=(\qbar_{\Delta a,k})_{k \in [0, \cdots,n_k]}$ is $r$–critical for a family of discrete actions $\bar{\mathfrak{B}}_{d,\Delta a}=\sum \bar{L}_d$ if and only if

\begin{equation}
\frac{\partial  \bar{L}_d}{\partial \qbar_{\Delta a,k}}(\qbar_{\Delta a,k},\qbar_{\Delta a,k+1}) + \frac{\partial  \bar{L}_d}{\partial \qbar_{\Delta a,k}}(\qbar_{\Delta a,k-1},\qbar_{\Delta a,k})=O(\Delta a^{r+1}) .
\end{equation}

The $r-$critical families of curves for the meshed modified Lagrangian $\bar{\mathfrak{L}}_{\rm mesh}$ satisfy

\begin{equation}
\frac{\partial \bar{\mathfrak{L}}_{mesh}}{\partial {\q}}  - \frac{d}{da} \left( \frac{\partial \bar{\mathfrak L}_{\rm mesh}}{\partial \q'} \right) = O(\Delta a^{r+1}).
\end{equation}

Based on Proposition $6$ from [20], the modified equation can be written in the following form

\begin{equation}
\epsilon_0(\qbar, \qbar', \qbar'') + \Delta a \epsilon_1(\qbar, \qbar', \qbar'') + \Delta a^2{\epsilon}_2(\qbar, \qbar', \qbar'',\qbar^{(3)}) + \cdots + \Delta a^r {\epsilon}_r(\qbar, \qbar', \cdots ,\qbar^{(r+1)})= O(\Delta a^{r+1}).
\label{eq:recursion}
\end{equation}

where for $r \geq 1$, $\epsilon_r$ depends on $\qbar, \qbar', \cdots, \qbar^{(r+1)}$ but not on higher derivatives of $\qbar$. Using $\bar {\mathfrak L}_{\rm mesh}([\qbar],\Delta a)=\bar L(\qbar, \qbar') + O(\Delta a)$, it is easy to show that $\epsilon_0(\qbar, \qbar',\qbar'')$ is simply the continuous Euler-Lagrange equation (<ref>). For a consistent discretization $\bar L_{\rm disc}$ and sufficiently small $\Delta a$, we can solve $\epsilon_0(\qbar, \qbar', \qbar'') + \Delta a \epsilon_1(\qbar, \qbar', \qbar'')=O(\Delta a^2)$ for $\q''$, i.e.

\begin{equation}
\q''=\pmb F_{1}(\q,\q',[t],\Delta a) + O(\Delta a^{2}).
\end{equation}

Following the recursion proposed in Section 4.6 of [20], we can show that the modified equation (<ref>) can be solved for $\q''$ in terms of $(\q,\q',[t],\Delta a)$

\begin{equation}
\q''=\pmb F_{r}(\q,\q',[t],\Delta a) + O(\Delta a^{r+1}).
\end{equation}

It is important to note that we can only solve the modified equation for $\q$ and not for $\qbar$.  In other words, we can solve (<ref>) for $\q''$ in terms of $\q,\q'$ and arbitrary number of derivatives of $t$.  However, the redundant nature of (<ref>)  prevents us from solving the modified equations for $\qbar=[\q,t]$ jointly for $\qbar''$ in terms of $\qbar$ and $\qbar'$.

For the fixed time-step variational integrator considered in  [20], the author used the properties of admissible curves to write the modified equation as a second order differential equation with $O(\Delta a^{r+1})$ defect. This was achieved by a recursive process in Section 4.6 of [20] where the leading order terms in the modified equation are solved for $\ddot{\q}$. This result doesn't carry over to our transformed setting due to the redundant nature of the continuous Euler-Lagrange equation (<ref>) corresponding to the time variable.

Unlike the classical Lagrangian defined on the tangent bundle, the meshed modified Lagrangian $\bar{\mathfrak L}_{\rm mesh}$ also depends on higher derivatives of $\qbar(a)$. From the meshed modified Lagrangian $\bar {\mathfrak L}_{\rm mesh}$, we define the classical modified Lagrangian $\bar{\mathfrak L}_{\rm mod}$

\begin{equation}
\bar{\mathfrak L}_{\rm mod}(\q,\q',[t],\Delta a)=\bar{\mathfrak L}_{\rm mesh}([\qbar],\Delta a)\bigg|_{\q''=\pmb \psi(\q, \q',[t],\Delta a), \q'''= \frac{d}{da} \pmb \psi( \q, \q',[t],\Delta a), \cdots  },
\end{equation}

where $\q''=\pmb \psi(\q,\q',[t],\Delta a)$ is the modified equation. The $r-$th truncation of the modified Lagrangian is denoted by $\bar {\mathfrak L}_{\rm mod,r}$,

\begin{equation}
\bar{\mathfrak L}_{\rm mod,r}({\q}, \q',[t],\Delta a)=\tau_r\left( \bar{\mathfrak L}_{\rm mod}( \q,\q',[t],\Delta a) \right)=\tau_r\left(\bar{\mathfrak L}_{\rm mesh}([\qbar],\Delta a)\bigg|_{\q^{(r)}=F^{(j)}_{r-2}( \q, \q',[t],\Delta a)}  \right).
\end{equation}

From the definition of $r-$th truncation of  modified Lagrangian, it follows that $\bar {\mathfrak L}_{\rm mod,r}({\q},{\q}',[t], \Delta a)=\bar {\mathfrak L}_{\rm mesh}( [\qbar] , \Delta a) + O(\Delta a^{r+1})$ for meshed $r-$critical families of curves. In fact, Vermeeren [20] showed that the meshed modified Lagrangian $\bar {\mathfrak L}_{\rm mesh}$ and truncated modified Lagrangian $\bar{\mathfrak L}_{\rm mod,r}$ have the same $r-$critical families of curves.

For a discretized Lagrangian $\bar{L}_{\rm disc}$ that is a consistent discretization of a regular Lagrangian $\bar L(\qbar(a),\qbar'(a))$, the $r-$th truncation of the Euler-Lagrange equation of $\bar{\mathfrak L}_{\rm mod,r}({\q}(a),{\q}'(a),[t(a)],\Delta a)$, solved for ${\q}''(a)$, is the $r-$th truncation of the modified equation.

This result follows directly from applying backward stability theorem from [20] to the transformed Lagrangian system defined in Section <ref>.

Although our theorem is based on the application of ideas from [20] to the transformed setting, the final result for the energy-preserving algorithm has two key differences due to the redundant nature of the energy equation.

* The state vector $\qbar$ in the transformed setting contains both $\qbar$ and $t$, the modified equation is only obtained and solved for $\q''$. This is due to the fact that the modified equation corresponding to the $t$ variable is redundant.

* The modified Lagrangian $\bar{\mathfrak L}_{\rm mod,r}$ in the transformed setting depends on $\q,\q',\Delta a$ and arbitrary high derivatives of $t$. This again is related to the redundant nature of the energy equation which prevents us from solving the modified equations jointly for $\qbar''$.

§.§ Mechanical System Example

We consider a one degree of freedom mechanical system with a separable Lagrangian to illustrate the backward stability result shown in Section <ref>. Consider a Lagrangian system $L(q,\dot{q})=\frac{1}{2}m\dot{q}^2 - V(q)$ where $V(q)$ is the potential energy of the mechanical system. We employ the time transformation (<ref>) to obtain the transformed Lagrangian

\begin{equation}
\bar L(\bar q(a), \bar q'(a))=t'(a)L\left(q(a), \frac{q'(a)}{t'(a)}\right)=t'(a) \left( \frac{1}{2}m\left( \frac{q'(a)}{t'(a)} \right)^2 - V(q(a)) \right).
\end{equation}

We use the midpoint rule to obtain the discretized Lagrangian $\bar L_{disc}$

\begin{equation}
\bar L_{\rm disc}(\bar q_k, \bar q_{k+1},\Delta a)=\frac{t_{k+1}-t_k}{\Delta a} \left( \frac{1}{2}m\left( \frac{q_{k+1}-q_k}{t_{k+1}-t_k} \right)^2 - V\left(\frac{q_{k+1}+q_k}{2}\right) \right).
\end{equation}

Using the power series expansions for $q_{k\pm1}(a)=q(a) \pm \frac{\Delta a}{2}q'(a) + \frac{\Delta a^2}{8}q''(a) \pm \frac{\Delta a^3}{24}q'''(a) + O(\Delta a^4) $ and $t_{k\pm1}(a)=t(a) \pm \frac{\Delta a}{2}t'(a) + \frac{\Delta a^2}{8}t''(a) \pm \frac{\Delta a^3}{24}t'''(a) + O(\Delta a^4) $, we obtain

\begin{align}\nonumber
\bar{\mathfrak L}_{\rm disc}([\bar q],\Delta a)&=\left( t' + \frac{\Delta a^2}{24}t''' + \cdots \right) \left( \frac{1}{2}m\left( \frac{ q' + \frac{\Delta a^2}{24}q''' + \cdots }{t' + \frac{\Delta a^2}{24}t''' + \cdots } \right)^2  - V\left(q + \frac{\Delta a^2}{8}q'' + \cdots \right) \right) \\
%=\left( t' + \frac{\Delta a^2}{24}t''' + O(\Delta a^4)\right)  \left( \frac{1}{2}m\left( \frac{ q' + \frac{\Delta a^2}{24}q''' + O(\Delta a^4)}{t' + \frac{\Delta a^2}{24}t''' + O(\Delta a^4)} \right)^2  - \left(V(q) + \frac{\Delta a^2}{8}V'(q)q''(a) + O(\Delta a^4)\right) \right) \\
&=t'\left( \frac{1}{2}m \left( \frac{q'}{t'} \right)^2 -V \right)  + \frac{\Delta a^2}{24} \left( \frac{mq'q'''}{t'}- \frac{m(q')^2t'''}{2(t')^2} - 3q'' t'V_q - t'''V \right) + O(\Delta a^4),
\end{align}

where $V$ with no arguments is $V(q(a))$ and $V_q=\frac{\partial V}{\partial q}(q(a))$. From $\bar{\mathfrak{L}}_{\rm disc}([\bar q],\Delta a)$, we calculate the meshed modified Lagrangian in the transformed setting

\begin{align}\nonumber
\bar{\mathfrak L}_{\rm mesh}([\bar q],\Delta a)&= \bar{\mathfrak L}_{\rm disc}([\bar q],\Delta a) - \frac{\Delta a^2}{24} \frac{d^2}{da^2} \bar{\mathfrak L}_{\rm disc}([\bar q],\Delta a)  + O(\Delta a^4) \\
%=t'(a)\left( \frac{1}{2}m \left( \frac{q'(a)}{t'(a)} \right)^2 -V(q) \right)  + \frac{\Delta a^2}{24} \left(  \frac{mq'(a)q'''(a)}{t'(a)} -  \frac{mq'(a)^2t'''(a)}{2t'(a)^2} - 3V_q(q)q''(a) t'(a) - V(q)t'''(a) \right) \\
%- \frac{\Delta a^2}{24} \frac{d^2}{da^2} \left( t'(a)\left( \frac{1}{2}m \left( \frac{q'(a)}{t'(a)} \right)^2 -V(q) \right) \right) + O(\Delta a^4)\\
%=t'(a)\left( \frac{1}{2}m \left( \frac{q'(a)}{t'(a)} \right)^2 -V(q) \right)  + \frac{\Delta a^2}{24} \left(  \frac{mq'(a)q'''(a)}{t'(a)} -  \frac{mq'(a)^2t'''(a)}{2t'(a)^2} - 3V_q(q)q''(a) t'(a) - V(q)t'''(a) \right) \\
%- \frac{\Delta a^2}{24}  \left( \frac{mq''^2}{t'} +  \frac{mq'q'''}{t'} - \frac{2mq'q''t''}{t'^2}   - \frac{mq'^2t'''}{2t'^2}+ \frac{mq'^2t''^2}{t'^3}- t'''V - 2t''V_qq' - t'V_{qq}q'^2 - t'V_qq''\right)  + O(\Delta a^4) \\
&=t'\left( \frac{1}{2}m \left( \frac{q'}{t'} \right)^2 -V \right)  + \frac{\Delta a^2}{24} \left(-\frac{m(q'')^2}{t'}  + \frac{2mq'q''t''}{(t')^2} - \frac{m(q')^2(t'')^2}{(t')^3} \right) \nonumber \\
&+ \frac{\Delta a^2}{24} \left( 2q't''V_q + (q')^2t'V_{qq} - 2q'' t'V_q \right) + O(\Delta a^4).
\label{eq:Lmesh}
\end{align}

The modified Euler-Lagrange and modified energy equation up to second order in the transformed setting are obtained from

\begin{align}\nonumber
O(\Delta a^4)&= \frac{\partial \bar {\mathfrak L}_{\rm mesh}}{\partial q} - \frac{d}{da}\left(\frac{\partial \bar{\mathfrak L}_{\rm mesh}}{\partial q'}\right) \\
%&= \left( -t'V' + \frac{\Delta a^2}{24} \left( 2V'' q' t'' + V^{(3)}(q')^2t' - 2V''q''t'\right) \right)\\
%&- \frac{d}{da} \left( \frac{mq'}{t'} + \frac{\Delta a^2}{24} \left(  \frac{2mq''t''}{(t')^2} - \frac{2mq'(t'')^2}{(t')^3} + 2V't'' + 2V''q't' \right) \right) \\
%&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V' + \frac{\Delta a^2}{24} \left( 2V'' q' t'' + V^{(3)}(q')^2t' - 2V''q''t' \right) \\
%&-   \frac{\Delta a^2}{24} \left(  \frac{2mq'''t''}{(t')^2} +  \frac{2mq''t'''}{(t')^2} -  \frac{4mq''(t'')^2}{(t')^3} - \frac{2mq''(t'')^2}{(t')^3} - \frac{4mq't''t'''}{(t')^3}+ \frac{6mq'(t'')^3}{(t')^4} \right) \\
%&-   \frac{\Delta a^2}{24} \left( 2V''q't'' + 2V't''' + 2V'''(q')^2t' + 2V''q''t'  + 2V''q't'' \right) \\
&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V_q + \frac{\Delta a^2}{24} \left(- (q')^2t'V_{qqq}-2t'''V_q -2q't''V_{qq}   - 4q''t'V_{qq} \right) \nonumber \\
&-   \frac{\Delta a^2}{24} \left(  \frac{2mq'''t''}{(t')^2} +  \frac{2mq''t'''}{(t')^2} -  \frac{6mq''(t'')^2}{(t')^3} - \frac{4mq't''t'''}{(t')^3}+ \frac{6mq'(t'')^3}{(t')^4} \right), \\\nonumber
O(\Delta a^4)&= \frac{\partial \bar {\mathfrak L}_{mesh}}{\partial t} - \frac{d}{da}\left(\frac{\partial \bar{\mathfrak L}_{mesh}}{\partial t'}\right)\\
&=\left(\frac{q'}{t'} \right)\left(-\frac{mq''}{t'} + \frac{mq't''}{(t')^2} - t'V_q \right) + O(\Delta a^2),
\end{align}

where $V_{qq}=\frac{\partial^2 V}{\partial q^2}(q(a))$ and $V_{qqq}=\frac{\partial^3 V}{\partial q^3}(q(a))$.
Due to the redundant nature of the leading order in modified energy equation, the modified equations for the transformed system can not be solved recursively for $\bar q''=[q''(a), t''(a)]$. This is a consequence of the redundancy in the Euler-Lagrange equations in the transformed setting. We solve recursively for $q''$ in terms of $(q,q')$ and $(t',t'',t''')$. In the leading order we have $q''=\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)$. Substituting this in the second-order terms in the modified Euler-Lagrange equation, we get

\begin{align} \nonumber
&= \frac{\partial \bar {\mathfrak L}_{\rm mesh}}{\partial q} - \frac{d}{da}\left(\frac{\partial \bar{\mathfrak L}_{\rm mesh}}{\partial q'}\right) + O(\Delta a^4) \\ \nonumber
&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V_q + \frac{\Delta a^2}{24} \left(- (q')^2t'V_{qqq}-2t'''V_q -2q't''V_{qq}   - 4\left( \frac{q't''}{t'} - \frac{(t')^2V_q}{m} \right)t'V_{qq} \right) \\ \nonumber
&-   \frac{\Delta a^2}{24} \left(  \frac{2m\frac{d}{da}\left( \frac{q't''}{t'} - \frac{(t')^2V_q}{m} \right)t''}{(t')^2} +  \frac{2m\left( \frac{q't''}{t'} - \frac{(t')^2V_q}{m} \right)t'''}{(t')^2} -  \frac{6m\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)(t'')^2}{(t')^3} \right) \\ \nonumber
&+   \frac{\Delta a^2}{24} \left(  \frac{4mq't''t'''}{(t')^3}- \frac{6mq'(t'')^3}{(t')^4} \right) + O(\Delta a^4)  \\
%&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V' + \frac{\Delta a^2}{24} \left(- V'''(q')^2t'-2V't''' -2V''q't''   - 4V''\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)t' \right) \\
%&-   \frac{\Delta a^2}{24} \left(  \frac{2m\left( \frac{q''t''}{t'}+ \frac{q't'''}{t'}-\frac{q't''^2}{t'^2} - \frac{2t't''V'}{m}-\frac{t'^2V''q'}{m} \right)t''}{(t')^2} -  2 V' t''' +  \frac{6 V'(t'')^2}{t'} - \frac{2mq't''t'''}{(t')^3} \right) \\
%&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V' + \frac{\Delta a^2}{24} \left(- V'''(q')^2t' -4V''q't'' +  \frac{4V''t'^3V'}{m}  \right) \\
%&-   \frac{\Delta a^2}{24} \left(  \frac{2m\left( \frac{\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)t''}{t'}-\frac{q't''^2}{t'^2} \right)t''}{(t')^2} +  \frac{2 V'(t'')^2}{t'}  \right) \\
&= -\frac{mq''}{t'}+\frac{mq't''}{(t')^2}  -t'V_q + \frac{\Delta a^2}{24} \left( \frac{4(t')^3V_qV_{qq}}{m} -4q't''V_{qq}- (q')^2t'V_{qqq}  \right) + O(\Delta a^4),
\end{align}

which can be re-written as

\begin{equation}
q''= \frac{q't''}{t'}  -\frac{(t')^2V_q}{m} + \frac{\Delta a^2}{24m} \left( \frac{4(t')^4V_qV_{qq}}{m} -4q't't''V_{qq}- (q')^2(t')^2V_{qqq}\right) + O(\Delta a^4).
\label{eq:EL_mesh}
\end{equation}

We also obtain the modified Lagrangian up to second order

\begin{align} \nonumber
\bar{\mathfrak L}_{\rm mod}& =\bar{\mathfrak L}_{\rm mesh}([\bar q],\Delta a)\bigg|_{q''=\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)} \\ \nonumber
&=t'\left( \frac{1}{2}m \left( \frac{q'}{t'} \right)^2 -V \right)  + \frac{\Delta a^2}{24} \left( \frac{2mq'\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)t''}{(t')^2}-\frac{m\left( \frac{q't''}{t'} - \frac{(t')^2V_q}{m} \right)^2}{t'}   - \frac{m(q')^2(t'')^2}{(t')^3} \right) \\ \nonumber
&+ \frac{\Delta a^2}{24} \left( 2q't''V_q + (q')^2t'V_{qq} - 2\left( \frac{q't''}{t'} - \frac{(t')^2V_q}{m} \right) t'V_q \right) + O(\Delta a^4)\\
%=t'(a)\left( \frac{1}{2}m \left( \frac{q'(a)}{t'(a)} \right)^2 -V(q) \right)  + \frac{\Delta a^2}{24} \left(-\frac{m\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)^2}{t'(a)}  + \frac{2mq'(a)\left( \frac{q't''}{t'} - \frac{t'^2V_q}{m} \right)t''(a)}{t'(a)^2} - \frac{mq'(a)^2t''(a)^2}{t'(a)^3} \right) \\
%+ \frac{\Delta a^2}{24} \left( V_{qq}(q)q'(a)^2t'(a) + \frac{2t'^3V^2_q}{m} \right) + O(\Delta a^4)\\
&=t'\left( \frac{1}{2}m \left( \frac{q'}{t'} \right)^2 -V \right)  + \frac{\Delta a^2}{24} \left( \frac{t'^3(V_q)^2}{m} + q'^2t'V_{qq} \right) +  O(\Delta a^4).
\end{align}

Unlike the meshed modified Lagrangian (<ref>) which depends on higher derivatives of $q$ and $t$, the modified Lagrangian in this case only depends on $\bar{q}=(q,t)$ and first derivative $\bar{q}'=(q',t')$.
The Euler-Lagrange equation of $\bar{\mathfrak L}_{\rm mod,3}$ is

\begin{align}\nonumber
0&= \frac{\partial \bar {\mathfrak L}_{\rm mod,3}}{\partial q} - \frac{d}{da}\left(\frac{\partial \bar{\mathfrak L}_{\rm mod,3}}{\partial q'}\right)\\\nonumber
&=-t'V_q + \frac{\Delta a^2}{24} \left( \frac{2t'^3V_qV_{qq}}{m} + V_{qqq}q'^2t' \right) - \frac{d}{da}\left( \frac{mq'}{t'} + \frac{\Delta a^2}{24} \left( 2V_{qq}q't' \right) \right) \\
&=-\frac{mq''}{t'}+\frac{mq't''}{(t')^2} -t'V_q + \frac{\Delta a^2}{24} \left( \frac{2t'^3V_qV_{qq}}{m}   - 2q'V_{qq} t'' - V_{qqq}q'^2t' -2V''q''t'  \right).
\label{eq:EL_mod}
\end{align}

It is important to note that the Euler-Lagrange equation for $\bar{\mathfrak L}_{\rm mod,3}$ does not contain an error term. However, when we solve the Euler-Lagrange equation (<ref>) for $q''$

\begin{align}\nonumber
%q''\left( \frac{m}{t'} + \frac{\Delta a^2 2V''t'}{24} \right)& = \frac{mq't''}{(t')^2} -t'V' + \frac{\Delta a^2}{24} \left( \frac{2t'^3V'V''}{m}   - 2q'V''t'' - V'''q'^2t' + 2V't''' \right)\\
q''&= \frac{\frac{mq't''}{(t')^2} -t'V_q + \frac{\Delta a^2}{24} \left( \frac{2(t')^3V_qV_{qq}}{m}   - 2q't''V_{qq} - q'^2t'V_{qqq} \right)}{\left( \frac{m}{t'} + \frac{\Delta a^2 t'V_{qq}}{12} \right)}\\\nonumber
%&= \frac{\frac{q't''}{t'} -\frac{t'^2V'}{m} + \frac{\Delta a^2}{24} \left( \frac{2t'^4V'V''}{m^2}   - \frac{2q'V''t''t'}{m} - \frac{V'''q'^2t'^2}{m}\right)}{\left( 1 + \frac{\Delta a^2 2V''t'^2}{24m} \right)}\\
&=\left[ \frac{q't''}{t'} -\frac{(t')^2V_q}{m} + \frac{\Delta a^2}{24m} \left( \frac{2(t')^4V_qV_{qq}}{m}   - 2q't't''V_{qq} - (q't')^2V_{qqq} \right) \right]\left( 1 - \Delta a^2  \frac{(t')^2V_{qq}}{12m} + O(\Delta a^4) \right)\\
%&=\frac{q't''}{t'} -\frac{t'^2V'}{m} + \frac{\Delta a^2}{24m} \left( \frac{2t'^4V'V''}{m}   - 2q'V''t''t' - V'''q'^2t'^2  \right) -  \frac{\Delta a^2}{24m} \left( 2V''q't''t' - \frac{2t'^4V'V''}{m}  \right)\\
&=\frac{q't''}{t'} -\frac{(t')^2V_q}{m} + \frac{\Delta a^2}{24m} \left( \frac{4(t')^4V_qV_{qq}}{m}   - 4q't't''V_{qq} - (q')^2(t')^2V_{qqq} \right) + O(\Delta a^4),
\end{align}

we again get the modified equation (<ref>). Thus, the modified equation for the variational integrator based on the midpoint rule, solved for $q''$, is $O(\Delta a^4)$ close to the Euler-Lagrange equation of the truncated modified Lagrangian $\bar{\mathfrak L}_{\rm mod,3}$.

§ CONCLUSIONS

We have presented a numerical and theoretical assessment of energy-preserving, adaptive time-step variational integrators. We first compared the time adaptation and energy performance of the energy-preserving adaptive algorithm with the adaptive variational integrators for Kepler's two-body problem. We then obtained a conservative bound on the discrete energy error using the discrete energy equation and implemented the energy-preserving adaptive algorithm using VPA to demonstrate exact discrete energy conservation up to machine precision.
Motivated by these numerical results, we proposed a time transformation approach to investigate the numerical stability of the adaptive algorithm. Using the equivalence of trajectories, we derived the energy-preserving, adaptive time-step variational integrators as fixed time-step variational integrators in the transformed setting. We then applied results from [20] in the transformed setting to show that the $r-$th truncation of the modified equation for the variational integrator in the transformed setting, solved for $\q''$, is $O(\Delta a^{r+1})$ close to the Euler-Lagrange equation of the $r-$truncated modified Lagrangian $\bar{\mathfrak L}_{\rm mod,r}$. Finally, we constructed a modified Lagrangian for a midpoint rule discretization of a simple mechanical system.
For future work, we would like to extend the backward stability results from this paper to the energy-preserving, adaptive time-step Lie group variational integrators. We would also like to study the connections between modified equations of adaptive variational integrators based on a prescribed monitor function and energy-preserving, adaptive time-step variational integrators.

Declaration of interest: None.

Funding: Harsh Sharma, Mayuresh Patil and Craig Woolsey were supported in part by the National Science Foundation under Grant No. 1826152. Jeff Borggaard was supported in part by the National Science Foundation under grant DMS-1819110.

[1]
Lew AJ, Mata P.
A brief introduction to variational integrators.
In: Structure-preserving integrators in nonlinear structural dynamics
and flexible multibody dynamics. Springer; 2016. p. 201–91.

[2]
Sharma H, Patil M, Woolsey C.
A review of structure-preserving numerical methods for engineering
Computer Methods in Applied Mechanics and Engineering.

[3]
Marsden JE, West M.
Discrete mechanics and variational integrators.
Acta Numer. 2001;10:357–514.

[4]
Wendlandt JM, Marsden JE.
Mechanical integrators derived from a discrete variational principle.
Physica D: Nonlinear Phenomena. 1997;106(3-4):223–246.

[5]
Leok M.
Generalized Galerkin variational integrators; 2005.
arXiv:math/0508360 [math.NA].

[6]
Kane C, Marsden JE, Ortiz M, West M.
Variational integrators and the Newmark algorithm for conservative
and dissipative mechanical systems.
International Journal for numerical methods in engineering.

[7]
Kane C, Marsden JE, Ortiz M.
Symplectic-energy-momentum preserving variational integrators.
J Math Phys. 1999;40(7):3353–71.

[8]
Guibout VM, Bloch AM.
Discrete variational principles and Hamilton-Jacobi theory for
mechanical systems and optimal control problems.
arXiv preprint math/0409296. 2004;.

[9]
Sharma H, Patil M, Woolsey C.
Energy-preserving variational integrators for forced Lagrangian
Communications in Nonlinear Science and Numerical Simulation.

[10]
Sharma H, Lee T.
Energy-preserving, adaptive time-step Lie group variational
integrators for the attitude dynamics of a rigid body.
In: 2019 American Control Conference (ACC). IEEE; 2019. p.

[11]
Sharma H, Patil M, Woolsey C.
Energy-preserving, adaptive time-step Lie group variational
integrators for rigid body motion in SE (3).
In: 2019 IEEE 58th Conference on Decision and Control (CDC). IEEE;
2019. p. 8079–8084.

[12]
Shibberu Y.
Is symplectic-energy-momentum integration well-posed?; 2006.

[13]
Shibberu Y.
How to regularize a symplectic-energy-momentum integrator.
arXiv preprint math/0507483. 2005;.

[14]
Wilkinson JH.
Error analysis of floating-point computation.
Numerische Mathematik. 1960;2(1):319–340.

[15]
Eirola T.
Aspects of backward error analysis of numerical ODEs.
Journal of Computational and Applied Mathematics.

[16]
Hairer E, Lubich C.
The life-span of backward error analysis for numerical integrators.
Numerische Mathematik. 1997;76(4):441–462.

[17]
Sanz-Serna J.
Backward error analysis of symplectic integrators.
Fields Inst Comm. 1996;10:193–205.

[18]
Reich S.
Backward error analysis for numerical integrators.
SIAM Journal on Numerical Analysis. 1999;36(5):1549–1570.

[19]
Hairer E.
Variable time step integration with symplectic methods.
Applied Numerical Mathematics. 1997;25(2-3):219–227.

[20]
Vermeeren M.
Modified equations for variational integrators.
Numerische Mathematik. 2017;137(4):1001–1037.

[21]
Marsden JE, Pekarsky S, Shkoller S, West M.
Variational methods, multisymplectic geometry and continuum
J Geom Phys. 2001;38(3):253–84.

[22]
Goldstein H.
Classical Mechanics.
Reading, MA: Addison-Wesley; 1980. 672 pages.

[23]
Schmitt J, Leok M.
Adaptive variational integrators.
arXiv preprint arXiv:170901975. 2017;.

[24]
Hairer E, Lubich C, Wanner G.
Geometric Numerical Integration: Structure-Preserving Algorithms for
Ordinary Differential Equations. vol. 31.
Springer Science & Business Media; 2006.