Fast matrix splitting preconditioners for higher dimensional spatial fractional diffusion equations

Journal Pre-proof Fast matrix splitting preconditioners for higher dimensional spatial fractional diffusion equations

Zhong-Zhi Bai, Kang-Ya Lu

PII:

S0021-9991(19)30822-8

DOI:

https://doi.org/10.1016/j.jcp.2019.109117

Reference:

YJCPH 109117

To appear in:

Journal of Computational Physics

Received date:

30 March 2018

Revised date:

26 April 2019

Accepted date:

10 November 2019

Please cite this article as: Z.-Z. Bai, K.-Y. Lu, Fast matrix splitting preconditioners for higher dimensional spatial fractional diffusion equations, J. Comput. Phys. (2019), 109117, doi: https://doi.org/10.1016/j.jcp.2019.109117.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • Structured fast preconditioning methods for higher dimensional spatial fractional diffusion equations. • Theoretical analyses about the properties of the preconditioned methods. • Numerical advantages of the preconditioned Krylov subspace iteration methods over the state-of-the-art methods such as multigrid.

Fast Matrix Splitting Preconditioners for Higher Dimensional Spatial Fractional Diffusion Equations Zhong-Zhi Bai†

‡ §

and Kang-Ya Lu‡

∗

§

‡

State Key Laboratory of Scientific/Engineering Computing Institute of Computational Mathematics and Scientific/Engineering Computing Academy of Mathematics and Systems Science Chinese Academy of Sciences, P.O. Box 2719, Beijing 100190, P.R. China Emails: [email protected], [email protected] §

School of Mathematical Sciences University of Chinese Academy of Sciences Beijing 100049, P.R. China

Abstract The discretizations of two- and three-dimensional spatial fractional diffusion equations with the shifted finite-difference formulas of the Gr¨ unwald-Letnikov type can result in discrete linear systems whose coefficient matrices are of the form D + T , where D is a nonnegative diagonal matrix and T is a block-Toeplitz with Toeplitz-block matrix or a block-Toeplitz with each block being block-Toeplitz with Toeplitz-block matrix. For these discrete spatial fractional diffusion matrices, we construct diagonal and block-circulant with circulant-block splitting preconditioner for the two-dimensional case, and diagonal and block-circulant with each block being block-circulant with circulant-block splitting preconditioner for the threedimensional case, to further accelerate the convergence rates of Krylov subspace iteration methods, and we analyze the eigenvalue distributions for the corresponding preconditioned matrices. Theoretical results show that except for a small number of outliners the eigenvalues of the preconditioned matrices are located within a complex disk centered at 1 with the radius being exactly less than 1, and numerical experiments demonstrate that these structured preconditioners can significantly improve the convergence behavior of the Krylov subspace iteration methods. Moreover, this approach is superior to the geometric multigrid method and the preconditioned conjugate gradient methods incorporated with the approximate inverse circulant-plus-diagonal preconditioners in both iteration counts and computing times. Keywords: spatial fractional diffusion equations, shifted finite-difference discretization, block Toeplitz-like matrix, block circulant-like matrix, preconditioning, eigenvalue distribution. AMS(MOS) Subject Classifications: 65F08, 65F10, 65M22, 65N06, 65N22, 65Z05; CR: G1.3. ∗ †

Supported by The National Natural Science Foundation (No. 11671393, and No. 11911530082), P.R. China. The Corresponding Author

1

2

1

Z.-Z. Bai and K-Y. Lu

Introduction

Because fractional derivatives can precisely describe the memory and hereditary properties of various materials and processes, they have been extensively applied in many different areas such as anomalous diffusion [24, 30], groundwater contaminant transport [7, 8, 12], turbulent flow [9, 28], and image processing [13, 16]. In particular, the spatial fractional diffusion equations can adequately and accurately characterize the anomalous diffusion process, in which the overall motions of the particles are strongly interdependent and the mean-square displacement usually has the fractional power law dependence [6, 19]. In this paper, we consider two- and threedimensional spatial fractional diffusion equations of the form ⎧  r  βs  ⎪ ∂u(x,t) ∂ u(x,t) ∂ βs u(x,t) ⎪ d(x, t) = f (x, t), − + ⎪ β β s s ∂t ⎨ s=1

∂ + xs

∂ − xs

u(x, t)|xs =0 = u(x, t)|xs =1 = 0, s ∈ {1, 2, . . . , r}, ⎪ ⎪ ⎪ ⎩u(x, 0) = u (x), (x, t) ∈ (0, 1)r × (0, 1), 0

(1.1)

where x = (x1 , x2 , . . . , xr )∗ ∈ Rr is the vectorized r-dimensional spatial variable, t ∈ R is the βs temporal variable, d(x, t) is a given nonnegative function, f (x, t) is the source term, ∂ u(x,t) βs and

∂ βs u(x,t) s ∂ − xβ s

∂ + xs

are the left and the right Riemann-Liouville fractional derivatives with respect to

the variable xs , respectively, with βs ∈ (1, 2) being the orders of the derivatives and r being the dimension of the problem. Here (·)∗ indicates the transpose of a vector, and R represents the real number field. For more details, we refer to [23, 27]. Discretizing the spatial fractional diffusion equation in (1.1) by making use of the implicit finite-difference scheme, we can obtain a discrete system of linear equations of the coefficient matrix D + T , where D is a nonnegative diagonal matrix, and T is a block-Toeplitz with Toeplitzblock (BTTB) matrix for the two-dimensional (2D) case (i.e., r = 2), and is a block-Toeplitz with BTTB (BTBTTB) matrix for the three-dimensional (3D) case (i.e., r = 3). In the light of properties of Toeplitz and BTTB matrices [10], a matrix-vector multiplication with respect to the matrix D + T in either 2D or 3D case can be accomplished in O(n1 n2 ) or O(n1 n2 n3 ) computer storage, and in O(n1 n2 log(n1 n2 )) or O(n1 n2 n3 log(n1 n2 n3 )) computational complexity, respectively. Here, ns is the number of inner grid points in the xs direction, for s = 1, 2, . . . , r. In particular, for the one-dimensional (1D) spatial fractional diffusion equation in (1.1), in order to solve the corresponding discrete linear system effectively, in [20] Ng and Pan constructed an approximate inverse circulant-plus-diagonal (AICD) preconditioner to accelerate the convergence speed of the corresponding preconditioned conjugate gradient (CG) method, and showed by numerical experiments that the AICD-preconditioned CG (AICD-CG) method can be also effective for solving the 2D equations. However, the theoretical analysis for the 1D case and the numerical results for both 1D and 2D cases have verified that this approach is feasible and efficient only if the elements of the diagonal matrix D differ from each other mildly and smoothly. We refer to [14, 21] and the references therein for more and recent results related to this topic. In order to overcome this limitation, also for the 1D spatial fractional diffusion equation, in [5] Bai, Lu and Pan constructed a diagonal and circulant splitting (DCS) preconditioner to further accelerate the convergence speeds of the corresponding preconditioned GMRES and BiCGSTAB methods. The DCS preconditioning matrix is built upon the diagonal and Toeplitz splitting (DTS) iteration method, with the involved Toeplitz matrix being replaced by certain circulant

3

Fast Matrix Splitting Preconditioners

matrix; see also [3, 4]. It was proved in [5] that for the 1D equation the eigenvalues of the DCS-preconditioned matrix are clustered around 1, especially when the discretization stepsize is small; and the convergence property of the DCS-preconditioned GMRES and BiCGSTAB (in short, DCS-GMRES and DCS-BiCGSTAB) methods is much better than that of the AICD-CG method in iteration counts and computing times, and they also show a convergence behavior independent of the meshsize even for the one-dimensional spatial fractional diffusion equations of discontinuous or big-jump coefficients. In this paper, we are going to extend the DCS preconditioner for the 1D equation and develop the corresponding convergence theory to 2D and 3D spatial fractional diffusion equations in (1.1). To this end, we construct block diagonal and circulant-based splitting (BDCS) preconditioners or, in specific, diagonal and block-circulant with circulant-block (BCCB) splitting preconditioner for the 2D case, and diagonal and block-circulant with each block being BCCB (BCBCCB) splitting preconditioner for the 3D case, to further accelerate the convergence rates of the Krylov subspace iteration methods, and we analyze the eigenvalue distributions for the preconditioned matrices corresponding to the BDCS preconditioners, that is, the diagonal and BCCB splitting (DBCCBS) preconditioner and the diagonal and BCBCCB splitting (DBCBCCBS) preconditioner. Theoretical results show that except for a small number of outliners the eigenvalues of the preconditioned matrices are located within a complex disk centered at 1 with the radius being exactly less than 1, and numerical experiments demonstrate that these structured BDCS preconditioners can significantly improve the convergence behavior of the Krylov subspace iteration methods. Moreover, in terms of the computing times, these new approaches are superior to the geometric multigrid (GMG) method; and in both iteration counts and computing times, it outperforms the 2D and 3D extensions of the AICD-CG method. The remainder of the paper is organized as follows. In Section 2, we describe the discrete linear systems for both 2D and 3D spatial fractional diffusion equations in (1.1). In Section 3, we construct and analyze the BDCS preconditioners. The numerical results are reported and discussed in Section 4. Finally, in Section 5, we end this paper with a few conclusions and remarks.

2

The Discrete Linear Systems

Let (·)∗ and · denote the transpose and the Euclidean norm in the real space of either a vector or a matrix, λ(·) indicate any eigenvalue of a square matrix, ⊗ represent the Kronecker product symbol, · and · stand for the floor and the ceiling of the corresponding real number, and | · | represent the absolute value of the corresponding real or complex number. In addition, we write the identity matrix of suitable dimension as I, and the rank of a matrix as rank(·). Define a spatial partition xi1 = ih1 (i = 0, 1, . . . , n1 + 1), xj2 = jh2 (j = 0, 1, . . . , n2 + 1) and xk3 = kh3 (k = 0, 1, . . . , n3 + 1), and a temporal partition t =  Δt ( = 0, 1, . . . , m). In addition, we denote by dij = d(xi1 , xj2 , t ) for the 2D equation and dijk = d(xi1 , xj2 , xk3 , t ) for the 3D equation. By making use of the backward Euler formula for the temporal derivative and the shifted finite-difference formulas of Gr¨ unwald-Letnikov type [17, 18] for the spatial fractional derivatives, corresponding to an implicit finite-difference discretization of the spatial fractional diffusion equation in (1.1) we can obtain the discrete linear system (D + T )u = Δt f  + D u−1 ,

 = 1, 2, . . . , m,

(2.1)

4

Z.-Z. Bai and K-Y. Lu

at the -th temporal level, where the elements in the unknown vector u is arrayed in the natural lexicographic order; and for the 2D or the 3D case, D = diag(d11 , . . . , d1n2 , . . . , dn1 1 , . . . , dn1 n2 ) or D = diag(d111 , . . . , d11n3 , . . . , d1n2 1 , . . . , d1n2 n3 , . . . , dn1 n2 1 , . . . , dn1 n2 n3 ) is a nonnegative diagonal matrix, T =

Δt hβ1 1

Δt

(Tβ1 + Tβ∗1 ) ⊗ I +

hβ2 2

I ⊗ (Tβ2 + Tβ∗2 )

(2.2)

or T =

Δt hβ1 1

(Tβ1 + Tβ∗1 ) ⊗ I ⊗ I +

Δt hβ2 2

is a BTTB or BTBTTB matrix, with ⎛ (βs ) g1 ⎜ (βs ) ⎜ g ⎜ 2 ⎜ .. ⎜ . T βs = − ⎜ ⎜ .. ⎜ . ⎜ ⎜ (βs ) ⎝ gns −1 (β ) gnss

I ⊗ (Tβ2 + Tβ∗2 ) ⊗ I +

(βs )

g1

g0

g2 s .. . .. .

g1 s .. . .. .

··· .. . .. . .. . .. .

(β )

···

···

g0

0

(βs )

(βs )

(β )

(β )

gnss−1

Δt hβ3 3

I ⊗ I ⊗ (Tβ3 + Tβ∗3 )

0 .. . .. . .. .

0 0 .. .

(β ) g1 s (β ) g2 s

(β ) g0 s (β ) g1 s

0

(2.3)

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

being Toeplitz matrices of lower Hessenberg form. We refer to [5] and the references therein for more details about the matrix Tβs and its properties. From [5] we know that both Tβs and Tβ∗s (s = 1, 2, . . . , r) are strictly diagonally dominant M -matrices, and their eigenvalues can be bounded as æo,s hβs s < |λ(Tβs )| < 2βs − æo,s hβs s

and

æo,s hβs s < |λ(Tβ∗s )| < 2βs − æo,s hβs s ,

respectively, where, for s = 1, 2, . . . , r, æo,s =

eβs +1 (4 − βs )9/2 s Γ(−βs )

(2.4)

e13/12 59/2 2βs β

are positive constants dependent on βs . Hence, T is a symmetric positive definite M -matrix for either the 2D or the 3D equation, and the eigenvalues of T satisfy 2

r  s=1

æo,s Δt < λ(T ) < 4

r  s=1

s βs h−β s

Δt − 2

r 

æo,s Δt,

(2.5)

s=1

which  r indicates  that the condition number of the matrix T in the Euclidean norm is of the order  −βs O . hs s=1

5

Fast Matrix Splitting Preconditioners

3

The Block DCS Preconditioning

In this section, based on certain circulant approximations of Toeplitz matrices, by making use of the Kronecker product structures of the matrix T in (2.2) and (2.3) we construct structured fast preconditioners for the sequence of discrete linear systems in (2.1) and analyze eigenvalue properties of the corresponding preconditioned matrices. For simplicity, we drop the superscripts of vectors and matrices in the discrete linear systems in (2.1) associated with the time, and denote by A=D+T

and

b = Δt f + Du.

The DTS preconditioner proposed in [5] for the discretized 1D equation, when straightforwardly extended to the higher dimensional cases, is of the form MC (α) =

1 (αI + D)(αI + C), 2α

(3.1)

where α is a given positive constant, C is a BCCB matrix defined as C=

Δt hβ1 1

(Cβ1 + Cβ∗1 ) ⊗ I +

Δt hβ2 2

I ⊗ (Cβ2 + Cβ∗2 )

(3.2)

corresponding to the BTTB matrix T in (2.2) for the 2D case, and C is a BCBCCB matrix defined as C=

Δt hβ1 1

(Cβ1 + Cβ∗1 ) ⊗ I ⊗ I +

Δt hβ2 2

I ⊗ (Cβ2 + Cβ∗2 ) ⊗ I +

Δt hβ3 3

I ⊗ I ⊗ (Cβ3 + Cβ∗3 )

(3.3)

corresponding to the BTBTTB matrix T in (2.3) for the 3D case, with Cβs (s = 1, 2, . . . , r) being the Strang’s circulant approximations of Tβs ; see [29, 11]. In brief, we call the preconditioners MC (α) in (3.1) specified with both matrices in (3.2) and (3.3) as the BDCS preconditioners. We claim that for either the 2D or the 3D case the matrix C is a good approximation to the matrix T , as the difference between them is only a sum of a small-norm matrix with a low-rank matrix. This fact is precisely demonstrated in the following lemma. Lemma 3.1 For the spatial fractional diffusion equation in (1.1), assume that max hs ≤ 1≤s≤r

1 6.

For s = 1, 2, . . . , r, let the positive constants æo,s be defined in (2.4), and æs =

eβs +1 (4 − βs )9/2 , e13/12 59/2 βs Γ(−βs )

θs =

3βs +1 . βs (3 − βs )βs +1 Γ(−βs )

Then, for the matrices T and C defined in (2.2) and (3.2) for r = 2, or in (2.3) and (3.3) for r = 3, if a positive constant  and positive integer ko satisfy         1 θs Δt 1/βs βs −βs − 1, (3.4) max 2 θs Δt ≤  < min θs hs Δt and ko = max 1≤s≤r 1≤s≤r 1≤s≤r hs  there exist two matrices E ∈ Rn1 n2 ×n1 n2 and F ∈ Rn1 n2 ×n1 n2 for r = 2, or E ∈ Rn1 n2 n3 ×n1 n2 n3 and F ∈ Rn1 n2 n3 ×n1 n2 n3 for r = 3, such that T − C = E + F,

6

Z.-Z. Bai and K-Y. Lu

with E and F satisfying rank(E) ≤ 2ko

r 

ns

s=1

and E ≤

r  

 s Δt, βs − æo,s hβs s h−β s

F  < r.

s=1

Proof: We only demonstrate the case for r = 2, as the case for r = 3 can be verified analogously. Denote by Ts =

Δt hβs s

(Tβs + Tβ∗s )

and

Cs =

Δt hβs s

(Cβs + Cβ∗s ).

Then, according to the proof of Lemma 4.3 in [5], for given  and ko satisfying (3.4), there exist symmetric matrices Es ∈ Rn1 n2 ×n1 n2 and Fs ∈ Rn1 n2 ×n1 n2 , satisfying rank(Es ) ≤ 2ko and   s Δt, Fs  < , (3.5) Es  ≤ βs − æo,s hβs s h−β s such that Ts − Cs = Es + Fs . Let E = E1 ⊗ I + I ⊗ E2

and

F = F1 ⊗ I + I ⊗ F 2 .

Then we have rank(E) ≤ 2ko (n1 + n2 ) and T − C = (T1 − C1 ) ⊗ I + I ⊗ (T2 − C2 ) = (E1 + F1 ) ⊗ I + I ⊗ (E2 + F2 ) = E + F. By making use of (3.5), we can straightforwardly obtain the bounds for the matrices E and F . 2 Based on this lemma, we can prove the following theorem, which describes clustering property for the eigenvalues of the BDCS-preconditioned matrix MC (α)−1 A. Theorem 3.1 For the spatial fractional diffusion equation in (1.1), represent by dmin , tmin and dmax , tmax lower and upper bounds, respectively, for the eigenvalues of the matrices D and T . Assume that max hs ≤ 16 . Denote by 1≤s≤r

  r r   s βs h−β Δt − 2 æ Δt 2 α+4 s o,s s=1 s=1   , ηr =  r r   α+2 æo,s Δt α+2 æs Δt s=1

s=1

7

Fast Matrix Splitting Preconditioners

where, for any s ∈ {1, 2, . . . , r}, æo,s , æs and θs are positive constants defined in (2.4) and Lemma 3.1. Then, for a given positive constant  and a prescribed positive integer ko satisfying (3.4), there exist two matrices E(α) ∈ Rn1 n2 ×n1 n2 and F (α) ∈ Rn1 n2 ×n1 n2 for r = 2, or E(α) ∈ r  ns and Rn1 n2 n3 ×n1 n2 n3 and F (α) ∈ Rn1 n2 n3 ×n1 n2 n3 for r = 3, satisfying rank(E(α)) ≤ 2ko s=1

E(α) ≤

r 

s (βs − æo,s hβs s )h−β s ηr Δt,

F (α) < rηr ,

(3.6)

s=1

such that MC (α)−1 A = M (α)−1 A + E(α) + F (α), where M (α) =

1 (αI + D)(αI + T ), 2α

(3.7)

and MC (α) is defined in (3.1). Moreover, all eigenvalues of the matrix M (α)−1 A are located within a complex disk centered at 1 with the radius being given by =

|α − d| |α − λ| max < 1. α + d tmin ≤λ≤tmax α + λ

max

dmin ≤d≤dmax

(3.8)

Proof: In a similar fashion, we only demonstrate the case for r = 2. In accordance with Lemma 3.1 and the definitions of the matrices M (α) and MC (α) in (3.7) and (3.1), we have MC (α)−1 A = MC (α)−1 M (α) M (α)−1 A = (αI + C)−1 (αI + T )M (α)−1 A   = I + (αI + C)−1 (T − C) M (α)−1 A = M (α)−1 A + (αI + C)−1 E M (α)−1 A + (αI + C)−1 F M (α)−1 A. From [5] we know that the eigenvalues of the matrix C can be bounded as 2

2 

æs Δt < λ(C) < 4

s=1

2 

s βs h−β s

Δt − 2

s=1

2  s=1

Hence, it holds that (αI + C)−1  ≤

1 . α + 2(æ1 + æ2 ) Δt

By introducing the matrix N (α) =

1 (αI − D)(αI − T ), 2α

we see that A = M (α) − N (α),

æs Δt.

8

Z.-Z. Bai and K-Y. Lu

so that M (α)−1 A = I − M (α)−1 N (α)   = (αI + T )−1 I − (αI + D)−1 (αI − D)(αI − T )(αI + T )−1 (αI + T ). By making use of (2.5) and (3.8), we obtain M (α)−1 A ≤ 2(αI + T )−1 αI + T    −β2 1 + β h ) Δt − 2(æ + æ ) Δt 2 α + 4(β1 h−β 2 o,1 o,2 1 2 ≤ . α + 2(æo,1 + æo,2 ) Δt Thereby, it straightforwardly follows that (αI + C)−1 M (α)−1 A ≤ η2 . Now, with the notation E(α) = (αI + C)−1 E M (α)−1 A

and

F (α) = (αI + C)−1 F M (α)−1 A,

we know that rank(E(α)) = rank(E) ≤ 2ko (n1 + n2 ). Furthermore, in the light of Lemma 3.1, we can derive the bounds for E(α) and F (α) in (3.6) immediately. 2 According to Theorem 3.1 and (3.8), we see that with respect to both 2D and 3D cases, except for a number of outliners, the eigenvalues of BDCS-preconditioned matrix MC (α)−1 A are located within a complex disk centered at 1 with the radius being less than , especially when s Δt is sufficiently small. This implies that the BDCS preconditioners can significantly max h−β s 1≤s≤r

accelerate the convergence speeds of the Krylov subspace iteration methods [1, 2, 25].

4

Numerical Results

In this section, we use examples about the 2D and 3D spatial fractional diffusion equations to show the feasibility and effectiveness of the BDCS-preconditioned Krylov subspace iteration methods such as GMRES and BiCGSTAB (in short, BDCS-GMRES and BDCS-BiCGSTAB; see, e.g., [25, 26]) in terms of the number of iteration steps (denoted as “IT”) and the computing time in seconds (denoted as “CPU”), when these iteration methods are employed to solve the discrete linear systems in (2.1) at the initial temporal level with  = 0, that is, the linear systems of the form Au = b,

with

A = D1 + T,

u = u1

and

b = Δt f 1 + D1 u0 .

(4.1)

The reason that we adopt GMRES and BiCGSTAB as linear iteration solvers is that the BDCS preconditioners are nonsymmetric and, hence, are only applicable to precondition the Krylov subspace methods designed for solving nonsymmetric linear systems. In addition, as comparisons, we also solve the linear systems in (4.1) by making use of the GMG method implemented

Fast Matrix Splitting Preconditioners

9

in [15, 22], and the AICD-CG methods extended from the AICD-CG method in [20] for solving the discrete linear systems with respect to the 1D spatial fractional diffusion equations, in which we take two interpolating points, correspondingly resulting in the iteration methods named as AICD-CG(2). In our implementations, all iterations are started from the zero vector, terminated once the relative residual error at the current iterate, say, u(k) , satisfies b − Au(k)  ≤ 10−6 × b, or the maximal number 10, 000 of iteration steps is exceeded. The latter case is denoted as “—” in the numerical tables. In the computations, we take h := h1 = h2 = h3 = Δt with the notation n := n1 = n2 = n3 , and set the parameters α involved in the BDCS preconditioners to be the experimentally computed optimal values that minimize the total numbers of iteration steps of the corresponding preconditioned GMRES and BiCGSTAB methods. In addition, all experiments are carried out using MATLAB (version R2014a) on a personal computer with 3.60 GHz central processing unit (Intel(R) Core(TM) i7-4790 CPU @3.60GHz), 8.00 GB memory, and Windows 7 operating system.

4.1

The Results for 2D Equations

Example 4.1 We consider the 2D fractional diffusion equation (1.1), in which the coefficient function d(x, t) := d(x1 , x2 , t) and the source term f (x, t) := f (x1 , x2 , t) are, respectively, set to be d(x1 , x2 , t) = | sin(2πx1 )| cosh(8x2 + π) and f (x1 , x2 , t) = 2tx1 x2 (1 − x1 )(1 − x2 ) d(x1 , x2 , t) −(t2 + 0.01) [x2 (1 − x2 ) g(x1 , β1 ) + x1 (1 − x1 ) g(x2 , β2 )] ,

(4.2)

with g(ς, β) =

2ς 2−β + 2(1 − ς)2−β ς 1−β + (1 − ς)1−β − , Γ(2 − β) Γ(3 − β)

(4.3)

where Γ(·) is the Gamma function. In addition, the initial condition u0 (x) := u0 (x1 , x2 ) is taken as u0 (x1 , x2 ) = 0.01x1 x2 (1 − x1 )(1 − x2 ).

(4.4)

Then, this equation has the exact solution u(x, t) := u(x1 , x2 , t) given by u(x1 , x2 , t) = (t2 + 0.01)x1 x2 (1 − x1 )(1 − x2 ). We test three choices of the orders β1 and β2 of the spatial fractional derivatives: Case (a) (β1 , β2 ) = (1.1, 1.1), Case (b) (β1 , β2 ) = (1.2, 1.8), and Case (c) (β1 , β2 ) = (1.9, 1.9).

(4.5)

10

Z.-Z. Bai and K-Y. Lu

In Tables 2-4, we report the iteration counts and computing times of CG, AICD-CG(2), BDCS-GMRES, BDCS-BiCGSTAB and GMG with respect to the above three choices of the spatial fractional derivative orders β1 and β2 , among them in BDCS-GMRES and BDCSBiCGSTAB we adopt the parameters α in Table 1. From Tables 2-4, we see that both BDCS-GMRES and BDCS-BiCGSTAB methods significantly outperform the CG and AICD-CG(2) methods either in iteration steps or in CPU times. Besides, BDCS-GMRES is superior to GMG in CPU time for Cases (a) and (b), while BDCSBiCGSTAB is superior to GMG also in CPU time for all three cases. Moreover, the iteration steps and the computing times of BDCS-GMRES are about twice as much as those of BDCSBiCGSTAB as the grids are refined. Hence, the BDCS-BiCGSTAB method is the most effective for solving this discrete 2D spatial fractional diffusion equation. We also point out that AICDCG(2) requires about one third iteration steps of CG, but costs much more CPU time when the grid number n is increasing, which indicates that the computational execution of the AICD preconditioner is time-consuming; and as GMG takes fewer iteration steps but costs much more CPU times, the computing time in each iteration step of GMG is very large. In addition, the iteration steps of both CG and AICD-CG(2) are increasing linearly as n becomes large, but decreasing linearly for every fixed n as both β1 and β2 grow up. The first of these phenomena mainly results from the more and more ill-conditioning of the BTTB matrix T when the meshsize becomes smaller, while the second of which stems from the more and more well-conditioning of the linear system in (4.1) with either (2.2) or (2.3) caused by the displacement matrix D when β1 and β2 become larger; see Figure 1. As AICD-CG(2) takes much more iteration steps than BDCS-GMRES and BDCS-BiCGSTAB, we observe that the BDCS preconditioner can significantly improve the conditioning of the linear system (4.1) with either (2.2) or (2.3), through tightly clustering the eigenvalues of the BDCS-preconditioned matrices; see Figures 1 and 2. Moreover, the BDCS-GMRES and BDCS-BiCGSTAB methods show h-independent convergence behavior for Case (a), see Table 2; and linear convergence behavior for both Cases (b) and (c), see Tables 3 and 4. Because the iteration steps of BDCS-GMRES and BDCS-BiCGSTAB are increasing as β1 and β2 become large, it is clear that the BDCS preconditioner can successfully remove the influence of the displacement matrix D and greatly reduce the influence of the BTTB matrix T on the conditioning of the matrix A; see Figure 2. Case

Method

(a)

BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB

(b) (c)

127 3.0 6.0 50.0 48.0 100.0 84.0

255 4.0 4.0 50.0 45.0 100.0 72.0

n 511 2.0 3.7 40.0 47.0 100.0 108.0

1023 2.0 2.0 50.0 49.1 130.0 147.0

2047 2.0 1.9 47.0 51.8 120.0 145.8

Table 1: The experimentally computed optimal values of the parameter α for Example 4.1.

Example 4.2 We consider the 2D fractional diffusion equation (1.1), in which the coefficient

11

Fast Matrix Splitting Preconditioners

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 666 3.56 267 2.83 11 0.12 5.5 0.10 2 0.15

255 788 17.30 307 12.66 10 0.36 5.5 0.38 2 0.65

n 511 907 141.22 354 95.38 10 2.39 5.0 1.85 2 4.31

1023 1031 548.23 399 562.45 10 8.50 5.0 6.38 2 24.96

2047 1159 3476.70 442 4988.13 9 43.87 4.5 32.38 2 182.83

Table 2: The iteration steps and computing times for Case (a) of Example 4.1.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 451 2.47 185 1.97 21 0.24 13.0 0.21 2 0.32

255 555 12.15 222 9.06 26 1.06 16.5 0.99 3 2.07

n 511 675 105.72 266 70.22 32 9.06 18.0 6.48 3 13.01

1023 808 410.49 317 444.96 38 53.73 23.5 28.91 4 80.80

2047 961 2888.13 371 4172.62 44 324.28 25.5 178.48 5 564.89

Table 3: The iteration steps and computing times for Case (b) of Example 4.1.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 295 1.67 110 1.27 21 0.26 12.0 0.20 3 0.32

255 360 9.20 131 5.97 28 1.32 16.0 0.98 3 1.41

n 511 430 67.09 157 43.32 36 12.62 20.0 7.12 3 8.63

1023 518 275.48 183 256.16 44 69.96 28.0 34.35 3 44.23

2047 635 1918.63 211 2383.90 54 440.02 33.5 234.06 3 291.64

Table 4: The iteration steps and computing times for Case (c) of Example 4.1.

12

Z.-Z. Bai and K-Y. Lu

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

−0.2

−0.2

−0.2

−0.4

−0.4

−0.4

−0.6

−0.6

−0.6

−0.8

−0.8

−0.8

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

Figure 1: The eigenvalues of the matrix A for Example 4.1 when n = 127: Case (a) (left), Case (b) (middle), Case (c) (right).

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

−0.2

−0.2

−0.2

−0.4

−0.4

−0.4

−0.6

−0.6

−0.6

−0.8

−0.8

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

−0.8 0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

Figure 2: The eigenvalues of the BDCS-preconditioned matrix MC (α)−1 A for Example 4.1 when n = 127: Case (a) (left), Case (b) (middle), Case (c) (right).

13

Fast Matrix Splitting Preconditioners

function d(x, t) := d(x1 , x2 , t) is set to be ⎧ ⎨

1 2 2+ 2

x1 x2 , d(x1 , x2 , t) = x2 sin(πx1 ) e ⎩1.2 − sin(πx1 ) sin(πx2 ),

(x1 , x2 ) ∈ (0.5, 1] × (0.5, 1], otherwise,

and the source term f (x, t) := f (x1 , x2 , t) and the initial condition u0 (x) := u0 (x1 , x2 ) are set to be the same as in (4.2)-(4.3) and (4.4), respectively. Then, this equation has the exact solution u(x, t) := u(x1 , x2 , t) given by (4.5). Note that for this example the coefficient function d(x, t) is discontinuous with big jumps at line segments (x1 , x2 ) ∈ 0.5 × [0.5, 1] and (x1 , x2 ) ∈ [0.5, 1] × 0.5. We also test three choices of the orders β1 and β2 of the spatial fractional derivatives: Case (a) (β1 , β2 ) = (1.1, 1.1), Case (b) (β1 , β2 ) = (1.2, 1.8), and Case (c) (β1 , β2 ) = (1.9, 1.9). In Tables 6-8, we report the iteration counts and computing times of CG, AICD-CG(2), BDCS-GMRES, BDCS-BiCGSTAB and GMG with respect to the above three choices of the spatial fractional derivative orders β1 and β2 , among them in BDCS-GMRES and BDCSBiCGSTAB we adopt the parameters α in Table 5. From the results, we observe that even for the discontinuous coefficient case, the BDCSGMRES and BDCS-BiCGSTAB methods still significantly outperform the CG and AICDCG(2) methods either in iteration steps or in CPU times, and BDCS-BiCGSTAB is superior to GMG in CPU time again. Moreover, both BDCS-GMRES and BDCS-BiCGSTAB show hindependent convergence behavior for Case (a), see Table 6; and linear convergence behavior for both Cases (b) and (c), see Tables 7 and 8. In addition, BDCS-BiCGSTAB takes about half of the iteration steps of BDCS-GMRES, and consumes about half, or even one third, of the CPU times of BDCS-GMRES for Case (a), or for Cases (b) and (c). Case

Method

(a)

BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB

(b) (c)

127 1.7 1.8 19.0 34.0 50.0 56.0

255 2.0 1.9 23.0 24.9 50.0 58.0

n 511 1.7 1.4 25.0 40.0 57.0 54.7

1023 1.7 1.4 30.0 36.0 64.0 80.8

2047 1.7 1.8 43.0 40.0 73.0 77.0

Table 5: The experimentally computed optimal values of the parameter α for Example 4.2.

14

Z.-Z. Bai and K-Y. Lu

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 1305 7.37 389 4.41 16 0.18 9.5 0.16 2 0.15

255 1743 46.35 564 25.32 16 0.62 9.0 0.54 2 0.64

n 511 1023 2035 2202 316.65 1145.92 840 1298 231.77 1815.43 17 17 4.10 15.50 10.0 10.0 3.62 12.38 2 2 4.26 24.89

2047 2277 6905.74 1988 22385.03 17 93.09 9.0 62.52 2 174.09

Table 6: The iteration steps and computing times for Case (a) of Example 4.2.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 1017 5.71 199 2.26 36 0.43 21.5 0.33 3 0.40

255 1403 35.93 298 13.76 47 2.79 27.0 1.54 3 1.72

n 511 1023 1739 1991 271.69 1057.67 475 736 131.76 1027.00 61 76 27.26 158.69 36.5 44.0 12.93 53.56 3 4 10.38 72.11

2047 2166 6553.69 1089 12287.18 88 925.92 53.5 366.71 4 392.05

Table 7: The iteration steps and computing times for Case (b) of Example 4.2.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

127 719 4.03 90 1.01 37 0.48 24.0 0.36 5 0.85

255 1059 26.99 139 6.24 50 2.80 33.5 1.99 6 4.55

n 511 1400 217.73 197 54.26 68 32.32 45.0 15.90 7 31.43

1023 1628 865.74 280 391.17 90 223.36 58.5 71.09 9 187.96

2047 1880 5677.83 410 4635.95 117 1486.84 83.5 576.55 11 1230.73

Table 8: The iteration steps and computing times for Case (c) of Example 4.2.

15

Fast Matrix Splitting Preconditioners

4.2

The Results for 3D Equation

Example 4.3 We consider the 3D fractional diffusion equation (1.1), in which the coefficient function d(x, t) := d(x1 , x2 , x3 , t) and the source term f (x, t) := f (x1 , x2 , x3 , t) are, respectively, set to be 1

d(x1 , x2 , x3 , t) = e x1 +0.1 | sin(2πx2 )| cosh(5x3 ) and f (x1 , x2 , x3 , t) = 2tx1 x2 x3 (1 − x1 )(1 − x2 )(1 − x3 ) d(x1 , x2 , x3 , t) − (t2 + 0.01) × [x2 x3 (1 − x2 )(1 − x3 ) g(x1 , β1 ) +x1 x3 (1 − x1 )(1 − x3 ) g(x2 , β2 ) + x1 x2 (1 − x1 )(1 − x2 ) g(x3 , β3 )] , with the function g(ς, β) being defined in (4.3). In addition, the initial condition u0 (x) := u0 (x1 , x2 , x3 ) is taken as u0 (x1 , x2 , x3 ) = 0.01x1 x2 x3 (1 − x1 )(1 − x2 )(1 − x3 ). Then, this equation has the exact solution u(x, t) := u(x1 , x2 , x3 , t) given by u(x1 , x2 , x3 , t) = (t2 + 0.01)x1 x2 x3 (1 − x1 )(1 − x2 )(1 − x3 ). We test three choices of the orders β1 and β2 of the spatial fractional derivatives: Case (a) (β1 , β2 , β3 ) = (1.1, 1.1, 1.1), Case (b) (β1 , β2 , β3 ) = (1.2, 1.5, 1.8), and Case (c) (β1 , β2 , β3 ) = (1.9, 1.9, 1.9). In Tables 10-12, we report the iteration counts and computing times of CG, AICD-CG(2), BDCS-GMRES, BDCS-BiCGSTAB and GMG with respect to the above three choices of the spatial fractional derivative orders β1 , β2 and β3 , among them in BDCS-GMRES and BDCSBiCGSTAB we adopt the parameters α in Table 9. We first see that the CG method fails to converge for Case (a) when n = 255, and the BDCSGMRES method fails to compute an approximate solution due to insufficient computer memory for Case (c) when n = 255. For the convergent situations, we observe from Tables 10-12 that the BDCS-BiCGSTAB method significantly outperforms the CG, AICD-CG(2) and BDCS-GMRES methods either in iteration steps or in CPU times, and it also dramatically outperforms the GMG method in CPU times. Hence, BDCS-BiCGSTAB is the most effective for solving this discrete 3D spatial fractional diffusion equation. In addition, these results show again that the AICD preconditioner can only slightly improve the conditioning of the linear systems in (4.1) with (2.3), while the BDCS preconditioner can successfully eliminate the influence of the displacement matrix D and greatly reduce the influence of the BTBTTB matrix T on the conditioning of the matrix A; see Figures 3 and 4. Besides, the computing times of GMG increase rapidly when the grids are refined, while its iteration steps keep constant. This implies that GMG should cost huge CPU time in each iteration step. Therefore, BDCS-BiCGSTAB is the most economic and high-effective method, although it takes more iteration steps than GMG.

16

Z.-Z. Bai and K-Y. Lu

Case

Method

(a)

BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB BDCS-GMRES BDCS-BiCGSTAB

(b) (c)

15 5.0 6.0 18.0 18.0 39.0 36.0

31 3.0 4.0 22.0 24.0 50.0 75.0

n 63 3.0 4.0 25.0 27.0 60.0 64.0

127 3.0 3.0 25.0 30.0 60.0 60.0

255 2.0 2.0 27.0 31.0 — 63.0

Table 9: The experimentally computed optimal values of the parameter α for Example 4.3.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

n 15 414 1.55 173 1.37 11 0.15 6.5 0.13 3 0.64

31 1449 35.20 484 29.58 13 0.89 7.5 0.70 3 0.70

63 3455 597.28 1188 394.47 14 6.31 7.5 4.75 3 8.72

127 6113 9779.05 2424 7340.75 14 56.17 8.0 42.49 3 195.19

255 — — 3857 105226.07 14 2473.32 8.0 393.44 3 11350.15

Table 10: The iteration steps and computing times for Case (a) of Example 4.3.

Method CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

Index IT CPU IT CPU IT CPU IT CPU IT CPU

n 15 279 1.06 113 0.99 13 0.18 7.0 0.13 3 0.41

31 924 21.90 315 18.63 16 1.06 9.0 0.96 3 0.91

63 2116 390.38 682 234.40 20 9.62 11.0 6.42 3 10.99

127 3662 5803.93 1160 3514.51 25 131.87 15.0 79.99 3 234.29

255 5269 90331.24 1767 46967.60 32 93675.17 18.5 903.70 3 12394.48

Table 11: The iteration steps and computing times for Case (b) of Example 4.3.

17

Fast Matrix Splitting Preconditioners

Method

Index

CG AICD-CG(2) BDCS-GMRES BDCS-BiCGSTAB GMG

15 207 0.89 67 0.56 15 0.16 9.0 0.19 4 0.34

IT CPU IT CPU IT CPU IT CPU IT CPU

31 673 15.63 160 9.19 20 1.29 12.0 1.08 4 0.99

n 63 1578 290.24 311 110.22 27 14.70 16.0 9.10 5 16.22

127 2764 4355.01 530 1591.68 37 236.06 20.5 107.63 4 262.96

255 3996 68101.13 817 21445.51 — — 28.0 1367.20 4 15616.59

Table 12: The iteration steps and computing times for Case (c) of Example 4.3.

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

−0.2

−0.2

−0.2

−0.4

−0.4

−0.4

−0.6

−0.6

−0.6

−0.8

−0.8

−0.8

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

Figure 3: The eigenvalues of the matrix A for Example 4.3 when n = 15: Case (a) (left), Case (b) (middle), Case (c) (right).

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

−0.2

−0.2

−0.2

−0.4

−0.4

−0.4

−0.6

−0.6

−0.6

−0.8

−0.8

−0.8

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

−1 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

Figure 4: The eigenvalues of the BDCS-preconditioned matrix MC (α)−1 A for Example 4.3 when n = 15: Case (a) (left), Case (b) (middle), Case (c) (right).

18

5

Z.-Z. Bai and K-Y. Lu

Concluding Remarks

The BDCS preconditioners are higher order extensions of the DCS preconditioner constructed in [5] for solving the discrete linear systems with respect to the 1D spatial fractional diffusion equations in (1.1). They can be effectively executed by the fast Fourier transform, and can tightly cluster the eigenvalues of the corresponding preconditioned matrices at around 1. Moreover, according to the numerical results, the BDCS-BiCGSTAB method outperforms the GMG method in terms of CPU time. However, the iteration steps of both BDCS-GMRES and BDCSBiCGSTAB may increase linearly for large values of the derivative orders. Hence, constructing faster and more robust circulant-based preconditioners for the discrete linear systems with respect to the 2D and 3D spatial fractional diffusion equations in (1.1) is still a challenge problem that deserves in-depth study. Moreover, in view of the almost constant iteration steps of the GMG method, constructing effective smoothers to further improve the computing efficiency of GMG is also a worthwhile topic that merits further discussion.

References [1] Z.-Z. Bai, Sharp error bounds of some Krylov subspace methods for non-Hermitian linear systems, Appl. Math. Comput., 109(2000), 273-285. [2] Z.-Z. Bai, Motivations and realizations of Krylov subspace methods for large sparse linear systems, J. Comput. Appl. Math., 283(2015), 71-78. [3] Z.-Z. Bai, G.H. Golub and M.K. Ng, Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, SIAM J. Matrix Anal. Appl., 24(2003), 603-626. [4] Z.-Z. Bai, G.H. Golub and M.K. Ng, On inexact Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, Linear Algebra Appl., 428(2008), 413-440. [5] Z.-Z. Bai, K.-Y. Lu and J.-Y. Pan, Diagonal and Toeplitz splitting iteration methods for diagonal-plus-Toeplitz linear systems from spatial fractional diffusion equations, Numer. Linear Algebra Appl., 24(2017), e2093, pp. 1-15. [6] E. Barkai, R. Metzler and J. Klafter, From continuous time random walks to the fractional Fokker-Planck equation, Phys. Rev. E, 61(2000), 132-138. [7] D.A. Benson, S.W. Wheatcraft and M.M. Meerschaert, Application of a fractional advection-dispersion equation, Water Resour. Res., 36(2000), 1403-1412. [8] D.A. Benson, S.W. Wheatcraft and M.M. Meerschaert, The fractional-order governing equation of l´evy motion, Water Resour. Res., 36(2000), 1413-1423. [9] B.A. Carreras, V.E. Lynch and G.M. Zaslavsky, Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence model, Phys. Plasmas, 8(2001), 5096-5103. [10] R.H. Chan and X.-Q. Jin, An Introduction to Iterative Toeplitz Solvers, SIAM, Philadelphia, PA, 2007.

Fast Matrix Splitting Preconditioners

19

[11] R.H. Chan and G. Strang, Toeplitz equations by conjugate gradients with circulant preconditioner, SIAM J. Sci. Statist. Comput., 10(1989), 104-119. [12] J.W. Kirchner, X.-H. Feng and C. Neal, Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 403(2000), 524-527. [13] W.S. Kuklinski, Utilization of fractal image models in medical image processing, Fractals, 2(1994), 363-369. [14] S.-L. Lei and H.-W. Sun, A circulant preconditioner for fractional diffusion equations, J. Comput. Phys., 242(2013), 715-725. [15] X.-L. Lin, M.K. Ng and H.-W. Sun, A multigrid method for linear systems arising from time-dependent two-dimensional space-fractional diffusion equations, J. Comput. Phys., 336(2017), 69-86. [16] B. Mathieu, P. Melchior, A. Oustaloup and Ch. Ceyral, Fractional differentiation for edge detection, Signal Processing, 83(2003), 2421-2432. [17] M.M. Meerschaert and C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, J. Comput. Appl. Math., 172(2004), 65-77. [18] M.M. Meerschaert and C. Tadjeran, Finite difference approximations for two-sided spacefractional partial differential equations, Appl. Numer. Math., 56(2006), 80-90. [19] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep., 339(2000), 1-77. [20] M.K. Ng and J.-Y. Pan, Approximate inverse circulant-plus-diagonal preconditioners for Toeplitz-plus-diagonal matrices, SIAM J. Sci. Comput., 32(2010), 1442-1464. [21] J.-Y. Pan, R.-H. Ke, M.K. Ng and H.-W. Sun, Preconditioning techniques for diagonaltimes-Toeplitz matrices in fractional diffusion equations, SIAM J. Sci. Comput., 36(2014), A2698-A2719. [22] H.-K. Pang and H.-W. Sun, Multigrid method for fractional diffusion equations, J. Comput. Phys., 231(2012), 693-703. [23] I. Podlubny, Fractional Differential Equations, Math. Sci. Engrg. 198, Academic Press, San Diego, 1999. [24] K. Razminia, A. Razminia and D. Baleanu, Investigation of the fractional diffusion equation based on generalized integral quadrature technique, Appl. Math. Model., 39(2015), 86-98. [25] Y. Saad, Iterative Methods for Sparse Linear Systems, Second Edition, SIAM, Philadelphia, PA, 2003. [26] Y. Saad and M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 7(1986), 856-869. [27] S.G. Samko, A.A. Kilbas and O.I. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach Science Publishers, Yverdon, 1993.

20

Z.-Z. Bai and K-Y. Lu

[28] M.F. Shlesinger, B.J. West and J. Klafter, L´evy dynamics of enhanced diffusion: application to turbulence, Phys. Rev. Lett., 58(1987), 1100-1103. [29] G. Strang, A proposal for Toeplitz matrix calculations, Stud. Appl. Math., 74(1986), 171-176. [30] G.M. Zaslavsky, Chaos, fractional kinetics, and anomalous transport, 371(2002), 461-580.

Phys. Rep.,

Declaration of interests ‫ ܈‬The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. տThe authors declare the following financial interests/personal relationships which may be considered as potential competing interests: