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

## Abstract

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.

###### keywords:
fractional diffusion equations, Riemann-Liouville derivative, weighted average methods, von Neumann stability analysis
1

## 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]

 ∂u∂t(x,t)=d(x)∂αu∂xα(x,t)+p(x,t) (1)

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

 u(x,0)=f(x), x∈IR (2)

and to the boundary condition

 u(x,t)=0as|x|→∞. (3)

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

 ∂αu∂xα(x,t)=1Γ(n−α)∂n∂xn∫xau(ξ,t)(x−ξ)n−α−1dξ,(n−1<α

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,

 ∂αu∂xα(x,t)=limΔx→01Δxα[x−aΔx]∑k=0(−1)k(αk)u(x−kΔx,t).(α>0) (5)

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,

 ∂αu∂xα(x,t)=1Γ(2−α)∂2∂x2∫x−∞u(ξ,t)(x−ξ)1−αdξ,1<α<2. (6)

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

 Iα(x)=∫x−∞u(ξ,t)(x−ξ)1−αdξ. (7)

First, we do the following approximation at

 ∂2∂x2Iα(xj)≃1Δx2[Iα(xj−1)−2Iα(xj)+Iα(xj+1)].

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

 Iα(xj)=∫xj−∞sj(ξ)(xj−ξ)1−αdξ. (8)

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

 sj(ξ)=j∑k=−∞u(xk,t)sj,k(ξ), (9)

with , in each interval , for , given by

 sj,k(ξ)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ξ−xk−1xk−xk−1,xk−1≤ξ≤xkxk+1−ξxk+1−xk,xk≤ξ≤xk+10otherwise, (10)

and for ,

 sj,j(ξ)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ξ−xj−1xj−xj−1,xj−1≤ξ≤xj0otherwise. (11)

From (8) and (9),

 Iα(xj)=j∑k=−∞u(xk,t)∫xk+1xk−1sj,k(ξ)(xj−ξ)1−αdξ. (12)

We have that

 ∫xk+1xk−1sj,k(ξ)(xj−ξ)1−αdξ = ∫xkxk−1ξ−xk−1Δx(xj−ξ)1−α+∫xk+1xkxk+1−ξΔx(xj−ξ)1−α (13) = Δx2−α(2−α)(3−α)aj,k,

where the are such that,

 aj,k=⎧⎪⎨⎪⎩(j−k+1)3−α−2(j−k)3−α+(j−k−1)3−α, k≤j−11, k=j. (14)

Therefore,

 Iα(xj)=Δx2−α(2−α)(3−α)j∑k=−∞u(xk,t)aj,k, (15)

and an approximation for , is given by,

 1Δx2[Iα(xj−1)−2Iα(xj)+Iα(xj+1)] (16)

that is,

 Δx−α(2−α)(3−α)[j−1∑k=−∞u(xk,t)aj−1,k−2j∑k=−∞u(xk,t)aj,k+j+1∑k=−∞u(xk,t)aj+1,k].

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

We define the fractional operator as

 δαUnj=1Γ(4−α){j+1∑k=−∞qj,kUnk}, (17)

where

 qj,k = aj−1,k−2aj,k+aj+1,k,k≤j−1 qj,j = −2aj,j+aj+1,j qj,j+1 = aj+1,j+1. (18)

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

We can also write the fractional operator (17) as

 δαUnj=1Γ(4−α)∞∑m=−1qj,j−mUnj−m. (19)

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

 μαj=djΔtΔxα.

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

 Un+1j−UnjΔt=djΔxαδαUnj+pnj, (20)
 Un+1j−UnjΔt=djΔxαδαUn+1j+pn+1j, (21)

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

 Un+1j−Unj=μαj{(1−τ)δαUnj+τδαUn+1j}+τΔtpn+1j+(1−τ)Δtpnj, (22)

where .

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

 δαUnj=Unj+1−2Unj+Unj−1.

We have the following numerical method

 (1−τμαjδα)Un+1j=(1+(1−τ)μαjδα)Unj+Δtpn+τj, (23)

where

 pn+τj=τpnj+(1−τ)pn+1j.

## 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 ,

 u(ξ)−sj,k(ξ)=−1r!3∑r=2u(r)(ξ)lk,r(ξ)−14!u(4)(ηk)lk,r(ξ), ηk∈[xk−1,xk],

where

 |lk,r(ξ)|≤Δxr.

Proof: For ,

 u(ξ)−sj,k(ξ)=u(ξ)−xk−ξΔxu(xk−1)−ξ−xk−1Δxu(xk).

Using Taylor expansions, we obtain

 u(ξ)−sj,k(ξ) = −1r!3∑r=2u(r)(ξ)lk,r(ξ)−14!u(4)(ηk)lk,r(ξ),

where are functions which depend on and , given by

 lk,r(ξ) = xk−ξΔx(xk−ξ−Δx)r−ξ−xk+ΔxΔx(xk−ξ)r (24) = (xk−ξ)r+p−1∑r=0(rp)(xk−ξ)p+1(−1)r−pΔxr−p−1. (25)

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

 ∂αu∂xα(xj)−δαuΔxα(xj)=ϵ1(xj)+ϵ2(xj),

where

 |ϵ1(xj)|≤C1Δx2|ϵ2(xj)|≤C2Δx2,

and and are independent of .

Proof: It is straightforward to prove that we have

 ∂αu∂xα(xj) = 1Γ(2−α)∂2∂x2Iα(xj) = 1Γ(2−α)1Δx2[Iα(xj−1)−2Iα(xj)+Iα(xj+1)]+ϵ1(xj),

where .

Let us define the error , such that,

 Iα(xj−1)−2Iα(xj)+Iα(xj+1)=Iα(xj−1)−2Iα(xj)+Iα(xj+1)+ES(xj).

We have

 ∂αu∂xα(xj) = 1Γ(2−α)1Δx2[Iα(xj−1)−2Iα(xj)+Iα(xj+1)] +1Γ(2−α)1Δx2ES(xj)+ϵ1(xj),

that is

 ∂αu∂xα(xj)=δαuΔxα(xj)+ϵ1(xj)+ϵ2(xj),

where

 ϵ2(xj)=1Γ(2−α)1Δx2ES(xj).

We are now going to compute the error . We have

 ES(xj) = j−1∑k=−∞ ∫xkxk−1(u(ξ)−sj−1,k(ξ))(xj−1−ξ)1−αdξ −2j∑k=−∞ ∫xkxk−1(u(ξ)−sj,k(ξ))(xj−ξ)1−αdξ +j+1∑k=−∞ ∫xkxk−1(u(ξ)−sj+1,k(ξ))(xj+1−ξ)1−αdξ.

Taking in consideration the previous lemma, let us denote

 ES(xj)=−4∑r=21r!Er(xj), (26)

where are defined as follows. For and ,

 Er(xj) = j−1∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ)(xj−1−ξ)1−αdξ (27) −2j∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ)(xj−ξ)1−αdξ +j+1∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ)(xj+1−ξ)1−αdξ,

and for

 Er(xj) = j−1∑k=−∞ u(4)(ηk)∫xkxk−1lk,r(ξ)(xj−1−ξ)1−αdξ (28) −2j∑k=−∞u(4)(ηk) ∫xkxk−1lk,r(ξ)(xj−ξ)1−αdξ +j+1∑k=−∞u(4)(ηk) ∫xkxk−1lk,r(ξ)(xj+1−ξ)1−αdξ.

For by changing variables, we obtain

 Er(xj) = j∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ−Δx)(xj−ξ)1−αdξ −2j∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ)(xj−ξ)1−αdξ +j∑k=−∞ ∫xkxk−1lk,r(ξ)u(r)(ξ+Δx)(xj−ξ)1−αdξ,

that is,

 Er(xj) = j∑k=−∞ ∫xkxk−1lk,r(ξ)[u(r)(ξ+Δx)−2u(r)(ξ)+u(r)(ξ−Δx)](xj−ξ)1−αdξ.

Let such that , for . For we have

 E2(xj) = j∑k=−∞ ∫xkxk−1lk,2(ξ)[u(r)(ξ+Δx)−2u(r)(ξ)+u(r)(ξ−Δx)](xj−ξ)1−αdξ = Δx22j∑k=Na+1u(4)(ξk)cj,k,2,ξk∈[xk−1,xk]

where

 cj,k,2=∫xkxk−1lk,r(ξ)(xj−ξ)1−αdξ

Since, by Lemma 1,

 |cj,k,2|≤Δx2∫xkxk−1(xj−ξ)1−αdξ

and

 ∫xjxa(xj−ξ)1−αdξ=12−α(xj−xa)2−α

we have

 |E2(xj)|≤Δx42(2−α)||u(4)||∞(xj−xa)2−α. (29)

For ,

 E3(xj)=j∑k=Na+1Δx(u(4)(ξk1)−u(4)(ξk2))cj,k,3, ξk1,ξk2∈[xk−1,xk]

and

 |cj,k,3|≤Δx3∫xkxk−1(xj−ξ)1−αdξ.

We have

 |E3(xj)|≤2Δx4(2−α)||u(4)||∞(xj−xa)2−α. (30)

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

 j−1∑k=Na+1u(4)(ηk)∫xkxk−1lk,4(ξ)(xj−1−ξ)1−αdξ ≤ Δx4||u(4)||∞j−1∑k=Na+1∫xkxk−1(xj−1−ξ)1−αdξ = Δx42−α||u(4)||∞(xj−1−xa)2−α.

Therefore, since for , we have

 j−1∑k=Na+1u(4)(ηk)∫xkxk−1lk,4(ξ)(xj−1−ξ)1−αdξ≤Δx42−α||u(4)||∞((xj−xa)2−α+Δx2−α).

Similarly, for the second integral we have

 j∑k=Na+1u(4)(ηk)∫xkxk−1lk,4(ξ)(xj−ξ)1−αdξ≤Δx42−α||u(4)||∞(xj−xa)2−α

and for the third integral

 j+1∑k=Na+1u(4)(ηk)∫xkxk−1lk,4(ξ)(xj+1−ξ)1−αdξ≤Δx42−α||u(4)||∞((xj−xa)2−α+Δx2−α).

Finally, we have

 |E4(xj)|≤3Δx42−α||u(4)||∞(xj−xa)2−α+2Δx6−α2−α||u(4)||∞. (31)

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

 Tnj=un+1j−unjΔt−djΔxα(τδαun+1j+(1−τ)δαunj)−pn+τj.

We have that

 un+1j−unjΔt=∂u(xj,tn)∂t+Δt2∂2u(xj,tn)∂t2+O(Δt2), (32)

and using the previous theorem we have

 Tnj = ∂u(xj,tn)∂t+Δt2∂2u(xj,tn)∂t2+O(Δt2)−τ(dj∂αu(xj,tn+1)∂xα+O(Δx2)) −(1−τ)(dj∂αu(xj,tn)∂xα+O(Δx2))−pn+τj.

Therefore

 Tnj = ∂u(xj,tn)∂t+Δt2∂2u(xj,tn)∂t2+−(1−τ)∂u(xj,tn)∂t−τ∂u(xj,tn+1)∂t +O(Δt2)+O(Δx2)

Finally,

 Tnj = (12−τ)Δt∂2u(xj,tn)∂t2+O(Δt2)+O(Δx2).

### 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

 Enj=Unj−unj (33)

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 (