Time-adaptive Lagrangian Variational Integrators
for Accelerated Optimization on Manifolds

Given a $n$-dimensional manifold $\mathcal{Q}$, a Lagrangian is a function $L:T\mathcal{Q}\rightarrow\mathbb{R}$. The corresponding action integral $\mathcal{S}$ is the functional

$$\mathcal{S}(q)=\int_{0}^{T}{L(q,\dot{q})dt},$$

(2.1)

over the space of smooth curves $q:[0,T]\rightarrow\mathcal{Q}$. Hamilton’s variational principle states that $\delta\mathcal{S}=0$ where the variation $\delta\mathcal{S}$ is induced by an infinitesimal variation $\delta q$ of the trajectory $q$ that vanishes at the endpoints. Given local coordinates $(q^{1},\ldots,q^{n})$ on the manifold $\mathcal{Q}$, Hamilton’s variational principle can be shown to be equivalent to the Euler–Lagrange equations,

$$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{q}^{k}}\right)=\frac{\partial L}{\partial q^{k}},\qquad\text{for }k=1,\ldots,n.$$

(2.2)

The Legendre transform $\mathbb{F}L:T\mathcal{Q}\rightarrow T^{*}\mathcal{Q}$ of $L$ is defined fiberwise by $\mathbb{F}L:(q^{i},\dot{q}^{i})\mapsto\left(q^{i},\frac{\partial L}{\partial\dot{q}^{i}}\right)$, and we say that a Lagrangian $L$ is regular or nondegenerate if the Hessian matrix $\frac{\partial^{2}L}{\partial\dot{q}^{2}}$ is invertible for every $q$ and $\dot{q}$, and hyperregular if the Legendre transform $\mathbb{F}L$ is a diffeomorphism. A hyperregular Lagrangian on $T\mathcal{Q}$ induces a Hamiltonian system on $T^{*}\mathcal{Q}$ via

$$H(q,p)=\langle\mathbb{F}L(q,\dot{q}),\dot{q}\rangle-L(q,\dot{q})=\sum_{j=1}^{n}{p_{j}\dot{q}^{j}}-L(q,\dot{q})\bigg{|}_{p_{i}=\frac{\partial L}{\partial\dot{q}^{i}}},$$

(2.3)

where $p_{i}=\frac{\partial L}{\partial\dot{q}^{i}}\in T^{*}\mathcal{Q}$ is the conjugate momentum of $q^{i}$. A Hamiltonian $H$ is called hyperregular if $\mathbb{F}H:T^{*}\mathcal{Q}\rightarrow T\mathcal{Q}$ defined by $\mathbb{F}H(\alpha)\cdot\beta=\frac{d}{ds}\big{|}_{s=0}H(\alpha+s\beta)$, is a diffeomorphism. Hyperregularity of the Hamiltonian $H$ implies invertibility of the Hessian matrix $\frac{\partial^{2}H}{\partial p^{2}}$ and thus nondegeneracy of $H$. Theorem 7.4.3 in [48] states that hyperregular Lagrangians and hyperregular Hamiltonians correspond in a bijective manner. We can also define a Hamiltonian variational principle on the Hamiltonian side in momentum phase space which is equivalent to Hamilton’s equations,

$$\dot{p}_{k}=-\frac{\partial H}{\partial q^{k}}(p,q),\qquad\dot{q}^{k}=\frac{\partial H}{\partial p_{k}}(p,q),\qquad\text{for }k=1,\ldots,n.$$

(2.4)

These equations can also be shown to be equivalent to the Euler–Lagrange equations (2.2), provided that the Lagrangian is hyperregular.

Variational integrators are derived by discretizing Hamilton’s principle, instead of discretizing Hamilton’s equations directly. As a result, variational integrators are symplectic, preserve many invariants and momentum maps, and have excellent long-time near-energy preservation (see [49]).

Traditionally, variational integrators have been designed based on the Type I generating function commonly known as the discrete Lagrangian, $L_{d}:Q\times Q\rightarrow\mathbb{R}$. The exact discrete Lagrangian that generates the time-$h$ flow of Hamilton’s equations can be represented both in a variational form and in a boundary-value form. The latter is given by

$$L_{d}^{E}(q_{0},q_{1};h)=\int_{0}^{h}L(q(t),\dot{q}(t))dt,$$

(2.5)

where $q(0)=q_{0},$ $q(h)=q_{1},$ and $q$ satisfies the Euler–Lagrange equations over the time interval $[0,h]$. A variational integrator is defined by constructing an approximation $L_{d}:Q\times Q\rightarrow\mathbb{R}$ to $L_{d}^{E}$, and then applying the discrete Euler–Lagrange equations,

$$p_{k}=-D_{1}L_{d}(q_{k},q_{k+1}),\qquad p_{k+1}=D_{2}L_{d}(q_{k},q_{k+1}),$$

(2.6)

where $D_{i}$ denotes a partial derivative with respect to the $i$-th argument, and these equations implicitly define the integrator map $\tilde{F}_{L_{d}}:(q_{k},p_{k})\mapsto(q_{k+1},p_{k+1})$. The error analysis is greatly simplified via Theorem 2.3.1 of [49], which 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}$ to order $r$, i.e.,

$$L_{d}(q_{0},q_{1};h)=L_{d}^{E}(q_{0},q_{1};h)+\mathcal{O}(h^{r+1}),$$

(2.7)

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, has order of accuracy $r$. Many other properties of the integrator, such as momentum conservation properties of the method, can be determined by analyzing the associated discrete Lagrangian, as opposed to analyzing the integrator directly.

Variational integrators have been extended to the framework of Type II and Type III generating functions, commonly referred to as discrete Hamiltonians (see [34; 46; 57]). Hamiltonian variational integrators are derived by discretizing Hamilton’s phase space principle. The boundary-value formulation of the exact Type II generating function of the time-$h$ flow of Hamilton’s equations is given by the exact discrete right Hamiltonian,

$$H_{d}^{+,E}(q_{0},p_{1};h)=p_{1}^{\top}q_{1}-\int_{0}^{h}\left[p(t)^{\top}\dot{q}(t)-H(q(t),p(t))\right]dt,$$

(2.8)

where $(q,p)$ satisfies Hamilton’s equations with boundary conditions $q(0)=q_{0}$ and $p(h)=p_{1}$. A Type II Hamiltonian variational integrator is constructed by using a discrete right Hamiltonian $H_{d}^{+}$ approximating $H_{d}^{+,E}$, and applying the discrete right Hamilton’s equations,

$$p_{0}=D_{1}H_{d}^{+}(q_{0},p_{1}),\qquad q_{1}=D_{2}H_{d}^{+}(q_{0},p_{1}),$$

(2.9)

which implicitly defines the integrator, $\tilde{F}_{H_{d}^{+}}:(q_{0},p_{0})\mapsto(q_{1},p_{1})$.

Theorem 2.3.1 of [49], which simplified the error analysis for Lagrangian variational integrators, has an analogue for Hamiltonian variational integrators. Theorem 2.2 in [57] states that if a discrete right Hamiltonian $H^{+}_{d}$ approximates the exact discrete right Hamiltonian $H_{d}^{+,E}$ to order $r$, i.e.,

$$H^{+}_{d}(q_{0},p_{1};h)=H_{d}^{+,E}(q_{0},p_{1};h)+\mathcal{O}(h^{r+1}),$$

(2.10)

then the discrete right Hamilton’s map $\tilde{F}_{H^{+}_{d}}:(q_{k},p_{k})\mapsto(q_{k+1},p_{k+1})$, viewed as a one-step method, is order $r$ accurate. Note that discrete left Hamiltonians and corresponding discrete left Hamilton’s maps can also be constructed in the Type III case (see [46; 20]).

Examples of variational integrators include Galerkin variational integrators  [46; 49], Prolongation-Collocation variational integrators [45], and Taylor variational integrators [58]. In many cases, the Type I and Type II/III approaches will produce equivalent integrators. This equivalence has been established in [58] for Taylor variational integrators provided the Lagrangian is hyperregular, and in [46] for generalized Galerkin variational integrators constructed using the same choices of basis functions and numerical quadrature formula provided the Hamiltonian is hyperregular. However, Hamiltonian and Lagrangian variational integrators are not always equivalent. In particular, it was shown in [57] that even when the Hamiltonian and Lagrangian integrators are analytically equivalent, they might still have different numerical properties because of numerical conditioning issues. Even more to the point, Lagrangian variational integrators cannot always be constructed when the underlying Hamiltonian is degenerate. This is particularly relevant in variational accelerated optimization since the time-adaptive Hamiltonian framework for accelerated optimization presented in [20] relies on a degenerate Hamiltonian which has no associated Lagrangian description. We will thus not be able to exploit the usual correspondence between Hamiltonian and Lagrangian dynamics and will have to come up with a different strategy to allow time-adaptivity on the Lagrangian side.

We now describe the construction of Taylor variational integrators as introduced in [58] as we will use them in our numerical experiments. A discrete approximate Lagrangian or Hamiltonian is constructed by approximating the flow map and the trajectory associated with the boundary values using a Taylor method, and approximating the integral by a quadrature rule. The Taylor variational integrator is generated by the implicit discrete Euler–Lagrange equations associated with the discrete Lagrangian or by the Hamilton’s equations associated with the discrete Hamiltonian. More explicitly, we first construct $\rho$-order and $(\rho+1)$-order Taylor methods $\Psi_{h}^{(\rho)}$ and $\Psi_{h}^{(\rho+1)}$ approximating the exact time-$h$ flow map $\Phi_{h}:TQ\rightarrow TQ$ corresponding to the Euler–Lagrange equation in the Type I case or the exact time-$h$ flow map $\Phi_{h}:T^{*}Q\rightarrow T^{*}Q$ corresponding to Hamilton’s equation in the Type II case. Let $\pi_{Q}:(q,p)\mapsto q$ and $\pi_{T^{*}Q}:(q,p)\mapsto p$. Given a quadrature rule of order $s$ with weights and nodes $(b_{i},c_{i})$ for $i=1,...,m$, the Taylor variational integrators are then constructed as follows:

1. (i)

   Approximate $p(0)=p_{0}$ by the solution $\tilde{p}_{0}$ of the problem $p_{1}=\pi_{T^{*}Q}\circ\Psi_{h}^{(\rho)}(q_{0},\tilde{p}_{0}).$

2. (ii)

   Generate approximations $(q_{c_{i}},p_{c_{i}})\approx(q(c_{i}h),p(c_{i}h))$ via $(q_{c_{i}},p_{c_{i}})=\Psi_{c_{i}h}^{(\rho)}(q_{0},\tilde{p}_{0}).$

3. (iii)

   Approximate $q_{1}$ via $\tilde{q}_{1}=\pi_{Q}\circ\Psi_{h}^{(\rho+1)}(q_{0},\tilde{p}_{0}).$

4. (iv)

   Use the continuous Legendre transform to obtain $\dot{q}_{c_{i}}=\frac{\partial H}{\partial p_{c_{i}}}$.

5. (v)

   Apply the quadrature rule to obtain the associated discrete right Hamiltonian,

   $$H_{d}^{+}(q_{0},p_{1};h)=p_{1}^{\top}\tilde{q}_{1}-h\sum_{i=1}^{m}{b_{i}\left[p_{c_{i}}^{\top}\dot{q}_{c_{i}}-H(q_{c_{i}},p_{c_{i}})\right]}.$$

6. (vi)

   The variational integrator is then defined by the implicit discrete right Hamilton’s equations,

   $$q_{1}=D_{2}H_{d}^{+}(q_{0},p_{1}),\qquad p_{0}=D_{1}H_{d}^{+}(q_{0},p_{1}).$$

Note that analogous constructions and error analysis results have been derived in [20; 58] for discrete left Hamiltonians in the Type III case.

Symplectic integrators form a class of geometric numerical integrators of interest since, when applied to conservative Hamiltonian systems, they yield discrete approximations of the flow that preserve the symplectic 2-form [26], which results in the preservation of many qualitative aspects of the underlying system. In particular, symplectic integrators exhibit excellent long-time near-energy preservation. However, when symplectic integrators were first used in combination with variable time-steps, the near-energy preservation was lost and the integrators performed poorly [7; 24]. Backward error analysis provided justification both for the excellent long-time near-energy preservation of symplectic integrators and for the poor performance experienced when using variable time-steps (see Chapter IX of [26]): it shows that symplectic integrators can be associated with a modified Hamiltonian in the form of a formal power series in terms of the time-step. The use of a variable time-step results in a different modified Hamiltonian at every iteration, which is the source of the poor energy conservation. The Poincaré transformation is one way to incorporate variable time-steps in geometric integrators without losing the nice conservation properties associated with these integrators.

Given a Hamiltonian $H(q,t,p)$, consider a desired transformation of time $t\mapsto\tau$ described by the monitor function $g(q,t,p)$ via

$$\frac{dt}{d\tau}=g(q,t,p).$$

(2.11)

The time $t$ shall be referred to as the physical time, while $\tau$ will be referred to as the fictive time, and we will denote derivatives with respect to $t$ and $\tau$ by dots and apostrophes, respectively. A new Hamiltonian system is constructed using the Poincaré transformation,

$$\bar{H}(\bar{q},\bar{p})=g(q,\mathfrak{q},p)\left(H(q,\mathfrak{q},p)+\mathfrak{p}\right),$$

(2.12)

in the extended phase space defined by $\bar{q}=\left[\begin{smallmatrix}q\\ \mathfrak{q}\end{smallmatrix}\right]\in\bar{\mathcal{Q}}$ and $\bar{p}=\left[\begin{smallmatrix}p\\ \mathfrak{p}\end{smallmatrix}\right]$ where $\mathfrak{p}$ is the conjugate momentum for $\mathfrak{q}=t$ with $\mathfrak{p}(0)=-H(q(0),0,p(0))$. The corresponding equations of motion in the extended phase space are then given by

$$\bar{q}^{\prime}=\frac{\partial\bar{H}}{\partial\bar{p}},\qquad\qquad\bar{p}^{\prime}=-\frac{\partial\bar{H}}{\partial\bar{q}}.$$

(2.13)

Suppose $(\bar{Q}(\tau),\bar{P}(\tau))$ are solutions to these extended equations of motion, and let $(q(t),p(t))$ solve Hamilton’s equations for the original Hamiltonian $H$. Then

$$\bar{H}(\bar{Q}(\tau),\bar{P}(\tau))=\bar{H}(\bar{Q}(0),\bar{P}(0))=0.$$

(2.14)

Therefore, the components $(Q(\tau),P(\tau))$ in the original phase space of the augmented solutions $(\bar{Q}(\tau),\bar{P}(\tau))$ satisfy

$$H(Q(\tau),\tau,P(\tau))=-\mathfrak{p}(\tau),\qquad H(Q(0),0,P(0))=-\mathfrak{p}(0)=H(q(0),0,p(0)).$$

(2.15)

Then, $(Q(\tau),P(\tau))$ and $(q(t),p(t))$ both satisfy Hamilton’s equations for the original Hamiltonian $H$ with the same initial values, so they must be the same. Note that the Hessian is given by

$$\frac{\partial^{2}\bar{H}}{\partial\bar{p}^{2}}=\begin{bmatrix}\frac{\partial H}{\partial p}\nabla_{p}g(\bar{q},p)^{\top}+g(\bar{q},p)\frac{\partial^{2}H}{\partial p^{2}}+\nabla_{p}g(\bar{q},p)\frac{\partial H}{\partial p}^{\top}\ \ \nabla_{p}g(\bar{q},p)\\ \nabla_{p}g(\bar{q},p)^{\top}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad 0\end{bmatrix},$$

(2.16)

which will be singular in many cases. The degeneracy of the Hamiltonian $\bar{H}$ implies that there is no corresponding Type I Lagrangian formulation. This approach works seamlessly with the existing methods and theorems for Hamiltonian variational integrators, but where the system under consideration is the transformed Hamiltonian system resulting from the Poincaré transformation. We can use a symplectic integrator with constant time-step in fictive time $\tau$ on the Poincaré transformed system, which will have the effect of integrating the original system with the desired variable time-step in physical time $t$ via the relation $\frac{dt}{d\tau}=g(q,t,p)$.

We now depart from the traditional reverse-engineering strategy for the Poincaré transformation and present a new way to think about the Poincaré transformed Hamiltonian by deriving it from a variational principle. This simple derivation gives additional insight into the transformation mechanism and provides natural candidates for time-adaptivity on the Lagrangian side and for more general frameworks.

As before, we work in the extended space $(q,\mathfrak{q},p,\mathfrak{p})$ where $\mathfrak{q}=t$ and $\mathfrak{p}$ is the corresponding conjugate momentum, and consider a time transformation $t\rightarrow\tau$ given by

$$\frac{dt}{d\tau}=g(q,t,p).$$

(3.1)

We define an extended action functional $\mathfrak{S}:C^{2}([0,T],T^{*}\bar{\mathcal{Q}})\rightarrow\mathbb{R}$ by

	$\displaystyle\mathfrak{S}(\bar{q}(\cdot),\bar{p}(\cdot))$	$\displaystyle=\bar{p}(T)\bar{q}(T)-\int_{0}^{T}{\left[\bar{p}(t)\dot{\bar{q}}(t)-H(q(t),t,p(t))-\mathfrak{p}(t)\right]dt}$		(3.2)
		$\displaystyle=\bar{p}(T)\bar{q}(T)-\int_{\tau(t=0)}^{\tau(t=T)}{\left[\bar{p}(\tau)\frac{d\tau}{dt}\bar{q}^{\prime}(\tau)-H(q(\tau),\mathfrak{q}(\tau),p(\tau))-\mathfrak{p}(\tau)\right]\frac{dt}{d\tau}d\tau}$		(3.3)
		$\displaystyle=\bar{p}(T)\bar{q}(T)-\int_{\tau(t=0)}^{\tau(t=T)}{\left\{\bar{p}(\tau)\bar{q}^{\prime}(\tau)-\frac{dt}{d\tau}\left[H(q(\tau),\mathfrak{q}(\tau),p(\tau))+\mathfrak{p}(\tau)\right]\right\}d\tau},$		(3.4)

where we have performed a change of variables in the integral. Then,

$$\mathfrak{S}(\bar{q}(\cdot),\bar{p}(\cdot))=\bar{p}(T)\bar{q}(T)-\int_{\tau(t=0)}^{\tau(t=T)}{\left\{\bar{p}(\tau)\bar{q}^{\prime}(\tau)-g(q(\tau),\mathfrak{q}(\tau),p(\tau))\left[H(q(\tau),\mathfrak{q}(\tau),p(\tau))+\mathfrak{p}(\tau)\right]\right\}d\tau}.$$

(3.5)

Computing the variation of $\mathfrak{S}$ yields

	$\displaystyle\delta\mathfrak{S}$	$\displaystyle=\bar{q}(T)\delta\bar{p}(T)+\bar{p}(T)\delta\bar{q}(T)-\int_{\tau(t=0)}^{\tau(t=T)}{\left[\mathfrak{q}^{\prime}\delta\mathfrak{p}+\mathfrak{p}\delta\mathfrak{q}^{\prime}-\left(g\frac{\partial H}{\partial\mathfrak{q}}+\frac{\partial g}{\partial\mathfrak{q}}(H+\mathfrak{q})\right)\delta\mathfrak{q}-g\delta\mathfrak{p}\right]d\tau}$
		$\displaystyle\qquad\qquad-\int_{\tau(t=0)}^{\tau(t=T)}{\left[q^{\prime}\delta p+p\delta q^{\prime}-\left(g\frac{\partial H}{\partial q}+\frac{\partial g}{\partial q}(H+\mathfrak{p})\right)\delta q-\left(g\frac{\partial H}{\partial p}+\frac{\partial g}{\partial p}(H+\mathfrak{p})\right)\delta p\right]d\tau},$

and using integration by parts and the boundary conditions $\delta\bar{q}(0)=0$ and $\delta\bar{p}(T)=0,$ gives

	$\displaystyle\delta\mathfrak{S}$	$\displaystyle=\int_{\tau(t=0)}^{\tau(t=T)}{\left[p^{\prime}+g\frac{\partial H}{\partial q}+\frac{\partial g}{\partial q}(H+\mathfrak{p})\right]\delta qd\tau}+\int_{\tau(t=0)}^{\tau(t=T)}{\left[g\frac{\partial H}{\partial p}+\frac{\partial g}{\partial p}(H+\mathfrak{p})-q^{\prime}\right]\delta pd\tau}$
		$\displaystyle\qquad\qquad+\int_{\tau(t=0)}^{\tau(t=T)}{\left[\mathfrak{p}^{\prime}+g\frac{\partial H}{\partial\mathfrak{q}}+\frac{\partial g}{\partial\mathfrak{q}}(H+\mathfrak{p})\right]\delta\mathfrak{q}d\tau}-\int_{\tau(t=0)}^{\tau(t=T)}{\left[\mathfrak{q}^{\prime}-g\right]\delta\mathfrak{p}d\tau}.$

Thus, the condition that $\mathfrak{S}(\bar{q}(\cdot),\bar{p}(\cdot))$ is stationary with respect to the boundary conditions $\delta\bar{q}(0)=0$ and $\delta\bar{p}(T)=0$ is equivalent to $(\bar{q}(\cdot),\bar{p}(\cdot))$ satisfying Hamilton’s canonical equations corresponding to the Poincaré transformed Hamiltonian,

	$\displaystyle\mathfrak{q}^{\prime}=g(q,\mathfrak{q},p),$		(3.6)
	$\displaystyle q^{\prime}=g(q,\mathfrak{q},p)\frac{\partial H}{\partial p}(q,\mathfrak{q},p)+\frac{\partial g}{\partial p}(q,\mathfrak{q},p)\left[H(q,\mathfrak{q},p)+\mathfrak{p}\right],$		(3.7)
	$\displaystyle p^{\prime}=-g(q,\mathfrak{q},p)\frac{\partial H}{\partial q}(q,\mathfrak{q},p)-\frac{\partial g}{\partial q}(q,\mathfrak{q},p)\left[H(q,\mathfrak{q},p)+\mathfrak{p}\right],$		(3.8)
	$\displaystyle\mathfrak{p}^{\prime}=-g(q,\mathfrak{q},p)\frac{\partial H}{\partial\mathfrak{q}}(q,\mathfrak{q},p)-\frac{\partial g}{\partial\mathfrak{q}}(q,\mathfrak{q},p)\left[H(q,\mathfrak{q},p)+\mathfrak{p}\right].$		(3.9)

An alternative way to reach the same conclusion is by interpreting equation (3.5) as the usual Type II action functional for the modified Hamiltonian,

$$g(q(\tau),\mathfrak{q}(\tau),p(\tau))\left[H(q(\tau),\mathfrak{q}(\tau),p(\tau))+\mathfrak{p}(\tau)\right],$$

(3.10)

which coincides with the Poincaré transformed Hamiltonian.

We will now derive a mechanism for time-adaptivity on the Lagrangian side by mimicking the derivation of the Poincaré Hamiltonian. We will work in the extended space $\bar{q}=(q,\mathfrak{q},\lambda)^{\top}\in\bar{\mathcal{Q}}$ where $\mathfrak{q}=t$ and $\lambda$ is a Lagrange multiplier used to enforce the time rescaling $\frac{dt}{d\tau}=g(t).$ Consider the action functional $\mathfrak{S}:C^{2}([0,T],T\bar{\mathcal{Q}})\rightarrow\mathbb{R}$ given by

	$\displaystyle\mathfrak{S}(\bar{q}(\cdot),\dot{\bar{q}}(\cdot))$	$\displaystyle=\int_{0}^{T}{\left[L(q(t),\dot{q}(t),\mathfrak{q}(t))-\lambda(t)\left(\frac{d\mathfrak{q}}{d\tau}-g(\mathfrak{q}(t))\right)\right]dt}$		(3.11)
		$\displaystyle=\int_{\tau(t=0)}^{\tau(t=T)}{\left[\frac{dt}{d\tau}L\left(q(\tau),\frac{d\tau}{dt}q^{\prime}(\tau),\mathfrak{q}(\tau)\right)-\lambda(\tau)\frac{dt}{d\tau}\left(\frac{d\mathfrak{q}}{d\tau}-g(\mathfrak{q}(\tau))\right)\right]d\tau}$		(3.12)
		$\displaystyle=\int_{\tau(t=0)}^{\tau(t=T)}{\left[\mathfrak{q}^{\prime}(\tau)L\left(q(\tau),\frac{d\tau}{dt}q^{\prime}(\tau),\mathfrak{q}(\tau)\right)-\lambda(\tau)\mathfrak{q}^{\prime}(\tau)\left[\mathfrak{q}^{\prime}(\tau)-g(\mathfrak{q}(\tau))\right]\right]d\tau},$		(3.13)

where, as before, we have performed a change of variables in the integral. This is the usual Type I action functional for the extended autonomous Lagrangian,

$$\bar{L}(\bar{q}(\tau),\bar{q}^{\prime}(\tau))=\mathfrak{q}^{\prime}(\tau)L\left(q(\tau),\frac{d\tau}{dt}q^{\prime}(\tau),\mathfrak{q}(\tau)\right)-\lambda(\tau)\mathfrak{q}^{\prime}(\tau)\left[\mathfrak{q}^{\prime}(\tau)-g(\mathfrak{q}(\tau))\right].$$

(3.14)

We now introduce a discrete variational formulation of these continuous Lagrangian mechanics. Suppose we are given a partition $0=\tau_{0}<\tau_{1}<\ldots<\tau_{N}=\mathcal{T}$ of the interval $[0,\mathcal{T}]$, and a discrete curve in $\mathcal{Q}\times\mathbb{R}\times\mathbb{R}$ denoted by $\{(q_{k},\mathfrak{q}_{k},\lambda_{k})\}_{k=0}^{N}$ such that $q_{k}\approx q(\tau_{k})$, $\mathfrak{q}_{k}\approx\mathfrak{q}(\tau_{k})$, and $\lambda_{k}\approx\lambda(\tau_{k})$. Consider the discrete action functional,

$$\bar{\mathfrak{S}}_{d}\left(\{(q_{k},\mathfrak{q}_{k},\lambda_{k})\}_{k=0}^{N}\right)=\sum_{k=0}^{N-1}{\left[L_{d}(q_{k},\mathfrak{q}_{k},q_{k+1},\mathfrak{q}_{k+1})-\lambda_{k}\frac{\mathfrak{q}_{k+1}-\mathfrak{q}_{k}}{\tau_{k+1}-\tau_{k}}+\lambda_{k}g(\mathfrak{q}_{k})\right]\frac{\mathfrak{q}_{k+1}-\mathfrak{q}_{k}}{\tau_{k+1}-\tau_{k}}},$$

(3.16)

where $L_{d}(q_{k},\mathfrak{q}_{k},q_{k+1},\mathfrak{q}_{k+1})$ is obtained by approximating the exact discrete Lagrangian, which is related to Jacobi’s solution of the Hamilton–Jacobi equation and is the generating function for the exact time-$h$ flow map. It is given by the extremum of the action integral from $\tau_{k}$ to $\tau_{k+1}$ over twice continuously differentiable curves $(q,\mathfrak{q})\in\mathcal{Q}\times\mathbb{R}$ satisfying the boundary conditions $(q(\tau_{k}),\mathfrak{q}(\tau_{k}))=(q_{k},\mathfrak{q}_{k}),$ and $(q(\tau_{k+1}),\mathfrak{q}(\tau_{k+1}))=(q_{k+1},\mathfrak{q}_{k+1})$:

$$L_{d}(q_{k},\mathfrak{q}_{k},q_{k+1},\mathfrak{q}_{k+1})\approx\operatorname*{ext}_{\begin{subarray}{c}(q,\mathfrak{q})\in C^{2}([\tau_{k},\tau_{k+1}],\mathcal{Q}\times\mathbb{R})\\ (q,\mathfrak{q})(\tau_{k})=(q_{k},\mathfrak{q}_{k}),\text{ }(q,\mathfrak{q})(\tau_{k+1})=(q_{k+1},\mathfrak{q}_{k+1})\end{subarray}}\text{ }\int_{\tau_{k}}^{\tau_{k+1}}{L\left(q,\frac{q^{\prime}}{g(\mathfrak{q})},\mathfrak{q}\right)d\tau}.$$

(3.17)

In practice, we can obtain an approximation by replacing the integral with a quadrature rule, and extremizing over a finite-dimensional function space instead of $C^{2}([\tau_{k},\tau_{k+1}],\mathcal{Q}\times\mathbb{R})$. This discrete functional $\bar{\mathfrak{S}}_{d}$ is a discrete analogue of the action functional $\bar{\mathfrak{S}}:C^{2}([0,T],\mathcal{Q}\times\mathbb{R}\times\mathbb{R})\rightarrow\mathbb{R}$ given by

	$\displaystyle\bar{\mathfrak{S}}(q(\cdot),\mathfrak{q}(\cdot),\lambda(\cdot))$	$\displaystyle=\int_{0}^{\mathcal{T}}{\bar{L}\left(q(\tau),\mathfrak{q}(\tau),\lambda(\tau),q^{\prime}(\tau),\mathfrak{q}^{\prime}(\tau),\lambda^{\prime}(\tau)\right)d\tau}$		(3.18)
		$\displaystyle=\int_{0}^{\mathcal{T}}{\left[L\left(q,\frac{q^{\prime}}{g(\mathfrak{q})},\mathfrak{q}\right)-\lambda\mathfrak{q}^{\prime}+\lambda g(\mathfrak{q})\right]\mathfrak{q}^{\prime}d\tau}.$		(3.19)

We can derive the following result which relates a discrete Type I variational principle to a set of discrete Euler–Lagrange equations: