1 Introduction

Fractional calculus is an extension of classical integration and differentiation theory to any real or complex order (Samko et al. 1993; Anatoly et al. 2006; Podlubny 1999; Oldham and Spanier 1974). A major feature enjoyed by fractional calculus is nonlocality which is used widely to model numerous phenomena in engineering, physics (including classical and quantum mechanics), biology, chemistry, etc. Baleanu (2008); Miller and Ross (1993); Silva et al. (2008); Agrawal and Baleanu (2007); Hilfer (2000b, 2000a).

The fractional variational principle was first proposed in Agrawal (2002). It is a tool to obtain a fractional equation of motion by defining Lagrangian actions depending on fractional derivatives. A discrete version of this approach was introduced in Bastos et al. (2011); Bourdin et al. (2013).

Fractional derivatives are successfully used to derive dissipative systems in the framework of variational principles (Riewe 1996; Cresson 2007). In this context, a new variational approach, the so-called restricted Hamilton’s principle, was developed in Jiménez and Ober-Blöbaum (2021). Contrary to the classical Hamilton’s principle, the restricted one gives only sufficient conditions for the extremals of the fractional variational problem, leading to equation that describes dynamics of fractionally damped systems, the so-called restricted fractional Euler–Lagrange equation:

$$\begin{aligned} \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{x}}\right) -\frac{\partial L}{\partial x}=-\mu \,D^{\alpha }_{-}x,\quad \mu >0,\quad \alpha \ge 0 \end{aligned}$$
(1)

where \(L(t,x,\dot{x})\) is a Lagrangian and \(D_-^{\alpha }\) is the fractional derivative operator. In particular, equation (1) with full order derivatives models the dynamics of forced systems (with a force linear in \({\dot{x}}\), such as linear dissipation). It can also describe the motion of a rigid plate immersed in a Newtonian fluid with fractional derivatives of order 3/2 (Torvik and Bagley 1984).

As a particular case, with \(\alpha =1\), there exist several approaches to derive equation (1) (Riewe 1996; Galley 2013; de Diego and de Almagro 2018; Cresson 2007, 2013; Jiménez and Ober-Blöbaum 2018, 2021).

From a numerical point of view, dealing with Lagrangian variational problems permits us to construct variational integrators (Marsden and West 2001; Hairer et al. 2006) via a discrete calculus of variations which are numerical schemes preserving the structural characteristics of the solution of the system. This approach is extended to a large class of systems including forced systems as well as nonholonomic ones (Monforte 2002; McLachlan and Perlmutter 2006; Marsden and West 2001).

Following (Jiménez and Ober-Blöbaum 2018), a discrete version of (1), also called Fractional Variational Integrator (FVI), is obtained from a discrete restricted Hamilton’s principle and fractional difference approximation. The (FVI) is mainly a combination of two different discretization schemes: one for the conservative part and the other for the fractional part. The authors use the midpoint rule for approximating a Lagrangian in the conservative part and the so-called Grünwald–Letnikov approach to approximate the \(\alpha \)-fractional derivatives in the fractional part. As stated in Podlubny (1999), this approximation is of first order. Then, the integrator (FVI) with \(\alpha =1\) has been proven to be of order one Jiménez and Ober-Blöbaum (2019).

Another type for approximating fractional integrals and derivatives is known as convolution quadrature (CQ). Such method was introduced by Lubich in Lubich (1986, 1988) and based on a specific quadrature rule for evaluating the convolution integral. The main difference between these methods and other numerical methods is that the weights of CQ are computed by the Laplace transform of the convolution kernel and multistep methods. As the left fractional integral is a particular convolution integral, the quadrature weights are obtained from the fractional order power of the rational polynomial of the generating functions of linear multistep methods (LMMs). In particular, the use of the classical backward difference formulas (BDFs), i.e., subclass of LMMs, is widely adopted for high accuracy (Lubich 1985, 1986).

In Diethelm and Ford (2002), the authors consider a particular type of fractional differential equations, called the Bagley–Torvik equation (Torvik and Bagley 1984) of the form \(Ax''(t)+B D_-^{3/2} x(t) +C x(t) = f(t)\). This equation can be recast as a system of four fractional differential equations of order 1/2, and then, Lubich’s correction term (Lubich 1986, 1988) has to be applied to obtain high-order numerical methods.

In the continuous setting, the integration by parts and semigroup properties are crucial in the restricted Hamilton’s principle for constructing (1). Thus, preserving such properties at the discrete level is then essential which can be fulfilled by CQ.

As we deal with higher-order approximation, it is also natural to apply higher-order techniques, also known as high-order variational integrators or Galerkin variational methods, for the conservative part. These techniques are applied to the action integral associated to a Lagrangian L in order to construct numerical schemes of arbitrarily high order. It is based on interpolating the trajectories and choosing a high-order quadrature for the approximation of the integral, see (Ober-Blöbaum and Saake 2015; Ober-Blöbaum 2016; Marsden and West 2001; Leok and Shingel 2012) and references therein.

Main contributions. In this paper, we follow the same strategy as in the continuous case in Jiménez and Ober-Blöbaum (2018). However, we use different discretization schemes, in particular, for the fractional part resulting in new fractional variational integration schemes. The main contributions of this paper are as follows:

  • We demonstrate how to use CQ based on fractional (BDF) methods for the approximation of fractional Lagrangians.

  • We combine those CQ-based approximations with Galerkin variational methods for the conservative part (including the midpoint rule) to derive high-order variational factional schemes for (1).

  • We prove the approximation order for a specific example with quadratic Lagrangian.

  • We verify numerically that our methods have second-order accuracy applied to the coupled damped oscillator and the scalar Bagley–Torvik equations including an 1/2-order derivative and an 3/2-order derivative.

Structure of the paper. In Sects. 2–3, we introduce and summarize the underlying mathematical theories our work is based on. In Sect. 2, we give a brief exposition of Lagrangian variational integrators that will be used throughout the work. Section 3 contains some necessary preliminaries of fractional calculus, including the notion of CQ as in §3.2 and, in particular, Lubich’s fractional linear multistep methods. We also discuss the main issue of using CQ with numerical experiments as illustration. After recalling a continuous restricted Hamilton’s principle in Sect. 4, we present our new contribution as a generalization of the one obtained by Jiménez and Ober-Blöbaum (2021) for deriving FVIs associated with (1) by means of the CQ in the framework of the discrete restricted Hamilton’s principle. We examine the accuracy, both, analytically for a special case and numerically for the coupled damped oscillator and the Bagley–Torvik equations. Section 5 treats the higher-order FVIs case by applying high-order techniques as presented in Ober-Blöbaum and Saake (2015); Ober-Blöbaum (2016) mixed with CQ.

2 Higher-Order Discrete Variational Mechanics

In this section, we recall the construction of variational integrators that will be used in this work. We refer to Marsden and West (2001); Leok and Shingel (2012); Ober-Blöbaum and Saake (2015) and references therein for more details.

2.1 Hamilton’s Principle and Euler–Lagrange Equations

Consider a mechanical system defined on the d-dimensional configuration manifold Q (later on we will particularize on \({\mathbb {R}}^d\), but in this section it can be considered as a general smooth manifold) with corresponding tangent bundle TQ. Let \(q(t)\in Q\) and \(\dot{q}(t)\in T_{q(t)}Q\), \(t\in [0,T]\subset {\mathbb {R}}\), \(0<T\), be the time-dependent configuration and velocity of the system. The action \({\mathscr {S}}:C^2([0,T],Q)\rightarrow {\mathbb {R}}\) of a mechanical system is defined as the time integral of the Lagrangian, i.e., 

$$\begin{aligned} {\mathscr {S}}(q)=\int _0^TL(q(t),\dot{q}(t))\,dt, \end{aligned}$$
(2)

where the \(C^2\) Lagrangian function \(L:TQ\rightarrow {\mathbb {R}}\) consists of kinetic minus potential energy. Hamilton’s principle seeks curves q, with fixed initial and final values q(0) and q(T), which are extremals of the action, i.e.,  satisfying

$$\begin{aligned} \delta \, {\mathscr {S}}(q)=0, \end{aligned}$$

for arbitrary variations \(\delta q\in T_qC^2([0,T],Q).\) A necessary and sufficient condition for the extremals is the so-called Euler–Lagrange equation

$$\begin{aligned} \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}}\right) -\frac{\partial L}{\partial q}=0, \end{aligned}$$
(3)

which is a second-order differential equation and describes the dynamics of conservative systems. See (Abraham and Marsden 1978) for more details.

2.2 Discrete Hamilton’s Principle and Discrete Euler–Lagrange Equations

The discretization of the objects described in the previous subsection is based on Marsden and West (2001); Moser and Veselov (1991). Let us consider a time grid \(t_k=\left\{ k\,h\,| k=0,\ldots , N\right\} \), where \(h\in {\mathbb {R}}_+\) is the time step and \(h\,N=T\). We replace the configuration q(t) by a discrete sequence \(q_d\equiv \left\{ q_k\right\} _{0:N}\in Q^{N+1}\) where \(q_k\) will be an approximation of \(q(t_k)\). The discrete Lagrangian \(L_d:Q\times Q\rightarrow {\mathbb {R}}\) will be an approximation of the action (2) in one time step \([t_k,t_{k+1}]\) based on two neighboring configurations \(q_k\) and \(q_{k+1}\), i.e., 

$$\begin{aligned} L_d(q_k,q_{k+1})\approx \int _{t_k}^{t_{k+1}}L(q(t),\dot{q}(t))\,dt. \end{aligned}$$
(4)

Furthermore, the discrete action sum \({\mathscr {S}}_d:Q^{N+1}\rightarrow {\mathbb {R}}\) is defined by

$$\begin{aligned} {\mathscr {S}}_d(q_d)=\sum _{k=0}^{N-1}L_d(q_k,q_{k+1}). \end{aligned}$$

The discrete Hamilton’s principle seeks extremals of the action \({\mathscr {S}}_d\) with fixed endpoints \(q_0,\,q_N\), i.e., 

$$\begin{aligned} \delta {\mathscr {S}}_d(q_d)=0, \end{aligned}$$

for arbitrary \(\delta q_k\in T_{q_k}Q\). A necessary and sufficient condition for the extremals is the discrete Euler–Lagrange equations

$$\begin{aligned} D_1L_d(q_k,q_{k+1})+D_2L_d(q_{k-1},q_{k})=0,\quad \quad k=1,\ldots , N-1, \end{aligned}$$
(5)

where \(D_i\) is the derivative with respect to the i-th argument. Given that the matrix \(D_{12}L_d(q_k,q_{k+1})\) is regular, equation (5) provides a discrete iteration scheme for (3) that determines \(q_{k+1}\) for given \(q_{k}\) and \(q_{k-1}\). This iteration scheme, which can be represented by the map \(F_{L_d}:Q\times Q\rightarrow Q\times Q\), \((q_{k-1},q_k)\mapsto (q_k,q_{k+1})\), is called variational integrator, and has interesting preservation properties, such as symplecticity and momentum preservation under the action of a symmetry (Marsden and West 2001; Moser and Veselov 1991).

2.3 Higher-Order Approximations of the Action

Considering only two neighboring configurations \(q_k\) and \(q_{k+1}\) in (4) limits the approximation order to \(O(h^2)\). With the aim of increasing this order, a well-known approach is to take into account inner nodes in between \([t_k, t_{k+1}]\) (Marsden and West 2001; Ober-Blöbaum and Saake 2015; Hall and Leok 2015). This higher-order approximation procedure consists of two steps: (1) the approximation of the space of trajectories and (2) the approximation of the integral of the Lagrangian by appropriate quadrature rules. (1) Trajectories space The space \({\mathscr {C}}([t_k,t_{k+1}],Q)=\left\{ q:[t_k,t_{k+1}]\rightarrow Q\,|\, q(t_k)=q_k,\,q(t_{k+1})=q_{k+1} \right\} ,\) will be approximated by \({\mathscr {C}}^s([t_k,t_{k+1}],Q)\subset {\mathscr {C}}([t_k,t_{k+1}],Q)\), where \({\mathscr {C}}^s([t_k,t_{k+1}],Q)\) denotes the space of polynomials of degree s. Given \(s+1\) control points \(0=d_0<d_1<\cdots<d_{s-1}<d_s=1\) and \(s+1\) configurations \(q_k=(q_k^0,q_k^1,\ldots ,q_k^{s-1},q_k^s)\), with \(q_k^0=q_k\) and \(q_k^{s}=q_{k+1}\), then \(q_d(t;k)\in {\mathscr {C}}^s([t_k,t_{k+1}],Q)\) can be defined by

$$\begin{aligned} q_d(t;k)=\sum _{\nu =0}^sq_k^{\nu }\,\ell _{\nu }\left( \frac{t}{h}\right) , \end{aligned}$$
(6)

where \(\ell _{\nu }(\tau )\) are Lagrange polynomials of degree s such that \(\ell _{\nu }(d_i)= \delta _{\nu i}\) (here \(\delta \) is the Kronecker symbol), and therefore \(q_d(h\,d_i;k)=q_k^i\) according to (6). Moreover, the time derivative of (6) is

$$\begin{aligned} \dot{q}_d(t;k)=\frac{1}{h}\sum _{\nu =0}^sq_k^{\nu }\,{\dot{\ell }}_{\nu }\left( \frac{t}{h}\right) . \end{aligned}$$

(2) Quadrature for the action integral For the approximation of the action integral (2), first we replace the curves q(t) and \(\dot{q}(t)\) by their polynomial counterparts \(q_d(t;k)\), \(\dot{q}_d(t;k)\) in the interval \([k\,h,\,(k+1)\,h]\), i.e., 

$$\begin{aligned} \int _{kh}^{(k+1)h}L(q_d(t;k),\dot{q}_d(t;k))\,dt,\quad k=0,\ldots , N-1. \end{aligned}$$

Next, in the same time interval, a quadrature rule \((b_i,c_i)_{i=1}^r\) is applied, with \(c_i\in [0,1]\). This defines the discrete Lagrangian:

$$\begin{aligned} L_d(q_k)\equiv L_d(q_k^0, \ldots ,q_k^s):=h\sum _{i=1}^rb_iL(q_d(c_i\,h;k),\dot{q}_d(c_i\,h;k)); \end{aligned}$$
(7)

it is important to remark that \(L_d\) depends on \(s+1\) variables. Naturally, the choice of quadrature should be adapted to the desired order of approximation with respect to the continuous action in order that now can be arbitrarily high.

The construction of higher-order variational integrators can be summarized in Fig. 1.

Fig. 1
figure 1

Polynomial interpolation principles. On each subinterval \([t_k,t_{k+1}]\), the trajectory interpolated by a polynomial passing through the points \(\{q_k^i\}_{i=0}^{s}\) associated with the control points \(\{d_\nu h\}\). The evaluations are made by the quadrature points for the \(\{c_ih\}\), indicated by cross points

Full size image

Now, following §2.2, the action sum can be defined from (7) by

$$\begin{aligned} {\mathscr {S}}_d(q_d)=\sum _{k=0}^{N-1}L_d(q_k), \end{aligned}$$

and the discrete Hamilton’s principle can be applied in order to obtain a necessary and sufficient condition for its extremals, i.e.,  the discrete Euler–Lagrange equations. These are:

$$\begin{aligned} \begin{aligned} D_{s+1}L_d(q_{k-1}^0,\ldots ,q_{k-1}^s)+D_1L_d(q_{k}^0,\ldots ,q_{k}^s)=0,\\ D_iL_d(q_{k}^0,\ldots ,q_{k}^s)=0,\quad \forall \,i=2,\ldots ,s; \end{aligned} \end{aligned}$$

for \(k=1,\ldots ,N-1,\) where the transition condition, namely \(q_{k-1}^s=q_k^0\), has to be taken into account. See (Marsden and West 2001; Ober-Blöbaum and Saake 2015) for further details.

Higher-order variational integrators (also denoted as Galerkin variational integrators) have been extensively studied in Hall and Leok (2015); Leok and Shingel (2012); Ober-Blöbaum and Saake (2015); Campos (2014) for conservative systems and in Campos et al. (2015) for optimal control problems. Obviously, the convergence order of higher-order variational integrators is limited by the order of the function space approximation and the order of the quadrature rule. In Hall and Leok (2015), a lower error bound is provided by the approximation order of the finite-dimensional function space (e.g., convergence order s is reached by using the space of polynomials of degree s). Numerical convergence studies in Ober-Blöbaum and Saake (2015) indicate that Galerkin variational integrators based on the Lobatto and Gauss quadrature rules are of order \(\min {(2s,u)}\), where s is the degree of the polynomial and u the order of the quadrature rule. A general proof of this superconvergence result is provided in Ober-Blöbaum and Vermeeren (2021) using backward error analysis in the context of the calculus of variations. For particular classes that are equivalent to so-called (modified) symplectic Runge–Kutta methods, a proof of a superconvergence has been made in Ober-Blöbaum (2016).

3 Fractional Integrals, Fractional Derivatives and Their Approximations

3.1 Fractional Integrals and Fractional Derivatives

Let us start by giving a brief overview concerning the fractional operators. We refer to Samko et al. (1993); Oldham and Spanier (1974); Anatoly et al. (2006) for more details.

3.1.1 Fractional Integrals

The Riemann–Liouville \(\alpha \)-fractional integrals, Re \(\alpha >0\), for \(f:[0,T]\rightarrow {\mathbb {R}}\)Footnote 1 are defined by

$$\begin{aligned} J^{\alpha }_{-}f(t)&=\frac{1}{\Gamma (\alpha )}\int _0^t(t-\tau )^{\alpha -1}f(\tau )\,d\tau ,\,\,\,t\in (0,T],\end{aligned}$$
(8a)
$$\begin{aligned} J^{\alpha }_{+}f(t)&=\frac{1}{\Gamma (\alpha )}\int _t^T(\tau -t)^{\alpha -1}f(\tau )\,d\tau ,\,\,\,t\in [0,T), \end{aligned}$$
(8b)

where \(\Gamma \) is the Euler Gamma function and we set \(J^{0}_{-}f=J^{0}_{+}f=f\). The fractional integrals satisfy the so-called semigroup property (Samko et al. 1993, Theorem 2.5, p.46), i.e., 

$$\begin{aligned} J_{\sigma }^{\alpha }\,J_{\sigma }^{\beta }f(t)=J_{\sigma }^{\alpha +\beta }f(t),\quad \sigma \in \left\{ -,+\right\} . \end{aligned}$$
(9)

Let \(m\in {\mathbb {N}}\). If the function f is continuously differentiable in [0, T], then it can be continued analytically to Re \(\alpha <0\) via

$$\begin{aligned} J^{\alpha }_{-} f(t)=\frac{d^m}{dt^m}J^{m+\alpha }_{-}f(t), \quad J^{\alpha }_{+} f(t)=\frac{d^m}{dt^m}J^{m+\alpha }_{+}f(t) \quad \text{ for } \text{ Re } \ \alpha >-m. \end{aligned}$$
(10)

In this case, the fractional integral is called fractional derivative, and can be denoted by \(D^{-\alpha }\) (note that the real part of \(\alpha \) is now negative).

In particular, from (10), and restricting ourselves to Re \(\alpha \in [0,1]\) (which will be the range of interest in this work), we can establish the following definition of Riemann–Liouville fractional derivatives, which is usually found in the literature:

3.1.2 Riemann–Liouville Fractional Derivatives

Let Re \(\alpha \in [0,1]\). The left and right Riemann–Liouville fractional derivatives are, respectively, defined by

$$\begin{aligned} \begin{aligned} D^{\alpha }_{-}f(t)&=\,\,\,\,\frac{d}{dt}\,J^{1-\alpha }_{-}f(t)=\quad \frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}\int _0^t(t-\tau )^{-\alpha }f(\tau )\,d\tau ,\,\,\,t\in (0,T], \\ D^{\alpha }_{+}f(t)&=-\frac{d}{dt}\,J^{1-\alpha }_{+}f(t)=-\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}\int _t^T(\tau -t)^{-\alpha }f(\tau )\,d\tau ,\,\,\,t\in [0,T), \end{aligned} \end{aligned}$$
(11)

provided that \(f\in AC([0,T],{\mathbb {R}})\) which is a very simple sufficient condition for the existence. It is easy to see that \(D^{0}_{-}f=D^{0}_{+}f=f\), whereas it can be proven that

$$\begin{aligned} D^{1}_{-}f=-D^{1}_{+}f=df/dt. \end{aligned}$$
(12)

The last relationships follow easily from the definitions (11) (first equality) and \(J^{0}_{-}f=J^{0}_{+}f=f\), but the latter are not trivial from the definitions (8).Footnote 2

Other relevant properties of fractional derivatives are

$$\begin{aligned} \int _0^Tf(t)D^{\alpha }_{\lambda }g(t)\, dt&=\int _0^T g(t)\big [D^{\alpha }_{-\lambda }f(t)\big ] \,dt,\quad \sigma \in \left\{ -,+\right\} , \end{aligned}$$
(13a)
$$\begin{aligned} D_{\lambda }^{\alpha }D_{\lambda }^{\beta }&=D_{\lambda }^{\alpha +\beta },\quad 0\le \alpha ,\beta \le 1/2. \end{aligned}$$
(13b)

Property (13a) is called “asymmetric integration by parts” of the fractional derivatives, whereas (13b) is called again the “semigroup property.”

Remark 3.1

The left fractional derivative in the sense of Caputo, also based on left Riemann–Liouville fractional integral (Anatoly et al. (2006), p. 73) is given by,

$$\begin{aligned} ^{*}D_-^\alpha f(t)=J_-^{1-\alpha }\dot{f}(t). \end{aligned}$$

In the case \(0< \alpha < 1\) with some regularity on f, we have the following relation

$$\begin{aligned} ^{*}D_-^\alpha f(t)=D_-^\alpha f(t) - \frac{t^{-\alpha }}{\Gamma (1-\alpha )}f(0). \end{aligned}$$

and it is important to note that, contrary to the Caputo derivative, the Riemann–Liouville derivative of a function f when \(f(0)\ne 0\) could lead to singularity at \(t=0\). That is why the Caputo derivatives are more useful in applications, otherwise we should impose \(f(0)=0\).

3.2 Discrete Convolution

For the time being, let us focus on the retarded fractional integral in (8), i.e.,  \(J^{\alpha }_{-}f(t)\). It is easy to see that it is a particular convolution integral

$$\begin{aligned} J^{\alpha }_{-}f(t)=(\kappa ^{(\alpha )}* f)(t):=\int _0^t\kappa ^{(\alpha )}(t-\tau )f(\tau )d\tau \end{aligned}$$
(14)

where the convolution kernel is given by

$$\begin{aligned} \kappa ^{(\alpha )}(t)=\frac{t^{\alpha -1}}{\Gamma (\alpha )}. \end{aligned}$$
(15)

Here, \((\alpha )\) must not be understood as a power in the left hand side, it is just a superscript denoting the dependence of the kernel on the parameter \(\alpha \).Footnote 3 The Laplace transform of this convolution kernel is given by

$$\begin{aligned} K^{(\alpha )}(s):={\mathscr {L}}(\kappa ^{(\alpha )})(s)=\int _0^{\infty }\frac{t^{\alpha -1}}{\Gamma (\alpha )}e^{-st}dt=s^{-\alpha }. \end{aligned}$$
(18)

The theory of discrete convolutions is developed in Lubich (1986, 1988, 1994) by Ch. Lubich (indeed, the topic of Lubich (1986) is the discretization of fractional integrals).

As in §2.2, let us consider a time grid \(t_k=\left\{ k\,h\,| k=0,\ldots , N\right\} \), where \(h\in {\mathbb {R}}_+\) is the time step and \(h\,N=T\). Moreover, consider a discrete series \(\left\{ f_{k}\right\} _{0:N}\in ({\mathbb {R}}^d)^{N+1}\), where \(f_k\) shall be an approximation of \(f(t_k)\).

Now, define the discrete convolution in the following as an approximation of the continuous convolution \((\kappa * f)(t_k)\) in (16), namely

$$\begin{aligned} (\kappa * f)(t_k)\approx (\omega *f)(t_k):=\sum _{n=0}^{\infty }\omega _nf_{k-n}=^{1}\sum _{n=0}^{k}\omega _nf_{k-n}, \end{aligned}$$
(16)

observe that the series is truncated after \(=^{1}\) since \(f_k\) is not defined for \(k<0\)Footnote 4, with a general convolution kernel \(\kappa \). The CQ weights \(\omega _n\) are defined as the coefficients of the generating power series

$$\begin{aligned} K\left( \frac{\gamma (z)}{h}\right) :=\sum _{n=0}^{\infty }\omega _nz^n,\quad \quad |z|\quad \text{ small }. \end{aligned}$$
(19)

Here, K(s) is the Laplace transform of the kernel \(\kappa \) and the so-called characteristic function \(\gamma (z)=\sum _{n=0}^{\infty }\gamma _nz^n\) is the quotient of the generating polynomials \((\rho ,\sigma )\) of a linear multistep method (LMM)

$$\begin{aligned} \rho _0y_k+\rho _1 y_{k-1}+\ldots +\rho _ny_{k-n}=h(\sigma _0f_k+\sigma _1 f_{k-1}+\ldots +\sigma _nf_{k-n}) \end{aligned}$$

for the differential equation \(y^{\prime }=f(y)\), i.e., 

$$\begin{aligned} \gamma (z)=\frac{\rho (z)}{\sigma (z)}=\frac{\rho _0+\rho _1 z+\ldots +\rho _nz^n}{\sigma _0+\sigma _1 z+\ldots +\sigma _nz^n}, \end{aligned}$$
(20)

where we assume that \(\rho _0/\sigma _0>0\), so that (19) is well defined at least for sufficiently small h. Note that if we define \(z_{-}\) to be the discrete backward operator, i.e., 

$$\begin{aligned} z_{-}\,f_k=f_{k-1}, \end{aligned}$$

the LMM can be defined as \(\rho (z_{-})\,y_k=h\,\sigma (z_{-})\,f_k\), where \(\left\{ y_k\right\} _{0:N}\in ({\mathbb {R}}^d)^{N+1}\) is also a discrete series.

More importantly for the purposes of this article, the discrete convolution approximating the fractional integral (8a),(14) can be redefined, according to (19), by

$$\begin{aligned} {\mathscr {J}}_{-}^{\alpha }f_k:=K^{(\alpha )}\left( \frac{\gamma (z_{-})}{h}\right) \,f_k=\sum _{n=0}^{k}\omega _n^{(\alpha )}f_{k-n}, \end{aligned}$$
(21)

where \(K^{(\alpha )}(\gamma (z)/h)\), with \(K^{(\alpha )}\) given in (17), has to be understood as an operator acting on \(\left\{ f_k\right\} _{0:N}\). On the other hand, considering the discrete forward operator

$$\begin{aligned} z_{+}\,f_k=f_{k+1}, \end{aligned}$$

the discrete convolution approximating the fractional integral (8b) will be

$$\begin{aligned} {\mathscr {J}}_{+}^{\alpha }f_k:=K^{(\alpha )}\left( \frac{\gamma (z_{+})}{h}\right) \,f_k=\sum _{n=0}^{N-k}\omega _n^{(\alpha )}f_{k+n}, \end{aligned}$$
(22)

where again the series (19) gets truncated because \(\left\{ f_k\right\} _{0:N}\) is undefined for \(k>N\).Footnote 5 It is interesting to note that, following these definitions, the convolution weights \(\omega ^{(\alpha )}_n\) are the same for \({\mathscr {J}}^{\alpha }_{-}\) and \({\mathscr {J}}^{\alpha }_{+}\).

Now, we take into consideration the convergence order of \({\mathscr {J}}_{\lambda }^{\alpha }\) with respect to \(J^{\alpha }_{\lambda }\) following (Lubich 1986); as we will see, this order is sensitive to the convergence order of the multistep method \((\rho , \sigma )\). In particular, for the generating function \(\gamma =\rho /\sigma \), which gives rise to \(\omega _n^{(\alpha )}\) according to (21), we say that the corresponding multistep method is convergent of order p if and only if

$$\begin{aligned} {\widetilde{\gamma }}_n\quad \text{ are } \text{ bounded },\,\,\text{(stability) } \end{aligned}$$
(23a)
$$\begin{aligned} h\,{\widetilde{\gamma }} (\,e^{-h}\,)=1+O(h^p),\quad \text{ as }\,\,h\rightarrow 0,\,\text{(order } \text{ p } \text{ consistency) }, \end{aligned}$$
(23b)

where \({\widetilde{\gamma }}_n\) are coefficients of the power series \(\gamma ^{-1}(z):=1/\gamma (z)\). Now, we introduce the notion of convergence of the quadrature \(\omega _n^{(\alpha )}\).

Definition 3.1

A CQ determined by the coefficients \(\omega _n^{(\alpha )}\) is convergent of order p (to \(J^{\alpha }_{\lambda }\)) if

$$\begin{aligned} J^{\alpha }_{\lambda }t^{\beta -1}-{\mathscr {J}}^{\alpha }_{\lambda }t^{\beta -1}=O(h^{\beta })+O(h^p), \end{aligned}$$
(24)

for all \(\beta \in {\mathbb {C}}\), \(\beta \ne 0,-1,-2,\ldots .\)

Theorem 3.1

(Theorem 2.6 in Lubich (1986)) Let \((\rho ,\sigma )\) denote an implicit linear multistep method which is convergent of order p (23) and assume that the zeros of \(\sigma \) have absolute value less than 1. Then, \({\mathscr {J}}^{\alpha }_{\lambda }\) (21), (22) are convergent of order p (Definition 3.1) to \(J^{\alpha }_{\lambda }\) (8).

In Lubich (1986) is also established, under the condition of p-convergence in Definition 3.1, that for functions \(f(t)=t^{\beta -1}g(t)\), g(t) smooth, there always exists a starting quadrature \(\varpi _{k,n}\)Footnote 6 which is defined by

$$\begin{aligned} \widetilde{{\mathscr {J}}}_{-}^{\alpha }f_k:=\sum _{n=0}^{k}\omega _{n}^{(\alpha )}f_{k-n}+h^{\alpha }\sum _{n=0}^s\varpi _{k,n}f_n, \end{aligned}$$
(25)

with s fixed, such that

$$\begin{aligned} J^{\alpha }_{\lambda }f-\widetilde{{\mathscr {J}}}_{\lambda }^{\alpha }f=O(h^p), \end{aligned}$$
(26)

uniformly for \(t\in [0,T]\). In other words, to achieve the order p convergence of the underlying LMM, the extra term in (26) should be introduced in order to eliminate low order terms in the error bound in (24).

It is important to remark that due to the presence of the extra \(\varpi \)-initial terms in (25) the convolution structure is violated, and therefore the properties proved in Lemma 4.1 are no longer true for \(\widetilde{{\mathscr {J}}}_{\lambda }^{\alpha }f\). In Lubich’s own words: “The convolution structure is only violated by the few correction terms of the starting quadrature which will be necessary for high order schemes.” Indeed, this is also relevant in terms of the convergence order, since from Definition 3.1, Theorem 3.1 and (26) we conclude that for a function \(f(t)=t^{\beta -1}g(t)\), with g(t) smooth, a multistep method \((\rho ,\sigma )\) is p-convergent would generate the following convergence bound:

$$\begin{aligned} J^{\alpha }_{\lambda }f-{\mathscr {J}}_{\lambda }^{\alpha }f=O(h^{\beta })+O(h^p). \end{aligned}$$
(27)

Thus, we expect the saturation of the convergence order at \(\min (\beta ,p)\).

Remark 3.2

In a similar way, the approximation for the Riemann–Liouville derivatives by CQ can be formally done by using \(-\alpha \) rather than \(\alpha \) in (21) and (22). Moreover, for a non-integer \(\alpha \), the standard relation between the Riemann–Liouville and Caputo derivatives allows one to use the CQ for approximating Caputo derivatives (Diethelm 2019)

$$\begin{aligned} ^{*}{\mathscr {J}}_{-}^{-\alpha }f_k:={\mathscr {J}}_{-}^{-\alpha }f_k-\sum _{n=0}^{\lceil \alpha \rceil -1}\frac{t_k^{n-\alpha }}{\Gamma (n+1-\alpha )}f^{(n)}(0). \end{aligned}$$

A disadvantage of using the Caputo definitions is that the associated CQ may destroy the integration by parts and the semigroup properties which we need in the next section. In particular, \(^{*}{\mathscr {J}}_{-}^{-\alpha }f_k\) and \({\mathscr {J}}_{-}^{-\alpha }f_k\) coincide if and only if \(f^{(n)}(0)=0\), \(n=0,1,\cdots ,\lceil \alpha \rceil -1\).

We illustrate the saturation issue (27) for the function \(t\mapsto t^{\beta -1}\sin (t)\) with \(\beta = 1,3,5\) using the Caputo derivative of order \(\alpha =1/2\) which is equivalent to the Riemann–Liouville derivative in this case. In Fig. 2, the error is defined by means of the maximum norm

$$\begin{aligned} e_h= \max _{0\le k\le N}\,\big |J^{\alpha }_{-}f(t_k)-{\mathscr {J}}_{-}^{\alpha }f_k\big |, \end{aligned}$$
(28)

versus h is plotted on Log–Log scale. Of our interested, the characteristic function \(\gamma (z_{-})\) of the classical backward differentiation formula (BDF) up to order 6, that is,

$$\begin{aligned} \gamma (z)=\sum _{k=1}^p \frac{1}{k}(1-z)^k:=\gamma _p(z). \end{aligned}$$
(29)
Fig. 2
figure 2

Log–Log plot the error Err(h) versus h corresponds to the Caputo fractional derivative of order 1/2 (\(\alpha =-1/2\) in (28)). As expected, the convergence starts to saturate at \(p=1\) (left), \(p=3\) (middle) and \(p=5\) (right)

Full size image

4 Restricted Variational Principle for the Dynamics of Lagrangian Systems Subject to Fractional Damping

In Jiménez and Ober-Blöbaum (2018, 2021) a restricted variational principle (both continuous and discrete) in order to obtain the dynamics of a Lagrangian system subject to fractional damping is delivered. As in other previous approaches treating dissipative systems in a variational fashion (Agrawal 2002; Bateman 1931; Cresson and Inizan 2012; Galley 2013; Riewe 1996), it is based on the doubling of variables and the introduction of their fractional derivatives (which in this work will be considered as fractional integrals with negative \(\alpha \) index, as explained in §3) in the state space of the relevant Lagrangians.

4.1 Continuous Setting

Let us consider the \(AC^2\) curves \(x,y:[0,T]\rightarrow {\mathbb {R}}^d\) and a \(C^2\)-Lagrangian function \(L:{\textsf{T}}{\mathbb {R}}^d\rightarrow {\mathbb {R}}\) (we use the natural identification \({\textsf{T}}{\mathbb {R}}^d\cong {\mathbb {R}}^d\times {\mathbb {R}}^d\)). Define a \(C^2\)-Lagrangian functionFootnote 7

$$\begin{aligned}&L^{(\alpha ,\beta )}:{\mathbb {R}}^d\times {\mathbb {R}}^d\times {\mathbb {R}}^d\times {\mathbb {R}}^d\times {\mathbb {R}}^d\times {\mathbb {R}}^d \rightarrow {\mathbb {R}}\nonumber \\&(x,\,y,\,\dot{x},\,\dot{y},\, J^{-\alpha }_{-}x,\, J^{-\beta }_{+}y) \mapsto L^{(\alpha ,\beta )}(x,y,\dot{x},\dot{y}, J^{-\alpha }_{-}x,J^{-\beta }_{+}y)\nonumber \\&L^{(\alpha ,\beta )}(x,\,y,\,\dot{x},\,\dot{y},\, J^{-\alpha }_{-}x,\, J^{-\beta }_{+}y)=L(x,\dot{x})+L(y,\dot{y})-\mu \,J^{-\alpha }_{-}x\,J^{-\beta }_{+}y, \end{aligned}$$
(30)

where \(\alpha ,\beta \in [0,1/2]\) and \(\mu \in {\mathbb {R}}_+\). Given this particular Lagrangian, we define the relevant action \({\mathscr {S}}: AC^2([0,T])\times AC^2([0,T])\rightarrow {\mathbb {R}}\) by \({\mathscr {S}}(x,y)={\mathscr {S}}^{^\text {cons}}(x,y)+{\mathscr {S}}^{^\text {frac}}(x,y)\) with

$$\begin{aligned} \begin{aligned} {\mathscr {S}}^{^\text {cons}}(x,y)&=\int _0^T(L(x(t),\dot{x}(t))+L(y(t),\dot{y}(t))\,)\, dt\\ {\mathscr {S}}^{^\text {frac}}(x,y)&=-\mu \,\int _0^TJ^{-\alpha }_{-}x(t)\,J^{-\beta }_{+}y(t)\,dt, \end{aligned} \end{aligned}$$
(31)

where cons goes after “conservative” and frac after “fractional.” Moreover, let us define restricted varied curves by means of

$$\begin{aligned} x_{\epsilon }(t)=x(t)+\epsilon \, \delta x(t),\,\,\,y_{\epsilon }(t)=y(t)+\epsilon \, \delta x(t), \end{aligned}$$
(32)

with \(\epsilon \in {\mathbb {R}}_{+}\) and an \(AC^2\) \(\,\,\delta x:[0,T]\rightarrow {\mathbb {R}}^d\) such that \(\delta x(0)=\delta x(T)=0\) (observe that the variations are equal for both curves, which is the base of the restriction). With these elements, and assuming fixed endpoints \(x(0),\,x(T),\,y(0),\,y(T)\), we can establish the restricted variational principle:

Theorem 4.1

The equations

$$\begin{aligned} \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{x}}\right) -\frac{\partial L}{\partial x}=-\mu \,J^{-(\alpha +\beta )}_{-}x,\end{aligned}$$
(33a)
$$\begin{aligned} \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{y}}\right) -\frac{\partial L}{\partial y}=-\mu \,J^{-(\alpha +\beta )}_{+}y, \end{aligned}$$
(33b)

are sufficient conditions for the extremals of \({\mathscr {S}}(x,y)\) (31) under restricted variations (32).

It is important to remark that in the proof of this theorem (see (Jiménez and Ober-Blöbaum 2021)), the use of the asymmetric integration by parts (13a) and the semigroup property (13b) of the fractional derivatives is crucial. In addition, it is also proven in (Jiménez and Ober-Blöbaum (2021), Proposition 3.2) that under even parity of L in the velocity variable, then (33b) reduces to (33a) in reversed time, i.e.,  \(y(t)=x(T-t)\). Finally, it is easy to see that the dynamics (33), say the Lagrangian dynamics subject to fractional damping, reduces to the usual linear damped dynamics when \(\alpha =\beta =1/2\), according to (12).

Remark 4.1

The statements in the traditional variational principles or the unrestricted fractional variational principles (Jiménez and Ober-Blöbaum (2021), Remark 3.2) are made with “sufficient and necessary conditions”; however, in the fractional variational principle under restricted variations the “necessary” part of the implication is eliminated, see (Jiménez and Ober-Blöbaum (2021), Theorem 3.1).

In the following, we give several intermediary lemmas which are obtained by generalizing the results presented in Jiménez and Ober-Blöbaum (2021).

4.2 Discrete Setting Based on CQ

In Jiménez and Ober-Blöbaum (2018, 2021), the author applied the discretization procedure described in §2.2 and the Grünwald–Letnikov approximation for the fractional derivative to derive the so-called fractional variational integrators. In the following, we propose a generalization of this process using CQ. For that, let us consider two discrete series \(x_d=\left\{ x_k\right\} _{0:N}\) and \(y_d=\left\{ y_k\right\} _{0:N}\), as well as two particular discretizations of (11), i.e., 

$$\begin{aligned} {\mathscr {J}}_{-}^{-\alpha }x_k=\sum _{n=0}^k\omega _n^{(-\alpha )}x_{k-n}, \qquad {\mathscr {J}}_{+}^{-\beta }y_k=\sum _{n=0}^{N-k}\omega _n^{(-\beta )}y_{k+n}, \end{aligned}$$
(34)

where the weights \(\omega _n^{(-\alpha )}\) are the coefficients of the generating power series of \(K^{(-\alpha )}(\gamma (z)/h)\) with \(K^{(\alpha )}\) defined in (17), namely

$$\begin{aligned} K^{(-\alpha )}\left( \frac{\gamma (z)}{h}\right) =\left( \frac{\gamma (z)}{h}\right) ^\alpha =\sum _{k=0}^\infty \omega _n^{(-\alpha )} z^n. \end{aligned}$$

For concreteness, we choose \(\gamma _p(z)\) as in (29), the generating function of the classical BDF methods. The Grünwald weights used in Jiménez and Ober-Blöbaum (2021) coincide with a particular choice of \(\gamma _p(z)\), i.e.,  \(\gamma _1(z)=1-z\) and the notations \({\mathscr {J}}_{-}^{-\alpha }x_k\) and \({\mathscr {J}}_{+}^{-\beta }y_k\) have been used as \(\Delta ^{\alpha }_{-}x_k\) and \(\Delta ^{\beta }_{+}y_k\), respectively.

Now we are in situation to show the proof of some properties relevant for future results. The following lemma is a generalization of the one given in Jiménez and Ober-Blöbaum (2021) by means of CQ.

Lemma 4.1

Consider two discrete series \(\left\{ f_k\right\} _{\tiny 0:N}, \left\{ g_k\right\} _{\tiny 0:N}\). Then, the following properties hold true:

  1. 1.

    The semigroup property of the discrete convolution: \({\mathscr {J}}^{\alpha }_{\lambda }\,{\mathscr {J}}^{\beta }_{\lambda }\,f_k={\mathscr {J}}^{\alpha +\beta }_{\lambda }\,f_k.\)

  2. 2.

    The asymmetric integration by parts:

    $$\begin{aligned} \sum _{k=0}^{N}g_k\,({\mathscr {J}}_{-}^{\alpha }f_k)=\sum _{k=0}^N({\mathscr {J}}_{+}^{\alpha }g_k)\,f_k. \end{aligned}$$
    (35)

Proof

(1) From the equalities (21) and (22) and as \(K^{(\alpha )}(s)=s^{-\alpha }\) (17), it follows that:

$$\begin{aligned} {\mathscr {J}}^{\alpha }_{\lambda }\,{\mathscr {J}}^{\beta }_{\lambda }\,f_k&=K^{(\alpha )}\left( \frac{\gamma (z_{\lambda })}{h}\right) K^{(\beta )}\left( \frac{\gamma (z_{\lambda })}{h}\right) \,f_k\\&=\left( \frac{\gamma (z_{\lambda })}{h}\right) ^{-\alpha }\left( \frac{\gamma (z_{\lambda })}{h}\right) ^{-\beta }\,f_k\\&=\left( \frac{\gamma (z_{\lambda })}{h}\right) ^{-\alpha -\beta } =K^{\alpha +\beta }\left( \frac{\gamma (z_{\lambda })}{h}\right) ={\mathscr {J}}^{\alpha +\beta }_{\lambda }\,f_k. \end{aligned}$$

(2) See (Bourdin et al. 2013) for more details:

$$\begin{aligned} \begin{aligned} \sum _{k=0}^Ng_k\,({\mathscr {J}}_{-}^{\alpha }f_k)&=^{_1}\sum _{k=0}^N\sum _{n=0}^k\omega _n^{(\alpha )}\,g_k\,f_{k-n}=^{_2}\sum _{n=0}^N\sum _{k=n}^N\omega _n^{(\alpha )}\,g_k\,f_{k-n}\\&=^{_3}\sum _{n=0}^N\sum _{k=0}^{N-n}\omega _n^{(\alpha )}\,g_{k+n}\,f_{k}=^{_4}\sum _{k=0}^N\sum _{n=0}^{N-k}\omega _n^{(\alpha )}\,g_{k+n}\,f_{k}=^{_6}\sum _{k=0}^{N}({\mathscr {J}}_{+}^{\alpha }g_k)\,f_k. \end{aligned} \end{aligned}$$

In \(=^{_1}\), the definition (21) (second equality) is used. To prove \(=^{_2}\) is enough to notice that, for fixed \(j=0, \ldots ,N,\) the elements \(a_{i}:=\omega _i^{(\alpha )}\,g_j\,f_{j-i}\), \(i=0, \ldots ,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, \ldots ,N\) and \(i=j, \ldots ,N\), account for the transposed matrix; therefore, their total sums are equal. In \(=^{_3}\), the sum index is rearranged. In \(=^{_4}\), equivalent arguments to \(=^{_2}\) can be used. Finally, in \(=^{_6}\) the definition (22) (second equality) is used; this concludes the proof. \(\square \)

The discrete action, counterpart of (31), is then

$$\begin{aligned} \begin{aligned}&{\mathscr {S}}_d(x_d,y_d)={\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)+{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d), \\&{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)=\sum _{k=0}^{N-1}(L_d(x_k,x_{k+1})+L_d(y_k,y_{k+1})),\\&{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d)=-\mu \,h\,\,\sum _{k=0}^{N}{\mathscr {J}}_{-}^{-\alpha }x_k \,{\mathscr {J}}_{+}^{-\beta }y_k. \end{aligned} \end{aligned}$$
(36)

Taking the discrete equivalent of the restricted varied curves, i.e., 

$$\begin{aligned} x_d^{\epsilon }=x_d+\epsilon \, \left\{ \delta x_k\right\} _{0:N},\,\,\,y_d^{\epsilon }=y_d+\epsilon \, \left\{ \delta x_k\right\} _{0:N}, \end{aligned}$$
(37)

such that \(\delta x_0=\delta x_N=0\), we can establish the discrete counterpart of Theorem 4.1:

Theorem 4.2

The equations

$$\begin{aligned} D_1L_d(x_k,x_{k+1})&+D_2L_d(x_{k-1},x_{k})=\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k},\end{aligned}$$
(38a)
$$\begin{aligned} D_1L_d(y_k,y_{k+1})&+D_2L_d(y_{k-1},y_{k})=\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_{k}, \end{aligned}$$
(38b)

for \(k=1,\ldots ,N-1\), are sufficient conditions for the extremals of \({\mathscr {S}}_d(x_d,y_d)\) (36) under restricted variations (37).

See (Jiménez and Ober-Blöbaum 2021) for the proof. As in the continuous side, the semigroup property and asymmetric integration by parts of the discrete fractional derivatives/integrals, properties (1) and (2) in Lemma 4.1, respectively, are crucial in the proof of this theorem. It is also true that (38b) reduces to (38a) in discrete reversed time, i.e.,  \(y_k=x_{N-k}\) ( (Jiménez and Ober-Blöbaum 2021), Proposition 4.2). Naturally, the equations (38) provide a discrete iteration scheme for the fractional damped dynamics (33), with \(\gamma _1(z)=1-z\), delivering a \(p=1\) convergent integrator; see (Jiménez and Ober-Blöbaum 2021) for further details.

4.2.1 Approximation Order

In a traditional variational integrator, the accuracy is assessed by variational error analysis (Marsden and West 2001); in other words, it is sufficient to construct a discrete Lagrangian that approximates the exact discrete Lagrangian with order p to obtain a variational integrator of order p. The question is then can we determine the accuracy of a fractional variational integrator from fractional variational error analysis? In general, this analysis concept is not directly applicable to a fractional Lagrangian since it is defined on a different space (including fractional derivatives) compared to the classical Lagrangian considered in Marsden and West (2001). The extension to the fractional case is not immediately clear and still open. In Diethelm and Ford (2002), a different approach is applied for a very particular class of fractional differential equations, the Bagley–Torvik equation with Caputo derivative of order 3/2 by reformulating the equation to a system of fractional differential equations of order 1/2. In that paper, the modified CQ method (25) is used which makes it possible to obtain numerical schemes of higher order. However, modified CQ methods destroy the structures, i.e., the integration by parts and the semigroup properties which we need in our approach. For those reasons, we focus on an error analysis for a particular case of fractional Lagrangian: Let us consider a Lagrangian associated to linearly damped systems

$$\begin{aligned} L(x,\dot{x})=\frac{1}{2}\dot{x}^\top \dot{x}-U(x).\end{aligned}$$
(39)

The Lagrangian defined in (30) reads

$$\begin{aligned} L^{(\alpha ,\beta )}(x,y,\dot{x},\dot{y}, J_-^{-\alpha } x,J_+^{-\beta } y)=L(x,\dot{x})+L(y,\dot{y})-\mu J_-^{-\alpha } xJ_+^{-\beta } y. \end{aligned}$$

The restricted fractional Euler–Lagrange equation associated to (33a) can be written as fractional Hamiltonian system

$$\begin{aligned} \begin{aligned} \dot{x}&= p\\ \dot{p}&= -\nabla U(x)-\mu D_-^{(2\alpha )}x. \end{aligned} \end{aligned}$$
(40)

With \(\alpha =1/2\) in (40), we recover the classical forced system

$$\begin{aligned} \begin{aligned} \dot{x}&= p\\ \dot{p}&= -\nabla U(x)-\mu \dot{x}. \end{aligned} \end{aligned}$$
(41)

To compare the first order Hamiltonian system (41) with the FVI (33a), we first have to introduce the discrete momenta based on the discrete Legendre transform. Now introduce the discrete Legendre transform for \(p_k^-\) and \(p_k^+\) as

$$\begin{aligned} p_{k}^-&=-D_1L_d(x_k,x_{k+1}) +\kappa \mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k}\\ p_{k+1}^{+}&=D_2L_d(x_k,x_{k+1}) - (1-\kappa )\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k+1}, \end{aligned}$$

where \(\kappa \in [0,1]\) is a real parameter. This allows to rewrite equivalently the fractional DEL equations as matching condition for the momenta \(p_k:=p_k^+ = p_k^-\) and hence a unique momentum \(p_k\) is defined at every grid point along a solution of the FVI.

With our particular Lagrangian (39), the FVI can then equivalently be written as

$$\begin{aligned} p_{k}&=\frac{x_{k+1}-x_k}{h} + \frac{h}{2} \nabla U\left( x_{k+1/2}\right) +\kappa \mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k} \end{aligned}$$
(42a)
$$\begin{aligned} p_{k+1}&=\frac{x_{k+1}-x_k}{h} - \frac{h}{2} \nabla U\left( x_{k+1/2}\right) -(1-\kappa ) \mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k+1} \end{aligned}$$
(42b)

and solved for \((x_{k+1},p_{k+1})\). In particular, for \(\alpha =\beta =1/2\), the CQ associated with BDF2 becomes

$$\begin{aligned} {\mathscr {J}}^{-(\alpha +\beta )}_{-}x_{k}=\frac{1}{h}\left( \frac{3}{2}x_k-2x_{k-1}+\frac{1}{2}x_{k-2}\right) . \end{aligned}$$

Theorem 4.3

Let \(\alpha = \beta = \kappa = 1/2\) with CQ method associated to BDF2. The local truncation order of the FVI for the system (42), with respect to the continuous dynamics (41), is two.

Proof for BDF2

With the notation \(x_k:= x(t_k )\), \(p_k:=p(t_k)=\dot{x}(t_k)\), we perform the Taylor expressions for \(x(t_{k+1})\) and \(p(t_{k+1})\) in terms of (41), namely:

$$\begin{aligned} x({t_{k+1}})&= x_k+ h p_k - \frac{h^2}{2} ( \rho p_k + \nabla U(x_k)) + \frac{h^3}{6}\rho \nabla U(x_k) \nonumber \\&\quad + \frac{h^3}{6} p_k\left( \rho ^2-\Delta U(x_k)\right) +O(h^4)\end{aligned}$$
(43)
$$\begin{aligned} p({t_{k+1}})&=p_k- h\left( \rho p_k + \nabla U(x_k)\right) + \frac{h^2}{2}\rho \nabla U(x_k) \nonumber \\&\quad + \frac{h^2}{2} p_k\left( \rho ^2-\Delta U(x_k)\right) +O(h^3) \end{aligned}$$
(44)

On the other hand, solving (42a)+(42b)\(= p_{k+1}+p_{k}\) w.r.t \(x_{k+1}\) leads to

$$\begin{aligned} x_{k+1}&= \frac{1}{4 - 3 h \rho + 3 \kappa h \rho }\Big [2 h (p_k + p_{k+1}) + (4 - 4 h \rho +\kappa h\rho ) x_k \\&\quad - \kappa h \rho x_{k-2} + h \rho x_{k-1} + 3 \kappa h \rho x_{k-1}\Big ]\\&= x_k+ h p_k - \left( \kappa \rho p_k + \frac{1}{2}\nabla U(x_k)\right) h^2 - \frac{h^3}{4} (1 - 2 \kappa ) \rho \nabla U(x_k) \\&\quad +\frac{\rho ^2}{8}h^3 ((1 - 5 \kappa + 6 \kappa ^2)p_k -\frac{h^3}{4} \Delta U(x_k)) +O(h^4) \end{aligned}$$

Comparing the previous equation with (43) and \(\kappa =\frac{1}{2}\) yields \(\Vert x(t_{k+1})-x_{k+1}\Vert =O(h^3)\). Furthermore, solving (42b)-(42a)\(= p_{k+1}-p_{k}\) w.r.t \(p_{k+1}\) leads to

$$\begin{aligned} p_{k+1}&= p_k + \frac{\rho }{2} (4 - 7 \kappa ) x_k - \frac{\kappa \rho }{2} x_{k-2} - (1 - 5 a)\frac{\rho }{2}x_{k-1} - 3 (1 - \kappa )\frac{\rho }{2} x_{k+1}\\&\qquad - h \nabla U(x_{k+1/2})\\&= p_k -(\rho p_k +\nabla U(x_k))h + (1 - \kappa ) \rho (\rho p_k + \nabla U(x_k))h^2 \\&\qquad -\frac{1}{2} p_k \Delta U(x_k)h^2 +O(h^3) \end{aligned}$$

From this last expression and (43) with \(\kappa =\frac{1}{2}\), we obtain that \(\Vert p(t_{k+1})-p_{k+1}\Vert =O(h^3)\). This completes the proof. \(\square \)

In the following section, we will apply the order two midpoint variational integrator for the conservative part (Marsden and West 2001) and BDFCQ for the fractional one. Then, we will denote by FVI-BDFCQ the scheme (38a) and write it simply FVI when no confusion can arise. All the experiments are carried out in Julia Version 1.9.3.

4.3 Numerical Experiment

Let us consider the Lagrangian associated to the undamped free coupled oscillators system of the form Pribõtkin and Tomasiello (2023)

The dynamics of fractionally forced oscillators can be deduced from (33a) with \(\alpha =\beta \) as follows:

$$\begin{aligned} m_1\ddot{x} +\alpha _1x_1 + \beta _1x_1^3 + \varepsilon (x_1-x_2)=- \mu _1\,J_-^{-(2\alpha )} x_1, \end{aligned}$$
(45a)
$$\begin{aligned} m_2\ddot{x} +\alpha _2x_2 + \beta _2x_2^3 + \varepsilon (x_2-x_1) =- \mu _2\,J_-^{-(2\alpha )}x_2. \end{aligned}$$
(45b)

In this experiment, we fix \(N=2\) with the parameters \(m_1=m_2=1,\alpha _1=0.01,\alpha _2=0.02,\beta _1=\beta _2=0,\varepsilon =0.5,\mu _1=\mu _2=3\sqrt{\alpha _1/2}, x_{01}=0.8, x_{02}=-0.4,\dot{x}_{01}=\dot{x}_{02}=0\).

Consider the problem on [0, 16] and choose \(\alpha =1/2\) so that the fractional operator \(D_-^{(2\alpha )}\) coincides with the usual operator d/dt. In this test run, we take the initial values \(x(0)=0,\ \dot{x}(0)=1.2\) and \(\mu =0.2\). Firstly, we plot the numerical solution obtained by the FVI-BDF1CQ, explicit and implicit Euler integrators on the interval [0, 16] with the stepsize \(h = 0.125\) in Fig. 3 (left). The corresponding results of energy dissipation and absolute errors are presented in Fig. 3 (right) and Fig. 4 (left), respectively. Secondly, we compute global errors as the maximum norm between the numerical solution and the exact solution, i.e.,

$$\begin{aligned} \max |x(t_k) - x_k|,\quad \forall k. \end{aligned}$$

The results are presented in Fig. 4 (right) as the global errors (in logarithmic scale) against stepsizes on [0, 16] with \(h =16/2^i,\ i=4,\ldots ,11\).

Fig. 3
figure 3

Damped harmonic oscillator (45) (\(\alpha =1/2\)). Left: Exact solution vs FVI-BDF1CQ method for \(h=0.125\). Right: Energy behavior for \(h=0.125\)

Full size image

The main property of a dissipative system is that the energy is always dissipated with time, and as we can see in Fig. 3 (right), FVI-BDF1CQ can preserve the dissipation structure of the damped harmonic oscillator which confirm that FVI-BDF1CQ gives good numerical behavior.

Fig. 4
figure 4

Damped harmonic oscillator (45). Left: absolute errors for \(h=0.125\). Right: Log–Log plot of the global error presented on \(t\in [0,16]\) for \(h =16/2^i,\ i=4,\ldots ,11\)

Full size image

We can confirm from Fig. 4 (right) that the order of FVI-BDF1CQ is one and this result has been discussed in Jiménez and Ober-Blöbaum (2021). However, we observe that the second-order convergence both FVI-BDF2CQ and FVI-BDF3CQ which is natural since the midpoint integrator being used is of second order.

We also consider another example. Let us choose a Lagrangian of the forced harmonic oscillator problem defined by \(L(t,x,\dot{x})=\dot{x}^2/2-x^2/2+x\, f(t)\) with a nonvanishing function f. So that equation (33a), again for \(\alpha =\beta \), reads

$$\begin{aligned} \ddot{x}+\mu \, D_-^{(2\alpha )} x+x=f(t),\quad t\in [0,1]. \end{aligned}$$
(46)

For a non-integer \(2\alpha \), this problem is known as the Bagley–Torvik equation, which can be used to describe, for example, the dynamics of a rigid plate immersed in a Newtonian fluid when \(\alpha =3/4\) (see (Podlubny 1999; Torvik and Bagley 1984)). Due to mathematical complexity, the analytic solutions of such equation are very few and are restricted to the one-dimensional case. In particular, with the initial conditions \(x(0)=\dot{x}(0)=0\) and \(\mu =1\), the Bagley–Torvik equation is exactly solvable (see (Jena and Chakraverty 2019; Ford and Connolly 2009)) by considering,

$$\begin{aligned} f(t)&= t^3 + 6 \, t + \frac{3.2}{\Gamma (1/2)}\, t^2 \sqrt{t},\quad \alpha =\frac{1}{4} , \end{aligned}$$
(47a)
$$\begin{aligned} f(t)&= \frac{15}{4} \sqrt{t} + \frac{15}{8} \sqrt{\pi } \, t + t^2 \sqrt{t},\quad \alpha = \frac{3}{4}. \end{aligned}$$
(47b)

where the analytic solutions are given, respectively, by \(x(t)=t^3\) (resp. \(=t^{\frac{5}{2}}\)). We numerically solve the Bagley–Torvik problem (46) using FVI on [0, 1] for \(\alpha = \frac{1}{4},\ \frac{3}{4}\). The global errors (on logarithmic scale) are presented in Figs. 5 for \(h = 1/2^i,\ i=1,\ldots ,8\).

Fig. 5
figure 5

Bagley–Torvik equation (46). Log–Log plot of the global errors on \(t\in [0,1]\) for \(h = 1/2^i,\ i=1,\ldots ,8\). Left: case (47a). Right: case (47b)

Full size image

Again, it can be observed from Fig. 5 that FVI-BDF1CQ leading to a convergence of order one. A convergence of order 2 is obtained for FVI-BDF2CQ and FVI-BDF3CQ (left plot), whereas a convergence of order cannot reach two for FVI-BDF2CQ and FVI-BDF3CQ (right plot).

We summarize the convergence order of (38a) for equations (45) and (46) in Table 1 where we consider the midpoint integrator for the conservative part.

Table 1 Convergence order of (38a) for equations (45) and (46)
Full size table

5 Higher-Order Fractional Variational Integrators Based on Convolution Quadrature

Now, we establish a particular discretization of the action (31), where we choose a higher-order approximation with quadrature rule \((b_i,c_i)_{i=1}^r\) (§2.3) for the conservative part \({\mathscr {S}}^{^\text {cons}}\) and CQ (21), (22) for the fractional integrals involved in \({\mathscr {S}}^{^\text {frac}}\) instead of the order 1 (34). For that, we take into account two discrete series \(x_d=\left\{ x_k\right\} _{0:N}\in ({\mathbb {R}}^d)^{N+1}\), \(y_d=\left\{ y_k\right\} _{0:N}\in ({\mathbb {R}}^d)^{N+1}\) and \(s+1\) inner nodes \(\left\{ x_k^{\nu }\right\} ^{0:s}\in ({\mathbb {R}}^d)^{s+1}\) in each interval \([k,k+1]\) such that \(x_k^s=x_{k+1}^0\) (equiv. for y). Namely

$$\begin{aligned} \begin{aligned}&{\mathscr {S}}_d(x_d,y_d)={\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)+{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d), \\&{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)=\sum _{k=0}^{N-1}(L_d(x_k)+L_d(y_k)),\\&\quad {\mathscr {S}}^{^\text {frac}}_d(x_d,y_d)=-\mu \,h\,\,\sum _{k=0}^{N}{\mathscr {J}}^{-\alpha }_{-}x_k\,{\mathscr {J}}^{-\beta }_{+}y_k,\\&L_d(x_k)=h\sum _{i=1}^rb_iL(x_d(c_i\,h;k),\,\dot{x}_d(c_i\,h;k)),\\&\,\,\,\,\,\, L_d(y_k)=h\sum _{i=1}^rb_iL(y_d(c_i\,h;k),\dot{y}_d(c_i\,h;k)), \end{aligned} \end{aligned}$$
(48)

where the definition (6) applies for \(x_d(t;k)\) and \(y_d(t;k)\) just by \(Q={\mathbb {R}}^d\). Now, considering restricted varied curves

$$\begin{aligned} x_d^{\epsilon }=\left\{ x_k^{\nu }\right\} _{0:N-1}^{0:s}+\epsilon \, \left\{ \delta x_k^{\nu }\right\} _{0:N-1}^{0:s},\,\,\,y_d^{\epsilon }=\left\{ y_k^{\nu }\right\} _{0:N-1}^{0:s}+\epsilon \, \left\{ \delta x_k^{\nu }\right\} _{0:N-1}^{0:s}, \end{aligned}$$
(49)

such that \(\delta x_0=\delta x_0^0=0\) and \(\delta x_N=\delta x_{N-1}^s=0\), we establish the following result (it is important to recall that, from now on, we shall consider the variation operator as \(\delta \equiv d/d\epsilon |_{\epsilon =0}\), applied over the “varied” quantities). Now, we set the following useful lemma.

Lemma 5.1

According to the definitions (21), (22) and considering varied curves (49), we have

$$\begin{aligned} \delta \, {\mathscr {J}}^{-\alpha }_{\lambda }x_k = {\mathscr {J}}^{-\alpha }_{\lambda }\,\delta x_k. \end{aligned}$$

Equivalently for y.

Proof

We pick \(\sigma =-\), the proof for \(+\) is equivalent. It is important to note that in the CQ (21) the inner nodes are not involved, and consequently from (49) we only take into consideration the main nodes, i.e.,  \(x_d^{\epsilon }=\left\{ x_k\right\} _{0:N}+\epsilon \, \left\{ \delta x_k\right\} _{0:N}\).

$$\begin{aligned} \begin{aligned} \delta \, {\mathscr {J}}^{-\alpha }_{-}x_k&=\frac{d}{d\epsilon }\Big |_{\epsilon =0}\, {\mathscr {J}}^{-\alpha }_{-}x_d^{\epsilon }=\frac{d}{d\epsilon }\Big |_{\epsilon =0} \sum _{n=0}^{k}\omega _n^{(-\alpha )}(x_{k-n}+\epsilon \,\delta x_{k-n})\\&=\sum _{n=0}^{k}\omega _n^{(-\alpha )}\delta x_{k-n}={\mathscr {J}}^{-\alpha }_{-}\,\delta x_k. \end{aligned} \end{aligned}$$

\(\square \)

Theorem 5.1

The equations

$$\begin{aligned}&D_{s+1}L_d(x_{k-1}^0, \ldots ,x_{k-1}^s)+ D_1L_d(x_k^0, \ldots ,x_k^s)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_k^0=0,\nonumber \\&\qquad k=1,\ldots , N-1,\end{aligned}$$
(50a)
$$\begin{aligned}&D_iL_d(x_k^0, \ldots ,x_k^s)=0,\,\hspace{1.3cm}\,\,\quad k=0,\ldots , N-1,\,\quad i=2,\ldots ,s,\end{aligned}$$
(50b)
$$\begin{aligned}&D_{s+1}L_d(y_{k-1}^0, \ldots ,y_{k-1}^s)+ D_1L_d(y_k^0, \ldots ,y_k^s)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_k^0=0,\nonumber \\ \qquad k=1,\ldots , N-1,\end{aligned}$$
(50c)
$$\begin{aligned}&D_iL_d(y_k^0, \ldots ,y_k^s)=0,\,\hspace{1.3cm}\,\,\quad k=0,\ldots , N-1,\qquad i=2,\ldots ,s, \end{aligned}$$
(50d)

are sufficient conditions for the extremals of (48) under restricted variations (49).

Proof

From (48), we have that

$$\begin{aligned} \delta {\mathscr {S}}_d(x_d,y_d)=\delta \,{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)+\delta \,{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d). \end{aligned}$$

Let start to simplify \(\delta \,{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)\)

$$\begin{aligned} \begin{aligned} \delta \,{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)=\sum _{k=0}^{N-1}\left( \frac{\partial L_d(x_k)}{\partial x_k^{\nu }}+\frac{L_d(y_k)}{\partial y_k^{\nu }}\right) \,\delta x_k^{\nu }, \end{aligned} \end{aligned}$$

where the summation over \(\nu \) is understood and we have employed the restricted variations (49). Concerning the term \(\delta \,{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d)\), we use the notation

$$\begin{aligned} \frac{\partial L_d(x_k)}{\partial x_k^{\nu }}\delta x_k^{\nu }=D_iL_d(x_k)\,\delta x_k^{\nu _i}, \end{aligned}$$
(51)

where on the left hand side \(\nu =0, \ldots ,s\) and on the right hand side \(\nu _1=0,\,\nu _2=1,\, \ldots ,\nu _{s+1}=s\) (in other words \(D_i=\partial /\partial x_k^{\nu _i}\)); this way, it is highlighted that \(L_d\) is a function of \(s+1\) variables (equiv. for y). Thus, we have

$$\begin{aligned} \begin{aligned} \delta \,{\mathscr {S}}^{^\text {cons}}_d(x_d,y_d)=&\sum _{k=0}^{N-1}\sum _{i=1}^{s+1}\left( D_iL_d(x_k)+D_iL_d(y_k)\right) \,\delta x_k^{\nu _i}.\\ \delta \,{\mathscr {S}}^{^\text {frac}}_d(x_d,y_d)=^1&-\mu \,h\,\,\sum _{k=0}^{N}{\mathscr {J}}^{-\alpha }_{-}x_k\,{\mathscr {J}}^{-\beta }_{+}\delta \, x_k -\mu \,h\,\,\sum _{k=0}^{N}{\mathscr {J}}^{-\alpha }_{-}\delta \, x_k\,{\mathscr {J}}^{-\beta }_{+}y_k\\ =^2&-\mu \,h\,\,\sum _{k=0}^{N}{\mathscr {J}}^{-\beta }_{-}{\mathscr {J}}^{-\alpha }_{-}x_k\,\delta \, x_k -\mu \,h\,\,\sum _{k=0}^{N}\delta \, x_k\,{\mathscr {J}}^{-\alpha }_{+}{\mathscr {J}}^{-\beta }_{+}y_k\\ =^3&-\mu \,h\,\,\sum _{k=1}^{N-1}\left( {\mathscr {J}}^{-(\beta +\alpha )}_{-}x_k^0+{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_k^0\right) \,\delta \, x_k^0. \end{aligned} \end{aligned}$$

In \(=^1\), we have employed the Leibnitz rule of the derivative and Lemma 5.1. In \(=^2\), we have employed the asymmetric integration by parts, i.e.,  property (2) in Lemma 4.1. Finally, in \(=^3\) we have rearranged terms, employed the semigroup property (1) in Lemma 4.1, taken into account that \(x_k=x_k^0\), \(y_k=y_k^0\) and \(\delta \, x_k=\delta \,x_k^0\) in terms of the inner nodes and taken also into account that \(\delta x_0=\delta x_N=0\), such that the terms \(k=0\) and \(k=N\) vanish. Putting everything together, we have

$$\begin{aligned} \delta \, {\mathscr {S}}_d(x_d,y_d)&=\sum _{k=0}^{N-1}\sum _{i=1}^{s+1}\left( D_iL_d(x_k)+D_iL_d(y_k)\right) \,\delta x_k^{\nu _i}\\&\quad -\mu \,h\,\,\sum _{k=1}^{N-1}\left( {\mathscr {J}}^{-(\alpha +\beta )}_{-}x_k^0+{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_k^0\right) \,\delta \, x_k^0\\&=\sum _{k=0}^{N-1}\sum _{i=2}^{s+1} D_iL_d(x_k)\,\delta x_k^{\nu _i}+\sum _{k=1}^{N-1}\left( D_1L_d(x_k)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_k^0\right) \delta \, x_k^0\\&\quad +\sum _{k=0}^{N-1}\sum _{i=2}^{s+1} D_iL_d(y_k)\,\delta x_k^{\nu _i}+\sum _{k=1}^{N-1}\left( D_1L_d(y_k)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_k^0\right) \delta \, x_k^0\\ \end{aligned}$$

where it is taken into account that \(\delta x_0=\delta x_0^0=0\) and that \(D_{s+1}L_d(x_{k-1})=D_1L_d(x_k).\) Now, given that \(\delta x_k^{\nu _i}\) are arbitrary for \(k=0,\ldots ,N-1\), \(i=1,\ldots ,s+1\) (except \(\delta x_0\)), we see from the last equality that

$$\begin{aligned} \begin{aligned}&D_iL_d(x_k)=0,\,\hspace{5cm}\,\,\quad k=0,\ldots , N-1,\,\quad i=2,\ldots ,s,\\&D_{s+1}L_d(x_{k-1})+ D_1L_d(x_k)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_k=0,\,\,\,\,\,\, k=1,\ldots , N-1,\\&D_iL_d(y_k)=0,\,\hspace{5cm}\,\,\quad k=0,\ldots , N-1,\,\quad i=2,\ldots ,s,\\&D_{s+1}L_d(y_{k-1})+ D_1L_d(y_k)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{+}y_k=0,\,\,\,\,\,\,\, k=1,\ldots , N-1, \end{aligned} \end{aligned}$$

is a sufficient condition for \(\delta \, {\mathscr {S}}_d(x_d,y_d)=0\) and the claim holds. \(\square \)

Remark 5.1

The partial derivatives of \(L_d(x_k)\) (equiv. y) are completely determined by the quadrature rule \((b_i,c_i)_{i=1}^r\) and the Lagrangian function \(L(q,\dot{q})\), according to (48) and (6). Namely

$$\begin{aligned} \begin{aligned} \frac{\partial L_d(x_k)}{\partial x_k^{\nu }}=&h\sum _{i=1}^rb_i\Bigg (\frac{\partial L}{\partial q}(x_d(c_i\,h;k),\,\dot{x}_d(c_i\,h;k))\frac{\partial x_d}{\partial x_k^{\nu }}\\&\qquad +\frac{\partial L}{\partial \dot{q}}(x_d(c_i\,h;k),\,\dot{x}_d(c_i\,h;k))\frac{\partial \dot{x}_d}{\partial x_k^{\nu }}\Bigg )\\ =&h\sum _{i=1}^rb_i\Bigg (\frac{\partial L}{\partial q}(x_d(c_i\,h;k),\,\dot{x}_d(c_i\,h;k))\,\ell _{\nu }(c_ih)\\&\qquad +\frac{\partial L}{\partial \dot{q}}(x_d(c_i\,h;k),\,\dot{x}_d(c_i\,h;k))\,\frac{1}{h}\dot{\ell }_{\nu }(c_ih)\Bigg ). \end{aligned} \end{aligned}$$

Naturally, equations (50a),(50b) can be employed as a discrete iteration scheme for the dynamics (33a), the same way (50c),(50d) can be used for (33b). We shall focus on the x-part, since y is equivalent.

The equations

$$\begin{aligned} \begin{aligned} p_{x_0}:=&-D_1L_d(x_0^0, \ldots ,x_0^s),\\ 0=&\,\,\,\,\,\,\,D_iL_d(x_k^0, \ldots ,x_k^s),\quad \forall \,i=2,\ldots ,s,\,\hspace{4cm}\,\,\,\quad k=0,\ldots , N-1,\\ 0=&\,\,\,\,\,\,\,D_{s+1}L_d(x_{k-1}^0, \ldots ,x_{k-1}^s)+ D_1L_d(x_k^0, \ldots ,x_k^s)-\mu \,h\,{\mathscr {J}}^{-(\alpha +\beta )}_{-}x_k^0,\,\,\,\,\,\,\,\,\,\,\,\, k=1,\ldots , N-1, \end{aligned} \end{aligned}$$

(where we include the initial momentum \(p_{x_0}:=-D_1L_d(x_0^0, \ldots ,x_0^s)\) as a definitionFootnote 8) conform a discrete iteration scheme

$$\begin{aligned} \begin{aligned} x_0^0&\,\,\mapsto \,\, (x_0^0, \ldots ,x_0^s=x_1^0),\\ (x_{k-1}^0, \ldots ,x_{k-1}^s=x_k^0)&\,\,\mapsto \,\,(x_{k}^0, \ldots ,x_{k}^s=x_{k+1}^0),\quad k=1, \ldots ,N-1, \end{aligned} \end{aligned}$$

that can be represented as an algorithm:

Algorithm 1
figure a

Higher-order fractional variational integrator (with CQ)

Full size image

It is important to remark that at each step a nonlinear system of s algebraic equations is solved in order to obtain the s unknowns \((x_k^1, \ldots ,x_k^s=x_{k+1}^0)\), even in the initialization step.

5.1 Numerical Experiment

We will employ a variational integrators of order 4 for the conservative part based on two points Gauss quadrature and a polynomial degree 2 (see (Ober-Blöbaum and Saake 2015)) and BDFCQ for the fractional one. To simplify notation, we continue to write FVI-BDFCQ (or only FVI when it is convenient) for equations 50a and 50b and the numerical solution can be computed using Algorithm 1.

Let us consider the previous examples as in §4.3. As expected, we notice in Figs. 6 and 7 that FVI-BDF1CQ and FVI-BDF2CQ are of first and second order, respectively. From the numerical point of view, one would expect a higher accuracy (order 3 using BDF3CQ mixed with third-order variational integrator) which we do not get. One possible reason might be mixing of integrators (VIs and BDFCQ). In particular, BDFCQ depends only the main nodes but not on the inner nodes which are considered in the conservative part. Another possible reason might be saturation effects coming from CQ as described in §3.2.

Table 2 Convergence order of (38a) for equations (45) and (46)
Full size table
Fig. 6
figure 6

Damped harmonic oscillator (45). Left: absolute errors for \(h=0.125\). Right: Log–Log plot of the global error presented on \(t\in [0,16]\) for \(h =16/2^i,\ i=4,\ldots ,11\)

Full size image
Fig. 7
figure 7

Bagley–Torvik equation (46). Log–Log plot of the global error on \(t\in [0,1]\) for \(h = 1/2^i,\ i=1,\ldots ,8\). Left: case (47a). Right: case (47b)

Full size image

As we have seen in §3.2, the main issue of using BDFCQ for certain class of solution functions is that one cannot achieve a high accuracy, see the saturation effects in Figs. 2. However, a correction term should be added as in (25) to recover the accuracy order as the one of the underlying BDF methods. We apply BDF3CQ with a correction term in the discrete Euler–Lagrange equation in Algorithm 1 for equation (46) when \(\alpha =3/4\) and as we observe in Fig. 8, the third accuracy order is almost achieved. However, this phenomenon does not work with the previous studied cases, which seems related to the accumulation of errors.

Fig. 8
figure 8

Bagley–Torvik equation (46), case (47b), i.e.,  \(\alpha =3/4\). Log–Log plot of the global error on \(t\in [0,1]\) for \(h = 1/2^i,\ i=1,\ldots ,8\)

Full size image

We summarize the convergence order of (38a) for equations (45) and (46) in Table 2 where we consider an integrator of order 4 for the conservative part.

6 Conclusions

Several methods have been proposed to deal with nonconservative mechanical systems, e.g., a restricted Hamilton’s principle (Jiménez and Ober-Blöbaum 2018, 2021), an approach introducing fractional derivatives. The main motivation of this approach is to derive the dynamics of fractionally damped systems (1) as direct as it is in the case of conservative systems, i.e., using a purely variational way.

We have developed a discrete version of the restricted Hamilton’s principle using two different procedures. One for conservative part (Marsden and West 2001; Ober-Blöbaum and Saake 2015; Hall and Leok 2015) and the other for the fractional part using the CQ method, which is particularly suitable (Lubich 1986, 1988) for approximating fractional derivatives.

This work centers around increasing the accuracy the numerical scheme associated to (1). Here, we have focused on implementing the FVIs and test numerically their accuracy using two mechanical systems, the coupled damped oscillator and the Bagley–Torvik problems. We notice that for FVI based on BDFCQ, it can only achieve the second-order accuracy even for a higher-order FVI (see Figs. 6 and 7) which is due to the fact that saturation effects are also a part of the problem. Another problem of using BDFCQ is that the inner nodes used in a higher-order approximation for the conservative action (31) are not taken into account in BDFCQ for the fractional one. Thus, a way to handle this is to apply the high-order Runge–Kutta CQ (RKCQ) (Lubich and Ostermann 1993; Banjai and Lubich 2011) for the fractional part, which will be a future work.