An explicit form for higher order approximations of fractional derivatives

Applied Numerical Mathematics 143 (2019) 51–60

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

An explicit form for higher order approximations of fractional derivatives ✩ W.A. Gunarathna a , H.M. Nasir b,∗ , W.B. Daundasekera c a b c

Department of Mathematics, General Sir John Kotelawala Defence University, Sri Lanka FracDiff Research Group, Department of Mathematics, Sultan Qaboos University, P. O. Box: 36, Al-Khoudh 123, Muscat, Oman Department of Mathematics, University of Peradeniya, Peradeniya, Sri Lanka

a r t i c l e

i n f o

Article history: Received 19 September 2018 Received in revised form 2 March 2019 Accepted 28 March 2019 Available online 2 April 2019 Keywords: Fractional derivative Grünwald approximation Generating function Vandermonde system Finite difference formula

a b s t r a c t An explicit form for coefficients of shifted Grünwald type approximations of fractional derivatives is presented. This form directly gives approximations for any order of accuracy with any desired shift leading to efficient and automatic computations of the coefficients. To achieve this, we consider generating functions in the form of power of a polynomial. Then, an equivalent characterization for consistency and order of accuracy established on a general generating function is used to form a linear system of equations with Vandermonde matrix. This linear system is solved for the coefficients of the polynomial in the generating function. These generating functions completely characterize Grünwald type approximations with shifts and order of accuracy. Incidentally, the constructed generating functions happen to be a generalization of the previously known Lubich forms of generating functions without shift. We also present a formula to compute leading and some successive error terms from the coefficients. We further show that finite difference formulas for integer-order derivatives with desired shift and order of accuracy are some special cases of our explicit form. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Fractional derivative is usually approximated by the standard Grünwald approximation obtained from the GrünwaldLetnikov formula by simply dropping the limit in its definition. It is completely identified by the generating function W 1 ( z) = (1 − z)α whose Taylor series expansion provides coefficients for the infinite sum in the Grünwald approximation formula [17]. The Grünwald approximation is of first order accuracy with respect to the interval domain discretization step size h. Meerchaert et al. [13] introduced a shifted form for the Grünwald approximation to remedy some stability difficulties that occur when the unshifted Grünwald formula is used to approximate space fractional diffusion equations (SFDEs). This shifted Grünwald approximation is also shown to be of first order accuracy and is, therefore, useful only in first order approximation methods for differential equations involving space-fractional derivatives.

✩

This research was supported by Sultan Qaboos University Internal Grant No. IG/SCI/DOMAS/16/10. Corresponding author. E-mail addresses: [email protected] (W.A. Gunarathna), [email protected], [email protected] (H.M. Nasir), [email protected] (W.B. Daundasekera).

*

https://doi.org/10.1016/j.apnum.2019.03.017 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

52

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

Since then, higher order approximations for fractional derivatives were sought for space-fractional diffusion type equations. Nasir et al. [16] obtained a second order accuracy through a non-integer shift in the shifted Grünwald approximation, displaying super convergence. Convex combinations of various shifts of the shifted Grünwald approximations were used to obtain higher order approximations in [20,7,23,11,22,2]. All of the above approximations are based on the first order shifted Grünwald approximations with the generating function W 1 ( z). Earlier, Lubich [12] obtained some higher order approximations for the fractional derivatives without shift for orders up to 6 targeting the time-fractional initial value problems. However, numerical experiments show that using these higher order approximations in space-fractional differential equations displays the same stability difficulties as appeared in the first order case. Shifted forms of these higher order approximations reduce the orders to one, making them uninteresting as Chen and Deng [3,4] noted. In search of higher order shifted Grünwald type approximations, Nasir and Nafa [14,15] have recently constructed generating functions that retain higher order accuracy for Grünwald type formula with shift. They generalized the concept of shifted Grünwald formula for approximation by identifying it with coefficient generating functions. The generating functions for each order of approximation have to be constructed separately by using Taylor series expansions. This process is almost impractical for hand computations and time consuming even for computer applications with symbolic computations. It is, therefore, highly desirable to have formulations, that are easily and efficiently computable with some level of automation. In this paper, we present an explicit form for approximating generating functions for the coefficients of shifted Grünwald type approximations for any order with an arbitrary shift. The non-shift versions of these generating functions reduce to the Lubich generating functions and hence our formulation can be regarded as a generalization with shift of the Lubich approximations. Moreover, our formulation also gives a unified general form for finite difference approximations for any integer-order derivative. The rest of the paper is organized as follows. In Section 2, some preliminaries of earlier works are introduced. In Section 3, the explicit form for generating functions is established. Section 4 demonstrates that the explicit formula can be used to obtain finite difference approximations for integer-order derivatives for any desired order and shift. A brief discussion on stability for these approximations is given in Section 5. Conclusions are drawn in Section 6. 2. Prelimineries, definitions and notations The space of Lebesgue integrable functions defined in  ⊆ R is denoted by L 1 (). Let f (x) ∈ L 1 (R) and assume that it is sufficiently differentiable so that the following definitions hold. The left and right Riemann-Liouville fractional derivatives of real order α > 0 are defined respectively as

RL

Dα x− f (x) =

1

dn

(n − α ) dxn

x −∞

f (η) dη , (x − η)α +1−n

(1)

and RL

Dα x+ f (x) =

(−1)n dn (n − α ) dxn

∞ x

f (η) dη , (η − x)α +1−n

(2)

where n is an integer satisfying n − 1 < α ≤ n and (·) denotes the gamma function. The corresponding left/right (∓) Grünwald-Letnikov (GL) definition is given respectively by

 

GL

where

∞ 1  α Dα (−1)k x∓ f (x) = lim α k h→0 h

α  k

f (x ∓ kh),

(3)

k =0

=

(α +1) (α +1−k)k! and h is the subinterval size of a uniform discretization of the interval of interest.

It is known that, under sufficient smoothness of f with vanishing constant (namely, f (0) = 01 ), the Riemann-Liouville and Grünwald-Letnikov definitions are equivalent (see [17] and Theorem 2.25 in [5]). Hence, from now onwards, we denote α the left and right fractional differential operators as D α x− and D x+ respectively. 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 in (3):

1

This additional condition is stated from the observation that, for a non-zero constant C , we have

GL

Dα x∓ C = 0 while

RL

Dα x∓ C = 0.

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

53

Table 1 Lubich approximating generating functions. W 1 ( z) = (1 − z)α W 2 ( z) = W 3 ( z) = W 4 ( z) = W 5 ( z) = W 6 ( z) =



3

2   

− 2z + 12 z2

α α

11 6

− 3z + 32 z2 − 13 z3

25 12

− 4z + 3z2 − 43 z3 + 14 z4

137 60 49 20

− 5z + 5z2 −

− 6z +

15 2 z 2

−

10 3 z 3 20 3 z 3

α

+ 54 z4 − 15 z5 +

15 4 z 4

α

− 65 z5 + 16 z6

α

∞ 1  (α ) gk f (x ∓ kh),

δxα∓,h f (x) = α h (α )

where gk

(1 − z)α .

= (−1)

k

α 

k =0

k

are the coefficients of the Taylor series expansion of the approximating generating function W 1 ( z) =

When f (x) is defined in the intervals [a, b], it is zero-extended outside the interval to adopt these definitions of fractional derivatives and their approximations.

a  The sums are

xthen  restricted to finite terms up to N which grows to infinity as h → 0. Ideally, N is chosen to be N = x− and N = b− for the left and right fractional derivatives respectively to cover the h h sums up to the boundary of the domain intervals, where [ y ] denotes the integer part of y. The Grünwald approximations are of first order accuracy: α Dα x∓ f (x) = δx∓,h f (x) + O (h )

and display unstable solutions in the approximation of SFDE by implicit and Crank-Nicolson (CN) type schemes [13]. As a remedy, a shifted form of the Grünwald formula with shift r is used: N +r 1  (α ) gk f (x ∓ (k − r )h),

δxα∓,h,r f (x) = α h

k =0

where the upper limit of the summation has been adjusted to cover the shift r. Meerchaert et al. [13] showed that δxα∓,h,r f (x) is also of first order accuracy for any shift, and for r = 1, the unconditional stability is restored in implicit Euler and CN type schemes for SFDEs. As for higher order approximations, Nasir et al. [16] derived a second order approximation by a non-integer shift, displaying super convergence:

δx∓,h,α /2 f (x) = D αx∓ f (x) + O (h2 ). Chen et al. [20] used convex combinations of different shifted Grünwald forms to obtain order 2 approximations:

λ1 δx∓,h, p + λ2 δx∓,h,q = D αx∓ f (x) + O (h2 ) α −2q

α −2p

with λ1 = 2( p −q) and λ2 = 2( p −q) for ( p , q) = (1, 0), (1, −1). We note that all of the above higher order approximations were based on some manipulations on the first order Grünwald approximations with various shifts. Higher order approximations were also obtained earlier by Lubich [12] establishing a connection with the characteristic polynomials of multistep methods for ordinary differential equations. Specifically, if ρ ( z), σ ( z) are, the characteristic polynomials of a multistep method of a given order of convergence [8], then W ( z) =

σ (1/z) α gives the coefficients for ρ (1/z)

the Grünwald type approximation of the same order for the fractional derivative of order α . From the pbackward multistep methods, Lubich [12] derived higher order approximations of orders up to six in the form W p ( z) = ( j =1 1j (1 − z) j )α . The generating functions W p ( z) given in Table 1 are of order p for 1 ≤ p ≤ 6. These higher order approximations also suffer the same instability issues without shift when used in the numerical schemes mentioned above. Shifted forms of these approximations reduce the order to one, making them uninteresting [3,4, 14]. 2.1. Higher order approximations with shift A generalization of the Grünwald approximation was proposed by Nasir and Nafa [14,15] to construct higher order approximating generating functions that can be used with a shift without reducing the order of accuracy. The definitions and results in their work are given below.

54

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

Consider any real sequence

W ( z) =

∞ 



(α )

w k,r

with its generating function

(α )

w k,r zk .

(4)

k =0

(α )

For a sufficiently smooth function f (x), denote the left or right Grünwald type operator with shift r and weights w k,r as ∞ 1 

αx∓,h,r f (x) = α h

(α )

w k,r f (x ∓ (k − r )h).

(5)

k =0

Definition 1. A sequence

(α )



w k,r

(or it’s generating function W ( z)) is said to approximate the left and right fractional

derivative D α x∓ f (x) at x with shift r in the sense of Grünwald if α Dα x∓ f (x) = lim x∓,h,r f (x). h→0

Definition 2. A sequence

(α )



w k,r

(or the generating function W ( z)) is said to approximate the fractional derivatives D α x∓ f (x)

with shift r and order p ≥ 1 if α p Dα x∓ f (x) = x∓,h,r f (x) + O (h ).

(6)

The following theorem establishes an equivalent characterization of the generator W ( z) for an approximation of the (α ) fractional differential operators D x∓ with order p ≥ 1 and shift r. We denote for conciseness the integer-order differential operators as D k =

dk ,k dxk

= 0, 1, 2, · · · .

Theorem 1 (Theorem 1, [14]). Let α ∈ R+ and m, n ∈ Z+ satisfying n − 1 < α ≤ n < m. Let f (x) ∈ C m+n+1 (R) and D k f (x) ∈ L 1 (R) (α ) for 0 ≤ k ≤ m + n + 1. Let W ( z) be the generating function of a real sequence { w k,r } and

1 G r ( z) := α W (e −z )erz z

(7)

be analytic in | z| ≤ R for some R ≥ 1. Then, the generating function W ( z) approximates the left and right fractional differential operators with shift r and order p, 1 ≤ p ≤ m, if and only if

G r ( z) = 1 + O ( z p ). ∞ Moreover, if G r ( z) = 1 + k= p ak zk , we have

αx∓,h,r f (x) = D αx∓ f (x) +

m −1 

(8)

+α hk ak D kx∓ f (x) + O (hm ).

(9)

k= p

2.2. Vandermonde matrix and a variant We also need the following results on the determinant of Vandermonde matrix and its variant. Lemma 1. For a finite sequence of parameters x0 , x1 , · · · , xq , the determinant of the Vandermonde matrix V (x0 , x1 , · · · , xq ) of size (q + 1) is a multivariate polynomial of degree q in each variable and of total degree q(q + 1)/2, and is given by

  1 1 1   x0 x1 x2  2 2 2  | V | = | V (x0 , x1 , · · · , xq )| =  x0 x1 x2  .. .. ..  .  q .q .q x x x 0

1

2

··· ··· ··· .. . ···



1   xq   2 xq =

 ..  .  q x  q

 0≤i < j ≤q

(x j − xi ).

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

55

Lemma 2. Let a variant of the Vandermonde matrix be

⎡

1 ⎢ x20 ⎢ 3 ⎢ U 2 (x0 , x1 , · · · , xq ) = ⎢ x0

⎢ ⎣

.. .

.. .

q +1

q +1

x0

⎤

··· ··· ··· .. . ···

1 x21 x31 x1

1 xq2 ⎥ ⎥ xq3 ⎥ . ⎥

.. .

⎥ ⎦

q +1

xq

The determinant of U 2 , which is of size (q + 1), is given by

⎛

⎞



|U 2 | = ⎝

0≤i < j ≤q

⎞

⎛

q q ⎟ ⎜  ⎟ ⎜ (x j − xi )⎠ ⎜ xl ⎟ . ⎠ ⎝

(10)

m =0 l =0 l =m

Proof. It is clear that |U 2 | is a multi variable homogeneous polynomial of degree q + 1 in each variable and of total degree 2 + 3 + · · · + (q + 1) = q(q + 3)/2. Moreover, for i = j, by replacing xi with x j , we have |U 2 | = 0 and thus (xi − x j ) are factors of |U 2 |, and there are q(q + 1)/2 such factors as was in the determinant | V | of the Vandermonde matrix in Lemma 1. Hence,

|U 2 | =



(x j − xi ) P (x0 , x1 , · · · , xq ),

0≤i < j ≤q

where P is a homogeneous polynomial of degree one in each variable and of total degree q(q + 3)/2 − q(q + 1)/2 = q. Since there are q + 1 variables, x j , 0 ≤ j ≤ q, P has the form, with constant A,

⎛

⎞

q q ⎜ ⎟  ⎜ ⎟ xl ⎟ . ⎝ ⎠

P = A (x1 x2 · · · xq + x0 x2 · · · xq + · · · + x1 x2 · · · xq−1 ) = A ⎜

m =0 l =0 l =m

Equating one of the terms of |U 2 |, we see that A = 1.

2

3. An explicit form In this section, we establish our main result by deriving explicit expressions for the coefficients β j for an approximating generating function of the form

W p ,r ( z) = (β0 + β1 z + · · · + β p z p )α

(11)

with a shift r and an order p. The following theorem gives a linear system of equations for the unknown coefficients β j for the generating function to approximate fractional derivatives with an order and a shift. Theorem 2. With the assumptions of Theorem 1, the generating function of the form W p ,r ( z) in (11) approximates the left and right fractional derivatives D x∓ f (x) at x with shift r and order p ≥ 1 if and only if the coefficients β j , 0 ≤ j ≤ p, govern the linear system p  (λ − j )n β j = δ1,n ,

n = 0, 1 , · · · , p ,

(12)

j =0

where λ = r /α and δ1,n is the Kronecker delta having value of one for n = 1 and zero otherwise. Proof. In view of Theorem 1, W p ,r ( z) approximates the fractional derivatives with shift r and order p if and only if

⎛ 1

1

G r ( z) = α W p ,r (e −z )erz = α ⎝ z z

p  j =0

⎞α 1

β j e − jz ⎠ erz = α z

⎛ ⎞α p  ⎝ β j e (r /α − j )z ⎠ j =0

⎛ ⎞α ⎛ ⎞α  α p p ∞ ∞    1  1 1 λ z n n n − 1 = α⎝ β je j ⎠ = ⎝ βj bn z λ z ⎠ = z z n! j j =0

j =0

n =0

n =0

56

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

 =

b0 z

+ b1 +

∞ 

α n −1

= 1 + O ( z p ),

bn z

n =2

where p 1 

bn =

n!

λnj β j ,

n = 0, 1 , 2 , · · · ,

(13)

j =0

and

λ j = λ − j.

(14)

Since G r ( z) is analytic, it does not have any poles and hence we have b0 = 0. Next, since G r ( z) = 1 + O ( z ) has the constant term 1, we must have b1 = 1. For p = 1, these two conditions give the system (12) and the proof ends. For p > 1, G r ( z) reduces to p



G r ( z) = where X =

1+

∞ 

α

n −1

=: (1 + X )α ,

bn z

n =2

∞

n−1 . Taylor expansion of G r ( z) with respect to X gives n=2 bn z

G r ( z) = 1 + α X +

α (α − 1) 2!

X2 +

α (α − 1)(α − 2) 3!

X 3 + · · · = 1 + O ( z p ).

(15)

The term z appears in the term α X at index n = 2 only on the left hand side of (15). This gives b2 = 0. The same is true for bn , n = 2, 3, · · · , p by successively comparing the coefficients of zn−1 , n = 2, 3, · · · , p in (15). Altogether, we have bn = δ1,n , n = 0, 1, 2, · · · , p which yield the system (12). 2 The system (12) can be expressed in Vandermonde matrix form as

V (λ0 , λ1 , · · · , λ p )b = d,

(16)

where b and d are column vectors given by respectively b = [β0 , β1 , · · · , β p ] and d = [0, 1, 0, · · · , 0] denoting transpose. Solving the Vandermonde system (16), we have the following main result. T

T

with superscript T

Theorem 3. The generating function W p ,r ( z) = (β0 + β1 z + · · · + β p z p )α approximates the fractional derivatives with shift r and order p if and only if the coefficients β j are given by

⎞ ⎞⎛ p p p ⎟   ⎜ 1 ⎟⎜ ⎟ ⎟⎜ βj = −⎜ (λl )⎟ ⎜ ⎠ ⎝ ⎠ j −m ⎝ ⎛

⎛

m =0 m = j

⎜ = −⎜ ⎝ p

m =0 m = j

⎞⎛

(17)

m =0 l =0 m = j l =m, j

⎞

⎜ 1 ⎟ λm ⎟ ⎟⎜ ⎟, ⎠ j −m ⎝ λm ⎠ p

j = 0, 1 , 2 , · · · , p ,

(18)

m =0 m = j

and the leading error constant is given by

Ep =

α ( p + 1)!

p  (λ − j ) p +1 β j .

(19)

j =0

Proof. Employing Cramer’s rule to solve the system (16), we have

βj =

| V j (d)| , | V (λ0 , λ1 , · · · , λ p )|

j = 0, 1 , 2 , · · · ,

(20)

where V j (d) are the matrices obtained  from V (λ0 , λ1 , · · · , λ p ) by replacing  their j-th columns with vector b. 0 1   1 0     = 1, while | V (λ0 , λ1 )| = λ1 − λ0 = (λ − 1) − = −1 and | V 1 (d)| =  When p = 1, we have | V 0 (d)| =  1 λ1  λ0 1  (λ − 0) = −1. Hence, β0 = 1, β1 = −1.

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

57

When p > 1, | V j (d)| is given by, with d in the j-th column and by using the standard determinant rule selecting the j −th column as pivot,

  1 1   λ0 λ1  2 2  | V j (d)| =  λ0 λ1  .. ..  .  p .p λ λ 1 0

··· ··· ··· .. . ···

0 ··· 1 ··· 0 ···

.. .

.. . 0 ···

   1 1  2   λ λ2 λp  1  0  2  λ p  = (−1) j +1  λ30 λ31  . ..  ..  .  . .  .  p p   λ0 λ1p λp 1 

··· ··· ··· .. . ···

1

1

λ2j −1 λ2j +1 λ3j −1 λ3j +1 .. .. . . p p λ j −1 λ j +1

··· ··· ··· .. . ···



1 

λ2p  λ3p  . ..  .  p  λp 

Note that the last determinant is of the form U 2 of size p × p missing the column with the parameter λ j . Hence, we have from (10),

⎛ ⎜ | V j (d)| = (−1) j +1 ⎜ ⎝

 0≤m
⎞ ⎞⎛ p p ⎟ ⎟⎜ ⎜  ⎟ (λn − λm )⎟ λl ⎟ . ⎜ ⎠⎝ ⎠ m =0 l =0 m = j l =m, j

The determinant of the Vandermonde matrix V (λ0 , λ1 , · · · , λ p )) can be written as

|V | =



(λn − λm ) =

0≤m


(λn − λm )(−1) j

0≤m
p 

(λm − λ j ),

m =0 m = j

where we have partitioned the product to extract the product part in | V j (d)|. We thus have from (20) and with the observation that λm − λ j = j − m,

⎞ ⎞⎛ ⎞⎛ ⎞ ⎛ p p p p p   ⎟ ⎜  1 ⎟⎜ ⎜  λ −m⎟⎜  1 ⎟ ⎟ ⎟⎜ ⎟⎜ ⎟. βj = −⎜ λl ⎟ = − ⎜ ⎜ ⎠ ⎝ ⎝ ⎠ j −m ⎝ j −m ⎠⎝ λ −m⎠ ⎛

m =0 m = j

m =0 l =0 m = j l =m, j

m =0 m = j

m =0 m = j

The leading error constant E p is given by the coefficient of z p in (15) with X given by X = n n= p bn+1 z . This is in the term α X and, from (13), is given by

∞

E p = α b p +1 =

α ( p + 1)!

p 

∞

n−1 n= p +1 bn z

=

(λ − j ) p +1 β j . 2

j =0

Remark 1. Note that formulas (17) and (18) are compatible for the case p = 1 as well although we treated it separately in the proof. Remark 2. When f is sufficiently smooth, the error constants for higher order terms up to 2p − 1 can be obtained from in (15). In fact, we have the error constants for orders n = p , p + 1, · · · , 2p − 1,

E n = α b n +1 =

α (n + 1)!

αX

p  (λ − j )n+1 β j . j =0

The error terms corresponding to the error constants E n are then given from Theorem 1 as

T n = E n hn f (n+α ) (ζ ), where ζ is in the interval containing the grid points involved. Remark 3. When the leading error constant E p is zero, the difference form has super convergence. The order and the error term are determined from the next non-zero error constant E p +1 . We list the coefficients for approximating generating functions W p ,r ( z) with order p and shift r obtained for 1 ≤ p ≤ 6 in Table 2. Note that when there is no shift, i.e., r = 0, we have λ = αr = 0 and these generating functions in Table 2 reduce to those in Table 1 obtained by Lubich [12].

58

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

Table 2 Coefficients βk of W p ,r ( z) = (β0 + β1 z + · · · + β p z p )α for order p approximation, 1 ≤ p ≤ 6. p 1

βk , 0 ≤ k ≤ p β0 = 1 β1 = −1

2

β0 = −λ + 32 β1 = −2 + 2λ β2 = −λ + 12

3

β0 =

λ2

− 2λ +

2

2

11 6

β1 = − 3λ2 + 5λ − 3 β2 =

3λ2 2

− 4λ + 2

β3 = − λ2 + λ − 4

3

β3 = β4 = β0 = β1 = β2 = β3 = β4 = β5 = 6

5λ2 λ − 35 + 25 4 12 12 2λ 9λ2 26λ − + − 4 3 2 3 −λ3 + 6λ2 − 192λ + 3 2 2λ3 − 7λ2 + 143λ − 43 3 λ3 3λ2 λ − 6 + 4 − 11 + 14 12

β0 = − λ6 + β1 = β2 =

5

3 2 1 3

3

λ4

3 2 − λ2 + 178λ − 154λ + 137 24 60 2 5λ4 7λ3 − 24 + 3 − 718λ + 776λ − 5 3 2 5λ4 λ − 133λ + 594λ − 107 +5 12 6 2 λ4 − 512 + 4λ3 − 494λ + 13λ − 10 3 3 2 5λ4 λ − 116λ + 418λ − 61 + 54 24 12 4 3 2 − λ24 + λ3 − 7λ8 + 56λ − 15 5

3 2 7λ4 λ − 3536λ + 4916λ − 203 + 49 48 45 20 5λ4 31λ3 29λ2 87λ − 6 + 6 − 2 + 5 −6 20 5 4 3 2 λ λ λ − λ8 + 9548λ − 137 + 461 − 117 + 15 12 16 4 2 λ5 5λ4 121λ3 254λ 20 2 − + − 31 λ + − 6 2 9 9 3 5 4 3 2 λ λ − λ8 + 8548λ − 107 + 307 − 332λ + 15 12 16 4 λ5 2λ4 19λ3 13λ2 27λ − + − + − 65 20 3 6 2 5 3 2 λ5 λ4 λ − 120 + 548 − 1736λ + 1516λ − 137 + 16 180

λ β0 = − 120 +

β1 = β2 = β3 = β4 = β5 = β6 =

λ5

(α )

Once the polynomial coefficients β j are computed, the Grünwald weights w k,r can be computed by a recurrence formula [18,21] given in (21) below. Lemma 3. Let P ( z) = recurrence relation

∞

w 0 = β0α , w m =

j =0 β j z

j

, where β j ∈ R,

α > 0 and W (z) = ( P (z))α =

m 1 

m β0

( j (α + 1) − m) w m− j β j .

∞

m=0

w m zm . Then, the coefficients w m satisfy the

(21)

j =1

For the generating function W p ,r ( z), being the power of a polynomial of degree p, the upper limit of the summation in (21) goes up to min(m, p ) only. We note that formulas of similar nature to the explicit formula of this paper have appeared with some variations [9,1,6]. However, we point out that those formulas are restricted to finite difference formulas for integer-order derivatives only. For example, Aceto et al. [1] considered a general finite difference form for the first order derivative and called the generalized backward differentiation formula (GBDF) and used to analyse backward difference multi-step methods where the first order derivative is for the time domain with initial and final conditions. In our work, the explicit formula is extended to fractional derivatives with respect to space variable which mostly involve boundary conditions. To the best of our knowledge, this type of general explicit form for fractional derivatives with arbitrary shifts and orders is not found in the literature. 4. Finite difference formulas for integer-order derivatives We point out that the explicit formula established in Section 3 can also be used to compute coefficients for finite difference formulas for integer-order derivatives as special cases.

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

59

In fact, when α = n, an integer, the generating function W p ,r ( z) is an integer power of a polynomial and hence has a finite series expansion. Particularly, for α = 1, the finite difference approximation with order p are the p-step GBDF class mentioned in [1]. This means that we have finite difference formulas for a given order with a shift for any integer-order derivative. The coefficients for the finite difference formulas are obtained by expanding the integer power of the polynomial obtained from Theorem 3. We list here some of well known and not well-known finite difference formulas for demonstration. Let f be a function defined on an interval which has enough derivatives to employ the finite difference forms. For notational brevity, we denote f (x + kh) = f k . 1. For the first derivative df (x)/dx, approximation of order 2 with shift 1 is obtained from (17) with parameters, 2

1, p = 2 and r = 1. The generating function is then given by W ( z) = − z2 + 1 2 h 6

1 2

and we have f  (x) =

1 h

1 2

f −1 −

α=

1 f 2 1



+

f (3) (ζ ), where ζ ∈ (x − h, x + h). This is the well known central difference formula of order 2 for the first derivative. 3

2. For order 3 with shift 2, we have W ( z) = − z3 −

 1

r = 2, is f  (x) = h

− 16 f −2

+ f −1 −

3. For order 3 with shift 3/2, we have

α = 1, p = 3 and r = 3/2, is given by

z2 2

+ z − 16 and the finite difference formula at x, with = 1, p = 3 and 1 3 (4) − − 12 h f (ζ ), where ζ ∈ (x − 2h, x + h), which is a GBDF in [1]. 2 z3 1 W ( z) = 24 − 9z8 + 9z − 24 and the finite difference formula at f (x) = f 0 , with 8  1 1 9 1 3  f (x) = h − 24 f −3/2 + 8 f −1/2 − 98 f 1/2 + 24 f 3/2 + 0h3 + 640 h4 f (5) (ζ ), which is a

1 f 2 0

1 f 3 1



α

central difference formula with functional values at midpoints of the discretized points. The error constant in this case is zero for order 3. This means that it has a super-convergence with order 4. 4. For the second derivative d2 f (x)/dx2 , the approximating generating function of order 4 with shift 2 is given by W ( z) =

2 3 2 6 5 z4 z8 z7 + z2 − 3z2 + 5z + 14 = 144 − 12 + z2 − 59z + − 12 6 36 f (x) = f 0 is, with p = 4, α = 2, r = 2, is given by



f  (x) =

1 h2



1 f 16 −2

+

5 f 12 −1

−

1 f 18 0

−

9 f 4 1

+

73z4 24

73 f 24 2

9z3 4

−

z2 18

59 f 36 3

+

1 f 2 4

−

−

+

5z 12

−

+

1 16

1 f 12 5

and the finite difference formula at

+

1 f 144 6



+

1 4 h 10

f (6) (ζ ),

where ζ ∈ (x − 2h, x + 6h). Note that this finite difference formula is not compact as it does not use minimum number of function values. However, these type non-compact formulas are applicable in computing derivatives when there are more data of function values available. We mention that finite difference formulas for integer-order derivatives with given stencil points can be computed by solving their corresponding systems of linear equations for each case [19] which is inefficient for exact computations and prone to round off errors for numerical computations. In [10], a MATLAB code is given to compute the coefficients, also by solving a system of equations for each case. A recursive algorithm to compute finite difference coefficients is also given in [6] in which the coefficients for the approximation of a higher order derivative require the coefficients of the previous order derivatives with supplied stencil points. Our explicit form, therefore, gives a unified expression for the coefficients of difference approximations for fractional and integer-order derivatives with a readily computable formula which does not require any solving process for linear systems. 5. Stability As for the stability of numerical solutions using these difference approximations in steady state and time dependent diffusion problems, we point out that only the first and second order Grünwald type approximations have been proved to be stable with shift r = 1. Particularly, the first order Grünwald approximation is shown to be stable with shift 1 for any α > 0 [13] and the second order Grünwald type approximation is proved to be stable with shift one for 1 < α ≤ 2 [14]. Stability of other higher order Grünwald approximations remains an open problem, as for the present authors’ knowledge. However, numerical results show that these higher order approximations are stable for only a limited range of values for the fractional order α . As for the approximation of integer-order derivatives, numerical results for diffusion problems confirm the result in [1] that for space derivatives also, only one shifted finite difference approximation of a particular order displays stability. 6. Conclusion Higher order approximations for fractional derivatives in the Grünwald sense with shifts are considered. The shifted Grünwald type approximation is completely characterised by a corresponding generating function from which the coefficients are obtained. Construction of approximating generating functions in the form of power of polynomials is considered and explicit forms for the coefficients for such generating functions are derived with arbitrary shift and order of accuracy. The non-shift versions of these approximations are the Lubich approximation forms given as special cases.

60

W.A. Gunarathna et al. / Applied Numerical Mathematics 143 (2019) 51–60

Moreover, the leading and successive error terms for the difference approximations are also given in easily computable explicit forms. The formula for approximating generating functions also gives finite difference formulas for any integer-order derivatives unifying the difference formula for fractional and integer-order derivatives. These explicit generating functions might be a useful tool for constructing difference approximations for fractional and integer-order derivatives with their error terms for analysis and computations. References [1] L. Aceto, D. Trigiante, On the a-stable methods in the gbdf class, Nonlinear Anal., Real World Appl. 3 (2002) 9–23. [2] B. Baeumer, M. Kovács, H. Sankaranarayanan, Higher order Grünwald approximations of fractional derivatives and fractional powers of operators, Trans. Am. Math. Soc. 367 (2015) 813–834. [3] M. Chen, W. Deng, Fourth order difference approximations for space Riemann-Liouville derivatives based on weighted and shifted Lubich difference operators, Commun. Comput. Phys. 16 (2014) 516–540. [4] M. Chen, W. Deng, Fourth order accurate scheme for the space fractional diffusion equations, SIAM J. Numer. Anal. 52 (2014) 1418–1438. [5] K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition using Differential Operators of Caputo Type, Springer Science & Business Media, 2010. [6] B. Fornberg, Generation of finite difference formulas on arbitrarily spaced grids, Math. Comput. 51 (1988) 699–706. [7] Z. Hao, Z. Sun, W. Cao, A fourth-order approximation of fractional derivatives with its applications, J. Comput. Phys. 281 (2015) 787–805. [8] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley & Sons, Inc., Interscience Publishers Inc., 1962. [9] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [10] R.J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, vol. 98, Siam, 2007. [11] C. Li, W. Deng, A new family of difference schemes for space fractional advection diffusion equation, Adv. Appl. Math. Mech. 9 (2017) 282–306. [12] C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (1986) 704–719. [13] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, Am. J. Comput. Appl. Math. 172 (2004) 65–77. [14] H.M. Nasir, K. Nafa, A new second order approximation for fractional derivatives with applications, Sultan Qaboos University, J. Sci. [SQUJS] 23 (2018) 43–55. [15] H.M. Nasir, K. Nafa, Algebraic construction of a third order difference approximations for fractional derivatives and applications, ANZIAM J. 59 (2018) C231–C245. [16] H.M. Nasir, B.L.K. Gunawardana, H.M.N.P. Aberathna, A second order finite difference approximation for the fractional diffusion equation, Int. J. Appl. Phys. Math. 3 (2013) 237. [17] I. Podlubny, 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. [18] L.B. Rall, Automatic Differentiation: Techniques and Applications, Springer, 1981. [19] C. Taylor, Finite difference coefficients calculator, http://web.media.mit.edu/~crtaylor/calculator.html, 2017. (Accessed 30 September 2018). [20] W. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comput. 84 (2015) 1703–1727. [21] M. Weilbeer, Efficient Numerical Methods for Fractional Differential Equations and Their Analytical Background, Papierflieger, 2005. [22] Y. Yu, W. Deng, Y. Wu, J. Wu, Third order difference schemes (without using points outside of the domain) for one sided space tempered fractional partial differential equations, Appl. Numer. Math. 112 (2017) 126–145. [23] H. Zhou, W. Tian, W. Deng, Quasi-compact finite difference schemes for space fractional diffusion equations, J. Sci. Comput. 56 (2013) 45–66.