SQU Journal for Science, 2018, 23(1), 43-55 DOI: http://dx.doi.org/10.24200/squjs.vol23iss1pp43-55 Sultan Qaboos University 43 A New Second Order Approximation for Fractional Derivatives with Applications Haniffa M. Nasir* and Kamel Nafa Department of Mathematics and Statistics, College of Science, Sultan Qaboos University, P.O. Box 36, Al-Khod 123, Sultanate of Oman. *Email: nasirh@squ.edu.om ABSTRACT: We propose a generalized theory to construct higher order Grünwald type approximations for fractional derivatives. We use this generalization to simplify the proofs of orders for existing approximation forms for the fractional derivative. We also construct a set of higher order Grünwald type approximations for fractional derivatives in terms of a general real sequence and its generating function. From this, a second order approximation with shift is shown to be useful in approximating steady state problems and time dependent fractional diffusion problems. Stability and convergence for a Crank-Nicolson type scheme for this second order approximation are analyzed and are supported by numerical results. Keywords: Grünwald approximation, Generating function, Fractional diffusion equation, Steady state fractional equation, Crank-Nicolson scheme, Stability and convergence. تطبيقاتها معمن الدرجة الثانية للمشتقات الكسرية تقريب جديد كمال نافعو *حنيفة محمد ناصر درجاتاللمشتقات الكسرية. نستخدم هذا التعميم لتبسيط براهين ل يةلتقريبالد اجرنومن نوع درجة عالية بناء ةنظريلنقترح في هذا البحث تعميما :صالملخ متتالية استخدامكسرية بالمشتقات لل يةلتقريبالد اجرنونوع بدرجة عالية لمشتقات الكسرية. كما نقوم ببناء مجموعةلالمعروفة ات تقريبال المتعلقة بأنواع من االنتشار الكسريومسائل ثابتة مسائل الحالة اللتقريب من الدرجة الثانية مع االنسحاب مفيد في تقريب حقيقية عامة والدالة المولدة لها. وبهذا تم برهان أن ا . نكلسن لتقريبات من الدرجة الثانية ودعمهما بنتائج عددية-كرانك وتم تحليل االستقرار والتقارب لطريقة من نوع المرتبطة بالزمن. نكلسن، االستقرار والتقارب.-، معادلة االنتشار الكسري، الحالة الثابتة للمعادلة الكسرية، طريقة كرانكالدالة المولدة ،لداجرنوتقريبات :مفتاحيةالكلمات ال 1. Introduction ractional calculus has a history that goes back to L’Hospital, Leibniz and Euler [1,2]. A historical account of early works on fractional calculus can be found, for example, in [3]. Fractional integral and fractional derivative are extensions of the integer order integrals and derivatives to a real or complex order. Various definitions of fractional derivatives have been proposed in the past, among which the Riemann-Liouville, Grünwald-Letnikov and Caputo derivative are common and established. Each definition characterizes certain properties of the integer order derivatives. Recently, fractional calculus found its way into the application domain in science and engineering. The field of application includes, but is not limited to, oscillation phenomena [4], visco-elasticity [5], control theory [6] and transport problems [7]. Fractional derivatives are also found to be suitable to describe anomalous transport in an external field derived from the continuous time random walk [8], resulting in a fractional diffusion equation. The fractional diffusion equation involves fractional derivative either in time, in space or in both variables. Fractional derivative is approximated by the Grünwald approximation obtained from the equivalent Grünwald- Letnikov formula of the Riemann-Liouville fractional derivative. Numerical experience and theoretical justifications have shown that application of this approximation as it is in the space fractional diffusion equation results in unstable solutions when explicit, implicit and even when the Crank-Nicolson (CN) type schemes are used [9]. F H.M. NASIR and K. NAFA 44 The latter two schemes are popular for their unconditional stability for classical diffusion equations. This peculiar phenomenon for the implicit and CN type schemes is corrected and the stability is restored when a shifted form of the Grünwald approximation is used [9,10]. The Grünwald approximation is known to be of first order with the space discretization size h in the shifted and non-shifted form and is, therefore, useful only in first order schemes such as explicit Euler (forward) and implicit Euler (backward) for the fractional diffusion equation. Since the CN approximation scheme is of second order in time step  , Tadjeran et al. [11] used extrapolation improvement for the space discretization to obtain a second order accuracy. Subsequently, second order approximations for the space fractional derivatives were obtained through some manipulations on the Grünwald approximation. Nasir et al. [12] obtained a second order accuracy through a non-integer shift in the Grünwald approximation, displaying super convergence. Convex combinations of various shifts of the shifted Grünwald approximation were used to obtain higher order approximations in Chinese schools [13, 14, 15, 16, 17], some of which are unconditionally stable for the space fractional diffusion equation with CN type schemes. Zhao and Deng [18] extended the concept of super convergence to derive a series of higher order approximations. Earlier, Lubich [19] obtained some higher order approximations for the fractional derivative without shift for orders up to 6. Numerical experiments show that these approximations are also unstable for the fractional diffusion equation when the methods mentioned above are used. Shifted forms of these higher order approximations diminish the order to one, making them unusable as Chen and Deng [20,21] noted. In this paper, we construct a new second order Grünwald type approximation which can be used with shifts without reducing its order. We apply this second order approximation in steady state problems and time dependent fractional diffusion problems. A CN type scheme is devised for this approximation and a justification for stability and convergence is given. The rest of the paper is organized as follows. In Section 2, definitions and notations are introduced. In Section 3, the main results of generalization are presented. In Section 4, approximations of order one and two with shifts are constructed. In Section 5 the constructed second order approximation is applied to devise numerical schemes for steady state and time dependant space fractional diffusion equations. Stability and convergence for the scheme of fractional diffusion equations are analysed in Section 6. Supporting numerical results are presented in Section 7 and conclusions are drawn in Section 8. 2. Definitions and notations Denote by R),( 1 L , the space of Lebesgue integrable functions:    <|)(||=)(1 dxxffL . Let )()( 1 RLxf  and assume that it is sufficiently differentiable so that the following definitions hold. The left and right Riemann-Liouville fractional derivatives of order R are defined respectively as , )( )( )( 1 =)( 1       d x f dx d n xfD n x n n x RL     (1) and , )( )( )( 1)( =)( 1       d x f dx d n xfD nxn nn RL x       (2) where =n is an integer with nn <1  and )( denotes the Gamma function. The corresponding left and right Grünwald-Letnikov (GL) fractional derivatives are given respectively by ),(1)( 1 lim=)( 0= 0 khxf kh xfD k kh x GL              (3) ),(1)( 1 lim=)( 0= 0 khxf kh xfD k kh GL x              (4) where !)1( 1)( = kkk           . It is known that the Riemann-Liouville and Grünwald-Letnikov definitions are equivalent [22]. Hence, from now onwards, we denote the left and right fractional derivatives as )(xfDx   and )(xfDx   respectively. A NEW SECOND ORDER APPROXIMATION 45 The Fourier transform (FT) and inverse FT of an integrable function )()( 1 Lxf are given by dxexffxf xi     )(=)( ˆ=:)))(((F and )(=)(ˆ=)))((ˆ( 1 xfdefxf xi       F respectively. The FT is linear: ).(ˆ)(ˆ=)))(()((  gfxgxf F For a function f at a point R,  x , the FT is given by ).(ˆ=)))(((   fexf xi F If )(),( xfDxf x   and )()( 1  LxfDx  , the FTs of the left and right fractional derivatives are given by )(ˆ)(=)))(((   fixfD x F and )(ˆ)(=)))(((   fixfD x   F respectively [22]. For a fixed h , the Grünwald approximations for the fractional derivatives in (1) and (2) are obtained by simply dropping the limit in the GL definitions (3) and (4) as ( ) ( ) , , =0 =0 and 1 1 ( ) = ( ) ( ) = ( ) h k h k k k f x g f x kh f x g f x kh h h                respectively, where ( ) = ( 1) k k k g         . When )(xf is defined in the intervals ],[ ba , it is zero-extended outside the interval to adopt the definitions of fractional derivatives and their approximations. The sums are restricted to a finite up to N which grows to infinity as 0h . Often, N is chosen to be = x a N h      and = b x N h      for the left and right fractional derivatives respectively, to cover the sum up to the boundary of these domain intervals, where ][ y is the integer part of y . The left and right fractional derivatives, in this case, are denoted by )(xfDxa  and )(xfDbx  respectively. The Grünwald approximations are of first order accuracy and display unstable solutions in the approximation of fractional diffusion equation by implicit and CN type schemes [9]. As a remedy, a shifted form of Grünwald formula with shift r is used: ,>),)(( 1 =)( )( 0= , axhrkxfg h xf k rN k rh        (5) ,<),)(( 1 =)( )( 0= , bxhrkxfg h xf k rN k rh        (6) where the upper limits of the summations have been adjusted to cover the shift r . Meerchaert et al. [9] showed that for a shift 1=r , the shifted approximations )(1, xfh    are also first order approximations with unconditional stability restored in implicit Euler and CN type schemes for space fractional diffusion equations. As for higher order approximations, Nasir et al.[12] derived a second order approximation with a non-integer shift, /2= r , displaying super convergence. ).()(=)( 2 /2, hOxfDxf xah     (7) Tian et al. [13] used convex combinations of different shifted Grünwald forms to obtain two second order approximations: )()(=)()( 2 ,2,1 hOxfDxfxf xaqhph     (8) with )2( 2 = 1 qp q    and )2( 2 = 2 qp p    for 1)(1,(1,0),=),( qp . Hao et al. [14] obtained an quasi-compact order 4 approximation. All the above approximations are based on the shifted Grünwald approximations rh , . The weights )( k g are the coefficients of the Taylor series expansion of the power function  )(1 z . H.M. NASIR and K. NAFA 46 Higher order approximations were also obtained earlier by Lubich [19], establishing a connection with the characteristic polynomials of multistep methods for ordinary differential equations. Specifically, if )(),( zz  are the characteristic polynomial [23] of a multistep method of order p , then          )(1/ )(1/ z z gives the weights for the Grünwald type approximation of the same order for the fractional derivative of order  . From the backward multistep methods, Lubich [19] derived higher order approximations of up to order six in the form   ))(1 1 (=)( 1= j j z j zW  . The generating functions )(zW given in Table 1 are of order  for 61  . Table 1. Lubich approximation generating functions.    zzW 1=)( 1         2 2 2 1 2 2 3 =)( zzzW         32 3 3 1 2 3 3 6 11 =)( zzzzW         432 4 4 1 3 4 34 12 25 =)( zzzzzW         5432 5 5 1 4 5 3 10 55 60 137 =)( zzzzzzW         65432 6 6 1 5 6 4 15 3 20 2 15 6 20 49 =)( zzzzzzzW 3. Generalization of Grünwald approximation In this section, we generalize the concept of shifted Grunwald approximation to an arbitrary sequence of weights and analyse its properties. This generalization is then used to construct higher order approximations for fractional derivatives. For a sufficiently smooth function )(xf , denote the left Grünwald type operator with shift r and weights rkw , as ).)(( 1 =)( , 0= , hrkxfw h xf rk k rh       (9) Definition 1 A sequence  rkw , of real numbers is said to approximate the fractional derivative )(xfDx   at x with shift r in the sense of Grünwald if ).(lim=)( , 0 xfxfD rh h x      Definition 2 A sequence  rkw , of real numbers is said to approximate the fractional derivative )(xfDx   at x with shift r and order 1p if ).()(=)( , p rhx hOxfxfD    (10) We consider the generating function k rkk zwzW ,0= =)(   of the weights rkw , which will play a central role in constructing approximations. Remark 1 1. Analogous definitions hold for the right fractional derivative with the same weights rkw , . We denote the right Grünwald approximation as )(, xfrh    . A NEW SECOND ORDER APPROXIMATION 47 2. The property (10) is the consistency condition for the Grünwald type approximation operator  rh   , to the fractional differential operator. This consistency condition tells that the order of approximation is at least one. 3. For convenience, we do not distinguish the sequence of weights rkw , and its generating function )( zW in the above definitions, and in the discussions and results that follow. The following theorem establishes an equivalent characterization of the generator )( zW for an approximation of fractional differential operator with order 1p and shift r . Theorem 1 Let mnn ,<1   be a non-negative integer, )()( 1 R   nm Cxf and )()( 1 RLxfD k x   for 10  nmk . Then, the generating function )( zW of a real sequence }{ ,rkw approximates the left fractional differential operator for )(xf with order p and shift r , mp 1 , if and only if ).(1=)( 1 :=)( przz r zOeeW z zG    (11) Moreover, if l lplr zrazG )(1=)( =   , we have for the left fractional derivative ).()()()(=)(, mp xp p xrh hOxfDrahxfDxf     (12) Proof. Taking the FT of )(, xfrh    in (9), with the help of linearity, we obtain )))()((( 1 =)))((( , 0= ,    hrkxfw h xf rk k rh     FF )(ˆ 1 = )( , 0=    few h ihrk rk k    )(ˆ)( )( = , 0=      fiew ih e kih rk k rih    )(ˆ)()(= , 0=    fiew z e kz rk k rz          )(ˆ))((=)(ˆ))((=    fizGfieW z e r z rz  ),(ˆ)()(=)(ˆ)()(= 0=0=   fihrafizra ll l l l l l    where we have used hiz = . Taking the inverse FT, we have ).()()(=)( 1 0= , mll xal m l rh hOhxfDraxf       (13) Now, (11) holds if and only if 0=)(1,=)( 0 rara l , for 1,1,2,= pl  . Equation (12) holds immediately. □ Remark 2 The equivalent condition (11) for order p approximation in Theorem 1 holds for the right fractional derivative as well. In fact, the condition on the generating function )( zW for the right fractional derivative is ))((1=)( )( 1 :=)( przz r zOeeW z zG      which is equivalent to (11). One of the consequences of Theorem 1 is the following consistency condition. Corollary 1 If the generating function )( zW gives a consistent approximation of the left and right fractional differential operators, then , =0 and hence(1) = 0, , = 0. k r k W w   Proof. By Remark 1 and 2, the order p is at least one. When 0z , the condition (11) becomes ))((1=)( przz zOzeeW    . Take the limit as 0z . H.M. NASIR and K. NAFA 48 Using Theorem 1, one can check algebraically the following propositions for the generating function of approximation operators known previously in [12, 13, 19, 20]. Proposition 1 The approximation given by (7) is of second order. Proof. The coefficients )( k g have the generating function )(1 z . Since the shift of the approximation is = / 2 ,r  it is enough to check the function   )(1 1 =)( /2 zz ee z zG   . Taylor series expansion gives that )( 24 1=)( 42 2 zOzzG   which confirms the second order. □ Proposition 2 The approximation given by (8) is of second order accuracy. Proof. Since ph , and qh , have the same generating function  )(1 z with shifts p and q respectively, check the Taylor series of      )(1)(1=)( 21 zqzzpz ee z ee z zG   which gives )( 224448 1=)( 32 2 zOz pqqp zG           . □ Proposition 3 The Lubich generating functions ,6,1,2,=),( pzWp in Table 1 are of order p accurate without shift. Moreover, if a non-zero shift r is introduced to the approximation, the orders reduce to one. Proof. When there is no shift ( 0=r ), Taylor expansion gives )(1=)/( pz p zOzeW    . When 0r , the order is determined by ),(1=))())(1((1=)/( zOzOzOzeWe pz p rz    for 61  p reducing the orders to 1. □ 4. Construction of higher order approximation We construct higher order approximation generating functions for fractional derivatives with shifts by the use of Theorem 1. The importance of Theorem 1 is that the construction process is entirely confined to algebraic manipulation with the aid of Taylor series expansion. In this paper, we choose the form   )(=)( 2 210 p p zzzzW   (14) for the generators. Consider the Taylor series expansion .)(=)( 1 =)( 0= l l l rzz r zraeeW z zG     (15) For an order p approximation, we set the first p coefficients )(ra l to satisfy (11) in Theorem 1. That is, we impose conditions 0=)(1,=)( 0 rara k for = 1, 2, , 1k p  . First and second order shifted Grünwald approximations have been constructed by the above algebraic method and we have the following. Theorem 2 The first and second order shifted Grünwald approximations with shift r are given respectively by  )(1=)( 1, zzW r  and   . 2 12 2 2 3 ==)( 22 2102,                                z p r z r p r zzzW r (16) The generating function )(1, zW r corresponds to the shifted Grunwald formula )(, xfrh  in (5) and (6) for the left and right fractional derivatives respectively, and is independent of the shift r . )( 2, zW r gives the second order Grünwald type approximations given by (9) with coefficients rkw , obtained from the Taylor series expansion of )(2, zW r . Note that when there is no shift )(2,0 zW gives the Lubich approximation )(2 zW in Table 1. We denote the Grünwald approximation of order p by , 1, 2,( ), =p ru x p for left and right fractional derivatives respectively. A NEW SECOND ORDER APPROXIMATION 49 The Grünwald weights rkw , can be computed by a recurrence formula [24,25]. Lemma 1 Let j jj zzW   0= =)( , where Rj , 0> and m mm zwzW   0= =))((  . Then, the coefficients m w satisfy the recurrence form 1.,)1)(( 1 =,= 1=0 00   mwmj m ww jjm m j m     (17) For the generating function )(2, zW r , being the power of a polynomial of degree 2 , the upper limit of the summation in (17) goes up to only 2=,2)(min m for 2m . 5. Applications to fractional differential equations In this section, we apply the fourth order Grünwald type approximation derived in the previous section to fractional differential equations. 5.1 Steady state problems Consider the steady state problem ,),(=)( bxaxfxuDxa   (18) .=)(,=)( 21  buau Consider a uniform partition bxxxxa N =<<<<= 210  of the problem domain ],[ ba with subinterval size 1,0,1,2,=,= 1   Nixxh ii  . Problem (18) at i x is approximated by ,,0,1,2,=),(= 2 2, NihOfu iir    where )(= ii xuu and Nixff ii ,0,1,2,=),(=  . In order to keep the discrete values i u within the computational domain, we choose the shift 1=r . Neglecting the )( 2 hO remainder term, we get the approximation scheme ,0,1,2,=,=ˆ12, Nifu ii     (19) where i û are the solutions of (19) and hence the approximation of the exact solution i u for 1,1,2,= Ni  . Let T NN uuuuuU ],ˆ,,ˆ,ˆ,[= 1210   , T N ffffF ],,,,[= 210  with the boundary conditions incorporated in U as )(= 10 au  and )(= 2 bu N  . Then, the matrix formulation of (19) is given by ,=2,1 FUA where 2,1A is an 1)(1)(  NN Toeplitz matrix given by     ,0, 1, =),( 1,1 2,1 elsewhere jiw jiA ji where 1,1 jiw are the coefficients of the Taylor series expansion of the generating function )(2,1 zW . After imposing boundary conditions, we obtain the ready-to-solve second order scheme as ,ˆ=ˆˆ 002,1 NN uAuAFUA  where 2,1  is the reduced matrix of size 1)(1)(  NN obtained from 2,1A by deleting the first and last rows and columns, FU ˆ,ˆ are obtained from FU , respectively by deleting their first and last boundary entries. N AA , 0 are the first and last column vectors of the matrix 2,1A reduced at both ends as above. 5.2 Approximation of fractional diffusion equation We consider the numerical approximation of the space fractional diffusion equation defined in the domain ][0,],[ Tba  : ),,(),(),(= ),( 21 txftxuDKtxuDK t txu bxxa     (20) with the initial and boundary conditions H.M. NASIR and K. NAFA 50 ],[),(=,0)( 0 baxxsxu  ],[0,),(=),(),(=),( 21 Ttttbuttau  where ).( txu is the unknown function to be determined; 21 , KK are non-negative constant diffusion coefficients with 0 21  KK , i.e, not both are simultaneously zero, and ),( txf is a known source term. The boundary conditions are set as follows: If 0 1 K , then 0)( 1 t and if 0 2 K , then 0)( 2 t . We assume that the diffusion problem has a unique solution. The space domain ],[ ba is partitioned into a uniform mesh of size N with subintervals of length Nabh )/(=  , and the time domain ][0, T into a uniform partition of size M with subintervals of length MT/= . These two partitions form a uniform partition of the 2-D domain ][0,],[ Tba  with grid points ),( mi tx , where ihax i = and MmNimt m  0,0,=  . We use the following notations for conciseness: ),(= mi m i txuu , )( 2 1 = 11/2 mmm ttt   and ),(= 1/2 1/2   mi m i txff . We present the CN type scheme with the order 2 approximation in space using 2,1 . Using the approximations  12,212,12,1 ='   KK of order 2, the Crank-Nicolson type scheme at 1,0),0,(  MmNitx mi is ).()( 2 1 '= 221/21 2,1 1 hOfuu uu m i m i m i m i m i       (21) Let m U be the solution of (21) after neglecting the )( 22 hO  terms with Tm N m N mmmm uuuuuU ]ˆ,ˆ,,ˆ,ˆ,ˆ[= 1210   , where m i û becomes the approximation of the exact values m i u with mm uu 00 =ˆ and m N m N uu =ˆ . Then, (21) becomes 1.0,0,)ˆˆ(' 2 =)ˆˆ( 1/21 2,1 1   MmNifuuuu m i m i m i m i m i   Thus, the Crank-Nicolson type scheme in matrix form reads 1,,0)(=)( 1/211   MmFUUBUU mmmmm   (22) where Tm N mmm fffF ],,,[= 1/2 1 1/2 1 1/2 0 1/2     . The matrix B corresponding to the operator '2,1 is given by ).( 2 = 2,122,11 T AKAKB    Re-arranging for 1m U and m U , we have 1.,0,1,2,=,)(=)( 1/21   MmFUBIUBI mmm   (23) Let B̂ be the reduced matrix from B and 1/2ˆ mF be the reduced vector from 1/2m F as was in Section 5.1. After imposing the boundary conditions, equation (23) reduces to the ready-to-solve form 1,,0,1,2,=,ˆˆˆ)ˆ(=ˆ)ˆ( 1/21   MmbFUBIUBI mmmm  where )()(=ˆ 1 0 1 00 m N m NN mmm uuBuuBb   and N BB , 0 are the first( th 0 ) and last( th N ) column vectors of the matrix B reduced again as before, and I is the unit matrix of appropriate size. 6. Stability and convergence We establish the stability of the Crank-Nicolson type scheme (23). We closely follow the analysis in [14] and present some required results. Let 0}==,),,,,(=|{= 010 NiNh vvvvvvvvV R be the space of grid functions in the computational domain in the space interval ],[ ba . For any h Vvu , , define the discrete inner product and its corresponding norm respectively as A NEW SECOND ORDER APPROXIMATION 51 .),(=||a=),( 1 1= uuundvuhvu ii N i ||  Then, we have the following: Lemma 2 Let   0= }{ kk t be a double sided real sequence such that (i) 0 kk tt for 0,k (ii) 0 =   j N Nj t for 0N . Then, the Toeplitz matrix ][= jiN tT  of size 1N is negative definite for 0N . Proof. For 0=N , 0 0 t . This matrix of size 1 is negative definite. For any positive integer N and for any non- zero 1)( N -dimensional vector T N vvvvv ],,,,[= 210  , consider the quadratic form jiij N j N iN T vvtvTv  0=0== . Summing the terms diagonally, we have jkj kN j kk N k j N j jkj kN j k N Nk N T vvttvtvvtvTv        1 0=1= 2 0= 0 1 0== )(== )|||(1/2)(|)(|| 22 1 0=1= 2 0 jkj kN j kk N k vvttvt     || 0||=||)(|| 2 = 2 1= 2 0     |||||| vtvttvt k N Nk kk N k which completes the proof. □ Lemma 3 The matrices 2,1  and T A 2,1 ˆ and hence their corresponding operators  12 ,  and  12 ,  respectively are negative definite for 2.1  Proof. The matrix 2,1  of size 1)(1)(  NN has the form . 1 =ˆ 1,12,13,12,11,1 0,11,14,13,12,1 0,11,12,13,1 0,11,12,1 0,11,1 2,1                       wwwww wwwww wwww www ww h A NNN NNN    By virtue of Lemma 2, it is enough to prove that the coefficients ,1kw of the generator k kk zwzW ,10=2,1 =)(   satisfy the following properties for 21  . 1 0,1 0,1 2,1 ,1 ,1 =0 0, 0, = 0 for 3 and 0 for all 2. N m m m k k w w w t t w m w N           Let     )(= 1 2 12 2 1 2 3 =)( 2 210 2 2,1 zzzzzW                          . Then, 0 1 2 3 = 0    , 0 2 2= 1    and 0 1 2 1 = 2    for 21  . Now, 0= 00,1   w and 0= 1 1 01,1     w . Also,   20 2 1 2 002,10,1 21)( 2 1 =     ww  )(21)( 2 1 = 020 2 1 2 0     H.M. NASIR and K. NAFA 52 0, 2 2)1)(( 21)( 2 1 = 0 2 1 2 0                  0.62)(1 6 1 = 20 2 11 3 03,1     w We show that 0,1 mw for all 4m inductively.    0,122)12(3)2)((1 4! 1 = 2 2 2 02 2 10 4 1 4 04,1     w   0.60)204)3)((( 5! 2)1)(( = 2 2 2 02 2 10 4 1 1 5 0 5,1       w Assuming that k w is non negative for 13  mk , we have for 6m , 0))2(2)1(( 1 = 22,111,1 0 ,1     mmm wmwm m w for 21  , since 2,11,1,  mm ww , 1 ),2(2),1(  mm  and 2  are all non-positive for 6m . From Corollary 1, we have the consistency condition 0= 0= kk w  in general which is true for )(2,1 zW as well. Since 0,1=   kNk w for 3N , the last inequality follows from the consistency condition. □ Lemma 4 The approximation operator '2,1 is negative definite. Proof. For any h Vv  , since the diffusion coefficients 21 , KK are non-negative, we have 0),(),(=),'( 12,212,12,1   vvKvvKvv  . □ Theorem 3 If ),,,(= 121 m N mmm vvvv   be the solution of the problem 101,1,=' 1/2 2,1 1/2   MmNiSvv m i m i m it  (24) 0=0,=0 m M m vv .0),(= 0 0 Nixvv ii  Then, .|||||| 1 0= 0          |||||| l m l m Svv  Proof. Taking the inner product of (24) with 1/2m v , we have ).,(=),(),( 1/21/21/2 2 1/21/2   m i m i m i m i m i m it vSvvvv Since 0),'( 1/21/2 2,1   m i m i vv from Lemma 4, we get ).,(),( 1/21/21/2   m i m i m i m i vSvv   Since )( 1 = 11/2 m i mm vvv     and )( 2 1 = 11/2 mmm vvv   , we have         )( 2 1 ),( 1 =),( 111/21/2 mmmmmm vvvvvv    ),(),( 2 1 = 11 mmmm vvvv    ||||),()||||( 2 1 = 1/21/2221 ||||||||   mmmmmm vSvSvv   .|||||| 2 1 1 |||||| mmm vvS   The inequality in the last two lines reduces to A NEW SECOND ORDER APPROXIMATION 53 1.0,|||||| 1   MmSvv mmm ||||||  Summing this for the first m inequalities, we have .0,|||||| 1 0= 0 MmSvv l m l m    ||||||  From the above estimates, we have the following stability result. □ Theorem 4. The CN type difference scheme (23) of order 2 is unconditionally stable for 21  . For the convergence of the approximate solution from the CN type scheme, we have the following. Theorem 5. The approximate solutions of the CN type scheme (23) with the given initial and boundary conditions are convergent for 21  . Proof. Let mmm Uue ˆ=  be the error vector of the exact and approximate solutions mm Uu ˆ, of the diffusion problem (20) respectively. Then the error of the internal grid values m ê satisfy the system 1,01,1,=' 1/2 2,1 1/2   MmNiRee m i m i m it  0,=0,= 0 m M m ee 1,10,= 0  Nie i where m i R are the remainder terms in (21) with )(|| 22 hcR m i  || for some constatnt 0>c . Theorem 3 gives the estimate ).(|||| 22 1 0= hNcRe l m l m     |||| The convergence is then established as 0, h . □ 7. Numerical results We test the approximation scheme devised in Subsection 5.1 using the steady state test problem 1,0, )1( 1)(10 =)( 0     xx n n xuD n x   10=(1)0,=(0) uu with the exact solution n xxu 10=)( . We set 8=n and test for values of the parameter  with 1.1,1.5= and 1.9 . The number of grid subintervals N corresponding to the discretization size = (1 0) /h N was considered for values = 16, 32, ,1024N . The maximum error   PP Uu and the computed convergence orders are listed in Table 2. Table 2. Second order approximation with shift 1=r using )(2,1 zW . = 1.1 = 1.5 = 1.9 N   |||| Uu Order   |||| Uu Order   |||| Uu Order 16 4.8893e-01 —- 2.5141e-01 —- 1.3365e-01 —- 32 1.1592e-01 2.08 6.4851e-02 1.95 3.3951e-02 1.98 64 2.7227e-02 2.09 1.6450e-02 1.98 8.5491e-03 1.99 128 6.3685e-03 2.10 4.1396e-03 1.99 2.1446e-03 2.00 256 1.4873e-03 2.10 1.0383e-03 2.00 5.3703e-04 2.00 512 3.5020e-04 2.09 2.5997e-04 2.00 1.3437e-04 2.00 1024 8.7574e-05 2.00 6.5044e-05 2.00 3.3606e-05 2.00 This test confirms the theoretical justification of second order for the approximation )(2,1 zW . We consider the following test example for the fractional diffusion problem (20). H.M. NASIR and K. NAFA 54 Let ))(1( )1( 1)( =),,(        mm xx m m mxG and 55 0 )(1=)( xxxs  . 1 2 Diffusion coefficients: = 1, = 1.K K 0 Source function: ( , ) = ( ( ) ( , 5, ) 5 ( , 6, ) t f x t e s x G x G x      )).,10,(),9,(5),8,(10),7,(10  xGxGxGxG  0 Initial condition: ( , 0) = ( ).u x s x Boundary conditions: (0, ) = 0, (1, ) = 0.u t u t 0 Exact solution: ( , ) = ( ) . t u x t s x e  The test problem was applied to the CN type numerical scheme developed in Subsection 5.2. All computations were performed using Python Language with Scipy libraries [26] on an i7 notebook computer with 2.7Ghz speed and 12Gb memory and Windows operating system. The second order approximation )(2,1 zW was tested for = 1.1,1.5 and 1.9 . Table 3 lists the maximum error and the order of convergence for grid sizes = = 16, 32, ,1024N M . The partition subinterval sizes are then = 1/ M and = 1/h N for time and space, respectively. Table 3. Second order of convergence for CN type scheme with )(2,1 zW . 1.1= 1.5= 1.9= hN 1/=   |||| Uu Order   |||| Uu Order   |||| Uu Order 16 1.0544e-05 — 9.0719e-06 — 5.6905e-06 — 32 2.8172e-06 1.90 2.3208e-06 1.97 1.4309e-06 1.99 64 7.3008e-07 1.95 5.8863e-07 1.98 3.5731e-07 2.00 128 1.8606e-07 1.97 1.4836e-07 1.99 8.9332e-08 2.00 256 4.6984e-08 1.99 3.7252e-08 1.99 2.2338e-08 2.00 512 1.1806e-08 1.99 9.3341e-09 2.00 5.5852e-09 2.00 1024 2.9592e-09 2.00 2.3362e-09 2.00 1.3964e-09 2.00 The test results show that the second order approximation with )(2,1 zW is justified for its order of convergence and unconditional stability with the CN type scheme. 8. Conclusion A new second order Grünwald type approximation for fractional derivative with shift is constructed algebraically from a generalization of the Grünwald approximation for the left and right fractional derivatives in terms of generating functions. Numerical schemes for this approximation to solve steady state problems and time dependent fractional diffusion problems were devised with proof of stability and convergence. The approach of generating functions could be a useful tool for constructing difference approximation formulas for fractional derivatives. Acknowledgement The authors are grateful to the reviewers for their useful comments and suggestions. This research is supported by Sultan Qaboos University Internal Grant No. IG/SCI/DOMAS/16/10. References 1. Leibniz, G.W. "Letter from Hanover, Germany to G.F.A. L'Hospital, September 30, 1695", Mathematische Schriften, reprinted 1962, Hildesheim, Germany (Olns Verlag) 2, 301-302. 2. Euler, L. “On Transcendental Progressions. That is, Those Whose General Terms Cannot be Given Algebraically”, Commentarii academiae scientiarum Petropolitanae, 1738, 5, 36-57. 3. English Translation by S. Langton at:http://home.sandiego.edu/~langton/eg.pdf 4. Ross, B. The development of fractional calculus 1695–1900. Historia Mathematica, 1977, 4(1), 75–89. A NEW SECOND ORDER APPROXIMATION 55 5. Mainardi, F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos, Solitons and Fractals, 1996, 7(9), 1461–1477. 6. Bagley, R.L., and Torvik, P. A theoretical basis for the application of fractional calculus to viscoelasticity. Journal of Rheology, 1983, 27(3), 201–210. 7. Vinagre, B., Podlubny, I., Hernandez, A., and Feliu, V. Some approximations of fractional order operators used in control theory and applications. Fractional Calculus and Applied Analysis, 2000, 3, 231–248. 8. Metzler, R., and Klafter, J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. Journal of Physics A: Mathematical and General, 2004, 37(31), R161-R204. 9. Barkai, E., Metzler, R., and Klafter, J. From continuous time random walks to the fractional Fokker-Planck equation. Physical Review E, 61, 2000, 1, 132. 10. Meerschaert, M.M., and Tadjeran, C. Finite difference approximations for fractional advection–dispersion flow equations. Journal of Computational and Applied Mathematics, 2004, 172(1), 65–77. 11. Meerschaert, M.M., and Tadjeran, C. Finite difference approximations for two-sided space-fractional partial differential equations. Applied Numerical Mathematics, 2006, 56(1), 80–90. 12. Tadjeran, C., Meerschaert, M.M., and Scheffler, H.-P. A second-order accurate numerical approximation for the fractional diffusion equation. Journal of Computational Physics, 2006, 213(1), 205–213. 13. Nasir, H.M., Gunawardana, B.L.K., and Aberathna, 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. 14. Tian, W., Zhou, H., and Deng, W. A class of second order difference approximations for solving space fractional diffusion equations. Mathematics of Computation, 2015, 84(294), 1703–1727. 15. Hao, Z.P., Sun, Z.Z., and Cao, W.R. A fourth-order approximation of fractional derivatives with its applications. Journal of Computational Physics, 2015, 281, 787–805. 16. Zhou, H., Tian, W., and Deng, W. Quasi-compact finite difference schemes for space fractional diffusion equations. Journal of Scientific Computing, 2013, 56(1) 45–66. 17. Li, C., and Deng, W. A new family of difference schemes for space fractional advection diffusion equation. Advances in Applied Mathematics and Mechanics, 2017, 9(2), 282–306. 18. Yu, Y., Deng, W., Wu, Y., and Wu, J. Third order difference schemes (without using points outside of the domain) for one sided space tempered fractional partial differential equations. Applied Numerical Mathematics, 2017, 112, 126–145. 19. 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. 20. Lubich, C. Discretized fractional calculus. SIAM Journal on Mathematical Analysis, 1986, 17(3), 704–719. 21. Chen, M., and Deng, W. Fourth order difference approximations for space Riemann-Liouville derivatives based on weighted and shifted Lubich difference operators. Communications in Computational Physics 2014,16(02), 516–540. 22. Chen, M., and Deng, W. Fourth order accurate scheme for the space fractional diffusion equations. SIAM Journal on Numerical Analysis, 2014, 52(3), 1418–1438. 23. 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. 24. Henrici, P. Discrete variable methods in ordinary differential equations, SIAM Rev., 1962, 4(3), 261-262. 25. Rall, L.B. Automatic differentiation: Techniques and Applications, Lecture Notes in Computer Science, Springer- Verlag, Vol. 120, 1981. 26. Weilbeer, M. Efficient Numerical Methods for Fractional Differential Equations and Their Analytical Background. Papierflieger, 2005. 27. Walt, S.V.D., Colbert, S.C., and Varoquaux, G. The numpy array: a structure for efficient numerical computation. Computing in Science and Engineering, 2011, 13(2), 22–30. Received 25 May 2017 Accepted 11 September 2017