An adaptation of the Gear scheme for fractional derivatives

Abstract The Gear scheme is a three-level step algorithm, backward in time and second-order accurate for the approximation of classical time derivatives. In this contribution, the formal power of this scheme is proposed to approximate fractional derivative operators in the context of finite difference methods. Some numerical examples are presented and analysed in order to show the effectiveness of the present Gear scheme at the power α (Gα-scheme) when compared to the classical Grunwald–Letnikov approximation. In particular, for a fractional damped oscillator problem, the combined Gα-Newmark scheme is shown to be second-order accurate.


Introduction
The notion of fractional calculus appears in diverse fields of science and engineering. If the notion first appears during the 17th century, a precise definition has been proposed by Riemann-Liouville and more recently by Caputo [1].
The importance of fractional calculus for modeling viscoelastic materials has been recognized by the mechanical scientific community since the article of Bagley and Torvik [2]. The numerical approximation of such systems has been intensively studied since the work of Padovan [3]. For applications to complex mechanical systems, one refers to [4]. On the other side, the numerical community is interested in the approximation of fractional derivatives. One refers to the pioneering theoretical work of Lubich [5] and the state of the art proposed by Diethelm et al. [6]. Most applications use the discrete convolution formula proposed by Grü nwald-Letnikov [7,8]. Another direction could be autonomous systems in the context of diffusive representations [9][10][11].
In this work, we develop a numerical method based on the Gear scheme for the approximation of fractional derivatives. After a first tentative [12] for preliminary tests using semi-derivatives, we focus on (i) the analytic determination of the coefficients of the numerical scheme and (ii) a set of preliminary tests in order to derive orders of convergence. Then, we study the harmonic oscillator with fractional damping from analytical and numerical points of view.

Mathematical background
Let us introduce the Riemann-Liouville fractional integral operator of order a I a uðtÞ ¼ 1 CðaÞ where a is a real number such that 0 < a < 1 and C is the Gamma function. The corresponding differential operator arises from Eq. (1), such that which is the classical Riemann-Liouville fractional differential operator of order a. It is well known that this operator can be alternatively defined by using the Caputo [1] approach: Obviously, the equivalence of Eqs. (2) and (3) is valid only when the function u satisfies the requirement u(0) = 0. Furthermore, we assume throughout this study that u is a causal function, i.e. u(t) = 0 for t 6 0, except for the fundamental fractional differential equation studied in Section 4.2.
In this work, we are interested in numerical methods to approximate the Riemann-Liouville fractional derivative. For didactical purposes, some basic definitions are recalled before showing the generalization of two finite difference methods to fractional derivative operators.
Let u be a time dependent function known only in its discretized values u n at each time t n , where n is a positive integer. The function u n is approximated by u(t n ) with t n = nDt, where Dt, which is supposed to be fixed, is the time step. Furthermore, let us introduce a delay operator given by

Two renowned schemes for integer derivatives
Recall that the first-order derivative of the function u n can be approximated using the Euler backward formula where (Eu) n is the Euler backward operator applied to the function u at time t n , defined as follows: Analogously, we can approximate an integer derivative of u by using the Gear three-time step scheme [13] ðD 1 uÞ n % 1 Dt where G is the Gear operator defined by It should be stressed that these two schemes can be extended to the case of fractional derivatives, as we shall see below.

The Grü nwald-Letnikov scheme
The extension of the Euler backward scheme to fractional calculus generates the so-called Grü nwald-Letnikov (GL) scheme that can be represented, following the work of Oustaloup [14], by In Eq. (9), the term in brackets can be computed by using the Newton binomial formula By applying this expression to Eq. (9), one obtains the following fractional derivative operator: with the coefficients ðÀ1Þ j C j a given in terms of the Gamma function where A a jþ1 are the so-called GL-coefficients, with A a 1 ¼ 1 for any a. Then, the a-derivative of function u evaluated at time t n is approximated by the Grü nwald-Letnikov scheme [7,8]: Using the property nC(n) = C(n + 1), the GL-coefficients in Eq. (12) can be computed by the recurrence formula It is important to note that this scheme has been widely used to approximate fractional derivatives in the context of finite element method due to its simplicity when implemented numerically. Furthermore, Eq. (14) describes correctly the fading memory phenomena observed in the behavior of viscoelastic materials (see, for example, [2,4,[15][16][17]).

Outline of the G a -scheme
Based on the previous procedure used to achieve the Grü nwald-Letnikov fractional differential operator GL from the Euler backward formula, one introduces the fractional differential operator G a as which is directly obtained by evaluating the a-power of Eq. (8). The basic idea is to reproduce the formalism used to achieve GL from the Gear scheme. In this way, using (10) and (12), the Gear operator for fractional derivatives is written as Thus, the a-derivative of u at time t n can be approximated by where the coefficients A a jþ1 are given in (14) and the coefficients B j 'þ1 are calculated using the recurrence formula The sign of these coefficients is alternate following the evolution of j, while coefficients from Eq. (14) are always negative when j increases. Furthermore, one notes that B j 1 ¼ 1 for any j. The calculation of the G a -coefficients is a hard task due to cumulative numerical errors. In order to overcome such a difficulty, the method employed here consists of calculating these coefficients analytically using Symbolic Matlab Toolbox. For illustrative purposes, the reader is referred to Table 1 and Fig. 1, where the first 20 G a -coefficients are presented for three values of a: 1/3, 1/2, and 3/4. Such coefficients are related to Eq. (17) by the following expression: where g is a rational number. As one can observe, this way of writing the G a -scheme is similar to that one for the GLscheme (see Eq. (13)). Thence, from a numerical point of view, it is obvious that Eq. (19) represents a handle tool for approximating fractional derivatives when compared to Eq. (17). Moreover, even if the method used for computing the G a -coefficients avoids cumulative numerical errors, it is limited to handle only a few hundred terms due to numerical overflow. For example, only 2 9 = 512 terms are used in the calculations performed in this work (see Section 4). This would correspond to a fraction where numerator comprises 605, 392, 468 digits and denominator 609, 397, 473 digits, when a equals respectively to 1/3, 1/2, 3/4.

Results and analysis
In this section, three sets of numerical tests are performed in order to validate and analyse the G a -scheme. The analysis starts by ''static'' tests related to the power function, where one compares the proposed scheme with an exact solution proposed in literature. The second set of tests deals with the solution of a ''fundamental'' fractional differential equation whose solution is the analogous of the exponential function. Finally, in the last but not least set of tests, one presents a now classical problem in structural dynamics, which consists of a fractional damped oscillator submitted to a dependent-time excitation. Calculations carried out using a combined (GL, G a )-(Heun, Newmark) scheme are compared with a proposed formal series solution.

First elementary test cases
In this section, some numerical tests are performed in order to establish a rate of convergence for the G a -scheme. The function chosen for this purpose is the power function The basic idea is (i) to find the approximated solution for the following fractional derivative equation: using (GL, G a )-scheme, then (ii) to compare it with an exact solution f(t) such as, for the power function (20), its Riemann-Liouville fractional derivative is [18,19] Furthermore, the order of accuracy of the G a -scheme can be quantified by calculating L 2 -and L 1 -errors. For a fixed time step Dt = 1/2 m , these errors are respectively calculated by where m is a positive integer. Preliminary results related to the rate of convergence of the G a -scheme are presented in Table 2. These results are obtained from the error estimates in L 2 -and L 1 -norms, with three values of a: 1/3, 1/2 and 3/4 and various values of the power m. Moreover, the time discretization is bounded to m 6 9. All the calculations are also performed using the GL-scheme. It should be emphasized that the rate of convergence for both schemes depends strongly on the power m (evolution of columns) and weakly on the order of the fractional derivative a (evolution of rows). As for finite difference methods, the smoothness of the approximated function plays an essential role in the accuracy of the scheme. In order to show this, rows in Table 2 are split on four row-blocks as values of m increase. One focuses the analysis on second (1 6 m 6 7/4) and third (2 6 m 6 11/4) row-blocks that correspond to the accuracy range of GL-and G a -schemes. As expected, the rate of convergence for the GL-scheme does not exceed 2 (see underlined elements in second row-block). The rate is bounded between 1 and 2 for all values of a, while for G a -scheme, the rate is always superior to 2 (see underlined elements in third row-block). Moreover, the slope of error estimates on norm L 1 corresponds to m or, in terms of a, one has a rate of convergence at best of order a + 1 for GL and at best a + 2 for G a . It might be worth mentioning that such a similar phenomenon have been described and theoretically analysed and proved for a different numerical method by Diethelm et al. [20].
In order to better understand the results presented in Table 2, Fig. 2 shows the evolution of the rate of convergence derived from an error estimate in L 1 -norm in terms of m. Points A-C stand for the convergence stagnation of G a -scheme (see also Table 2 for corresponding superscripts). One observes that for same values of a and m, the convergence of the GLscheme stagnates before (points D-F). For example, for m = 9/4 and a = 1/2, GL-convergence has been unchanged since m = 3/2 while the G a -convergence goes on increasing.  Fig. 3), G a -scheme is poor accurate, i.e. it has same slope than GL-scheme in L 1 -norm. Such a result lead us to think that this is due to the non-regularity of the function t in origin. Fig. 4 presents the tests related to the semi-derivative of the function t 2 . In this case, the G a -convergence is larger or equal to 2, then being better accurate than the GL-scheme.
In Fig. 5, the semi-derivative of the function t 3 is analysed. As for the previous result, the G a -scheme is better accurate, noting that the rate of convergence of the GL-scheme is the same than in Fig. 4.
These results for a preliminary test case show that the accuracy of the G a -scheme depends on the regularity of the evaluated function in origin. The set of tests carried out in the following sections emphasizes such a numerical phenomenon.

The fundamental fractional differential equation
Consider the following elementary differential equation:  where s > 0 is a time type constant. The solution of (25), as well-known, can be explicited by using the exponential function. For fractional differential operators, the situation is analogous. Thus, the fundamental fractional differential equation can be written as The solution of this problem is given in terms of the Mittag-Leffler's fundamental ''exponential'' where E a (•) is the one-parameter Mittag-Leffler function [21], defined as Table 3 presents error estimates in L 2 -and L 1 -norms of GL-and G a -schemes for three values of a when solving Eq. (26). This test case shows a limitation of the present G a -scheme when the fractional derivative drives the governing equation, i.e., when the highest derivative term is the fractional one. In this case, both schemes are poor accurate, having the same order which equals at about to a. In Fig. 6(a), one presents exact and approximated (with GL-and G a -schemes) solutions of (26) for a = 1/2 and s = 1 (in a suitable unit system). As in previous examples (see Section 4.1), the accuracy is probably lost due to the non-regularity of the function at the origin as illustrated by the spurious points in the beginning. For this example, the error estimates in L 1norm of both schemes is at about a according to Fig. 6(b).

The fractional damped oscillator problem
Consider a fractional one-dof system submitted to a constant step load f for t > 0 with zero initial conditions. The damping is taken into account by introducing a fractional damping term or a spring-pot element in the formulation. The corresponding governing equation as well as the initial conditions are given below Table 3 Error estimates in L 2 -and L 1 -norms for various values of a where m and k are mass and stiffness constants; and cs a is a fractional damping constant with s the relaxation time and c the classical damping constant. The aim of this section is to solve the set of Eq. (29) with a direct time integration method (Newmark and Heun schemes) in conjunction with an approximation for the a-derivative D a u (G a -or GL-scheme). Furthermore, in order to validate such a combination, the approximated solution is compared to an exact solution proposed below. Finally, error estimates in L 1 -norm are performed (see formula (24)).

Exact solution
Let us introduce an original exact solution for the problem (29) such that a = p/q, where p and q are positive integers. One proposes to write this exact solution in terms of formal power series as where the corresponding S j coefficients are defined as follows: These coefficients have been obtained after some analytical calculations. More details will be presented in a coming paper. The above solution is used for the calculation of the order of convergence of our numerical scheme using a finite number of terms in the series (see the numerical application).

Algorithm
Concerning the numerical approximation of the a-derivative, one follows the strategy used by Galucio et al. in [4] with a combined GL-Newmark algorithm for structural dynamics. In other words, the displacement history arising from the aderivative approximation (damping term) is shifted to the right-hand side of Eq. (29). In this way, using (16), the governing equation in its discretized form is written as where the modified terms j and / arise from the approximation of the a-derivative: One notes that the stiffness term j is constant in time, depending only on the time step, which is supposed to be fixed. Concerning the modified loading /, it depends on the displacement history. Let us recall two classical predictor-corrector integrators: Newmark and Heun schemes. Both are one-time step and second-order accurate algorithms when used to solve second-order differential equations sufficiently smooth. Newmark is an implicit scheme while Heun is an explicit one. Below, we shall outline how to adapt both integrators to Eq. (32).

Update time step and return to 2
In all calculations performed in this work, we assume that m = k = s = f = 1 in a suitable unit system. In Table 4, as well in Figs. 7-9, one assumes that c = 1. Moreover, as one has seen in previous sections, three values of a are tested. The final time is chosen to be T = 15 for various values of time step. Moreover, 400 terms are retained in the series for the calculation of the exact solution (30).
In Table 4, error estimates in L 1 -norm are presented when integrating Eq. (29) with Newmark and Heun algorithms. These errors are measured by using Eq. (24) over two time discretizations, i.e. when the number of time steps is bounded  Fig. 7. (a) Exact and approximated solutions of (29) using Newmark scheme for a = 1/3 and Dt = T/2 6 ; (b) error estimates in L 1 -norm. between 2 7 and 2 9 . One notes that the combined G a -Newmark scheme keeps the second-order accuracy of both G a and Newmark schemes for any value of a. However, the use of a GL-Newmark algorithm decreases the order of accuracy to 1. Concerning the combination G a -or GL-Heun scheme, one might say that due to driving acceleration term in Eq. (29), the accuracy is not conserved neither for G a nor for GL. Figs. 7-9 show the evolution of the displacement for various values of a as well the error estimates on L 1 -norm, when using a combined G a -or GL-Newmark scheme. In Figs. 7-9(a), the exact solution of Eq. (29) and its corresponding numerical approximations (GL and G a methods) are presented, with a time discretization corresponding to 2 6 = 64 time steps. One can easily note that the solution obtained by using the G a -scheme is very close to the exact solution while that one obtained by the GL-method is overestimated.
Error estimates in L 1 -norms are presented in Figs. 7-9(b). In the three situations, the combined G a -Newmark scheme shows a better accuracy than the GL one. One recalls that the rates of convergence presented in Table 4 are computed with 7-9 meshes, otherwise the slopes are wrongly estimated.
It should be emphasized that the order of the fractional derivative does not affect the rate of convergence (according to Table 4 and Figs. 7-9) when using a Newmark integrator. The influence of a is observed in the mechanical behavior of the fractional damped oscillator by means of a damping factor. In other words, when a decreases, the damping and the time required to achieve the quasistatic time solution increase.
In order to show the influence of an added damping, the results presented below are computed for a fixed value of a = 1/ 2 and different values of the classical damping constant c. According to Table 5 (see also Table 4), the rate of convergence remains the same when using the Newmark scheme, while a slight difference is observed for the Heun scheme. As in previous results, using the combined G a -Newmark algorithm, the order of accuracy is 2.
For illustrative purposes, the responses of the oscillator computed with the G a -Newmark scheme are presented in Fig. 10(a) for three values of damping: c = 0.5, 1.0, 1.5 (see Eq. (29)). These results are obtained for a semi-derivative problem. Comparing to the exact solution (30), the corresponding error in L 1 -norm is plotted in Fig. 10(b). One notes that the rate of convergence remains the same for all c.

Conclusion
A numerical method based on the Gear scheme to approximate fractional derivatives is proposed. This G a -scheme is written in terms of a formal power series as the Grü nwald-Letnikov approximation. This means that the G a -scheme is no more expensive than the Grü nwald-Letnikov one. The numerical evaluation of G a -coefficients is delicate due to a bad conditioning of the recurrence formula. However, with the help of formal calculus, cumulative numerical errors are avoided.
Three sets of tests are performed. In a first time, numerical experiences based on power functions yield a L 1 error at best of order a + 1 for the Grü nwald-Letnikov scheme and a + 2 for the G a -scheme. Secondly, the fundamental fractional differential problem is treated. Results obtained show that the G a -scheme is poor accurate (of order a) due to the non-smoothness of the function at origin. Finally, a single degree-of-freedom oscillator with a fractional damping is analysed. In order to validate our numerical approach, an exact solution for the single dof problem submitted to a constant load is proposed. Two classic algorithms are used to integrate the governing equation: Newmark (average acceleration) and Heun (secondorder Runge-Kutta) schemes. The combined algorithm G a -Newmark seems to be an effective tool for dynamic problems since a two-order accuracy is obtained. Obviously, for real problems in structural dynamics, i.e. with many dof's, the effectiveness of the present method should be demonstrated.