New fractional linear multi-step methods of order four with improved stabilities

We present two implicit fractional linear multi-step methods (FLMM) of order four for fractional initial value problems. These FLMMs are of a new type that has not appeared before in the literature. The methods are obtained from the second order super-convergence of the Grünwald-Letnikov approximation of the fractional derivative at a non-integer shift point, taking advantage of the fact that the error coefficients of this super-convergence vanish not only at first order, but also at the third order terms. The weight coefficients of the methods are obtained from the Grünwald weights and hence computationally efficient compared with that of the fractional backward difference formula method of order four. The stability regions of the proposed methods are larger than that of the fractional Adams-Moulton method and the fractional backward difference formula method. Numerical results and illustrations are presented to justify results

Fractional calculus and fractional differential equations, despite their long history, have only recently gained places in science, engineering, artificial intelligence and many other fields.
Many numerical methods have been developed in the recent past for solving (1) approximately.We are interested in the numerical methods of type commonly known as fractional liner multi-step methods (FLMM).
Lubich (LUBICH, 1985) introduced a set of higher order FLMMs as convolution quadratures for the Volterra integral equation (VIE) obtained by reformulating (1) (See also eg.(DIETHELM, 2010)).The quadrature coefficients are obtained from the fractional order power of the rational polynomial of the generating functions of linear multi-step method (LMM) for ordinary differential equations (ODEs).As a particular subfamily of these FLMMs, the fractional backward difference formulas (FBDFs) were also proposed by Lubich in (LUBICH, 1986).Other forms of FLMM are the fractional trapezoidal method of order 2 and the fractional Adams methods.
In this work, we propose two fourth order implicit FLMMs of new type that does not come under the subfamilies of FLMMs listed in Section 2. The weight coefficients of the methods are obtained from the simple Grünwald weights and has an improved stability region compared to the previously known FLMMs of order four.

Prelimineries
For a sufficiently smooth function () defined for  ≥  0 , the left Riemann-Liouville (RL) fractional derivative of order  > 0 is defined by (see eg. (PODLUBNY, 1999)) where  = ⌈⌉ -the smallest integer larger than or equal to .
The left Caputo fractional derivative of order  > 0 is defined as where  () is the -th derivative of .
Often, for practical reasons, the integer ceiling  of the fractional order  is considered to be one or two.In this paper, we investigate the case of 0 <  ≤ 1 when  = 1.Further, there is no loss in generality in the assumptions  0 = 0 and (0) = 0.
In addition to the above two definitions, the Grünwald-Letnikov(GL) definition is useful for numerical approximations of fractional derivatives.
.The coefficients can be successively computed by the recurrence relation For theoretical purposes, the function () is zero extended for  < 0 and hence the infinite summation in the GL formulation (4).Practically, the upper limit of the sum is  = [/ℎ], where [•] is the integer part function.

Numerical approximations of fractional derivatives
For numerical approximation of the fractional derivative, the GL definition is commonly used by dropping the limit in (4) giving the Grunwald Approximation (GA) for a fixed step ℎ (OLDHAM; SPANIER, 1974).
where  is the shift parameter.
For an integer shift , the SGA is of order one consistency (MEERSCHAERT; TADJERAN, 2004).However, it is shown in (NASIR; GUNAWARDANA; ABEYRATHNA, 2013) that the SGA gives a second order approximation at a non-integer shift  = /2 displaying super convergence.
Some higher order Grünwald type approximations with shifts have been presented in (GU-NARATHNA; NASIR; DAUNDASEKERA, 2019) with the weight coefficients obtained from some generating functions in an explicit form according to the order and shift requirements.

Fractional linear multi-step methods
Among the several numerical methods to solve (1), we list the numerical methods that fall under the category of FLMM.
Lubich (LUBICH, 1986) presented and studied numerical approximation methods for the FIVP (1) through some convolution quadrature for the equivalent Volterra integral equation of the FIVP.
An analogous equivalent formulation for the FIVP is also given in (GALEONE; GARRAPPA, 2008) in the classical LMM form , are starting weights to compensate the reduction of order of convergence for certain class of solution functions having singular derivatives at the initial point.
This FLMM have some subclasses in the literature with generating functions of the following general forms: 1. Fractional Trapezoidal rule: The fractional trapezoidal method of order 2 (FT2) obtained from the trapezoidal rule for the ODE has the generating function It is the only method known so far in the form () = () ()

Fractional backward difference formula:
The fractional backward difference formula (FBDF) obtained from the BDF for ODE has the generating functions of the form () = (())  .
For orders 1 ≤  ≤ 6, a set of 6 FDBF methods have been obtained in (LUBICH, 1986) with polynomials corresponding to the generating polynomials of the BDF of order  given by

Rational approximation:
In (ACETO; MAGHERINI; NOVATI, 2015), a classical LMM type of approximation is proposed to obtain a class of FLMMs by rational approximations of the FBDF generating functions.These methods have generating functions in the rational polynomial form () = () () .

New FLMMs of order four
We present the main result of constructing two new FLMMs of order 4. The fractional derivative in (1a) is replaced by the super convergence approximation (8) of order two.This gives at  =   , with the error coefficient of order 2: where  2 ≡  2 () =  24 and  2 =  2 / 2 is the second differential operator.Note that the super-convergence approximation diminishes all the odd order terms of ℎ 1 , ℎ 3 , etc. Now, we perform two approximation on the super-convergence equation ( 8).First, since  − /2 is not integer for 0 <  ≤ 1, the point   − ( − /2)ℎ is not aligned with the discrete points of the computational domain {  ,  = 0, 1, ...,  }.
We approximate this non-aligned function value by an order 4 approximation using function values at the nodal points where  − =   − ℎ and Next, the differential operator  2 is approximated by the order 2 backward difference approximation given by which gives for the function  (  , (  )) =:  (  ) Substituting ( 11) and ( 12) in ( 8), dropping the error term  (ℎ 4 ), choosing ℎ = /,  ∈ N, and denoting   =  ℎ,   ≈ (  ) and   =  (  ,   ), for  = 0, 1, ..., , we obtain a new implicit FLMM approximation scheme of order 4: Again, the second differential operator  2 can also be approximated by the order 2 difference formula Substituting ( 11) and ( 14) in ( 8), we obtain another new implicit FLMM approximation scheme of order 4: The coefficients in the new FLMMs ( 13) and ( 15) are linear combinations of the Grünwald weights  ()  and thus does not involve any heavy computations.Theorem 1 The new FLMMs in (13) and (15) are consistent of order 4 and has the generating function of the form with () =  0 +  1  +  2  2 +  3  3 for both FLMMs, and To the best our knowledge, this form of generating functions has not been appeared in the previously known FLMMs in the literature.
The computational order of the method is computed by the formula where   , ℎ  are the Maximum error and the step size for the computational domain size   .Tables 1 and 2 list the results obtained in the computations for the proposed FLMM4.1 and FLMM4.2 respectively for  = 0.4, 0.6 and 0.6.

Stability regions and comparisons
For the analysis of stability of a FLMM, the analytical solution of the test problem    () = (), (0) =  0 is given by () =   (  ) 0 , where   (•) is the the Mittag-Leffler function The analytical solution () of the test problem is stable in the sense that it vanishes in the -angled region , where the angle /2 is measured from the positive real axis of the complex plane.The analytical unstable region is thus the infinite wedge with angle  complement to the analytical stability region Σ  .
For the numerical stability of FLMM, we have the following criteria: Definition 1 Let  be the numerical stability region of a FLMM.For an angle , define the wedge where  is measured from the negative real axis of the complex plane.The FLMM is said to be 1.()-stable if () ⊆ .We denote our new FLMMs in ( 13) and (15) as NFLMM4.1 and NFLMM4.2respectively.We compare the stability regions of previously established implicit FLMMs of order 4 with our new FLMMs.For this, we consider the Lubich's FBDF4 (LUBICH, 1986) and the FAM3(GALEONE; GARRAPPA, 2008) given by their respective generating functions , where  From the Dahlquist's second barrier for FLMM (see (GALEONE; GARRAPPA, 2008)), it is clear that order 4 methods are not A-stable.However, (/2)-stability and (0)-stability could be measuring tools for comparing the FLMM methods of orders higher than 2.
In Figure 1, the unstable regions of the proposed two FLMMs for various values of  are illustrated.Here, the straight lines in the figures depicts the A-stable boundaries of the analytical stability region where the left side of the lines are the analytical stability regions Σ  for various  values.Note that the unstable regions surpasses the A-stable boundaries for all values of .This confirms that the methods are not A-stable.As for the (/2)-stability, we see from the figure that the unstable regions are on the right side of the complex plane for many values of .However, there are values of  for which the unstable region peeks in to the left side as well, for example  = 1.
We numerically computed the intervals for  where the FLMMs are (/2)-stable.As we see, the NFLMM4.2 has the highest interval for  for (/2) stability followed by FBDF4, NFLMM4.1 and FAM3.The FAM3 has a far lower interval size in this sense.
Note also that the A(0)-stability of FAM3 is destroyed as  passes the (/2)-stability bound  * .The stability region for  >  * becomes bounded and falls in the left side of the complex plane (GALEONE; GARRAPPA, 2008) displaying only conditional stability.
As, for the computational efficiency, the weights of the NFLMMs have the simplest computational effort as they involve only a linear combinations of the Grünwald coefficients  ()  .Obviously, the weights of FBDF4 require computations by the Miller's formula with four previous weights.

Table 2 :
Computational order of the new FLMM4.2