Spectral analysis for preconditioning of multi-dimensional Riesz fractional diffusion equations
SSPECTRAL ANALYSIS FOR PRECONDITIONING OFMULTI-DIMENSIONAL RIESZ FRACTIONAL DIFFUSIONEQUATIONS ∗ XIN HUANG † , XUE-LEI LIN ‡ , MICHAEL K. NG § , AND
HAI-WEI SUN ¶ Abstract.
In this paper, we analyze the spectra of the preconditioned matrices arising from dis-cretized multi-dimensional Riesz spatial fractional diffusion equations. The finite difference methodis employed to approximate the multi-dimensional Riesz fractional derivatives, which will generatesymmetric positive definite ill-conditioned multi-level Toeplitz matrices. The preconditioned conju-gate gradient method with a preconditioner based on the sine transform is employed to solve theresulting linear system. Theoretically, we prove that the spectra of the preconditioned matrices areuniformly bounded in the open interval (1 / , /
2) and thus the preconditioned conjugate gradientmethod converges linearly. The proposed method can be extended to multi-level Toeplitz matricesgenerated by functions with zeros of fractional order. Our theoretical results fill in a vacancy inthe literature. Numerical examples are presented to demonstrate our new theoretical results in theliterature and show the convergence performance of the proposed preconditioner that is better thanother existing preconditioners.
Key words.
Riesz fractional derivative, multi-level Toeplitz matrix, sine transform based pre-conditioner, condition number, fractional order zero, preconditioned conjugate gradient method
AMS subject classifications.
1. Introduction.
In this paper, we study the preconditioning technique for thefollowing multi-dimensional Riesz fractional diffusion equations(1.1) − m (cid:88) i =1 d i ∂ α i u ( x ) ∂ | x i | α i = y ( x ) , x ∈ Ω = m (cid:89) i =1 [ a i , b i ] ⊂ R m , subject to the boundary condition u ( x ) = 0 , x ∈ ∂ Ω , where d i > i = 1 , . . . , m , x = ( x , . . . , x m ) ∈ R m , y ( x ) : R m (cid:55)→ R is the sourceterm, and ∂ αi u ( x ) ∂ | x i | αi is the Riesz fractional derivative of α i ∈ (1 ,
2) with respect to x i defined by(1.2) ∂ α i u ( x ) ∂ | x i | α i = c ( α i ) (cid:0) a i D α i x i u ( x ) + x i D α i b i u ( x ) (cid:1) , c ( α i ) = −
12 cos( α i π ) > . The above left and right Riemann-Liouville (RL) fractional derivatives are defined by(1.3) a i D α i x i u ( x ) = 1Γ(2 − α i ) ∂ ∂x i (cid:90) x i a i u ( x , x , . . . , x i − , ξ, x i +1 , . . . , x m )( x i − ξ ) α i − d ξ, ∗ Submitted to the editors DATE.
Funding:
This work is supported by research grants of the Science and Technology DevelopmentFund, Macau SAR (file no. 0118/2018/A3), MYRG2018-00015-FST from University of Macau, andthe HKRGC GRF 12306616, 12200317, 12300218, 12300519 and 17201020. † Department of Mathematics, University of Macau, Macao ([email protected]). ‡ Department of Mathematics, University of Macau, Macao ([email protected]). § Corresponding author. Department of Mathematics, The University of Hongkong, Hong Kong([email protected]). ¶ Department of Mathematics, University of Macau, Macao ([email protected]).1
This manuscript is for review purposes only. a r X i v : . [ m a t h . NA ] F e b XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN (1.4) x i D α i b i u ( x ) = 1Γ(2 − α i ) ∂ ∂x i (cid:90) b i x i u ( x , x , . . . , x i − , ξ, x i +1 , . . . , x m )( ξ − x i ) α i − d ξ, respectively, where Γ( · ) is the gamma function.Fractional calculus has received an increasing interest since its applications in-volve various fields including physics, chemistry, engineering; see [11, 17, 18, 28]. TheRiesz fractional derivative, which derives from the kinetic of chaotic dynamics [31],is generalized as one of the most popular fractional calculus. Recently, the study ofthe Riesz fractional derivative has been urgent and significant as it can be applied tolattice model with long-range interactions [12], nonlocal dynamics [36], and so on.After discretization by the finite difference method, the resulting coefficient matrixof the above multi-dimensional Riesz fractional derivative in (1.1) is a dense multi-level Toeplitz matrix; i.e., each block has a Toeplitz structure. It is interesting to notethat such multi-level Toeplitz matrix can be generated by a continuous real-valuedeven function which is nonnegative defined on the interval [ − π, π ] [26]. Moreover, thediagonal entries are the Fourier coefficients of the generating function with α -th orderzero at the origin (1 < α < n α , where n denotes the matrix size; see [6, 7, 33]. Therefore, theresulting linear system arising from (1.1) is ill-conditioned and thus the conjugategradient (CG) method for this system converges slowly.In order to speed up the convergence of the CG method, preconditioning tech-niques have been proposed and developed for ill-conditioned Toeplitz systems. For ex-amples, banded Toeplitz preconditioners [3,7] were proposed to handle ill-conditionedToeplitz linear systems where the generating functions of Toeplitz matrices have zerosof even order. When these banded Toeplitz preconditioners are applied to Toeplitzmatrices generated by functions with α -th order zero at the origin (1 < α < τ -preconditioners (which can be diagonalized by the discrete sinetransform matrix) [2, 4, 9, 34], circulant preconditioners [10, 24, 29, 30], and multigridmethods [3, 8, 15, 16, 32, 35]. These approaches can significantly speed up the con-vergence of iterative methods for solving Toeplitz systems where their generatingfunctions have an α -th order zero at the origin (1 < α < theoretically confirmed for solving such Toeplitzsystems. Lately, some efficient preconditioners are developed to solve linear systemsarising from fractional diffusion equations; see [1,13,14,21,23]. Nevertheless, from thetheoretical point of view, the spectra of these preconditioned matrices are not shown tobe bounded independent of the matrix sizes and hence the linear convergence cannotbe guaranteed.In order to tackle this theoretical problem, Noutsos, Serra and Vassalos [25] ex-ploited a multiple step preconditioning and applied to the case where the generatingfunctions of coefficient matrices have fractional order zeros. In their method, theyproposed a new τ -preconditioner that is constructed from the generating function ofthe given Toeplitz matrix. Numerical results were shown that their method workedvery well for the Toeplitz matrices whose generating functions have fractional orderzeros. Theoretically, they have proved that the largest eigenvalue of the precondi-tioned matrix has an upper bound independent of the matrix size. Nevertheless, it isstill unclear whether the smallest eigenvalue of the preconditioned matrix is boundedbelow away from zero; see the remark in [25]. This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE τ -preconditionerfor the ill-conditioned multi-level Toeplitz system arising from the discretized Rieszfractional derivatives. Theoretically, we prove that the spectra of the τ -preconditionedmatrices are uniformly bounded in the open interval (1 / , /
2) and thus the precon-ditioned CG (PCG) method converges linearly. Furthermore, the proposed methodcan be extended to multi-level Toeplitz matrices generated by functions with zeros offractional order. Similarly, we show that the spectra of these preconditioned multi-level Toeplitz matrices are uniformly bounded. Numerical examples are presented toverify our new theoretical results in the literature and show the good performance ofthe proposed preconditioner.The outline of the rest paper is as follows. In Section 2, multi-level Toeplitzmatrices are generated. In Section 3, the preconditioner is proposed and developed.In Section 4, the spectral analysis of the preconditioned matrix is discussed. InSection 5, a new preconditioning technique for multi-level Toeplitz matrices is studied.Numerical experiments are given in Section 6 to show the performance of the proposedpreconditioner. Finally, some concluding remarks are given in Section 7.
2. Multi-level Toeplitz matrices.
A matrix, whose entries are constant alongthe diagonals with the following form(2.1) T n = t t − . . . t − n t − n t t t − . . . t − n ... t t . . . ... t n − . . . . . . . . . t − t n − t n − . . . t t , is called a Toeplitz matrix. Assume that the diagonals { t k } n − k =1 − n of T n are the Fouriercoefficients of a function f ; i.e., t k = 12 π (cid:90) π − π f ( θ ) e − i kθ dθ, then the function f is called the generating function of T n . Generally, denote T n = T n ( f ) to emphasize that an n × n Toeplitz matrix T n is generated by f . Moreover, ifthe generating function is real and even, then T n is real symmetric for all n .Let θ = ( θ , θ , . . . , θ m ) ∈ [ − π, π ] m , f ( θ ) = f ( θ , . . . , θ m ) ∈ L ([ − π, π ] m ). Thematrix generated by function f ( θ ) is an m -level Toeplitz matrix with Toeplitz struc-ture on each level. Let n i for i = 1 , , . . . , m be positive integers. Denote N = m (cid:81) i =1 n i .Then an m -level Toeplitz matrix with the size N × N is a block Toeplitz matrix withToeplitz block, and its form is(2.2) T mN = T m − T m − − · · · T m − − n m T m − T m − . . . ...... . . . . . . T m − − T m − n m − · · · T m − T m − , where each block T m − j for j = 0 , . . . , n m − m − n · · · n m − ) × ( n · · · n m − ) with ( m − n · · · n m − ) × This manuscript is for review purposes only.
XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN ( n · · · n m − ), and so on. Let ( j , . . . , j m ) ∈ Z m be a multi-index. Then the coefficientsof T mN can be obtained by [27](2.3) t j ,j ,...,j m = 1(2 π ) m (cid:90) [ − π,π ] m f ( θ ) e − i m (cid:80) i =1 θ i j i d θ . Firstly, we consider theone dimensional case. Let n be the partition of the interval [ a, b ]. We define a uniformspatial partition h = b − an + 1 , η j = a + jh, for j = 0 , , . . . , n + 1 . Then the following shifted Gr¨unwald-Letnikov formula is exploited to approximatethe left- and right- RL fractional derivatives at grid point η j ,(2.4) a D αx u ( η j ) = 1 h α j +1 (cid:88) k =0 g ( α ) k u ( η j − k +1 ) + O ( h ) , (2.5) x D αb u ( η j ) = 1 h α n − j +2 (cid:88) k =0 g ( α ) k u ( η j + k − ) + O ( h ) , where the coefficients g ( α ) k are defined by(2.6) g ( α )0 = 1 , g ( α ) k = (cid:18) − α + 1 k (cid:19) g ( α ) k − , for k ≥ . It can be shown that the coefficients g ( α ) k defined above have the following properties. Lemma
For α ∈ (1 , , the coefficients g ( α ) k , k = 0 , , . . . , satisfy g ( α )0 = 1 , g ( α )1 = − α < , g ( α )2 > g ( α )3 > · · · > , ∞ (cid:88) k =0 g ( α ) k = 0 , n (cid:88) k =0 g ( α ) k < , for n ≥ . By applying (2.4) and (2.5) to the model equation (1.1), we obtain(2.7) − d c ( α ) h α (cid:32) j +1 (cid:88) k =0 g ( α ) k u ( η j − k +1 ) + n − j +2 (cid:88) k =0 g ( α ) k u ( η j + k − ) (cid:33) = y ( η j ) + O ( h ) . Defining u j as the numerical approximation of u ( η j ), setting y j = y ( η j ), and omittingthe small term O ( h ), the finite difference scheme for solving (1.1) is constructed asfollows(2.8) − d c ( α ) h α (cid:32) j +1 (cid:88) k =0 g ( α ) k u j − k +1 + n − j +2 (cid:88) k =0 g ( α ) k u j + k − (cid:33) = y j , for j = 1 , . . . , n. Let u = [ u , . . . , u n ] (cid:124) , y = [ y , . . . , y n ] (cid:124) . Then, the numerical scheme (2.8) can besimplified as the following matrix-vector form(2.9) A n u = y, This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE A n = wG ( α ) n , where w = dc ( α ) h α > G ( α ) n = − g ( α )1 g ( α )0 + g ( α )2 g ( α )3 . . . g ( α ) n − g ( α ) n g ( α )0 + g ( α )2 g ( α )1 g ( α )0 + g ( α )2 g ( α )3 . . . g ( α ) n − ... g ( α )0 + g ( α )2 g ( α )1 . . . . . . ...... . . . . . . . . . . . . g ( α )3 g ( α ) n − . . . . . . . . . 2 g ( α )1 g ( α )0 + g ( α )2 g ( α ) n g ( α ) n − · · · · · · g ( α )0 + g ( α )2 g ( α )1 . It is evident that the matrix G ( α ) n is a symmetric Toeplitz matrix, which is generated byan integrable real-value even function defined on the interval [ − π, π ]. In the following,the generating function of matrix G ( α ) n will be presented. Lemma
The generating function of matrix G ( α ) n is (2.11) g α ( θ ) = (cid:40) − α +1 (sin − θ ) α cos[ α ( π + θ ) − θ ] , θ ∈ [ − π, , − α +1 (sin θ ) α cos[ α ( π − θ ) + θ ] , θ ∈ [0 , π ] . By Lemma 2.2, the generating function of G ( α ) n possesses the following properties,which will be exploited to further study. Lemma
For any α ∈ (1 , , it holds that ≤ | θ | α g α ( θ ) ≤ π − πα ) , where θ ∈ [ − π, π ] . Proof.
We only consider the case that θ ∈ (0 , π ] as it is an even function. FromLemma 2.2, we have | θ | α g α ( θ ) = − θ α α +1 (sin θ ) α cos[ α ( π − θ ) + θ ] . Note that 1 < α <
2, it holds12 < θ α α +1 (sin θ ) α = 12 ( θ ) α (sin θ ) α ≤ π , and 1 ≤ − α ( π − θ ) + θ ] ≤ − πα ) . Similarly, we can derive the same conclusion for θ ∈ [ − π,
0) and the result is con-cluded.In the light of the definition of the fractional order zero defined in [14], we deducethat the generating function g α ( θ ) has a zero of order α at θ = 0. Now, we extend ourinvestigation to the multi-dimensional cases as in (1.1). To obtain the discretized formof multi-dimensional Riesz fractional diffusion equations, some notations are required.
This manuscript is for review purposes only.
XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
Denote I k be a k × k identity matrix. Let n − i = i − (cid:81) j =1 n j and n + i = m (cid:81) j = i +1 n j for i = 1 , , . . . , m , respectively. In particular, take n − = n + m = 1. Let h i = b i − a i n i +1 , for i = 1 , . . . , m , η ij = a i + jh i , and y j ,j ,...,j m = y ( η j , η j , . . . , η mj m ). Denote u = [ u , ,..., , . . . , u n , ,..., , u , ,..., , . . . , . . . , u n ,n ,...,n m ] (cid:124) ∈ R N and y = [ y , ,..., , . . . , y n , ,..., , y , ,..., , . . . , . . . , y n ,n ,...,n m ] (cid:124) ∈ R N . Analogously, using the shifted Gr¨unwald-Letnikov formula to discretize the multi-dimensional Riesz fractional diffusion equations (1.1), we obtain the matrix-vectorform of the resulting linear system as(2.12) Au = y , where(2.13) A = m (cid:88) i =1 I n − i ⊗ A n i ⊗ I n + i , in which A n i = w i G ( α i ) n i , w i = d i c ( α i ) h αii >
0, and G ( α i ) n i is defined as in (2.10). Moreover,the generating function of A is as follows,(2.14) f α ( θ ) = m (cid:88) i =1 w i g α i ( θ i ) , where α = ( α , . . . , α m ).Note that the generating function g α i ( θ i ) is nonnegative and hence so is f α ( θ ),which manifests that the coefficient matrices A n i and A are both symmetric positivedefinite. It is well-known that the CG method has been recognized as the most appro-priate method for solving symmetric positive definite Toeplitz systems. Actually, sincethe generating functions have zeros, the corresponding matrices are ill-conditioned,which will slow down the convergent rate of the CG method. In order to speed up theconvergent rate, the preconditioning technique is employed to solve the ill-conditionedsystem.
3. Sine transform based preconditioners.
Let T n be a given symmetricToeplitz matrix whose first column is [ t , t , . . . , t n − ] (cid:124) . Then, the natural τ matrix τ ( T n ) of T n can be determined by the Hankel correction [4](3.1) τ ( T n ) = T n − H n , where H n is a Hankel matrix whose entries are constants along the antidiagonals, inwhich the antidiagonals are given as[ t , t , . . . , t n − , , , , t n − , . . . , t , t ] (cid:124) . More preciously, the entries p ij of τ ( T n ) can be generalized as(3.2) p ij = t | j − i | − t i + j , i + j < n − ,t | j − i | , i + j = n − , n, n + 1 ,t | j − i | − t n +2 − ( i + j ) , otherwise . This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE τ matrix can be diagonalized by the sine transform matrix S n [5];i.e.,(3.3) τ ( T n ) = S n Λ n S n , where Λ n is a diagonal matrix holding all the eigenvalues of τ ( T n ) and the entries of S n are given by(3.4) [ S n ] j,k = (cid:114) n + 1 sin (cid:18) πjkn + 1 (cid:19) , ≤ k, j ≤ n. Obviously, S n is a symmetric orthogonal matrix, and the matrix-vector multiplication S n v for any vector v can be computed with only O ( n log n ) operations by the fast sinetransform. Likewise, the matrix-vector product τ ( T n ) − v = S n Λ − n S n v can be donein O ( n log n ) operations. Moreover, the eigenvalues of the τ matrix can be determinedby its first column. Thus, O ( n ) storage and O ( n log n ) computational complexity arerequired for saving and computing the eigenvalues of the τ matrix, respectively.Now we consider the preconditioner matrix of the linear system (2.9). Recallingthe form of the coefficient matrix, we obtain the preconditioner as(3.5) P n = τ ( A n ) = wτ ( G ( α ) n ) , where G ( α ) n is the symmetric positive definite Toeplitz matrix defined in (2.10). Forgeneral τ matrix, we obtain its eigenvalues by the following lemma. Lemma
Let T n be a symmetric Toeplitz matrix whose first columnis [ t , t , . . . , t n − ] (cid:124) . Then the eigenvalues of τ ( T n ) can be expressed as σ j = t + 2 n − (cid:88) k =1 t k cos( jζ k ) , j = 1 , . . . , n, where ζ k = πkn +1 , k = 1 , . . . , n − . Lemma 3.1 provides a methodology to compute the eigenvalues of the τ matrix.Taking advantage of Lemma 3.1, the preconditioner P n can be verified to be positivedefinite in the following lemma. Lemma
The preconditioner P n = wτ ( G ( α ) n ) defined in (3.5) is symmetricpositive definite.Proof. The first column of G ( α ) n in (2.10) is (cid:104) − g ( α )1 , − g ( α )0 − g ( α )2 , − g ( α )3 , . . . , − g ( α ) n (cid:105) (cid:124) . According to Lemma 3.1, the j th eigenvalue of τ ( G ( α ) n ) can be expressed as δ j = − g ( α )1 − (cid:16) g ( α )0 + g ( α )2 (cid:17) cos (cid:18) πjn + 1 (cid:19) − n (cid:88) k =3 g ( α ) k cos (cid:18) πj ( k − n + 1 (cid:19) ≥ − g ( α )1 − (cid:16) g ( α )0 + g ( α )2 (cid:17) − n (cid:88) k =3 g ( α ) k = − n (cid:88) k =0 g ( α ) k > . (by Lemma 2.1) This manuscript is for review purposes only.
XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
Therefore, we deduce that τ ( G ( α ) n ) is symmetric positive definite and conclude thatso is P n .Making use of Lemma 3.2, we derive the following corollary. Corollary
The preconditioner P n defined in (3.5) is invertible. Now we consider the τ -preconditioner for multi-level Toeplitz matrices. Recallthe coefficient matrix A = m (cid:88) i =1 I n − i ⊗ A n i ⊗ I n + i . The preconditioner matrix of the linear system (2.12) can be expressed as(3.6) P = m (cid:88) i =1 I n − i ⊗ τ ( A n i ) ⊗ I n + i , where τ ( A n i ) = w i τ ( G ( α i ) n i ) for i = 1 , . . . , m . Denote(3.7) A i = I n − i ⊗ A n i ⊗ I n + i , and the corresponding τ -matrix can be defined as τ ( A i ) = I n − i ⊗ τ ( A n i ) ⊗ I n + i . Lemma
The preconditioner P defined in (3.6) is symmetric positive definite.Proof. Let Λ n i be a diagonal matrix holding the eigenvalues of τ ( A n i ). Then τ ( A n i ) = S n i Λ n i S n i . Lemma 3.2 has shown that all the eigenvalues of τ ( G ( α i ) n i ) arepositive; that is, τ ( A n i ) for i = 1 , . . . , m are positive definite. By virtue of theproperties of Kronecker product, it is easy to obtain I n − i ⊗ τ ( A n i ) ⊗ I n + i = S ( I n − i ⊗ Λ n i ⊗ I n + i ) S , where S = m (cid:78) i =1 S n i , which manifests that all the eigenvalues of τ ( A i ) are positive. As P is the summation of symmetric positive definite matrices, it follows that P is alsosymmetric positive definite.
4. Spectral analysis for preconditioned matrices.
In this section, we dis-cuss the spectra of the preconditioned matrices.
In the following, we discuss the spectrum of τ ( G ( α ) n ) − H n , which is needed in the theoretical analysis.First of all, we focus on the matrix τ ( G ( α ) n ) = G ( α ) n − H n , where H n is the Hankelmatrix as in (3.1), and G ( α ) n is defined in (2.10). Denote the first column of G ( α ) n as(4.1) [ t , t , . . . , t n − ] (cid:124) = [ − g ( α )1 , − ( g ( α )0 + g ( α )2 ) , . . . , − g ( α ) n ] (cid:124) . By Lemma 2.1, we have
Lemma
Let t j , j = 0 , , . . . , n − , be defined in (4.1) . Then t > , t < t < · · · < t n − < , This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE and t + 2 n − (cid:88) j =1 t j > . Proof.
Using (2.6) and Lemma 2.1, the results are immediately concluded.Let h ij and p ij be the entries of H n and τ ( G ( α ) n ), respectively. From Lemma 4.1,we have(4.2) h ij = t i + j , i + j < n − , , i + j = n − , n, n + 1 ,t n +2 − ( i + j ) , otherwise , and(4.3) p ij = t | i − j | − h ij . It is obvious that h ij ≤ ≤ i, j ≤ n . By Lemma 4.1, we will show that p ij possesses the following properties. Lemma
Let p ij be the entries of τ ( G ( α ) n ) defined in (4.3) . For ≤ i, j ≤ n ,it holds that (4.4) p ii > , and p ij < , for i (cid:54) = j. Proof.
It is evident that p ii = t − h ii > i = 1 , . . . , n . On the otherhand, let 1 ≤ i < j ≤ n . It is easy to check that | i − j | < i + j, and | i − j | = | ( n − i ) − ( n − j ) | < n − ( i + j ) + 2 . By Lemma 4.1 and formulas (4.2)–(4.3), we derive that p ij = t | i − j | − t i + j ( or t | i − j | − t n − ( i + j )+2 ) < . The proof is complete.
Lemma
Let G ( α ) n be the Toeplitz matrix defined in (2.10) , and H n be thecorresponding Hankel matrix whose entries are defined in (4.2) . Then, the eigenvaluesof τ ( G ( α ) n ) − H n fall inside the open interval ( − / , / .Proof. Let λ be the eigenvalue of τ ( G ( α ) n ) − H n , and z = [ z , z , . . . , z n ] (cid:124) be thecorresponding eigenvector with the largest component having the magnitude 1; i.e.,max ≤ j ≤ n | z j | = 1. Then, we have H n z = λτ ( G ( α ) n ) z. For each i , it holds that n (cid:88) j =1 h ij z j = λ n (cid:88) j =1 p ij z j , which can be written as λp ii z i = n (cid:88) j =1 h ij z j − λ n (cid:88) j =1 ,j (cid:54) = i p ij z j . This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
Let | z k | = 1. Then, by the above formula, it follows | λ || p kk | ≤ n (cid:88) j =1 | h kj | + | λ | n (cid:88) j =1 ,j (cid:54) = k | p kj | . Accordingly, it is resulted that | λ | ≤ n (cid:80) j =1 | h kj || p kk | − n (cid:80) j =1 ,j (cid:54) = k | p kj | . By Lemma 4.1 and Lemma 4.2, we have | p kk | − n (cid:88) j =1 ,j (cid:54) = k | p kj | − n (cid:88) j =1 | h kj | =( t − h kk ) − n (cid:88) j =1 ,j (cid:54) = k ( h kj − t | k − j | ) + 2 n (cid:88) j =1 h kj (by (4.3) and (4.4))= t + n (cid:88) j =1 ,j (cid:54) = k t | k − j | + n (cid:88) j =1 h kj = t + k − (cid:88) j =1 t j + n − k (cid:88) j =1 t j + n − (cid:88) j = k +1 t j + n − (cid:88) j = n − k +2 t j ≥ t + 2 n − (cid:88) j =1 t j > , which manifests | λ | < /
2. Therefore, we derive that the eigenvalues of τ ( G ( α ) n ) − H n fall inside the interval ( − / , / τ ( G ( α ) n ) − H n are bounded, whichis the key to obtain the spectrum of τ ( G ( α ) n ) − G ( α ) n . Theorem
Let λ ( τ ( G ( α ) n ) − G ( α ) n ) be the eigenvalues of matrix τ ( G ( α ) n ) − G ( α ) n .Then the following inequality holds < λ ( τ ( G ( α ) n ) − G ( α ) n ) < . Proof.
Using I n + τ ( G ( α ) n ) − H n = τ ( G ( α ) n ) − ( τ ( G ( α ) n ) + H n ) = τ ( G ( α ) n ) − G ( α ) n , taking advantage of the conclusion obtained in Lemma 4.3, the proof is complete.Based on Theorem 4.4, the following corollary which gives an upper bound interms of the condition number of the preconditioned matrix can be achieved. Corollary
Let P n defined in (3.5) be the preconditioner of the linear system (2.9) . Then, the condition number of the preconditioned matrix P − n A n is less than .This manuscript is for review purposes only. PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE Proof.
It is easy to verify that P − n A n = τ ( G ( α ) n ) − G ( α ) n . By Theorem 4.4, we obtain κ = λ max ( P − n A n ) λ min ( P − n A n ) < . Corollary 4.5 manifests that the spectrum of the preconditioned matrix is uni-formly bounded independent of the matrix size, where the smallest eigenvalue is awayfrom 0. Therefore, we conclude that the τ -preconditioner is efficient for 1D Riesz frac-tional diffusion equations. Next, we extend our discussion on the multi-dimensionalcases. Recalling the coefficient matrix A defined in(2.13), the matrix A n i is the coefficient matrix corresponding to 1D case. In the lightof the results shown in previous subsection, we obtain the spectrum of the matrix τ ( A i ) − A i firstly. Lemma A i is a multi-level Toeplitz matrix defined in (3.7) . We then havefor each i = 1 , . . . , m , the eigenvalues of τ ( A i ) − A i satisfying < λ ( τ ( A i ) − A i ) < . Proof.
It is clear that τ ( A i ) − A i = (cid:16) I n − i ⊗ τ ( A n i ) ⊗ I n + i (cid:17) − (cid:16) I n − i ⊗ A n i ⊗ I n + i (cid:17) = (cid:16) I − n − i ⊗ τ ( A n i ) − ⊗ I − n + i (cid:17) (cid:16) I n − i ⊗ A n i ⊗ I n + i (cid:17) = I n − i ⊗ ( τ ( A n i ) − A n i ) ⊗ I n + i . Invoking Theorem 4.4, for each i , it holds that12 < λ ( τ ( A n i ) − A n i ) = λ ( τ ( G ( α i ) n i ) − G ( α i ) n i ) < . Utilizing the properties of Kronecker product, it is easy to check each i holding12 < λ ( τ ( A i ) − A i ) < . This lemma shows that the spectrum of τ ( A i ) − A i is bounded for each i , whichwill be exploited to prove our aim conclusion that the spectrum of the preconditionedmatrix P − A is uniformly bounded. Theorem
The spectrum of the preconditioned matrix P − A is uniformlybounded below by / and bounded above by / .Proof. Let z ∈ R N be any nonzero vector. By the Rayleigh quotients theorem(see Theorem 4.2.2 in [19]) and Lemma 4.6, for each i , it holds12 < λ min ( τ ( A i ) − A i ) ≤ z (cid:124) A i zz (cid:124) τ ( A i ) z ≤ λ max ( τ ( A i ) − A i ) < , This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN that is, 12 z (cid:124) τ ( A i ) z < z (cid:124) A i z < z (cid:124) τ ( A i ) z . Thus we have 12 z (cid:124) m (cid:88) i =1 τ ( A i ) z < z (cid:124) m (cid:88) i =1 A i z < z (cid:124) m (cid:88) i =1 τ ( A i ) z . It results that 12 < z (cid:124) m (cid:80) i =1 A i zz (cid:124) m (cid:80) i =1 τ ( A i ) z < , i.e., 12 < z (cid:124) Azz (cid:124) Pz < . Therefore, we have λ min ( P − A ) = min z z (cid:124) Azz (cid:124) Pz > / , λ max ( P − A ) = max z z (cid:124) Azz (cid:124) Pz < / . Theorem 4.7 indicates that the spectrum of the preconditioned matrix is uniformlybounded. Furthermore, applying the results of Theorem 4.7, we obtain the followingcorollary.
Corollary
Let A be the coefficient matrix of the multi-dimensional Rieszfractional diffusion equation defined in (2.13), P be the preconditioner defined in (3.6) .Then the condition number of the preconditioned matrix P − A is less than . Up to now, for multi-dimensional Riesz fractional diffusion equations, we haveproved that the spectrum of the preconditioned matrix is uniformly bounded. Sincethe smallest eigenvalue of the preconditioned matrix is bounded away from 0, whichmanifests that the PCG method converges linearly. Moreover, we deduce that thecondition number is less than 3, which reveals the number of iterations is independentof the size of the coefficient matrix. From the theoretical point of view, we have provedthe efficiency of the τ -preconditioner for ill-conditioned linear systems (2.9) and (2.12)arising from Riesz fractional derivative, whose generating function is with fractionalorder zeros.
5. Extension to ill-conditioned multi-level Toeplitz matrices.
In this sec-tion, we extend our discussion to general multi-level Toeplitz matrices whose gener-ating functions are with fractional order zeros at the origin.Let α = ( α , . . . , α m ) with each α i ∈ (1 , θ = ( θ , . . . , θ m ) ∈ [ − π, π ] m . Suppose p α ( θ ) ∈ L ([ − π, π ] m ) is a nonnegative integrable even function with p α ( θ ) = 0 where θ = (0 , , . . . , p α ( θ ) generates an m -level Toeplitz matrix B , which has the same structure as (2.2), where the Toeplitz block at i -th level iswith the matrix size n i × n i for i = 1 , . . . , m . Let N = m (cid:81) i =1 n i . We derive an N × N This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE Bu = b , where u ∈ R N is unknown, b ∈ R N is the right hand side.As known that the inverse of the τ matrix can be obtained in O ( N log N ) com-putational cost by the fast discrete sine transform, we then construct a new τ -preconditioner based on the discretized Riesz derivatives that can match the fractionalorder zero of the generating function.Denote(5.2) q α ( θ ) = m (cid:88) i =1 l i | θ i | α i , where l i > i = 1 , . . . , m are constants. We refer to the multi-level Toeplitzmatrix arising from q α ( θ ) as Q . For convenience of our investigation, the generatingfunction p α ( θ ) is required to satisfy the following assumption. Assumption There exists two positive constants c < c such that < c ≤ p α ( θ ) q α ( θ ) ≤ c . Note that this assumption indicates that the generating function p α ( θ ) on eachdirection i has an α i -th order zero at θ i = 0 for i = 1 , . . . , m . According to preciousanalysis, we know that system (5.1) is symmetric positive definite and ill-conditioned.Therefore, a positive definite preconditioner is indispensable. However, the positivedefiniteness of B cannot guarantee that so is the corresponding τ -preconditioner τ ( B )[34]. Therefore, the τ ( B ) cannot be directly applied to precondition B . Actually, thenumerical results shown in Table 5 of Section 6 exhibit the bad performance of thepreconditioner τ ( B ). Therefore, it is essential to design a preconditioner aiming atsuch kind of ill-conditioned systems arising from generating functions with fractionalorder zeros. In the following, we construct a new preconditioner based on the τ -preconditioner of the Riesz fractional derivative that is symmetric positive definite.Denote(5.3) g α ( θ ) = m (cid:88) i =1 l i g α i ( θ i ) , where g α i ( θ i ) is defined as (2.11). Then, we obtain the corresponding matrix as G = m (cid:88) i =1 I n − i ⊗ l i G ( α i ) n i ⊗ I n + i , where G ( α i ) n i is defined as (2.10).Taking advantage of the above auxiliary tools, the specific procedures of con-structing preconditioner P of B are depicted as follows.1. The matrix Q generated by q α ( θ ) defined in (5.2) is employed to approximate B , denoted as P (1) = Q .2. The matrix G generated by (5.3) is exploited to approximate Q , denoted as P (2) = GQ − .3. The matrix τ ( G ) is used to approximate G , denoted as P (3) = τ ( G ) G − . This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
4. Finally, the preconditioner of B is constructed by P = P (3) P (2) P (1) ; i.e.,(5.4) P = τ ( G ) . To obtain the spectrum of the preconditioned matrix τ ( G ) − B , the propertiesof matrices P (1) , P (2) , P (3) will be considered. Under Assumption 1, the followinglemma is achieved. Lemma
Let B , Q be the multi-level Toeplitz matrices generated by p α ( θ ) and q α ( θ ) , respectively. Then, for any nonzero vector z ∈ R N , we have z (cid:124) Bzz (cid:124) Qz ∈ [ c , c ] . Proof.
Note that the eigenvalues of matrix Q − B are subjected to the Grenander–Szeg¨o’s theorem [17]; i.e.,min p α ( θ ) q α ( θ ) ≤ λ min ( Q − B ) ≤ λ max ( Q − B ) ≤ max p α ( θ ) q α ( θ ) . Combining the Rayleigh quotients theorem with Assumption 1, we obtain0 < c ≤ λ min ( Q − B ) ≤ z (cid:124) Bzz (cid:124) Qz ≤ λ max ( Q − B ) ≤ c . The proof is complete.Similarly, we have the following lemma.
Lemma
The matrix G is generated by g α ( θ ) defined in (5.3) . For α i ∈ (1 , , i = 1 , . . . , m , it holds that ≤ λ ( G − Q ) ≤ c , where c = max i π − πα i / .Proof. It is easy to check thatmin i | θ i | α i g α i ( θ i ) ≤ q α ( θ ) g α ( θ ) = m (cid:80) i =1 l i | θ i | α i m (cid:80) i =1 l i g α i ( θ i ) ≤ max i | θ i | α i g α i ( θ i ) . From Lemma 2.3, we derive that for each i holds12 ≤ | θ i | α i g α i ( θ i ) ≤ π − πα i / . Then, it immediately arrives12 ≤ q α ( θ ) g α ( θ ) ≤ max i π − πα i /
2) = c . In view of the proof of Lemma 5.1, it suffices to show that12 ≤ min g α ( θ ) q α ( θ ) ≤ λ ( G − Q ) ≤ max g α ( θ ) q α ( θ ) ≤ c . This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE
Theorem
The spectrum of preconditioned matrix τ ( G ) − B is bounded belowand above by constants independent of the matrix size; i.e., (5.5) c ≤ λ ( τ ( G ) − B ) ≤ c c . Proof.
Note that by Theorem 4.7, for any nonzero vector z ∈ R N , it follows that12 < z (cid:124) Gzz (cid:124) τ ( G ) z < . Then, combining Lemma 5.1 and Lemma 5.2, it holds that c ≤ z (cid:124) Bzz (cid:124) τ ( G ) z = z (cid:124) Gzz (cid:124) τ ( G ) z · z (cid:124) Qzz (cid:124) Gz · z (cid:124) Bzz (cid:124) Qz ≤ c c . Invoking the Rayleigh quotients theorem again, we derive c ≤ λ ( τ ( G ) − B ) ≤ c c . This theorem implies the spectrum of the preconditioned matrix is bounded,and the smallest eigenvalue of the preconditioned matrix is bounded away from 0,which means the linear convergent rate if the PCG method is exploited to solve ill-conditioned multi-level Toeplitz systems. Moreover, we derive the following corollary.
Corollary
The condition number of the preconditioned matrix τ ( G ) − B isbounded by constant.Proof. From Theorem 5.3, we obtain κ ( τ ( G ) − B ) = λ max ( τ ( G ) − B ) λ min ( τ ( G ) − B ) ≤ c c c . This corollary indicates the condition number of the preconditioned matrix isbounded by a constant. Hence, the number of iterations of the PCG method withour preconditioner for solving the system (5.1) is expected to be independent of thematrix size.
6. Numerical results.
In this section, the numerical experiments are carriedout to examine the efficiency of the proposed methods. All results are performed viaMATLAB R2017a on a PC with the configuration: Intel(R) Core(TM)i7-8700 [email protected] 3.20GHz and 8 GB RAM.To exhibit the performance of the PCG method with the τ -preconditioner forsolving linear systems (2.9) and (2.12), we also implement other preconditioners ascomparisons, including the Strang circulant preconditioner [20] and the banded pre- This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN conditioner [21], where the banded preconditioner B n is defined as B n = b b . . . b k − b b a ...... b ... b k − b k − ... ... ...... ... ... b b k − . . . b b + ∗ b k ... 2 ∗ n − (cid:80) j = k b j with ‘ k ’ being the bandwidth of the preconditioner. In practical computation, wechoose k = 8. Besides, the multigrid methods [13, 35] are also presented in ournumerical experiments.In the following tables, ‘ τ pre ’, ‘ C pre ’, ‘ B pre ’ and ‘ N pre ’ represent the PCG methodwith the τ -preconditioner, the Strang circulant preconditioner, the banded precondi-tioner, and no preconditioner, respectively. ‘ M pre ’ denotes the multigrid precondi-tioned method with the banded preconditioner proposed in [13] and ‘ M GM ’ repre-sents the algebraic multigrid method shown in [35]. Besides, ‘ n ’ denotes the spatialgrid points, ‘Iter’ displays the number of iterations required for convergence by thosemethods, and ‘CPU(s)’ signifies the CPU time in second for solving the linear sys-tems. In particular, ‘ − ’ implies the CPU time over 10 and ‘ ∗ ’ represents the numberof iterations over 10 . For all tested methods, let u = 0 be the initial guess and thestopping criterion is chosen as (cid:107) r q (cid:107) (cid:107) r (cid:107) < − , where r q is the residual vector at the q -th iteration. Example Consider 1D Riesz fractional diffusion equations defined in [0 , .Take the diffusion coefficient d = 1 , and the source term is determined by the ex-act solution, which is given by u ( x ) = x (1 − x ) . The number of iterations by those iterative methods for solving Example 1 aredisplayed in Table 1. It is obvious that both the τ -preconditioner and the circu-lant preconditioner show good performance, while the banded preconditioner needsmore iterations. On the other hand, Figure 1 exhibits the spectral distribution of thecoefficient matrix and the preconditioned matrices with α = 1 . n = 2 − x -axis direction represent the logarithmic values of eigen-values. We observe that the spectrum of the coefficient matrix A is scattered onthe coordinate axis from the left figure of Figure 1 and hence more iterations arerequired by the CG method to converge. The comparisons of the spectral distribu-tion of the preconditioned matrices with different preconditioners are shown in theright one, where ‘ τ − A ’, ‘ B − A ’, and ‘ C − A ’ denote the preconditioned matrices withthe τ -preconditioner, the banded preconditioner, and the circulant preconditioner, re-spectively. Note that the τ -preconditioner has a highly clustered spectrum, whichillustrates the better performance of the τ -preconditioner.Moreover, we list the extreme eigenvalues of τ − A for α = 1 . / , / This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE Table 1
Comparisons of the iterations for solving Example for different α by the CG method, and thePCG methods with the τ -preconditioner, the circulant preconditioner, and the banded preconditioner. n + 1 α = 1 . α = 1 . α = 1 . τ pre C pre B pre N pre τ pre C pre B pre N pre τ pre C pre B pre N pre Fig. 1 . The spectral distribution of coefficient matrix A and preconditioned matrices P − A . Table 2
Extreme eigenvalues of preconditioned matrix τ − A of Example with α = 1 . . n 2 λ max λ min Example In this example, we consider D Riesz fractional diffusion equationsin Ω ∈ [0 , × [0 , , where the exact solution is u ( x , x ) = x (1 − x ) x (1 − x ) . Take the diffusion coefficients d = d = 1 . The source term is given by y = d cos( πα / x (1 − x ) ( y ( x , α ) + y (1 − x , α ))+ d cos( πα / x (1 − x ) ( y ( x , α ) + y (1 − x , α )) , where y ( x, α ) = 2Γ(3 − α ) x − α − − α ) x − α + 24Γ(5 − α ) x − α . In the 2D case, we exploit the multigrid preconditioned method and the algebraic
This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
Table 3
Comparisons for solving Example for different α by the CG method, and the PCG methodswith the τ -preconditioner, and the circulant preconditioner, and multigrid methods. ( α , α ) n + 1 τ pre C pre M pre MGM N pre
Iter CPU(s) Iter CPU(s) Iter CPU(s) Iter CPU(s) Iter CPU(s)(1.1,1.2) 2 multigrid method for comparisons, where the weight of the Jacobi iterative methodis chosen as w = 2 / n = n . From Table 3, note that both of the numberof iterations and the CPU time provided by the τ -preconditioner are much less thanthose by other methods. It demonstrates the superiority of the τ -preconditioner. Example In this example, we test 3D Riesz fractional diffusion equations. Con-sider
Ω = [0 , × [0 , × [0 , , d = d = d = 1 ,y = d cos( πα / x (1 − x ) x (1 − x ) ( y ( x , α ) + y (1 − x , α ))+ d cos( πα / x (1 − x ) x (1 − x ) ( y ( x , α ) + y (1 − x , α ))+ d cos( πα / x (1 − x ) x (1 − x ) ( y ( x , α ) + y (1 − x , α )) . The exact solution is u = x (1 − x ) x (1 − x ) x (1 − x ) . In this example, let n = n = n . It is evident that the τ -preconditioner is stillefficient for 3D Riesz fractional diffusion equations from the Table 4. In contrast withthe circulant preconditioner, the number of iterations by the τ -preconditioner is almostunchanged when the matrix size increases. Hence, τ -preconditioner is an excellent toolfor solving ill-conditioned multi-level Toeplitz systems arising from multi-dimensionalRiesz fractional diffusion equations. This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE Table 4
Comparisons for solving Example for different α by the CG method, and the PCG methodswith the τ -preconditioner and the circulant preconditioner. ( α , α , α ) n + 1 τ pre C pre N pre Iter CPU(s) Iter CPU(s) Iter CPU(s)(1.1,1.2,1.3) 2 Finally, we verify the effectiveness of the preconditioning of discretized Rieszfractional derivatives for handling the ill-conditioned multi-level Toeplitz matriceswhose generating functions are with fractional order zeros at the origin.
Example Consider a two-level Toeplitz matrix whose generating function isdefined by p α ( θ ) = p α ( θ ) + p α ( θ ) − p ( θ ) p ( θ ) , where (see [8]) p α i ( θ i ) = (cid:26) | θ i | α i , | θ i | < π , , | θ i | ≥ π . It is obvious that 0 < − π ≤ p α ( θ ) q α ( θ ) ≤ , where q α ( θ ) is defined in (5.2) with l = l = 1. The above inequality exemplifiesthat the generating function p α ( θ ) satisfies the Assumption 1.To test the efficiency of our proposed preconditioner defined in (5.4), the PCGmethods with the τ -preconditioner, the circulant preconditioner, and the case without This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN
Table 5
Comparisons for solving Example for different α by the CG method, and the PCG methodswith the τ -preconditioner, the circulant preconditioner, and our preconditioner, and the algebraicmultigrid method. ( α , α ) n + 1 R pre τ pre C pre N pre MGM
Iter CPU(s) Iter CPU(s) Iter CPU(s) Iter CPU(s) Iter CPU(s)(1.9,1.5) 2
26 0.42 13 0.23 26 0.31 79 0.36 9 0.222
26 0.92 13 0.44 36 1.17 134 2.95 10 0.672
27 4.97 16 2.58 45 4.81 225 15.69 12 2.362
27 15.41 19 12.77 69 35.14 376 136.80 14 11.752
27 46.48 25 50.16 192 278.50 566 572.09 18 58.382
27 188.77 42 302.30 239 1.44e+3 956 4.19e+3 21 329.052
27 659.05 57 1.36e+3 782 – * – 26 1.71e+3(1.9,1.7) 2
26 0.53 14 0.14 29 0.63 90 0.42 9 0.222
26 0.84 17 0.63 44 1.34 159 4.44 9 0.772
26 3.95 20 2.98 79 8.02 278 20.34 10 2.242
26 15.47 30 19.05 145 73.41 467 132.67 11 9.772
27 47.89 49 83.78 270 385.14 846 907.67 12 48.602
27 188.72 85 565.31 736 5.20e+3 * – 13 216.622
27 690.84 219 5.09e+3 * – * – 15 1.01e+3(1.9,1.9) 2
27 0.39 15 0.19 30 0.27 97 0.64 9 0.332
27 0.91 21 0.93 57 1.77 179 3.56 9 0.642
27 4.17 26 4.12 85 8.58 327 22.75 9 2.332
27 17.27 44 26.34 212 111.33 573 210.27 9 8.412
27 47.45 80 135.02 622 884.95 943 948.84 9 37.302
27 190.06 195 1.29e+3 * – * – 10 166.472
27 660.47 530 – * – * – 10 683.06 preconditioner are proposed as comparisons. Denote our preconditioner as ‘ R pre ’ inTable 5.From Table 5, we see that the number of iterations deriving from the circulantpreconditioner, the natural τ -preconditioner, and no preconditioner increase rapidlyas the matrices size increase, while that by our proposed preconditioner almost keepsconstant. In terms of the multigrid method, we observe that the method imple-ments well when the orders of the zeros are closed. Since the cost per iteration ofusing the algebraic multigrid method is about 8 / τ -preconditioner [35], in light of Table 5, we derive that the performance of the multigridmethod is as efficient as the proposed method. Nevertheless, the multigrid methodwill be inefficient provided that the orders of the zeros are quite difference; see theitem for ( α , α ) = (1 . , .
5) in Table 5. Moreover, the linearly convergent rate is stilla question for these cases by the algebraic multigrid method.
7. Conclusion remarks.
In this paper, we have studied the spectra of the τ -preconditioned matrices for the multi-level Toeplitz systems arising from the multi-dimensional Riesz spatial fractional diffusion equations. Theoretically, we have provedthat the spectra of the preconditioned matrices are bounded below by 1/2 and boundedabove by 3/2, and hence the condition numbers of the preconditioned matrices are allless than 3. Besides, we proposed a new preconditioner for ill-conditioned multi-levelToeplitz systems, which are generated by the generating functions with fractionalorder zeros at the origin. We have proved that the spectra of the proposed precondi-tioned matrices are bounded by constants which are independent of the matrices size. This manuscript is for review purposes only.
PECTRAL ANALYSIS FOR PRECONDITIONING OF MULTI-DIMENSIONAL RFDE τ -preconditioner with other methods to handle more gen-eral (non-symmetric) multi-level Toeplitz-like systems arising from multi-dimensionalfractional partial differential equations. REFERENCES[1]
N. Barakitis, S. E. Ekstr¨om and P. Vassalos , Preconditioners for fractional diffusionequations based on the spectral symbol , ArXiv preprint arXiv: 1912.13304 (2019).[2]
F. Di Benedetto , Preconditioning of block Toeplitz matrices by sine transforms , SIAM J.Sci. Comput., 18 (1997), pp. 499–515.[3]
F. Di Benedetto, G. Fiorentino and S. Serra Capizzano , CG preconditioning for Toeplitzmatrices , Comput. Math. Appl., 25 (1993), pp. 35–45.[4]
D. Bini and F. Benedetto , A new preconditioner for the parallel solution of positive definiteToeplitz systems , in Proceedings, 2nd SPAA Conference, Crete, Greece, July 1990, pp.220–223.[5]
D. Bini and M. Capovani , Spectral and computational properties of band symmetric Toeplitzmatrices , Linear Algebra Appl., 52/53 (1983), pp. 99–126.[6]
A. B¨ottcher and S. Grudsky , On the condition numbers of large semi-definite Toeplitzmatrices , Linear Algebra Appl., 279 (1998), pp. 285–301.[7]
R. H. Chan , Toeplitz preconditioner for Toeplitz system with nonnegative generating function ,IMA J. Numer. Anal., 11 (1991), pp. 333–345.[8]
R. H. Chan, Q. S. Chang and H. W. Sun , Multigrid method for ill-conditioned symmetricToeplitz systems , SIAM J. Sci. Comput., 19 (1998), pp. 516–529.[9]
R. H. Chan, M. K. Ng and C. K. Wong , Sine transform based preconditioners for symmetricToeplitz systems , Linear Algebra Appl., 232 (1996), pp. 237–259.[10]
R. H. Chan, A. M. Yip and M. K. Ng , The best circulant preconditioners for HermitianToeplitz systems , SIAM J. Numer. Anal., 38 (2000), pp. 876–896.[11]
W. K. Ching , Iterative methods for Queuing and Manufacturing systems , Springer-Verlag,London, 2001.[12]
H. F. Ding and Y. X. Zhang , New numerical methods for the Riesz space fractional partialdifferential equations , Comput. Math. Appl., 63 (2012), pp. 1135–1146.[13]
M. Donatelli, R. Krause, M. Mazza and K. Trotti , Multigrid preconditioners for aniso-tropic space-fractional diffusion equations , Adv Comput Math., 49 (2020).[14]
M. Donatelli, M. Mazza and S. Serra Capizzano , Spectral analysis and structure pre-serving preconditioners for fractional diffusion equations , J. Comput. Phys., 307 (2016),pp. 262–279.[15]
G. Fiorentino and S. Serra Capizzano , Multigrid methods for Toeplitz matrices , Calcolo.,28 (1991), pp. 283–305.[16]
G. Fiorentino and S. Serra Capizzano , Multigrid methods for symmetric positive definiteblock Toeplitz matrices with nonnegative generating functions , SIAM J. Sci. Comput., 17(1996), pp. 1068–1081.[17]
U. Grenander and G. Szeg¨o , Toeplitz Forms and Their Applications , 2nd ed., Chelsea, NewYork, 1984.[18]
R. Hilfer , Applications of Fractional Calculus in Physics , World Scientific, Singapore, 2000.[19]
R. A. Horn and C. R. Johnson , Matrix Analysis , Cambridge University Press, Cambridge,2012.[20]
S. L. Lei and H. W. Sun , A circulant preconditioner for fractional diffusion equations , J.Comput. Phys., 242 (2013, pp. 715–725.[21]
F. R. Lin, S. W. Yang and X. Q. Jin , Preconditioned iterative methods for fractional diffu-sion equation , J. Comput. Phys., 256 (2014), pp. 109–117.[22]
M. M. Meerschaert and C. Tadjeran , Finite difference approximations for fractionaladvection-dispersion flows equations , J. Comput. Appl. Math., 172 (2004), pp. 209–219.[23]
H. Moghaderi, M. Dehghan, M. Donatelli and M. Mazza , spectral analysis and multi-grid preconditioners for two-dimensional space-fractional diffusion equations , J. Comput.Phys., 350 (2017), pp. 992–1011.[24] D. Noutsos, S. Serra Capizzano and P. Vassalos , Matrix algebra preconditioners formultilevel Toeplitz systems do not insure optimal convergence rate , Theor Comput Sci.,315 (2004), pp. 557–579.
This manuscript is for review purposes only. XIN HUANG, XUE-LEI LIN, MICHAEL K. NG, AND HAI-WEI SUN[25]
D. Noutsos, S. Serra Capizzano and P. Vassalos , Essential spectral equivalence via mul-tiple step preconditioning and applications to ill conditioned Toeplitz matrices , LinearAlgebra Appl., 491 (2016), pp. 276–291.[26]
H. K. Pang and H. W. Sun , Fast numerical contour integral method for fractional diffusionequations , J. Sci. Comput., 66 (2016), pp. 41–66.[27]
J. Pestana , Preconditioners for symmetrized Toeplitz and multi-level Toeplitz matrices ,SIAM J. Matrix. Anal. Appl., 44 (2019), pp. 870–887.[28]
I. Podlubny , Fractional differential equations , Academic Press, New York, 1999.[29]
D. Potts and G. Steidl , Preconditioners for ill-conditioned Toeplitz Matrices , BIT., 39(1999), pp. 513–533.[30]
D. Potts and G. Steidl , Preconditioners for ill-conditioned Toeplitz systems constructedfrom positive kernels , SIAM J. Sci. Comput., 22 (2001), pp. 1741–1761.[31]
A. I. Saichev and G. M. Zaslavsky , Fractional kinetic equations: solutions and applications ,Chaos., 7 (1997), pp. 753–764.[32]
S. Serra Capizzano , New PCG based algorithms for the solution of Hermitian Toeplitzsystems , Calcolo., 32 (1995), pp. 153–176.[33]
S. Serra Capizzano , On the extreme eigenvalues of Hermitian (block) Toeplitz matrices ,Linear Algebra Appl., 270 (1998), pp. 109–129.[34]
S. Serra Capizzano , Superlinear PCG methods for symmetric Toeplitz systems , Math. Com-put., 68 (1999), pp. 793–803.[35]
H. W. Sun, X. Q. Jin and Q. S. Chang , Convergence of the multigrid method for ill-conditioned block Toeplitz systems , BIT., 41 (2001), pp. 179–190.[36]
V. E. Tarasov , Fractional Dynamics: Applications of Fractional Calculus to Dynamics ofParticles , Fields and Media Higher Education Press, Beijing, 2010., Fields and Media Higher Education Press, Beijing, 2010.