Contact variational integrators

TL;DR

This paper introduces geometric numerical integrators for contact flows based on a discretized variational principle, demonstrating their properties and comparing their performance to symplectic integrators using a damped harmonic oscillator example.

Contribution

It develops contact variational integrators, proving they are contact transformations and establishing their backward error analysis, extending variational methods to contact geometry.

Findings

The integrators are contact transformations.

They can be derived from a variational principle.

Performance comparison shows advantages over symplectic integrators.

Abstract

We present geometric numerical integrators for contact flows that stem from a discretization of Herglotz' variational principle. First we show that the resulting discrete map is a contact transformation and that any contact map can be derived from a variational principle. Then we discuss the backward error analysis of our variational integrators, including the construction of a modified Lagrangian. Throughout the paper we use the damped harmonic oscillator as a benchmark example to compare our integrators to their symplectic analogues.

Equations

$$\eta =\mathrm{d}z-\sum _i{p}_i\mathrm{d}{x}_i.$$

$${\mathcal{L}}_{X_H}\eta =f_H\eta \text{and}\eta (X_H)=-H,$$

$$\{\begin{array}{l}\dot{x}=\frac{\partial H}{\partial p}\\ \dot{p}=-\frac{\partial H}{\partial x}-p\frac{\partial H}{\partial z}\\ \dot{z}=p\frac{\partial H}{\partial p}-H,\end{array}$$

$$\frac{\mathrm{d}H}{\mathrm{d}t}=-H\frac{\partial H}{\partial z}.$$

$$\{\begin{array}{l}\dot{x}=p\\ \dot{p}=-V^{\prime }(x)-\alpha p\\ \dot{z}=p^2-H.\end{array}$$

$$S={\int }_0^T\mathcal{L}(t,x(t),\dot{x}(t))\mathrm{d}t$$

$$\dot{z}(t)=\mathcal{L}(t,x(t),\dot{x}(t),z(t)).$$

$$\frac{\partial \mathcal{L}}{\partial x}-\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \mathcal{L}}{\partial \dot{x}}+\frac{\partial \mathcal{L}}{\partial z}\frac{\partial \mathcal{L}}{\partial \dot{x}}=0,$$

$$\delta \dot{z}=\frac{\partial \mathcal{L}}{\partial x}\delta x+\frac{\partial \mathcal{L}}{\partial \dot{x}}\delta \dot{x}+\frac{\partial \mathcal{L}}{\partial z}\delta z.$$

$$A(t)=\frac{\partial \mathcal{L}}{\partial x}\delta x+\frac{\partial \mathcal{L}}{\partial \dot{x}}\delta \dot{x}\text{and}B(t)={\int }_0^t\frac{\partial \mathcal{L}}{\partial z}(\tau )\mathrm{d}\tau ,$$

$$\delta \dot{z}(t)=A(t)+\frac{\mathrm{d}B(t)}{\mathrm{d}t}\delta z$$

$$\delta z(t)=e^{B(t)}\left[{\int }_0^{t}A(\tau )e^{-B(\tau )}\mathrm{d}\tau +\delta z(0)\right].$$

$$\delta z(T)$$

$$\displaystyle\begin{split}&=e^{B(T)}\bigg{[}\int_{0}^{T}\left(\frac{\partial{\mathcal{L}}}{\partial{x}}-\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}+\frac{\partial{\mathcal{L}}}{\partial{z}}\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}\right)\delta x\,e^{-B(\tau)}\,\mathrm{d}\tau\\ &\hskip 56.9055pt+\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}(T)\delta x(T)\,e^{-B(T)}-\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}(0)\delta x(0)+\delta z(0)\bigg{]}.\end{split}$$

$$\ddot{x}$$

$$\dot{z}$$

$$\dot{E}=\frac{\partial \mathcal{L}}{\partial z}E$$

$$\frac{\partial \mathcal{L}}{\partial \dot{x}}(t)\dot{x}(t)-\dot{z}(t)=e^{B(t)}\left(\frac{\partial \mathcal{L}}{\partial \dot{x}}(0)\dot{x}(0)-\dot{z}(0)\right),$$

$$\mathrm{d}z-p\mathrm{d}x,$$

$$v=(\delta x(0),\delta \dot{x}(0),\delta z(0))\in T_{(x(0),\dot{x}(0),z(0))}M$$

$$F_{\ast }^{t}v=(\delta x(t),\delta \dot{x}(t),\delta z(t))\in T_{(x(t),\dot{x}(t),z(t))}M$$

$$\mathrm{d}z(F_{\ast }^{t}v)=e^{B(t)}\left[p(t)e^{-B(t)}\mathrm{d}x(F_{\ast }^{t}v)-p(0)\mathrm{d}x(v)+\mathrm{d}z(v)\right],$$

$${\left(F^t\right)}^{\ast }(\mathrm{d}z-p\mathrm{d}x)=e^{B(t)}(\mathrm{d}z-p\mathrm{d}x),$$

$$\exp (B(t))=\exp \left({\int }_0^t\frac{\partial \mathcal{L}}{\partial z}(\tau )\mathrm{d}\tau \right).$$

$$H(x,p,z)=p\dot{x}-\mathcal{L}(x,\dot{x},z),$$

$$\frac{\partial H}{\partial z}=-\frac{\partial \mathcal{L}}{\partial z},$$

$$\dot{H}=-\frac{\partial H}{\partial z}H.$$

$$\frac{\partial H}{\partial p}$$

$$\frac{\partial H}{\partial x}$$

$$=-\dot{p}-p\frac{\partial H}{\partial z}.$$

$$z_{j+1}-z_j=hL(x_j,x_{j+1},z_j,z_{j+1};h).$$

$$\frac{\partial z_{j+1}}{\partial x_j}=0,$$

Full text

Keywords: contact geometry, geometric integrators, Herglotz’ variational principle

MSC2010: 65D30, 34K28, 34A26

1 Introduction

The last few years have seen a rise in importance of the field of contact geometry. As the theory gets more relevant to scientific applications, there is an increasing demand for the development of numerical integrators preserving the contact structure.

Contact geometry appears in fluid dynamics [19, 38, 47], thermodynamics [7, 28, 29, 43, 53], statistical physics [9], statistics [5, 27], quantum mechanics [11, 14, 21, 33, 46], gravity [40], information geometry [3, 4, 26], shape dynamics [52], biology [8, 30], optimal control [37, 45], and integrable systems [6, 35, 36, 49, 50]. One of the applications that has recently attracted a lot of attention is the classical mechanics of dissipative systems [1, 2, 10, 17, 18, 23, 41], on which we will focus in this paper.

Contact geometry can be thought of as an odd-dimensional analogue of symplectic geometry. A contact manifold is a pair $(M,\xi)$ where $M$ is an $(2n+1)$-dimensional manifold and $\xi\subset{TM}$ is a contact structure, that is, a maximally non-integrable distribution of hyperplanes. Locally, such distribution is given by the kernel of a one form $\eta$ satisfying $\eta\wedge(\mathrm{d}\eta)^{n}\neq 0$ (see e.g. [24] for more details). The 1-form $\eta$ is called the contact form. Note that if we multiply $\eta$ by a non-vanishing function we obtain another 1-form giving rise to the same contact structure $\xi$. This means that in order to preserve $\xi$ we can act on the 1-form by conformal transformations. Thus we must keep in mind that $\eta$ is just a representative element in an equivalence class of 1-forms describing the same $\xi$.

Once we fix $\eta$, Darboux’s theorem for contact manifolds states that for any point on $M$ there exists a neighbourhood with local coordinates $(x_{1},\ldots,x_{n},p_{1},\ldots,p_{n},z)$ such that the contact 1-form can be written as

[TABLE]

Throughout the paper we will write $\mathrm{d}z-p\,\mathrm{d}x$ as a short form for $\eta$.

Moreover, given $\eta$, to any smooth function $H:M\rightarrow\mathbb{R}$ we can associate a contact Hamiltonian vector field $X_{H}$, defined by

[TABLE]

where $\mathcal{L}$ is the Lie derivative, $f_{H}=-R_{\eta}(H)$ and $R_{\eta}$ is the Reeb vector field corresponding to $\eta$ [24]. In Darboux coordinates the flow of $X_{H}$ is given by

[TABLE]

where $x=(x_{1},\ldots,x_{n})$, $p=(p_{1},\ldots,p_{n})$ and the standard scalar product is assumed where two vectors are multiplied. The flow of a contact Hamiltonian system preserves the contact structure, but it does not preserve the Hamiltonian. Instead we have

[TABLE]

For example a Hamiltonian of the form $H=\frac{1}{2}p^{2}+V(x)+\alpha z$ leads to the equations of motion of a damped mechanical system:

[TABLE]

Another remarkable similarity with standard symplectic Hamiltonian systems is the fact that contact Hamiltonian systems have an associated variational principle, which is due to Herglotz [25, 34] (see also [13, 55]), and a corresponding theory of generating functions [10].

Geometric integrators for contact Hamiltonian systems have been studied in [20] by exploiting their symplectification and the corresponding generating functions. However, a variational approach is missing. So far, [20] has received little attention, most likely because the authors did not discuss any physically relevant examples.

In this paper we present a natural way to develop numerical integrators for contact systems by exploiting Herglotz’ variational principle. Our result furnishes a variational scheme to integrate contact Hamiltonian systems in such a way that the contact structure is preserved. Furthermore, in analogy with the theory of symplectic numerical integrators, we find that modified equations for our method are again contact systems. This suggests that the numerical results will remain very close to a “nearby” contact system for very long times.

The purpose of this paper is to lay the theoretical groundwork of contact integrators and to show their promise with some simple numerical experiments. To keep the discussion direct and self-contained, the main results will be presented only for contact systems without an explicit time dependence, with a section showing how the method is easily extended to general contact systems by means of an explanatory example. Furthermore, some interesting discussions and developments are postponed to future works. These include the extension to a more general sub-Riemannian setting and the comparison with [20] and other known approaches to some physically relevant examples.

It is important to remark at this point that not every kind of dissipative system can be written with a contact Hamiltonian of the form above, as the geometry underlying the construction will enforce some structure. For more details on this, we refer the readers to the thorough investigation in [14, 18, 23]. Nevertheless, the contact Hamiltonian structure allows to describe a large number of physically relevant systems, including most of the ones from the literature cited at the beginning of this introduction.

The paper is structured as follows. In section 2 we set the stage introducing Herglotz’ variational principle and some of its relevant properties. In section 3 we develop the central idea of the paper defining a contact integrator obtained from a discretization of Herglotz’ variational principle. In section 4 we study the modified equations for the contact integrators introduced in section 3, showing that up to truncation errors, the numerical solutions are interpolated by contact systems. In section 5 we present an example to show how the ideas can be extended in a straightforward fashion to systems with an explicit time dependence. Finally, in section 6 we illustrate by numerical examples how our integrators perform on contact systems in comparison to symplectic integrators.

2 Herglotz’ variational principle

The usual variational principle in mechanics looks for a curve $x:[0,T]\rightarrow Q$ in configuration space $Q$, such that the action integral

[TABLE]

is critical with respect to variations of $x$ that vanish at the endpoints, where $\mathcal{L}:\mathbb{R}\times TQ\rightarrow\mathbb{R}$ is a given Lagrange function. Herglotz [34] generalized this variational principle by defining the action in terms of a differential equation instead of an integral.

Definition 1**.**

Let $\mathcal{L}:\mathbb{R}\times TQ\times\mathbb{R}\rightarrow\mathbb{R}$. Given a curve $x:[0,T]\rightarrow Q$, define the function $z:[0,T]\rightarrow\mathbb{R}$ by an initial condition $z(0)=z_{0}$ and the differential equation

[TABLE]

The curve $x$ is a solution to Herglotz’ variational principle with initial condition $z_{0}$ if every variation of $x$ that vanishes at the boundary of $[0,T]$ leaves the action $z(T)$ invariant.

If $\mathcal{L}$ does not depend on $z$, then the differential equation (2) is solved by the integral (1) and Herglotz’ variational principle reduces to the classical variational principle. A modern discussion of Herglotz’ variational principle can be found for example in [25].

Proposition 1**.**

A (sufficiently regular) curve $x$ is a solution to Herglotz’ variational principle if and only if it satisfies the generalized Euler-Lagrange equations

[TABLE]

where $z$ is given in terms of $x$ by Equation (2).

Note that the Euler-Lagrange equations are not linear in $\mathcal{L}$, hence they are not invariant under scaling of the Lagrangian $\mathcal{L}$.

Proof of Proposition 1.

Consider an arbitrary variation $\delta x$ of $x$, vanishing at the endpoints, and the corresponding induced variation $\delta z$ of $z$. Since the initial condition $z(0)=z_{0}$ is independent of $x$ we have $\delta z(0)=0$. From Equation (2) it follows that

[TABLE]

If we set

[TABLE]

this differential equation reads

[TABLE]

and its solution is

[TABLE]

Plugging in the expression for $A$ and noting that $\frac{\mathrm{d}B}{\mathrm{d}t}=\frac{\partial{\mathcal{L}}}{\partial{z}}$ we find

[TABLE]

The boundary terms vanish because $\delta x(0)=\delta x(T)=\delta z(0)=0$. Since $\delta x$ is otherwise arbitrary, the action $z(T)$ is critical if and only if Equation (3) holds. ∎

If the classical variational principle is satisfied on the interval $[0,T]$ it is automatically satisfied on any subinterval. For the Herglotz variational principle this property is not obvious from the definition, but it still follows from the generalized Euler-Lagrange equations.

Proposition 2**.**

If $x:[0,T]\rightarrow Q$ solves the Herglotz variational principle with initial condition $z_{0}$, then for any interval $[a,b]\subset[0,T]$, the restriction $x|_{[a,b]}$ solves the Herglotz variational principle with initial condition $z(a)$.

Proof.

If $x$ is critical on $[0,T]$ then the generalized Euler-Lagrange equations are satisfied everywhere on this interval. In particular, they hold on $[a,b]$, hence $x$ is critical on $[a,b]$. ∎

In the following we will assume that the Lagrangian is regular, i.e. $\left|\frac{\partial{{}^{2}\mathcal{L}}}{\partial{\dot{x}^{2}}}\right|\neq 0$. Then the generalized Euler-Lagrange equations can be written explicitly as a second order ODE. Together with the evolution of $z$ we find the system of ODEs

[TABLE]

An important aspect of the Herglotz variational principle is that the energy is not conserved (unless the Lagrangian is independent of $z$). Instead we find a differential equation governing its evolution.

Proposition 3**.**

If the Lagrangian does not explicitly depend on time, then the energy $E=\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}\dot{x}-\mathcal{L}$ satisfies the differential equation

[TABLE]

Proof.

Consider a uniform time-shift of the critical curve $x$ and the function $z$. Then $\delta x=\dot{x}$ and $\delta z=\dot{z}$. If the Lagrangian does not explicitly depend on time, it follows from Equation (4) that

[TABLE]

for any $t\in[0,T]$ because criticality on $[0,T]$ implies criticality on the subinterval $[0,t]$. It follows that $E=\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}\dot{x}-\dot{z}$ satisfies Equation (5). ∎

The usual argument that Lagrangian flows are symplectic, as presented for example in [42, Section 1.2], can be extended to show that flows of Herglotz’ variational principle are contact transformations.

Proposition 4**.**

Let $M=TQ\times\mathbb{R}$ with local coordinates $(x,\dot{x},z)$. The flow $F:\mathbb{R}\times M\rightarrow M:(t,x,\dot{x},z)\mapsto F^{t}(x,\dot{x},z)$ of the generalized Euler-Lagrange equations consists of contact transformations $F^{t}$ with respect to the 1-form

[TABLE]

where $p=\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}$.

Proof.

On solutions of the generalized Euler-Lagrange equations, the value of $z(t)$ is uniquely defined by the initial values $x(0)$, $\dot{x}(0)$, and $z(0)$. Any variation

[TABLE]

of the initial data induces a variation

[TABLE]

at the endpoint, where $F^{t}_{*}:TM\rightarrow TM$ denotes the pushforward of $F^{t}$.

Since we are working on solutions of the generalized Euler-Lagrange equations, the integrand in Equation (4) vanishes and only the boundary terms remain. They can be written as

[TABLE]

where $p=\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}$. It follows that

[TABLE]

where $\left(F^{t}\right)^{*}:T^{*}M\rightarrow T^{*}M$ denotes the pullback of $F^{t}$. Hence the flow consists of contact transformations with respect to the 1-form $\mathrm{d}z-p\,\mathrm{d}x$ with conformal factor

[TABLE]

∎

To close this section, let us briefly state the link of the Herglotz variational principle to the more common Hamiltonian formulation of contact dynamics. The contact Hamiltonian $H:T^{*}Q\times\mathbb{R}\rightarrow\mathbb{R}$ is defined by Legendre transformation

[TABLE]

where $\dot{x}$ is eliminated from the right hand side by $p=\frac{\partial{\mathcal{L}}}{\partial{\dot{x}}}$. Taking the partial derivative with respect to $z$ gives

[TABLE]

hence from Equation (5) it follows that

[TABLE]

Differentiating instead with respect to $p$ and $x$ gives us the contact Hamiltonian equations:

[TABLE]

3 Discrete Herglotz variational principle

As is standard in discrete mechanics, in what follows we replace $TQ$ by $Q^{2}$ (see e.g. [42]).

Definition 2**.**

Let $L:Q^{2}\times\mathbb{R}^{2}\times\mathbb{R}\rightarrow\mathbb{R}$ and $h>0$. Given a discrete curve $x=(x_{0},\ldots,x_{N})\in Q^{N+1}$, we define $z=(z_{0},\ldots,z_{N})\in\mathbb{R}^{N+1}$ by $z_{0}=0$ and

[TABLE]

The discrete curve $x$ is a solution to the discrete Herglotz variational principle if

[TABLE]

for all $j\in\{1,\ldots,N-1\}$.

Note that Equation (7) is a discrete version of Equation (2), and that Equation (8) means that for a critical discrete curve $x$, a variation of $x_{k}$ can affect $z_{k}$ but none of the other $z_{j}$. In particular, this implies that $z_{N}$ is critical with respect to variations of $x_{1},\ldots,x_{N-1}$. Most of the time we will consider a fixed step size $h$ and omit it from the notation of the discrete Lagrangian $L(x_{j},x_{j+1},z_{j},z_{j+1})$.

Theorem 1**.**

For a sufficiently small step size $h$, the discrete curve $x$ is a solution of the discrete Herglotz variational principle, with $z$ defined by Equation (7), if and only if it satisfies the discrete generalized Euler-Lagrange equations

[TABLE]

where $\operatorname{D}_{i}$ denotes the partial derivative with respect to the $i$-th entry.

Note that while in general the $x_{j}$ have several components, the $z_{j}$ are always scalar, hence $\operatorname{D}_{1}L$ and $\operatorname{D}_{2}L$ are vectors but $\operatorname{D}_{3}L$ and $\operatorname{D}_{4}L$ are scalars.

Equation (9) is equivalent to

[TABLE]

In the first line one recognizes the usual discrete Euler-Lagrange equations. The term in the second line is a discretization of $\frac{\partial{L}}{\partial{\dot{x}}}\frac{\partial{L}}{\partial{z}}$.

Proof of Theorem 1.

From Equation (7) it follows that

[TABLE]

On solutions we have

[TABLE]

where the derivative $\operatorname{D}_{3}L$ is omitted because $\frac{\partial{z_{j-1}}}{\partial{x_{j}}}=0$. It follows that

[TABLE]

hence for sufficiently small $h$, $\frac{\partial{z_{j+1}}}{\partial{x_{j}}}=0$ is equivalent to Equation (9). ∎

In analogy to the continuous case, we will always assume the non-degeneracy condition

[TABLE]

which guarantees that for sufficiently small $h$, Equation (9) can be solved for $x_{j+1}$.

Theorem 2**.**

The map $Q^{2}\times\mathbb{R}\mapsto Q^{2}\times\mathbb{R}:(x_{j-1},x_{j},z_{j-1})\mapsto(x_{j},x_{j+1},z_{j}),$ given by the generalized discrete Euler-Lagrange equations, induces a map

[TABLE]

where

[TABLE]

and

[TABLE]

The map $F$ is a contact transformation with respect to the 1-form $\mathrm{d}z-p\,\mathrm{d}x$.

Equations (12) and (13) define the discrete Legendre transforms for contact systems, compare [42, Section 1.5].

Proof.

First note that the second equality in Equation (11) follows from Equation (9) and the definitions (12) and (13).

To prove that $F$ is a contact transformation, we consider the case $j=2$. The general statement is obtained by shifting all indices by the same integer.

From

[TABLE]

it follows that

[TABLE]

hence, on solutions of the generalized Euler-Lagrange equations,

[TABLE]

It follows that

[TABLE]

Note that the conformal factor

[TABLE]

is consistent with the continuous $\exp\!\left(\int_{0}^{h}\frac{\partial{\mathcal{L}}}{\partial{z}}\mathrm{d}t\right)$, cf. Equation (6). We stress that Theorems 1 and 2 also apply to the case where $L$ does not depend on $z_{2}$, i.e. $\operatorname{D}_{4}L=0$.

A natural question to ask at this point is whether every contact transformation comes from a variational principle. Just like the symplectic counterpart to this question, it is answered in the affirmative using generating functions.

Remarkably, the following result is stronger than the literal inverse to Theorem 2, which said that any discrete Lagrangian $L(x_{1},x_{2},z_{1},z_{2})$ yields a contact transformation. We will show that every contact transformation can be obtained from a discrete Lagrangian $L(x_{1},x_{2},z_{1})$ that does not depend on the second instance of $z$.

We stress the importance of this result, since it implies that every contact integrator is variational.

Theorem 3**.**

Iterations of any contact transformation $(x_{0},p_{0},z_{0})\mapsto(x_{1},p_{1},z_{1})$ yield a discrete curve $x=(x_{0},\ldots,x_{N})$ that solves the discrete Herglotz variational principle for some discrete Lagrangian $L(x_{j},x_{j+1},z_{j})$.

Note that $L$ does not depend on $z_{j+1}$ in the statement of Theorem 3. Hence without loss of generality we can restrict our attention to Lagrangians depending only on the first of the $z$-coordinates, as we will do e.g. in Example 1.

Proof of Theorem 3.

As pointed out in [10], the coordinate $z_{1}$ of a contact transformation $(x_{0},p_{0},z_{0})\mapsto(x_{1},p_{1},z_{1})$ can be considered as a generating function. We have

[TABLE]

Writing $z_{1}=S(x_{0},x_{1},p_{0},p_{1},z_{0})$ we find

[TABLE]

It follows that $\frac{\partial{S}}{\partial{p_{0}}}=\frac{\partial{S}}{\partial{p_{1}}}=0$ and

[TABLE]

Note that $S$ does not depend on $p_{0}$ or $p_{1}$, hence from now on we will write $S(x_{0},x_{1},z_{0})$. Setting

[TABLE]

iterations of the contact map satisfy

[TABLE]

Furthermore, using Equation (14) we calculate that

[TABLE]

so the discrete curve $x$ obtained by iteration of the contact map satisfies the discrete Herglotz variational principle for $L$. ∎

Example 1**.**

The Lagrangian $\mathcal{L}=\frac{1}{2}\dot{x}^{2}-V(x)-\alpha z$ describes a mechanical system with Rayleigh dissipation (i.e. a friction force linear in the velocity). The generalized Euler-Lagrange equation is

[TABLE]

Note that $x$ need not be a scalar: the Lagrangian $\mathcal{L}=\frac{1}{2}|\dot{x}|^{2}-V(x)-\alpha z$ yields the analogous multi-component equation. This contrasts many other variational descriptions of the damped harmonic oscillator, which only apply to the scalar case [44, 15]. The same comment applies to the following discretization, which we write down for scalar $x$ but can easily be adapted to higher dimensions.

A discretization of the Lagrangian is

[TABLE]

Note that this Lagrangian depends only on $z_{j}$, not on $z_{j+1}$. Its discrete generalized Euler-Lagrange equations read

[TABLE]

The discrete momentum can be calculated as

[TABLE]

or

[TABLE]

We can implement the integrator explicitly in position-momentum formulation as

[TABLE]

Let us consider the damped harmonic oscillator, $V(x)=\frac{1}{2}x^{2}$. Its equation of motion is

[TABLE]

The above discrete Lagrangian then becomes

[TABLE]

and it discrete generalized Euler-Lagrange equations read

[TABLE]

The position-momentum formulation of the integrator gives

[TABLE]

Example 2**.**

For the theory of discrete contact systems by itself, it is sufficient to have the Lagrangian depend on $z_{j}$ but not on $z_{j+1}$, as we saw in Theorem 3. For the sake of a good numerical approximation, however, it is beneficial to relax this condition. Continuing the example of a damped mechanical system, we can take the discrete Lagrangian

[TABLE]

Note the difference with Example 1: now $L$ depends also on $z_{j+1}$. Its discrete generalized Euler-Lagrange equations read

[TABLE]

Equations (16) and (20) are related by a simple change in the parameter $\alpha$. This minor difference should not be dismissed, though, as the discrete Lagrangian (19) is a second order approximation of the continuous Lagrangian, compared to the first order approximation of Equation (15). What we mean by this will be clarified in the next section: see Example 4 and Example 5.

The discrete momentum for the Lagrangian (19) can be calculated as

[TABLE]

or

[TABLE]

The generalized Euler-Lagrange equations state that both formulas for the discrete momentum agree. On solutions, we have that

[TABLE]

We can implement the integrator explicitly in position-momentum formulation as

[TABLE]

Equation (20) is also the second order difference equation corresponding to the leapfrog (Störmer-Verlet) method

[TABLE]

Indeed, eliminating the momentum $\pi$ from this system we find

[TABLE]

hence

[TABLE]

which is equivalent to Equation (20).

The momenta at integer steps can be included in the leapfrog method by adding one internal stage:

[TABLE]

We find that

[TABLE]

hence when initialized with the same momentum, the difference between the result of our contact method and the leapfrog solution will be of order $h^{2}(\alpha+\alpha^{2})$. Whether $p_{0}$ or $\pi_{0}$ is a better approximation for the true initial momentum $\dot{x}$ depends on the initial conditions.

Example 3**.**

For a more general contact system, motivated by [52], consider the Lagrangian

[TABLE]

The equations of motion are

[TABLE]

or, in position-momentum formulation,

[TABLE]

Note that the Euler-Lagrange equation explicitly involves $z$ in this case, so it is not possible to solve the equations for $x$ and $p$ separately. This means that symplectic integrators cannot be applied (unless one adds a dummy variable to the systems in order to obtain an even-dimensional system once again). In comparison, in Examples 1 and 2 a symplectic integrator could be applied, but it would not respect the contact structure.

Consider the discrete Lagrangian

[TABLE]

The discrete momenta are

[TABLE]

and

[TABLE]

Hence we find an implicit contact integrator

[TABLE]

4 Backward error analysis

A central idea to explain the long-time behavior of symplectic integrators is the study of modified differential equations whose solutions interpolate the discrete solutions of a discrete system of equations. This idea, looking for a perturbed continuous system that exactly corresponds to the discretization, is an example of backward error analysis. It is a well-known and essential fact that if a symplectic integrator is applied to a Hamiltonian equation, then the resulting modified equation is Hamiltonian as well. Similarly, when a classical variational integrator is applied to a Lagrangian system, the resulting modified equation is Lagrangian [54]. Below we establish that an analogous result holds for contact variational integrators.

First let us have a look at the general form of the discrete generalized Euler-Lagrange equations.

Proposition 5**.**

Consider a continuous non-degenerate Lagrangian $\mathcal{L}(x,\dot{x},z)$ with generalized Euler-Lagrange equation $\ddot{x}=f(x,\dot{x},z)$ and a consistent discretization $L(x_{j},x_{j+1},z_{j},z_{j+1};h)$ of $\mathcal{L}$, by which we mean that for any smooth $x$ and $z$ there holds

[TABLE]

Then the discrete generalized Euler-Lagrange equation is a consistent discretization of the continuous generalized Euler-Lagrange equation, i.e. it takes the form

[TABLE]

where for any smooth $x$ and $z$

[TABLE]

We give a formal proof. A rigorous version of this argument is obtained by generalizing the corresponding proof in [54] to the case of the Herglotz variational principle. A different proof strategy can be found in [42, Section 2.3].

Proof of Proposition 5.

Let $x$ be a smooth curve interpolating solutions of the discrete generalized Euler-Lagrange equation (10). Then $\mathcal{E}_{1}+\mathcal{E}_{2}\mathcal{E}_{3}=0$, where

[TABLE]

[TABLE]

and

[TABLE]

We start by showing that $\mathcal{E}_{1}$ is a consistent discretization of the classical Euler-Lagrange equation. In terms of the Taylor expansions of the shifted variables, it becomes

[TABLE]

Since the Lagrangian is assumed to be a consistent discretization, we have

[TABLE]

hence

[TABLE]

In other words

[TABLE]

where each $\mathcal{L}$ is evaluated at $(x(t),\dot{x}(t),z(t))$. We can simplify this expression using the fact that

[TABLE]

and obtain

[TABLE]

A perfectly analogous computation yields that

[TABLE]

Finally, we compute

[TABLE]

It follows that

[TABLE]

The claimed result follows by isolating the term $\ddot{x}$ in the right hand side of this equation and the term $\frac{x(t-h)-2x(t)+x(t+h)}{h^{2}}$ in the left hand side. ∎

Now we turn our attention to the modified equations of the system

[TABLE]

That is, we look for differential equations whose solutions interpolate solutions of the difference equations. The precise definition of a modified equation is a bit more involved because of convergence issues.

Definition 3**.**

The system of modified equations for the difference system (22) is defined by the formal expressions

[TABLE]

such that for any $k\in\mathbb{N}$, every solution $(x,z)$ of the truncated differential equations

[TABLE]

satisfies the difference equations with a defect of order $k+1$, in the sense that

[TABLE]

Given a difference equation of the form (21) one can recursively compute the coefficients $\mathcal{L}_{i}$ and $f_{i}$ of the system of modified equations. Examples of such calculations can be found for example in [32, Chapter IX] and [54]. Note that in the leading order of the modified equations we recover the original differential equations. This is because we are dealing with consistent discretizations. The additional terms of the power series contain information about the integrator. In particular, the order of an integrator is the smallest $k>0$ such that the $h^{k}$-term in Equation (23) is non-zero.

If we look at the difference equation for $z$ by itself, i.e. with an arbitrary smooth curve $x$ instead of one interpolating a discrete solution, we get a different modified equation, depending on higher derivatives of $x$:

[TABLE]

such that for any $k\in\mathbb{N}$ and any smooth curve $x$, every solution $z$ of the truncated differential equation

[TABLE]

satisfies

[TABLE]

In the system of modified equations (23), $\mathcal{L}_{\mathrm{mod}}(x,\dot{x},z,h)$ is the Herglotz Lagrangian for $\ddot{x}=f_{\mathrm{mod}}(x,\dot{x},z;h)$ in the following sense:

Theorem 4**.**

A truncation after order $h^{k}$ of the power series $\mathcal{L}_{\mathrm{mod}}(x,\dot{x},z,h)$ yields as its generalized Euler-Lagrange equations $\ddot{x}=f_{\mathrm{mod}}(x,\dot{x},z;h)+\mathcal{O}(h^{k+1})$.

Sketch of proof.

Let $(x(t),z(t))$ be any solution to the truncated system of modified equations (24). By definition of a modified equation, the discrete curve $(x_{j},z_{j})_{j\in\mathbb{Z}}$ defined by $x_{j}=x(jh)$ and $z_{j}=z(jh)$ satisfies the discrete system (22) with a defect of order $\mathcal{O}(h^{k+1})$.

Now consider the action $z_{N}=\bar{z}(Nh)$ where $\bar{z}$ is a solution of the higher order Equation (26), with $x$ as before. The discrete Herglotz variational principle implies that $z_{N}$ is critical with respect to variations of $x(t)$, supported on the interval $(0,Nh)$, again up to a defect of order $\mathcal{O}(h^{k+1})$. This means that $x(t)$ solves the continuous Herglotz problem for $\mathcal{L}_{\mathrm{mod}}^{z}$ with the same defect.

We want to prove the same property with $\mathcal{L}_{\mathrm{mod}}$ instead of $\mathcal{L}_{\mathrm{mod}}^{z}$. If we define $z$ by Equation (24), independent of second or higher derivatives of $x$, then it will only interpolate a discrete solution if $x$ solves the modified equation $\ddot{x}=f_{\mathrm{mod}}(x,\dot{x},z,h)+\mathcal{O}(h^{k+1})$. Hence we do not have the freedom to take variations of $x$ as needed in the argument above.

To show that $\mathcal{L}_{\mathrm{mod}}(x,\dot{x},z,h)$ is nevertheless a Lagrangian for the modified equation, we need to show that replacing higher derivatives of $x$ in $\mathcal{L}_{\mathrm{mod}}^{z}(x,\dot{x},\ddot{x},\ldots,z,h)$ using the modified equation does not change its generalized Euler-Lagrange equations. Sufficient conditions for this are that

[TABLE]

where $x^{(\ell)}$ denotes the $\ell$-th derivative of $x$. These conditions can be obtained from the so-called meshed variational problem, as explained in [54]. Then it follows that $x(t)$ satisfies the Herglotz variational principle for the first order Lagrangian $\mathcal{L}_{\mathrm{mod}}(x,\dot{x},z;h)$. ∎

Example 4** (Example 1 continued).**

Let us calculate the first order approximation of the modified equations for our first discretization of the damped harmonic oscillator. Assume that $x$ is a solution to the modified equation. Then it satisfies Equation (18) when we replace $x_{j}$ by $x(t)$ and $x_{j\pm 1}$ by $x(t\pm h)$:

[TABLE]

A Taylor expansion gives

[TABLE]

where all instances of $x$ and its derivatives are evaluated at $t$. Since in the leading order of the modified equation we recover the original equation, we know that $\ddot{x}=-x-\alpha\dot{x}+\mathcal{O}(h)$, which we can use to simplify the right hand side. We obtain

[TABLE]

Using the same procedure we calculate the modified equation for $\frac{z_{j+1}-z_{j}}{h}=L$, with $L$ given by (17). In terms of the interpolating functions, the difference equation reads

[TABLE]

and its Taylor expansion is

[TABLE]

Solving this for $\dot{z}$ we find

[TABLE]

In the right hand side we replace $\ddot{z}$ using the leading order equation

[TABLE]

to find

[TABLE]

The right hand side of this modified equation should give us the modified Lagrangian. A simple calculation shows that the generalized Euler-Lagrange equation for

[TABLE]

is indeed Equation (27). Note that up to the $\mathcal{O}(h^{2})$ error term, the modified Lagrangian $\mathcal{L}_{\mathrm{mod}}$ is a rescaling of the original Lagrangian $\mathcal{L}$. Unlike for the classical variational principle, this does not imply that both Lagrangians have the same generalized Euler-Lagrange equations.

The second order term of the modified equations can be calculated by including one more term in the Taylor expansions and simplifying the right hand sides using the first order modified equations (27)–(28) instead of the leading order equations. We find

[TABLE]

and

[TABLE]

This process can be continued to recursively find the modified equations to any order.

Example 5** (Example 2 continued).**

We can repeat the above calculation to obtain a modified equation for the second order integrator too. We find

[TABLE]

which shows that the discretization of Example 2 is indeed a second order method. This is a consequence of the symmetry of the discretization:

[TABLE]

5 On the explicit time dependence

Even though in the previous sections we focused on Lagrangians that do not explicitly depend on time, going through the previous proofs and examples one can observe that there is no obstruction to considering explicitly time-dependent systems. In fact, when the system depends explicitly on time, the resulting flow yields a time-dependent contact transformation, in compete analogy to what happens with time-dependent canonical transformations in the symplectic case. We refer to [10] for more on the theory of time-dependent contact transformations and their generating functions. Therefore, with some efforts and modulo a slight complication of the notation in few instances, it is possible to extend all the previous results to such systems in a straightforward way, so that the resulting maps are discretizations of the corresponding time-dependent contact transformations.

How the explicitly time-dependent terms appear in the discrete Lagrangian and the resulting difference equation depend both on the choice of discretization and on the form of time-dependent terms in the continuous Lagrangian. In many cases, such as for external forcing, the time-dependence can be separated neatly and the final result will be elegant and readable. To illustrate the time-dependent case, we build upon the Lagrangian presented in Example 1 and Example 2, and consider a forced damped harmonic oscillator.

Example 6**.**

In general, the Lagrangian of a mechanical system with Raileigh dissipation and external forcing $f(t)$ looks like

[TABLE]

Indeed, the generalized Euler-Lagrange equation in this case is

[TABLE]

Here and in what follows, we emphasize the difference with the equations from Example 2 in boxes.

Let $t_{j}=t_{j-1}+h=t_{0}+jh$, then a natural discretization of the Lagrangian above is

[TABLE]

and the discrete generalized Euler-Lagrange equation reads

[TABLE]

The position-momentum formulation of the integrator is

[TABLE]

6 Numerical results

In this section we discuss the behaviour of our contact variational integrators in comparison with some classical fixed step methods. In what follows we consider the damped harmonic oscillator with and without forcing, integrated using:

• •

our contact variational integrators of both first and second order as presented in the previous examples Example 1 and Example 2 (respectively denoted “Contact (1st)” and “Contact (2nd)”),

• •

the symplectic second order Leapfrog (also known as Störmer-Verlet) [31, 32],

• •

the third order Ruth3 integrator [12, 48],

• •

a second order variational but non-contact (VNC) method for forced Lagrangian systems [16] obtained by a Verlet discretization of a Lagrangian in duplicated phase space [22],

[TABLE]

and

• •

the fourth order Runge-Kutta integrator (RK4) [39].

For the comparison with the damped oscillator with forcing, the symplectic methods are extended in a natural way by additionally adding the forcing term when evaluating the acceleration component in each step.

The error plots in Figures 1–6 show a regularised relative error computed as follows: if $x_{i}$ denotes the value of the exact solution at time $t_{i}$ and $x^{*}_{i}$ is the corresponding value of the approximate solution, $\mathrm{err}_{i}=\frac{10+x^{*}_{i}}{10+x_{i}}-1$.

The simulations have been performed in python, with support from the scipy, numpy and matplotlib libraries. The plots have been generated using matplotlib, with a style imported from the seaborn library. All code is released with an MIT license and available from GitHub and Zenodo [51].

In our numerical experiments the fourth order Runge-Kutta method shows an impressive level of accuracy, and at least for this example. The main reason for choosing our method over it is when guarantees of the geometric invariants are more important than the solution accuracy. In all cases under consideration, regardless of the size of the error, our first and second order contact integrators guarantee the conservation of the contact structure, unlike any of the other methods.

As already shown in Example 2, the leapfrog method and our second order contact method are equivalent, except for the initialization of the momentum. If $x_{j}=x(jh)$ for some smooth interpolating curve $x$, then we have

[TABLE]

hence if $x^{(3)}(0)$ and $\alpha\ddot{x}(0)$ have the same sign, then it is best to initialize with $p$, i.e. use the contact method. If they are of opposite sign, which is very likely for overdamped systems, the leapfrog method will be better.

Furthermore, in the limit $\alpha\to 0$, i.e. in the limit of the system becoming symplectic, both the integrators presented in Example 1 and 2 converge to the same symplectic leapfrog scheme. The same is true for the time-dependent case of Example 6. Thus for small values of $\alpha$ we a priori expect our integrator to be on par with the leapfrog integrator, and in general perform worse than the third order Ruth3 integrator. As can be seen from Figures 1 and 4, for $\alpha=0.01$ this is indeed the case: the contact integrators are performing very similarly and Ruth3 performs much better.

One interesting fact, however, is that already for $\alpha=0.1$ our method outperforms the third order Ruth3 method. We believe that the reason for this is that the Ruth3 method is only symplectic for separable Hamiltonians, i.e. if $\dot{p}=f(q)$ and $\dot{q}=g(p)$, whereas with damping, the acceleration depends also on $p$.

Similarly, even though the contact integrator is outperforming the VNC integrator in all the simulations, we see that for $\alpha=0.01$ the VNC integrator and our contact integrators have similar performances, but the contact integrator is performing much better when $\alpha$ increases. In this case, the lack of separability does not explain the difference in performance, and we are left to believe that it is the preservation of the contact structure that makes the difference.

7 Conclusions

In this work we have begun a thorough investigation of geometric numerical integrators for contact flows. Contrary to [20], our approach is variational: we discretize Herglotz’ variational principle and obtain the discrete version of the generalized Euler-Lagrange equations (see Theorem 1). Furthermore, in Theorems 2 and 3 we have proved that the discrete map thus obtained is contact and that any geometric integrator for contact flows must be of this (variational) type.

In Theorem 4 we presented a formal backward error analysis for contact variational integrators, showing that the numerical solutions are interpolated by contact flows.

Finally, we have considered the implementation of the first and second order contact integrators for the benchmark example of a damped harmonic oscillator, both with and without external forcing. Our numerical experiments show that contact variational integrators in general are comparable with both symplectic and variational but non-contact (VNC) methods for forced Lagrangian systems of the same order, but that in situations where the contact structure is more relevant (for instance, when the damping increases), they usually outperform them.

Motivated by the results of this work, we expect to extend our analysis in multiple directions. On the one hand, we plan to derive higher order analogues of the first and second order contact methods presented here. On the other hand, we would like to compare our approach with the purely Hamiltonian integrators put forward in [20] in a number of systems. With this future work in mind, we expect that the implementation of contact integrators will be beneficial for the study of a wide range of applications where dissipation plays a central role.

Acknowledgements

The authors would like to thank the organizers of the VI Iberoamerican Meeting on Geometry, Mechanics and Control, during which part of this work was initiated, and the NWO Visitor Travel Grant 040.11.698 that sponsored the visit of AB at the Bernoulli Institute. MV is funded by the SFB Transregio 109 “Discretization in Geometry and Dynamics”. AB acknowledges FORDECYT (project number 265667) for financial support. MS research is supported by the NWO project 613.009.10. The authors would also like to thank the anonymous referees for useful comments that improved the final version of this paper.