A weighted finite difference method for the fractional diffusion equation based on the Riemann-Liouville derivative

A weighted finite difference method for the fractional diffusion equation based on the Riemann-Liouville derivative


A one dimensional fractional diffusion model with the Riemann-Liouville fractional derivative is studied. First, a second order discretization for this derivative is presented and then an unconditionally stable weighted average finite difference method is derived. The stability of this scheme is established by von Neumann analysis. Some numerical results are shown, which demonstrate the efficiency and convergence of the method. Additionally, some physical properties of this fractional diffusion system are simulated, which further confirm the effectiveness of our method.

fractional diffusion equations, Riemann-Liouville derivative, weighted average methods, von Neumann stability analysis

1 Introduction

Recently, a large number of applied problems have been formulated on fractional differential equations and consequently considerable attention has been given to the solutions of those equations. Fractional space derivatives are used to model anomalous diffusion or dispersion, a phenomenon observed in many problems. There are some diffusion processes for which the Fick’s second law fails to describe the related transport behavior. This phenomenon is called anomalous diffusion, which is characterized by the nonlinear growth of the mean square displacement, of a diffusion particle over time. The anomalous diffusions differ according to the values of , where is the order of the fractional derivative. Some works providing an introduction to fractional calculus related to diffusion problems are, for instance, [2, 6, 11, 12, 28, 29]. In this work we will be interested in the anomalous diffusion, called supperdiffusion, for and experimental evidence of this type of diffusion is already reported in several works [1, 7, 13, 14].

Fractional derivatives are non-local opposed to the local behaviour of integer derivatives. Therefore, different challenges appear when we try to derive numerical methods for this type of equations. Numerical approaches to different types of fractional diffusion models are increasingly appearing in literature. We can found recent work on numerical solutions for the fractional diffusion equation describing superdiffusion [5, 9, 10, 19, 15, 27, 16] and also for several transport equations including this type of diffusion [18, 25, 31]. Some other works consider subdiffusion, which is represented by a time fractional derivative of positive order and less than one [3, 30]. However, the challenges for these equations are different from the ones that arise when we consider a space fractional derivative of order .

Numerical methods, for models with superdiffusion, have been obtained with mathematical techniques which do not necessarily consider a second order discretization for the fractional derivative to achieve second order accuracy. In this work, we present a second order approximation for the fractional Riemann-Liouville derivative of order , . This approach uses some of the tools described in [4, 8] and also applied in [26] to derive an approximation for the Caputo fractional derivative defined in bounded domains. Here, we consider the Riemann-Liouville fractional derivative in an unbounded domain and its discretization is represented by a series instead of a finite sum. We prove the order of consistency of this discretization is second order.

A weighted average finite difference -scheme is considered, for , which includes the Crank-Nicolson method () and the back forward Euler method (). The consistency and stability of the -scheme are established and we prove the -scheme is unconditionally stable. Also for we have second order accuracy in time and space as expected.

Consider the one-dimensional fractional diffusion equation [1, 7, 16]


on the domain , where and , subject to the initial condition


and to the boundary condition


The usual way of representing the fractional derivatives is by the Riemann-Liouville formula. The Riemann-Liouville fractional derivative of order , for , , is defined by


where is the Gamma function and , with denoting the integer part of .

The function under consideration, that is, which is solution of (1), should be such that the corresponding integral (4) converges. If the function vanishes at infinity, as assumed when we impose the boundary condition (3), we have absolute convergence of such integrals for a wide class of functions [24]. However, these functions do not necessarily need to vanish at infinity and we can found under which conditions these integrals converge in [24] (section 14.3). There are very complete works about the fractional calculus [17, 20, 21, 22, 24], where the theoretical properties of this type of derivative are studied in detail.

Another way to represent the fractional derivative is by the Grünwald-Letnikov formula, that is,


The Grünwald-Letnikov approximation is often used to numerically approximate the Riemann-Liouville derivative and it was the first algorithm to appear for approximating fractional derivatives [21, 22]. However, this approximation has consistency of order one and also very frequently numerical approximations based in this formula originate unstable numerical methods and henceforth a shifted Grünwald-Letnikov formula is used [16, 18].

The plan of the paper is as follows. In section 2 we derive a numerical approximation for the Riemann-Liouville derivative. The full discretisation of the fractional diffusion equation is given in section 3, where a weighted finite difference method in time is applied with the weight . In section 4 we prove the convergence of the numerical method by showing consistency and stability. In the fifth section we present numerical results which confirm the theoretical results and in the last section we give some conclusions.

2 The numerical method

In this section we present a numerical approximation for the Riemann-Liouville derivative and also the numerical method that gives an approximate solution to the fractional diffusion equation.

2.1 Approximation of the Riemann-Liouville derivative

Let us consider the Riemann-Liouville derivative [21, 22], that is,


We define the mesh points where denotes the uniform space step. For a fixed time , let us denote


First, we do the following approximation at

For each we need to calculate .

We compute these integrals by approximating , at a fixed instant , by a linear spline , whose nodes and knots are chosen at , , that is, an approximation to becomes defined by


The spline interpolates the points and is of the form [23]


with , in each interval , for , given by


and for ,


From (8) and (9),


We have that


where the are such that,




and an approximation for , is given by,


that is,

Let us assume there are approximations to the values , where and is the uniform time-step.

We define the fractional operator as




Therefore, an approximation of (6), for , can be given by

We can also write the fractional operator (17) as


Remark: Note that for and the coefficients (18) are such that , for . For , , , and for , , , .

Remark: The series (19) converges absolutely for each and for every bounded function , for a fixed . This result is a straightforward consequence of some results given in section 3 about the convergence of the series of the .

In this section we have considered a linear spline to approximate the integral representation of the Riemann-Liouville derivative with the purpose of obtaining a second order approximation. In the next section we describe the full discretisation of the differential equation.

2.2 Weighted average finite difference methods

We discretize the spatial -order derivative following the steps of the previous section. The discretization in time consists of the weighted average discretization.

We consider the time discretization . Additionally, let , . For the uniform space step and time step , let

From equation (1) we can arrive at the explicit Euler and implicit Euler numerical methods, respectively


Let (20) multiplies and (21) multiplies . We obtain the following weighted -scheme


where .

Note that for , the operator (17) is the central second order operator , that is,

We have the following numerical method



3 Convergence of the numerical scheme

In this section we prove the convergence of the numerical method by showing it is consistent and von Neumann stable. First, we start to study the consistency of the numerical method and lastly we present the stability results.

3.1 Consistency

In the beginning of this section, for the sake of clarity, we omit the variable and we denote the partial derivative of in of order by .

Lemma 1

Let . For ,


Proof: For ,

Using Taylor expansions, we obtain

where are functions which depend on and , given by


It is easy to conclude that , for .

Theorem 2

(Order of accuracy of the approximation for the fractional derivative): Let and such that , for , being a real constant. We have that


and and are independent of .

Proof: It is straightforward to prove that we have

where .

Let us define the error , such that,

We have

that is


We are now going to compute the error . We have

Taking in consideration the previous lemma, let us denote


where are defined as follows. For and ,


and for


For by changing variables, we obtain

that is,

Let such that , for . For we have


Since, by Lemma 1,


we have


For ,


We have


Finally for , we bound each integral of (28) separately. For the first integral we have

Therefore, since for , we have

Similarly, for the second integral we have

and for the third integral

Finally, we have


From (29), (30) and (31) it is easy to conclude that the error defined by (45) is of order and therefore the is of order .

Theorem 3

The truncation error of the weighted numerical method (23) is of order , where , for and , for .

Proof: Let be a solution to the fractional partial differential equation and satisfying the conditions of the previous theorem. Note that the truncation error for the numerical method (23) is given by

We have that


and using the previous theorem we have



3.2 Fourier decomposition of the error

In order to derive stability conditions for the finite difference schemes, we apply the von Neumann analysis or Fourier analysis. Fourier analysis assumes that we have a solution defined in the whole real line. It is also applied to problems defined in finite domains with periodic boundary conditions since the solution is seen as a periodic function in .

If is the exact solution , let


be the error at time level in mesh point . To apply the von Neumann analysis we also consider locally constant, and we denote by .

Considering the scheme (23) and inserting equation (