Analysis of Fractional Linear Multi-Step Methods of Order Four from Super-convergence

: We analyzed two implicit fractional linear multi-step methods of order four for solving fractional initial value problems. The methods are derived from the G r̈ unwald-Letnikov approximation of fractional derivative at a non-integer shift point with super-convergence. The weight coefficients of the methods are computed from fundamental G r̈ unwald weights, making them computationally efficient when compared with other known methods of order four. We also show that the stability regions are larger than those of the fractional Adams-Moulton and fractional backward difference formula methods. We present numerical results and illustrations to verify that the theoretical results obtained are indeed satisfied

( 0 ) =  0 , where   0  is the left Caputo fractional derivative operator defined in Section 2, (, ) is a function satisfying the Lipschitz condition in the second argument  which guarantees a unique solution to the problem [1,2].When  = 1, problem (1) with (2) becomes the classical initial value problem(IVP) with first order derivative.Many numerical schemes for approximately solving the FIVP (1) have been proposed in the recent past.The numerical methods referred to as fractional linear multi-step methods (FLMMs) are of particular interest.

C
The simplest and most highly investigated FLMM is the fractional Euler method (also known as the Grünwald-Letnikov method) obtained by approximating the fractional derivative   0  () in (1), after some modifications, by the Grunwald-Letnikov (GL) approximation (See Section 2).
Converting the FIVP (1) in the form of Volterra integral equation (VIE) of the second kind, Lubich [3][4][5] proposed a class of higher-order FLMMs for the VIE as convolution quadrature rules.The quadrature coefficients of the FLMMs are obtained from the fractional order power of some rational polynomials from the linear multistep methods(LMMs) for classical IVPs.Moreover, Lubich [5], suggested a set of implicit fractional backward difference formula (FBDF) methods as a subclass of these FLMMs.
Aceto [10] constructed another subclass of FLMMs by approximating Lubich's generating functions of the FBDFs by Pade approximations.However, in this class, the orders of the FLMMs are reduced compared to the source FLMMs.
The present authors proposed two new implicit FLMMs of order 4 with preliminary properties and tests presented in [11].The methods use the super convergence of the GL approximation.Earlier, the authors used super convergence to derive an FLMM of order 2 in [12].
This paper analyzes the two implicit FLMMs of order 4 presented in [11].As FLMMs of orders higher than two are not A-stable according to Dahlquist's second barrier for FIVPs [13], we analyze the stability of the methods through (/2)-stability and unconditional stability.We also show that the methods are better in stability than the FAM4 of order 4 and one of the methods is better than FBDF4 of order 4.
The computational costs of these methods have also been compared with other order 4 methods and show that the new FLMMs are computationally competitive with the FAM and simpler than that the FBDF4.This paper is organized as follows.Section 2 gathers the necessary definitions, theories, and results on factional calculus and numerical solutions of FIVPs.In Section 3, we give the main results by constructing the new FLMMs along with an algorithm to compute approximate solutions using these methods.In Section 4, we analyze the stability of the methods.Section 5 compares our methods with other known methods of order 4. In Section 6, we draw conclusions.

Preliminaries
The fundamental definitions of fractional derivatives in fractional calculus are typically presented as follows: Definition 2.1 Let () be a function defined in the interval domain [ 0 , ] and is sufficiently smooth to hold the following: 1.When  ∈  1 ([ 0 , ]), the Riemann-Liouville (RL) fractional integral of order  > 0 of () is defined as 2. The RL fractional derivative of order  > 0 is defined by where  = ⌈⌉ is the integer ceiling of , and Γ(⋅) denotes Euler's gamma function.

5.
A shifted form of GL fractional derivative is also available [14]: where  is the shift which is often taken to be an integer, but any real shift is valid.

Remark:
The fractional derivatives given above are called the left fractional derivatives as there are also their right counterparts.For details see [15,16] and the references therein.The fractional integrals and derivatives are related such that the RL and Caputo derivatives are two left inverses of the RL integral [1]: However, the two fractional derivatives in (4) and ( 5) are related by where Hence, the RL and Caputo fractional derivatives are equal under the homogeneous initial conditions  () ( 0 ) = 0,  = 0,1, … ,  − 1 [1,2].The GL fractional derivative in (6) and its shifted counterpart in (7) are also equivalent to the Caputo derivative under homogeneous conditions and are often utilized as tools for numerical approximations of fractional derivatives.

Approximation of fractional integrals and derivatives
To approximate the fractional integrals and derivatives, the involving domain [ 0 , ] is discretized into a computational domain with uniformly spaced discrete points   =  0 + ,  = 0,1, … ,  with a fixed step size  = ( −  0 )/.The fractional integral can be approximated by using a quadrature rule as the sum of weighted function values at the discrete points of the involved integrating domain.Common quadrature rules used in this sense are the rectangular and trapezoid rules [17][18][19].Lubich [5] introduced a convolution quadrature approximation formula for the fractional integral where the weights   are obtained from the power series expansion of the generating function (, ) being the pair of generating polynomials of the LMM for classical IVPs [1].The order of consistency for the FLMM is the same as that of the underlying LMM.As for the approximation of fractional derivatives, the fundamental approximation for the RL fractional derivative is obtained from the GL (or generally the shifted GL) definition by simply dropping the limit for a fixed step size .This gives an order one approximation  ,  () with an integer shift  [5,21].
where the initial value  0 has been subtracted from () to satisfy the homogeneous initial condition so that the different definitions coincide.However, for the particular non-integer shift  = /2, the above Grünwald approximation gives order 2 displaying super convergence [22].
Analogous to the convolution quadrature approximation for fractional integral, fractional derivatives can also be approximated by convolution quadrature formula in a similar form , where   are the coefficients of a generating function ().The order of consistency of an FLMM can be obtained from its generating function through the following theorem.Theorem 2.2 [14,23,24] Let () be the generating function of an approximation in the shifted form of the fractional derivative   0  (), where () is sufficiently smooth.The order of the shifted approximation with shift r is  if and only if ) where   ≡   (, ) are the coefficients of the series expansion of ():

FLMM scheme
The general form of an FLMM scheme for the FIVP in ( 1) and ( 2) is given by where   and   are the coefficients of the generating functions and   and   denote   ≈ (  ) and   = (  ,   ).
In numerical computations using the FLMMs of order more than one, the intended order  is achieved only for a certain class of functions, specifically for functions of the form () =   (),  ≥ , with () analytic [5].However, for  < , the order is reduced to (ℎ  ) only.To remedy this order reduction, an additional sum is introduced in (11) to have the approximation scheme The starting weights  , in ( 13) are to compensate for the reduced order of convergence.However, computing the starting weights poses many difficulties in practice.Since the starting weights do not affect the stability or convergence of the solution, we do not include them in the computation and analysis in the subsequent sections.For some developments on the starting weights, we refer to [20,26,27].

Stability
For the analysis of the stability of an FLMM, consider the linear test problem for which the analytical solution is () =   (  ) 0 , where is the Mittag-Leffler function.For analytical stability, the solution () of the test problem ( 14) is stable in the sense that the series of () converges in the region Σ  = { ∈ : |()| > /2}.The unstable region is then the infinite wedge with angle .(See Figure 1).3. (/2)-stable when the entire left half of the complex plane is included in . 4. Unconditionally stable if it is (0)-stable.That is, the negative real line is included in S. The stability region of an FLMM is also characterized by its generating function: Theorem 2.4 [5] The stability region of an FLMM with generating function () is given by where   = {(): || ≤ 1} is the unstable region.

New FLMMs of order 4
In this section, we give the construction of these methods.Denote by   (ℝ) the class of functions having continuous n th derivatives.

Lemma 3.2
The second derivative of a function () can be approximated by the backward difference forms of order 2 as Proof.Taylor series expansions.∎  in ( 9) is an even function since Hence, the odd order terms of the series expansion of () are zero.Moreover, expanding for the first few terms reveals  0 = 1,  2 =  24 .∎ From Theorem 2.2 again, we have from (10), Writing   0 +2 = We approximate ( − + /2) by ( 15 Dropping the error term, with the notations in (12), equation ( 22) gives an implicit FLMM scheme Again, approximating the second derivative in ( 20) by (18) in Lemma 3.2, with the same operations and notations, we get another implicit FLMM: For brevity of presentation, we call these FLMMs in ( 23) and ( 24) as NFLMM4.1 and NFLMM4.2respectively.

Analysis of linear stability
The generating functions of the new implicit schemes NFLMM4.1 and NFLMM4.2 are given by: respectively, where (),  4.1 () and  4.2 () are polynomials for which the coefficients are given by ( 26),( 27) and (28) repectively.The following elementary results on complex numbers is useful.

Noticing that
Proof.For the boundedness, we see that the numerator part of (; ), , where we have used the facts that  0 ,  2 > 0 and  1 ,  3 < 0. For the denominator part for NFLMM4.1, Hence, for 0 <  ≤ 1 and for || ≤ 1, Since  4. (; ) =  4.x (; ), we immediately see that the unstable regions are symmetry about the real axis.∎ Since the NFLMM4.xare of order 4, the Dahlquist barrier for FIVPs tells us that they are not A-stable [13] .
Therefore, we look for the  (   )-stable boundaries.

Comparisons
In this section, we compare the order 4 NFLMM4.xwith FAM3 and FBDF4 for their performances in terms of computations and stability.
The linear equation was solved by using the schems FAM3, FBDF4, NFLMM4.1 and NFLMM4.2 in ( 23) and (24) with different values of fractional orders  = 0. Since all the four methods are of order 4, the computational solutions for all choices of discretization are expected to be nearly the same.
As for the computational cost, the weights of the NFLMM4.xneed only a linear combination of the Grunwald coefficients   () , that have the simplest computational cost.The weights of FBDF4 obviously require computations using Miller's formula with four prior weights.
As for the memory requirement, the FAM3 requires keeping all the   values stored during the iteration to be used on the right side of scheme (25).The NFLMM4.x require only the last three or four values of   as in (23) and (24).

Comparison of stability
We compare the stability regions of the four implicit FLMMs.The generating functions FBDF4 and FAM3 are provided below:  4 () = (  Since the FLMM methods with orders greater than 2 are not A-stable, the  (  2 )-stable and (0)-stable could be used as comparison tools of those methods.As shown in Figure 3, the unstable region for many values of  is on the right side of the complex plane.However, there are some  values for which the unstable region also extends to the left side, such as  = 1.The intervals for  where the FLMMs are  ( )-stability are given in Table 4. )-stability the NFLMM4.2 has the highest interval for  followed by FBDF4, NFLMM4.1 and FAM3.The lower interval size is gained for FAM3.Also note that as  approaches the  (  2 )-stability bound  * , FAM3's A(0)stability vanishes.For  >  * the stability region becomes bounded and falls on the left complex plane, resulting in only conditional stability.

Conclusion
We analyzed and compared the new fractional linear multistep methods NFLMM4.1 and NFLMM4.2 with FAM3 and FBDF4 methods of order 4. We see that NFLMM4.2 is  (  2 )-stable over a wider fractional-order interval while FAM3 displays  (  2 )-stablity for a small range of  values.Furthermore, the proposed methods have lower computing cost and minimal storage compared to FAM3.

Figure 2 .
Figure 2. Stability regions of NFLMM4.x for some  values.In Figure 2, the unstable regions of the two methods NFLMM4.x,x=1,2 are shown for different fractional orders 0 <  ≤ 1.The straight lines in the figures represent the stability region boundaries of the methods, where the stability regions are shown on the left side of the lines.These lines also correlate to the analytical stability region's boundary Σ  .It is clear from the figure that the unstable regions surpass the A-stable boundaries for all values of .Thus, the methods are verified to be not A-stable.The regions in blue are the unstable regions for the threshold values  4, * which indicate  (  2 () = (, ()) = () + (), 0 ≤  ≤ 1, 0 <  ≤ 1, with the initial condition (0) = 0, where () = Γ(+1) Γ(+1−)  − − Γ() Γ(−)  −1− +   −  −1 with  = −1.

Table 1 .
4, 0.6 and 0.8.The problem is computed on the domain [0,1], with   = 2  ,  = 3,4, ...,11 as the number of subintervals and   = 1  as the step size.The maximum errors   for step size   are compared for the methods.The order of the method NFLMM4.1 is computed by the formula The orders of other methods are nearly the same and are not presented.The results obtained are listed in Table1, 2 and 3 respectively.Comparing Maximum errors for  = 0.4.