ShuyingZhai,DongweiGui,JianpingZhaoandXinlongFeng,?

1School of Mathematics Science,Huaqiao University,Quanzhou,Fujian 362011,P.R.China.

2Cele National Station of Observation&Research for Desert Grassland Ecosystem,Xinjiang Institute of Ecology and Geography,Chinese Academy of Sciences,Urumqi,Xinjiang 830011,P.R.China.

3College of Mathematics and Systems Science,Xinjiang University,Urumqi,Xinjiang 830046,P.R.China.

High Accuracy Spectral Method for the Space-Fractional Diffusion Equations

ShuyingZhai1,DongweiGui2,JianpingZhao3andXinlongFeng3,?

1School of Mathematics Science,Huaqiao University,Quanzhou,Fujian 362011,P.R.China.

2Cele National Station of Observation&Research for Desert Grassland Ecosystem,Xinjiang Institute of Ecology and Geography,Chinese Academy of Sciences,Urumqi,Xinjiang 830011,P.R.China.

3College of Mathematics and Systems Science,Xinjiang University,Urumqi,Xinjiang 830046,P.R.China.

.In this paper,a high order accurate spectral method is presented for the space-fractional diffusion equations.Based on Fourier spectral method in space and Chebyshev collocation method in time,three high order accuracy schemes are proposed.The main advantages of this method are that it yields a fully diagonal representation of the fractional operator,with increased accuracy and efficiency compared with low-order counterparts,and a completely straightforward extension to high spatial dimensions.Some numerical examples,including Allen-Cahn equation,are conducted to verify the effectiveness of this method.

AMS subject classifications:35K55,65M70,65L06,65L12

Chinese Library classifications:O241.82

Space-fractional diffusion equation,fractional Laplacian,Chebyshev collocation method,Fourier spectral method,implicit-explicit Runge-Kutta method.

1 Introduction

Fractional differential equations have been proved to be valuable tools in modeling of many phenomena in various fields.In water resources,fractional models provide a useful description of chemical and contaminant transport in heterogeneous aquifers[1,2].In transport dynamics,they have been used to describe transport dynamics in complexsystems which are governed by anomalous diffusion and non-exponential relaxation patterns[3].Moreover,they are also used in finance,engineering and physics(see[4–6]and references cited therein).

In this paper,we consider the following space fractional diffusion equation

with the homogeneous Dirichlet or homogeneous Neumann boundary conditions.Here K ＞0 is the conductivity or diffusion tensor,and(?Δ)α/2is the fractional Laplacian operator[7]with 1＜α＜2.The function f=f(x,t,u)denotes the nonlinear source term.

There are many numerical methods to discretize the fractional Laplacian operator of problem(1.1).However,fractional differential operator is non-local red operator,which generates computational and numerical difficulties that have not been encountered in the context of the classical second-order diffusion equations.For space-fractional diffusion equations,numerical methods often generate full coefficient matrices with complicated structures[8–11].In this paper we use Fourier spectral methods[12-14]to discretize the space-fractional derivative.This approach gives a full diagonal representation of the fractional operator and achieves spectral convergence regardless of the fractional power in the problem.Meanwhile,the application to high spatial dimensions is the same as the one-dimensional problem.For the temporal discretization,based on Chebyshev nodes[15,16],the second-order Crank-Nicolson(CN)method and third-order implicitexplicit(IMEX)Runge-Kutta method[17]are used on the Chebyshev grids,respectively.Numerical experiments in Section 3 show that the time accuracy using Chebyshev grids is more accurate than using uniform grids.

The outline of this paper is as follows.In Section 2,three collocation/spectral numerical schemes are given for the space fractional diffusion equation(1.1).In Section 3,three numerical examples are carried out to verify the high efficiency of the proposed method,including the space-fractional Allen-Cahn equation in two dimensions.Finally,conclusions are drawn in Section 4.

2 High-order accurate schemes

In this section,we present three numerical schemes to simulate the asymptotic behavior of solution for the space fractional diffusion equation(1.1).The proposed schemes are based on Fourier spectral method in space and the collocation technique in time.In order to simplify the notations and without lose of generality,we only present numerical schemes for the one-dimensional space-fractional diffusion equation.

2.1 Fourier spectral spatial discretization

This subsection starts to present high-order accurate spectral method approximating the initial boundary value problem(1.1).The symbol(?Δ)α/2has the usual meaning as a function of Laplacian(?Δ),which is defined in terms of its spectral decomposition.In order to illustrate the main idea of the proposed method,the following definition is adopted.

Definition 2.1.Suppose the one-dimensional Laplacian(?Δ)has a complete set of orthonormal eigenfunctions ?icorresponding to eigenvalues λion the bounded region[a,b],i.e.,(?Δ)?i=λi?i.Let

Then for any u∈Uα,the Laplacian(?Δ)α/2is de fined by

where λiand ?iwill depend on the speci fied boundary conditions:

(1)Homogeneous Dirichlet boundary condition

(2)Homogeneous Neumann boundary condition

From Definition 2.1,we know thatcan be used to approximate the exact

solution u(x),where N is a positive integer.

Meanwhile,by combined with Eq.(2.1),the i-th Fourier mode of Eq.(1.1)becomes

Remark 2.1.For the one-dimensional problems,it is widely assumed that fractional Laplacian operatoris equivalent to the Riesz fractional derivativeunder homogeneous Dirichlet boundary conditions[24].However,it is difficult to extend this method to Caputo or Riemann Liouville derivative.For both kinds of fractional derivatives,much effort has been devoted to develop high order spectral methods,such as least square spectral method[25],spectral collocation method[26,27].The research on these aspects will be reported in our future work.

2.2 Chebyshev collocation time discretization

It is well known that Chebyshev points are the best for Cauchy optimality and they are very powerful in high order polynomial approximation.Continuous functions defined on[-1,1]can be approximated very accurately by using the polynomial interpolation with enough Chebyshev points.Above all,the rate of convergence of a scheme can be accelerated using Chebyshev points.In this subsection,we use Chebyshev points to discrete the time variable,and three high order accuracy schemes are proposed based on CN method[21,22]and IMEX Runge-Kutta method with the non-uniform time step size.

Chebyshev polynomials are a well known family of orthogonal polynomials that have many applications[19,20].They are define don the interval[?1,1]and relatedre cursively by

where T0(z)=1 and T1(z)=z.

The Chebyshev nodes zmof degree M are the zeros of TM,namely

For practical use of the Chebyshev nodes on the time interval of interest t∈[0,T],it is necessary to shift these nodes by the following relationship:

Then the Chebyshev gridsT={tm|0≤m≤M+1}with t0=0,tM+1=T are given correspondingly.

Based on the time interval[tm,tm+1],a CN type finite difference scheme for Eq.(2.2)can be given as follows:

However,for the nonlinear case,we need use an iterative method to solve the resulting algebraic system.In order to overcome this drawback,the following linearization scheme is needed

where τ=max0≤m≤M?1(tm+1?tm).

Then,it follows from scheme S1 that

Scheme S2 can be seen as an IMEX scheme,which is a combination of second-order Adams-Bashforth scheme for the explicit term f(t,u)and CN scheme for the implicit term(?Δ)α/2u.

Note that both schemes are only second-order accuracy in time.In order to construct higher-order scheme,we will investigate the performance of IMEX Runge-Kutta scheme for the stiff and non stiff terms.The s-stage IMEX Runge-Kutta scheme from tm?1to tmcan be generally represented as:

with internal stages given by

where δτ=tm?tm?1,and g1is the implicitly treated part,while g2is the explicitly treated part.

IMEX Runge-Kutta schemes can be represented concisely by two Butcher tableaus[17],and notationscanbede fined similarly.Moreover,the coefficients c and e c are given respectively by the usual relation

where

Now we use a third-order accurate IMEX Runge-Kutta scheme

to solve Eq.(2.2).The scheme,which will be called S3,consists of applying an implicit discretization forand an explicit one for

Remark 2.2.Note that the two-dimensional case can be handled trivially in the previous formulations by simply replacing λrespectively,where

and the orthonormal eigenfunctions ?ijcorresponding to eigenvalues λi+λjin a rectangular region[a,b]2.

Remark 2.3.It is well known that the CN scheme is an unconditionally stable,implicit scheme with second-orderaccuracy in time[21,22],i.e.,scheme S1 in this paper is unconditionally stable.Nevertheless,scheme S2 is conditionally stable,which has a reasonable time step restriction for largerKand small space step[23].For scheme S3,we know from the discussion given in[17]that it is L-stable.

3 Numerical experiments

In this section,three numerical examples are presented to demonstrate the efficiency and accuracy of the proposed method.We compute the maximum norm errors

for one-and two-dimensional cases,respectively.We also compute the temporal convergence order

The order of accuracy is formally defined when the mesh size approaches to zero.Therefore,when M is relatively small,the numerical scheme may not achieve its formal order of accuracy.

Comparison with the related work[28]is presented to show the effectiveness of the proposed method.Meanwhile,numerical results on the uniform time step sizes are also provided.

Problem 1

In order to compare our schemes with Bueno-Orovib et al.’scheme[28](denoted by S0),we first consider a one-dimensional problem with homogeneous Dirichlet boundary conditions in their paper.The exact analytical solution and the corresponding force term inx∈(0,1)are given by

The data in Tables 1 and 2 show the maximum norm errors for the numerical solution withα=1.5,K=10,N=51 andT=1.From both tables we find that the numerical results on Chebyshev grids are much better than those on uniform grids in time.As predicted,all of them generate the correspondingly temporal convergence orders whenMis large enough.

Meanwhile,numerical results in Table 1 show that schemes S1-S3 are more accurate than S0 whenM≥1000,and scheme S1 can obtain the highest accuracy.Of course,this scheme need iteration.In fact,we find four or less iterations are sufficient to obtain high precise numerical solution in our computation.Moreover,although scheme S3 is thirdorder accurate in time,it does not show a selective superiority untilM=4000.Asa whole,scheme S1 would be the best choice.

Table 1:Numerical results of Problem 1 using Chebyshev collocation points in time at T=1,N=51 and α=1.5 withK=10.

Table 2:Numerical results of Problem 1 using uniform grids in time at T=1,N=51 and α=1.5 withK=10.

Problem 2

To better illustrate the efficiency of the proposed method,we extend Problem 1 to two dimensional case.The exact solution and the corresponding force term in(x,y)∈(0,1)2are given by

In this test,we fixα=1.5,K=10 andN=51.The numerical results atT=1 are presented in Tables 3 and 4.Again,the data in Table 3 are more accurate than those in Table 4,which further con firm that scheme S1 is the best choice.

Table 3:Numerical results of Problem 2 using Chebyshev collocation points in time at T=1,N=51 and α=1.5 withK=10.

Table 4:Numerical results of Problem 2 using uniform grids in time at T=1,N=51 and α=1.5 withK=10.

Problem 3

The Allen-Cahn equation represents a model for anti-phase domain coarsening in a binary mixture.The continuous problem has a decreasing total energy[29–31].Now we consider the following space-fractional Allen-Cahn equation

with homogeneous Neumann boundary conditions,and the initial conditions are

and

for one-and two-dimensional cases,respectively.

In this test,settingK=0.01,N=100 and τ=1.All numerical results are obtained by scheme S1.Fig.1(a-c)show the time evolution of the one-dimensional Allen-Cahn equation for varying α.Fig.1(a)shows that the initial datum evolves to an intermediate unstable equilibrium,followed by a rapid transition to stable state of u=±1.As the fractional power is decreased,Fig.1(b)shows the lifetime of the unstable interface is largely prolonged,eventually becoming fully stable due to the long-tailed influence of the fractional diffusion process(Fig.1(c)).Correspondingly,Fig.1(d-e)show the trend of energy evolution E(u),which can be written as

Figure 1:Numerical solution and corresponding energy of the one-dimensional space-fractional Allen-Cahn equation for varying α.

Figure 2:The contour plots of numerical solution of the two-dimensional space-fractional Allen-Cahn equation for varying α.

Figs. 2 and 3 show the contour plots and corresponding energy of the two-dimensional space-fractional Allen-Cahn equation for varying a. The same conclusions as the one-dimensional case can be obtained.

Figure 3:The energy of the two-dimensional space-fractional Allen-Cahn equation for varying α.

4 Conclusions

In this work,three numerical schemes for solving space-fractional diffusion equation are proposed based on Fourier spectral method in space and collocation method in time.Numerical experiments have shown that the numerical results on Chebyshev grids are more accurate than those on the uniform grids in time,and scheme S1 may be the best choice in this work.Meanwhile,although scheme S3 has third-order convergence in time,the precision advantage would not be displayed unless the Chebyshev nodes M is big enough.So scheme S3 is not applicable to long time behavior,such as Allen-Cahn equation.Moreover,the delay effects[18]of fractional operator are also con firmed.

The authors would like to thank the editor and referees for their valuable comments and suggestions which helped us to improve the results of this paper.This work is in parts supported by the Distinguished Young Scholars Fund of Xinjiang Province(No.2013711010),the Western Light Program of Chinese Academy of Sciences(No.XBBS201105),and the NSF of China(No.11271313,No.61163027,No.41471031).

[1]D.A.Benson,S.Wheatcraft and M.M.Meerschaert.Application of a fractional advection-dispersion equation.Water Resources Res.,36:1403–1412,2000.

[2]M.M.Meerschaert,D.A.Benson and S.W.Wheat craft.Subordinated advection-dispersion equation for contaminant transport.Water Resource Res.,37:1543–1550,2001.

[3]R.Metzler and J.Klafter.The random walk’s guide to anomalous diffusion:a fractional dynamics approach.Phys.Rep.,339:1–77,2000.

[4]W.H.Deng and C.P.Li.The evolution of chaotic dynamics for fractional unified system.Phys.Lett.A,372:401–407,2008.

[5]S.Y.Zhai,X.L.Feng and Z.F.Weng.New high-order compact ADI algorithms for 3D nonlinear time-fractional convection-diffusion equation.Math.Probl.Engrg.,http://dx.doi.org/10.1155/2013/246025.

[6]S. Y. Zhai, X. L. Feng and Y. N. He. An unconditionally stable compact ADI method for three dimensional time-fractional convection-diffusion equation. J. Comput. Phys., 269: 138–155,2014.

[7]Q.Q.Yang,I.Turner,F.W.Liu and M.Ili’c.Novel numerical methods for solving the time space fractional diffusion equation in 2D.SIAM J.Sci.Comp.,33:1159–1180,2011.

[8]M.M.Meerschaert and C.Tadjeran.Finite difference approximations for fractional advection-dispersion flow equations.J.Comput.Appl.Math.,172:65–77,2004.

[9]V.J.Ervin,N.Heuer and J.P.Roop.Numerical approximation of a time dependent,nonlinear,space-fractional diffusion equation.SIAM J.Numer.Anal.,45:572–591,2007.

[10]H.K.Pang and H.W.Sun.Multigrid method for fractional diffusion equations.J.Comput.Phys.,231(2):693–703,2012.

[11]H.Wang and N.Du.A super fast-preconditioned iterative method for steady-state space-fractional diffusion equations.J.Comput.Phys.,240:49–57,2013.

[12]Y.M.Lin and C.J.Xu.Finite difference/spectral approximations for the time-fractional diffusion equation.J.Comput.Phys.,225:1533–1552,2007.

[13]X.J.Li and C.J.Xu.A space-time spectral method for the time fractional diffusion equation.SIAM J.Numer.Anal.,47:2108–2131,2009.

[14]J.Shen,T.Tang and L.Wang.Spectral Methods:Algorithms,Analysis and Applications,volume 41 of Springer Series in Computational Mathematics.Springer,2011.

[15]J.C.Mason and D.C.Handscomb.Chebyshev Polynomials.Chapman and Hall/CRC,2003.[16]X.L.Feng,T.Tang and J.Yang.Long time numerical simulations for phase- field problems using p-adaptive spectral deferred correction methods.Preprint,2013.

[17]U.M.Ascher,S.J.Ruuth and R.J.Spiteri.Implicit-explicit Runge-Kutta methods for time dependent partial differential equations.Appl.Numer.Math.,25:151–167,1997.

[18]K.Seki,M.Wojcik and M.Tachiya.Fractional reaction-diffusion equation.J.Chem.Phys.,119(4):2165–2170,2003.

[19]J.Shen.efficient Spectral-Galerkin method II.Direct solvers of second-and fourth-Order equations using Chebyshev polynomials.SIAM J.Sci.Comput.,16(1):74–87,1995.

[20]K.Maleknejad,S.Sohrabi and Y.Rostami.Numerical solution of nonlinear Volterra integral equation of the second kind by using Chebyshev polynomials.Appl.Math.Cmput.,188:123–128,2007.

[21]J.C.Strikwerda.Finite Difference Schemes and Partial Differential Equations.Pacific Grove,CA:Brooks/Cole,1989.

[22]J.W.Thomas.Numerical Partial Differential Equations:Finite Difference Methods.New York:Springer-Verlag,1998.

[23]U.Ascher,S.Ruuth and B.Wetton.Implicit-explicit methods for time-dependent PDE’s.SlAM J.Numer.Anal.,32:797–823,1995.

[24]Q.Yang,F.Liu and I.Turner.Numerical methods for fractional partial differential equations with Riesz space fractional derivatives.Appl.Math.Model.,34:200–218,2010.

[25]A.R.Carella and C.A.Dorao.Least-Squares Spectral Method for the solution of a fractional advection-dispersion equation.J.Comput.Phys.,232:33–45,2013.

[26]M.M.Khader.On the numerical solutions for the fractional diffusion equation.Commun.Nonlinear Sci.Numer.Simul.,16:2535–2542,2011.

[27]M.Zayernouri and G.E.Karniadakis.Fractional Spectral Collocation Methods.SIAM J.Sci.Comput.,36(1):40–62,2014.

[28]A.Bueno-Orovio,D.Kay and K.Burrage.Fourier spectral methods for fractional-in-space reaction-diffusion equations.BIT Numer.Math.,DOI 10.1007/s10543-014-0484-2.

[29]X.B.Feng and A.Prohl.Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows.Numer.Math.,94:33–65,2003.

[30]X.L.Feng,T.Tang and J.Yang.Stabilized Crank-Nicolson/Adams-Bash forth schemes for phase field models.East Asian J.Appl.Math.,3:59–80,2013.

[31]X.L.Feng,H.L.Song,T.Tang and J.Yang.Nonlinear stability of the implicit-explicit methods for the Allen-Cahn equation.Inverse Problems&Imaging,7(3):679–695,2013.

6 January 2014;Accepted 13 July 2014

?Corresponding author.Email addresses:zhaishuying123456@163.com(S.Zhai),guidwei@163.com(D.Gui),zhaojianping@126.com(Z.Zhao),fxlmath@xju.edu.cn(X.Feng)