# A New Grünwald-Letnikov Derivative Derived from a Second-Order Scheme

Authors: B. A. Jacobs
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

A novel derivation of a second-order accurate Grünwald-Letnikov-type approximation to the fractional derivative of a function is presented. This scheme is shown to be second-order accurate under certain modifications to account for poor accuracy in approximating the asymptotic behavior near the lower limit of differentiation. Some example functions are chosen and numerical results are presented to illustrate the efficacy of this new method over some other popular choices for discretizing fractional derivatives.


---

## Body

## 1. Introduction

Fractional differential equations and fractional order derivatives have become increasingly popular in recent years. The use of fractional order derivatives has become abundant in modeling techniques as well as computational methods for the numerical solution of these models. Zeng et al. [1] derive a generalized fractional diffusion equation by addressing deficiencies in previously proposed models for anomalous diffusion. An application of modelling plumes moving at a fractional rate was investigated by Benson et al. [2] which showed a fractional advection-dispersion model to be appropriate as the mean squared displacement followed t
1
/
α. He [3] presented a more exact model for seepage flow by relaxing the continuity assumption and went on to solve the arising equation using the variational iteration method (VIM). Metzler et al. [4] provide a review of anomalous diffusion equations presented as well as the formulation of a general diffusion case overcoming previous inconsistencies. Moving past anomalous diffusion Gafiychuk et al. [5, 6] studied a fractional reaction-diffusion system modelling an activator-inhibitor dynamic. The application of fractional order models is various and with this increased complexity requires both analytic and numerical methods to find useful solutions.Meerschaert and Tadjeran [7] construct an implicit scheme for two-sided space-fractional partial differential equations based on two first-order approximations to the left- and right-shifted Grünwald-Letnikov fractional derivative. Scherer et al. [8] also derive an implicit scheme for the fractional heat equation in the Caputo sense, making use of the standard Grünwald-Letnikov scheme as well as the relationship between the Riemann-Liouville derivative and the Caputo derivative. There has been work conducted on comparing numerical methods with analytical methods; El-Sayed et al. [9], for example, compare a modified numerical method, wherein the authors apply a transform to homogenize the initial conditions, with the Adomian decomposition method. Another approach taken by Karatay et al. [10] derives an unconditionally stable implicit scheme based on a first-order approximation to the Caputo derivative. An interesting study is presented by Su et al. [11] which combines the first-order shifted Grünwald-Letnikov scheme with characteristic methods and a departure from the Riemann-Liouville and Caputo derivatives leads Chen et al. [12] to derive a conditionally stable explicit scheme and unconditionally stable implicit scheme for the Caputo-Riesz fractional diffusion equation that are superlinearly convergent.Tadjeran and Meerschaert [13] implement a second-order accurate Crank-Nicolson-type finite difference scheme for the space-fractional diffusion equation, making use of the right-shifted Grünwald-Letnikov fractional derivative. Diethelm et al. [14] propose a predictor-corrector-type approach for the numerical solution of both linear and nonlinear fractional differential equations including systems of equations. Mohebbi et al. [15] study a high-order accurate numerical method for obtaining the solution to the subdiffusion equation with a nonlinear source term. Cui [16] goes on to solve the two-dimensional time-fractional diffusion equation using a high-order compact alternating direction implicit method. Meerschaert et al. [17] solve the two-dimensional space-fractional dispersion equation using a right-shifted Grünwald-Letnikov discretization since the order of the fractional derivative is above 1 but less than 2. Lynch et al. [18] implement a finite difference approach for the space-fractional dispersion problem but implement discretization switching to maximize accuracy based on the order of derivative. Momani and Odibat [19] depart from the finite-difference approach and instead compare two analytic methods for the solution of linear fractional differential equations, namely, the variational iteration method and the Adomian decomposition method. These authors go on to explore the application of these methods to nonlinear partial differential equations of fractional order in [20].The present work derives a difference approximation in a novel way that produces a form similar to the known result presented by Oldham and Spanier in [21]. This text cites an error in the approximation of the fractional derivative of a constant and linear function of order N
-
1 that is proportional to the number of terms in the series approximation. A novel technique is coupled with this approximation to attain order N
-
2 accuracy. Although the result obtained is known, to the best of the author’s knowledge this derivation is novel.Another important technique used in the numerical application of fractional differential equations is the “short-memory” principle [22]. This technique truncates the approximation sum to avoid high computational costs. It is shown in Section 4 that the present method is highly accurate even with very few terms in the approximation sum.The following section presents some preliminary results used in the remainder of the paper. Section3 goes on to derive the new method of discretizing a fractional derivative with results on the performance of this technique reported in Section 4. We make some concluding remarks and discuss opportunities for further work in Section 5.


## 2. Preliminaries

This section presents some preliminary results that are put to use in subsequent sections. However, in the interest of brevity, we omit some useful properties but direct the interested reader to the work by Li and Deng [23]. Often for numerical simulation of fractional order problems the Grünwald-Letnikov derivative is used as a means of discretizing a fractional order derivative. The Grünwald-Letnikov discretization has been used for numerical schemes for fractional partial differential equations [13, 24, 25] and is defined as follows.Definition 1.
The Grünwald-Letnikov derivative of orderα
>
0 of a function u
(
t
) is (1)DGL0,tαut=limΔt→01Δtα∑k=0Nωkαut-kΔt,where (2)ωkα=-1kαkand N
Δ
t
=
t.However the limit definition is only practically useful for finite-difference implementations; hence the following definition is used.Definition 2.
The Grünwald-Letnikov derivative of orderα
>
0 of a function u
(
t
)
∈
C
m
[
0
,
t
] is [22] (3)D0,tαut=∑k=0m-1uk0t-a+kΓ-α+k+1+1m-α·∫0tt-τm-α-1umτdτ,where m
-
1
≤
α
<
m
∈
Z
+.Definition2 is used as a means of calculating the exact fractional derivative of a function. The exact solution obtained here is then used as a measure for accuracy when compared with the discrete approximations. It is prudent to mention that if u
(
t
) is sufficiently smooth, u
(
t
)
∈
C
m
[
0
,
t
], then the Grünwald-Letnikov derivative is equivalent to the popular Riemann-Liouville defined below.Definition 3.
The Riemann-Liouville derivative of fractional orderα
>
0 of a function u
(
t
) is [22] (4)DRL0,tαut=1Γm-αdmdtm∫0tt-τm-α-1uτdτ.The right-shifted Grünwald-Letnikov approximation introduced by Meerschaert and Tadjeran [25], as used in [13], leads to stable numerical schemes and is defined below.Definition 4.
The right-shifted Grünwald-Letnikov approximation of a functionu
(
t
) for 1
<
α
≤
2 is [25] (5)DRSGL0,tαut=limN→∞1Δtα∑k=0N-1kαkut-k-1Δt.Oldham and Spanier obtain a second-order Grünwald-Letnikov approximation simply by noting that the Riemann sum present in the derivative would converge faster if a simple modification was made.Definition 5.
The Grünwald-Letnikov derivative of orderα
>
0 of a function u
(
t
) with a central difference modification is [21] (6)DGL0,tαut=limΔt→01Δtα∑k=0N-1-1kαkut-kΔt+αΔt2,where N
Δ
t
=
t.The primary difference between Definition5 and the present work is the derivation and the choice of N, the upper limit of the Riemann sum. Although these discrepancies are slight we also address the issue of the asymptotic behavior of the fractional derivative as t approaches zero, a difficult behavior for a difference scheme to capture.


## 3. Method

In much the same way the Grünwald-Letnikov derivative is generalized from the first-order backward difference scheme, as described by Podlubny in [22], we derive a new approximate fractional derivative in the form of a difference scheme which is derived from the second-order central difference scheme. Yuste and Acedo [24] state that the standard Grünwald-Letnikov derivative is accurate to O
(
Δ
t
), while we derive a scheme based on the central-difference scheme; we hope to achieve an approximation that is O
(
Δ
t
2
). Podlubny [26] shows, by example, that the Grünwald-Letnikov derivative of a function f
(
t
)
=
1 or f
(
t
)
=
t
m is accurate to O
(
Δ
t
) and hence conclude that the Grünwald-Letnikov derivative of any function that can be written as a power series will be accurate to O
(
Δ
t
). The central difference approximation to the first-order derivative is given by (7)u′t=limΔt→0ut+Δt-ut-Δt2Δt+OΔt2or (8)u′t=limΔt→0ut+Δt/2-ut-Δt/2Δt+OΔt222.Similarly the second-order derivative can be approximated by (9)u′′t=limΔt→0ut+Δt-2ut+ut-ΔtΔt2+OΔt2.Following this trend we may write (10)unt=limΔt→01Δtn∑k=0n-1knkut-kΔt+nΔt2.Generalizing the order of the derivative to be a positive real value α and substituting into (10) we obtain a general formula: (11)DGL20,tαut=uαt=limΔt→01Δtα∑k=0N-1kαkut-kΔt+αΔt2,where the binomial coefficient is calculated by the use of the Gamma function that generalizes the factorial operator: (12)αk=α!k!α-k!=Γα+1Γk+1Γα-k+1.It is well known [26, 27] that the fractional differentiation operator is a nonlocal operator. This effect is seen in the discretization being affected by every point in the function and not just a neighborhood of the point t. This insight guides our choice of N as the upper limit of the summation. We choose N so that the last term in the series corresponds with the terminating point of the function u
(
0
), and in terms of fractional differential equations and fractional partial differential equations this would correspond with a boundary or initial condition: (13)N=tΔt+α2,where · represents the ceiling function. Due to the singularity of a fractional derivative D
G
L
0
,
t
α
u
(
t
) when t
=
0 it is convenient to select u
(
t
) such that u
0
=
0 to alleviate the asymptotic behavior. However we do not want to restrict ourselves only to functions that obey this property. Hence we define (14)gt=ut-u0,so that g
(
0
)
=
0 and then (15)DGL20,tαut=DGL20,tαgt+DGL20,tαu0,since the fractional derivative is a linear operator. We may now exploit Property  2.2 in [23] which yields (16)LDGL20,tαu0=sαLu0,which then leads to (17)L-1sαLu0=t-αu0Γ1-α.Finally (15) becomes (18)DGL20,tαut=DGL20,tαgt+t-αu0Γ1-α.While correcting for the asymptotic behavior at t
=
0 produces sufficiently accurate results, we may go one step further and substitute the exact solution of the linear component of u
(
t
). This is done to ensure that the accuracy of the scheme is at least O
(
N
-
2
). Let(19)ht=ut-u0-tu′0;then (20)DGL20,tαut=DGL20,tαgt+t-αu0Γ1-α+t1-αu′0Γ2-α.


### 3.1. Accuracy

As described by Oldham and Spanier in [21] the error in the above scheme for functions u
(
t
)
=
C, u
(
t
)
=
t, and u
(
t
)
=
t
2 is given by α
(
α
+
1
)
/
2
N, α
2
/
2
N, and α
(
α
-
1
)
(
α
-
2
)
/
24
N
2, respectively. Since we make the exact substitutions for the linear components we are able to increase the accuracy of the scheme to O
(
N
-
2
) for any analytic function.Theorem 6.
The error ofD
G
L
2
0
,
t
α
u
(
t
), denoted by ϵ
(
u
t
)
G
L
2
=
D
0
,
t
α
u
t
-
D
G
L
2
0
,
t
α
u
t is at least O
(
N
-
2
) accurate for the function u
(
t
)
=
t
p for all p
≥
2 given that the exact derivative of order α of t
p is (21)D0,tαtp=tp-αΓp+1Γp-α+1.Proof.
Letp
=
2. Then (see [21], Chapter 8) (while Oldham and Spanier cite the relative error it is easy to show that the more restrictive, absolute error is also O
(
N
-
2
). More information can be found in [22], Chapter 7)(22)ϵutGL2=αα-1α-224N2.Assume that the order of error is at least O
(
N
-
2
) for p
=
m; that is, (23)DGL20,tαtm=tm-αΓm+1Γm-α+1+ON-2.Now let p
=
m
+
1; then by the product rule for fractional derivatives (see [21], Chapter 5) we have (24)DGL20,tαtm+1=DGL20,tαtmt=tGL2D0,tαtm+αGL2D0,tα-1tm≈ttm-αΓm+1Γm-α+1+αtm-α+1Γm+1Γm-α+ON-2=tm-α+1Γm+1Γm-α+11+αm+1-α+ON-2.Then using the recursive property of the Gamma function x
Γ
(
x
)
=
Γ
(
x
+
1
) we have (25)DGL20,tαtm+1=tm+1-αΓm+2Γm-α+2+ON-2.Using Theorem6 we can deduce that the accuracy of D
G
L
2
0
,
t
α
u
(
t
) for any analytic function u
(
t
) is O
(
N
-
2
) provided that the linear and constant components of that function are dealt with accordingly.


## 3.1. Accuracy

As described by Oldham and Spanier in [21] the error in the above scheme for functions u
(
t
)
=
C, u
(
t
)
=
t, and u
(
t
)
=
t
2 is given by α
(
α
+
1
)
/
2
N, α
2
/
2
N, and α
(
α
-
1
)
(
α
-
2
)
/
24
N
2, respectively. Since we make the exact substitutions for the linear components we are able to increase the accuracy of the scheme to O
(
N
-
2
) for any analytic function.Theorem 6.
The error ofD
G
L
2
0
,
t
α
u
(
t
), denoted by ϵ
(
u
t
)
G
L
2
=
D
0
,
t
α
u
t
-
D
G
L
2
0
,
t
α
u
t is at least O
(
N
-
2
) accurate for the function u
(
t
)
=
t
p for all p
≥
2 given that the exact derivative of order α of t
p is (21)D0,tαtp=tp-αΓp+1Γp-α+1.Proof.
Letp
=
2. Then (see [21], Chapter 8) (while Oldham and Spanier cite the relative error it is easy to show that the more restrictive, absolute error is also O
(
N
-
2
). More information can be found in [22], Chapter 7)(22)ϵutGL2=αα-1α-224N2.Assume that the order of error is at least O
(
N
-
2
) for p
=
m; that is, (23)DGL20,tαtm=tm-αΓm+1Γm-α+1+ON-2.Now let p
=
m
+
1; then by the product rule for fractional derivatives (see [21], Chapter 5) we have (24)DGL20,tαtm+1=DGL20,tαtmt=tGL2D0,tαtm+αGL2D0,tα-1tm≈ttm-αΓm+1Γm-α+1+αtm-α+1Γm+1Γm-α+ON-2=tm-α+1Γm+1Γm-α+11+αm+1-α+ON-2.Then using the recursive property of the Gamma function x
Γ
(
x
)
=
Γ
(
x
+
1
) we have (25)DGL20,tαtm+1=tm+1-αΓm+2Γm-α+2+ON-2.Using Theorem6 we can deduce that the accuracy of D
G
L
2
0
,
t
α
u
(
t
) for any analytic function u
(
t
) is O
(
N
-
2
) provided that the linear and constant components of that function are dealt with accordingly.


## 4. Results

This section presents some numerical results for different functions. Each example shows the absolute error over the domain of the function, the consistency of each derivative ash tends to zero, the absolute error of the truncated series, and the absolute error over a range of α. We compare the aforementioned criteria for four different methods: the standard Grünwald-Letnikov approximation D
G
L
0
,
t
α, the second-order Grünwald-Letnikov approximation D
G
L
2
0
,
t
α, the second-order Grünwald-Letnikov approximation with constant correction described in Section 3  C
C
G
L
2
D
0
,
t
α, and finally the right-shifted Grünwald-Letnikov approximation D
R
S
G
L
0
,
t
α.


### 4.1. Example  1

Let(26)ut=t2,and let the order of derivative be α
=
0.75, with time step Δ
t
=
0.001. Since this choice of u
(
t
) satisfies u
(
0
)
=
0 the approximation does not suffer from poor accuracy near t
=
0 and hence the second-order Grünwald-Letnikov D
G
L
2
0
,
t
α and the second-order Grünwald-Letnikov with constant correction C
C
G
L
2
D
0
,
t
α approximations are similar.


### 4.2. Example  2

Let(27)ut=cos⁡t,the order of derivative α
=
1.5, and Δ
t
=
0.001. As mentioned by Oldham and Spanier [21] this method is robust for all α; this example shows highly accurate results for α
>
1.


### 4.3. Example  3

Let(28)ut=cos⁡t+sin⁡tt2,the order of derivative α
=
0.9, and again Δ
t
=
0.001. This function suffers from the asymptotic behavior of the fractional derivative near t
=
0 and hence the constant correction method enjoys a vastly superior accuracy over the other methods.


### 4.4. Example  4

Let(29)ut=e-tcos⁡t,the order of derivative α
=
2.1, and Δ
t
=
0.001.


### 4.5. Example  5

Let(30)ut=sin⁡t+cos⁡t,the order of derivative α
=
0.75, and Δ
t
=
0.001.The above results indicate that the present method enjoys a high accuracy. The nature of the fractional derivative near the lower limit of differentiation results in a high gradient. Difference schemes with fixed interval steps are unable to accurately capture very steep gradients, similar to those at the asymptote of the fractional derivative. The use of the constant correction allows one to avoid the asymptotic behavior neart
=
0 and the provable accuracy.Figures1, 5, 9, 13, and 17 show the increasing accuracy of the approximations as t deviates from 0, as well as the robustness of our method to the asymptotic behavior. The consistency of the methods is evident in Figures 2, 6, 10, 14, and 18. These figures indicate that as h tends to 0 the absolute error of the approximations also tends to 0. As mentioned in Section 1 the “short-memory” principle allows for the truncation of the approximating series in order to alleviate computational cost; our present method is highly accurate even for small values of N as is evident in Figures 3, 7, 11, 15, and 19. The final set of figures, Figures 4, 8, 12, 16, and 20, indicate the robustness of our method to varying values of α. It is important to note here the log scale of the graph.Figure 1
Absolute error over time domain for four different methods.Figure 2
Absolute error over varying values ofΔ
t for four different methods.Figure 3
Absolute error for increasingN for four different methods.Figure 4
Absolute error over varying values ofα for four different methods.Figure 5
Absolute error over time domain for four different methods.Figure 6
Absolute error over varying values ofΔ
t for four different methods.Figure 7
Absolute error for increasingN for four different methods.Figure 8
Absolute error over varying values ofα for four different methods.Figure 9
Absolute error over time domain for four different methods.Figure 10
Absolute error over varying values ofΔ
t for four different methods.Figure 11
Absolute error for increasingN for four different methods.Figure 12
Absolute error over varying values ofα for four different methods.Figure 13
Absolute error over time domain for four different methods.Figure 14
Absolute error over varying values ofΔ
t for four different methods.Figure 15
Absolute error for increasingN for four different methods.Figure 16
Absolute error over varying values ofα for four different methods.Figure 17
Absolute error over time domain for four different methods.Figure 18
Absolute error over varying values ofΔ
t for four different methods.Figure 19
Absolute error for increasingN for four different methods.Figure 20
Absolute error over varying values ofα for four different methods.


## 4.1. Example  1

Let(26)ut=t2,and let the order of derivative be α
=
0.75, with time step Δ
t
=
0.001. Since this choice of u
(
t
) satisfies u
(
0
)
=
0 the approximation does not suffer from poor accuracy near t
=
0 and hence the second-order Grünwald-Letnikov D
G
L
2
0
,
t
α and the second-order Grünwald-Letnikov with constant correction C
C
G
L
2
D
0
,
t
α approximations are similar.


## 4.2. Example  2

Let(27)ut=cos⁡t,the order of derivative α
=
1.5, and Δ
t
=
0.001. As mentioned by Oldham and Spanier [21] this method is robust for all α; this example shows highly accurate results for α
>
1.


## 4.3. Example  3

Let(28)ut=cos⁡t+sin⁡tt2,the order of derivative α
=
0.9, and again Δ
t
=
0.001. This function suffers from the asymptotic behavior of the fractional derivative near t
=
0 and hence the constant correction method enjoys a vastly superior accuracy over the other methods.


## 4.4. Example  4

Let(29)ut=e-tcos⁡t,the order of derivative α
=
2.1, and Δ
t
=
0.001.


## 4.5. Example  5

Let(30)ut=sin⁡t+cos⁡t,the order of derivative α
=
0.75, and Δ
t
=
0.001.The above results indicate that the present method enjoys a high accuracy. The nature of the fractional derivative near the lower limit of differentiation results in a high gradient. Difference schemes with fixed interval steps are unable to accurately capture very steep gradients, similar to those at the asymptote of the fractional derivative. The use of the constant correction allows one to avoid the asymptotic behavior neart
=
0 and the provable accuracy.Figures1, 5, 9, 13, and 17 show the increasing accuracy of the approximations as t deviates from 0, as well as the robustness of our method to the asymptotic behavior. The consistency of the methods is evident in Figures 2, 6, 10, 14, and 18. These figures indicate that as h tends to 0 the absolute error of the approximations also tends to 0. As mentioned in Section 1 the “short-memory” principle allows for the truncation of the approximating series in order to alleviate computational cost; our present method is highly accurate even for small values of N as is evident in Figures 3, 7, 11, 15, and 19. The final set of figures, Figures 4, 8, 12, 16, and 20, indicate the robustness of our method to varying values of α. It is important to note here the log scale of the graph.Figure 1
Absolute error over time domain for four different methods.Figure 2
Absolute error over varying values ofΔ
t for four different methods.Figure 3
Absolute error for increasingN for four different methods.Figure 4
Absolute error over varying values ofα for four different methods.Figure 5
Absolute error over time domain for four different methods.Figure 6
Absolute error over varying values ofΔ
t for four different methods.Figure 7
Absolute error for increasingN for four different methods.Figure 8
Absolute error over varying values ofα for four different methods.Figure 9
Absolute error over time domain for four different methods.Figure 10
Absolute error over varying values ofΔ
t for four different methods.Figure 11
Absolute error for increasingN for four different methods.Figure 12
Absolute error over varying values ofα for four different methods.Figure 13
Absolute error over time domain for four different methods.Figure 14
Absolute error over varying values ofΔ
t for four different methods.Figure 15
Absolute error for increasingN for four different methods.Figure 16
Absolute error over varying values ofα for four different methods.Figure 17
Absolute error over time domain for four different methods.Figure 18
Absolute error over varying values ofΔ
t for four different methods.Figure 19
Absolute error for increasingN for four different methods.Figure 20
Absolute error over varying values ofα for four different methods.


## 5. Conclusion

This work has shown a novel derivation of a second-order accurate scheme for approximating the fractional order derivative of a function. To the best of the author’s knowledge performing a linear transformation on the function to avoid asymptotic behavior near the origin is also novel. Section4 provides several examples with varying parameters and behaviors to illustrate the efficacy of the present method.Using a higher-order accurate approximation leads to the potential for higher-order accurate numerical schemes for fractional differential equations and fractional partial differential equations. Tadjeran and Meerschaert [13] implement the right-shifted Grünwald-Letnikov derivative for fractional derivatives of order 1
<
α
≤
2 which lead to an unconditionally stable numerical scheme. According to Oldham and Spanier [21] the present method is robust for all values ofalpha. Potential extensions to this work include constructing numerical schemes centered on this approximation and the fractional order derivative for nonlinear fractional differential equations as well as time-marching schemes for time or space fractional partial differential equations.The accuracy of the method has been proven to be at least ofO
(
N
-
2
), which suggests that the accuracy should be of order 1
0
-
6 for Δ
t. However in these examples the error is observed to be in the region of 1
0
-
8. Nevertheless the mathematical proof of accuracy as well as the compelling results presented substantiates the efforts toward this field.



---

*Source: 952057-2015-08-26.xml*