Analysis of New Type of Second-order Fractional Linear Multi-step Method with Improved Stability

: We present and investigate a new type of implicit fractional linear multi-step method of order two for fractional initial value problems. The method is obtained from the second-order superconvergence of the Grünwald-Letnikov approximation of the fractional derivative at a non-integer shift point. The method coincides with the backward difference method of order two for the classical initial value problem when the order of the derivative is one. The weight coefficients of the proposed method are obtained from the Grünwald weights and are hence computationally efficient compared with that of the fractional backward difference formula of order two. The stability properties are analyzed and it is shown that the stability region of the method is larger than that of the fractional Adams-Moulton method of order two and the fractional trapezoidal method. Numerical results and illustrations are presented to justify the analytical theories


Introduction
Consider the fractional initial value problem (FIVP) 0 ( ) = ( , ( )), ≥ 0 , 0 < ≤ 1, where 0 is the left Caputo fractional derivative operator of order defined in Section 2, ( , ) is a bounded function satisfying the Lipschitz condition in the second argument guaranteeing a unique solution to the problem [1]. There is no loss in considering the fractional order in the interval ∈ (0,1]. For, when 1 < ≤ = ⌈ ⌉, with appropriate initial conditions, the FIVP (1) can be formulated as a system of FIVP of order 0 < / ≤ 1, just as in the case of classical initial value problems (IVPs) with higher integer order derivatives [2]. Fractional calculus, despite its long history, has only recently gained a place in science, engineering, artificial intelligence, and many other fields [3][4][5][6].
In the recent past, many numerical methods have been developed for solving (1) approximately. We are interested in the numerical methods of the type commonly known as the fractional linear multi-step methods (FLMMs).
The basic numerical method of FLMM type of consistency order one for (1) is obtained from the Grünwald-Letnikov form for the fractional derivative [7,8]. The weight coefficients for this basic FLMM are the Grünwald weights obtained from the series of the generating function 1 ( ) = (1 − ) . Lubich [9] introduced a set of higher-order FLMMs as convolution quadrature methods for the Volterra integral equation (VIE) obtained by reformulating (1) (See also eg. [1]). The quadrature coefficients are obtained from the fractional-order power of a rational polynomial of the generating polynomials for the linear multi-step method (LMM) of classical initial value problems (IVPs). As a particular subfamily of these FLMMs, the fractional backward difference formulas (FBDFs) were also proposed by Lubich in [10]. Another particular form of FLMM type is the fractional trapezoid method of order 2.
Many researchers have utilized these formulations to construct variations of the FLMMs (see eg. [11] and the references therein). Galeone and Garrappa [12] studied some implicit FLMMs generalizing the Adams-Moulton methods for classical IVPs. Galeone and Garrappa [13] and Garrappa [14] have also investigated a set of explicit FLMMs generalizing the Adams-Bashforth methods. In this paper, we propose and analyze a new type of FLMM of order 2. The method is computationally efficient and has improved stability. We also present algorithms to solve linear and non-linear FIVPs using the proposed FLMM. We also compare the method with other known FLMMs of order 2 and show that the presented method outweighs the other methods in terms of stability and/or computational efficiency.
We classify the previously known FLMMs into subclasses based on the form of their generating functions which indicate that our FLMM presented here falls into a new subclass not encountered in the past literature. This paper is organized as follows. In Section 2, the preliminaries and previous relevant works are summarized. In Section 3, the new FLMM of order 2 is introduced along with a computational algorithm. Numerical examples for testing the method are given in Section 4. In Section 5, the stability of the method is analyzed. In Section 6, the new method is compared with other FLMMs and Section 7 draws some conclusions.

Preliminaries
For a sufficiently smooth function ( ) defined for ≥ 0 , the left Riemann-Liouville (RL) fractional derivative of order 0 < ≤ 1 is defined by (see eg. [20]) where Γ(⋅) is the Euler-Gamma function. The left Caputo fractional derivative of order > 0 is defined as In addition to the above two definitions, the Grünwald-Letnikov(GL) definition is useful for numerical approximations of fractional derivatives.
For theoretical purposes, the function ( ) is zero-extended for < 0 , hence the infinite summation in (4). Practically, the upper limit of the sum is = ⌈( − 0 )/ℎ⌉. The three definitions are equivalent under homogeneous initial conditions [8].

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.
The fundamental and widely investigated numerical approximation scheme is the Grünwald-Letnikov method (also called the fractional backward Euler method) obtained by replacing the fractional derivative operator in (1) with its GA operator ℎ of order one [8].
By choosing the discretization step ℎ appropriately to align the discrete points − ℎ with the endpoints of the problem domain [0, ] and assuming zero extension for the unknown function ( ) for < 0, the infinite sum in (6) is reduced to a finite sum. Dropping the first order error term, choosing ℎ = / , ∈ ℕ, and denoting = ℎ, ≈ ( ) and = ( , ), equation (6) gives the GL scheme ∑ =0 A shifted GL approximation of (6) is given by replacing by − and a shifted GL scheme is then given by where is the shift. The shifted scheme is also of the first order when the shift is an integer. However, at = /2, the scheme (8) gives super-convergence of order 2 [15].
However, this super-convergence scheme introduces an additional difficulty in dealing with values of the function at points not aligned with the grid points [15]. In Section 3, we modify this super-convergence scheme to construct a new scheme of order 2.
A general way to approximate the FIVP (1) is to replace the fractional derivative by the form with general weights as where the weights are to be determined for a desired order of consistency. Thus, Grunwald-type approximation schemes for the FIVP have the form expressed in conformance with the classical backward difference form (BDF) as where are coefficients in the expansion of some appropriate generating function ( ). In Lubich [10], the weights in (11) are chosen as the coefficients of the series expansion of the generating function where , are the generating polynomials of the LMM for classical IVP. It is enough to consider generating functions in the form where , , , and are polynomials, and = 1/ [10]. The generating function of an FLMM completely characterizes the approximation scheme, its stability, and the order of consistency through the following theorems.
Theorem 1 [16,17,18]. The order of an FLMM with generating function ( ) is if and only if Moreover, the approximation corresponding to ( ) satisfies, with denoting the RL fractional derivative, where ( ) is assumed to be sufficiently smooth.
Theorem 2 [10] The stability region of an FLMM with generating function ( ) is given by

Subclasses of FLMMs
We classify the FLMMs in the literature into subclasses. This classification suggests that the proposed FLMM in this paper belongs to a new subclass.

1.
The fractional trapezoid subclass: The fractional trapezoidal method of order 2 (FT2) obtained from the trapezoidal rule for the ODE has the generating function [10] It is the only method known so far in the form ( ) = ( ) with deg ( ) ≥ 1.

The FBDF subclass:
The fractional backward difference formula (FBDF) [13] obtained from the BDF for classical ODE has generating function of the form ( ) = ( ( )) . For orders 1 ≤ ≤ 6, a set of six FBDF methods have been obtained with polynomials corresponding to the generating polynomials of the BDF of order given by ( ) = ∑ =1 1 (1 − ) . The second order FBDF2, for example, is given by [10] 2 ( ) = ( 3. Fractional Adams subclass: The fractional Adams methods have the generating functions of the form , where the polynomial ( ) is determined to have order of consistency for the method [11][12][13][14]. When (0) = 0, the method is explicit and is called fractional Adams-Bashforth methods [13,14] while (0) ≠ 0 gives implicit methods which are called fractional Adams-Moulton methods (FAMs) [12]. The second order FAM method is given by the generating function.
4. Rational polynomial subclass: In [19], a classical LMM type of approximation is proposed to obtain a class of FLMMs by approximating the generating function of the FBDF methods by rational polynomials in the form . This approach, however, reduces the order of the methods and requires higher degree polynomials and , to achieve orders close to the order of FBDF considered.

A new fractional linear multi-step method
We present the main result of constructing an FLMM of order 2 which belongs to a new subclass. We need the following lemma from the Taylor series expansion.
With the notations in (7), we obtain the new implicit FLMM approximation scheme where the function values − and − +1 are properly aligned with the grid points in the computational domain.

Theorem 3 The generating function of the new implicit FLMM is given by
where ( ) = (1 + 2 ) − 2 . Moreover, the generating function satisfies confirming order 2 consistency.
Proof. The sum on the left side of (21) is manipulated, with 0 = 1 + /2 and 1 = − /2, as follows: where we have set −1 ( ) = 0. The weights which, by Theorem 1, confirms the order 2 consistency of the method.
The generating function for the new FLMM is of the form (23) and is different from those subclasses listed in subsection 2.2. Therefore, it can be considered to belong to a new subclass in the family of FLMMs. The notion of super-convergence and nodal alignment have been applied for space fractional diffusion equations in [15] and [20]. The fractional Adam-Moulton method (FAM1) [12] of order 2 derived from fractional Newton-Gregory functions and Taylor series expansion methods, can also be derived from the super-convergence of the Grünwald approximation in the form by replacing the right-hand side with the respective second-order approximations The generating function for the FAM1 is thus given by (18). Dimitrov et al. [21] formulated an order 2 scheme with super-convergence from asymptotic expansions of the superconvergence. However, the non-aligned points of super-convergence have not been re-aligned in their work. To the knowledge of the authors, the application of super-convergence of Grünwald approximation for time-fractional differential equations with re-alignments of the super-convergence point has not appeared before in the literature. The following theorem relates the proposed FLMM scheme (22) with the FBDF of order 2 by Lubich [10] and the FAM1 of order 2 by Galeone and Garrappa [12]. Proof. Immediate by substituting = 1 in (18) and (16) respectively.

Implementation
Here, we give two algorithms to compute the approximate solutions for the FIVP for linear and non-linear cases respectively using the new FLMM.

Numerical tests
We   Proof. For the boundedness of , we see that for | | ≤ 1,

Analysis of linear stability
For the symmetry about the real axis, we immediately see that ( ̅ ) = ( ). , 0] in the quadrants III and IV for 0 < ≤ 1 where sin( + 1) < 0.
Therefore, from the symmetry, the unstable region is contained in the wedge { : |arg( )| ≤ 2 } = ℂ\Σ meaning that the new FLMM is -stable.
The -stability confirms that, for 0 < ≤ 1, our new FLMM is ( /2)-stable and hence unconditionally stable. In Figure 1, the unstable regions and the A-stable tangent boundaries of the new FLMM given by the generating function (23) for fractional order values = 0.25,0.5,0.75,1 are shown. Note that the unstable region for a generating function ( ) is given by the graph (see also (15) ) = { ( ): | | ≤ 1}.

Comparison of stability regions
We compare the stability regions of previously established implicit FLMMs of order 2 with our new FLMM which we now denote by NFLMM2 for want of an abbreviation.
For this, we consider the Lubich's fractional backward difference method FBDF2 [10], the fractional Adams-Moulton method FAM1 [6] and the fractional Trapezoid rule [10,22] given by their respective generating functions in (17), (18) and (16). In Figure 2, the unstable regions for these FLMMs and our NFLMM2 are shaded for various values of . Note that the straight lines in the figures depict the boundary of the stability region of the FT2 method in which the left sides of the lines are the stability region which also corresponds to the boundary of the analytical stability regions Σ . The unstable regions of FT2 are not shaded for clarity.
The advantage of our NFLMM2, in terms of the unstable regions (UR), is that the UR of the NFLMM2 is smaller than that of the FAM1 and is very much closer to the UR of the FBDF2. Also, the UR of the FT2 is the largest among all the URs.
The computational costs of all the FLMMs of order 2 in the general form (11) are the same. Therefore, the efficiency of the methods is measured by the computational cost of the weights in (11). The weights of NFLMM2 have the simplest computational effort as they involve only a linear combination of the Grünwald weights ( ) given in (25) which can be recursively computed by (5) with only one previous weight.
In contrast, the weights of FBDF2 require computations using Miller's formula (see eg. [12] ) with two previous weights.
The weights of FAM1 require more effort as its generating function involves a rational function. However, in this case, the right-hand side of the FAM1 scheme has the form (26) with almost the same computational effort as NFLMM2. Nevertheless, the right-hand side of this scheme requires two coefficients and two values of the function ( , ) requiring an additional memory [6]. Finally, computing the weights for FT2 needs more effort as it requires the first coefficients of its generating function and FFT [12].

Conclusion
We proposed and analyzed a new FLMM of order two for FIVPs that falls under a new subclass of FLMM. The new FLMM is as -stable as the other known order two methods. However, the proposed method has a larger stability region than that of the FAM1 and FT2 methods. The computational cost of the NFLMM2 is better than that of the FBDF2 method, whereas the FAM1 requires an additional memory requirement in its iterations. Hence, the proposed method can be considered competitive with the other methods of order 2 in terms of stability and/or computational cost.

Conflict of interest
The authors declare no conflict of interest.