Determination of Coefficients of High-Order Schemes for Riemann-Liouville Derivative

Although there have existed some numerical algorithms for the fractional differential equations, developing high-order methods (i.e., with convergence order greater than or equal to 2) is just the beginning. Lubich has ever proposed the high-order schemes when he studied the fractional linear multistep methods, where he constructed the pth order schemes (p = 2, 3, 4, 5, 6) for the αth order Riemann-Liouville integral and αth order Riemann-Liouville derivative. In this paper, we study such a problem and develop recursion formulas to compute these coefficients in the higher-order schemes. The coefficients of higher-order schemes (p = 7,8, 9,10) are also obtained. We first find that these coefficients are oscillatory, which is similar to Runge's phenomenon. So, they are not suitable for numerical calculations. Finally, several numerical examples are implemented to testify the efficiency of the numerical schemes for p = 3,…, 6.

Although there are more than six kinds of fractional derivatives, the commonly used derivatives are Riemann-Liouville, Grünwald-Letnikov, and Caputo ones. In this paper, we focus on Riemann-Liouville derivative. Under suitable conditions, the Riemann-Liouville derivative can be discretized by the discrete form of the Grünwald-Letnikov one.
In the following, we give some basic definitions, notations, and properties of the fractional calculus [1,3,7,26].
The present paper is organized as follows. We divide Section 2 into three subsections. In Section 2.1, we introduce the recursion formulas of coefficients of orders 1 and 2, which are usually used, and some properties are given. In Section 2.2, coefficients of orders from 3 to 6 are presented for reference. And coefficients of orders from 7 to 10 are displayed in Section 2.3 which are oscillatory. Such a phenomenon is similar to Runge's phenomenon. So, the high-order (more than 6th-order) schemes seem not to be suitable for numerical calculations, which is the same as the case of ordinary differential equation. In Section 3, some numerical experiments are carried out to support the computational schemes. Section 4 concludes this paper.
In most situations, we naturally use the following formula to approximate the Riemann-Liouville derivative: where ( ) are the binomial coefficients. And they have the following recurrence relationships: Obviously, ( ) 1,ℓ (ℓ = 0, 1, . . . , ) are just the first + 1 coefficients of Taylor series of the expansion of the following function: We can see that formula (10) has only the first-order accuracy if ( +) = 0. Therefore, to seek high accurate numerical methods for the fractional derivatives is of great importance. In [16], Lubich firstly proposed numerical schemes of orders 2, 3, 4, 5, and 6 for fractional linear multistep formulas. Here, it must be mentioned that the fractional linear multistep method is different from the usual linear multistep method. The former is of varied steps. That is to say, the value of the th step relies on the preceding step values 0 , 1 , . . . , −1 , which means that the number of multisteps is increasing, while the latter is of fixed number of multisteps. Under the homogeneous initial value conditions, that is, ( ) ( +) = 0 ( = 0, 1, . . . , − 1), the fractional linear multistep scheme for the Riemann-Liouville derivative has the following form: If ≥ 2, the corresponding numerical methods are often called the high-order methods. The coefficients ( ) ,ℓ in the above equations are those of the Taylor series expansions of the corresponding generating functions ( ) ( ); that is, The Scientific World Journal 3 The corresponding generating functions for = 2, . . . , 6 are given as follows [16]: From the above introduction, the key question is how to compute the coefficients ( ) ,ℓ , = 2, . . . , 10. As far as we know, there are three methods to compute these coefficients, in which one way is to use the fast Fourier transform method [19].
Here, the values ( = 0, 1, . . .) in the above formula denote the Taylor expansion coefficients of the generating functions The third method is to use the expansion of series, where the coefficients are not given in the form of recurrence relationships but the exact expressions; see the Appendix. Such analytical expressions are useful for theoretical analysis, such as stability and convergence. Here, we introduce the second method.

Remark 7.
In [15], the high-order schemes for Caputo derivative were firstly derived. Here, one can get another way to construct the high-order numerical algorithms for Caputo derivatives. If the homogeneous initial value conditions are satisfied, one has the following numerical schemes due to Lemma 5: 2.3. Determination of ,ℓ ( = 7, 8,9,10). In this subsection, we present the recursion formulas of ,ℓ ( = 7, 8, 9, 10) for reference.

6
The Scientific World Journal The value of coefficient The Scientific World Journal

Solutions
The numerical solution The exact value of RL D 0,1 x 3 Figure 11: The exact value of RL 0,1 3 and numerical solution of = 9 for = 0.5 and step length of ℎ = 1/120. The exact value of RL D 0,1 x 3 Figure 12: The exact value of RL 0,1 3 and numerical solution of = 10 for = 0.5 and step length of ℎ = 1/120.
Next, we plot the figures of coefficients of ( ) ,ℓ ( = 7, 8,9,10) to show the evolutions with ℓ. From Figures 5, 6, 7, and 8, we can see that these coefficients are violently oscillatory such that the approximation behaves like Runge's phenomenon, which is similar to the case of ordinary differential equation. So, to seek high-order (≥7th-order) schemes by this form seems not to be appropriate.

Numerical Examples
In order to verify the reasonability of the coefficients for = 3, 4, 5, 6, we give the following two numerical examples. These numerical results show that the coefficients are efficient. Example 9. Let us consider a fractional ordinary differential equation with initial values RL −1 0, The exact solution of the above equation is given by At this moment, we use the numerical formula (14) with different order to solve this equation. The absolute error and numerical convergence order are listed in Tables 5, 6, 7, and 8.
From the numerical results presented in Tables 1-8, we can see that the coefficients of the fractional linear multistep method for = 3, 4, 5, 6 are efficient.
The coefficients of ( ) ,ℓ ( = 7, 8, 9, 10) are violently oscillatory which may not be suitable for numerical calculations. Now, we take an example to show this.
We use the scheme (14) with = 7, 8, 9, 10 to numerically compute RL 0, 3 . See Figures 9, 10, 11, and 12. From these figures, we can see that the results are not numerically stable, which can be regarded as fractional Runge's phenomenon. So, it is not necessary to derive more than 6th-order schemes for Riemann-Liouville derivatives.

Conclusion
In this paper, we propose recursion formulas to compute the coefficients of the fractional linear multistep schemes. The numerical experiments have been carried out to support the derived numerical schemes. Here, we should note that the th order ( ≥ 7) schemes are not stable.

Appendix
The Computations of Coefficients ( ) ,ℓ for Case = 3, 4, 5, 6 Here, we first introduce the Fourier transform method to compute the coefficients for = 3, 4, 5, 6. Letting = and substituting it into (15) yield That is, we need compute the following integrals: Obviously, it is not an easy task to obtain the explicit analytical expressions of the coefficients ( ) ,ℓ ( = 2, 3, 4, 5, 6) if ∉ Z. These analytical expressions are seemingly complicated but are very useful for theoretical analysis, such as stability and convergence analysis. In the following, we give an effective and interesting method to obtain the explicit analytical expressions of coefficients ( ) ,ℓ ( = 2, 3, 4, 5, 6). For case = 2, it is easy to get; see [27] for more details, For cases = 3, 4, 5, 6, the details are given as follows.
By some calculations, one has ( ) 3 ( ) = ( (A.5) Comparing the above equation and (14) with = 3 gives ( ) 3,ℓ = ( The Scientific World Journal It follows from the above equation and (14) with = 4 that ( ) 4,ℓ = ( By complex computations, one obtains      The main reason to list the exact expression of the coefficients is that they are helpful for stability and convergence analysis. The coefficients for = 7, 8, 9, 10 are omitted due to inefficiency.