A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations✩ Xin-Hui Shao, Zhen-Duo Zhang, Hai-Long Shen

∗

Department of Mathematics, College of Sciences, Northeastern University, Shenyang 110819, PR China

article

info

a b s t r a c t The spatial fractional diffusion equation can be discretized by employing the implicit finite difference scheme with the shifted Grünwald formula and the given discretized linear systems have a diagonal-plus-Toeplitz structure. In this paper, based on trigonometric transformation splitting (TTS), we study efficient iterative method called GTTS has two parameters. As a focus, we give the simple and effective optimal forms of these two parameters respectively. We carry out some numerical experiments to illustrate the effectiveness and accuracy of this new algorithm. © 2019 Elsevier Ltd. All rights reserved.

Article history: Received 21 October 2018 Received in revised form 7 July 2019 Accepted 6 October 2019 Available online xxxx Keywords: Spatial fractional diffusion equation Toeplitz linear system Two parameters Trigonometric transformation splitting

1. Introduction We consider spatial fractional diffusion equations of the form

⎧ ∂ u(x, t) ∂ β u(x, t) ∂ β u(x, t) ⎪ ⎪ − − = f (x, t), d(x, t) ⎪ ⎪ ∂t ∂ + xβ ∂− xβ ⎪ ⎨ u(xL , t) = u(xR , t) = 0, ⎪ ⎪ ⎪u(x, 0) = u0 (x), ⎪ ⎪ ⎩ (x, t) ∈ (xL , xR ) × (0, +∞), where d(x, t) is a prescribed nonnegative function, ∂ β u(x,t) ∂− xβ

∂ β u(x,t) ∂ + xβ

(1)

is the left Riemann–Liouville fractional derivative [1,2] and

is the right one, with the order being β ∈ (1, 2), that is,

∂ β u(x, t) 1 d2 = ∂+ xβ Γ (2 − β ) dx2

∫

x

(x − ξ )1−β u(ξ , t)dζ , xL

∫ x ∂ β u(x, t) 1 d2 = (x − ξ )1−β u(ξ , t)dζ , ∂− xβ Γ (2 − β ) dx2 xL where Γ (·) is the Gamma function. In recent years, the fractional diffusion equation (FDE) which tends to be better for description of various character and processes that the traditional integer order one has broad applications and has gained ✩ Project supported by the Natural Science Foundation of Liaoning Province, PR China (No. 20170540323). ∗ Corresponding author. E-mail address: [email protected] (H.-L. Shen). https://doi.org/10.1016/j.camwa.2019.10.003 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

2

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

considerable popularity and emphasis. The spatial FDEs (1), which is a kind of FDEs, can provide an accurate description for anomalous diffusion process which cannot be modeled accurately by traditional one. Eq. (1) has a broad application in many areas involving biology [3,4], chemistry [5,6], physics [7,8], finance [9,10], and so on. Hence, considerable numerical solutions as the major methods for FDEs have been developed [11–16]. However, certain salient features of the fractional differential operator lead to unconditional instability [17]. Besides, most numerical methods for FDEs generate a full coefficient matrix, which require computational cost of O(n3 ) and storage of O(n2 ) where n is the number of grid points [18]. Because the closed-form expression of these equations is usually not available, the numerical solution is a common and effective approximate solution for Eq. (1). There are a lot of methods about solving the discretized FDEs recently [19–25]. Through a proper finite difference discretization, the equation results in a linear system with the coefficient matrix D + T ∈ Rn×n , where D is a diagonal matrix with nonnegative elements and T is a symmetric positive definite Toeplitz matrix. This kind of linear systems cannot be solved by common Toeplitz solvers due to the addition of the diagonal matrix. Therefore, although the pure Toeplitz linear system can be solved in O(n2 ) or even O(n log2 n) and O(n log n), see the previous studies [26], the above diagonal-plus-Toeplitz linear system can only be solved directly via the Cholesky factorization in O(n3 ) [27]. The Krylov subspace iteration methods, represented by conjugate gradient (CG), are a kind of common and effective solution methods for solving these linear systems [28]. The farther studies provide some high-quality preconditioning methods, such as the preconditioned conjugate gradient (PCG), the preconditioned Krylov subspace iteration methods [27]. They have better convergence rate. However, in the meantime, these methods are confronted with deterioration of convergence rate [29]. Recently, the trigonometric transformation splitting (TTS) method is proposed by Liu, Wu, Qin and Zhang [30] to solve symmetric Toeplitz systems. Compared with other two-step splitting iteration methods, such as the classic symmetric Gauss–Seidel method (SGS), the positive definite and skew-symmetric splitting iterative method (PSS) in Bai et al. (2005) and the CSCS iterative method only for Toeplitz matrix system in Ng (2003), TTS method has smaller spectral radius and converges faster. Splitting the Toeplitz matrix by TTS is the forepart motivation of this paper. However, in our studies, we split the matrix D and generate a new parameter, which is the most important feature of this paper. This paper is organized as follows. In Section 2, we obtain the discretized linear system by finite difference discretization of Eq. (1). In Section 3, we present the GTTS iteration method and establish its convergence theory. As the most important part, in Section 4, we consider the two optimal parameters that have a huge impact on convergence rate. In Section 5, we give the numerical results. At last, we end this paper with a few conclusions and remarks. 2. The finite difference discretization of FDE According to [29], we can rewrite Eqs. (1) as follows (Dℓ+1 + T )uℓ+1 = ∆tf ℓ+1 + Dℓ+1 uℓ , ℓ = 0, 1, . . . , m − 1,

(2)

where 1 ℓ+1 ℓ+1 Dℓ+1 = diag(dℓ+ 1 , d2 , . . . , dn ),

uℓ = [uℓ1 , uℓ2 , . . . , uℓn ]∗ , f ℓ = [f1ℓ , f2ℓ , . . . , fnℓ ]∗ , is a diagonal matrix at the (ℓ + 1)th temporal level, and T =

∆t hβ

(Tβ + Tβ∗ )

(3)

with

⎛

(β)

g1

⎜ ⎜g (β) ⎜ 2 ⎜. ⎜. ⎜. Tβ = − ⎜ ⎜.. ⎜. ⎜ ⎜ (β) ⎝gn−1 (β)

gn

(β)

g0

0

.

..

.

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

(β) gn−1

···

···

(β)

g1

(β)

g2

.. ..

.

(β)

g0

(β)

g1

..

.

⎞

0

0

g1

⎟ ⎟ ⎟ .. ⎟ ⎟ . ⎟ ⎟. ⎟ 0 ⎟ ⎟ (β) ⎟ g0 ⎠

g2

g1

.. .. ..

. . .

(β) (β)

0

(β)

The process above and the analysis of the eigenvalue bounds are clearly described in [29]. We can see that T has only positive eigenvalues, and it is also a real symmetric positive definite matrix. Actually, the matrix T in Eq. (2) can be further expressed as the sum of two matrices: T = TC + TS . Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

3

The process is elaborated below. For the specific theory, see [30]. We denote the first version of the discrete cosine transform (DTC) matrix of order n by CnI , the first version of discrete sine transform (DST) matrix of order n by SnI :

√ SnI m,k

[ ]

=

2 n+1

√ [

CnI +2 m,k

]

=

sin

2 n+1

π (m + 1)(k + 1) , m, k = 0, 1, . . . , n − 1, n+1

εk εm cos

π mk n+1

, m, k = 0, 1, . . . , n + 1,

where

⎧ ⎨1,

εj =

if j ̸ = 0 and j ̸ = n + 1,

1

⎩√

2

if j = 0 or j = n + 1.

Let

√

√

Gn+2 = diag(1/ 2, 1, . . . , 1, 1/ 2),

√

2

(Cˆ n )m,k =

cos

π mk , m, k = 0, 1, . . . , n + 1. n+1

n+1 Then the discrete cosine transform matrix is

√ CnI +2

=

2 n+1

Gn+2 (cos(

π mk

))nm+,k1=0 Gn+2 .

n+1

It can possess the block structure as follows

⎡ √ CnI +2 =

1

⎢2 ⎢ ⎢ 1 ⎢√ e n+1 ⎢ ⎢ 2 ⎣1 2

2

1

√ eT √2 n+1 ˆ

Cn

2 1

√ fT

1

⎤

2 1

⎥ ⎥ ⎥ √ f⎥ ⎥, 2 ⎥ 1 ⎦ 2

2

where e = (1, 1, . . . , 1) and f = (−1, 1, . . . , (−1)n )T are the vectors of dimension n. Let αn+2 = [a0 , a1 , . . . , an+1 ]T ∈ Rn+2 , and it used to represent the elements of the above real symmetric positive definite matrix T . We define λ = (λ0 , . . . , λn+1 )T with T

λ=

√

2(n + 1)Gn+2 CnI +2 Gn+2 αn+2 ,

and we can get the following conclusion. Theorem 2.1 ([31,32]). A real symmetric Toeplitz matrix An admits a splitting An = TC + TS , where TC =

1 2

(Cˆ n ΛCˆ n + R2 ) and TS =

with Λ = diag(λ1 , . . . , λn ) and R2 =

1 2

(Sˆn ΛSˆn + R2 ),

λ0 n+1

eeT +

λn+1 T ff . n+1

Obviously, if we have no further information, we can set an = an+1 = 0 or choose an and an+1 such that λ0 = λn+1 = 0 and An = 12 (Cˆ n ΛCˆ n + Sˆn ΛSˆn ). According to Theorem 2.1, the TTS method with the matrix An , as the coefficient matrix for real linear systems given by Liu and Wu [30] is introduced here. The following lemma shows the fact that TC and TS all are positive definite. Lemma 2.1 ([30]). If the generating function of the Toeplitz matrix A is a real positive even function, then the matrices TC and TS are positive definite for sufficient large n. In fact, the matrix A conforming to Lemma 2.1 is a real symmetric Toeplitz matrix. The TTS iteration. [30] Given an initial guess u(0) ∈ Rn . For k = 0, 1, 2, . . . until {u(k) } converges, compute

{

1

(α I + TC )u(k+ 2 ) = (α I − TS )u(k) + b 1

(α I + TS )u(k+1) = (α I − TC )u(k+ 2 ) + b, where α is a prescribed positive constant. Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

4

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

3. The GTTS iteration methods Based on the discussion in the previous section, we consider GTTS matrix splitting iteration methods for solving the following system Au = b with A ∈ Rn×n , b ∈ Rn , where A = D + T , D ∈ Rn×n is a diagonal matrix with nonnegative elements and T ∈ Rn×n is a real symmetric positive definite Toeplitz matrix. According to Theorem 2.1, the matrix A can be further expressed as 1 1 A = D + TC + TS = ( D + TC ) + ( D + TS ), 2 2 where TC and TS are defined as in Theorem 2.1. The generalized TTS (GTTS) scheme is based on the splitting above. We can construct the fixed-equations with a given positive constant α for Eq. (2)

⎧ ⎪ ⎨(α I + ( 1 D + TC ))u = (α I − ( 1 D + TS ))u + b 2

2

⎪ ⎩(α I + ( 1 D + T ))u = (α I − ( 1 D + T ))u + b. S C 2

2

Note that we split the matrix D equally. But we find that the different splitting forms of D have a great influence on the convergence speed from numerical experiments. So we can introduce a real parameter ω ∈ [0, 1] and we have A = D + TC + TS = (ωD + TC ) + [(1 − ω)D + TS ], and (α I + (ωD + TC ))u = (α I − (1 − ω)D − TS )u + b (α I + (1 − ω)D + TS )u = (α I − (ωD + TC ))u + b.

{

Based on the above discussion, we can propose the following GTTS method. The GTTS iteration. Given an initial guess u(0) ∈ Rn . For k = 0, 1, 2, . . . until {u(k) } converges, compute

{

1

(α I + (ωD + TC ))u(k+ 2 ) = (α I − (1 − ω)D − TS )u(k) + b 1

(α I + (1 − ω)D + TS )u(k+1) = (α I − (ωD + TC ))u(k+ 2 ) + b,

(4)

where α and ω ∈ [0, 1] are the prescribed real positive constants. As shown in [30], there are some algorithms which can fast be calculated the above matrices and vectors by DCT and DST. By straightforward computations, elimination of 1 u(k+ 2 ) from (4) yields Mu(k+1) = Nu(k) + b, where M=

1 2α

(α I + ωD + TC )(α I + (1 − ω)D + TS )

and N =

1 2α

(α I − ωD − TC )(α I − (1 − ω)D − TS ).

Let M1 = ωD + TC and M2 = (1 − ω)D + TS . Then we can obtain the iteration matrix L = M −1 N = (α I + M2 )−1 (α I + M1 )−1 (α I − M1 )(α I − M2 ). As the two-step splitting iteration method, the convergence of the GTTS iterative method is an easy consequence of a classical result , which shows that it is convergent if and only if the spectral radius ρ (L) < 1 [33]. The following convergence result indicates that the GTTS iteration method is the unconditional convergence. Theorem 3.1. Let A = D + TC + TS = (ωD + TC ) + [(1 − ω)D + TS ] = M1 + M2 . The matrix D ∈ Rn×n is a diagonal matrix with nonnegative elements, TC , TS are defined as in Theorem 2.1, and ω ∈ [0, 1] is a prescribed real positive constant. Then, for any positive constant α and any initial guess u(0) ∈ Rn , the alternating iteration (4) converges unconditionally to the unique solution of Ax = b. Moreover, then the spectral radius ρ (L) of the GTTS iteration matrix L is bounded by

σ =

{ max

θmin ≤θ ≤θmax

|α − θ| α+θ

}

{ max

ηmin ≤η≤ηmax

|α − η| α+η

}

< 1.

(5)

where θ is an arbitrary eigenvalue of M1 , η is an arbitrary eigenvalue of M2 . Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Proof. The iteration matrix is similar to Lˆ = (α I + M1 )−1 (α I − M1 )(α I − M2 )(α I + M2 )−1 . We denote that ∥·∥ is the Euclidean norm of a matrix. Clearly,

ˆ ρ (L) = ρ (L)   ≤ (α I + M1 )−1 (α I − M1 )(α I − M2 )(α I + M2 )−1 

   ≤ (α I + M1 )−1 (α I − M1 ) (α I − M2 )(α I + M2 )−1  } { } { |α − η| |α − θ| max = max θ ∈sp(M1 ) α{+ θ η∈}sp(M1 ) α + { η } |α − θ| |α − η| ≤ max max θmin ≤θ ≤θmax α + θ ηmin ≤η≤ηmax α + η = σ, Because D is a diagonal matrix with nonnegative elements, TC and TS are also positive definite matrices, according to Lemma 2.1, we have σ < 1. □ The main potential advantage of the GTTS scheme than the classical methods for solving FDEs is that D is split into two parts and assigned to M1 and M2 respectively. There are two parameters including classical α and ω, and the different splitting forms of D have a great influence on the convergence speed. So the GTTS iteration matrix is more flexible for solving FDEs. 4. The optimal α and ω As we mentioned, the GTTS has two parameters and how to choose them has a huge impact on convergence rate. First, we consider the optimal parameter α . It has been shown in [29] that the choice of α leads to the minimization of convergence bound. Lemma 4.1 ([29]). Let A = D + T ∈ Rn×n , where D ∈ Rn×n is a diagonal matrix with nonnegative elements, T ∈ Rn×n is a Hermitian positive definite Toeplitz matrix. Denote that dmax and dmin the largest and the smallest elements of D, tmax and tmin the largest and the smallest elements of T . For the two-step implicit iteration of the splitting A = D + T , the spectral radius of the iteration matrix is bounded by

σ (α ) =

{ max dmin ≤d≤dmax

{

}

|α − d| α+d

max

tmin ≤λ≤tmax

|α − λ| α+λ

}

< 1,

where α is prescribed positive constant. And, with d∗ =

√

√ dmin dmax and t∗ =

tmin tmax ,

the optimal parameter α∗ that minimizes the function σ (α ) is determined in the following cases: (a) If d∗ = t∗ , then it holds that α∗ = d∗ = t∗ and

√ √ √ √ dmax − dmin tmax − tmin σ (α ) = √ . √ √ √ dmax + dmin tmax + tmin

(b) If d∗ < t∗ , then (b1) for dmin ≥ tmin , or for tmax > dmax √, tmin > dmin +tmax dmax tmax and ddmax + ≥ t d t min

min

min min

it holds that α∗ = d∗ and

σ (α ) =

√ √ dmax − dmin ; √ √ tmax + d∗ dmax + dmin

tmax − d∗

(b2) for dmax ≥ tmax , or for tmax > dmax √, tmin > dmin and

dmin +tmin dmax +tmax

≥

dmin tmin , dmax tmax

Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

6

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

it holds that α∗ = t∗ and

σ (α ) =

t∗ − dmin t∗ + dmin

√ √ tmax − tmin ; √ √ tmax + tmin

(c) If t∗ < d∗ , the conclusion is similar to (b) and just swap t and d in (b). Actually, the lemma is also true when D and T all are Hermitian positive definite matrices and the new estimate

σ (α ) =

{ max dmin ≤λ1 ≤dmax

|α − λ1 | α + λ1

}

{ max

tmin ≤λ2 ≤tmax

|α − λ2 | α + λ2

}

is a little bit different and the other conclusions and proof processes are exactly the same, see [30]. Hence, the optimal parameter α∗ of GTTS iteration method satisfies the following theorem. Theorem 4.1. The minimal value of σ (α ) in (5) for iteration GTTS method over all positive α is attained at

α∗ = θ∗ =

√

θmax θmin or α∗ = η∗ =

√

ηmax ηmin ,

where θ and η are defined as in Theorem 3.1. The exact cases and proofs are exactly as same as the discussion in Lemma 4.1 (replace d and t to θ and η, respectively). Now, the optimal parameter of α is clear and we consider optimal ω. Before the discussion, we emphasize that we try to find an estimate of ω that can effectively accelerate the convergence of GTTS rather than finding the actual optimal value of ω is hard to get. Note the splitting A = D + T = D + TC + TS of GTTS for Eq. (1), where T is determinate for a certain n as (3) such that TC and TS are uniquely determined by n as well. Moreover, D is a diagonal matrix with nonnegative elements. Because of the special structures of these matrices and actual situation of Eq. (1), we can make some reasonable assumptions so that we can estimate ω pithily. We remark that if the eigenvalues of the matrices TC and TS contain in Ω = [τmin , τmax ], then the fact is demonstrated as follows. Theorem 4.2. Solving Eq. (1) by GTTS methods defined in (5), if the matrices dimension n is sufficiently large and dmin ≫ τmax , for all real ω ∈ [0, 1], then

√ ω∗ = √

√ dmax

√

dmax +

dmin

dmin or ω∗ = √ , √ dmax + dmin

can all be an available estimate, and σ (ω∗ ) is the minimal value of an upper bound of spectral radius of the iteration matrix, where d is the element of D define in (4). Proof. Now we find the estimate of the minimal value of

σ (ω ) =

{ max

θmin ≤θ ≤θmax

|α − θ| α+θ

}

|α − η| α+η

{ max

ηmin ≤η≤ηmax

}

,

where θ ∈ sp(ωD + TC ) and η ∈ sp((1 −ω)D + TS ) are defined in Theorem 3.1(where sp(·) refers to the set of all eigenvalues of a matrix). x According to the monotonically decreasing property of the one variable function φ (x) = α− with respect to x ∈ [0, ∞], α+x we obtain

σ (ω) = max

{(

|α − θmax | |α − θmin | , α + θmax α + θmin

)

( ×

|α − ηmax | |α − ηmin | , α + ηmax α + ηmin

)}

.

Now we expound a fact that the spectral distributions of TC and TS are very close and as n increases their difference is becoming smaller. From the standpoint of spectral distributions, TC and TS may be considered as the same matrix when n is sufficiently large. Note that

θ ∈ sp(ωD + TC ) and η ∈ sp((1 − ω)D + TS ) and D usually has different elements actually so that we use the following approximation to simplify the discussion

{

θmax ≈ ωdmax + τmax θmin ≈ ωdmin + τmin ,

{

ηmax ≈ (1 − ω)dmax + τmax ηmin ≈ (1 − ω)dmin + τmin .

Because D is diagonal, TC and TS are similar, the approximation is available and efficient. Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Therefore, we obtain

) |α − ωdmax − τmax | |α − ωdmin − τmin | σ (ω) = max , α + ωdmax + τmax α + ωdmin + τmin )} ( |α − (1 − ω)dmax − τmax | |α − (1 − ω)dmin − τmin | , . × α + (1 − ω)dmax + τmax α + (1 − ω)dmin + τmin {(

Because the left and right structures of the multiplication are same, we actually consider the minimization of

|α − ωdmax − τmax | |α − (1 − ω)dmax − τmax | |α − ωdmin − τmin | |α − (1 − ω)dmin − τmin | , σ (ω) = max α + ωdmax + τmax α + (1 − ω)dmax + τmax α + ωdmin + τmin α + (1 − ω)dmin + τmin

{

} (6)

In order to get a brief expression, we re-denote the formula above as 1 ⃝} 2 . σ (ω) = max {⃝,

Thus, with the conclusion of Theorem 4.1, the optimal ω that minimizes the function (6) is determined as follows: 1 , then (a) If σ (ω) = ⃝ √ 1 into a block form corresponding (6) as following (a1) for α = θ∗ = (ωdmax + τmax )(ωdmin + τmin ), we can rewrite ⃝ {⏐ a ⏐ ⏐ c ⏐} ⏐ ⏐ ⏐ ⏐ 1 = ⏐ ⏐·⏐ ⏐ . ⃝

b d By concrete computations, verifiably, we obtain the result that only the block c may be equal to zero, a, b, d > 0. In other words, the root of c = 0 is the optimal ω attained at. However, it is a hard task to solve the equations. Recall that we make a hypothesis dmin ≫ τmax and it usually works because τ is really small and almost can be seen as zero. We can ignore τ and obtain the optimal

√

ω∗ = √

dmax

dmin +

√

dmax

.

The above arguments √ readily show that statements in other cases are valid. (a2) for α = η∗ = [(1 − ω)dmax + τmax ][(1 − ω)dmin + τmin ], we obtain the optimal in same way as (a1)

√

ω∗ = √

dmin dmin +

√

dmax

.

2 , same as (a), (b) If σ (ω) = ⃝ (b1) for α = θ∗ , then

√

ω∗ = √

dmin dmin +

√

dmax

.

(b2) for α = η∗ , then

√

ω∗ = √

dmax

dmin +

√

dmax

.

1 and block ⃝ 2 are It is not a coincidence that these two results of ω∗ are obtained in any case. Actually, the block ⃝ symmetrical and have the same ‘‘status’’ for the upper bound σ in some sense. Note that the sum of these two optimal values of ω∗ is 1. In practice, both of them have almost the same effect as well. As we mentioned, the estimate of ω that can effectively accelerate the convergence of GTTS rather than the actual optimal value of ω is more hard to get. In the actual operations of choosing optimal α and ω, we use Theorem 4.2 to get a ω∗ , then we obtain α∗ by Theorem 4.1.

5. Numerical examples In this section, we have solved the FDE (1) by GTTS method and propose two examples to discuss that GTTS methods are more effective than HSS [34], as well as some Krylov subspace iteration methods such as CG [28]. For reference, in the first example, we compare GTTS with DTS [29] method, which is one of the inspiration sources of GTTS method. We use the initial guess u(0) = 0 as the zero vector at each time step, and each iteration process is terminated when the current residual satisfies b − Au(k)  ≤ 10−5 ∥b∥. All numerical experiments are performed in MATLAB (version 8.3.0(R2014a)) with machine precision 10−16 on an Intel Core (CPU 2.50 GHz, 4G RAM) windows 7 system. Actually, in implementations, the discretized linear system that needs to be solved is as follows A = D1 + T and b = △tf 1 + D1 u0 where we take ∆t = h and f ≡ 0. Example 5.1 ([29]). Consider the spatial FDE (1) on domain (0, 1) × (0, +∞) in which u0 = x2 (1 − x), d(x, t) =

1 . x2 (1−x)2

Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

8

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx Table 5.1 The numerical results for Example 5.1 with different choices of parameters for β = 1.2.

β = 1.2 method

Index

N 2^6

2^7

2^8

2^9

2^10

2^11

2^12

2^13

HSS

IT CPU IT CPU IT CPU IT CPU

26 0.0021 22 0.0047 160 0.0033 11 0.0016

37 0.0026 33 0.0057 271 0.0073 14 0.0022

52 0.0030 47 0.0072 425 0.0105 17 0.0028

74 0.0061 70 0.0140 746 0.0732 21 0.0032

104 0.0527 100 0.1090 1222 0.4834 26 0.0145

147 0.3912 147 0.6953 1989 11.925 32 0.0876

208 1.7827 213 3.4041 – – 38 0.3722

293 8.4364 320 17.9210 – – 48 1.3358

CG DTS GTTS

Table 5.2 The numerical results for Example 5.1 with different choices of parameters for β = 1.8.

β = 1.8 method

Index

N 2^6

2^7

2^8

2^9

2^10

2^11

2^12

2^13

HSS

IT CPU IT CPU IT CPU

– – 21 0.0049 2133 0.0339

– – 32 0.0066 18 0.0018

96 0.0048 48 0.0072 13 0.0020

88 0.0077 68 0.0141 17 0.0040

92 0.0467 101 0.1046 21 0.0208

105 0.2916 149 0.7224 27 0.0753

151 1.3066 213 3.3421 33 0.2968

225 6.5225 318 17.7199 41 1.1496

CG GTTS

Fig. 5.1. The curves of ρ (LHSS (α )) and ρ (LGTTS (α )), and the theoretical optimal parameters α of HSS and GTTS are denoted by pints ‘o’ and ‘∗’ respectively.

In Tables 5.1 and 5.2, we report the number of iteration steps and the CPU time of the methods used in Example 5.1, and observe that GTTS is still outstanding, especially when the matrix is larger. The theoretical optimal parameters α are shown in Fig. 5.1 as well, and their values are really close to the experimental optimal α values. Note that the span of α is large in Fig. 5.1, we zoom in the area near the point of theoretical optimal parameters α as Fig. 5.2 and thereby we can observe more clearly. As shown in Fig. 5.2 that the theoretical optimal parameters α almost reach the experimental one, it is a powerful evidence that our estimated parameters are valid. Now we consider the parameters ω and we plot the function of ρ (LGTTS (ω)) in Fig. 5.3, when β = 1.2, N = 29 and optimal α (ω) according to ω. Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Fig. 5.2. The area, from Fig. 5.1, near the theoretical optimal parameters α .

Fig. 5.3. The curves of ρ (LGTTS (ω)), and the theoretical exact optimal parameters ω of GTTS are denoted by pints ‘∗’ and the other kind of optimal parameters are denoted by ‘+’.

In Section 4, we get two possible parameters ω,

√

ω∗ = √

√

dmin dmin +

√

dmax

and ω∗ = √

dmax

dmin +

√

dmax

denoted by ‘∗’ and ‘+’ respectively, and the first one should be chosen in these cases. Obviously, the two parameters add up to equal 1, hence they look symmetrical in Fig. 5.3. Actually, they have almost the same effect on the optimization of the spectral radius ρ (LGTTS (ω)), although the former has a slight advantage. In the proof of the optimal parameters ω, we ignore some factors for simplicity and make a hypothesis that N is large enough, which generate the error between theoretical and experimental optimal parameters. Nevertheless, the error is acceptable. In the next experiment we will see this fact according to the hypothesis that the larger the N, the smaller the error. Example 5.2 ([29]). Consider the spatial FDE (1) on domain (0, 1) × (0, +∞) in which u0 = x(1 − x), d(x, t) =

512 . x3 (1+8x)3

In Table 5.3, we report the number of iteration steps and the CPU time of the methods used in Example 5.2, and observe that GTTS is still outstanding, especially when the matrix is larger. Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

10

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx Table 5.3 The numerical results for Example 5.2 with different choices of parameters for β = 1.2.

β = 1.2 method

Index

N 2^6

2^7

2^8

2^9

2^10

2^11

2^12

2^13

HSS

IT CPU IT CPU IT CPU

95 0.0047 63 0.0076 18 0.0019

124 0.0060 127 0.0106 22 0.0031

161 0.0082 255 0.0206 26 0.0028

210 0.0160 458 0.0652 31 0.0148

284 0.1361 706 0.6968 38 0.0203

393 0.9797 2047 9.5356 46 0.1229

553 4.6689 4096 62.8932 55 0.4660

780 22.3664 8191 450.5419 67 1.9482

CG GTTS

Fig. 5.4. The curves of ρ (LGTTS (ω)), three functions of N = 27 , 29 , 210 according to ω denoted by ‘-’, ‘–’, ‘-.’ respectively, and the respective theoretical optimal parameters ω are denoted by pints ‘∗’.

Fig. 5.5. The area, from Fig. 5.4, near the theoretical optimal parameters ω.

Of course, one major reason that GTTS exhibits preferable numerical behavior when matrix is larger is that, at the same time, theoretical optimal parameters ω are closer to the experimental one, see Fig. 5.5, which is zoomed to show the optimal points of Fig. 5.4. The consequence is expected. Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

11

6. Concluding remarks In this paper, we propose a GTTS method for solving diagonal-plus-Toeplitz linear systems from spatial FDEs. We show the conditions that the GTTS methods converge to the unique solution. Then we analyze and derive an upper bound of the contraction factor of the GTTS. As the part we highlighted in this paper, we consider the optimizations of two parameters in GTTS, which have a huge impact on convergence rate. Based on reasonable assumptions, we have concise optimal parameters, which are both easy to get computationally. At last, numerical experiments show that the GTTS method is more effective than the HSS and spectral radius of iteration matrix is much lower. Compared with the Krylov subspace iteration methods such as CG, the GTTS is more effective, especially when the matrix is larger. In a word, the GTTS has a satisfactory experimental result. However, as shown in the previous sections, in the proof of the optimal parameters ω, we ignore some factors to simplify the process and make some hypotheses, which make a error between theoretical and experimental optimal parameters ω. Moreover, we are not considering the preconditioner of GTTS. It is important to solve those two problems however they are both very challenging. Hence, we believe there should be more in-depth study from the viewpoint of theory and computation in these areas. References [1] S.-G. Samko, A.-A. Kilbas, O.-I. Marichev, Fractional integrals and derivatives, Yverdon-Les-Bains, Gordon and Breach Science Publishers, Yverdon, Switzerland, 1993. [2] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations to Methods of their Solution and Some of their Applications, Elsevier, 1998. [3] R.-L. Magin, Fractional Calculus in Bioengineering, Begell House, Redding, 2006. [4] R.-K. Upadhyay, Dynamics of fractional order modified Morris-Lecar neural model, Netw. Biol. 5 (3) (2015) 113–136. [5] J.-W. Kirchner, X. Feng, C. Neal, Fractal stream chemistry and its implications for contaminant transport in catchments, Nature 403 (6769) (2000) 524–527. [6] M. Kreer, A. Kızılersü, A.-W. Thomas, Fractional Poisson processes and their representation by infinite systems of ordinary differential equations, Statist. Probab. Lett. 84 (2014) 27–32. [7] I.-M. Sokolov, J. Klafter, A. Blumen, Fractional kinetics, Phys. Today 55 (11) (2002) 48–54. [8] B.-A. Carreras, V.-E. Lynch, G.-M. Zaslavsky, Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence model, Phys. Plasmas 8 (12) (2001) 5096–5103. [9] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial data: an empirical study, Physica A 314 (1–4) (2002) 749–755. [10] L. Sabatelli, S. Keating, J. Dudley, et al., Waiting time distributions in financial markets, Eur. Phys. J. B 27 (2) (2002) 273–275. [11] S. Shen, F. Liu, V. Anh, et al., The fundamental solution and numerical solution of the Riesz fractional advection–dispersion equation, IMA J. Appl. Math. 73 (6) (2008) 850–872. [12] G.-H. Gao, Z.-Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (3) (2011) 586–595. [13] F. Liu, P. Zhuang, K. Burrage, Numerical methods and analysis for a class of fractional advection–dispersion models, Comput. Math. Appl. 64 (10) (2012) 2990–3007. [14] B. Baeumer, M. Kovács, M.-M. Meerschaert, Numerical solutions for fractional reaction–diffusion equations, Comput. Math. Appl. 55 (10) (2008) 2212–2226. [15] W.-H. Deng, Finite element method for the space and time fractional Fokker–Planck equation, SIAM J. Numer. Anal. 47 (1) (2008) 204–226. [16] C. Tadjeran, M.-M. Meerschaert, H.-P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (1) (2006) 205–213. [17] M.-M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Numer. Math. 56 (1) (2006) 80–90. [18] K.-X. Wang, H. Wang, A fast characteristic finite difference method for fractional advection–diffusion equations, Adv. Water Resour. 34 (7) (2011) 810–816. [19] M.-M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (1) (2004) 65–77. [20] I. Podlubny, A. Chechkin, T. Skovranek, et al., Matrix approach to discrete fractional calculus II: Partial fractional differential equations, J. Comput. Phys. 228 (8) (2009) 3137–3153. [21] H. Wang, K.-X. Wang, T. Sircar, A direct O (N log2 N) finite difference method for fractional diffusion equations, J. Comput. Phys. 229 (21) (2010) 8095–8104. [22] H. Wang, K.-X. Wang, An O (N log2N) alternating-direction finite difference method for two-dimensional fractional diffusion equations, J. Comput. Phys. 230 (21) (2011) 7830–7839. [23] S.-L. Lei, H.-W. Sun, A circulant preconditioner for fractional diffusion equations, J. Comput. Phys. 242 (2013) 715–725. [24] J.-Y. Pan, R.-H. Ke, M.-K. Ng, et al., Preconditioning techniques for diagonal-times-Toeplitz matrices in fractional diffusion equations, SIAM J. Sci. Comput. 36 (6) (2014) A2698–A2719. [25] F.-R. Lin, S.-W. Yang, X.-Q. Jin, Preconditioned iterative methods for fractional diffusion equation, J. Comput. Phys. 256 (2014) 109–117. [26] Pedro Alonso, Manuel F. Dolz, et Vida, Block pivoting implementation of a symmetric Toeplitz solver, J. Parallel Distrib. Comput. 74 (5) (2014) 2392–2399. [27] H.-G. Golub, C.-F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore and London, 1996. [28] A.-J. Wathen, Preconditioning, Acta Numer. 24 (2015) 329–376. [29] Z.-Z. Bai, K.-Y. Lu, 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 (4) (2017) 2093. [30] Z.-Y. Liu, N.-C. Wu, X.-R. Qin, et al., Trigonometric transform splitting methods for real symmetric Toeplitz systems, Comput. Math. Appl. 75 (8) (2018) 2782–2794. [31] G. Heinig, K. Rost, Representations of Toeplitz-plus-Hankel matrices using trigonometric transformations with application to fast matrix–vector multiplication, Linear Algebra Appl. 275 (1998) 225–248.

Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.

12

X.-H. Shao, Z.-D. Zhang and H.-L. Shen / Computers and Mathematics with Applications xxx (xxxx) xxx

[32] T. Huckle, Iterative methods for ill-conditioned Toeplitz matrices, Calcolo 33 (3–4) (1996) 177–190. [33] D.-A. Benson, S.-W. Wheatcraft, M.-M. Meerschaert, Application of a fractional advection–dispersion equation, Water Resour. Res. 36 (6) (2000) 1403–1412. [34] Z.-Z. Bai, G.-H. Golub, M.-K. Ng, Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, SIAM J. Matrix Anal. Appl. 24 (3) (2003) 603–626.

Please cite this article as: X.-H. Shao, Z.-D. Zhang and H.-L. Shen, A generalization of trigonometric transform splitting methods for spatial fractional diffusion equations, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.003.