Fractional damping through restricted calculus of variations

TL;DR

This paper introduces a new variational approach to fractional damping in mechanical systems, establishing a restricted Hamilton's principle and developing numerical integrators called Fractional Variational Integrators (FVIs) that outperform traditional methods in energy tracking.

Contribution

It proposes a novel restricted Hamilton's principle incorporating fractional derivatives, and develops FVIs for better numerical simulation of fractional damping systems.

Findings

FVIs have local truncation order 1.

FVIs outperform explicit and implicit Euler schemes in energy tracking.

The restricted Hamilton's principle accurately captures fractional damping dynamics.

Abstract

We deliver a novel approach towards the variational description of Lagrangian mechanical systems subject to fractional damping by establishing a restricted Hamilton's principle. Fractional damping is a particular instance of non-local (in time) damping, which is ubiquitous in mechanical engineering applications. The restricted Hamilton's principle relies on including fractional derivatives to the state space, the doubling of curves (which implies an extra mirror system) and the restriction of the class of varied curves. We will obtain the correct dynamics, and will show rigorously that the extra mirror dynamics is nothing but the main one in reversed time; thus, the restricted Hamilton's principle is not adding extra physics to the original system. The price to pay, on the other hand, is that the fractional damped dynamics is only a sufficient condition for the extremals of the action. In addition, we proceed to discretise the new principle. This discretisation provides a set of numerical integrators for the continuous dynamics that we denote Fractional Variational Integrators (FVIs). The discrete dynamics is obtained upon the same ingredients, say doubling of discrete curves and restriction of the discrete variations. We display the performance of the FVIs, which have local truncation order 1, in two examples. As other integrators with variational origin, for instance those generated by the discrete Lagrange-d'Alembert principle, they show a superior performance tracking the dissipative energy, in opposition to direct (order 1) discretisations of the dissipative equations, such as explicit and implicit Euler schemes.

Equations

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right)-\frac{\partial L}{\partial x}=-\rho D_{-}^{\alpha }{D}_{-}^{\alpha }x,$$

$$D_1{L}_d(x_k,x_{k+1})+D_2{L}_d(x_{k-1},x_k)=h\rho {\Delta }_{-}^{\alpha }{\Delta }_{-}^{\alpha }{x}_k,$$

$${}_{{}_{RL}}{D}_{-}^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\int }_a^t(t-\tau {)}^{-\alpha }f(\tau )d\tau ,{}_{{}_{RL}}{D}_{+}^{\alpha }f(t)=-\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\int }_t^b(\tau -t{)}^{-\alpha }f(\tau )d\tau ,$$

$${}_{{}_C}{D}_{-}^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_a^t(t-\tau {)}^{-\alpha }\dot{f}(\tau )d\tau ,{}_{{}_C}{D}_{+}^{\alpha }f(t)=-\frac{1}{\Gamma (1-\alpha )}{\int }_t^b(\tau -t{)}^{-\alpha }\dot{f}(\tau )d\tau ,$$

$${}_{{}_{RL}}{D}_{-}^{\alpha }f(t)=\frac{-1}{\Gamma (1-\alpha )}\frac{f(a)}{(t-a{)}^{\alpha }}{+}_{{}_C}{D}_{-}^{\alpha }f(t),{}_{{}_{RL}}{D}_{+}^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}\frac{f(b)}{(b-t{)}^{\alpha }}{+}_{{}_C}{D}_{+}^{\alpha }f(t).$$

$${\int }_a^{b}f(t)D_{\sigma }^{\alpha }g(t)dt$$

$$D_{\sigma }^{\alpha }{D}_{\sigma }^{\beta }$$

$$D_{-}^{\alpha }=d/dt,$$

$$\delta {\int }_a^{b}L(q(t),\dot{q}(t))dt=0,$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)-\frac{\partial L}{\partial q}=0.$$

$$\mathbb{F}L:TQ\to T^{\ast }Q;(q,\dot{q})\mapsto \left(q,p=\frac{\partial L}{\partial \dot{q}}\right).$$

$$H(q,p):=⟨p,\dot{q}⟩-L(q,\dot{q}),$$

$$\delta {\int }_a^b\left\{⟨p(t),\dot{q}(t)⟩-H(q(t),p(t))\right\}dt=0.$$

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

$$E(q,\dot{q}):=\Big{<}\frac{\partial L}{\partial\dot{q}},\dot{q}\Big{>}-L(q,\dot{q}),$$

$$f_L:TQ\to T^{\ast }Q.$$

$$\delta {\int }_a^{b}L(q,\dot{q})dt+{\int }_a^b⟨f_L(q,\dot{q}),\delta q⟩dt=0,$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)-\frac{\partial L}{\partial q}=f_L(q,\dot{q}).$$

$$\delta {\int }_a^b\left\{⟨p,\dot{q}⟩-H(q,p)\right\}dt+{\int }_a^b⟨f_H(q,p),\delta q⟩dt=0,$$

$$\dot{q}=\frac{\partial H}{\partial p},\dot{p}=-\frac{\partial H}{\partial q}+f_H(q,p).$$

$$S_d({\gamma }_d)=\sum _{k=0}^{N-1}{L}_d(q_k,q_{k+1}).$$

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

$${\mathbb{F}}^{-}{L}_d,{\mathbb{F}}^{+}{L}_d:Q\times Q\to T^{\ast }Q,$$

$${\mathbb{F}}^{-}{L}_d(q_k,q_{k+1})$$

$${\mathbb{F}}^{+}{L}_d(q_k,q_{k+1})$$

$$p_k^{-}=p_k^{+},$$

$${\tilde{F}}_{L_d}={\mathbb{F}}^{+}{L}_d\circ ({\mathbb{F}}^{-}{L}_d{)}^{-1}={\mathbb{F}}^{+}{L}_d\circ F_{L_d}\circ ({\mathbb{F}}^{+}{L}_d{)}^{-1}={\mathbb{F}}^{-}{L}_d\circ F_{L_d}\circ ({\mathbb{F}}^{-}{L}_d{)}^{-1};$$

$$f_{L_d}^{-},f_{L_d}^{+}:Q\times Q\to T^{\ast }Q$$

$$⟨f_{L_d}^{-}(q_k,q_{k+1}),\delta q_k⟩+⟨f_{L_d}^{+}(q_k,q_{k+1}),\delta q_{k+1}⟩\simeq {\int }_{t_k}^{t_k+h}⟨f_L(q,\dot{q}),\delta q⟩dt.$$

$$\delta\sum_{k=0}^{N-1}L_{d}(q_{k},q_{k+1})+\sum_{k=0}^{N-1}\Big{[}\langle f_{L_{d}}^{-}(q_{k},q_{k+1}),\delta q_{k}\rangle+\langle f_{L_{d}}^{+}(q_{k},q_{k+1}),\delta q_{k+1}\rangle\Big{]}=0.$$

$$D_1{L}_d(q_k,q_{k+1})+D_2{L}_d(q_{k-1},q_k)+f_{L_d}^{-}(q_k,q_{k+1})+f_{L_d}^{+}(q_{k-1},q_k)=0,k=1,...,N-1;$$

$${\mathbb{F}}^{-}{L}_d^{\text{mbox}f}(q_k,q_{k+1})$$

Full text

Fernando Jiménez and Sina Ober-Blöbaum

**Department of Engineering Science, University of Oxford

Parks Road, Oxford OXI 3PJ, UK

**  e-mail: fernando.jimenez@eng.ox.ac.uk

e-mail: sina.ober-blobaum@eng.ox.ac.uk

1. Introduction

The problem of obtaining a variational description of mechanical systems subject to external forces has been present in the literature for long time. Concretely, this article is concerned with the variational nature of the dynamical equations of a Lagrangian system subject to what we call fractional damping, namely:

[TABLE]

where $x$ is the dynamical variable, $L(x,\dot{x})$ is a Lagrangian function, $\rho\in\mathbb{R}_{+}$ and $D^{\alpha}_{-}$, with $\alpha\in[0,1]$, is the retarded fractional derivative defined below (3) (as usual, the dot notation represents the time derivative). The right hand side is non-local in time, and therefore the previous equation represents a particular example of non-locally damped mechanical system (we shall focus on Lagrangians given by kinetic minus potential energy), which are ubiquitous in mechanical engineering applications (see [3, 35, 38] and references therein). In addition, (1) also involves the linear damping case, since $D^{2\alpha}_{-}x=\dot{x}$ when $\alpha=1/2$, where we have used (5b),(5c). This case is a pardigmatic example of non-Lagrangian/Hamiltonian system, i.e.  its dynamics cannot be obtained from the Hamilton’s principle [2] given a Lagrangian/Hamiltonian function, as proven in [8]. A remarkable approach towards the variational modelling of external forces, due to its phenomenological versatility, is the Lagrange-d’Alembert principle [10], where the variation of the action is set equal to the work done by external forces under virtual displacements. As happens with the usual Hamilton’s principle, the Lagrange-d’Alembert’s can be performed in Lagrangian and Hamiltonian fashions, both related to each other by the Legendre transformation [2]. It is important to note, however, that Lagrange-d’Alembert principle is not variational in the pure sense on the word, circumstance that we try to avoid in our apprach.

There are previous approaches to generate a purely variational principle by including fractional derivatives into the state space of the considered Lagrangian functions; more concretely, we base our work on [14, 34]. These references delivered promising results, but they present some drawbacks. Say: in [34] Riewe did not take into account the asymmetric integration by parts of the fractional derivatives (5a); Cresson and collaborators, in [14, 15], designed the so-called asymmetric embeddings to surpass that issue, but such objects are unclear from the point of view of calculus of variations. This last approach implies the doubling of curves in the state space; i.e. it is necessary to add an extra $y$-mirror system, which is natural when treating externally forced system in a variational way, see [7, 16, 24] and references therein.

Taking as well into account a doubled space of curves, in this work we establish a novel approach to surpass the asymmetry issue while obtaining the correct fractional damped dynamics (1), embodied in a particular restriction of the usual calculus of variations. Namely, we will apply the Hamilton’s principle over a properly designed (in a geometric way) Lagrangian function, but we shall restrict the class of varied curves. We denote this principle by restricted Hamilton’s principle. Out of this principle, we shall obtain both (1) and the dynamics of the $y$-mirror system. However, we shall rigorously prove, both in Lagrangian and Hamiltonian settings, that the latter is nothing but the former in reversed time, which implies that the restricted Hamilton’s principle is not adding extra physics to the studied system. As it is made clear below, the price to pay for applying this new principle is that the dynamical equations are not anymore necessary and sufficient conditions for the extremals of the action (as in the classical Hamilton’s principle), but only sufficient. Finally, given that the linear damping case is also included, we elucidate the connection between the restricted Hamilton’s and Lagrange-d’Alembert principles.

During the last years, the discretisation of variational principles (Hamilton’s, Lagrage-d’Alembert and others [23, 29]) has been of high interest in numerical integration theory, since such discretisations generate numerical integrators that approximate the systems’ dynamics faithfully both from dynamical and geometrical perspectives, presenting as well a superior behaviour in the integration of the energy in the long-term. Thus, we proceed to discretise the restricted Hamilton’s principle, which leads to the discrete counterpart of (2), i.e.:

[TABLE]

where $x_{k}$ is the discrete dynamical variable, $L_{d}(x_{k},x_{k+1})$ the discrete Lagrangian and $\Delta_{-}^{\alpha}x_{k}$ is a given approximation of $D^{\alpha}_{-}x(t_{k})$ for a time grid $t_{k}$ with time step $h$ (other approaches to discretise fractional mechanical problems can be found in [6, 11, 13, 15] and references therein). The elements of the discrete principle are analogous to the continuous’ ones: the discrete dynamical equations are sufficient conditions for the extremals of the discrete action, and we have an extra $y$-mirror system (which again is proven to be the $x$-system in reversed discrete time). The initialisation of (2) as a numerical integrator of (1) (which we will call Fractional Variational Integator, FVI) requires an initial condition which is based on the discrete Legendre transform linking the Lagrangian and Hamiltonian approaches. With that aim, we develop the restricted Hamilton’s principle also in a Hamiltonian fashion (both in continuous and discrete scenarios), where, as in the usual discrete mechanics [23], the momentum mathching condition will be a crucial element when constructing numerical integrators approximating the Hamiltonian version of (1). As in the continuous side, we elucidate the connection between the restricted Hamilton’s and Lagrange-d’Alembert principles in the discrete side. Finally, we test the generated integrators with respect to well-known examples. With that aim, we use a particular benchmark approximation of the solution of inhomogeneous fractional-differential equations [32], based on a matrix discretisation of the fractional-differential operators and, afterwards, the resolution of the discrete dynamics as a matricial-algebraic equation [33].

The paper is organised as follows: In §2 we introduce the basics on fractional derivatives (since they will be used as a tool, we do not put emphasis on the mathematical specifics, and refer the interested reader to the proper literature); moreover we present the variational description of Lagrangian/Hamiltonian systems (Hamilton’s principle) both in the continuous and discrete settings, as well as externally forced systems (Lagrange-d’Alembert principle). §3 is devoted to develop the continuous restricted Hamilton’s principle, which is stated in its Lagrangian version in Theorem 3.1. The relationship between the $x$-system and the $y$-system in reversed time is given in Proposition 3.2. The preservation of a particular presymplectic form by (1) is established in Proposition 3.3, for which the existence and uniqueness of solutions of (1) is necessary, whose proof is developed in Appendix. The Hamiltonian version of the restricted Hamilton’s principle is given in Theorem 3.2. §4 accounts for the discrete restricted Hamilton’s principle, established in Proposition 4.1 and Theorem 4.1. The relationship between the $x$-system and the $y$-system in reversed discrete time is given in Proposition 4.2. The connection between the discrete restricted Hamilton’s principle and the discrete Lagrange-d’Alembert’s is given in Corollary 4.2, for which Lemma 4.2 is essential. The definition of the discrete Legendre transform, needed for the initialisation of the FVIs, is established in Definition 4.3, whereas the momentum matching condition and the Hamiltonian version of the FVIs is given in Proposition 4.3. Finally, in §5 we display the performance of the FVIs in two examples.

2. Preliminaries

2.1. Fractional derivatives

Let $\alpha\in[0,1]\subset\mathbb{R}$ and $f:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}$, $[a,b]\subset\mathbb{R}$, a smooth function. The fractional derivatives are defined by

[TABLE]

for $t\in$ [a,b], and $\Gamma(z)$ the gamma function, in their Riemann-Liouville and Caputo expressions. These two kinds of fractional derivatives are related to each other, indeed it can be shown that:

[TABLE]

In this work, the function $f$ will represent the dynamical variable of a mechanical system. Thus, it is interesting to remark that we can always set $f(a)=0$, i.e. the system is at the origin of coordinates at initial time, and consequently both retarded Riemman-Liouville and Caputo versions are equivalent, according to the first equation in (4)111This is particularly apparent for $\alpha=0$. Namely: ${}_{{}_{RL}}D^{0}_{-}f(t)=f(t)$, whereas ${}_{{}_{C}}D^{0}_{-}f(t)=f(t)-f(a)$., for dynamical purposes. We will assume henceforth that $D_{\pm}^{\alpha}$ is determined by the Riemann-Liouville expression (3a), unless otherwise stated. Further relevant properties are:

[TABLE]

For the proof of these properties and more details on fractional derivatives, we refer to [36].

2.2. Continuous Lagrangian and Hamiltonian description of mechanics

In this subsection we shall consider the configuration space of the studied systems as a finite dimensional smooth manifold $Q$. Moreover, $TQ$ and $T^{*}Q$ will denote its tangent and cotangent bundles, locally represented by coordinates $(q,\dot{q})$ and $(q,p)$, respectively. For more details on the geometric formulation of mechanics we refer to [2].

2.2.1. Conservative systems

Given a Lagrangian function $L:TQ\rightarrow\mathbb{R}$, the associated action functional in the time interval $t\in[a,b]\subset\mathbb{R}$ for a smooth curve $q:[a,b]\rightarrow Q$ is defined by $S(q)=\int_{a}^{b}L(q(t),\dot{q}(t))\,dt$. Through Hamilton’s principle, i.e. the true evolution of the system $q(t)$ with fixed endpoints $q(a)$ and $q(b)$ will satisfy

[TABLE]

we obtain the Euler-Lagrange equations via calculus of variations:

[TABLE]

Define the Legendre transformation:

[TABLE]

If (8) is a global diffeomorphism we say that it is hyperregular, and we call the Lagrangian function hyperregular. Under the assumption of hyperregularity (which we will take throughout the article), via the Legendre transformation we can define the Hamiltonian function $H:T^{*}Q\rightarrow\mathbb{R}$:

[TABLE]

where $\langle\cdot,\cdot\rangle:T^{*}Q\times TQ\rightarrow\mathbb{R}$ is the natural pairing. From the definition of the Hamiltonian function (9) it follows that $L(q,\dot{q})=\langle p,\dot{q}\rangle-H(q,p)$. Furthermore, from (6) we can write the stationary condition of the action functional in a Hamiltonian version, i.e.

[TABLE]

Again, using calculus of variations we obtain the Hamilton equations:

[TABLE]

We shall consider that the physical energy of the system is given by the Hamiltonian function (9). It is easy to check that $dH(q,p)/dt=0$ under (11), showing that the system is conservative. Equivalently, the Lagrangian energy

[TABLE]

is invariant under (7), i.e. $dE(q,\dot{q})/dt=0.$

2.2.2. Forced systems

First we model the external forces (which might include damping, dragging, etc.) through the mapping:

[TABLE]

The forced dynamics is provided by the Lagrange-d’Alembert principle [5, 10]: the true evolution of the system $q(t)$ between fixed points $q(a)$ and $q(b)$ will satisfy

[TABLE]

where $\delta q\in TQ$, which provides the forced Euler-Lagrange equations:

[TABLE]

Now, the Lagrangian energy of the system (12) is not preserved by (15). In particular dE(q,\dot{q})/dt=\big{<}f_{L}(q,\dot{q}),\dot{q}\big{>}, showing that this kind of systems is not conservative.

The dual version of the Lagrange-d’Alembert principle (14) is naturally obtained through the Legendre transformation (8). The dual external forces $f_{H}:T^{*}Q\rightarrow T^{*}Q$ are defined by $f_{H}:=f_{L}\circ\left(\mathbb{F}L\right)^{-1}$ (we recall that we are assuming $L$ hyperregular), while the dynamics is established by

[TABLE]

yielding the forced Hamilton equations:

[TABLE]

Obviously, the Hamiltonian function (9) is not preserved under (17). In particular dH(q,p)/dt=\big{<}f_{H}(q,p),\frac{\partial H}{\partial p}\big{>}.

2.3. Discrete Lagrangian and Hamiltonian description of mechanics

2.3.1. Conservative systems

The construction of the discrete version of mechanics relies on the substitution of $TQ$ by the Cartesian product $Q\times Q$ (note that these two spaces contain the same amount of information at local level) [23, 27]. The continuous curves $q(t)$ will be replaced by discrete ones, say $\gamma_{d}=\left\{q_{k}\right\}_{0:N}:=\left\{q_{0},q_{1},...,q_{N}\right\}\in Q^{N+1}$, where $N\in\mathbb{N}$ and the power $N+1$ indicates the Cartesian product of $N+1$ copies of $Q$. Given an increasing sequence of times $\left\{t_{k}=a+hk\,|\,k=0,...,N\right\}\subset\mathbb{R}$, with $h=(b-a)/N$, the points in $\gamma_{d}$ will be considered as an approximation of the continuous curve at time $t_{k}$, i.e. $q_{k}\simeq q(t_{k}).$ Defining the discrete Lagrangian $L_{d}:Q\times Q\rightarrow\mathbb{R}$ as an approximation of the action integral in one time step, say $L_{d}(q_{k},q_{k+1},h)\simeq\int_{t_{k}}^{t_{k}+h}L(q(t),\dot{q}(t))\,dt$ (we shall omit the $h$ dependence of the discrete Lagrangian unless needed), we can establish the so called discrete action sum:

[TABLE]

Applying the Hamilton’s principle over (18), i.e. considering variations of $\gamma_{d}$ with fixed endpoints $q_{0}=q(a),\,q_{N}=q(b)$ and extremizing $S_{d}$, we obtain the discrete Euler-Lagrange equations

[TABLE]

where $D_{1}$ and $D_{2}$ denote the partial derivative with respect to the first and second variables, respectively. If $L_{d}$ is regular, i.e. the matrix \big{[}D_{12}L_{d}\big{]} is invertible, the equations (19) define a discrete Lagrangian flow $F_{L_{d}}:Q\times Q\rightarrow Q\times Q$; $(q_{k},q_{k+1})\mapsto(q_{k+1},q_{k+2})$, which is normally called variational integrator of the continuous dynamics provided by the Euler-Lagrange equations (7) (indistinctly, we shall call the equations (19) also variational integrator). Moreover, (19) are a discretisation in finite differences of (7).

In order to establish the Hamiltonian picture we need to introduce the discrete Legendre transforms. From $L_{d}$, two of them can be defined:

[TABLE]

in particular

[TABLE]

We observe that the momentum matching condition, i.e.

[TABLE]

provides the discrete Euler-Lagrange equations (19) according to (20) (based on this, we shall refer indistinctly to the discrete Legendre transform as momentum matching). Under the regularity of $L_{d}$, both discrete Legendre transforms are invertible and the discrete Hamiltonian flow $\tilde{F}_{L_{d}}:T^{*}Q\rightarrow T^{*}Q$; $(q_{k},p_{k})\mapsto(q_{k+1},p_{k+1})$ can be defined by any of the following identities:

[TABLE]

see [23] for the proof. At the Hamiltonian level, the map $\tilde{F}_{L_{d}}$ is called variational integrator of the continuous dynamics provided by the Hamilton equations (11). Moreover, the discrete equations provided by (20) are a discretisation in finite differences of (11).

Remark 2.1**.**

Another advantage of the Hamiltonian version of variational integrators is that it provides a natural way of initialisating the numerical scheme. As showed in the previous discussion, two points $q_{0},q_{1}$ are necessary in order to start (19) and establish the discrete flow $F_{L_{d}}$. On the other hand, a mechanical problem involves as initial data $q(a)=q_{0}$, $\dot{q}(a)=v_{0}$ and $p(a)=p_{0}$; thus an extra step, providing $q_{1}$, is in order. This step is naturally determined by (20a), leading to the algorithm (where the regularity of $L_{d}$ is assumed):

Algorithm 1**.**

Variational Integrator Scheme

• 1:Initial data: $N,\,h,\,q_{0},\,p_{0}.$

2:solve for $q_{1}$ from $p_{0}=-D_{1}L_{d}(q_{0},q_{1}).$

3:Initial points: $q_{0},\,q_{1}.$

4:for $k=1:N-1$ do solve for $q_{k+1}$ from $D_{1}L_{d}(q_{k},q_{k+1})+D_{2}L_{d}(q_{k-1},q_{k})=0$

5:end for

6:Output: $(q_{2},...,q_{N}).$

end

Other definitions of $q_{1}$, different from Step 2, can lead to the instability of the discrete flow $F_{L_{d}}$. $\diamond$

A crucial feature of variational integrators is their symplecticity. If $\Omega_{T^{*}Q}$ is the canonical symplectic form on $T^{*}Q$ (which, according to Darboux theorem, can be locally written as $\Omega_{T^{*}Q}=dq\wedge dp$), define $\Omega_{Q\times Q}:=(\mathbb{F}^{-}L_{d})^{*}\Omega_{T^{*}Q}=(\mathbb{F}^{+}L_{d})^{*}\Omega_{T^{*}Q}$. Thus, the symplecticity of the variational integrators imply $F_{L_{d}}^{*}\Omega_{Q\times Q}=\Omega_{Q\times Q}$ [23, 27], which furthermore imply that the energy cannot be conserved at the same time [17]. However, symplectic integrators have proven to present a stable energy behaviour even in long-term simulations [37], behaviour that can be explained in terms of Backward Error Analysis [18, 19].

2.3.2. Forced systems

As discrete version of the external forces (13) we consider the maps:

[TABLE]

such that

[TABLE]

Note that the previous equation implies that $f_{L_{d}}^{-}(q_{k},q_{k+1})\in T^{*}_{q_{k}}Q$ and $f_{L_{d}}^{+}(q_{k},q_{k+1})\in T^{*}_{q_{k+1}}Q$. The discrete Lagrange-d’Alembert principle [23, 29] provides discrete curves between fixed $q_{0},q_{N}$ satisfying the critical condition

[TABLE]

These curves are given by the forced discrete Euler-Lagrange equations

[TABLE]

they are a discretisation in finite differences of (15) and, under the regularity of the matrix \big{[}D_{12}L_{d}(q_{k},q_{k+1})+D_{2}f_{L_{d}}^{-}(q_{k},q_{k+1})\big{]}, provide a forced discrete Lagrangian map $F_{L_{d}}^{\tiny\mbox{f}}:Q\times Q\rightarrow Q\times Q$ approximating their continuous solution. In the forced case, the discrete Legendre transformation is defined by

[TABLE]

The momentum matching condition (21) reproduces the forced discrete Euler-Lagrange equations (23). Moreover, (24) provide a discretisation in finite differences of (17), whereas (22) (for $F_{L_{d}}^{\tiny\mbox{f}}$ and $\mathbb{F}^{\pm}L_{d}^{\tiny\mbox{f}}$) yields an approximation $\tilde{F}_{L_{d}}^{\tiny\mbox{f}}:T^{*}Q\rightarrow T^{*}Q$ of their continuous flow.

Remark 2.2**.**

Establishing the forced variational flow $F_{L_{d}}^{\tiny\mbox{f}}:Q\times Q\rightarrow Q\times Q$ requires as well two initial points $q_{0},q_{1}$ as in the free case discussed in Remark 2.1. According to (24a) the algorithm for the forced variational integrator is given by:

Algorithm 2**.**

Forced Variational Integrator Scheme

• 1:Initial data: $N,\,h,\,q_{0},\,p_{0}.$

2:solve for $q_{1}$ from $p_{0}=-D_{1}L_{d}(q_{0},q_{1})-f^{-}_{L_{d}}(q_{0},q_{1}).$

3:Initial points: $q_{0},\,q_{1}.$

4:for $k=1:N-1$ do solve for $q_{k+1}$ from $D_{1}L_{d}(q_{k},q_{k+1})+D_{2}L_{d}(q_{k-1},q_{k})+f_{L_{d}}^{-}(q_{k},q_{k+1})+f_{L_{d}}^{+}(q_{k-1},q_{k})=0$

5:end for

6:Output: $(q_{2},...,q_{N}).$

end

$\diamond$

3. continuous restricted Hamilton’s principle

3.1. Fractional state space

Consider a smooth curve $\gamma:[a,b]\subset\mathbb{R}\rightarrow\mathbb{R}^{d}$ for $d\in\mathbb{N}$. The local representation of the curve is given by $\gamma(t)=(x^{1}(t),...,x^{d}(t))$, $t\in[a,b]$. For the set of all smooth curves $C^{\infty}([a,b],\mathbb{R}^{d})$, let us define the fractional tangent vector of the curve $\gamma$ by means of the following mapping

[TABLE]

where the fractional tangent vector is defined by

[TABLE]

and $D^{\alpha}_{\sigma}$ represents the $\alpha$-fractional derivatives (3), $\sigma=\left\{-,+\right\}$. In the sequel we shall omit the $t$-dependence in curves and fractional tangent vectors; furthermore, we will omit as well the coordinate superindex $i$, refering to the local expression of $X_{\sigma}^{\alpha}\gamma$ just as $D^{\alpha}_{\sigma}x$.

Remark 3.1**.**

Observe that we are choosing $\mathbb{R}^{d}$ as configuration space instead of a $d$-dimensional smooth manifold $Q$, as in §2.2. We do so because the definition of $D^{\alpha}_{\sigma}\gamma(t)$ for $\gamma(t)\subset Q$ is not unique and depends on the particular set of charts employed to cover the manifold. $\diamond$

Definition 3.1**.**

We define the fractional tangent space of the curve $\gamma$ as

[TABLE]

Proposition 3.1**.**

$V^{\alpha}_{\sigma}\mathbb{R}^{d}$* is a vector space with dim $=d$.*

Proof.

The vector structure for curves is defined pointwise, i.e. for $\gamma,\delta\in C^{\infty}([a,b],\mathbb{R}^{d})$ and $\lambda\in\mathbb{R}$

[TABLE]

Noting that $D^{\alpha}_{\sigma}$ is a linear operator according to (3), the vector structure for $V^{\alpha}_{\sigma}\mathbb{R}^{d}$ is defined by

[TABLE]

Considering that for any curve $\gamma$ the canonical basis of $\mathbb{R}^{d}$ is a linearly independent generating system with $d$ elements for $X^{\alpha}_{\sigma}\gamma$ with coordinates $D^{\alpha}_{\sigma}x^{i}$, we conclude that $V^{\alpha}_{\sigma}\mathbb{R}^{d}$ is a vector space with dim $=d$. ∎

Now, we enlarge our space of cuves in the following way: given $\gamma_{x},\gamma_{y}\in C^{\infty}([a,b],\mathbb{R}^{d})$ it is straightforward to define the curve $\tilde{\gamma}\in C^{\infty}([a,b],\mathbb{R}^{d}\times\mathbb{R}^{d})$ by $\tilde{\gamma}(t):=(\gamma_{x}(t),\gamma_{y}(t))$ for $t\in[a,b]$, which locally we shall denote $\tilde{\gamma}=(x,y)$ (which we will use indistinctly). Furthermore, for the space of curves $C^{\infty}([a,b],\mathbb{R}^{d}\times\mathbb{R}^{d})$ we establish the vector bundle $\mathbb{T}^{\alpha}\mathbb{R}^{d}$, with $V^{\alpha}_{-}\mathbb{R}^{d}\times V^{\alpha}_{+}\mathbb{R}^{d}$ (whose vector structure is straightforward after Definiton 3.1 and the Cartesian product; see [2] for further details) as the fiber over the base space $\mathbb{R}^{d}\times\mathbb{R}^{d}$. More particuarly, for a given curve $\tilde{\gamma}$, $\Omega_{\tilde{\gamma}}^{\alpha}\in\mathbb{T}^{\alpha}\mathbb{R}^{d}$ is defined by $\Omega_{\tilde{\gamma}}^{\alpha}:=((\gamma_{x},\gamma_{y}),(X^{\alpha}_{-}\gamma_{x},X^{\alpha}_{+}\gamma_{y}))$, $t\in[a,b],$ and has natural coordinates $(x,y,D^{\alpha}_{-}x,D^{\alpha}_{+}y)$ (henceforth, we shall use indistinctly $\Omega_{\tilde{\gamma}}^{\alpha}$ and its coordinates). The bundle projection $\tau^{\alpha}:\mathbb{T}^{\alpha}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$ is given by $\tau^{\alpha}(\Omega_{\tilde{\gamma}}^{\alpha})=(x,y)$; furthermore, it is apparent that $\mathbb{T}^{\alpha}\mathbb{R}^{d}\cong\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d}$. On the other hand, $\mathbb{T}\mathbb{R}^{d}$ is defined as the vector bundle with $V^{1}_{-}\mathbb{R}^{d}\times V^{1}_{-}\mathbb{R}^{d}$ as fiber over $\mathbb{R}^{d}\times\mathbb{R}^{d}$, with elements $\Xi_{\tilde{\gamma}}:=((\gamma_{x},\gamma_{y}),(X^{1}_{-}\gamma_{x},X^{1}_{-}\gamma_{y}))=((\gamma_{x},\gamma_{y}),(D^{1}_{-}\gamma_{x},D^{1}_{-}\gamma_{y}))=((\gamma_{x},\gamma_{y}),(\dot{\gamma}_{x},\dot{\gamma}_{y}))$, $t\in[a,b]$, and local coordinates $((x,y),(\dot{x},\dot{y})).$ The bundle projection is given by $\tau:\mathbb{T}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$, $\tau(\Xi_{\tilde{\gamma}})=(x,y).$ With these elements we define the fractional state space.

Definition 3.2**.**

Consider the vector bundles $\mathbb{T}^{\alpha}\mathbb{R}^{d}$ and $\mathbb{T}\mathbb{R}^{d}$. We define the fractional state space as the bundle product of them over $\mathbb{R}^{d}\times\mathbb{R}^{d}$, i.e.

[TABLE]

Thus, $\Sigma_{(x,y)}:=(\gamma_{x},\gamma_{y},X_{-}\gamma_{x},X_{-}\gamma_{y},X_{-}^{\alpha}\gamma_{x},X_{+}^{\alpha}\gamma_{y})\in\mathfrak{T}\mathbb{R}^{d}$, $t\in[a,b]$, is locally described by

[TABLE]

*The bundle projection $\mathcal{T}:\mathfrak{T}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$ is defined by $\mathcal{T}(\Sigma_{(x,y)})=(x,y)$. *

For more details on bundle products we refer to [2]. The construction of the dual bundle

[TABLE]

which we will denote fractional phase space, follows straightforwardly from the dual bundles $\mathbb{T}^{\alpha*}\mathbb{R}^{d}$ and $\mathbb{T}^{*}\mathbb{R}^{d}$. For $P_{(x,y)}\in\mathfrak{T}^{*}\mathbb{R}^{d}$, we fix local coordinates

[TABLE]

The bundle projection $\mathcal{P}:\mathfrak{T}^{*}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$ is locally given by $\mathcal{P}(P_{(x,y)})=(x,y)$. It is easy to see that both $\mathfrak{T}\mathbb{R}^{d}$ and $\mathfrak{T}^{*}\mathbb{R}^{d}$ are locally equivalent to the Cartesian product of 6 copies of $\mathbb{R}^{d}$.

3.2. Fractional dynamics

In order to establish the fractional dynamics we are going to consider a subclass of curves $C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})\subset C^{\infty}([a,b],\mathbb{R}^{d}\times\mathbb{R}^{d})$, in particular those $\tilde{\gamma}=(\gamma_{x},\gamma_{y})$ such that $\gamma_{x}(a)=x_{a},\,\gamma_{x}(b)=x_{b}$ and $\gamma_{y}(a)=y_{a},\,\gamma_{y}(b)=y_{b}$ for $x_{a},\,x_{b},\,y_{a},\,y_{b}\in\mathbb{R}^{d}$; i.e. those curves with fixed endpoints. Let us define de action sum $S:C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})\rightarrow\mathbb{R}$ for a Lagrangian function $\mathcal{L}:\mathfrak{T}\mathbb{R}^{d}\rightarrow\mathbb{R}$ (which henceforth we shall consider $C^{2}$) by

[TABLE]

and a particular set of varied curves and variations, namely:

Definition 3.3**.**

Define the set of restricted varied curves as $\Gamma_{(\eta,\epsilon)}(t):=\tilde{\gamma}(t)+\epsilon\eta(t)$, $\epsilon\in\mathbb{R}$, where $\Gamma_{(\eta,\epsilon)}(t)\in C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})$ and $\eta(t):=(\delta\gamma_{x}(t),\delta\gamma_{x}(t))$, $t\in[a,b]$, is defined such that $\delta\gamma_{x}\in C^{\infty}([a,b],\mathbb{R}^{d})$ with $\delta\gamma_{x}(a)=\delta\gamma_{x}(b)=0$. As it is easy to see, an unrestricted variation would be defined by $(\delta\gamma_{x}(t),\delta\gamma_{y}(t))$, with $\delta\gamma_{x}\neq\delta\gamma_{y}$. We impose $\delta\gamma_{x}(t)=\delta\gamma_{y}(t)$ to the variations, which locally is expressed by $\delta x=\delta y,$ vanishing at the endpoints.

Note that the introduction of $D^{\alpha}_{+}y(t)$ in the action (27) breaks its causality, since it depends on future times $[t,b]$. Thus, apparently it is not suitable for the physical description of the $x$ and $y$ systems. We will see in the following discussion, however, that their dynamics can be decoupled thanks to the restriction of the variations presented in Definition 3.3. Moreover, the choice of mechanical Lagrangians and the definition of the fractional derivatives (3a) help to surpass the causality issue.

We have already all the ingredients to establish the restricted Hamilton’s principle:

Theorem 3.1**.**

A sufficient conditions for a curve $\tilde{\gamma}\in C^{\infty}([a,b],\mathbb{R}^{d}\times\mathbb{R}^{d})$, subject to the restricted variations in Definition 3.3, to be an extremal of the action $S:C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})\rightarrow\mathbb{R}$ (27) is the so-called restricted fractional Euler-Lagrange equations:

[TABLE]

Proof.

To find the extremals of $S$ for restricted varied curves $\Gamma_{(\eta,\epsilon)}(t)$ we impose the usual critical condition, i.e. \delta S:=\frac{d}{d\epsilon}S(\Gamma_{(\eta,\epsilon)})\big{|}_{\epsilon=0}=0. Using that

[TABLE]

where

[TABLE]

according to $\Gamma_{(\eta,\epsilon)}$ in Definition 3.3, we obtain:

[TABLE]

where in the second equality we have integration by parts with respect to the total and fractional derivatives (5a). According to the vanishing endpoint conditions, all the border terms are equal to zero, leading to

[TABLE]

Finally, from this last expression and considering arbitrary $\delta x$, it is easy to see that the restricted fractional Euler-Lagrange equations (30) are a sufficient condition for $\delta S=0$; and the claim holds. ∎

Remark 3.2**.**

For unrestricted variations $\delta x\neq\delta y$, both arbitrary, and fixed endpoint conditions, the necessary and sufficient conditions for the extremals of (27) are the (unrestricted) fractional Euler-Lagrange equations:

[TABLE]

Other derivations of fractional Euler-Lagrange equations can be found in [4]. Moreover, from the linearity of the fractional derivatives (3), it is easy to prove that the restricted Euler-Lagrange equations (28) are invariant under linear constant change of variables, i.e. $x=\Lambda\,z$, $y=\Lambda\,\bar{z}$, with $\Lambda\in\mathbb{R}^{d\times d}$ a regular constant matrix. If furthermore we pick the Caputo definition of the fractional derivatives (3b), the equations are invariant under affine change of variables $x=\Lambda\,z+\zeta$, $y=\Lambda\,\bar{z}+\bar{\zeta}$; for constant $\zeta,\bar{\zeta}\in\mathbb{R}^{d}$ (see [20] for more details). $\diamond$

Define now the Lagrangian $\mathcal{L}:\mathfrak{T}\mathbb{R}^{d}\rightarrow\mathbb{R}$ by

[TABLE]

where $L:T\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a $C^{2}$ function, $\rho=\,$diag$(\rho_{1},\cdots,\rho_{d})\in\mathbb{R}^{d\times d}$ and we represent the matrix product $\mu^{T}\,M\,\nu$, for $\mu,\nu\in\mathbb{R}^{d}$ and $M\in\mathbb{R}^{d\times d}$, by $\mu\,M\,\nu$. In this case, the restricted fractional Euler-Lagrange equations (28) read:

[TABLE]

which are the usual Euler-Lagrange equations for a Lagrangian system plus a fractional damping term (1). The following proposition provides an interpretation of the $y$-mirror system.

Proposition 3.2**.**

If we set $y(t):=x(a+b-t)$, which implies that $\tilde{\gamma}\in C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})$ with $\gamma_{y}(a)=x_{b},\,\gamma_{y}(b)=x_{a}$, and we consider an even Lagrangian function $L:T\mathbb{R}^{d}\rightarrow\mathbb{R}$ in the second variable, i.e. $L(z,-\dot{z})=L(z,\dot{z})$, then (32b) is (32a) in reversed time.

Proof.

Defining the reverse time by $\tilde{t}:=a+b-t$ for $t\in[a,b]$, it is easy to see that $\dot{y}(t)=-x^{\prime}(\tilde{t})$ by applying the chain rule, where we consider $x^{\prime}(\tilde{t})=dx(\tilde{t})/d\tilde{t}$. Then

[TABLE]

where in $=^{{}_{1}}$ we have used the parity $L(x,-x^{\prime})=L(x,x^{\prime})$, and in $=^{{}_{2}}$, $=^{{}_{3}}$, the chain rule and the redefinition of the time $\tilde{t}=a+b-t$ and parity again, respectively. On the other hand, using (5b) and according to (3a) and $y(t)=x(a+b-t)$:

[TABLE]

where in $=^{{}_{1}}$ we have used the change of variables $\tilde{\tau}=a+b-\tau$, in $=^{{}_{2}}$ the redefinition of time and, in $=^{{}_{3}}$, the definition of $D^{\alpha}_{-}$ in (3a). Plugging all these elements in (32b) we obtain

[TABLE]

for $\tilde{t}:b\rightarrow a$, and the claim holds. ∎

Thus, our conclusion is that, for the (mechanical) Lagrangians of interest in this work, i.e.

[TABLE]

with $m=\,$diag $(m_{1},\cdots,m_{d})\in\mathbb{R}^{d\times d}_{+}$ and $U:\mathbb{R}^{d}\rightarrow\mathbb{R}$ a smooth function, we can consider the $y$-mirror system as the $x$-system in reversed time, as long as $\gamma_{y}(t)$ is bounded by $y(a)=x_{b}$ and $y(b)=x_{a}$. Thus, doubling the space of curves as exposed in §3.1 allows us to surpass the issue raised by the asymetric integration by parts (5a) without adding extra unnecessary (physical) information to the dynamics. In addition, the non-causal terms in the action (27) due to the presence of the advanced fractional derivative become causal in reversed time since $D^{\alpha}_{+}y(t)=D^{\alpha}_{-}x(\tilde{t})$, as we just proved. Finally, we observe that the restricted varied curves introduced in Definition 3.3 acquire a particular shape when we set $y(t)=x(a+b-t)$, as we show in Figure 1.

Remark 3.3**.**

Note that the (unrestricted) fractional Euler-Lagrange equations (30) for the Lagrangians (31) read

[TABLE]

which are a coupled system of fractional differential equations with meaningless dynamics. $\diamond$

In addition, $D^{\alpha}_{-}D^{\alpha}_{-}=D^{2\alpha}_{-}$ (5b); therefore $\alpha=1/2$ implies $D^{2\alpha}_{-}=d/dt$ (5c). According to this, it is easy to see that, given $\alpha=1/2$, then equation (32a) is equivalent to the forced Euler-Lagrange equations (15) when $f_{L}(x,\dot{x})=-\rho\,\dot{x}$. This establishes the relationship between the restricted Hamilton’s principle and the Lagrange-d’Alembert principle, expressed in Figure 2.

Picking the mechanical Lagrangian (33), equations (32) read

[TABLE]

which are the dynamical equations of a Lagrangian system with $2\alpha$ fractional damping, and will be considered in the examples, §5.

3.3. Preservation of presymplectic form

Before proving the preservation of a given presymplectic form in $\mathcal{T}\mathbb{R}^{d}$ we point out that the equations (32) define a continuous flow $F_{t}:\mathcal{T}\mathbb{R}^{d}\rightarrow\mathcal{T}\mathbb{R}^{d}$, with $F_{0}=$Id. The existence of this flow is ensured by the proof of existence and uniqueness of the cited equations in Appendix. In particular, it is proven that given initial conditions $(x_{0},\dot{x}_{0})$, the existence and uniqueness of smooth $(x(t),\dot{x}(t))$ is ensured, ensuring furthermore the existence and uniqueness of $D_{-}^{\alpha}x(t)$ according to definition (3) (equivalently for the mirror $y$-system). We shall ignore issues related to global versus local flows, which are easily dealt with by restricting the domains of the flows. Given initial conditions $X_{0}^{\alpha}:=(x_{0},y_{0},\dot{x}_{0},\dot{y}_{0},0,y^{\alpha}_{0})$222Observe that, according to (3), the only available value of $D^{\alpha}_{-}x(t_{0})$ is 0., the flow provides $F_{t}(X_{0}^{\alpha})=X^{\alpha}(t)$, with $X^{\alpha}(t):=(x(t),y(t),\dot{x}(t),\dot{y}(t),D^{\alpha}_{-}x(t),D^{\alpha}_{+}y(t)).$

Now, define the one-form and two-form on $\mathcal{T}\mathbb{R}^{d}$, with $L:T\mathbb{R}^{d}\rightarrow\mathbb{R}$ a $C^{2}$ function:

[TABLE]

It is straightforward to see that $\Omega$ is presymplectic on $\mathcal{T}\mathbb{R}^{d}$, since the fractional variables are absent, therefore it is degenerate, and $\mbox{d}\Omega=\mbox{d}\mbox{d}\Theta=0.$ We can establish the following preservation result.

Proposition 3.3**.**

The flow generated by the restricted fractional Euler-Lagrange equations when we set the Lagrangian function (31), i.e.  equations (32), preserves the presymplectic form $\Omega$ (35b). In other words, $F_{t}^{*}\Omega=\Omega.$

Proof.

The action (27), in our fractional context, can be seen as a function $S:\mathcal{T}\mathbb{R}^{d}\rightarrow\mathbb{R}$, and its variation $\delta S=\langle\mbox{d}S\,,\,\delta X^{\alpha}_{0}\rangle$, where we set $\delta X_{0}^{\alpha}:=(\delta x_{0},\delta y_{0},\delta\dot{x}_{0},\delta\dot{y}_{0},\delta x^{\alpha}_{0},\delta y^{\alpha}_{0})$ $\in T\mathcal{T}\mathbb{R}^{d}$ arbitrary. According to (3.2), when we pick the Lagrangian (31) and along the solutions of (32), we have that

[TABLE]

where we have used arbitrary initial and final times $t_{0},t$, respectively. In $=^{{}_{1}}$ we have employed the definition of $\Theta$ (35a); in $=^{{}_{2}}$ we have used the flow $F_{t}:\mathcal{T}\mathbb{R}^{d}\rightarrow\mathcal{T}\mathbb{R}^{d}$ introduced above. It follows straightforwardly that

[TABLE]

Now, taking the differential in both sides, considering that the differential and the pull-back commute and the definition of $\Omega$ (35b), we obtain $(F_{t})^{*}\mbox{d}\Theta-\mbox{d}\Theta=0$ and then $F_{t}^{*}\Omega=\Omega,$ from which the result follows. ∎

3.4. The Legendre transformation

Let us define the fractional Legendre transform $\mathcal{F}\mathcal{L}:\mathfrak{T}\mathbb{R}^{d}\rightarrow\mathfrak{T}^{*}\mathbb{R}^{d}$ as the fiber derivative for a Lagrangian function $\mathcal{L}:\mathfrak{T}\mathbb{R}^{d}\rightarrow\mathbb{R}$, i.e.

[TABLE]

where $D_{(x,y)}$ denotes the partial derivative in the fiber $\mathcal{T}^{-1}((x,y))$. Locally we have

[TABLE]

It is easy to check that $\mathcal{F}\mathcal{L}$ is fiber preserving [2]. Moreover, we will say that $\mathcal{F}\mathcal{L}$ is regular if it is a diffeomorphism, and furthermore we will call $\mathcal{L}$ regular if that is the case (which we will assume throughout the article). Hence, we define the Hamiltonian function $\mathcal{H}:\mathfrak{T}^{*}\mathbb{R}^{d}\rightarrow\mathbb{R}$ by

[TABLE]

where the coordinates of $P_{(x,y)}:=\mathcal{F}\mathcal{L}(\Sigma_{(x,y)})$ are given in (26) and $\langle\cdot,\cdot\rangle:\mathfrak{T}_{(x,y)}^{*}\mathbb{R}^{d}\times\mathfrak{T}_{(x,y)}\mathbb{R}^{d}\rightarrow\mathbb{R}$ denotes the natural pairing. Employing these elements, we can establish the Hamiltonian version of the restricted Hamilton’s principle.

Theorem 3.2**.**

A sufficient condition for a curve $\tilde{\gamma}\in C^{\infty}([a,b],\mathbb{R}^{d}\times\mathbb{R}^{d})$, subject to the restricted variations in Definition 3.3, to be an extremal of the action $S:C^{\infty}(\tilde{\gamma}^{{}_{(a,b)}};\mathbb{R}^{d}\times\mathbb{R}^{d})\rightarrow\mathbb{R}$ (27) is the so-called restricted fractional Hamilton equations:

[TABLE]

Proof.

Setting the action (27) in its Hamiltonian form, i.e.

[TABLE]

and imposing the critical condition with restricted variations, i.e. \delta S=\frac{d}{d\epsilon}S(\Gamma_{(\eta,\epsilon)})\big{|}_{\epsilon=0}=0, after applying fractional and total integration by parts we arrive to

[TABLE]

which leads to

[TABLE]

since the boundary terms vanish. From this last expression is straightforward to see that (39) is a sufficient condition for the extremals given arbitrary variations $\delta x,\,\delta p_{x},\,\delta p_{y},\,\delta p^{\alpha}_{x},\,\delta p^{\alpha}_{y}$. ∎

Remark 3.4**.**

For unrestricted variations $\delta x\neq\delta y$, both arbitrary, and fixed endpoint conditions, the necessary and sufficient conditions for the extremals of (27) are the (unrestricted) fractional Hamilton equations:

[TABLE]

$\diamond$

From (31) and (38) we get the Hamiltonian

[TABLE]

where $H(z,p_{z}):=\langle p_{z},\,\dot{z}\rangle-L(z,\dot{z})$. For the particular Hamiltonian (40), the restricted Hamilton equations (39) read:

[TABLE]

We see that the $x$ and $y$ dynamics are also coupled in principle. On the other hand, in the next result we show that the relationship between $x$ and $y$ systems is consistent with the Lagrangian case under inversion of time.

Proposition 3.4**.**

If we set $y(t):=x(a+b-t)$, $p_{y}(t):=-p_{x}(a+b-t)$ and $p_{y}^{\alpha}(t):=p_{x}(a+b-t)$, and we consider an even Hamiltonian function in the second variable, i.e. $H(z,-p_{z})=H(z,p_{z})$, then (41b) is (41a) in reversed time.

Proof.

Given that $\dot{y}(t)=-x^{\prime}(\tilde{t})$ and

[TABLE]

where in $=^{{}_{1}}$ we use the hypotheses, in $=^{{}_{2}}$ the chain rule and in $=^{{}_{3}}$ the parity of the Hamiltonian; we see that

[TABLE]

From the proof of Proposition 3.2, $D^{\beta}_{+}y(t)=D^{\beta}_{-}x(\tilde{t})$, and therefore it is easy to see that

[TABLE]

Finally, considering

[TABLE]

where again we use the chain rule and the parity of the Hamiltonain, and that $\dot{p}_{y}(t)=p_{x}^{\prime}(\tilde{t})$ according to the hypotheses, we get

[TABLE]

and the claim holds. ∎

The mechanical Hamiltonian

[TABLE]

presents even parity in the second variable, and therefore the previous result applies. The restricted fractional Hamilton equations (41) in this case, after replacing the fractional relationship in the dynamical equation for $p_{x}$ and $p_{y}$, read

[TABLE]

where we recognize the usual Hamilton equations for mechanical systems plus a fractional damping term. We observe that, after getting rid of the pure fractional equation, the dynamics of $x$ and $y$ are again decoupled. Furthermore, when $\alpha=1/2$ it is easy to see that the $x$-equations (analogous arguments can be applied to the $y$-equations) are equivalent to the forced Hamilton equations (17) when $f_{H}(x,p_{x})=-\rho\,m^{-1}p_{x}$. Indeed, when $\alpha=1/2$:

[TABLE]

which are the usual Hamilton equations for mechanical systems with linear damping.

Remark 3.5**.**

We observe in the right hand side of $=^{{}_{2}}$ that the damping force is not actually defined in $T^{*}\mathbb{R}^{d}$, but in $T\mathbb{R}^{d}$, i.e. $f_{L}(x,\dot{x})=-\rho\,\dot{x}$, and that we can relate it to a pure “cotangent” force thanks to $=^{{}_{1}}$, say $f_{H}(x,p_{x})=-\rho\,m^{-1}p_{x}$. This is a specific phenomenon of our approach, and differs from the usual description of forced systems §2.2.2. However, we observe as well that both dynamics, i.e.

[TABLE]

define the same subspace in $TT^{*}\mathbb{R}^{d}$, whose natural coordinates are $(x,p_{x},\dot{x},\dot{p}_{x}).$ $\diamond$

Remark 3.6**.**

Define the one-form on $\mathfrak{T}^{*}\mathbb{R}^{d}$ given by $\Theta_{H}:=(\mathcal{F}\mathcal{L})^{*}\Theta$, with $\Theta$ defined in (35a) and $\mathcal{F}\mathcal{L}$ in (37). Using the Lagrangian (31) and the coordinates (26), it is easy to see that $\Theta_{H}=(p_{x}+p_{y})\,dx$. Using analogous arguments than in Proposition 3.3, it is easy to see that the flow generated by equations (41) preserves the presymplectic form $\Omega_{H}:=\mbox{d}\Theta_{H}=dx\wedge dp_{x}+dx\wedge dp_{y}$. $\diamond$

4. Discrete restricted Hamilton’s principle

4.1. Discrete Lagrangian dynamics

Let us consider the increasing sequence of times $\{t_{k}=a+hk\,|\,k=0,...,N\}\subset\mathbb{R}$ where $h$ is the fixed time step determined by $h=(b-a)/N$. Define a discrete curve as a collection of points in $\mathbb{R}^{d}$ i.e. $\gamma_{d}:=\left\{x_{0},x_{1},...,x_{N-1},x_{N}\right\}=\left\{x_{k}\right\}_{0:N}\in\mathbb{R}^{(N+1)d}$ (here $\mathbb{R}^{(N+1)d}$ denotes $(\mathbb{R}^{d})^{(N+1)}:=\mathbb{R}^{d}\times\cdots\times\mathbb{R}^{d}$, the Cartesian product of $(N+1)$ copies of $\mathbb{R}^{d}$). As usual, we will consider these points as an approximation of the continuous curve at time $t_{k}$, namely $x_{k}\simeq x(t_{k}).$ Given $\left\{z_{k}\right\}_{0:N}$ (later on we shall particularise in $\left\{x_{k}\right\}_{0:N}$ and $\left\{y_{k}\right\}_{0:N}$) define the following sequences:

[TABLE]

where

[TABLE]

The discrete series $\Delta_{-}^{\alpha}x_{k}$ (resp. $\Delta_{+}^{\alpha}y_{k}$) are an approximation of $D^{\alpha}_{-}x(t_{k})$ (resp. $D^{\alpha}_{+}y(t_{k})$). For more details on the approximation of fractional derivatives we refer to ([15], Chapter 5).

Remark 4.1**.**

We observe that (45b), (45c) are well-defined for $k=0$ and $k=N$. Namely, straightforward computations lead to $\Delta_{-}^{\alpha}x_{0}=\alpha_{0}x_{0}/h^{\alpha}$ and $\Delta_{+}^{\alpha}y_{N}=\alpha_{0}y_{N}/h^{\alpha}$. $\diamond$

Given two sequences $\left\{F_{k}\right\}_{\tiny 0:N},\left\{G_{k}\right\}_{\tiny 0:N}$, the discrete fractional derivatives (45b), (45c), obey the following discrete integration by parts relationships:

Lemma 4.1**.**

Consider $\left\{F_{k}\right\}_{\tiny 0:N},\left\{G_{k}\right\}_{\tiny 0:N}$, with $F_{0}=F_{N}=G_{0}=G_{N}=0$. Then

[TABLE]

Proof.

(47a) (see [11] for more details):

[TABLE]

In $=^{{}_{1}}$ the definition (45b) is used. In $=^{{}_{2}}$ $F_{0}=0$ is taken into account. To prove $=^{{}_{3}}$ is enough to notice that, for a fixed $j=0,...,N,$ the elements $a_{i}:=\alpha_{i}G_{j}F_{j-i}$, $i=0,...,j$, on the left hand side, disposed in columns, form an upper diagonal $(N+1)\times(N+1)$, whereas the same elements on the right hand side, for $j=0,...,N$ and $i=j,...,N$, account for the transposed matrix; therefore their total sums are equal. In $=^{{}_{4}}$ the sum index is rearranged. In $=^{{}_{5}}$ equivalent arguments to $=^{{}_{3}}$ can be used. In $=^{{}_{6}}$ $F_{N}=0$ is taken into account. Finally, in $=^{{}_{7}}$ the definition (45c) is used.

(47b):

[TABLE]

In $=^{{}_{1}}$ we have arranged the sum index. In $=^{{}_{2}}$ we have used (47a) and, finally, in $=^{{}_{3}}$ we have used $F_{0}=0.$ ∎

Remark 4.2**.**

By simple inspection it is easy to check that the discrete asymetric integration by parts does not hold for out of phase indices, for instance $\sum_{k=0}^{N-1}(\Delta^{\alpha}_{+}G_{k})F_{k+1}\neq\sum_{k=1}^{N-1}G_{k}(\Delta^{\alpha}_{-}F_{k+1})$ and similar cases. $\diamond$

Next, according to §2.3, we shall consider the discrete Lagrangian as an approximation of the continuous action (27), i.e.

[TABLE]

for $a+h<t_{k}<b-h$. As mentioned above, we shall use (45b),(45c) as discrete counterparts for the fractional derivatives. We see that these discrete series imply the whole discrete past for $x$ and the whole discrete future for $y$ at time $t_{k}$. Thus, it is clear that the approximation of the action in the discrete time interval $[k,k+1]$ would depend on $x_{(0,k+1)}:=(x_{0},...,x_{k+1})\in\mathbb{R}^{d(k+2)}$ and $y_{(k,N)}:=(y_{k},...,y_{N})\in\mathbb{R}^{d(N+1-k)}$, accounting for the discrete version of the non-locality of the fractional derivatives. Moreover, it is explicit that the approximation of the action shall depend as well on the interval $[k,k+1]$ and therefore on $k$. According to this, we establish the function

[TABLE]

where we remark the correspondence of the first $k+2$ entries with the discrete $x$-path, and the last $N+1-k$ with the $y$-path. If we define the sets of discrete curves as

[TABLE]

then the discrete action sum is defined naturally by $S_{d}:C_{d}^{x}\times C_{d}^{y}\rightarrow\mathbb{R}$:

[TABLE]

where $\tilde{\gamma}_{d}:=(\gamma_{d}^{x},\gamma_{d}^{y})\in C_{d}^{x}\times C_{d}^{y}$. Next, we introduce the discrete restricted variations.

Definition 4.1**.**

Given a discrete curve $\tilde{\gamma}_{d}=(\gamma_{d}^{x},\gamma_{d}^{y})\in C_{d}^{x}\times C_{d}^{y}$, we define the set of varied discrete curves by

[TABLE]

where $\delta\gamma_{d}^{x}:=\left\{\delta x_{k}\right\}_{0:N}$, $\delta\gamma_{d}^{y}:=\left\{\delta y_{k}\right\}_{0:N}$ are the discrete variations, defined such that

[TABLE]

We define the set of restricted varied discrete curves, by

[TABLE]

where we establish $\eta_{d}=\delta\gamma_{d}^{x}=\delta\gamma_{d}^{y}$. In other words, we set $\delta x_{k}=\delta y_{k}$ for $k=1,...,N-1$.

In the following proposition we establish the necessary and sufficient condition for the extremals of (49) under restricted variations (discrete restricted Hamilton’s principle).

Proposition 4.1**.**

Given the discrete action (49) subject to restricted variations (52) and endpoint conditions (51), the necessary and sufficient conditions for the extremals are

[TABLE]

for $k=1,...,N-1$, where $D_{x_{k}}:=\partial/\partial x_{k}$ and $D_{y_{k}}:=\partial/\partial y_{k}$.

Proof.

The condition for the extremals is \frac{d}{d\epsilon}S_{d}((\Gamma^{x}_{(\epsilon,\eta)},\Gamma^{y}_{(\epsilon,\eta)}))\Big{|}_{\epsilon=0}=0, which, as it is easy to see, is equivalent to

[TABLE]

given that $\delta x_{k}=\delta y_{k}$. Considering the variations with respect to $x$ and $y$ separately, we obtain

[TABLE]

where $D_{x_{i}}$ represents the partial derivative with respect to $x_{i}$, and

[TABLE]

Adding up the columns, we arrive to

[TABLE]

Equating the last expression to 0, considering $\delta x_{0}=\delta x_{N}=0$ and arbitrary variations $\delta x_{k}$, we arrive directly to (53). ∎

Remark 4.3**.**

Observe that the discrete Lagrangian problem established by (48) and (49) is of the higher-order type, i.e. the discrete Lagrangian depends on multiple copies (more than two) of the configuration manifold (see [9, 12] for more details). As striking difference with respect to the kind of problem in these references, we note that the number of copies of $\mathbb{R}^{d}_{X}$ and $\mathbb{R}^{d}_{Y}$ the discrete Lagrangian (48) depends on is not fixed and is determined by $k$. This circumstance prevents in general the definition of a discrete flow in the sense expressed in §2.3, where $F_{L_{d}}:Q\times Q\rightarrow Q\times Q$ can be obtained from (19) under regularity conditions. This can be seen as well noticing that all variables are present at the same time in (53) for any $k$, which makes mandatory to solve them simultaneously in order to obtain the sequences $\left\{x_{k}\right\}_{0:N}$, $\left\{y_{k}\right\}_{0:N}$. However, we will see below that particular choices of $\mathcal{L}_{d}^{k}$ lead to actual discrete flows. $\diamond$

Remark 4.4**.**

For unrestricted variations (50), i.e. $\delta x_{k}\neq\delta y_{k}$ and both arbitrary, and endpoint conditions (51), the necessary and sufficient conditions for the extremals of (49) read

[TABLE]

Next, let us pick the discrete Lagrangian

[TABLE]

where $L_{d}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a particular discretisation of the continuous action integral defined for the Lagrangian $L:T\mathbb{R}^{d}\rightarrow\mathbb{R}$ in (31), the discrete fractional derivatives are defined in (45b), (45c) and $\rho\in\mathbb{R}^{d\times d}$. Furthermore, we set $\mathcal{L}_{d}^{N}:=0$ since $L_{d}$ is not defined in such a case. Naturally, this is the discrete counterpart of the continuous Lagrangian (31).

Theorem 4.1**.**

Given the discrete Lagrangian (54), restricted variations (52) and endpoint conditions (51), a sufficient condition for the extremals of (49) is

[TABLE]

for $k=1,...,N-1$, and $D_{i}$ denotes the partial derivative with respect to the $i$-th variable.

Proof.

To prove the claim it is enough to show that equations (55) are a sufficient condition of (53) to be satisfied.

First, we see that

[TABLE]

In (56a) we use that $\sum_{i=0}^{k-2}\partial\mathcal{L}_{d}^{i}/\partial x_{k}=0$ because there is no dependence of $\mathcal{L}_{d}^{i}$ on $x_{k}$ in the range $[0,k-2]$, and that $\mathcal{L}_{d}^{N}=0$ in the last equality (which implies that we have already picked the Lagrangian (54)). Equivalent arguments lead to $\sum_{i=k+1}^{N-1}\partial\mathcal{L}_{d}^{i}/\partial y_{k}=0$ in order to prove (56b).

Recalling that $k=1,...,N-1$, from (56a) we get:

[TABLE]

where $\delta_{ij}$ represents the Kronecker delta. In the previous computation, right hand side of $=^{{}_{1}}$, we have taken into account that

[TABLE]

given that $\partial(\Delta_{-}^{\alpha}x_{0}\,\rho\,\Delta_{+}^{\alpha}y_{0})/\partial x_{k}=0$ for $k=1,...,N-1$. Moreover, in the right hand side of $=^{{}_{2}}$ we have used (47a), $\partial(x_{0}\,\rho\,\Delta_{+}^{\alpha}\Delta_{+}^{\alpha}y_{0})/\partial x_{k}=0$ and $\partial x_{i}/\partial x_{k}=\delta_{ij}$ for $k$ in the range $[1,N-1]$.

Using equivalent arguments, it can be proven from (56b) that:

[TABLE]

where again (47a) and the range $k=1,...,N-1$ are used.

Adding both terms together and equating the sum to 0, which accounts for (53), we obtain:

[TABLE]

Swapping the fractional terms and equating each block to 0, we directly obtain (55) as a sufficient condition for (53) to hold. ∎

We note that equations (55) are the natural discrete version of (32).

Remark 4.5**.**

Observe that a term $h\Delta_{-}^{\alpha}x_{k+1}\Delta_{+}^{\alpha}y_{k+1}$ is admissible in (54) according to the definition (48), leading to a term $h\sum_{k=0}^{N-1}\Delta_{-}^{\alpha}x_{k+1}\Delta_{+}^{\alpha}y_{k+1}$ in the action sum (49). However, it provides the same discrete dynamics as $\Delta_{-}^{\alpha}x_{k}\Delta_{+}^{\alpha}y_{k}$, according to (47b), which makes it redundant. Further terms, as those described in Remark 4.2, are meaningless since the asymetric integration by parts is not defined for them. $\diamond$

In the next result we prove that the $x$ and $y$ dynamics in (55) are also related under inversion of time at a discrete level.

Proposition 4.2**.**

Given $y_{k}:=x_{N-k}$, then (55b) is (55a) in reversed discrete time.

Proof.

We define the reversed discrete time as $\tilde{k}:=N-k$, such that $\tilde{k}=N,...,0$, for $k=0,...,N$. Given that, we observe that, under inversion of time:

[TABLE]

On the other hand:

[TABLE]

Multiplying $\Delta_{-}^{\alpha}\Delta_{-}^{\alpha}x_{\tilde{k}}$ by $h\,\rho$, adding it to $D_{1}L_{d}(x_{\tilde{k}},x_{\tilde{k}+1})\,+\,D_{2}L_{d}(x_{\tilde{k}-1},x_{\tilde{k}})$ and equating the sum to 0, the claim holds. ∎

As advanced in Remark 4.3, in spite all variables are present at the same time in (53), which a priori makes necessary to solve them all simultaneously, we can find particular and meaningful Lagrangian functions (54) such that the discrete dynamics (55a) (we restrict ourselves to the $x$-system since we have just proved that $y$ can be interpreted as $x$ in reversed time) is provided by a discrete flow. That is established in the following algorithm, accounting for the definition of the FVIs:

Algorithm 3**.**

Fractional Variational Integrator Scheme

• 1:Initial data: $N,\,h,\,\alpha,\,\rho,\,x_{0},\,p_{x_{0}}.$

2:solve for $x_{1}$ from $p_{x_{0}}=I(x_{0},x_{1}).$

3:Initial points: $x_{0},\,x_{1}.$

4:for $k=1:N-1$ do solve for $x_{k+1}$ from $D_{1}L_{d}(x_{k},x_{k+1})+D_{2}L_{d}(x_{k-1},x_{k})-h\,\rho\,\Delta_{-}^{\alpha}\Delta_{-}^{\alpha}x_{k}=0$

5:end for

6:Output: $(x_{2},...,x_{N}).$

end

Observe that the initialisation condition in Step 2, i.e. $p_{x_{0}}=I(x_{0},x_{1})$, has to be properly determined. This will be object of discussion when defining the fractional discrete Legendre transformation in §4.3.

4.2. Discrete mechanical Lagrangian

Now, let us pick the discrete Lagrangian

[TABLE]

where $S^{\kappa}z_{k}$ is given in (45a), as the usual discretisation of the mechanical Lagrangian (33). In this case, (55) read

[TABLE]

where we have divided both sides by $h$. We observe that these equations are a discretisation in finite differences of (34).

According to what happens in the continous case, i.e. $D^{2\alpha}_{-}=d/dt$ when $\alpha=1/2$, we expect a particular discretisation of the total time derivative from the term $\Delta_{-}^{\alpha}\Delta_{-}^{\alpha}x_{k}$. This is proven in the following result.

Lemma 4.2**.**

For $\alpha=1/2$ and $k=1,...,N-1$,

[TABLE]

Proof.

According to (46) we have that $\alpha_{0}=1$ and $\alpha_{1}=-1/2$ for $\alpha=1/2$, leading, after expanding the summations, to

[TABLE]

In this expansion, the value of the sum when $k=0$, explicited in Remark 4.1, has been taken into account. The claim automatically holds for $k=1$. For $k\geq 2$, arranging the sum indices we see that the previous expression can be rewritten as

[TABLE]

where, for a fixed $k$, we set $r=k-1$ and

[TABLE]

(it is apparent that $j$ is not a power but a superindex). For a fixed $k=\tilde{k}$, (59) acquires the form

[TABLE]

According to this, it is enough to prove that $\beta_{0}^{j}=0$ for any $j$, for which we proceed by induction. From (60) and (46), it follows that $\beta^{2}_{0}=2\alpha_{2}+\alpha_{1}\alpha_{1}$, which vanishes for $\alpha=1/2$. Taking this as the first induction step, it is enough to prove that $\beta^{j+1}_{0}=0$ assuming that $\beta^{j}_{0}=0$. This is shown next:

[TABLE]

where we have set $r=j+1.$ Hence the claim follows. ∎

Using similar arguments, one can prove that

[TABLE]

It follows straightforwardly that

[TABLE]

showing that $\Delta_{-}^{{}_{1/2}}\Delta_{-}^{{}_{1/2}}x_{k}$, $\Delta_{+}^{{}_{1/2}}\Delta_{+}^{{}_{1/2}}y_{k}$ are the backward and forward (up to a minus sign) difference operators, respectively; thus order one approximations of the velocity. This leads to the following corollary of Theorem 4.1.

Corollary 4.2**.**

*Given $\alpha=1/2$, then (55a) is equivalent to the forced discrete Euler-Lagrange equations (23) when $f^{-}_{L_{d}}(x_{k},x_{k+1})=0$, $f^{+}_{L_{d}}(x_{k},x_{k+1})=-\rho\,(x_{k+1}-x_{k}).$ *

Proof.

The result follows straightforwardly by replacing (61) (first relationship) in (55a) and comparing with (23) for $f^{-}_{L_{d}}(x_{k},x_{k+1})=0$, $f^{+}_{L_{d}}(x_{k},x_{k+1})=-\rho\,(x_{k+1}-x_{k}).$ ∎

An equivalent result can be obtained for the $y$-mirror system setting $f_{L_{d}}^{-}(y_{k},y_{k+1})=\rho(y_{k+1}-y_{k})$ and $f_{L_{d}}^{+}(y_{k},y_{k+1})=0$. This discussion makes explicit the relationship between the discretisation of the Lagrange-d’Alembert principle (§2.3.2, (23)) and the restricted Hamilton’s principle developed in this work, that we display in Figure 3 (where we omit the $y$-system for sake of simplicity).

4.3. Discrete Legendre transformation

The main guidelines to construct the discrete Legendre transformation in the fractional case are the following:

As in the usual discrete mechanics §2.3, we want to reproduce the discrete Lagrangian dynamics (55) through momentum matching.

• 2.

We seek for a fair discretisation of the fractional Hamilton equations in the case of mechanical Hamiltonians (41).

• 3.

We intend to obtain initialisation condition (Step 2) for Algorithm 3.

According to this, we provide the following definition of the discrete Legendre transformation. Previously, we introduce the intermediate variables:

[TABLE]

and replace them into the Lagrangian (54) in order to establish the following definition

[TABLE]

For the sake of generality, we allow the presence of $x_{k+1}^{\alpha}=\Delta_{-}^{\alpha}x_{k+1}$ and $y_{k+1}^{\alpha}=\Delta_{+}^{\alpha}y_{k+1}$, which are admissible as discussed in Remark 4.5.

Definition 4.3**.**

Given the discrete Lagrangian (62), we define the discrete Legendre transformation by:

[TABLE]

where we consider the row matrix $[\Delta_{-}^{\alpha}\,\,\Delta_{+}^{\alpha}]$ in the sense of operators.

Proposition 4.3**.**

Given the discrete Legendre transformation in Definition (4.3), the particular discrete Lagrangian $\tilde{\mathcal{L}}_{d}^{k}=L_{d}(x_{k},x_{k+1})+L_{d}(y_{k},y_{k+1})-h\,x_{k+1}^{\alpha}\,\,\rho\,\,y_{k+1}^{\alpha}$ and the intermediate variables $x_{k}^{\alpha}:=\Delta^{\alpha}_{-}x_{k},\,\,y_{k}^{\alpha}:=\Delta^{\alpha}_{+}y_{k},$ the following statements are true:

The momentum matching condition $p_{x_{k}}^{-}=p_{x_{k}}^{+}$, $p_{y_{k}}^{-}=p_{y_{k}}^{+}$, is equivalent to the discrete Lagrangian dynamics (55).

• 2.

Given the mechanical Hamiltonian (42) and $L_{d}(z_{k},z_{k+1})=\frac{1}{2h}(z_{k+1}-z_{k})\,m\,(z_{k+1}-z_{k})-h\,U(z_{k+1})$, i.e. we pick $\kappa=0$ in (45a),(57) , then (63) provides a discretisation of (41).

• 3.

Under the hypotheses of Statement 2, when $\alpha=1/2$ (63) provides a discretisation of (44), i.e. the Hamiltonian dynamics with linear damping.

Proof.

Statement 1: from (63a), (63b) and $\tilde{\mathcal{L}}_{d}^{k}=L_{d}(x_{k},x_{k+1})+L_{d}(y_{k},y_{k+1})-h\,x_{k+1}^{\alpha}\rho y_{k+1}^{\alpha}$ it follows:

[TABLE]

Now setting the momentum matching condition $p_{x_{k}}^{-}=p_{x_{k}}^{+}$, $p_{y_{k}}^{-}=p_{y_{k}}^{+}$ it is straightforward to obtain (55).

Statement 2: under the hypotheses, we have that (41) read

[TABLE]

On the other hand, from (63) we get

[TABLE]

(Observe that, in the last equation, we have taken a $k$-step forward in (63c)). From this, splitting $x$ and $y$ sides and rearranging terms, we obtain

[TABLE]

which are a natural discretisation of (64).

Statement 3: we focus on the $x$-system in (65). Replacing the fractional equation into the third one, we obtain

[TABLE]

where, in $=^{{}_{1}}$ we have employed Lemma 4.2 and in $=^{{}_{2}}$ we have employed the discrete $x$-dynamics, i.e. $x_{k+1}=x_{k}+h\,m^{-1}p_{x_{k}}$. Naturally, the previous equations are a discretisation of (44). ∎

Furthermore, the discrete Legendre transformation in Definition 4.3 provides a initialisation step $p_{x_{0}}=I(x_{0},x_{1})$ (Step 2) for Algorithm 3. Namely, the $x$-part of (63a) reads

[TABLE]

which, for $k=0$, only involves $p_{x_{0}},x_{0}$ and $x_{1}$. For the particular $\tilde{\mathcal{L}}_{d}^{k}$ in the theorem above, the initial condition reads $p_{x_{0}}=-D_{1}L_{d}(x_{0},x_{1})$.

Remark 4.6**.**

The matrix $\begin{bmatrix}0&1\\ 1&0\end{bmatrix}$ in second term of the right hand side of (63a) and (63b) obeys to the necessity of decoupling $x$ and $y$ dynamics at the discrete level, which we achieved by restricting the variations and setting the critical conditions as (only) sufficient in the Lagrangian side, as shown in Proposition 4.1 and Remark 4.4. In other words, it can be considered as a discrete Hamiltonian consequence of the restricted Hamilton’s principle. $\diamond$

Remark 4.7**.**

It is interesting to note that the result remains the same for any $\tilde{\mathcal{L}}_{d}^{k}=L_{d}(x_{k},x_{k+1})+L_{d}(y_{k},y_{k+1})-\tilde{\kappa}\,h\,x_{k}^{\alpha}\,\,\rho\,\,y_{k}^{\alpha}-(1-\tilde{\kappa})\,h\,x_{k+1}^{\alpha}\,\,\rho\,\,y_{k+1}^{\alpha}$, with $\tilde{\kappa}\in[0,1]$, which is a way of rephrasing Remark 4.5. However, the presence of $x^{\alpha}_{k}$ turns the initial condition (66) meaningless from a physical point of view, which makes convenient setting $\tilde{\kappa}=0$. In that case, the pick of $L_{d}(x_{k},x_{k+1})$, which implies a particular choice of $\kappa$ in (45a) and (57), leads to different discretisations of (64) and (44). We remark that the chosen one ($\kappa=0$) preserves the semi-implicitness of variables $x,p_{x}$ of the symplectic-Euler methods [37] for (44); say: the final integrator is explicit in the variable $p_{x}$ and implicit in the variable $x$. $\diamond$

5. Numerical simulations

As a first test example, we employ the linearly damped harmonic oscillator with potential function $U(x)=\mathfrak{c}x^{2}/2$ and dynamical equation (exactly solvable):

[TABLE]

with $m=1$, $\mathfrak{c}=1$, $\rho=0.2$, $x(0)=1$ and $p_{x}(0)=\dot{x}(0)=0.5$ in the simulations.

In Figure 4, we show the outcome of Algorithm 3 with initial condition $p_{x_{0}}=-D_{1}L_{d}(x_{0},x_{1})$ according to (66), for several $\alpha$’s, where we choose $\kappa=1/2$ in (57) since it provides the midpoint rule for the potential and it is where the maximum local truncation order (namely 2) is achieved in usual low-order variational integrators [23]. We observe that the FVI approximates properly the solution of (67) when $\alpha=1/2$, which is natural since that is the case when (34a) $\rightarrow$ (67) (in other words, $D^{{}_{1/2}}_{-}D^{{}_{1/2}}_{-}=d/dt$). Moreover, according to Corollary 4.2, that is also the case when the FVI is equivalent to the Forced Variational Integrator (coming out of the discrete Lagrange-d’Alembert principle), Algorithm 2, when $f_{L_{d}}^{+}(x_{k},x_{k+1})=-\rho\,(x_{k+1}-x_{k}).$ This theoretical agreement is numerically tested (and shown) up to machine rounding error in Figure 5 (Lower-Left plot), for the different implementations of both Algorithms 3 and 2. We also show the comparison of the FVI to implicit and explicit Euler integrators for (44), choosing a smaller $h$ for the latter, which is necessary to obtain stable simulations for the explicit Euler scheme. In particular, we compare the $x$-trajectories in Figure 5 (Upper plots) and the energy in Figure 6 (Left and Middle plots), where naturally, we define the continuous and discrete energies by $E(t)=p(t)^{2}/2m+\mathfrak{c}\,x(t)^{2}/2$ and $E_{k}=p_{k}^{2}/2m+\mathfrak{c}\,x_{k}^{2}/2$ for $k=0,...,N$; respectively. While implicit and explicit Euler artificially gains respectively looses energy, the FVI respects the energy decay due to the dissipation much better and is very close to the exact solution.

We finally do a numerical convergence study by investigating the global error in both $x$ and $p_{x}$ variables (Lower-Right plot in Figure 5) as well as for the energy (Right plot in Figure 6). Here, the global error is defined as

[TABLE]

and equivalently for any other quantity. For all quantities the convergence is of order $O(h^{0.94})$, i.e. we obtain a convergence rate of approximately $1$.

As for $\alpha\neq 1/2$, we employ the following test example:

[TABLE]

and initial conditions $x(0)=\dot{x}(0)=0$. This equation corresponds to (34a) when $U(x)=x^{2}/2$, $m=\rho=1$, $\alpha=3/4$ and we add the inhomogeneous external force $F(t)$ in the right hand side (which can be easily carried out by equating the virtual work of such a force and the variation of the actions (27) and (49), in continuous and discrete scenarios, respectively, under restricted calculus of variations). Equation (68) has an exact solution333Particularly, it is a so-called inhomogeneous sequential fractional differential equation, i.e.

$\left(D^{n/q}_{-}+a_{1}D^{(n-1)/q}_{-}+\cdots+a_{n}D^{0}_{-}\right)\,x(t)=F(t),$

with $(n,q)=(4,2)$, $a_{1}=a_{4}=1$, $a_{2}=a_{3}=0$ [26]., but it is of difficult implementation. For that reason, we employ the benchmark numerical solution designed for MatLab in [32]. We display the performance of the FVI versus the benchmark solution (with a much smaller time step) (Left plot in Figure 7) and observe a global convergence of order $O(h^{0.97})$ (Right plot in Figure 7), i.e. a convergence rate of approximately $1$ as for the $\alpha=1/2$ case.

6. Conclusions

We have developed a restricted Hamilton’s principle providing the dynamics of Lagrangian systems subject to fractional damping (1) (continuous setting), (2) (discrete setting), as sufficient conditions for the extremals of the action. The discrete dynamics (2) is the result of the discretisation of the mentioned restricted Hamilton’s principle (instead of the discretisation of the equations (1) themselves), following the spirit of discrete mechanics and variational integrators [23]. As it is well-know, the variational principles and preservation properties (symplecticity, Noether’s theorem [2]) of the generated dynamics are closely related; we find a particular example in our approach, say the preservation of the presymplectic form $\Omega$ (35b), as proven in Proposition 3.3. Nevertheless, the dynamical importance of $\Omega$ has to be clarified in future work. This two-form is defined on the fractional state space, Definition 3.2, which is designed to accomodate the fractional derivatives. It is a vector bundle over the real space, particularly $\mathbb{R}^{d}\times\mathbb{R}^{d}$, which is necessary due to the unclear unique definition of fractional derivatives (3) on a general smooth manifold $Q$. From the geometrical perspective, an interesting challenge for future work is to carry out this generalisation.

The discretisation of the restricted fractional principle leads to the discrete equations (2), and to what we denote Fractional Variational Integrators FVIs, via the Algorithm 3. When $\alpha=1/2$, the theoretical local truncation order of FVIs is 1 [21], which is consistent with the observed global convergence order, i.e. $O(h^{0.94})$, in Figure 5, Lower-Right plot. However, this global convergence and the local truncation order 1 represents an improvement from what is expected from the order theorems in [23, 31], since, as proven in [25], $\Delta^{\alpha}_{-}x_{k}$ (45b) is only a consistent approximation of $D^{\alpha}_{-}x(t_{k})$ (3a), and thus one would expect a slower global convergence. This is an interesting phenomenom to explore in future works. Moreover, the performance of FVI is proven to be superior to other order 1 methods, such as both Euler’s. This is particularly evident with respect to the energy, as shown in Figure 6, accounting for another example of the superior performance of integrators with variational origin in this aspect. In the context of this work, this is furthermore explained thanks to the relationship between the FVI and the Forced Variational Integrators, Algorithm 2, obtained from the discrete Lagrange-d’Alembert principle, which is established in Corollary 4.2. Naturally, we also test the FVI for $\alpha\neq 1/2$.

We notice that the chosen discretisations limit the local truncation order of the FVIs and their global convergence. Thus, as a natural extension of this work, we intend to carry out the application of higher-order techniques, as introduced in [28, 30].

Appendix: Existence and uniqueness of solutions of fractional differential equations (32)

We study the existence and uniqueness of solutions of (32). It has been proven, Proposition 3.2, that $x$ and $y$ systems are equivalent in reversed directions of time. Thus, we focus on the $x$-system and set $d=1$ for simplicity. We are considering $L:T\mathbb{R}\rightarrow\mathbb{R}$ a $C^{2}$ function; furthermore we shall consider $\left(\partial^{2}L/\partial\dot{x}\partial\dot{x}\right)^{-1}$ smooth. All in all, (32a), can be expressed as

[TABLE]

with initial condition $x(t_{0})=x_{0},\,v(t_{0})=v_{0}$, $\beta=2\alpha$ and $f,\bar{\rho}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ given by

[TABLE]

with $\rho\in\mathbb{R}_{+}$, adding up for the vector field $F:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$

[TABLE]

In proving the existence and uniqueness of solutions of (69), we shall take a local approach; in particular we consider the set $\mathcal{B}=I_{t}\times I_{x}\times I_{v}\subset\mathbb{R}^{3}$, with $I_{t}=[t_{0}-\delta_{t},t_{0}+\delta_{t}],\,I_{x}=[x_{0}-\delta_{x},x_{0}+\delta_{x}],\,I_{v}=[v_{0}-\delta_{v},v_{0}+\delta_{v}]$, $I_{t}\subset[a,b]$ and $\delta_{t},\delta_{x},\delta_{v}\in\mathbb{R}_{+}$. We consider $\mathbb{R}^{d}$ as a Banach space equipped with the norm

[TABLE]

Given that, we define the cylinder $\mathcal{C}=I_{t}\times\mathcal{D}$, with $\mathcal{D}=\left\{(x,v)\in\mathbb{R}^{2}\,|\,||(x,v)-(x_{0},v_{0})||\leq\mathfrak{b}\right\}$, $\mathfrak{b}=||(\delta_{x},\delta_{v})||.$ We establish the following hypotheses:

• H1.

$||f(x,v)-f(\tilde{x},\tilde{v})||\leq K||(x,v)-(\tilde{x},\tilde{v})||$ in $\mathcal{D}$ for $K\in\mathbb{R}_{+}$.

• H2.

$\bar{\rho}$ is continuous in $\mathcal{D}$.

Note that these two hypotheses follow directly from the assumptions over $L$ and $\left(\partial^{2}L/\partial\dot{x}\partial\dot{x}\right)^{-1}$, say they are $C^{2}$ and smooth, respectively. Given this, our strategy is to prove that the vector field (70) satisfy the required conditions to apply both Peano and Picard-Lindelöf theorems [22], ensuring the existence and uniqueness of solutions of (69) in $\tilde{\mathcal{C}}=\tilde{I}_{t}\times\mathcal{D}$, with $\tilde{I}_{t}=[t_{0}-\tilde{\delta}_{t},t_{0}+\tilde{\delta}_{t}]$, $\tilde{\delta}_{t}=\,$min$\left\{\delta_{t}\,,\,\mathfrak{b}/M\right\}$, where $M\in\mathbb{R}_{+}$ is constructed in the proof below. Moreover, this solution $(x(t),v(t))$ can be set smooth by establishing the proper Banach space of functions when applying the aforementioned theorems.

Proposition 6.1**.**

Given the hypotheses H1, H2, the following is true: $F:\mathcal{D}\rightarrow\mathbb{R}^{2}$

1. (1)

is bounded. 2. (2)

is Lipschitz continuous.

Proof.

(1) Let $(x,v)\in\mathcal{D}$, then $||F(x,v)||=$max $\left\{|v|,|f(x,v)+\bar{\rho}(x,v)D^{\beta}_{-}x|\right\}$. On the one hand, $|v|\leq|v_{0}+\mathfrak{b}|<\infty$. On the other, by H1 $f$ is continuous in $\mathcal{D}$ and therefore $|f(x,v)|\leq M_{f}<\infty$. By H2, $\bar{\rho}$ is also continuous; thus $|\bar{\rho}(x,v)|\leq M_{\bar{\rho}}<\infty$. Now, we shall use the Caputo definition (3b) of the fractional derivative in (70), just by using the relationship (4) and setting $x_{0}=0.$ With that, in $I_{t}\times\mathcal{D}$ we have

[TABLE]

where, according to $\beta=2\alpha$, we consider $\alpha\in[0,1/2]$, since $\beta$ must be less of equal to 1444Observe that this limits the original range of $\alpha$, i.e.  $\alpha\in[0,1]$. Allowing that, the displayed proof of existence and uniqueness is not valid as long as we admit $\delta_{t}\rightarrow 0$, which does not mean that such existence and uniqueness of solutions may be proven in a different way.. From this, it follows that $||F(x,v)||\leq M<\infty,$ where $M=\left\{|v_{0}+\mathfrak{b}|\,\,\mbox{or}\,\,M_{f}+\frac{M_{\bar{\rho}}|v_{0}+\mathfrak{b}|}{\Gamma(1-\beta)}\delta_{t}^{1-\beta}\right\}$, depending on whether the max $\left\{|v|,|f(x,v)+\bar{\rho}(x,v)D^{\beta}_{-}x|\right\}$ is achieved in the first or second entry.

(2) We have that

[TABLE]

which is equal to the maximum of the absolute value of both entries. On the one hand

[TABLE]

by construction. On the other

[TABLE]

where we have used H1, H2 and (72). Thus, it follows that

[TABLE]

where $\tilde{K}=\left\{1\,\mbox{or}\,K+\frac{M_{\bar{\rho}}\delta_{t}^{1-\beta}}{\Gamma(1-\beta)}\right\}$, depending on whether the maximum on the right hand side of (71) is achieved on the first or second entry. In both cases, $\tilde{K}\in\mathbb{R}_{+}.$ ∎

Acknowledgments: This work has been funded by the EPSRC project: ‘’Fractional Variational Integration and Optimal Control”; ref: EP/P020402/1. FJ thanks Farhang Haddad Farshi for his assistance in the design of Figure 1.