SQU Journal for Science, 2022, 27(2),107-118 DOI:10.53539/squjs.vol27iss2pp107-118 Sultan Qaboos University 107 Analysis of New Type of Second-order Fractional Linear Multi-step Method with Improved Stability H.M. Nasir* and Khadija Al-Hassani Department of Mathematics, Sultan Qaboos University, P.O. Box: 36, Al Khod, Muscat, Sultanate of Oman. *Email: nasirh@squ.edu.om. ABSTRACT: 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. Keywords: Grünwald approximation; Generating functions; Fractional Adams-Moulton methods; Backward difference method; Super-convergence; Stability regions. محّسنتحليل نوع جديد من الطريقة الكسرية الخطية متعددة الخطوات من الدرجة الثانية مع استقرار حنيفة محمد ناصر و خديجة الحسني متعددة الخطوات من الدرجة الثانية لمشاكل القيمة األولية الكسرية. يتم الخطية الكسرية نوع جديد من الطريقة الضمنية في هذا البحث نقدم :صلخمال وقد تبين ان هذه غير صحيحة. ذات قيمة إزاحة عند نقطة هالكسري هللمشتق جرنوالد الدرجة الثانية لتقريب ذوالفائق الحصول على الطريقة من التقارب واحدًا. يتم الحصول على هالقيمة األولية الكالسيكية عندما يكون ترتيب المشتق ئلهالعكسي من الرتبة الثانية لمستتطابق مع طريقة االختالف الطريقة ، وبالتالي فهي فعالة من الناحية الحسابية مقارنة بمعادلة الفرق الجزئي المتخلف من الدرجة الثانية. تم دجرنوال معامالت الوزن للطريقة المقترحة من أوزان الكسرية من الدرجة الثانية وطريقة شبه مولتن -أدمز طريقةمنطقة االستقرارلأكبر من المقترحة يل خصائص الثبات وتبين أن منطقة االستقرار للطريقةتحل . عرض النتائج العددية والرسوم التوضيحية لتبرير النظريات التحليلية وقد تمالمنحرف الجزئية. .االستقرار، منطقة التقارب الفائق ، الكسرية ، طريقة االختالف العكسيمولتن -ق أدمزطر، المولدة الدالة جرنوالد، تقريبات :مفتاحيةالكلمات ال H.M. NASIR and KHADIJA AL-HASSANI 108 1. Introduction Consider the fractional initial value problem (FIVP) 𝑡0 𝐶 𝐷𝑡 𝛽 𝑦(𝑡) = 𝑓(𝑡, 𝑦(𝑡)), 𝑡 ≥ 𝑡0, 0 < 𝛽 ≤ 1, (1) 𝑦(𝑡0) = 𝑦0, 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-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. 2. 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]) 𝑡0 𝑅𝐿 𝐷𝑡 𝛽 𝑦(𝑡) = 1 Γ(1−𝛽) 𝑑 𝑑𝑡 ∫ 𝑡 𝑡0 𝑦(𝜏) (𝑡−𝜏)𝛽 𝑑𝜏, 0 < 𝛽 ≤ 1, (2) where Γ(⋅) is the Euler-Gamma function. The left Caputo fractional derivative of order 𝛽 > 0 is defined as 𝑡0 𝐶 𝐷𝑡 𝛽 𝑦(𝑡) = 1 Γ(1−𝛽) ∫ 𝑡 𝑡0 𝑦′(𝜏) (𝑡−𝜏)𝛽 𝑑𝜏, 0 < 𝛽 ≤ 1. (3) In addition to the above two definitions, the Grünwald-Letnikov(GL) definition is useful for numerical approximations of fractional derivatives. 𝑡0 𝐺𝐿 𝐷𝑡 𝛽 𝑦(𝑡) = lim ℎ→0 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑦(𝑡 − 𝑘ℎ), (4) where 𝑔𝑘 (𝛽) = (−1)𝑘 Γ(𝛽+1) Γ(𝛽−𝑘+1)𝑘! are the Grünwald weights and are the coefficients of the series expansion of the Grünwald generating function ANALYSIS OF A NEW TYPE OF SECOND-ORDER FRACTIONAL LINEAR MULTI-STEP METHOD 109 𝜔1(𝑧) = (1 − 𝑧) 𝛽 = ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑧𝑘 . The coefficients 𝑔𝑘 (𝛽) can be successively computed by the recurrence relation 𝑔0 (𝛽) = 1, 𝑔𝑘 (𝛽) = (1 − 𝛽+1 𝑘 ) 𝑔𝑘−1 (𝛽) , 𝑘 = 1,2, . .. . (5) 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]. 2.1 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]. 𝛿ℎ 𝛽 𝑦(𝑡): = 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑦(𝑡 − 𝑘ℎ) = 𝑓(𝑡, 𝑦) + 𝑂(ℎ). (6) 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 𝑓𝑛 = 𝑓(𝑡𝑛, 𝑦𝑛 ), (7) equation (6) gives the GL scheme ∑ 𝑛𝑘=0 𝑔𝑘 (𝛽) 𝑦𝑛−𝑘 = ℎ 𝛽 𝑓𝑛, 𝑛 = 1,2, . . . , 𝑁. A shifted GL approximation of (6) is given by replacing 𝑘 by 𝑘 − 𝑟 and a shifted GL scheme is then given by ∑ 𝑛𝑘=0 𝑔𝑘 (𝛽) 𝑦𝑛−𝑘+𝑟 = ℎ 𝛽 𝑓𝑛, 𝑛 = 1,2, . . . , 𝑁, (8) 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]. 𝛿ℎ,𝛽/2 𝛽 𝑦(𝑡): = 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑦(𝑡 − (𝑘 − 𝛽/2)ℎ) = 𝑓(𝑡, 𝑦) + 𝑂(ℎ2). (9) 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 Ωℎ 𝛽 𝑦(𝑡): = 1 ℎ𝛽 ∑ ∞𝑘=0 𝑤𝑘 𝑦(𝑡 − 𝑘ℎ), (10) 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 ∑ 𝑛𝑘=0 𝑤𝑘 𝑦𝑛−𝑘 = ℎ 𝛽 𝑓𝑛, 𝑛 = 1,2, . . ., (11) 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 𝑤(𝑧) = ( 𝜌(1/𝑧) 𝜎(1/𝑧) ) 𝛽 , (12) where 𝜌, 𝜎 are the generating polynomials of the LMM for classical IVP. It is enough to consider generating functions in the form 𝛿(𝜉) = ( 𝑎(𝜉) 𝑏(𝜉) ) 𝛽 𝑝(𝜉) 𝑞(𝜉) , (13) H.M. NASIR and KHADIJA AL-HASSANI 110 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 1 𝑥𝛽 𝛿(𝑒 −𝑥 ) = 1 + 𝑂(𝑥 𝑝). (14) Moreover, the approximation corresponding to 𝛿(𝜉) satisfies, with 𝐷𝑡 𝛽 denoting the RL fractional derivative, 𝛿ℎ 𝛽 𝑦(𝑡) = 𝐷𝑡 𝛽 𝑦(𝑡) + ℎ𝑝 𝑎𝑝(𝛽)𝐷𝑡 𝛽+𝑝 𝑦(𝑡) + ℎ𝑝+1𝑎𝑝+1(𝛽)𝐷𝑡 𝛽+𝑝+1 𝑦(𝑡)+. . ., where 𝑦(𝑡) is assumed to be sufficiently smooth. Theorem 2 [10] The stability region of an FLMM with generating function 𝛿(𝜉) is given by 𝑆 = {𝛿(𝜉): |𝜉| > 1}. (15) 2.2. 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] 𝛿𝐹𝑇2(𝜉) = (2 1−𝜉 1+𝜉 ) 𝛽 . (16) It is the only method known so far in the form 𝛿(𝜉) = ( 𝑎(𝜉) 𝑏(𝜉) ) 𝛽 with deg𝑏(𝜉) ≥ 1. 2. 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 2 − 2𝜉 + 1 2 𝜉2) 𝛽 . (17) 3. Fractional Adams subclass: The fractional Adams methods have the generating functions of the form 𝛿(𝜉) = (1−𝜉)𝛽 𝑞(𝜉) , where the polynomial 𝑞(𝜉) is determined to have order 𝑝 of consistency for the method [11— 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. 𝛿𝐹𝐴𝑀1(𝜉) = (1−𝜉)𝛽 (1− 𝛽 2 )+ 𝛽 2 𝜉 . (18) 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. 3. 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. ANALYSIS OF A NEW TYPE OF SECOND-ORDER FRACTIONAL LINEAR MULTI-STEP METHOD 111 Lemma 1 Let 𝑦(𝑡) ∈ 𝐶1[𝑎, 𝑏] and 𝑦′′(𝑡) exists. Then, for 𝜇 ∈ (𝑎, 𝑏) and ℎ > 0, we have 𝑦(𝑡 + 𝜇ℎ) = (1 + 𝜇)𝑦(𝑡) − 𝜇𝑦(𝑡 − ℎ) + 𝑂(ℎ2). (19) The fractional derivative of the FIVP (1) (assuming with no loss 𝑡0 = 0 and 𝑦(𝑡0) = 0) is replaced by the approximation with super convergence (9) of order 2. This gives, at 𝑡 = 𝑡𝑛, 𝛿ℎ,𝛽/2 𝛽 𝑦(𝑡𝑛) = 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑦(𝑡𝑛−𝑘+𝛽/2) = 𝑓(𝑡𝑛, 𝑦(𝑡𝑛)) + 𝑂(ℎ 2). (20) Since 𝛽/2 is not an integer for 0 < 𝛽 ≤ 1, the point 𝑡𝑛−𝑘+𝛽/2 in (20) is not aligned with the discrete points in the computational domain {𝑡𝑚, 𝑚 = 0,1, . . . , 𝑁}. Using Lemma 1 with 𝑡 = 𝑡𝑛−𝑘 and 𝜇 = 𝛽/2, we replace 𝑦(𝑡𝑛−𝑘+𝛽/2) by (19). 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) [(1 + 𝛽 2 ) 𝑦(𝑡𝑛−𝑘 ) − 𝛽 2 𝑦(𝑡𝑛−𝑘−1)] = 𝑓(𝑡𝑛, 𝑦(𝑡𝑛)) + 𝑂(ℎ 2). (21) With the notations in (7), we obtain the new implicit FLMM approximation scheme Δℎ,𝛽/2 𝛽 𝑦𝑛 : = ∑ 𝑛 𝑘=0 𝑔𝑘 (𝛽) [(1 + 𝛽 2 ) 𝑦𝑛−𝑘 − 𝛽 2 𝑦𝑛−𝑘−1] = ℎ 𝛽 𝑓𝑛, 𝑛 = 1,2, ⋯, (22) 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 𝛿(𝜉) = (1 − 𝜉)𝛽 𝑝(𝜉), (23) where 𝑝(𝜉) = (1 + 𝛽 2 ) − 𝛽 2 𝜉. Moreover, the generating function satisfies 1 𝑥𝛽 𝛿(𝑒 −𝑥 ) = 1 + 𝑂(𝑥 2) confirming order 2 consistency. Proof. The sum on the left side of (21) is manipulated, with 𝑝0 = 1 + 𝛽/2 and 𝑝1 = −𝛽/2, as follows: ∑ ∞𝑘=0 𝑔𝑘 (𝛽) (𝑝0𝑦𝑛−𝑘 + 𝑝1𝑦𝑛−𝑘−1) = 𝑝0 ∑ ∞ 𝑘=0 𝑔𝑘 (𝛽) 𝑦𝑛−𝑘 + 𝑝1 ∑ ∞ 𝑘=0 𝑔𝑘 (𝛽) 𝑦𝑛−𝑘−1 = 𝑝0 ∑ ∞ 𝑘=0 𝑔𝑘 (𝛽) 𝑦𝑛−𝑘 + 𝑝1 ∑ ∞ 𝑘=1 𝑔𝑘−1 (𝛽) 𝑦𝑛−𝑘 = ∑ ∞ 𝑘=0 (𝑝0𝑔𝑘 (𝛽) + 𝑝1𝑔𝑘−1 (𝛽) )𝑦𝑛−𝑘 , (24) where we have set 𝑔−1 (𝛽) = 0. The weights 𝑤𝑘 = 𝑝0𝑔𝑘 (𝛽) + 𝑝1𝑔𝑘−1 (𝛽) , 𝑘 = 0,1, . .. (25) are the coefficients of the generating function 𝛿(𝜉) = 𝑝0(1 − 𝜉) 𝛽 + 𝑝1𝜉(1 − 𝜉) 𝛽 = (1 − 𝜉)𝛽 (𝑝0 + 𝑝1 𝜉). Moreover, we have 1 𝑥𝛽 𝛿(𝑒 −𝑥 ) = 1 − 𝛽(3𝛽+5) 24 𝑥 2 + 𝑂(𝑥 3) 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 1 ℎ𝛽 ∑ ∞𝑘=0 𝑔𝑘 (𝛽) 𝑦(𝑡𝑛−𝑘 ) = 𝑓𝑛−𝛽/2 + 𝑂(ℎ 2) by replacing the right-hand side with the respective second-order approximations 𝑓𝑛−𝛽/2 = (1 − 𝛽 2 ) 𝑓𝑛 + 𝛽 2 𝑓𝑛−1 + 𝑂(ℎ 2). (26) H.M. NASIR and KHADIJA AL-HASSANI 112 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 super- convergence. 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]. Theorem 4 When 𝛽 → 1, 1. the new FLMM converges to the BDF2 method of order 2 for the classical ODE by generating polynomial 𝛿(𝜉) = 3 2 − 2𝜉 + 1 2 𝜉2. 2. the FAM1 converges to the fractional trapezoid method of order 2 with generating function 𝛿(𝜉) = 2 1−𝜉 1+𝜉 . Proof. Immediate by substituting 𝛽 = 1 in (18) and (16) respectively. 3.1. 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. For brevity, we use the following notations: For a sequence 𝑎 = {𝑎𝑘 }, the vector slice [𝑎𝑖 , 𝑎𝑖 + 1, … , 𝑎𝑗 ] is denoted by 𝑎𝑖:𝑗 . The convolution of two vectors 𝑎, 𝑏 of size 𝑛 + 1 is denoted by 𝑎 ∗ 𝑏 = ∑ 𝑛 𝑘=0 𝑎𝑘 𝑏𝑛−𝑘. We reformulate the new FLMM scheme (22) with (24) and (25) as ∑ 𝑛𝑘=0 𝑤𝑘 𝑦𝑛−𝑘 = 𝑤0:𝑛 ∗ 𝑦0:𝑛 = 𝑤0𝑦𝑛 + 𝑤1:𝑛 ∗ 𝑦0:𝑛−1 = ℎ 𝛽 𝑓𝑛. (27) In the case of linear FIVP, we have 𝑓(𝑡, 𝑦) = 𝜆𝑦(𝑡) + 𝑠(𝑡) for some constant 𝜆 and function 𝑠(𝑡). The scheme (27) for this case, with 𝑠𝑛 = 𝑠(𝑡𝑛), is then 𝑤0𝑦𝑛 + 𝑤1:𝑛 ∗ 𝑦0:𝑛−1 = ℎ 𝛽 (𝜆𝑦𝑛 + 𝑠𝑛 ) which gives 𝑦𝑛 = 1 𝑤0 − 𝜆ℎ 𝛽 [ℎ𝛽 𝑠𝑛 − 𝑤1:𝑛 ∗ 𝑦0:𝑛−1], 𝑛 = 1,2, . .. . Algorithm 1 is given for the linear FIVP. Algorithm 1 (For linear FIVP) 1. Define 𝑠(𝑡) (for 𝑓(𝑡, 𝑦) = 𝜆𝑦 + 𝑠(𝑡)). 2. Input 𝛽, 𝜆, ℎ, and 𝑁. 3. Define 𝑔 with 𝑔0 = 1, 𝑔𝑛 = (1 − 𝛽+1 𝑛 ) 𝑔𝑛−1, 𝑛 = 1,2, … , 𝑁. 4. Define 𝑝 = [1 + 𝛽/2, −𝛽/2]. 𝑤0:𝑁 == {𝑝 ∗ 𝑔𝑘−1:𝑘 , 𝑘 = 0,1, … , 𝑁}. 5. Define array 𝒚 = {𝑦𝑘 , 𝑘 = 0,1, … , 𝑁} and set 𝑦0 = 0. 6. For 𝑛 = 1,2, … , 𝑁 𝑦𝑛 = 1 𝑤0 − 𝜆ℎ 𝛽 [ℎ𝛽 𝑠(𝑡𝑛) − 𝑤1:𝑛 ∗ 𝑦0:𝑛−1]. 7. Return 𝒚. For non-linear FIVP, the non-linear equation (27) in 𝑦𝑛 needs to be solved for the unknown 𝑦𝑛 . The Newton- Raphson method numerically solves this with an initial seed 𝑦𝑛,0 = 𝑦𝑛−1. Algorithm 2 is given for the non-linear FIVP. Algorithm 2 (For non-linear FIVP) 1. Define 𝑓(𝑡, 𝑦), 𝑓𝑦 (𝑡, 𝑦). 2. Input 𝛽, ℎ, and 𝑁. 3. Define 𝑔 with 𝑔0 = 1, 𝑔𝑛 = (1 − 𝛽+1 𝑛 ) 𝑔𝑛−1, 𝑛 = 1,2, … , 𝑁. 4. Define 𝑝 = [1 + 𝛽/2, −𝛽/2] . 𝑤0:𝑁 = {𝑝 ∗ 𝑔𝑘−1:𝑘 , 𝑘 = 0,1, … , 𝑁}. 5. Define array 𝒚 = {𝑦𝑘 , 𝑘 = 0,1, … , 𝑁} and set 𝑦0 = 0. 6. 𝑇𝑜𝑙𝑒𝑟𝑒𝑛𝑐𝑒 = 10−15, 𝐸𝑟𝑟𝑜𝑟 = 108 . 7. For 𝑛 = 1,2, … , 𝑁 𝑥0 ← 𝑦𝑛−1. 𝑐𝑛 = 𝑤1:𝑛 ∗ 𝑦0:𝑛−1. 𝐸𝑟𝑟𝑜𝑟 > 𝑇𝑜𝑙𝑒𝑟𝑎𝑛𝑐𝑒 𝐹 = 𝑤0𝑥0 − ℎ 𝛽 𝑓(𝑡𝑛, 𝑥0) + 𝑐𝑛 . 𝐽𝐹 = 𝑤0 − ℎ 𝛽 𝑓𝑦(𝑡𝑛, 𝑥0). 𝑦𝑛 = 𝑥0 − 𝐹 𝐽𝐹 . 𝐸𝑟𝑟𝑜𝑟 = |𝑦𝑛 − 𝑥0|. 𝑥0 ← 𝑦𝑛 . 8. return 𝒚. ANALYSIS OF A NEW TYPE OF SECOND-ORDER FRACTIONAL LINEAR MULTI-STEP METHOD 113 Table 1. Computational order of the new FLMM for example 1. 𝛽 = 0.4 𝛽 = 0.8 𝛽 = 1.0 𝑀 Max. Error Order Max Error Order Max Error Order 8 6.533e-03 – 1.803e-02 – 2.538e-02 – 16 1.882e-03 1.79544 5.319e-03 1.76094 7.569e-03 1.74543 32 5.052e-04 1.89743 1.449e-03 1.87619 2.078e-03 1.86499 64 1.309e-04 1.94874 3.783e-04 1.93741 5.448e-04 1.93116 128 3.330e-05 1.97438 9.665e-05 1.96858 1.395e-04 1.96533 256 8.400e-06 1.98720 2.443e-05 1.98427 3.530e-05 1.98261 512 2.109e-06 1.99360 6.140e-06 1.99213 8.879e-06 1.99130 1024 5.285e-07 1.99680 1.539e-06 1.99606 2.227e-06 1.99564 2048 1.323e-07 1.99840 3.853e-07 1.99803 5.575e-07 1.99782 4096 3.309e-08 1.99920 9.640e-08 1.99902 1.395e-07 1.99891 Table 2. Computational order of the new FLMM for example 2. 𝛽 = 0.4 𝛽 = 0.8 𝛽 = 1.0 𝑀 Max. Error Order Max Error Order Max Error Order 8 1.698e-01 – 7.835e-02 – 6.985e-02 – 16 2.779e-02 2.61128 1.978e-02 1.98599 1.769e-02 1.98155 32 6.648e-03 2.06349 5.060e-03 1.96667 4.466e-03 1.98563 64 1.663e-03 1.99866 1.286e-03 1.97645 1.122e-03 1.99286 128 4.186e-04 1.99047 3.245e-04 1.98628 2.812e-04 1.99660 256 1.052e-04 1.99271 8.155e-05 1.99260 7.037e-05 1.99836 512 2.638e-05 1.99566 2.044e-05 1.99616 1.760e-05 1.99920 1024 6.605e-06 1.99764 5.117e-06 1.99804 4.402e-06 1.99960 2048 1.653e-06 1.99877 1.280e-06 1.99901 1.101e-06 1.99980 4096 4.133e-07 1.99938 3.202e-07 1.99950 2.752e-07 1.99990 4. Numerical tests We used the new FLMM to compute approximate solutions of the FIVP (1) with a linear and a non-linear source function 𝑓(𝑡, 𝑦) in the time interval [0,1]. Example 1: 𝐶 𝐷0 𝛽 𝑦(𝑡) = Γ(𝑚 + 1) Γ(𝑚 + 1 − 𝛾) 𝑡 𝑚−𝛾 − Γ(𝑚) Γ(𝑚 − 𝛾) 𝑡 𝑚−1−𝛾 + 𝜆𝑦(𝑡) + 𝑡 𝑚 − 𝑡 𝑚−1, 𝑦(0) = 0, where 𝜆 = −1 and we set 𝑚 = 5. The exact solution is given by 𝑦(𝑡) = 𝑡 𝑚 − 𝑡 𝑚−1. Example 2: 𝐶 𝐷0 𝛽 𝑦(𝑡) = Γ(2𝛽 + 5) Γ(𝛽 + 5) 𝑡𝛽+4 − 240 Γ(6 − 𝛽) 𝑡 (5−𝛽) + (𝑡 2𝛽+4 − 2𝑡 5)2 − 𝑦(𝑡)2, 𝑦(0) = 0, with exact solution 𝑦(𝑡) = 𝑡 2𝛽+4 − 2𝑡 5. The problems are solved for fractional order values 𝛽 = 0.4,0.8 and 1.0. The computational domain for both problems are {𝑡𝑛 = 𝑛/𝑀, 𝑛 = 0,1, ⋯ , 𝑀} and step size ℎ = 1/𝑀, where 𝑀 is the number of subintervals of the problem domain [0,1]. The problems were solved for 𝑀 = 2𝑗 , 𝑗 = 3,4, . . . ,12. The computational order of the new FLMM method is computed by the formula 𝑝𝑗+1 = log(𝐸𝑗+1/𝐸𝑗 )/log(ℎ𝑗+1/ℎ𝑗 ) where 𝐸𝑗 , ℎ𝑗 are the merror and the step size for 𝑀 = 2 𝑗. H.M. NASIR and KHADIJA AL-HASSANI 114 Tables 1 and 2 list the maximum errors and computational orders for Examples 1 and 2 respectively for various grid sizes 𝑀 = 8,16, . . . ,4096. The computational orders confirm the theoretical order 2 of the new FLMM. 5. Analysis of linear stability For the analysis of the stability of an FLMM, we have the following preparations. The analytical solution to the test problem 𝐶 𝐷𝑡 𝛽 𝑦(𝑡) = 𝜆𝑦(𝑡), 𝑦(0) = 𝑦0 is given by 𝑦(𝑡) = 𝐸𝛽 (𝜆𝑡 𝛽 )𝑦0 , where 𝐸𝛽 (⋅) is the Mittag-Leffler function 𝐸𝛽 (𝑥) = ∑ ∞ 𝑘=0 𝑥 𝑘 Γ(𝛽𝑘 + 1) . The analytical solution 𝑦(𝑡) of the test problem is stable in the sense that it vanishes in the 𝛽𝜋-angled region Σ𝛽 = {𝜉 ∈ ℂ: |arg(𝜉)| > 𝛽𝜋 2 }, where the angle 𝛽𝜋/2 is measured from the positive real axis. The analytical unstable region is thus the infinite wedge {𝜉 ∈ ℂ: |arg(𝜉)| ≤ 𝛽𝜋 2 } = ℂ\Σ𝛽 . For the numerical stability of FLMM, we have the following criteria: Definition 1 Let 𝑆 be the numerical stability region of an FLMM. For an angle, 𝛼, define the sector 𝑆(𝛼) = {𝜉: |𝜋 − arg(𝜉)| < 𝛼}, where the angle 𝛼 is measured from the negative real axis. The FLMM is said to be 1. 𝐴(𝛼)-stable if 𝑆(𝛼) ⊆ 𝑆. 2. 𝐴-stable if it is 𝐴(𝜋 − 𝛽𝜋/2)-stable. That is, Σ𝛽 ⊆ 𝑆. 3. unconditionally stable if it is 𝐴(0)-stable. That is if the negative real line (−∞, 0) ⊆ 𝑆. We analyze the stability of the new FLMM through its stability region 𝑆 = {𝛿(𝜉) = (1 − 𝜉)𝛽 𝑝(𝜉): |𝜉| > 1} = ℂ\𝑆 𝑐 , where 𝑆 𝑐 = {𝛿(𝜉) = (1 − 𝜉)𝛽 𝑝(𝜉): |𝜉| ≤ 1} is the unstable region. Theorem 5 The unstable region 𝑆 𝑐 is bounded and symmetric about the real axis. Moreover, For 0 < 𝛽 ≤ 1, if the imaginary part of 𝜉, ℑ(𝜉) > 0, then the real part ℜ(𝛿(𝜉)) > 0 and ℑ(𝛿(𝜉)) < 0. Proof. For the boundedness of 𝑆 𝑐, we see that for |𝜉| ≤ 1, |𝛿(𝜉)| ≤ (1 + |𝜉|)𝛽 [(1 + 𝛽 2 ) + 𝛽 2 |𝜉|] ≤ 2𝛽 (1 + 𝛽) < ∞. For the symmetry about the real axis, we immediately see that 𝛿(𝜉)̅ = 𝛿(𝜉). For 𝜉 = 𝑒 𝑖𝜃 , we have 1 − 𝜉 = (𝑒 − 𝑖𝜃 2 − 𝑒 𝑖𝜃 2 ) 𝑒 𝑖𝜃 2 = 2𝑖sin 𝜃 2 𝑒 𝑖𝜃 2 = 2sin 𝜃 2 𝑒 𝑖( 𝜃 2 − 𝜋 2 ) =: 𝑏𝑒 𝑖𝜙 , where 𝜙 ≡ 𝜙(𝜃) = 𝜃 2 − 𝜋 2 and 𝑏 ≡ 𝑏(𝜃) = 2sin 𝜃 2 = 2cos𝜙 > 0 for 0 < 𝜃 < 2𝜋. Now, writing 𝛿(𝜉) = (1 − 𝜉)𝛽 + 𝛽/2(1 − 𝜉)𝛽+1, we have for the real part of 𝛿(𝜉), ℜ(𝛿(𝜉)) = 𝑏𝛽 [cos𝛽𝜙 + 𝛽cos𝜙cos(𝛽 + 1)𝜙] =: 𝑏𝛽 𝑔(𝜃) (28) where, with some trigonometric manipulations, 𝑔(𝜃) = (1 + 𝛽 2 ) cos𝛽𝜙 + 𝛽 2 cos(𝛽 + 2)𝜙. ANALYSIS OF A NEW TYPE OF SECOND-ORDER FRACTIONAL LINEAR MULTI-STEP METHOD 115 Now, 𝑔′(𝜃) = −𝛽 (1 + 𝛽 2 ) [sin𝛽𝜙 + sin(𝛽 + 2)𝜙]𝜙′ = − 1 2 𝛽(2 + 𝛽)sin(𝛽 + 1)𝜙cos𝜙 > 0, because, for 0 < 𝜃 < 𝜋, 𝜙 ∈ (−𝜋/2,0) where cos𝜙 > 0 and (𝛽 + 1)𝜙 ∈ [− (𝛽+1)𝜋 2 , 0] in the quadrants III and IV for 0 < 𝛽 ≤ 1 where sin(𝛽 + 1)𝜙 < 0. Hence, 𝑔(𝜃) is increasing with 𝑔(0) = cos(𝛽𝜋/2) > 0. Thus, 𝑔(𝜃) > 0 for 0 < 𝜃 < 𝜋. It then follows from the symmetry that ℜ(𝛿(𝜉)) = ℜ(𝛿(𝜉)̅) > 0 . For the imaginary part of 𝛿(𝜉), ℑ(𝛿(𝜉)) = 𝑏𝛽 [sin𝛽𝜙 + 𝛽cos𝜙sin(𝛽 + 1)𝜙] =: 𝑏𝛽 ℎ(𝜃) < 0, (29) because, when 0 < 𝜃 < 𝜋, we see that 𝜙 and 𝛽𝜙 are in the quadrant IV where cos𝜙 > 0 and sin𝛽𝜙 < 0, and (𝛽 + 1)𝜙 is in the quadrants III and IV where sin(𝛽 + 1)𝜙 < 0. This gives ℑ(𝛿(𝜉(𝜃)) < 0 for 0 < 𝜃 < 𝜋. When 𝜋 < 𝜃 < 2𝜋, we again see that 𝜙 and 𝛽𝜙 are in quadrant I where cos𝜙 > 0 and sin𝛽𝜙 > 0, and (𝛽 + 1)𝜙 is in the quadrants I and II where sin(𝛽 + 1)𝜙 > 0. This gives ℑ(𝛿(𝜉(𝜃)) > 0 for 𝜋 < 𝜃 < 2𝜋. and the proof is completed. Theorem 5 tells us that the new FLMM is 𝐴( 𝜋 2 )-stable for 0 < 𝛽 ≤ 1. We have a stronger result. Theorem 6 The FLMM in (22) is 𝐴-stable for 0 < 𝛽 ≤ 1. Proof. From (28) and (29), the tangent at 𝜃 ∈ [0, 𝜋] on the stability region boundary {𝛿(𝜉): |𝜉| = 1} is ℎ(𝜃)/𝑔(𝜃) with its derivative given by 𝑑 𝑑𝜃 ℎ(𝜃) 𝑔(𝜃) = 𝛽(𝛽 + 1)(𝛽 + 2)cos2𝜙 2(𝑔(𝜃))2 > 0. Thus, the tangent is monotonically increasing in [0, 𝜋] with minimum at 𝜃 = 0, (ℎ/𝑔)(0) = −tan(𝛽𝜋/2). 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. Figure 1. Unstable regions and 𝐴-stable tangent boundaries for the new FLMM in (23). 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) ) H.M. NASIR and KHADIJA AL-HASSANI 116 𝑆 𝑐 = {𝛿(𝜉): |𝜉| ≤ 1}. 6. 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). (a) 𝛽 = 0.25 (b) 𝛽 = 0.50 (c) 𝛽 = 0.75 (d) 𝛽 = 0.90 Figure 2. Comparing the unstable regions of FAM1, NFLMM2 and FBDF2 for various 𝛽. 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. We note this from the observation that the unstable regions (see also the figures in Figure 2) satisfy the ordering 𝛿𝐹𝐵𝐷𝐹2(−1) < 𝛿𝑁𝐹𝐿𝑀𝑀2(−1) < 𝛿𝐹𝐴𝑀1(−1) < 𝛿𝐹𝑇2(−1) = +∞. 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. ANALYSIS OF A NEW TYPE OF SECOND-ORDER FRACTIONAL LINEAR MULTI-STEP METHOD 117 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]. 7. 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. Acknowledgment The authors are grateful for the reviewers' constructive comments and suggestions to improve the presentation of the paper and to Dr. Kamel Nafa for the proofreading and corrections. References 1. Diethelm, K. The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type. Springer Science and Business Media, 2010. 2. Brandibur, O., Garrappa, R., and Kaslik, E. Stability of systems of fractional-order differential equations Caputoaputo derivatives. Mathematics 2021, 9(8), 914-933. https://doi.org/10.3390/math9080914. 3. Jia, J., Huang, X., Li, Y., Cao, J., and Alsaedi, A. Global stabilization of fractional-order memristor-based neural networks with time delay. IEEE transactions on neural networks and learning systems, 2019, 31(3), 997-1009. 4. Metzler, R., and Klafter, J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics reports, 2000, 339(1), 1-77. 5. Metzler, R., and Nonnenmacher, T. F. Space-and time-fractional diffusion and wave equations, fractional Fokker-Planck equations, and physical motivation. Chemical Physics, 2002, 284(1-2), 67-90. 6. Povstenko, Y. Z. Fractional heat conduction equation and associated thermal stress. Journal of Thermal Stresses, 2004, 28(1), 83-102. 7. Oldham, K., and Spanier, J. The fractional calculus theory and applications of differentiation and integration to arbitrary order. Elsevier, 1974. 8. Podlubny, I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, vol. 198. Academic Press, 1998. 9. Lubich, C. Fractional linear multistep methods for Abel-Volterra integral equations of the second kind. Mathematics of computation, 1985, 45(172), 463-469. 10. Lubich, C. Discretized fractional calculus. SIAM Journal on Mathematical Analysis, 1986, 17(3), 704-719. 11. Galeone, L., and Garrappa, R. On multistep methods for differential equations of fractional order. Mediterranean Journal of Mathematics, 2006, 3(3), 565-580. 12. Galeone, L., and Garrappa, R. Fractional Adams–Moulton methods. Mathematics and Computers in Simulation, 2008, 79(4), 1358-1367. 13. Galeone, L., and Garrappa, R. Explicit methods for fractional differential equations and their stability properties. Journal of Computational and Applied Mathematics, 2009, 228(2), 548-560. 14. Garrappa, R. On some explicit Adams multistep methods for fractional differential equations. Journal of computational and applied mathematics, 2009, 229(2), 392-399. 15. Nasir, H.M., Gunawardana, B.L.K., and Abeyrathna, H.M.N.P. A second order finite difference approximation for the fractional diffusion equation. International Journal of Applied Physics and Mathematics, 2013, 3(4), 237-243. 16. Gunarathna, W.A., Nasir, H.M., and Daundasekera, W.B. An explicit form for higher order approximations of fractional derivatives. Applied Numerical Mathematics, 2019, 143, 51-60. https://doi.org/10.3390/math9080914 H.M. NASIR and KHADIJA AL-HASSANI 118 17. Nasir, H.M., and Nafa, K. Algebraic construction of a third order difference approximation for fractional derivatives and applications. Australia and New Zealand Industrial and Applied Mathematics Journal, Engineering Mathematics and Application Conference special edition 2017, 2018, 59, C254-C245. 18. Nasir, H.M., and Nafa, K. A new second order approximation for fractional derivatives with applications. SQU Journal of Science, 2018, 23, 1, 43-55. 19. Aceto, L., Magherini, C., and Novati, P. On the construction and properties of m-step methods for FDEs. SIAM Journal on Scientific Computing, 2015, 37(2), A653-A675. 20. Zhao, L., and Deng, W. A series of high-order quasi-compact schemes for space fractional diffusion equations based on the superconvergent approximations for fractional derivatives. Numerical Methods for Partial Differential Equations, 2015, 31(5), 1345-1381. 21. Dimitrov, Y., Miryanov, R., and Todorov, V. Asymptotic expansions and approximations for the Caputo derivative. Computational and Applied Mathematics, 2018, 37, 4, 5476-5499. 22. Garrappa, R. Trapezoidal methods for fractional differential equations: Theoretical and computational aspects. Mathematics and Computers in Simulation, 2015, 110, 96-112. Received 29 April 2022 Accepted 5 September 2022