Finite Difference Approximation Method for a Space Fractional Convection–Diffusion Equation with Variable Coefficients

Abstract

Space non-integer order convection–diffusion descriptions are generalized form of integer order convection–diffusion problems expressing super diffusive and convective transport processes. In this article, we propose finite difference approximation for space fractional convection–diffusion model having space variable coefficients on the given bounded domain over time and space. It is shown that the Crank–Nicolson difference scheme based on the right shifted Grünwald–Letnikov difference formula is unconditionally stable and it is also of second order consistency both in temporal and spatial terms with extrapolation to the limit approach. Numerical experiments are tested to verify the efficiency of our theoretical analysis and confirm order of convergence.

1. Introduction

Fractional differential equations (FDE) have attracted the attention of many researchers and scientists due to their importance in different fields of study such as viscoelasticity, fluid mechanics, physics, biology, engineering, and flows in porous media (see [1,2,3,4,5,6] and the references cited therein). As different experiments and implementations have shown, non-integer space derivatives have been used to develop anomalous diffusion to which a particle spreads at a rate inconsistent with the integer Brownian motion problem in the direction of both time and space. When non-integer order is replaced by the second order derivative in a diffusion equation, it acts to enhance the process which we call super-diffusion [7,8,9,10,11,12]. Laboratory experiments and field-scale tracer dispersion breakthrough curves (BTCs) are suitable for exhibiting early time arrivals that are not captured by the integer order derivatives and these non-Fickian phenomena can be controlled by non-classical order convection–diffusion and dispersion equations (FCDE) as it was explained in [13]. To increase the number of applications, there should be significant interest in constructing numerical schemes to solve a well known space fractional convection–diffusion model that has space variable coefficients. In most cases, non-integer order differential problems have no exact solution, so various iterative and numerical approximations [3,9,14] must be pointed out in advance. In general, these kinds of approaches have become important in finding the approximate solutions of fractional differential equations, so extensive numerical methods have been developed for space fractional convection–diffusion equations such as the spectral method [15], finite volume method [16,17], finite difference method [2,9,14,18,19,20,21,22,23,24,25,26], finite element method [27,28,29,30] and collocation method [31,32].

When the discretization of domain over the region (which belongs to the geometry) is not complex, finite difference approximations are easier and faster than other methods (see [16,33] for further details) to get numerical solutions. In [34], the author used an unconditional stable difference method for time−space fractional convection–diffusion problems with space variable coefficients with first order convergence both in time and space. The Crank–Nicolson finite difference method for one-sided space fractional diffusion equations using an extrapolation method to get second order convergence was studied in [23]. In [9], the explicit and implicit finite difference methods are discussed for a one-sided space fractional convection–diffusion equation with first order convergence in both time and space. A first-order implicit finite difference discretization method for a two-sided space fractional diffusion equation (SFDE) is also applied in [10]. Recently, an unconditionally stable second order accurate difference method for a two-sided time−space fractional convection–diffusion equation was constructed in [35] using the weighted and Shifted Grünwald–Letnikov difference approximation. It is not suitable to apply the weighted combined with shifted Grünwald–Letnikov difference approximation for one-sided Riemann–Liouville fractional derivative to have second order accurate in space. To deal with such issues, it is important to develop a numerical scheme that leads to evaluate a one-sided space fractional convection–diffusion problem. Thus, the main focus of our study is to have temporal and spatial second order convergence estimates for one-sided space fractional convection–diffusion equations based on a stable finite difference method and using spatial extrapolation to the limit approach. The scheme has been treated using the Crank–Nicholson method with the novel Shifted Grünwald–Letnikov difference approximation and the algorithm has been examined both theoretically and experimentally.

Let us consider space-fractional convection–diffusion equation with variable coefficients:

∂u(x,t)∂t+c(x)∂u(x,t)∂x=d(x)^∂αu(x,t)∂^xα+p(x,t),x∈(L,R),t∈(0,T],α∈(1,2];

with the given initial condition:

u(x,0)=g(x),L≤x≤R,

and homogeneous Dirichlet boundary conditions:

u(L,t)=0,u(R,t)=0,0≤t≤T,

wherec(x),d(x)andg(x)are continuous functions on[L,R]andp(x,t)is continuous function on[L,R]×[0,T]. Hereu(x,t)is the concentration,d(x)>0is the variable diffusion coefficient,c(x)>0is the fluid variable velocity which means the system is evolving in space due to a velocity field andp(x,t)is sink term so that the fluid transport is from left to right. For the case of integer order (α=2 ), Equation (1) gives to the classical convection–diffusion equation(CDE). In this study, we have only considered the fractional derivative case which describes a physical meaning in [36] and it involves only a left-sided fractional order derivative. We have assumed that this one-dimensional space fractional convection–diffusion problem has sufficiently smooth and unique enough solutions.

The structure of this paper is arranged as follows. In Section 2, we introduce some preliminary remarks, lemmas and definitions and we show the formulation of the new Crank–Nicolson with right Shifted Grünwald–Letnikov difference scheme in Section 3. In Section 4, we describe the unconditional stability using Gerschgorin Theorem and convergence order analysis of the scheme. In Section 5, numerical tests are implemented to show the relevance of our theoretical study and the conclusions are put in Section 6.

2. Preliminary Remarks

Definition 1.

The Riemann fractional derivative operator_D∗αwith order α is written as:

_D∗αu(x)=1Γ(r−α)^drd^xr_∫Lxu(t)^{(x−t)α−r+1}dt,α>0

wherer−1<α<r,r∈N,t>0.

Definition 2.

The left hand side and the right hand side fractional order derivatives, respectively, in Equation (1) are the Riemann–Liouville fractional derivatives with order α which are given by:

(_D+αu)(x)=1Γ(r−α)^drd^xr_∫Lx ^{(x−s)r−α−1}u(s)ds(_D−αu)(x)=^(−1)rΓ(r−α)^drd^xr_∫xR ^{(s−x)r−α−1}u(s)ds

forr−1<α<r,x∈ℜ.

Definition 3

([3]). Let u be given on ℜ. The standard Grünwald–Letnikov estimate for1<α≤2with positive order α is defined by the formula,

^Dαu(x,t)≈1^hα∑k=0_Nxωk(α)u(x−kh,t),

we also define the Grünwald–Letnikov difference operator as:

^h−α(_Δhαu)(x,t)≈∑k=0_Nxωk(α)u(x−kh,t),h>0,x∈ℜ,

where

_ωk(α)=α(α−1)…(α−k+1)k!,

is called Grünwald–Letnikov coefficient which is the Taylor series expansionω(z)=^(1−z)αwhich is the generating function. We can expressed the coefficients by the following recursive relations.

_ω0(α)=1,_ωk(α)=(1−α+1k)_ωk−1(α),k=1,2,….

Lemma 1

([37]). Assume that1<α≤2, then Grünwald–Letnikov coefficients_ωk(α)satisfy:

_ω0(α)=1,_ω1(α)=−α<0,_ω2(α)=α(α−1)2>01≥_ω2(α)≥_ω3(α)≥…≥0,∑k=0∞_ωk(α)=0,∑k=0_Nxωk(α)<0,_Nx≥1.

The Shifted Grünwald–Letnikov difference operator expression is suitable for our purpose because, it allows us to estimate_D∗αu(x) , which is defined in Equation (2), numerically in an accurate way. According to [14], right shifted Grünwald–Letnikov difference operator with p shifts for^αthorder Left R-L fractional derivative ofu(x,t),x∈[L,R]atx=_xmcan be expressed as:

_D∗αu(x,t)≈1^hα∑k=0_xm−Lh+p_ωk(α)u(x−(k−p)h,t)

where

_xm=L+mh,h=R−L_Nx,m=0,1,2,…_Nx.

Lemma 2

([38,39]). Letu∈^C2n(ℜ)that has a finite degree of smoothness with(_D+αu)(x)which is approximated by^h−α_Δhαu(x)possesses an asymptotic expansion in integer powers of the step-length h, then an expansion in even powers of h for the Shifted operator can be written in the form:

_Δh,pαu(x)=∑j=0∞^−1j_jαux+αh2−jh,h>0.

Lemma 3

([39]). Letu∈^Cn+3(ℜ)all derivative of u up to the ordern+4belong to^L1(ℜ) . Then the Fourier transform of the Grünwald–Letnikov difference operator defined in Equation (5), is

ϕ^(x)=_∫ℜϕ(t)^eixtdt.

Theorem 1.

Letu∈^C2n+3(ℜ)with all derivatives of u up to order2n+3belong to^L1(ℜ). Forp≥0define the shifted Grünwald–Letnikov operator:

(_Δh,pα)u(x)=∑k=0∞_ωk(α)ux−k−ph,

with_ωk(α)=^(−1)2kα_a2k=α_a2k. Then, ifL=−∞ in Equation (2), for any computable coefficient_a2k, which is independent ofh,uand x, we have

^h−α_Δh,pαu(x)=_D+αu(x)+∑k=1n−1_b2kD+α+2ku(x)^h2k+O(^h2n)

uniformly inx∈ℜ.

Proof of Theorem 1.

We closely follow the result described in [9,10] for the unshifted Grünwald–Letnikov formula and also in [23] for the shifted Grünwald–Letnikov formula. We can see that with the Riemann-Lebesgue lemma, the assumptions on u indicates for real positive constant_C1and from the condition which is imposed on u, we have

u˜(t)≤_C1 ^{(1+|t|)−2n−3}.

From Lemma 3 for allt∈ℜthe Fourier transform foru(x)of the Grünwald–Letnikov approximation is

u˜(t)=_∫ℜu(x)^eixtdx.

From the definition of Fourier transform, we have observed that for a constanta∈ℜ, we have:

F([u(x−a)])(t)=^eiatu˜(t).

The function

^{1−e−zzα ezp}=_ωα,p(z),

have the Taylor expansion

_ωα,p(z)=∑k=0∞_a2k ^z2k,

where_a2k=^(−1)2kα_a2k=α_a2k, converges absolutely forz≤1since the function_ωα,p(z)is bounded on ℜ. The shifted Grünwald difference approximation(_Δh,p)u(x)∈^L1(ℜ).

Thus, we have

F(^h−α _Δh,pαu)(t)=^{h−α e−itph}∑k=0∞_2kα^eikthu˜(t)=^{h−α e−itph 1−eitphα}u˜(t)=^{(−it)α 1−eith−ithα e−itph}u˜(t)=^(−it)α _ωα,p(−ith)u˜(t)

since_ωα,p(−ith)is analytic around the origin, we express it as an even power expansions

_ωα,p(z)=∑k=0∞_a2k ^z2k

which absolutely convergent for all|z|≤R. For this a bounded function_ωα,p(z)on ℜ, there exist a real positive constant_C2which satisfy:

^{1−eix−ixα}−∑k=0n−1_a2k ^−ix2k≤_C2 ^x2n

is bounded uniformly inx∈ℜ. For any value|x|≤R, we have

(_ωα,p(−ix)−∑k=0n−1_a2k ^(−ix)2k=∑k=n∞_a2k ^(−ix)2k≤^x2n∑k=n∞_a2kα^x2(k−n)≤_C3 ^x2n

which is bounded on ℜ. For the other case|x|>Ralso, we have

_ωα,p−ix=^{1−eix−ixα eipx}≤^{2α Rα}<_C4 ^x2n

where_C4=^{2α Rα+2n}<∞and also

∑k=0n−1_a2k ^(−ix)2k≤^x2n∑k=0n−1|(_a2kα)|^x2(k−n)≤_C5 ^x2n

with_C5=_{∑k=0n−1a2kα}^R2k−2n<∞. Now, we set that

_C2=max∑k=0∞_a2k^R2k−2n,^{2α Rα+2n}+∑k=0n−1_a2k^R2k−2n

since

∑k=0∞_a2k^R2k−2n=∑k=0n−1_a2k^R2k−2n+∑k=n∞_a2k^R2k−2n

_C2=^{2α Rα+2n}+∑k=0n−1_a2k^R2k−2n

Then, this implies that Equation (15) holds for allx∈ℜ. From Equation (17), we can write

F(^h−α _Δh,pαu)(t)=∑k=0n−1_a2k ^{−itα+2k h2k}u˜(t)+φ˜(t,h)

where

φ˜(k,h)=^−itα_ωα,p−ith−∑k=0n−1_a2k ^−ith2ku˜(t)

since

^−itα+2ku˜(t)=_D+α+2ku˜(t).

Therefore, we have

^(−it)α+2ku˜(t)∈^L1(ℜ).

Moreover, we see that

φ˜(t,h)∈^L1(ℜ),

and with the conditions imposed on u, we can say that(1+^x2n+3u˜(t)is bounded on ℜ.

Thus,^t2α−3|u˜(t)|∈^L1(R). This implies that,

φ˜(t,h)≤C^{h2n 1+t2α−3}

fork∈ℜwithC=_{C1 C2}. Therefore using the Fourier inversion transform, we have

^h−α_Δh,pαu(x)=_D+αu(x)+∑k=1n−1_a2kD+α+2ku(x)^h2k+φ(x,h),

where

φ(x,h)=C_∫R ^e−itxφ˜(t,h)dt≤C_∫Rφ˜(t,h)dt≤C^h2n.

At last, we have

^h−α_Δh,pαu(x)=_D+αu(x)+∑k=1n−1_a2k(_D+α+2ku)(x)^h2k+O(^h2n).

□

Remark 1.

From Equation (10), it can be seen that forp=α/2, the error takes its minimum value and a second order convergence is achieved. We need the grid points_xm−(k−p)hto find an optimal positive integer p that makesp−α/2 is minimum. It is numerically proved in [3] that for the value0<α≤1,p=0is acceptable; while for1<α≤2,p=1 is optimal.

Remark 2.

Theorem 1 is the base of Extrapolation to the limit. Therefore one can apply it the Shifted Grünwald–Letnikov difference operator to obtain the convergence rate with arbitrary high order^hk,k=1,2,3,…,nsuch that

^h−α(^q−α _Δqh,pαu)(x)−q(_Δh,pαu)(x)1−q,0<q<1

(q is fixed) converges to(_D+αu)(x)+O(^h2).

3. Problem Formulation of the Scheme

Consider the following one-dimensional space fractional convection–diffusion problem:

∂u(x,t)∂t=−c(x)∂u(x,t)∂x+d(x)^∂αu(x,t)∂^xα+p(x,t),(x,t)∈(L,R)×(0,T]u(x,0)=g(x),x∈[L,R]u(L,t)=0,u(R,t)=0,t∈[0,T]

which is based on shifted Grünwald–Letnikov difference method with1<α≤2on a finite domainL<x<R.

Crank–Nicolson Scheme for Time and Shifted Grünwald Difference Scheme for Space Discretization

We partition the finite interval[L,R]with a uniform mesh in the space size steph=(R−L)/_Nxand the time stepτ=T/_Nt, in which_Nx,_Ntare non-negative integers and the set of grid size points is symbolized by_xm=mhand_tn=nτfor0≤m≤_Nx,0≤n≤_Nt. Set_tn+1/2=(_tn+1+_tn)/2with0≤n≤_Nt−1.

We use the following notations:

_umn=u(_xm,_tn),_pmn+1/2=p(_xm,_tn+1/2),_{δt umn}=_umn+1−_umnτ,_cm=c(_xm),_dm=d(_xm).

Applying the C-N technique for the time discretization of Equation (20) gives to

_{δt umn}=−_cm4h_um+1n+1−_um−1n+1+_um+1n−_um−1n+_dm2^hα∑z=01∑k=0_Nx−1_{ωk(α)um−k+1n+z}=_pmn+1/2+O(^τ2).

In space discretization we have used the central finite difference method for the convection term and the Shifted Grünwald–Letnikov operator for the space fractional derivative with the approach of spatial Extrapolation to the limit, respectively.

See the full discretization of the scheme:

_umn+1−_umnτ=−_cmum+1n−_um−1n+_um+1n+1−_um−1n+14h+_dm2^hα∑z=01∑k=0m+1_{ωk(α) um−k+1n+z}+_pmn+_pmn+12.

Multiplying Equation (22) byτthe discretization equation, we have

_umn+1−_umn=−_cmτ4h(_um+1n−_um−1n+_um+1n+1−_um−1n+1)+_dmτ2^hα∑z=01∑k=0m+1_{ωk(α) um−k+1n+z}+τ_pmn+1/2

The above equation is used to predict the values ofu(x,t)at timen+1, so all the values of u at time n are assumed to be known. For simplification

_μm=_cmτh,_ηm=_dmτ^hα, then we have

1−_ηm2_ω1αumn+1+−_μm2−_ηm2_{ω2αum−1n+1}+−_μm2−_ηm2_ω0αum+1n+1−_ηm2∑k=3m+1_{ωkα um−k+1n+1}=1+_ηm2_ω1αumn+_μm2+_ηm2_ω2αum−1n+_ηm2_ω0α+_μm2_um+1n+_ηm2∑k=3m+1_{ωk(α) um−k+1n}+τ_pmn+12.

Both the convection and diffusion variable coefficients are(_Nx−1)×(_Nx−1)diagonal matrices which are defined by

_μm=τ2hdiag_C1,_C2,_C3,…_CNx−1,_ηm=τ^hαdiag_d1,_d2,_d3,…_dNx−1.

These discretization together with Dirichlet boundary conditions which results in a linear system of equations for which the coefficient matrix is the sum of lower triangular and upper-diagonal matrices. The above discretization can be re-arranged to yield:

1−_ηm2_ω1αumn+1+(−_μm2−_ηm2_ω2α)_um−1n+1+−_μm2−_ηm2_ω0α)_um+1n+1−_ηm2(∑k=3m+1_{ωkα um−k+1n+1}=(1+_ηm2_ω1α)_umn+_μm2+_ηm2_ω2αum−1n+(_ηm2_ω0α+_μm2)_um+1n+_ηm2(∑k=3m+1_{ωkα um−k+1n})+τ(_Pmn+12).

Denoting_Umnas the numerical approximation of_umn , we can construct the C-N scheme for Equation (20)

1−_ηm2_ω1αUmn+1+−_μm2−_ηm2_{ω2αUm−1n+1}+−_μm2−_ηm2_ω0αUm+1n+1−_ηm2∑k=3m+1_{ωkα Um−k+1n+1}=1+_ηm2_ω1αUmn+_μm2+_ηm2_ω2αUm−1n+_ηm2_ω0α+_μm2_Um+1n+_ηm2∑k=3m+1_{ωkα Um−k+1n}+τ(_Pmn+12).

I is the_Nx−1×_Nt−1identity matrix with_Am,nas the matrix coefficients. These coefficients, form=1,2,3,…,_Nx−1,n=1,2,…,_Nt−1are given by:

_{Am ,n}=0,n≥m+2−_μm2−_ηm2_ω0(α),n=m+1(1−_ηm2_ω1(α)),n=m(−_ηm2_ω2(α)−_μm2),n=m−1−_ηm2_{ωm−n+1(α)}n≤m−1.

The finite difference scheme (24) and (26) defines a linear system of equations as

(I+A)^Un+1=(I−A)^Un+τ(_pmn+12)^Un+1=^{[_u1n+1,_u2n+1,…,_uNx−1n+1]⊤Un}+τ_Pmn+12=^{[0,τ_p1n+12,τ_p2n+12,…,τ_pNx−1n+12+(_ηNx−12+_μNx−12),0]⊤}.

Theorem 2.

Suppose that1<α≤2, the coefficient matrix defined in Equations (24)–(27), then the diagonal matrix and the coefficient matrix satisfy:

_Am,m>∑n=0,m≠1_Nx−1|_Am,n|,m=1,2,3,…,_Nx−1.

Proof of Theorem 2.As we have seen from the coefficient matrix defined in Equation (27),

_Am,m+1=_μm2−_ηm2_ω0(α)=_μm2−_ηm2<0

_Am,m−1=−_ηm2_ω2(α)−_μm2=−_ηm2(^α2−α2), but from Lemma 1,^α2−α2>0for1<α≤2mean that−_ηm2(^α2−α2)<0.

Whenn<m−1, we have,−_ηm2_{ωm−n+1(α)}<0and whenn=m,_Am,m=1−_ηm2_ω1(α)=1+_ηm2α>0. This implies that_{∑n=0,m≠1Nx−1}|_Am,n|<_Am,m.

Therefore, the diagonal matrix is strictly dominant. □

4. Theoretical Analysis of Finite Difference Scheme In general for analyzing convergence and stability, we consider the following description.

Let_χh=ν:ν=_νm:_xm=mhm=0Nx ,_ν0=_νNx =0be the grid function.

For anyν=_νm∈_χh, we define our point-wise maximum norm as

_||ν||∞=ma_x1≤m≤Nx|_νm|,

and the discrete^L2-norm

ν=h∑m=1_Nx−1_νm2.

4.1. Boundedness of the Fractional Scheme The Classical Crank–Nicolson scheme combines the stability of an implicit finite difference method with its accuracy which produce second order convergence in both space and time.

Theorem 3.

Crank–Nicolson scheme for solving space fractional convection–diffusion equations given by the following problem:

∂u(x,t)∂t+c(x)∂u(x,t)∂x=d(x)^∂αu(x,t)∂^xα+p(x,t).

which is based on shifted Grünwald–Letnikov difference approximation scheme is bounded for1<α≤2.

Proof of Theorem 3.

Consider C-N scheme for the space-fractional convection–diffusion problem for1<α≤2

_umn+1−_umnτ=−_cmumn−_um−1n+_umn+1−_um−1n+12h+_dm2^hα∑j=01∑k=0m+1_{ωk(α) um−k+1n+j}+_pmn+_pmn+12.

Here, we have shown the convergence and boundedness of the scheme by taking the smaller time-step in terms of Lax–Richtmyer stability analysis that uses a weaker bound (see [40]). Our matrix A has an eigenvalues ofλthat have positive real parts, and, we also have found a strictly dominant matrix. These eigenvalues which are centered in the disks at each diagonal entries as:

_Am,m=(1−_ηm2_ω1α)=1+α_ηm2.

with_μm=_cmτh,_ηm=τ_dm^hα . From the Gerschgorin Theorem in [41], the radius of the matrix can be expressed as

_{∑n=0,m≠1NxAm,n22}=^{−_ηm2−_μm2∑n=0m+1_{ωm−n+1(α)}2}≤^{−_ηm2−_μm22 ∑n=0m+1_{ωm−n+1(α)}2}.

Since from the Grünwald coefficients we have_{ωm−n+1(α)}≤_ω1(α)and_ω1(α)=−α, we have that:

_{∑n=0,m≠1Nx(Am,n)22}≤_{∑n=0,m≠1Nx(Am,m)22}≤_Am,m22≤_{(−ηm2−μm2)22 ω1(α)22}≤_1+ηm2α22.

For a bounded ratio of time-stepτand space-step h withnτ≤T, we have

_^Am,mn2≤^1+_ηm2αn/2.

From the relation of Parseval’s Theorem, [40]

_Am,m2≤^1+_ηm2αn/2≤^eαT/2.

which shows that the scheme is bounded. □

4.2. Stability Analysis

Theorem 4.

Let_Umnbe the numerical approximation of the exact solution_umn , then the C-N finite difference scheme (28) is unconditionally stable.

Proof of Theorem 4.

Consider the matrix coefficient of the difference approximation for the problem (20) can be written as described above

(I+A)^Un+1=(I−A)^Un+τ_pmn+1/2.

Let^en=_e1n,_e2n,_e3n,…,_eNx−1n, and take the relation between the error^en+1in^Un+1and the error^enin^Unwhich is given by the linear system

^en+1=^(I+A)−1(I−A)^en.

First of all, we must show that the (non-real valued) eigenvalues of the coefficient matrices A have positive real parts. For_ω1(α)=−αwith fractional order1<α<2andk≠1; we have_ωk(α)>0. In addition to this,−_ω1α=α≥_{∑k=0,k≠1N ωkα}for the valueN>1 . As stated in Gerschgorin Theorem ([41], pp. 136–139), the eigenvalues of the given matrix A are inside the disks centered at each diagonal entry.

_Am,m=(1−_ηm2_ω1(α))=1+_ηm2α>0,

with radius

_rm=∑n=0,m≠1_Nx|_Am,n|=_ηm2∑n=0m+1_{ωm−n+1(α)}<(1+_ηm2).

These Gerschgorin disks are belong to the right half of the complex plane. Thus, the eigenvalue of the coefficient matrix A has positive real part which implies that A has an eigenvalueλif and only if(I−A)has an eigenvalue(1−λ)if and only if^(I+A)−1(I−A)has an eigenvalue1−λ1+λ. From the first part of this sentence, we have seen that all the eigenvalues of the matrix given by(I+A)have a radius larger than unity which implies the matrix is invertible. Now we can see from the above description the real part ofλis non-negative which we can conclude that(1−λ)(1+λ)<1.

Thus, the spectral radius of the system matrix^(I+A)−1(I−A)is strictly less than unity which implies that the difference scheme is unconditionally stable. □

4.3. Convergence Analysis

First of all we have given the Truncation error of the C-N scheme. It is obvious to conclude that:

u(_xm,_tn+1)−u(_xm,_tn)τ=^{∂u(x,t)∂tn+1/2}+O(^τ2)._{c(x)∂u(x,t)∂x+d(x)^∂αu(x,t)∂^xαmn+1/2}=12_cm∂u(_xm,_tn+1)∂x+_dm^∂αu(_xm,_tn+1∂^xα+12_cm∂u(_xm,_tn)∂x+_dm^∂αu(_xm,_tn∂^xα+O(^τ2).

c(_xm)∂u(x,t)∂x≈u(_xm+1,_tn+1)−u(_xm−1,_tn+1)2h+O(^h2).

From the above Extrapolation to the limit Theorem forn=1, we got

^∂αu(x,t)∂^xα≈∑k=0m+1_{gk(α) um−k+1}+O(^h2).

Therefore the local truncation error of (20) is given by_Tmn+1=O(^τ2+τh)

Theorem 5.

Let_umn be the exact solution of problem (20), and_Umnbe the solution of the finite difference scheme (26), then for all1≤n≤_Nt , we have the estimate

_umn−Umn∞≤c(^τ2+h)

where_umn−Umn∞=ma_{x1≤m≤Nxumn}−_Umn, c is a non-negative constant independent of h and τ with||.||stands for the discrete^L2-norm.

Proof of Theorem 5.

Denote^en=_umn−_Umnwhere^en=(_e1n,_e2n,…,_eNx−1n). We have^e0=0, we have from Equations (26) and (27) ifn=0,

_Rm1=−_μm2−_ηm2_{ω0(α)em−11}+1+_ηm2α_e1m+−_μm2−_ηm2_ω2(α)em+11−_ηm2∑k=3_{Nxωk(α) em−n+11}.

ifn>0,

_Rmn+1=−_μm2−_ηm2_{ω0(α)em−1n+1}+1+_ηm2α_en+1m+−_μm2−_ηm2_{ω2(α)em+1n+1}−_ηm2∑k=3_{Nxωk(α) em−n+1n+1}.

where_Rmn+1≤c(^τ2+h),m=1,2,…,_Nx−1,n=1,2,3,…,_Nt−1, c is non-negative constant independent of h andτ.

We can use the mathematical induction to prove the Theorem. Letn=1and assume|_ej|=ma_{x1≤m≤Nx−1}|_em1|, we have the following expression.

||^e1 _||∞=_ej1≤−_μj2−_ηj2_{ω0(α)ej−11}+1+_ηj2α_e1j+−_μj2−_ηj2_ω2(α)ej+11−_ηj2∑k=3_{Nxωk(α)ej−n+11}≤−_μj2−_ηj2_{ω0(α)ej−11}+1+_ηj2α_e1j+−_μj2−_ηj2_ω2(α)ej+11−_ηj2∑k=3_{Nxωk(α) ej−n+11}=_Rj1≤c(^τ2+h).

Suppose that ifn≤r,||^er _||∞≤c(^τ2+^h2)hold and assumen=r+1, let_ejr+1=ma_{x1≤m≤Nx−1emr+1}, notice that from Lemma 1, we have_{∑k=0Nx ωk(α)}<0,m=1,2,…,_Nx. Therefore,

_^er+1∞=_ejr+1≤−_μj2−_ηj2_{ω0(α)ej−1r+1}+1+_ηj2α_er+1j+−_μj2−_ηj2_{ω2(α)ej+1r+1}−_ηj2∑k=3_{Nxωk(α)ej−n+1r+1}≤−_μj2−_ηj2_{ω0(α)ej−1r+1}+1+_ηj2α_er+1j+−_μj2−_ηj2_{ω2(α)ej+1r+1}−_ηj2∑k=3_{Nxωk(α) ej−n+1r+1}=_Rjr+1≤c(^τ2+h)

which completes the proof. □

Remark 3.The Crank–Nicolson scheme, for classical convection–diffusion equation, provides stable C-N finite difference method that is second order convergence in time and space. Also a study based on C-N finite difference method with the spatial extrapolation to the limit method, see Theorem 1, is used to get temporal and spatial second order for one-sided SFCDEs with space variable coefficients.

5. Numerical Tests Problem test 1

1. Consider the space-fractional diffusion type of problem:

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

with initial condition

u(x,0)=(^x2−^x3);0≤x≤1

homogeneous Dirichlet boundary condition

u(0,t)=0=u(1,t)

with variable diffusion coefficient,

d(x)=Γ(1.2)^xα,

and source term

p(x,t)=(6^x3−3^x2)^e−t

The exact solution is

u(x,t)=(^x2−^x3)^e−t

All numerical experiments are implemented using Theorem 1 and C-N scheme with the space domain,0<x<1and time domain,0<t<T . Figure 1 shows the maximum error produced by C-N scheme for large enough time domain and numerical solution is close enough to the exact solution using C-N scheme withα=1.5 in Figure 2. The maximum error and second order convergence for the fractional diffusion and fractional convection–diffusion equation with variable coefficients are given in Table 1, Table 2 and Table 3.

Problem test 2

2. Consider the space-fractional convection–diffusion type of equation with variable coefficients:

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

with initial condition

u(x,0)=(^xα−x);0≤x≤1

homogeneous Dirichlet boundary condition

u(0,t)=0=u(1,t)

with variable convection–diffusion coefficients respectively,

c(x)=^x15,d(x)=^x1100,

and source term

p(x,t)=^e−2t(2(x−^xα)−Γ(α)+Γ(α+1)Γ(α)^xα−1−1)

The exact solution is

u(x,t)=^e−2t(^xα−x)

Figure 3 and Figure 4 show the numerical and exact solutions for fractional diffusion and fractional convection–diffusion problems with large enough time domain in example 1 and 2, respectively.The exact and numerical solution of fractional convection–diffusion equation by C-N scheme is also given in Figure 5. In Table 4, the maximum error and first order convergence in space is obtained using C-N scheme without extrapolation to the limit approach by fixing the time step.

Problem Test 3

3. Consider the space-fractional convection–diffusion type of equation with variable coefficients:

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

with initial condition

u(x,0)=^x2(1−x)

homogeneous Dirichlet boundary condition

u(0,t)=0=u(1,t)

with variable convection–diffusion coefficients respectively,

c(x)=^x0.6,d(x)=Γ(2.8)^x3/4

and the forcing function

p(x,t)=2^x2(1−x)^t1.3/Γ(2.3)+0.3^{x1.8 e−t}

The exact solution is

u(x,t)=^x2(1−x)^e−t

Problem test 4

4. Consider the space-fractional convection–diffusion equation with variable coefficients:

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

with initial condition

u(x,0)=^xα(1−x)

homogeneous Dirichlet boundary condition

u(0,t)=0=u(1,t)

with variable convection–diffusion coefficients respectively,

c(x)=^x3/5,d(x)=^x3/4

and the forcing function

p(x,t)=2^xα(1−x)^t1.3/Γ(2.3)+0.3^{x1.8 e−t}

The exact solution is

u(x,t)=^xα(1−x)^e−t

Problem test 4 is experimented with the grid size reduction extrapolation approach stated in [23]. We have smooth enough numerical and exact solutions by using C-N scheme in Figure 6, and Table 5 shows the maximum error with the error rate is given for space fractional convection–diffusion problem with a grid size reduction extrapolation method.

6. Conclusions The one dimension space fractional diffusion and fractional convection–diffusion problem with space variable coefficients is solved by the fractional C-N scheme based on the Extrapolation to the limit approach of right shifted Grünwald–Letnikov approximation. The fractional C-N method, for the fractional diffusion problem and fractional convection–diffusion equation with space variable coefficients, is consistent and unconditionally stable with second order convergence. Numerical examples confirmed that the C-N method is suitable for the space fractional convection–diffusion problem even for a large value of time domain.

Figure 1. The Maximum error by C-N scheme at(T=10,Max-Error=6.5276e-07),(T=20,Max-Error=1.7244e-08),α=1.5left to right, respectively, for Example 1.

Figure 2. The exact (left) and numerical (right) solution by C-N scheme atT=1,α=1.5,τ=0.01=hfor Example 1.

Figure 3. Numerical and exact solution by C-N scheme atα=1.5,τ=h=0.01,with(T=10,T=30,T=40)left to right-down respectively, for Example 1.

Figure 4. The exact (left) and numerical (right) solution by C-N scheme for the FCDE at(h=τ=0.005,α=1.5,(t=5,max-error=4.0657e-05)for Example 2.

Figure 5. The exact (left) and numerical (right) solution by C-N scheme for the FCDE at(h=τ=0.01,(t=2,max-Error=4.2158e-04),α=1.75)for Example 2.

Figure 6. The exact (left) and numerical (right) solution by C-N scheme at(h=τ=0.0025,(t=0.1,max-Error=1.4e-03,α=1.1)for Example 4.

		α=1.25		α=1.5		α=1.8
Δt	Δx	Max-Error	Order	Max-Error	Order	Max-Error	Order
1/50	1/50	4.9807e−04	–	4.0046e−04	–	1.4048e−04	–
1/100	1/100	1.0660e−04	2.2241	8.8946e−05	2.1707	3.6848e−05	1.9307
1/200	1/200	2.4413e−05	2.1265	2.0643e−05	2.1073	9.4393e−06	1.9648
1/400	1/400	5.8239e−06	2.0676	4.9592e−06	2.0575	2.3887e−06	1.9825
1/800	1/800	1.4211e−06	2.0350	1.2146e−06	2.0296	6.0078e−07	1.9913

		T=1		T=5
Δt	Δx	Max-Error	Order	Max-Error	Order
1/50	1/50	1.4048e−04	–	2.5297e−05	–
1/100	1/100	3.6848e−05	1.9307	7.4748e−06	1.7589
1/200	1/200	9.4393e−06	1.9648	2.0122e−06	1.8933
1/400	1/400	2.3887e−06	1.9825	4.9017e−07	2.0374
1/800	1/800	6.0078e−07	1.9913	1.0620e−07	2.2065

Δt	Δx	Max-Error	Order
1/50	1/50	2.6e−03	–
1/100	1/100	7.695e−04	1.7563
1/150	1/150	2.144e−04	1.8436
1/200	1/200	5.688e−05	1.9143

		α=1.35		α=1.5		α=1.75
Δx	Max-Error	Order	Max-Error	Order	Max-Error	Order
1/50	4.5e−03	–	2.8e−03	–	1.7e−03	–
1/100	2.7e−03	0.7370	1.6e−03	0.8074	8.9641–04	0.97224
1/200	1.6e−03	0.7549	8.6405e−04	0.8889	4.6491e−04	0.8981
1/400	9.5896e−04	0.7385	4.7955e−04	0.8494	2.4086e−04	0.9488
1/800	5.7034e−04	0.7496	2.6609e−04	0.8498	1.2473e−04	0.9494

		α=1.25		α=1.55
Δt	Δx	Max-Error	Error-Rate	Max-Error	Error-Rate
1/50	1/50	1.91e−02	–	1.52e−02	–
1/100	1/100	9.9e−03	1.93	7.9e−03	1.9
1/200	1/200	5.2e−03	1.90	4.3e−03	1.84
1/400	1/400	2.8e−03	1.86	2.4e−03	1.79
1/800	1/800	1.6e−03	1.75	1.4e−03	1.7

Author Contributions

The authors contributed equally to the writing and approved the final manuscript of this paper. All authors have read and agreed to the published version of the manuscript.

Funding

This research was financially supported by the National Key Research (Grant No. 2017YFB0305601) and Development Program of China (Grant No. 2017YFB0701700).

Acknowledgments

The authors would like to thank the editor and the anonymous reviewers for their helpful comments for revising the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Diethelm, K. The Analysis of Fractional Differential Equations; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010.

2. Fomin, S.; Chugunov, V.; Hashida, T. Application of fractional differential equations for modeling the anomalous diffusion of contaminant from fracture into porous rock matrix with bordering alteration zone. Trans. Porous Med. 2010, 81, 187-205.

3. Guo, B.; Pu, X.; Huang, F. Fractional Partial Differential Equations and Their Numerical Solutions; World Scientific: Singapore, 2015.

4. Podlubny, I. Fractional Differential Equations Academic; Academic Press: New York, NY, USA, 1999.

5. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: New York, NY, USA, 1998; Volume 198.

6. Salman, W.; Gavriilidis, A.; Angeli, P. A model for predicting axial mixing during gas-liquid Taylor flow in microchannels at low Bodenstein number. Chem. Eng. J. 2004, 101, 391-396.

7. Berkowitz, B.; Cortis, A.; Dror, I.; Scher, H. Laboratory experiments on dispersive transport across interfaces. Water Resour. Res. 2009, 45.

8. Cortis, A.; Berkowitz, B. Computing "anomalous" contaminant transport in porous media. Ground Water 2005, 43, 947-950.

9. Meerschaert, M.M.; Tadjeran, C. Finite difference approximations for fractional advection-dispersion flow equations. J. Comput. Appl. Math. 2004, 172, 65-77.

10. Meerschaert, M.M.; Tadjeran, C. Finite difference approximations for two-sided space fractional partial differential equations. Appl. Numer. Math. 2006, 56, 80-90.

11. Ren, J.; Sun, Z.Z.; Zhao, X. Compact difference scheme for the fractional sub-diffusion equation with Neumann boundary conditions. J. Comput. Phys. 2013, 232, 456-467.

12. Wang, K.; Wang, H. A fast characteristic finite difference method for fractional advection-diffusion equations. Adv. Water. Resour. 2011, 34, 810-816.

13. Singh, J.; Swroop, R.; Kumar, D. A computational approach for fractional convection-diffusion equation via integral transforms. Ain Shams Eng. J. 2016, 9, 1019-1028.

14. Liu, F.; Zhuang, P.; Burrage, K. Numerical methods and analysis for a class of fractional advection-dispersion models. Comput. Math. Appl. 2012, 64, 2990-3007.

15. Tian, W.; Deng, W.; Wu, Y. Polynomial spectral collocation method for space fractional advection-diffusion equation. Numer. Methods Partial Differ. Equ. 2014, 30, 514-535.

16. Hejazi, H.; Moroney, T.; Liu, F. Comparison of finite difference and finite volume methods for solving the space fractional advection-dispersion equation with variable coefficients. ANZIAM J. 2013, 54, 557-573.

17. Liu, F.; Zhuang, P.; Turner, I.; Burrage, K.; Anh, V. A new fractional finite volume method for solving the fractional diffusion equation. Appl. Math. Model. 2014, 38, 15-16.

18. Li, C.; Zeng, F. Numerical Methods for Fractional Calculus; Chapman and Hall/CRC: Boca Raton, FL, USA, 2015.

19. Li, C.; Zeng, F. Finite difference methods for fractional differential equations. Int. J. Bifurcat. Chaos 2012, 22, 1230014.

20. Luchko, Y. Multi-dimensional fractional wave equation and some properties of its fundamental solution. Ind. Math. 2014, 6.

21. Lynch, V.E.; Carreras, B.A.; del-Castillo-Negrete, D.; Ferreira-Mejias, K.M.; Hicks, H.R. Numerical methods for the solution of partial differential equations of fractional order. J. Comput. Phys. 2003, 192, 406-421.

22. Sweilam, N.H.; Khader, M.M.; Mahdy, A.M.S. Crank-Nicolson finite difference method for solving time-fractional diffusion equation. J. Fract. Calc. Appl. 2012, 2, 1-9.

23. Tadjeran, C.; Meerschaert, M.M.; Scheffler, H.P. A second-order accurate numerical approximation for the fractional diffusion equation. J. Comput. Phys. 2006, 201, 205-213.

24. Wang, Y.M. A compact finite difference method for a class of time fractional convection-diffusion-wave equations with variable coefficients. Num. Algor. 2015, 70, 625-651.

25. Wang, Y.M.; Wang, T. A compact Alternative difference implicit method and its extrapolation for time fractional sub-diffusion equations with nonhomogeneous Neumann boundary conditions. Comput. Math. Appl. 2018, 75, 721-739.

26. Zhou, F.; Xu, X. The third kind Chebyshev wavelets collocation method for solving the time-fractional convection-diffusion equations with variable coefficients. Appl. Math. Comput. 2016, 280, 11-29.

27. Jin, B.; Lazarov, R.; Zhou, Z. A Petrov-Galerkin finite element method for fractional convection-diffusion equations. SIAM J. Numer. Anal. 2016, 54, 481-503.

28. Gao, F.; Yuan, Y.; Du, N. An upwind finite volume element method for nonlinear convection-diffusion problem. AJCM 2011, 1, 264.

29. Aboelenen, T. A direct discontinuous Galerkin method for fractional convection-diffusion and Schrödinger-type equations. Eur. Phys. J. Plus 2018, 133, 316.

30. Xu, Q.; Hesthaven, J.S. Discontinuous Galerkin method for fractional convection-diffusion equations. SIAM J. Numer. Anal. 2014, 52, 405-423.

31. Bhrawy, A.H.; Baleanu, D. A spectral Legendre-Gauss-Lobatto collocation method for a space-fractional advection-diffusion equations with variable coefficients. Rep. Math. Phys. 2013, 72, 219-233.

32. Saadatmandi, A.; Dehghan, M.; Azizi, M.R. The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients. Commun. Nonlinear. Sci. Numer. Simulat. 2012, 17, 4125-4136.

33. Pang, G.; Chen, W.; Sze, K.Y. A comparative study of finite element and finite difference methods for two-dimensional space-fractional advection-dispersion equation. Adv. Appl. Math. Mech. 2016, 8, 166-186.

34. Zhang, Y. A finite difference method for fractional partial differential equation. Appl. Math. Comput. 2009, 2015, 524-529.

35. Gu, X.M.; Huang, T.Z.; Ji, C.C. Carpentieri B, Alikhanov AA, Fast iterative method with a second-order implicit difference scheme for time-space fractional convection-diffusion equation. J. Sci. Comput. 2017, 72, 957-985.

36. Benson, D.A.; Wheatcraft, S.; Meerschaert, M.M. Application of a fractional advection-dispersion equation. Water Resour. Res. 2000, 36, 1403-1412.

37. Tian, W.; Zhou, H.; Deng, W. A class of second order difference approximations for solving space fractional diffusion equations. Math. Comp. 2015, 84, 1703-1727.

38. Henrici, P. Elements of Numerical Analysis; John Wiley and Sons: New York, NY, USA, 1984.

39. Tuan, V.K.; Gorenflo, R. Extrapolation to the Limit for Numerical Fractional Differentiation. Z. Angew. Math. Mech. 1995, 75, 646-648.

40. LeVeque, R.J. Finite Difference Methods for Ordinary and Partial Differential Equations; SIAM: Philadelphia, PA, USA, 2007; Volume 98.

41. Isaacson, E.; Keller, H.B. Analysis of Numerical Methods; Courier Corporation: Chelmsford, MA, USA, 2012.