Analytical and Numerical Solutions of Riesz Space Fractional Advection-Dispersion Equations with Delay

Mahdi Saedshoar Heris Mohammad Javidi and Bashir Ahmad

Abstract: In this paper,we propose numerical methods for the Riesz space fractional advection-dispersion equations with delay (RFADED).We utilize the fractional backward differential formulas method of second order (FBDF2)and weighted shifted Grünwald difference (WSGD)operators to approximate the Riesz fractional derivative and present the finite difference method for the RFADED.Firstly,the FBDF2 and the shifted Grünwald methods are introduced.Secondly,based on the FBDF2 method and the WSGD operators,the finite difference method is applied to the problem.We also show that our numerical schemes are conditionally stable and convergent with the accuracy of O(κ + h2)and O(κ2 + h2)respectively.Thirdly we find the analytical solution for RFDED in terms Mittag-Leffler type functions.Finally,some numerical examples are given to show the efficacy of the numerical methods and the results are found to be in complete agreement with the analytical solution.

Keywords: Riesz fractional derivative,shifted Grünwald difference operators,fractional advection-dispersion equation,delay differential equations,FBDF method.

1 Introduction

Fractional calculus finds its applications in diverse areas of science,engineering,economics and finance [Bagley and Calico (1991);Weaver Jr,Timoshenko and Young (1990);Marks and Hall (1981);Simo and Woafo (2016);Yang (2019);Cattani (2018);Yang,Abdel-Aty and Cattani (2019);Feng (2017)].In most of the cases,fractional differential equations (FDEs)cannot be solved exactly,so one needs to seek approximate and numerical techniques to solve these equations.Various numerical methods for solving FDEs are discussed in Galeone et al.[Galeone and Garrappa (2006);Garrappa (2015);Gorenflo (1997);Lubich (1986);Sulaiman,Yavuz,Bulut et al.(2019);Yang,Han,Li et al.(2016)].However,there are fewer works dealing with numerical methods for delay fractional differential equation [Heris and Javidi (2017b,a,2018a,2019);Javidi and Heris(2019)].In recent years,some authors investigated fractional order partial differential equations[Momani,Odibat and Erturk(2007);Momani and Odibat(2008a);Meerschaert and Tadjeran (2006)],and delay fractional partial differential equations [Zubik-Kowal(2000);Jackiewicz and Zubik-Kowal (2006);Tanthanuch (2012)].Such equations appear in mathematical modeling of several phenomena occurring in biology,medicine,population ecology,control systems and climate models[Wu(2012)].For details on numerical method for fractional partial differential equations,for instance,see Tadjeran et al.[Tadjeran,Meerschaert and Scheffler(2006);Liu,Zeng and Li(2015);Ding and Li(2013)].

The fractional advection-dispersion equation (FADE)is an important tool of groundwater hydrology to deal with the transport of passive tracers carried by fluid flow in a porous medium [Momani and Odibat (2008b);Liu,Anh,Turner et al.(2003);Huang and Liu(2005)].Meerschaert et al.[Meerschaert and Tadjeran (2004)] presented numerical methods for solving one-dimensional fractional advection-diffusion equation involving a Riemann-Liouville fractional derivative on a finite domain.Liu et al.[Liu,Anh and Turner(2004)]transformed the space fractional advection-diffusion equation into a system of ordinary differential equations and solved the resulting equations by using backward differentiation formulas.Shen et al.[Shen,Liu and Anh(2008)]discussed the fundamental solution and discrete random walk model for a time space fractional diffusion equation of distributed order.Numerical approximations and solution techniques for the space time Riesz-Caputo fractional advection-diffusion equation were studied in Shen et al.[Shen,Liu and Anh(2011)].For more details and examples,we refer the reader to the articles[Ding,Li and Chen (2015);Sousa (2012);Yang,Liu and Turner (2010);Wu,Baleanu and Xie(2016);Wu,Baleanu,Deng et al.(2015)].

In this paper,we focus on designing a numerical method for solving the following Riesz fractional advection-dispersion equation(RFADE)with time delay:

subject to the initial and boundary conditions:

where 0＜α ＜1,1＜β ≤2,Kα ≥0,Kβ ＞0 andτ ＞0.The Riesz space fractional operator on a finite domain[0,L]is defined as[Yang,Liu and Turner(2010)]

where

Partial differential equations with delay are more complicated and less studied in the literature.The solutions of PDEs are different from the ones for PDEs with delay.If we take exact solution of our model as initial delay solution,then the solution of our model without delay shifts to the right as

Thus,we can change initial delay solution to control dynamic properties of solutions.

Partial differential equations with delay have recently been studied by many authors and important aspects such as,existence and stability of solutions for these equations,are presented[Wu(2012)].In general,the exact solution of these equations cannot be obtained.Moreover,it is difficult to study the long-term dynamic properties of these equations.So one resorts to numerical simulation of such equations.In order to investigate the long-term dynamic properties,partial differential equations with delay are considered.As a typical example in the delay field with derivatives of integer order,we takeKα=0,β=2 in Eq.(1).This equation was considered as the reaction-diffusion equation with delay in Huang et al.[Huang and Vandewalle (2012)].For the delay field with fractional derivatives,we takeKα= 0 in Eq.(1)and the resulting equation is known as the Riesz space fractional diffusion equation with delay and nonlinear source term,see Yang[Yang(2018)].

We first obtain analytical solution for the problem at hand.Then we apply FBDF2 method for 0＜α ＜1 and shifted Grünwald difference operators for 1＜β ≤2 to approximate the Riesz space fractional derivative.Furthermore,we propose the finite difference method for the RFADED.We also show that the schemes for allhsmaller than

are stable and convergent with the accuracy of O(κ+h2)and O(κ2+h2)respectively.

The paper is organized as follows.Section 2 contains some definitions.In Section 3,analytical solution for the given problem is obtained.FBDF2 method is presented in Section 4,while Section 5 deals with Shifted Grünwald method.In Section 6,numerical methods for the RFADED are presented.The paper concludes with numerical simulations and results.

2 Preliminaries

In this section,we recall some definitions related to our work.Let C(J,R)denote the Banach space of all continuous functions from J=[0,T]into R endowed with the norm

and Cn(J,R)denotes the class of all real valued functions possessing derivatives upto ordernon J = [0,T],T＞0.

Definition 2.1[Kilbas,Srivastava and Trujillo (2006)].The fractional integral of orderα＞0 of the functionf∈C(J,R)is defined as

Definition 2.2[Kilbas,Srivastava and Trujillo (2006)].The Riemann-Liouville fractional derivative of orderα＞0 of the functionf∈Cn(J,R)is defined as

Definition 2.3[Kilbas,Srivastava and Trujillo (2006)].The Caputo fractional derivative of orderα＞0 of the functionf∈Cn(J,R)is defined as

Definition 2.4[ermk,Hornek and Kisela (2016)].The generalized delay exponential function (of Mittag-Leffler type)is defined by

whereλ ∈C,α,β,τ ∈R andm ∈Z andH(z)is the Heaviside step function.Ifλ ∈C,α,β,τ ∈R andm ∈Z,then the Laplace transform ofis

3 Analytical solution of problem

In this section,we derive the analytical solution of the RFADED(1-2).

Lemma 3.1.[Ilic,Liu,Turner et al.(2005)] Suppose that the one-dimensional Laplacian(-Δ)supplemented with Dirichlet boundary condition atx=0 andx=Lhas a complete set of orthonormal eigenfunctionsφnassociated with eigenvaluesλ2non a boundary region Ω = [0,L],that is,(-Δ)φn=λ2nφnleads to the eigenvaluesand the corresponding eigenfunctions

In the initial and boundary conditions (2),it is assumed thatμ1(t)andμ2(t)are nonzero smooth functions with first order continuous derivatives.We firstly transform the nonhomogeneous condition into a homogeneous one.Let

where

Substituting(11)into(1)leads to the the following problem with homogeneous boundary conditions satisfied byW(x,t):

where

Assume that the solution of(13)has the form:

Substituting(15)into(13),we obtain the Sturm-Liouville problem:

and

whereλ ＞0 andν(x,t)is the initial function arising from (2).By Lemma 3.1,the Sturm-Liouville problem(16)has eigenvalues and corresponding eigenfunctions

Therefore we set

Inserting(19)into(13)leads to

where

Taking Laplace transform of(20),we get

where

Writing

we obtain

Similarly

Let

with

Then

Thus we have

which yields the analytical solution of(1)-(2)given by

4 Fractional backward Differential Formulas of second order(FBDF2 method)

We consider the initial value problem

wherefis a sufficiently smooth function.We now introduce the FBDF method of second order(FBDF2)for(32)[Heris and Javidi(2018b)].For 0＜α ＜1,we have

where

with the coefficients given by

Lemma 4.1.[Heris and Javidi(2018b)]For 0＜α ＜1,the coefficients?j(α)satisfy

5 Shifted Grünwald method

The standard Grünwald-Letnikov difference formula was often unstable for 1＜β ≤2 and for time dependent problems [Meerschaert and Tadjeran (2004)].On the other hand,the Shifted Grünwald difference operators formula is stable and is given by

Observe that

wherep,q ∈Zand

Lemma 5.1.[Tian,Zhou and Deng(2015)]For 1＜γ ≤2,the coefficientssatisfy

Theorem 5.1.[Tian,Zhou and Deng (2015)] Let 1＜γ ≤2 andu ∈L1(R),-∞Dxγu,xD+γ∞uand their Fourier transforms belong toL1(R)and let the weighted and shifted Grünwald difference operators be defnied by

Then

wherepandqare integers and

Remark 5.1.In case of a well defined functionu(x)on the bounded interval [a,b] withu(a)=0 oru(b)=0,it can be extended to zero forx ＜aorx ＞b.Then the left and right Riemann-Liouville fractional derivatives ofuat each pointx ∈(a,b)can be approximated by the WSGD operators with second order accuracy as

where

Remark 5.2.For 1＜γ ≤2 and (p,q)= (1,0),Eq.(42)on the domain [0,L] can be written as

where

Lemma 5.2.[Tian,Zhou and Deng(2015)]For 1＜γ ≤2,the coefficientssatisfy

6 Numerical methods

In this section,we consider two cases.In case 1,we approximate the time derivative with order one and the Riesz space fractional derivative with order two.In case 2,we approximate the Riesz space fractional derivative and derive the Crank-Nicolson scheme for the equation.We partition the interval[0,L]into an uniform mesh with the space step sizeh=L/Mand the time step sizet=T/N,whereM,Nare two positive integers.The set of grid points are denoted byxi=ihandtj=jκfori=1,...,Mandj=n+1,···.

Case 1.We apply numerical method to Eq.(1)as follows.

Letu(xi,tj)=uji,f(xi,tj)=fji.Then

wherei= 1,...,M,j=n+ 1,···andR1is a truncation error term.Takingλα=cακKαh-αandλβ=cβκKβh-β,we have

Introducing

and

Eq.(47)takes the form:

where

Case 2.Here we approximate the Riesz space fractional derivative and derive the Crank-Nicolson scheme for the Eq.(1).Lettingwe have

wherei=1,...,M,j=n+1,...andR2is a truncation error term.Fixingand,we write

whereAandBare defined by(48).Thus Eq.(51)simplifies to the following form:

where

and

6.1 Stability of methods

Lemma 6.1.[Thomas (2013)] LetAbe a positive definite matrix of orderm-1.Then,for any parameterν ≥0,the following inequalities hold:

Theorem 6.1.ForDis a strictly diagonally dominant matrix.

Proof.We have

whereλα ＞0 for 0＜α ＜1 andλβ ＜0 for 1＜β ＜2.According to Lemma 4.1,whenk ≥4 and according to Lemma 5.2,whenk ≥4.Fork= 2,3,ifthenλα＜0 andλα＜0.Therefore,Di,j ＜0 whenj ＞i+1 orj ＜i-1.According to the Lemmas 4.1 and 5.2,we havewhich leads to

ThereforeDi,i ＞0.ForDi,i+1andDi,i-1,we have

Sinceλα ＞0 andλβ ＜0,therefore

For a giveni,we can write

Therefore

Therefore,Di,j ＜0 whenj ＞i+1 orj ＜i-1.Then relation(54)is valid for＜α ＜1,1＜β ≤2.Thus the matrixDis strictly diagonally dominant matrix.

Lemma 6.2.For,the matrixDis symmetric and positive definite.

Proof.In view of Eq.(49),the matrixDis clearly symmetric.Letυ0be an eigenvalue of the matrixD.Then it follows by the Ger?gorin circles Theorem[Varga(2010)]that

or

Then,by Theorem 6.1,we have

which shows thatDis positive definite.

Remark 6.1.Foris strictly diagonally dominant and symmetric positive definite.

Theorem 6.2.The first numerical method(50)foris stable.

Proof.LetUjbe the numerical solution andujbe an exact solution.Since the matrix(I+D)is invertible,we have

where

By Lemma 6.1,we have

Thus the numerical method(50)is stable.

Theorem 6.3.For,the second numerical method(53)is stable.

Proof.LetUjbe the numerical solution andujbe an exact solution.Since the matrix(I+)is invertible,we have

where

Then

which,by Lemma 6.1,yields

This shows that numerical method(53)is stable.

6.2 Convergence of methods

This subsection is concerned with convergence of the methods presented in Section 5.For the first method,we can write

Thus the local truncation error of(47)will be of the form:

Theorem 6.4.LetUjbe the numerical solution andujbe an exact solution of(50).Then for,we have

whereCis a positive constant.

Proof.It is easy to find that

and

Lettingeji=Uij -ujiand using(59)and(60),we obtain

which can equivalently be written in matrix-vector form as

where

Writing

we get

which,on iterating and using the given initial condition,yields

By Lemma 6.1,we can write

Therefore we have

For the second method,we can write

Thus the local truncation error of(51)will be of the following form:

Theorem 6.5.LetUjbe the numerical solution andujbe an exact solution of(53).Then,for,we have

whereCis a positive constant.

Proof.Obviously

and

Let us seteji=Uij -ujiand use(64)and(65)to obtain

which can be expressed in matrix-vector form as

where

Let us take

so that

which,on iterating and using the given initial condition,becomes

On the other hand,it follows by Lemma 6.1 that

In consequence,we obtain

7 Test examples

In this section,we show the efficacy and accuracy of the proposed methods with the aid of examples.

Example 7.1.We consider the following RFDED

subject to the initial condition:

where 1＜β ≤2 and

The exact solution of the problem(67)-(68)isu(x,t)=x2(1-x)2e-t.Numerical results for the problem(67)-(68)are given by Tabs.1-2 and shown in Fig.1.

Table 1:The absolute errors and the convergence orders of the first method (50)for example 7.1(67-68)

Figure 1:The numerical approximation and exact solution by the second method for example 7.1(67-68)(RFDED),for β =1.8,when T=1

Figure 2:The numerical approximation and exact solution by the second method for example 7.2(69-70)(RFDED),for α=0.1 and β =1.8,when T=1

Example 7.2.Consider the RFDED

subject to the initial condition:

where 0＜α ＜1,1＜β ≤2 and

Table 2:The absolute errors and the convergence orders of the second method (53)for example 7.1(67-68)

Figure 3:The numerical approximation by the second method for example 7.2(69-70)(RFDED),for various α=0.1,0.5,0.9,0.99,when T=1 and β =1.01

Figure 4:The numerical approximation by the second method for example 7.2(69-70)(RFDED),for various β =1.2,1.5,1.9,2,when T=1 and α=0.99

Table 3:The absolute errors and the convergence orders of the second method (53)for example 7.2(69-70)

The exact solution of this problem isu(x,t)=cos(t)x2(1-x)2.Numerical results for the given problem are given by Tab.3 and Figs.2-4.In Tab.3,scheme(53)forh ＜0.0952 is stable and convergent forα=0.8,β=1.2.

8 Conclusion

In this paper,we applied the FBDF method of second order and the shifted Grünwald method for solving the Riesz space fractional advection-dispersion equations with delay.The FBDF2 and the shifted Grünwald methods are introduced.Furthermore we find the analytical solution for RFDED in terms of t Mittag-Leffler type functions.The approximation of solution for the Riesz space fractional advection-dispersion equations with delay relies on the FBDF2 method and the WSGD operators,and is obtained by applying the finite difference method.It is shown that the schemes forh ＜are stable and convergent with the accuracy of O(κ+h2)and O(κ2+h2)respectively.Numerical methods presented in this paper are illustrated with the help of the examples.The obtained results clearly demonstrate that our methods are efficient and produce accurate results.

Acknowledgement:The authors would like to express special thanks to the referees for carefully reading,constructive comments and valuable remarks which significantly improved the quality of this paper.