1 Introduction

Replacement of the second order derivative in the conventional advection-diffusion equation(ADE)by a derivative of orderα∈(1,2)yields the fractional diffusion term and thus a fractional advection-diffusion equation(FADE).The superdiffusive flow feature of FADE is suitable for describing contaminant transport in groundwater[Benson,Wheatcraft and Meerschaert(2000);Benson,Schumer,Meerschaert and Wheatcraft(2001);Zhang,Benson and Reeves(2009)].For 2D and 3D FADE,Meerschaert et al[Meerschaert,Benson and Baumer(1999)]introduced a directionally-weighted fractional diffusion operatorwhich is an extension of the well-known fractional Laplacian withMdesignating a directional probability measure[Chen and Holm(2004);Pang,Chen,and Sze(2013);Duan(2005);Yang,Liu and Turner(2010);Katzav and Adda-bedia(2008)].The inaccessibility of closed-form analytical solutions of FADE fosters the development of the relevant numerical solution techniques.

Although numerical methods for time-fractional derivative equations are extensively investigated[Zeng,Li,Liu and Turner(2013);Li,Chen and Ye(2011);Wang,Meng,Ma and Wu(2013);Chen,Sun,Li and Fu(2013)],numerical solution of space-fractional derivative equations,especially in high-dimensional case,is less reported.Finite difference method(FDM)[Yang,Liu,and Turner(2010)], finite element method(FEM)[Huang,Huang and Zhan(2008);Zheng,Li and Zhao(2010)], finite volume method[Zhang,Craw ford,Deeks,Stutter,Bengough and Young(2005)],matrix transform technique(MTT)[Ilic,Liu,rner and Anh(2005);Ilic,Liu,Turner and Anh(2006)]and spectral methods[Hanert(2010);Li and Xu(2010)]have been employed to solve 1D FADE.In solving FADEs in high spatial dimensions,FDM s[Meerschaert,Scheffler and Tadjeran(2006);Tadjeran and Meerschaert(2007);Wang and Wang(2011)]are commonly used.Unfortunately,they are limited to the fractional differential operators in coordinate directions whilst the generalized operatorinvolves the fractional derivative along all directions.

To the best of our know ledge,FEM[Roop(2006)],vector Grünwald formula(VGF)[Meerschaert and Mortenson(2004)]and MTT[Yang,Turner,Liu,and Ilic(2011)]are the available methods for approximatingBased on the variational statements in[Ervin and Roop(2006)],Roop presented a finite element scheme for fractional advection-diffusion equation withand the convergence rate was proven to be equal to or exceed two,i.e.,O(hd)whered≥2 andhis the grid or element size.Despite the relatively high accuracy,the spatial nonlocality of spacefractional derivative leads to element stiffness matrices which couple all nodal dof.Accordingly,the global stiffness matrix is non-sparse even though the trial functions are local.These significantly increase the assembling and solution time of FEM.The VGF is actually the multi-variate extension of the standard one-variate Grünwald-Letnikov formula[Yang,Liu and Turner(2010)]for fractional derivative and the method has the first-order accuracy,i.e.,O(h).MTT is widely employed for approximating fractional Laplacian(-?)α/2but may not be a good choice for approximatingas they have different forms of definitions:the former is defined by the spectral decomposition while the latter is defined by a vector integral.

Pitfalls of the three afore-mentioned methods motivateus to seek more accurate and efficient alternatives Differential quadrature method(DQM)[Shu(2000);Chen,Zhong and Shu(1998)]is well-known for its high accuracy and low computational effort.In the current literature,DQM has only been applied to 1D diffusion equations with time fractional derivative and the DQ formula was con fined to approximating the spatial integer-order derivatives[Mokhtazri(2011)].In this paper we consider the 1D FADE with space-fractional derivative and apply DQ formulas to both spatial integer-order and fractional derivatives.Furthermore we shall apply the differential cubature method(DCM)to solve 2D FADE and go further in developing an improved DCM with better efficiency.

In the integration quadrature method,a weighted sum of the integrands at the quadrature points is taken to approximate an integral.Similarly,in the DQM,the derivative of a given function at each sampling point in a coordinate direction is expressed in terms of a weighted sum of the functional values of all the sampling points along the same direction[Bert and Malik(1996)].Lagrange interpolants are most commonly used as the test functions of DQM to derive the weighting coefficients yet the Fourier and radial basis functions can also be employed[Shu,Ding and Yeo(2003);Shu and Chew(1997)].This paper only considers the polynomial test functions.The distribution of sampling points and the accuracy of weighting coefficients computed are critical in the accuracy of the DQM.When the sampling points are Gauss-Chebyshev or Gauss-Legendre-Lobatto points,DQM is equivalent to Chebyshev or Legendre collocation method[Lui(2011)].Though explicit expressions of the DQ weighting coefficients for integer-order derivatives are available[Quan and Chang(1989);Shu and Richards(1992)],those for fractional derivatives are difficult,if not impossible,to derive.Here,they are obtained numerically by matrix inversion.For a large number of sampling points,the matrix becomes ill-conditioned and the accuracy as well as the convergence of the DQM may be adversely affected It will be seen,however,from the numerical tests that matrix inversion remains to be an acceptable expedient.

Unlike DQM which employs only the sampling points along a coordinate direction,DCM uses the weighted sum of functional values at the sampling points in both coordinate and non-coordinate directions to approximate the derivatives.This method is originally developed for numerical quadrature of multi-dimensional functions and subsequently used to solve high-dimensional differential equations.A comparative study of DQM and DCM for solving conventional advection-diffusion-reaction problems can be found in[Malik and Civan(1995)].

The rest of paper is organized as follows.In Section 2,the DQ,DC and improved DC formulas for fractional directional derivatives are presented.Note that a weighted sum of fractional directional derivatives constitutes the directionally-weighted fractional diffusion operator(see Section 3 below).Section 3 briefly reviews the definition of,followed by’s DQ and DC approximations.Section 4 presents some preliminary numerical results for steady-state 1D and 2D FADEs.Section 5 gives some conclusions

2 Approximation for fractional derivatives in a given direction

In this section,DQ and DC formulas are employed to approximate fractional directional derivatives of one-and two-variate functions,respectively.An improved DC formula is developed for a fast computation of weighting coefficients of DCM

2.1 DQ formula

In DQ formula for integer-order derivative,then-th order derivative is approximated by a linear weighted sum of functional values at sampling points,i.e.

wherexiis the sampling point,andis the weighting coefficient.The sampling points can be equally or unequally spaced[Bert and Malik(1996)]:

Now,consider the fractional-order counterpart of(1),i.e.,forα∈(n-1,n]andx∈[ab],

In particular,θ=0 anddefine respectively the left and right Caputo fractional differential operators[Podlubny(1999);Samko,Kilbas and Marichey(1993)],i.e.

It is trivial that[Kilbas,Srivastava and Trujillo(2006)]andIn compliance with the fractional differential order in FADE,αfalls into the interval(1,2].There are two common approaches to derivein(1):

·Numerical approach—numerically solve a Vandermonde system of linear equations formed by replacinguin(1)with monomials;

·Analytical approach—use the Lagrange interpolation ofuto derive the explicit formulas.

Using the former approach,the convergence of the DQM is adversely affected by the inherently ill-conditioned Vandermonde matrix as the number of sampling points goes up.The latter approach does not suffer the same pitfall and therefore is always employed for derivingwherenis an integer To obtainfor fractional-orderαin(3),we first consider the analytical approach.The Lagrange interpolation ofugives

where

Differentiation of(5)leads to

Comparison of(7)and(3)gives the explicit form of the weighting coefficients:

For the trivial caseα=2 andθ=0,the explicit expression of(8)has been given in[Shu and Richards(1992)].Since the fractional Leibniz differential rule is complicated[Podlubny(1999)],the deduction of(8)forα∈(1,2)is cumbersome.

A seem ingly alternative for determining(8)is the termwise differentiation onlk(x)which reads

whereis the coefficient before the monomialxp-1in the expansion oflk(x).Ascan be determined easily,the remaining work is to derive the coefficients.This can be done by the Matlab functionpoly()which returns the coefficients of a polynomial with zeros{xi}.Despite this convenience,our numerical results on 1D problems show that the solution accuracy of DQM based on(9)drops substantially when the number of equispaced sampling points grows,i.e.Nxexceeds 18(or larger).The convergence is deteriorated by the numerical difficulty of a very largetimes a very smallxp-1for a largeNxin(9).

Noting the limitation of the direct analytical approach,we consider the numerical approach The weighting coefficients are determined by the following linear systems(α∈(1,2))

Using the basic properties of Caputo fractional derivative the lefthand side of the above equations can be explicitly determined as follows[Podlubny(1999)]

From our numerical tests,the weighting coefficients of(3)obtained by solving the linear systems lead to higher convergence than evaluating(9)directly.Thus,the numerical approach will be adopted in the rest of the paper to solve the weighting coefficients.

2.2 DC formula

To consider a rectangular domain Ω=[a,b]×[c,e],(4)is extended to two dimensions as:

in which the fractional directional integral is defined as[Samko,Kilbas,and-Marichey(1993);Pang,Chen,and Sze(2013)],forσ∈(0,1),

Here,is the second-order directional derivative.The upper integration limit of the integral in(15),i.e.d,is the “backward distance”of internal node(x,y)to the boundary of the problem domain Ω(see Fig.1).This quantity is determined by the node location(xy),the direction vectorθ,and the shape of the computational domain.

Figure 1:The “backward distance”d of the internal node(x,y)to the boundary of domain Ω in the direction θ .

The DC formula of fractional directional derivative can be written as

and the corresponding DC weighting coefficients are determined through solving the following linear system:for

The subscriptsi,jin(16)and(17)are both global subscripts which correspond to a list of grid(or sampling)points(xkyk)sorted by a given order.The lefthand side of(17)can be computed by the Gauss-Lobatto-type quadrature rules discussed in our previous work[Pang,Chen and Sze(2013)].

2.3 Improved DC formula

To obtain the weighting coefficientsin(17),one needs to solveNxNylinear systems of dimensionNxNy.Here,we introduce a more efficient procedure which only needs to solve one linear system of dimensionNxorNy.

Figure 2:The DQ sampling points(“·”)contributing to the approximation of fractional directional derivative evaluated at point P

For simplicity,we consider a rectangular domain modeled by a 4×4 grid as shown in Fig.2 and the backward distance of pointPstarts from the pointEon the lower boundary in directionθ.A longEP,there areNDQsampling points(ˉxk,ˉyk)s which will be termed as the DQ sampling points.The fractional directional derivative at pointPis first approximated by a weighted linear sum of functional valuesˉuks at the DQ sampling points:agrid,the highest order monomial term isA long a straight line such asEP,xandyare linearly related and the highest order term can be expressed aswhich can be exactly interpolated atNx+Ny-1 DQ sampling points.It thus requiresNDQ=7 for the grid as shown in Fig.2 wherecan be obtained bywhereand

whereis the DQ weighting coefficients.Using Lagrange interpolation for the reference DQ weighting coefficientis solved from(10)by lettinga=0,Nx=NDQand introducing the normalized coordinate?xk∈[0,1]for the reference sampling points.The ordinates of the reference DQ sampling points alongEPare related to the normalized coordinate?xk∈[0,1]throughˉyk=yE+(yP-yE)?xiwhereyEandyPare the ordinates of pointEand pointP,respectively,as seen in Fig.3.

Figure 3:Auxiliary points“□”onto which the DQ weighting coefficients are distributed.“·”denotes the reference sampling point.

Conventionally,one can relate the values ofuat the DQ sampling points to those at the grid points by the 2D Lagrange interpolation which involvesNxNyterms.A more efficient procedure proposed here is to use a two-step 1D interpolation.Considering thei-th DQ sampling point with ordinateˉyiin Fig.3 as an example,the auxiliary points are the points with the same ordinateon theNxvertical grid lines.In the first step,the values of theuat auxiliary points are first obtained from those at theNygrid points on the same vertical grid line by the 1DNy-point Lagrange interpolation.In the second step,the values ofuat theNxauxiliary points are interpolated to thei-th DQ sampling point using the 1DNx-point Lagrange interpolation.

Finally,the DQ weighting coefficients at sampling points on the line PE are distributed onto all the grid points thru the two-step Lagrange interpolation,thus leading to the modified DC weighting coefficientswhich are similar toin(16).In fact,becauseandcan both be represented in terms of the fractional directional derivative of Lagrange interpolationli(x)lj(y),these two types of weighting coefficients are identical in theory,and therefore the quadrature solutions of DCM and improved DCM should have been completely identical.However,numerical errors enter differently toandthe pertinent quadrature solutions are different especially when the number of grid points is large.This point will be illustrated in Section/Example 4.2.

Fig.4 shows the auxiliary points when the starting point of the backward distance is on different portion of the domain boundary.In subsequent descriptions,the improved DCM will be abbreviated as iDCM.

Figure 4:Selections of auxiliary points for other cases of

3 DQ(-DC)formulas for fractional diffusion term

The 2D steady-state FADE can be written as[Meerschaert,Benson and Baumer(1999);Roop(2006);Pang,Chen and Sze(2013)]

with fractional diffusion operatorconstant diffusion factorc,constant medium velocityb,source termfand the rectangular domain Ω.in R2is defined as[Meerschaert,Benson and Baumer(1999)]

with the probability measureM(dθ)which can be discrete or continuous[Roop(2006);Pang,Chen and Sze(2013)]:

The directional differential operatorhere belongs to the Riemann-Liouville type.For computational convenience,we consider its Caputo counterpartAccordingly,we replace the fractional diffusion operator in(19)with the following operator

For a one-dimensional operator,(21)reduces to

and from(3)its DQ formula is given by

For a two-dimensional continuous operator in(22),numerical quadrature can be used to approximate(23),i.e.

whereθkandwkare the quadrature point and weight,respectively.The corresponding DC and improved DC formulas are

wherei=1,2,...NxNy.Similar formulas can also be found for discrete operator in(21).For convenience,we use compound trapezoidal rule here by lettingandθk=2kπ/NwhereNis the number of integration intervals with respect toθ

4 Numerical results and discussions

4.1 One-dimensional advection-diffusion problem

with exact solution

(a)u(x)=x8x5+3x2+5;

(b)u(x)=excos(x)+x3

Here,is prescribed to beand the Dirichlet boundary condition is taken;the source termf(x)is obtained by the10-point Gauss-Jacobi-Lobatto quadrature as given in[Pang,Chen and Sze(2013)].The first-order derivative in(28)is approximated by conventional DQ formula,i.e.

where the weighting coefficients are given by[Shu and Richards(1992)]:

We use DQ formula(25)to approximatechoose five test pointsz=0.1,0.3,0.5,0.7 and 0.9,and define the normalized error as

Figure 5:Convergence of the DQ solutions of the Example 4.1(a)

Same as in integer-order derivative problems,it can be seen that DQM can achieve the spectral accuracy for fractional derivative problem and unequally spaced sampling points can delay the numerical errors.This is because the Vandermenda matrix system for determining DQ weighting coefficients generated by non-uniform sampling points is less ill-conditioned than that generated by uniform sampling points.Despite this advantage,the machine precision～10-15for double precision starts to lose atNx≈25.On the other hand,for a numerical solution method having spectral accuracy,a less dense grid is often adequate for most problems..

Figure 6:Convergence of the DQ solutions of Example 4.1(b)

4.2 Two-dimensional homogenous diffusion problem

with exact solution

(a)u(xy)=x2(1-x)2y2(1-y)2

(b)u(xy)=x3cos(x2+y2)+ex+2y

The Dirichlet boundary condition is taken.To compute the value offat each grid point,integration quadrature is implemented.The backward distancedin(33)can be derived by[Roop(2006)]:

The integral with respect toθis partitioned into five sub-integrals according to(34).High-order Gauss-Jacobi-type quadrature rules[Pang,Chen and Sze(2013)]are used to evaluate each sub-integral for a highly accurate source termf.

Here,81 test points with coordinates(xiyj)=(i/10,j/10),i,j=1,2,...,9 are employed to quantify the normalized error which is defined as

where the approximant?uijare derived by 2D Lagrange interpolation at the gird points.

We first consider the case(a)with exact solutionx2(1-x)2y2(1-y)2.Roop’s FEM in[Roop(2006)]is taken as a comparison.Four-node rectangular elements and uniform ly spaced nodes are employed In iDCM,if the grid points are equally(unequally)distributed,the DQ sampling points would also be equally(unequally)distributed.Other combinations of grid point and DQ sampling point distributions are possible but the numerical solutions are plagued by numerical errors more readily as the grid density increases.Fig.7 compares the accuracy of FEM,DCM and iDCM.LetN=64 in(26)for DCM and iDCM.Fig.8 compares the CPU time(using an Intel i5 dual-core CPU and 2G RAM)for computing the weighting coefficient matrix of(improved)DCM with the CPU time for generating the global stiffness matrix in FEM.Here,Nx=Ny=Mwhich vary from 4 to 20.Note that the inherent ill-conditioning of the system matrix in(17)simply allows a less dense grid(i.e.,a smallNxandNy)when one needs acceptable solution accuracy.

Figure 7:Convergence of the FEM,DCM and iDCM(NDQ=12)solutions of Example 4.2(a)

Figure 8:CPU time for generating discretized matrices of FEM,DCM and iDCM(NDQ=12)in Example 4.2(a)

It can be seen from Fig.7 that DCM and iDCM at smallNpossess the spectral convergence which is much higher than that of FEM.Moreover,the non-uniform grid yields more stable solutions.As discussed in Section 2.3,for smallNxorNy,DCM and iDCM yield identical solutions;while for largeNxorNy,their solutions exhibit obvious distinction.The computational efficiency of the improved DCM over DCM can be seen in Fig.8

It has been discussed in Section 2.3 that for aNx×Nygrid,the number of DQ sampling points,i.e.,NDQ,in the iDCM is sensible to takeNx+Ny-1 based on the interpolation consideration.On the other hand,ifNDQis too large,the Vandermenda matrix may be ill-conditioned and adversely affect the accuracy and convergence.With reference to Example 4.2(a),the exact solution can be exactly interpolated by a 5×5 grid whenNDQ≥9.This can be seen by an obvious increase of the accuracy forNx(=Ny)≥5 andNDQ≥9 as shown in Fig.9.

Figure 9:Convergence of the improved DCM for different NDQ with unequally spaced gird points in Example 4.2(a)

Next,we consider case(b)with exact solutionx3cos(x2+y2)+ex+2y.Since the finite element scheme for the 2D FADE(19)is only developed for a problem with homogenous Dirichlet boundary condition[Roop(2006)],we will not consider the FEM solution for(b)which does not vanish identically over the domain boundary.The comparison of the solutions from the DCM and the iDCM(NDQ=12)forN=64 in(26)is given in Fig.10.It can be seen that accurate and stable solutions can still be secured by the iDCM with a unequispaced grid.

Figure 10:Convergence of the DCM and iDCM(NDQ=12)solutions of Example 4.2(b)

Now discuss the influence ofNDQon the solution accuracy of the iDCM.Fig.11 displays the variation of the normalized error of the iDCM with theNDQ.The observations show that the influence ofNDQcan only be clearly seen for a largerNxorNyand an increasingNDQwill not obviously improve the solution accuracy.We conjecture that two factors,i.e.,(a)the inherently ill-conditioning of the linear system that determines the DQ weighting coefficients and(b)the numerical errors accumulated in the two-step 1D Lagrange interpolation deteriorate the solution convergence.From the interpolation consideration in Section 2.3,a possible way to find an optimalNDQis to letNDQvary with the gird size,namely,NDQ=Nx+Ny-1.The corresponding convergence curve is shown in Fig.12,from which one can see the poor accuracy for a largerNxorNy.Comparison of Figs.11 and 12,however,indicates that using a constantNDQis usually more suitable than using a variable one,but for a smallerNxorNy,both cases are available.Considering the influence of constantNDQas shown in Fig.11,we suggest a constantNDQ∈[7 m in{18,Nx+Ny-1}]forNxNy≥4

Figure 11:Convergence of iDCM for constant NDQ with unequally spaced grid points in Example 4.2(b)

Figure 12:Convergence of iDCM for variable NDQ=Nx+Ny-1 with unequally spaced grid points in Example 4.2(b).

It is worth noting that the integration quadrature accuracy in(26)also plays a crucial role in the final solution accuracy.This will be true for both the case(a)and case(b).For instance,the variation of normalized errorErrwith the interval numberNin compound trapezoidal rule for case(b)whereNx=Ny=10 andNDQ=16∈[7,18]is illustrated in Fig.13.The solution accuracy will gradually increase with the increasingN.Consideration of both the accuracy and the CPU time leads to an intermediate but acceptableN.Here,we tookN=64 just as mentioned before

Thru numerical tests,it is interesting to find that using higher-order quadrature rules such as compound Simpson and even the(piecewise)Gauss quadrature rules will not obviously improve the convergence in Fig.13.The is also the reason why we use compound trapezoidal rule instead of high-order rules to discretize the angle integral in(26)and(27).

Figure 13:Variation of the normalized error with the number of the integration interval in compound trapezoidal rule using nonuniform grid for Example 4.2(b)

5 Conclusions

This paper presents differential quadrature and cubature methods for numerical solution of steady-state advection-diffusion equations with directionally-weighted fractional diffusion operator.The present methods feature spectral convergence and enjoy lower computational effort compared to the finite element method.An improved version of the original differential cubature method is proposed to reduce the computational cost considerably On the other hand,the weighting coefficients derived from solving the inherently ill-conditioned linear system fail these methods in large-scale problems.A partial antidote to this problem is to deduce explicit expressions for DQ/DC weighting coefficients that take the similar form to(30),which is under study now.

Acknow ledgement:The research is financially supported by the Hong Kong Research Grant Council(GRF grant HKU 7173 09E),the 111 project under grant B12032,the National Basic Research Program of China(973 Project No.2010CB832702),and the National Science Funds for Distinguished Young Scholars(11125208).

Benson,D.A.;Schumer,R.;M eerschaert,M.M.;Wheatcraft,S.W.(2001):Fractional diffusion,Levy motions,and the MADE tracer tests.Transp.Porous Media,vol.42,pp.211-240.

Benson,D.A.;Wheatcraft,S.;Meerschaert,M.M.(2000):Application of a fractional advection-diffusion equation.Water Resour.Res.,vol.36,pp.132-138.

Bert,C.W.;Malik,M.(1996):Differential quadrature method in computational mechanics:A review.Appl.Mech.Rev.,vol.49,pp.1-28.

Chen,W.;Holm,S.(2004):Fractional Laplacian time-space models for linear and nonlinear lossy media exhibiting arbitrary frequency power-law dependency.J.Acoust.Soc.Am.,vol.115,pp.1424-1430.

Chen,Y.;Sun,L.;Li,X.;Fu,X.(2013):Numerical solution of nonlinear fractional integral differential equations by using the second kind Chebyshev wavelets.CMES:Computer Modeling in Engineering&Sciences,vol.90,pp.359-378.

Chen,W.;Zhong,T.;Shu,C.(1998):A Lyapunov formulation for efficient solution of the Poisson and convection-diffusion equations by the differential quadrature method,J.Comput.Phys.,vol.141,pp.78-84.

Duan,J.S.(2005):Time-and space-fractional partial differential equations.J.Math.Phys.,vol.46,pp.13504-13511.

Ervin,V.J.;Roop,J.P.(2006):Variational solution of fractional advection diffusion equation on bounded domains in Rd.Numer.Meth.Part Differ.Equ.,vol.23,pp.256-281.

Hanert,E.(2010):A comparison of three Eulerian numerical methods for fractional-order transport models.Environmental Fluid Mechanics,vol.10,pp.7-20.

Huang,Q.Z;Huang,G.H.;Zhan,H.B.(2008):A finite element solution for the fractional advection-diffusion equation.Advances in Water Resources,vol.31,pp.1578-1589.

Ilic,M.;Liu,F.W.;Turner,I.;Anh,V.(2005):Numerical approximation of a fractional-in-space diffusion equation,I.Fractional Calculus&Applied Analysis,vol.8,pp.323-341.

Ilic,M.;Liu,F.W.;Turner,I.;Anh,V.(2006):Numerical approximation of a fractional-in-space diffusion equation(II)—with nonhomogeneous boundary conditions.Fractional Calculus&Applied Analysis,vol.9,pp.333-349.

Katzav,E.;Adda-bedia,M.(2008):The spectrum of the fractional Laplacian and First-Passage-Time statistics.Europhysics Letters,vol.83,pp.30006.

Kilbas,A.A.;Srivastava,H.M.;Trujillo,J.J.(2006):Theory and applications of fractional differential equations.Elsevier B.V..

Li,C.P.;Chen,A.;Ye,J.J.(2011):Numerical approaches to fractional calculus and fractional ordinary differential equation.Journal of Computational Physics,vol.230,pp.3352-3368.

Li,X.J.;Xu,C.J.(2010):Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and spectral method approximation.Commun.Comput.Phys.,vol.8,pp.1016-1051.

Lui,S.H.(2011):Numerical analysis of partial differential equations.WELEY.

Malik,M.;Civan,F.(1995):A comparative study of differential quadrature and cubature methods vis-à-vis some conventional techniques in context of convection diffusion-reaction problems.Chemical Engineering Science,vol.50,pp.531-547.

Meerschaert,M.M.;Benson,D.A.;Baumer,B.(1999):Multidimensional advection and fractional diffusion.Phys.Rev.E,vol.59,pp.5026-5028.

Meerschaert,M.M.;Mortenson,J.(2004):Vector Grünwald formula for fractional derivatives.Fract.Calc.Appl.Anal.,vol.7,pp.61-81.

Meerschaert,M.M.;Scheffler,H.P.;Tad jeran,C.(2006):Finite difference methods for two-dimensional fractional diffusion equation.J.Comput.Phys.,vol.211,pp.249-261.

Mokhtari,R.(2011):GDQM for solving the fractional diffusion equation,Communications in Fractional Calculus,vol.2,pp.1-13.

Pang,G.F.;Chen,W.;Sze,K.Y.(2013):Gauss-Jacobi-type quadrature rules for fractional directional integrals.Computation and Mathematics with Applications,vol.66,pp.597-607.

Podlubny,I.(1999):Fractional differential equations:An introduction to Fractional Derivatives,Fractional Differential Equations,to methods of their Solutionand some of their Applications.ACADEM IC PRESS.

Quan,J.R.;Chang,C.T.(1989):New insights in solving distributed system equations by the quadrature method-I.Analysis.Computers Chem.Engng.,vol.13,pp.779-788.

Roop,J.P.(2006):Computational aspects of FEM approximation of fractional advection diffusion equations on boundary domains in R2.J.Comput.Appl.Math.,vol.193,pp.243-268.

Samko,S.G.;Kilbas,A.A.;Marichev,O.I.(1993):Fractional integrals and derivatives:Theory and Applications.Gordon and Breach,Newark,N.J..

Shu,C.(2000):Differential quadrature and its application in engineering.Springer-Verlag London Limited.

Shu,C.;Chew,Y.T.(1997):Fourier expansion-based differential quadrature and its application to Helmholtz eigenvalue problems.Commun.Numer.Meth.Engrg.,vol.13,pp.643-653.

Shu,C.;Ding,H.;Yeo,K.S.(2003):Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations.Comput.Methods Appl.Mech.Engrg.,vol.192,pp.941-954.

Shu,C.;Richards,B.E.(1992):Application of generalized differential quadrature to solving two-dimensional incompressible Navier-Stokes equations.International Journal for Numerical Methods in Fluids,vol.15,pp.791-798.

Shu,C.;Yao,Q.;Yeo,K.S.(2002):Block-marching in time with DQ discretization:an efficient method for time-dependent problem.Comput.Methods Appl.Mech.Engrg.,vol.191,pp.4587-4597.

Tad jeran,C.;Meerschaert,M.M.(2007):A second-order accurate numerical method for the two-dimensional fractional diffusion equation.J.Comput.Phys.,vol.220,pp.813-823.

Wang,L.;M eng,Z.;M a,Y.;Wu,Z.(2013):Numerical solution of fractional partial differential equations using Haar wavelets.CMES:Computer Modeling in Engineering&Sciences,vol.91,pp.269-287.

Wang,H.;Wang,K.X.(2011):An O(Nlog2N)alternating-direction finite difference method for two-dimensional fractional diffusion equations.J.Comput.Phys.,vol.230,pp.7830-7839.

Yang,Q.Q.;Liu,F.W.;Turner,I.(2010):Numerical methods for fractional partial differential equations with Riesz space fractional derivatives.Applied Mathematical Modeling,vol.34,pp.200-218.

Yang,Q.Q.;Turner,I.;Liu,F.W.;Ilic,M.(2011):Novel numerical methods for solving the time-space fractional diffusion equation in two dimensionsSIAM J.SCI.COMPUT.,vol.33,pp.1159-1180.

Zhang,Y.;Benson,D.A.;Reeves,D.M.(2009):Time and space nonlocalities underlying fractional-derivatives models:Distinction and literature review of field applicationsAdvances in Water Resources,vol.32,pp.561-581.

Zhang,X.;Craw ford,J.W.;Deeks,L.K.;Stutter,M.I.;Bengough,A.G.;Young,I.M.(2005):A mass balance based numerical method for the fractional advection diffusion equation:theory and application.Water Resour.Res.,vol.41,pp.W 07029.

Zeng,F.H.;Li,C.P.;Liu,F.W.;Turner,I.(2013):The use of finite difference/element approaches for solving the time-fractional subdiffusion equation.SIAM J.Sci.Comput.,vol.35,pp.A2976-A3000.

Zheng,Y.Y.;Li,C.P.;Zhao,Z.G.(2010):A note on the finite element method for the space-fractional advection diffusion equation.Computers and Mathematics with Applications,vol.59,pp.1718-1726.