Numerical approximation of the space-time Caputo-Fabrizio fractional derivative and application to groundwater pollution equation

Abstract

Recently, Caputo and Fabrizio proposed a new derivative with fractional order without singular kernel. The derivative has several interesting properties that are useful for modeling in many branches of sciences. For instance, the derivative is able to describe substance heterogeneities and configurations with different scales. In order to accommodate researchers dealing with numerical analysis, we propose a numerical approximation in time and space. We modified the advection dispersion equation by replacing the time derivative with the new fractional derivative. We solve numerically the modified equation using the proposed numerical approximation. The stability and convergence analysis of the numerical scheme were presented together with some simulations.

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 DOI 10.1186/s13662-016-0871-x

R E S E A R C H Open Access

http://crossmark.crossref.org/dialog/?doi=10.1186/s13662-016-0871-x&domain=pdf

Web End = Numerical approximation of the space-time Caputo-Fabrizio fractional derivative and application to groundwater pollution equation

Abdon Atangana1 and Rubayyi T Alqahtani2*

*Correspondence: mailto:[email protected]

Web End [email protected]

2Department of Mathematics and Statistics, College of Science, Al-Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, 11566, Saudi ArabiaFull list of author information is available at the end of the article

Abstract

Recently, Caputo and Fabrizio proposed a new derivative with fractional order without singular kernel. The derivative has several interesting properties that are useful for modeling in many branches of sciences. For instance, the derivative is able to describe substance heterogeneities and congurations with dierent scales. In order to accommodate researchers dealing with numerical analysis, we propose a numerical approximation in time and space. We modied the advection dispersion equation by replacing the time derivative with the new fractional derivative. We solve numerically the modied equation using the proposed numerical approximation. The stability and convergence analysis of the numerical scheme were presented together with some simulations.

Keywords: Caputo-Fabrizio derivative; numerical approximation; advection diusion equation; stability analysis

1 Introduction

In the last decade, many physical problems have been modeled using the concept of noninteger-order derivative. These derivatives of fractional order range from Riemann-Liouville via Caputo to Caputo-Fabrizio [, ]. We can nd in the literature many analytical approaches to deal with dierential equations with fractional equations []. Most of these techniques are dealing with linear fractional dierential equations. However, most fractional dierential equations describing real-world problems are highly complicated and cannot sometime be handled via analytical methods. In order to solve these problems in many cases, researchers rely on the use of numerical methods because these problems have initial conditions, boundary condition, and source terms that turn hard to nd an analytical solution.

Several numerical approaches in connection with derivatives of fractional order describing real-world problems alter essentially in the many in which the derivative of fractional order is tailored []. Approximation representation of a derivative of fractional order has a highly complicated formula compared to those of integer order because fractional derivatives are nonlocal, and therefore the calculation at a particular point requires knowl-

2016 Atangana and Alqahtani. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 2 of 13

edge of the function further out of the region close to that point. Accordingly, nite difference approximations of derivatives of fractional order engage a quantity of points that alters according to how faraway we are from the border line [].

The most recent derivative of fractional order was proposed by Caputo and Fabrizio [], who demonstrated that the new-fangled derivative encompasses extra encouraging properties in comparison with the old version. They demonstrated, for example, that it can depict substance heterogeneities and congurations with dierent scales, which obviously cannot be overseeing with the prominent local theories and also the well-known fractional derivative. An additional application is in the investigation of the macroscopic behaviors of some materials, associated with nonlocal communications between atoms, which are recognized to be important of the properties of material. We present the denition of the Caputo fractional derivative.

Denition The Caputo derivative of fractional order old editor of a function f is given as

C Dx[parenleftbig]f (x)[parenrightbig] =

(n )

[integraldisplay] x

(x t)n

dndtn [parenleftbig]f (t)[parenrightbig] dt, n <

n. ()

Denition Let f H(a, b), b > a, [, ]. Then the new Caputo fractional derivative is dened as

CF Dt[parenleftbig]f (t)[parenrightbig] = M(

)

[bracketrightbigg] dx, ()

where M() is a normalization function such that M() = M() = []. However, if the function does not belong to H(a, b) then, the derivative can be redened as

CF Dt[parenleftbig]f (t)[parenrightbig] =

M()

[integraldisplay] ta f (x)

exp[bracketleftbigg] t x

[bracketrightbigg] dx. ()

The aim of this paper is to propose a numerical approximation of the space and time Caputo-Fabrizio derivative of fractional order that will be used by researchers in the eld of numerical analysis.

2 Caputo-Fabrizio approximations

In this section, we derive a numerical approximation based upon the denition of newly proposed derivative of fractional order [],

CF Dt[parenleftbig]f (t)[parenrightbig] = M(

)

[integraldisplay] t a

[parenleftbig]f (t) f (x)[parenrightbig] exp[bracketleftbigg] t x

[bracketrightbigg] dx. ()

For some positive integer N, the grid size in time for nite dierence technique is dened by

k =

N .

The grid points in the time interval [, T] are labeled tn = nk, n = , , , . . . , TN. The value of the function f at the grid point is fi = f (ti).

[integraldisplay] t

f (x)

exp[bracketleftbigg] t x

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 3 of 13

A discrete approximation to the Caputo-Fabrizio derivative of fractional order can be obtained by simple quadrature formula as follows:

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

[bracketrightbigg] dx. ()

This equation can be modied using the rst-order approximation to

CF Dt[parenleftbig]f (tj)[parenrightbig] = M(

)

[integraldisplay] tn

f (x)

exp[bracketleftbigg] tn x

n

[summationdisplay]

j=

[integraldisplay] jk

(j)k

[parenleftbigg]f k+ f k

t + O(

t)[parenrightbigg] exp[bracketleftbigg] tj x

[bracketrightbigg] dx. ()

Before integration we obtain the following expression

M()

n

j=

[parenleftbigg]f j+ f j

[summationdisplay] t + O(

t)[parenrightbigg] [integraldisplay]

jk

(j)k

exp[bracketleftbigg] tn x

[bracketrightbigg] dx, ()

CF Dt[parenleftbig]f (tj)[parenrightbig] = M(

)

n

j=

[parenleftbigg]f j+ f j

[summationdisplay] t + O(

t)[parenrightbigg]dj,k,

where

dj,k = exp[bracketleftbigg] k

(n j + )

[bracketrightbigg] exp[bracketleftbigg] k

(n j)

[bracketrightbigg]. ()

We nally have that

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

n

[summationdisplay]

j=

[parenleftbigg]f j+ f j t

[parenrightbigg]dj,k + M(

)

n

[summationdisplay]

j=

dj,kO( t). ()

Theorem Let f (x) be a function in C[a, b], and let the order of the fractional derivative be < . Then the rst-order approximation of the Caputo-Fabrizio derivative at a point tn is

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

n

[summationdisplay]

j=

[parenleftbigg]f j+ f j

t

[parenrightbigg]dj,k + O[parenleftbig]( t)[parenrightbig]. ()

Proof From equation () we have

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

n

[summationdisplay]

j=

[parenleftbigg]f j+ f j

t

[parenrightbigg]dj,k

+ M(

)

n

j=

[summationdisplay] [parenleftbigg]exp[bracketleftbigg] k

(n j + )

[bracketrightbigg] exp[bracketleftbigg] k

(n j)

[bracketrightbigg][parenrightbigg]O( t).

However,

n

j=

[summationdisplay] [parenleftbigg]exp[bracketleftbigg] k

(n j + )

[bracketrightbigg] exp[bracketleftbigg] k

(n j)

[bracketrightbigg][parenrightbigg] = exp[bracketleftbigg] k

(n)

[bracketrightbigg] . ()

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 4 of 13

Now the approximation of the exponential function can be obtained as

exp[bracketleftbigg] k

(n)

[bracketrightbigg]

k (n). ()

Then replacing the above in equation (), we obtain

n

j=

[summationdisplay] [parenleftbigg]exp[bracketleftbigg] k

(n j + )

[bracketrightbigg] exp[bracketleftbigg] k

(n j)

[bracketrightbigg][parenrightbigg]

k (n). ()

Then equation () becomes

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

n

[summationdisplay]

j=

[parenleftbigg]f j+ f j t

[parenrightbigg]dj,k + M(

)k (n)O(

t). ()

We therefore obtain the requested result

[parenrightbigg]dj,k + O[parenleftbig]( t)[parenrightbig]. ()

This completes the proof.

We now conclude that the rst-order approximation method for the computation of the Caputo-Fabrizio derivative of time fractional order is given as

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

CF Dt[parenleftbig]f (tn)[parenrightbig] = M(

)

n

[summationdisplay]

j=

[parenleftbigg]f j+ f j t

[parenleftbigg]f j+ f j

t

[parenrightbigg]dj,k. ()

We next propose the rst order for the space fractional order.

For some positive integer N, the grid sizes in time for nite dierence technique is dened by

i =

M .

The grid points in the time interval [, X] are labeled xi = mi, m = , , , . . . , XM.

The value of the function f at the grid point is f ki = f (xi, tk). We have

CF Dt[parenleftbig]f (xm, ti)[parenrightbig] = M(

) ( )

n

[summationdisplay]

j=

[integraldisplay] xm

yf (y, ti)

exp[bracketleftbigg] (xm y)

( )

[bracketrightbigg] dy. ()

Now employing the Crank-Nicolson approximation for the rst-order derivative, the above equation is converted to

CF Dx[parenleftbig]f (xm, tk)[parenrightbig] = M(

) ( )

[integraldisplay] xm

[parenleftbigg](f k+i+ f k+i) (f ki+ f ki)

x + O(

t)[parenrightbigg]

exp[bracketleftbigg] (xm y)

( )

[bracketrightbigg] dy. ()

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 5 of 13

The latter equation can be converted to

CF Dx[parenleftbig]f (xm, ti)[parenrightbig] = M(

) ( )

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl) x + O(i)

[bracerightbigg]

[integraldisplay] li

(l)i

exp[bracketleftbigg] (im y)

( )

[bracketrightbigg] dy, ()

where the integral part is given as

[integraldisplay] li

(l)i

exp[bracketleftbigg] (im y)

( )

[bracketrightbigg] dy

= (

)

[braceleftbigg]erf[bracketleftbigg](mi li)

[bracketrightbigg] erf[bracketleftbigg](mi li + i)

[bracketrightbigg][bracerightbigg], ()

so that equation () becomes

CF Dx[parenleftbig]f (xm, tk)[parenrightbig]

= M(

)

[bracerightbigg]

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl) x + O(i)

( )

[braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg]. ()

From the above we obtain

CF Dx[parenleftbig]f (xm, tk)[parenrightbig]

= M(

)

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl) x

( )

[bracketrightbigg][bracerightbigg][bracerightbigg]

[braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

+ O(i)(

)

[bracketrightbigg][bracerightbigg]. ()

Theorem Let f (x, t) be a function in C([a, b] [, T]), and let the order of the fractional derivative be < . Then the rst-order approximation of the Caputo-Fabrizio derivative at a point (xm, tn) is

CF Dx[parenleftbig]f (xm, tk)[parenrightbig] = M(

)

m

l=

[summationdisplay] [braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl)

x

[bracerightbigg]di,l + R(, i, l), ()

where

di,l = [braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg], [vextenddouble][vextenddouble]R(, i, l)[vextenddouble][vextenddouble] < M.

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 6 of 13

Proof From equation () we have that

CF Dx[parenleftbig]f (xm, tk)[parenrightbig]

= M(

)

[braceleftbigg]erf[bracketleftbigg](m l)

[bracketrightbigg]

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl) x

i

erf[bracketleftbigg](m l + )

[bracketrightbigg][bracerightbigg][bracerightbigg]

i

+ O(i)M(

)

m

l=

[summationdisplay] [braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg].

We put

R(, i, l) = O(i)M(

)

m

l=

[summationdisplay] [braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg]. ()

Then taking the norm to both sides, we have

[vextenddouble][vextenddouble]R(, i, l)[vextenddouble][vextenddouble] =

[vextenddouble][vextenddouble][vextenddouble][vextenddouble][vextenddouble]

O(i)M(

)

m

l=

[summationdisplay] [braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg][vextenddouble][vextenddouble][vextenddouble][vextenddouble][vextenddouble]

, ()

[vextenddouble][vextenddouble]R(, i, l)[vextenddouble][vextenddouble] =

[vextenddouble][vextenddouble][vextenddouble][vextenddouble]O(i)

[bracketrightbigg][parenrightbigg][vextenddouble][vextenddouble][vextenddouble][vextenddouble].

This completes the proof.

Then, the rst-order approximation method for the computation of Caputo-Fabrizio derivative of space fractional order is given as

CF Dx[parenleftbig]f (xm, ti)[parenrightbig] = M(

)

M()

[parenleftbigg]erf[bracketleftbigg]m

i

m

[summationdisplay]

l=

[braceleftbigg](f k+l+ f k+l) (f kl+ f kl)

x

[braceleftbigg]erf[bracketleftbigg](m l)

i

[bracketrightbigg] erf[bracketleftbigg](m l + )

i

[bracketrightbigg][bracerightbigg][bracerightbigg]. ()

3 Application to some well-known equations

In this section, we present a numerical solution of the time fractional advection diusion equation in heterogeneous medium. The fractional derivative used here is of the Caputo-Fabrizio type.

The reason of this modication is that the fractional derivatives are recollection operational which recurrently distinguish indulgence of force or damage in the passable as in the case of inelastic media or reconsideration of the porosity in the thinning out in permeable media and supplementary in comprehensive they are in traditional values throughout the subsequent theory of hydrology. They are accredited not merely for the motivation that they match appropriately a variety of noticeable actuality, nevertheless, additionally for the motive that they own the well-designed alongside with scrupulous property that although the order of dierentiation is integer, they match by means of the traditional derivative of that order. On the other hand, this chattel is not pertinent to the eect they characterize in

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 7 of 13

the physical observable fact and conjectures if using other dierential operators, probably simpler nevertheless devoid of this property, one may get similar responds of fractional order derivative. Therefore, in order to well replicate the ow of the particles via porous media in dierent scale in the medium, we replace the ordinary derivative in time with the scale time derivative proposed by Caputo and Fabrizio. The equation under consideration here is

CF Dt[parenleftbig]P(x, t)[parenrightbig] +

u c

P(x, t)

P(x, t)

c . ()

In equation (), four terms represents transient, advection, and source terms, respectively, P(x, t) is the particle, heat, pollution, or other physical quantities; c is the specic of heat, particles, or other physical quantities; is the porosity that is the ratio of the liquid volume to the total volume of the medium via which the ow is taken place, , are the mass density and thermal conductivity, respectively, and, nally, Q(x, t) is the source term. Now substituting equation () into (), we obtain

M()

j+

[summationdisplay]

k=

x

x =

Q(x, t)

Pki Pki

[parenleftbigg]erf[parenleftbigg](j k)

k

[parenrightbigg] erf[parenleftbigg](j k + )

[parenrightbigg][parenrightbigg]

k

= hc [braceleftbig][parenleftbig]Pj+i+ Pj+i + Pj+i[parenrightbig] + [parenleftbig]Pji+ Pji + Pji[parenrightbig][bracerightbig]

u c

[braceleftbig][parenleftbig]Pj+

i+ Pj+i[parenrightbig] + [parenleftbig]Pji+ Pji[parenrightbig][bracerightbig] + Qj+i + Qji

c . ()

The above equation can be converted to

M()

[parenleftBigg]Pj+i Pji

+

j

[summationdisplay]

k=

Pj+ki Pjki

[parenrightBigg]

dk,j

= hc [braceleftbig][parenleftbig]Pj+i+ Pj+i + Pj+i[parenrightbig] + [parenleftbig]Pji+ Pji + Pji[parenrightbig][bracerightbig]

uhc [braceleftbig][parenleftbig]Pj+i+ Pj+i[parenrightbig] + [parenleftbig]Pji+ Pji[parenrightbig][bracerightbig] +

Qj+i + Qji

c . ()

For simplicity, let us put

a = M(

) , b =

u hc .

hc , c =

Rearranging, we obtain the following recursive formula:

[parenleftbig]ad

k,j + b[parenrightbig]Pj+

i = [parenleftbig]adk,j b[parenrightbig]Pji + a

j

k=

[summationdisplay] [parenleftbig]Pj+k

i Pjki[parenrightbig]dk,j

+ b[braceleftbig][parenleftbig]Pj+

i+ + Pj+i[parenrightbig] + [parenleftbig]Pji+ + Pji[parenrightbig][bracerightbig] c[braceleftbig][parenleftbig]Pj+i+ Pj+i[parenrightbig] + [parenleftbig]Pji+ Pji[parenrightbig][bracerightbig]

+ Qj+i + Qji

c . ()

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 8 of 13

4 Stability analysis of the numerical scheme

We present in this section the stability analysis of the Crank-Nicolson scheme for time fractional advection diusion equation. For this, we let eji = Pji pji with pji the approximate solution at the point (xi, tj) (i = , , . . . , N, j = , , , . . . , M) and, as usual, ej = [ej, . . . , ejN]T.

The error committed while solving the time fractional advection diusion equation with the Crank-Nicolson scheme is

[parenleftbig]ad

k,j + b[parenrightbig]ej+

i = [parenleftbig]adk,j b[parenrightbig]eji + a

j

k=

[summationdisplay] [parenleftbig]ej+k

i ejki[parenrightbig]dk,j

+ b[braceleftbig][parenleftbig]ej+

i+ + ej+i[parenrightbig] + [parenleftbig]eji+ + eji[parenrightbig][bracerightbig] c[braceleftbig][parenleftbig]ej+i+ ej+i[parenrightbig] + [parenleftbig]eji+ eji[parenrightbig][bracerightbig]

+ Qj+i + Qji

c . ()

Here, we assume that

eji = f (j) exp(ij), ()

where is the real spatial wave number []. However, substituting equation () into equation (), we obtain, for j = ,

[parenleftbigg]adk, + b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f () = [parenleftbigg]adk, b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f (), ()

and for j > , we have

[parenleftbigg]adk,j + b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f (j) = [parenleftbigg]adk,j b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f (j )

a

j

[summationdisplay]

l=

f (j l)dk,j + f (j + )dk,. ()

Theorem Assume that f (k) satises equations () and (). Then, for all k > ,

[vextendsingle][vextendsingle]f (j)[vextendsingle][vextendsingle] [vextendsingle][vextendsingle]f ()[vextendsingle][vextendsingle]. ()

Proof We shall prove this theorem by employing the recursive method on the natural number j. Then, when j = , we have equation (), and we reformulate it as follows:

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

f () f ()

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

= [vextendsingle][vextendsingle][vextendsingle][vextendsingle]

(adk, b sin(i)) (adk, + b sin(i))

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

. ()

This implies

[vextendsingle]f ()[vextendsingle][vextendsingle] [vextendsingle][vextendsingle]f ()[vextendsingle][vextendsingle].

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 9 of 13

The property is veried for j = . Let us assume that this property is also satised for any j . We shall verify that the property holds also for j + :

[parenleftbigg]adk,j + b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f (j + ) = [parenleftbigg]adk,j b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg]f (j) a

j

[summationdisplay]

l=

f (j l)dk,l. ()

Now taking the norms of both sides of equation (), we obtain

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]ad

k,j + b sin[parenleftbigg]

i

[parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[vextendsingle][vextendsingle]f (j + )[vextendsingle][vextendsingle]

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[parenleftbigg]adk,j b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[vextendsingle][vextendsingle]f (j)[vextendsingle][vextendsingle] +

j

l=

[summationdisplay] [vextendsingle][vextendsingle]f (j l)[vextendsingle][vextendsingle]dk,l. ()

Nonetheless, we recall that the property holds up to j. Thus, we transform the above equation into

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]ad

k,j + b sin[parenleftbigg]

i

[parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[vextendsingle][vextendsingle]f (j + )[vextendsingle][vextendsingle]

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[parenleftbigg]adk,j b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[vextendsingle][vextendsingle]f ()[vextendsingle][vextendsingle] +

j

l=

[summationdisplay] [vextendsingle][vextendsingle]f ()[vextendsingle][vextendsingle]dk,l.

Rearranging, we obtain

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]ad

k,j + b sin[parenleftbigg]

i

[parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[vextendsingle][vextendsingle]f (j + )[vextendsingle][vextendsingle]

[braceleftbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

[parenleftbigg]adk,j b sin[parenleftbigg]

i

[parenrightbigg][parenrightbigg][vextendsingle][vextendsingle][vextendsingle][vextendsingle]

+ [braceleftbigg]erf[bracketleftbigg] m

i

[bracketrightbigg] erf[bracketleftbigg]

i

[bracketrightbigg][bracerightbigg][bracerightbigg][vextendsingle][vextendsingle]f ()[vextendsingle][vextendsingle]. ()

It is important to recall that

[vextendsingle][vextendsingle]erf[x][vextendsingle][vextendsingle] , erf

[bracketleftbigg] mi

[bracketrightbigg] erf[bracketleftbigg]

[bracketrightbigg]

. ()

Therefore,

|f (j + )|

i

|f ()|

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

|(adk,j b sin(i))| + {erf[mi] erf[ i ]} adk,j + b sin(i)

[vextendsingle][vextendsingle][vextendsingle][vextendsingle]

.

Then,

|f (j + )|

. ()

The property also holds for j + . According to the inductive technique, the property is satised for any natural number. This completes the proof of Theorem . Theorem shows that the Crank-Nicolson scheme is stable for the advection diusion equation with the time fractional Caputo-Fabrizio derivative.

5 Convergence analysis of the numerical solution

Let us suppose that, at the point (xi, tj), the exact solution of our considered equation is P(xi, tj) (i = , , , . . . , N; j = , , , , . . . , M). We assume that the dierence between the exact solution and the approximate solution at that particular point is provided by

|f ()|

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 10 of 13

ji = P(xi, tj) Pji. The transpose matrix associated with the matrix ji = P(xi, tj) Pji (i = , , , , . . . , N; j = , , , , , . . . , M) is (j, j, . . . , jN)T. However, the row is zero because it represents the initial condition. The recursive relation in connection with the Crank-

Nicolson scheme for the time fractional advection diusion equation is given as

[parenleftbig]ad

k, + b[parenrightbig]i + (c b)[parenleftbig]i+ + i[parenrightbig] Qi Qi

c = Ri for j = , ()

[parenleftbig]ad

i,j + b[parenrightbig]j+

i [parenleftbig]adi,j + b[parenrightbig]ji + (c b)[parenleftbig]j+i+ + j+i[parenrightbig] + (c b)[parenleftbig]ji+ + ji[parenrightbig] Qj+i Qji

c

= a

j

l=

[summationdisplay] b(xi, tjl)di,l b[braceleftbig][parenleftbig]P(xi+, tj+) P(xi, tj+) + P(xi, xj+)[parenrightbig]

+ [parenleftbig]P(xi+, tj) P(xi, tj) + P(xi, xj)[parenrightbig][bracerightbig] + c[braceleftbig][parenleftbig]P(xi+, tj+) + P(xi, xj+)[parenrightbig]

+ [parenleftbig]P(xi+, tj) + P(xi, xj)[parenrightbig][bracerightbig] Qj+i Qji

c . ()

Now using the full approximation and considering the remaining terms, we obtain the following relation:

Rj+i D[parenleftbig]k + h[parenrightbig]. ()

Theorem The Crank-Nicolson scheme for the advection diusion equation with time fractional Caputo-Fabrizio derivative converges, and there exists a positive constant D such that

[vextenddouble][vextenddouble]P(xi, ti) Pji[vextenddouble][vextenddouble] D[parenleftbig]k + h[parenrightbig] for all (i = , , , . . . , M, j = , , , . . . , N). ()

6 Numerical simulation for different values of alpha

In this section, using the new numerical scheme, we present the numerical simulation of the advection diusion equation with the Caputo-Fabrizio derivative of fractional order for dierent values of alpha. We chose in this case Q(x, t) = sin[x + ], P(x, ) = cos[x],

P(, t) = cos[t], P(x, ) = , u = ., c = ., = , = .. The numerical simulations are depicted in Figures , , , and . It is worth noting that each gure represents the ow at scale alpha. It is very important to realize that fractional dierentiation is able to control the variabilities of the plume movement within the geological formations. The pollution does not only move within a homogeneous medium but also within heterogeneous one; therefore, the plume paths cannot be predicted by the classical advection dispersion equation. In these gures, we can see that the proportionally of the density of pollution within the geological formation is not the same everywhere due to the heterogeneity, and

j

[summationdisplay]

l=

jlidi,l + Rj+i for j > .

The remainder term of the approximation for using the Crank-Nicolson scheme to solve the modied advection diusion equation is given in this case as

Rj+i = aP(xi, tj+)di,l + a

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 11 of 13

Figure 1 Numerical simulation at scale alpha = 0.15.

Figure 2 Numerical simulation at scale alpha = 0.55.

this is better described via the concept of fractional dierentiation with nonsingular kernel.

7 Conclusions

We have proposed in this work the numerical approximation of the newly proposed derivative of fractional order in order to t this derivative in the scope of numerical investigations. The new derivative is easy to use even numerically and display important characteristics that cannot be observed in the commonly used fractional derivatives. In order to test the possible application of the new numerical approximation of the new Caputo-Fabrizio derivative of fractional order, we presented a model of advection diusion equation with the time fractional of the new derivative. We solved this equation numerically using the Crank-Nicolson technique. We showed the stability analysis together with some numerical simulations for dierent values of alpha.

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 12 of 13

Figure 3 Numerical simulation at scale alpha = 0.85.

Figure 4 Numerical simulation at scale alpha = 0.95.

Competing interests

The authors declare that they have no competing interests.

Authors contributions

The authors have equal contributions to each part of this paper. All the authors read and approved the nal manuscript.

Author details

1Faculty of Natural and Agricultural Sciences, Institute for Groundwater Studies, Bloemfontein, 9300, South Africa.

2Department of Mathematics and Statistics, College of Science, Al-Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, 11566, Saudi Arabia.

Acknowledgements

We acknowledge the Editorial Board and the referees for their eorts and constructive criticism, which have improved the manuscript.

Received: 17 April 2016 Accepted: 24 May 2016

References

1. Caputo, M, Fabrizio, M: A new denition of fractional derivative without singular kernel. Prog. Fract. Dier. Appl. 1(2), 73-85 (2015)

Atangana and Alqahtani Advances in Dierence Equations (2016) 2016:156 Page 13 of 13

2. Losada, J, Nieto, JJ: Properties of a new fractional derivative without singular kernel. Prog. Fract. Dier. Appl. 1(2), 87-92 (2015)

3. Tian, Y, Wang, JR: On some Hermite-Hadamard type inequalities for convex functions via Hadamard fractional integrals. Prog. Fract. Dier. Appl. 1(2), 103-110 (2015)

4. Atangana, A, Nieto, JJ: Numerical solution for the model of RLC circuit via the fractional derivative without singular kernel. Adv. Mech. Eng. 7(10), 1-7 (2015)

5. Diethelm, K, Ford, NJ, Freed, AD: Detailed error analysis for a fractional Adams method. Numer. Algorithms 36(1), 31-52 (2004)

6. Li, CP, Tao, CX: On the fractional Adams method. Comput. Math. Appl. 58(8), 1573-1588 (2009)7. Podlubny, I, Chechkin, A, Skovranek, T, Chen, Y, Vinagre Jara, BM: Matrix approach to discrete fractional calculus II: partial fractional dierential equations. J. Comput. Phys. 228(8), 3137-3153 (2009)

8. Atangana, A, Secer, A: A note on fractional order derivatives and table of fractional derivatives of some special functions. Abstr. Appl. Anal. 2013, Article ID 279681 (2013)

9. Kilbas, AA, Srivastava, HM, Trujillo, JJ: Theory and Applications of Fractional Dierential Equations. North-Holland Mathematics Studies, vol. 204. Elsevier, Amsterdam (2006)

10. Ashyralyev, A: A note on fractional derivatives and fractional powers of operators. J. Math. Anal. Appl. 357(1), 232-236 (2009)

11. Hanert, E: On the numerical solution of space-time fractional diusion models. Comput. Fluids 46(1), 33-39 (2011)12. Zhuang, P, Liu, F, Anh, V, Turner, I: Numerical methods for the variable-order fractional advection-diusion equation with a nonlinear source term. SIAM J. Numer. Anal. 47(3), 1760-1781 (2009)

13. Lin, R, Liu, F, Anh, V, Turner, I: Stability and convergence of a new explicit nite-dierence approximation for the variable-order nonlinear fractional diusion equation. Appl. Math. Comput. 212(2), 435-445 (2009)

14. Crank, J, Nicolson, P: A practical method for numerical evaluation of solutions of partial dierential equations of the heat conduction type. Proc. Camb. Philos. Soc. 43(1), 50-67 (1947)

15. Atangana, A, Oukouomi Noutchie, SC: Stability and convergence of a time-fractional variable order Hantush equation for a deformable aquifer. Abstr. Appl. Anal. 2013, Article ID 691060 (2013)

16. Akyz-Dasclolu, A: Chebyshev polynomial solutions of systems of linear integral equations. Appl. Math. Comput. 151(1), 221-232 (2004)

17. El Tom, MEA: Application of spline functions to systems of Volterra integral equations of the rst and second kinds. IMA J. Appl. Math. 17(3), 295-310 (1976)

18. Saeed, RK, Ahmed, C: Approximate solution for the system of non-linear Volterra integral equations of the second kind by using block-by-block method. Aust. J. Basic Appl. Sci. 2(1), 114-124 (2008)

19. Maleknejad, K, Karami, M: Numerical solution of non-linear Fredholm integral equations by using multiwavelets in the Petrov-Galerkin method. Appl. Math. Comput. 168(1), 102-110 (2005)

20. Meerschaert, MM, Tadjeran, C: Finite dierence approximations for fractional advection-dispersion ow equations.J. Comput. Appl. Math. 172(1), 65-77 (2004)21. Doungmo Goufo, EF: Application of the Caputo-Fabrizio fractional derivative without singular kernel to Korteweg-de Vries-Bergers equation. Math. Model. Anal. 21(2), 188-198 (2016)

22. Doungmo Goufo, EF, Pene, MK, Mwambakana, JN: Duplication in a model of rock fracture with fractional derivative without singular kernel. Open Math. 13(1), 839-846 (2015)

23. Gmez-Aguilar, JF, Lpez-Lpez, MG, Alvarado-Martnez, VM, Reyes-Reyes, J, Adam-Medina, M: Modeling diusive transport with a fractional derivative without singular kernel. Phys. A, Stat. Mech. Appl. 447, 467-481 (2016)24. Gmez-Aguilar, JF, Crdova-Fraga, T, Escalante-Martnez, JE, Caldern-Ramn, C, Escobar-Jimnez, RF: Electrical circuits described by a fractional derivative with regular kernel. Rev. Mex. Fis. 62(2), 144-154 (2016)