Elsevier

Journal of Computational Physics

Volume 256, 1 January 2014, Pages 109-117
Journal of Computational Physics

Preconditioned iterative methods for fractional diffusion equation

https://doi.org/10.1016/j.jcp.2013.07.040Get rights and content

Abstract

In this paper, we are concerned with numerical methods for the solution of initial–boundary value problems of anomalous diffusion equations of order α(1,2). The classical Crank–Nicholson method is used to discretize the fractional diffusion equation and then the spatial extrapolation is used to obtain temporally and spatially second-order accurate numerical estimates. Two preconditioned iterative methods, namely, the preconditioned generalized minimal residual (preconditioned GMRES) method and the preconditioned conjugate gradient for normal residual (preconditioned CGNR) method, are proposed to solve relevant linear systems. Numerical experiments are given to illustrate the efficiency of the methods.

Introduction

Fractional differential equations emerge from numerous topics such as turbulent flow [7], [27], chaotic dynamics of classical conservative systems [39], groundwater contaminant transport [4], [5], and applications in physics [28], finance [25], biology [18], and image processing [2]. Since there are very few fractional differential equations whose analytical closed-form solutions are available, the research on numerical methods for the solution of the equation has attracted many researchers, and a lot of research results have been obtained; see [3], [6], [11], [12], [13], [15], [16], [17], [19], [20], [21], [22], [23], [29], [30], [31], [32], [33], [34], [35] and references therein. The existence and uniqueness of the solution also have been studied, e.g., Wang and Yang proved the existence and uniqueness of the solutions to variable-coefficient space-fractional diffusion equations [38].

In this paper, we consider the following initial–boundary value problem for a one-dimensional space fractional diffusion equation (FDE):u(x,t)t=d+(x,t)αu(x,t)+xα+d(x,t)αu(x,t)xα+f(x,t),x(xL,xR),t(0,T],u(xL,t)=u(xR,t)=0,0tT,u(x,0)=u0(x),x[xL,xR], where 1<α<2, f(x,t) is the source term, and d±(x,t)0. In this paper, we use the Grünwald–Letnikov form [24] to define the left-sided and the right-sided fractional derivatives as follows:αu(x,t)+xαlimh0+1hαk=0(xxL)/hgk(α)u(xkh,t),αu(x,t)xαlimh0+1hαk=0(xRx)/hgk(α)u(x+kh,t), where a denotes the largest integer that is not larger than a, g0(α)=1, andgk(α)=(1)kk!α(α1)(αk+1),k=1,2,.

We now briefly review recent developments of numerical methods for the solution of FDEs. In 2004 and 2006, Meerschaert and Tadjeran proposed a shifted Grünwald discretization to approximate the FDE with a left-sided fractional derivative and the FDE with two-sided fractional derivatives, respectively [20], [21]. In 2006, Tadjeran, Meerschaert, and Scheffler proposed a second-order accurate numerical scheme by using the classical Crank–Nicholson method and the Richardson extrapolation scheme [32]. The above methods were proved to be unconditionally stable. Nevertheless, the above-mentioned methods, as well as other finite difference and finite element methods [11], [12], [13], yield full coefficient matrices, which result in O(N2) storage and O(N3) complexity.

Recently, Wang, Wang, and Sircar [36] discovered that the coefficient matrix of Meerschaert–Tadjeranʼs method [21] possesses a Toeplitz-like structure: it can be written as a sum of diagonal-multiply-Toeplitz matrices. Thus the storage requirement is reduced from O(N2) to O(N), and the matrix–vector multiplications can be done with O(NlogN) complexity by using a fast algorithm based on the fast Fourier transform (FFT) [9], [10], [14]. By using the Toeplitz-like structure, Wang et al. proposed a fast finite difference scheme by approximating the coefficient matrix with a banded matrix of bandwidth O(logN). As a result, the computational cost per time step is only O(Nlog2N), which is much lower than other methods. However, the unconditional stability of their scheme has not been proved.

Most recently, iterative solvers for the discretized system of the FDE by Meerschaert–Tadjeranʼs method have been proposed and studied [23], [37]. Wang and Wang, two authors of [36], proposed the conjugate gradient normal residual (CGNR) method [37] to solve the discretized system. Their numerical results show that the CGNR method converges fast if the diffusion coefficients d±(x,t) are very small (in that case the coefficient matrix is well-conditioned). However, if the diffusion coefficients are not sufficiently small, the resulting system becomes ill-conditioned and hence the CGNR method converges very slowly; see for instance, [23] and numerical results in Section 4. Pang and Sun developed an efficient multigrid method by following the idea in [8] and using the Toeplitz-like structure of the coefficient matrices, which only requires O(NlogN) operations at each time step [23]. Fast algorithms for two-dimensional space fractional diffusion equations also have been studied; see [34], [35] and references therein.

In this paper, we use the well-known Crank–Nicholson method to discretize the FDE and then use the spatial extrapolation method to obtain temporally and spatially second-order accurate numerical estimates. We propose preconditioned conjugate gradient normal residual (preconditioned CGNR) method and preconditioned generalized minimal residual (preconditioned GMRES) method for the discretized systems, with preconditioners being constructed based on banded matrices of fixed bandwidth. As the multigrid method [23], the cost per iteration of the above preconditioned iterative methods is O(NlogN) operations.

The paper is organized as follows. In Section 2, we discretize the FDE (1) by the Crank–Nicholson method and give the structure of the coefficient matrix. In Section 3, we propose the preconditioned CGNR method and the preconditioned GMRES for the discretized system of (1). In Section 4 we present numerical results to illustrate the efficiency of the proposed algorithms.

Section snippets

The Crank–Nicholson finite difference scheme

Let N and M be positive integers, and let h=(xRxL)/N and Δt=T/M be the sizes of spatial grid and time step, respectively. We define a spatial and temporal partition:xn=xL+nh,n=0,1,,N;tm=mΔt,m=0,1,,M. Letun(m)=u(xn,tm),d±,n(m)=d±(xn,tm),fn(m+1/2)=f(xn,(tm+tm+1)/2). The Crank–Nicholson finite difference scheme (CN scheme) is given byun(m+1)un(m)Δt12(d+,n(m+1)hαj=0n+1gj(α)unj+1(m+1)+d,n(m+1)hαj=0Nn+1gj(α)un+j1(m+1))12(d+,n(m)hαj=0n+1gj(α)unj+1(m)+d,n(m)hαj=0Nn+1gj(α)un+j1(m))=fn(m

Preconditioned CGNR method and preconditioned GMRES method

Since the coefficient matrix IA(m+1) in (5) is not symmetric positive definite, we consider solving it by using the preconditioned CGNR method or the preconditioned GMRES method.

We require the following proposition [20], [21], [36] about the properties of gj(α).

Proposition 1

Let 1<α<2 and gj(α) be defined in (2). Then we have{g0(α)=1,g1(α)=α<0,g2(α)>g3(α)>>0,j=0gj(α)=0,j=0ngj(α)<0,for n1,gj(α)=O(j(α+1)).

We see from Proposition 1 that|g1(α)|>j=0,j1ngj(α), which implies that IA(m+1) is strongly

Numerical results

In this section, we employ the preconditioned CGNR method and the preconditioned GMRES method proposed in the previous section to solve the linear system (5). For all methods, we choose the initial guess asv0={u(0)[u0(x1),,u0(xN1)]T,m=0,u˜(m+1)2u(m)u(m1),m>0. We set u˜(m+1) as an initial guess since it is a good approximation of u(m+1) [36]. The stopping criterion isrj2b(m+1)2<107, where rj is the residual vector after j iterations. In our tests, we set the half bandwidth of the

Concluding remarks

We present two preconditioned iterative methods for the solution of initial–boundary value problems of anomalous diffusion of order α(1,2) in one space-dimension. We propose banded matrices as preconditioners for the discretization linear systems by exploiting the off-diagonal decay property of the coefficient matrices. Numerical results show that the proposed methods are very efficient.

We would like to make comments on related recent works and possible extension of our methods.

  • 1.

    In [34], a

Acknowledgements

We thank the anonymous referees for providing valuable comments and suggestions.

References (39)

Cited by (110)

  • An efficient multigrid method with preconditioned smoother for two-dimensional anisotropic space-fractional diffusion equations

    2022, Computers and Mathematics with Applications
    Citation Excerpt :

    The discretization on uniform mesh leads to a Toeplitz-like coefficient matrix [15] which brings out the design of iterative solvers for FDEs. For one-dimensional (1D) FDEs, we list some numerical methods such as the multigrid method [16], band preconditioning [17], circulant preconditioning [18], and the tridiagonal preconditioning [19]. In terms of the 2D case, we outline the splitting preconditioner presented in [20], the multigrid method employing banded-splitting smoother in [21], and the extrapolation cascadic multigrid (EXCMG) method investigated in [22].

  • On τ matrix-based approximate inverse preconditioning technique for diagonal-plus-Toeplitz linear systems from spatial fractional diffusion equations

    2022, Journal of Computational and Applied Mathematics
    Citation Excerpt :

    People use the FDEs to provide an adequate and accurate description for these anomalous diffusion, which include modeling chaotic dynamic of classical conservation systems [5], groundwater contaminant transport [6], turbulent flow [7], applications in biology [8], finance [9], image processing [10], physics [11] and flow in human meniscus [12–14]. As closed-form analytical solutions of fractional differential equations are usually unavailable especially in the presence of variable coefficients, hence, a lot of useful numerical approximations have been developed for these fractional derivatives, see, for instance, [15–20]. Owing to the Toeplitz-like structure of the resulting discretized systems, one may consider the circulant preconditioners for solving such systems.

View all citing articles on Scopus

The research was supported by NSF of China No. 11271238 and the Guangdong Provincial NSF under contract No. 10151503101000023 and the research grant UL020/08-Y5/MAT/JXQ01/FST from University of Macau.

View full text