Fully Finite Element Adaptive Algebraic Multigrid Method for Time-Space Caputo-Riesz Fractional Diffusion Equations
Xiaoqiang Yue Weiping Bu Shi Shu Menghuan Liu Shuai Wang
math.NA

0in 6.0in 8.0in










[5]

Journal of Computational and Applied Mathematics






























elsarticle-num


numbers,sort   compress

TheoremTheorem
ExampleExample
LemmaLemma
RemarkRemark
CorollaryCorollary
DefinitionDefinition
AlgorithmAlgorithm
AssumptionAssumption
Proof[1][Proof.]

  #1








mymainaddress]Xiaoqiang Yue

mymainaddress]Weiping Bu

mysecondaryaddress]Shi Shumycorrespondingauthor
[mycorrespondingauthor]Corresponding author
shushi@xtu.edu.cn

mymainaddress]Menghuan Liu

mymainaddress]Shuai Wang

[mymainaddress]School of Mathematics and Computational Science, Xiangtan University, Hunan 411105, P.R. China
[mysecondaryaddress]Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Hunan 411105, P.R. China


The paper aims to establish a fully discrete finite element (FE) scheme and provide cost-effective solutions for 
one-dimensional time-space Caputo-Riesz fractional diffusion equations on a bounded domain Ω.
Firstly, we construct a fully discrete scheme of the linear FE method in both temporal and spatial directions,
derive many characterizations on the coefficient matrix
and numerically verify that the fully FE approximation possesses the saturation error order under L^2(Ω) norm.
Secondly, we theoretically prove the estimation 1+𝒪(τ^α h^-2β) on the condition number of the coefficient matrix,
in which τ and h respectively denote time and space step sizes. Finally, on the grounds of the estimation and fast Fourier transform,
we develop and analyze an adaptive algebraic multigrid (AMG) method with low algorithmic complexity,
reveal a reference formula to measure the strength-of-connection tolerance
which severely affect the robustness of AMG methods in handling fractional diffusion equations,
and illustrate the well robustness and high efficiency of the proposed algorithm compared with the classical AMG,
conjugate gradient and Jacobi iterative methods.



Caputo-Riesz fractional diffusion equation, fully time-space FE scheme, condition number estimation, algorithmic complexity, adaptive AMG method
[2010] 35R1165F1065F1565N55








§ INTRODUCTION


In recent years, there has been an explosion of research interest in numerical solutions for fractional differential equations,
mainly due to the following two aspects: (i) the huge majority can't be solved analytically,
(ii) the analytical solution (if luckily derived)
always involve certain infinite series which sharply drives up the costs of its evaluation.
Various numerical methods have been proposed to approximate
more accurately and faster, such as finite difference (FD) method <cit.>,
finite element (FE) method <cit.>,
finite volume <cit.> method and spectral (element) method <cit.>.
An essential challenge against standard differential equations
lies in the presence of the fractional differential operator,
which gives rise to nonlocality (space fractional, nearly dense or full coefficient matrix) or
memory-requirement (time fractional, the entire time history of evaluations) issue, resulting in a vast computational cost.

Preconditioned Krylov subspace methods are regarded as one of the potential solutions to the aforementioned challenge.
Numerous preconditioners with various Krylov-subspace methods have been constructed respectively for one- and two-dimensional,
linear and nonlinear space-fractional diffusion equations (SFDE) <cit.>.
Multigrid method has been proven to be a superior solver
and preconditioner for ill-conditioned Toeplitz systems as well as SFDE.
Pang and Sun propose an efficient and robust geometric multigrid (GMG)
with fast Fourier transform (FFT) for one-dimensional SFDE by an implicit FD scheme <cit.>.
Bu et al. employ the GMG to one-dimensional multi-term time-fractional advection-diffusion equations
via a fully discrete scheme by FD method in temporal and FE method in spatial directions <cit.>.
Jiang and Xu construct optimal GMG for two-dimensional SFDE to get FE approximations <cit.>.
Chen et al. make the first attempt to present an algebraic multigrid (AMG) method with line smoothers to
the fractional Laplacian through localizing it into a nonuniform elliptic equation <cit.>.
Zhao et al. invoke GMG for one-dimensional Riesz SFDE by an adaptive FE scheme using hierarchical matrices <cit.>.
From the survey of references, in spite of quite a number of contributions to numerical methods and preconditioners,
there are no calculations taking into account of fully discrete FE schemes and AMG methods
for time-space Caputo-Riesz fractional diffusion equations.

In this paper, we are concerned with the following time-space Caputo-Riesz fractional diffusion equation (CR-FDE)

    _0^CD^α_tu(x,t)=∂^2βu(x,t)/∂|x|^2β+f(x,t),  
      t∈ I=(0,T], x∈Ω=(a,b) 
       u(x,t)=0,  t∈ I, x∈∂Ω
       u(x,0)=ψ_0(x),  x∈Ω

with orders α∈(0,1) and β∈(1/2,1), the Caputo and Riesz fractional derivatives are respectively defined by

    _0^CD^α_tu=1/Γ(1-α)∫_0^t(t-s)^-α∂ u/∂ sds, ∂^2βu/∂|x|^2β=-1/2cos(βπ)(_xD_L^2βu
      + _xD_R^2βu),

where

    _xD_L^2βu = 1/Γ(2-2β)∂^2/∂ x^2∫_a^x(x-s)^1-2βuds, _xD_R^2βu = 1/Γ(2-2β)∂^2/∂ x^2∫_x^b(s-x)^1-2βuds.


The remainder of this paper proceeds as follows. A fully discrete FE method of (<ref>)-(<ref>) is developed in Section 2.
Section 3 comes up with the theoretical estimation and verification experiments on the condition number of the coefficient matrix.
The classical AMG method is introduced in Section 4 followed by its uniform convergence analysis
and the construction of an adaptive AMG method. Section 5 reports and analyzes numerical results to show the benefits.
We close in Section 6 with some concluding remarks.



§ FULLY DISCRETE FINITE ELEMENT SCHEME FOR THE CR-FDE


For simplicity, following <cit.>, we will use the symbols ≲, ≳ and ≃ throughout the paper.
u_1≲ v_1 means u_1≤ C_1v_1, u_2≳ v_2 means u_2≥ c_2v_2 while u_3≃ v_3 means
c_3v_3≤ u_3≤ C_3v_3, where C_1, c_2, c_3 and C_3 are generic positive constants independent of variables, time and space step sizes.



 §.§ Reminder about fractional calculus


In this subsection, we briefly introduce some fractional derivative spaces and several auxiliary results.
Here the L^2 inner product and norm are denoted by

    (u,v)_L^2(Ω)=∫_Ωuvdx, u_L^2(Ω)=(u,u)^1/2_L^2(Ω).


 (Left and right fractional derivative spaces)
  For constant μ>0, define norms

    u_J_L^μ(Ω) := (u^2_L^2(Ω)+_xD_L^μ u^2_L^2(Ω))^1/2, u_J_R^μ(Ω) := (u^2_L^2(Ω)+_xD_R^μ u^2_L^2(Ω))^1/2,

  and let J_L,0^μ(Ω) and J_R,0^μ(Ω) be closures of C_0^∞(Ω)
  under ·_J_L^μ(Ω) and ·_J_R^μ(Ω), respectively.


 (Fractional Sobolev space)
  For constant μ>0, define the norm

    u_H^μ(Ω) := (u^2_L^2(Ω)+|ξ|^μũ^2_L^2(Ω_ξ))^1/2,

  and let H_0^μ(Ω) be the closure of C_0^∞(Ω) under ·_H^μ(Ω),
  where ũ is the Fourier transform of u.



  The equivalence between (<ref>) and the general definition of the norm has been established in <cit.>,
  which implies the reasonability of Definition <ref>.


 (see <cit.>, Proposition 1)
  If constant μ∈(0,1), u,v∈ J_L,0^2μ(Ω) (or J_R,0^2μ(Ω)), then

    (_xD_L^2μu,v)_L^2(Ω)=(_xD_L^μ u,_xD_R^μ v)_L^2(Ω), 
      (_xD_R^2μu,v)_L^2(Ω)=(_xD_R^μ u,_xD_L^μ v)_L^2(Ω).



 (see <cit.>, Lemma 2.4)
  For constant μ > 0, we have

    (_xD_L^μ u, _xD_R^μ u)_L^2(Ω) = cos(πμ)_xD_L^μ u^2_L^2(Ω).



 (Fractional Poincaré-Friedrichs inequality, see <cit.>, Theorem 2.10)
  For u∈ J_L,0^μ(Ω), we have

    u_L^2(Ω)≲_xD_L^μ u_L^2(Ω).





 §.§ Derivation of the fully discrete scheme


By Lemma <ref>, we get the variational (weak) formulation of (<ref>)-(<ref>):
given f∈ L^2(Ω,I), ϕ_0∈ L^2(Ω) and Q_t:=Ω×(0,t), to find u∈ℋ subject to u(x,0)=ψ_0(x) and

    (_0^CD^α_σ u,v)_Q_t+B^t_Ω(u,v)=(f,v)_Q_t, ∀ v∈ℋ^*,

where ℋ:=H_0^β(Ω)× H^1(I), ℋ^*:=H_0^β(Ω)× L^2(I), and

    (_0^CD^α_σ u,v)_Q_t=∫_0^t(_0^CD^α_σ u,v)_L^2(Ω)dσ,
       (f,v)_Q_t=∫_0^t(f,v)_L^2(Ω)dσ, 
       B^t_Ω(u,v)=∫_0^t1/2cos(βπ)[(_xD^β_Lu,_xD^β_Rv)_L^2(Ω)
      +(_xD^β_Ru,_xD^β_Lv)_L^2(Ω)]dσ.





In order to acquire numerical solutions of u, we firstly make a (possibly nonuniform) temporal discretization by points 0=t_0<t_1<⋯<t_N=T,
and a uniform spatial discretization by points x_i = a+ih (i=0,1,⋯,M), where h=(b-a)/M represents the space step size. Let

    I_j=(t_j-1,t_j), Ĩ_j=(0,t_j), j=1,2,⋯,N; Ω_h = {Ω_l:Ω_l=(x_l-1,x_l), l=1,2,⋯,M}.


We observe that it is convenient to form the FE spaces in tensor products

    𝒱_n=𝒱^β_h(Ω_h)×𝒱_τ(Ĩ_n), 𝒱^*_n=𝒱^β_h(Ω_h)×𝒱^*_τ(I_n),

where

    𝒱^β_h(Ω_h)={w_h∈ H_0^β(Ω)∩ C(Ω̅): 
      w_h(x)|_Ω_l∈𝒫_1(Ω_l), l=1,⋯,M},
       𝒱_τ(Ĩ_n) = {v_τ∈𝒞(Ĩ_n): v_τ(0)=1, 
      v_τ(t)|_I_j∈𝒫_1(I_j), j=1,⋯, n},
       𝒱^*_τ(I_n) = {v_τ∈ L^2(I_n): 
      v_τ(t)|_I_n∈𝒫_0(I_n)},

and 𝒫_k denotes the set of all polynomials of degree ≤ k.


  Apparently, for a given u_hτ(x,t)∈𝒱_n, we have ∂ u_hτ/∂ t∈𝒱^*_n,
  where ∂ u_hτ/∂ t is obtained by differentiating u_hτ
  with respect to t on each subinterval I_j (j=1,2,⋯,N).


We obtain a fully discrete FE scheme in temporal and spatial directions of problem (<ref>):
given Q_n:=Ω_h× I_n, to find u_hτ∈𝒱_n such that u_hτ(x,0)=ψ_0,I(x) and

    (_0^CD^α_t u_hτ,v_hτ)_Q_n+B^n_Ω(u_hτ,v_hτ)=
      (f,v_hτ)_Q_n, ∀ v_hτ∈𝒱^*_n,

where ψ_0,I(x)∈𝒱_n satisfying ψ_0,I(x_i)=ψ_0(x_i) (i=0,1,⋯,M), and

    (_0^CD^α_tu_hτ,v_hτ)_Q_n=∫_t_n-1^t_n (_0^CD^α_t u_hτ, v_hτ)_L^2(Ω)dt, 
      (f,v_hτ)_Q_n=∫_t_n-1^t_n(f,v_hτ)_L^2(Ω)dt,
        B^n_Ω(u_hτ,v_hτ)
      = ∫_t_n-1^t_n1/2cos(βπ)[(_xD^β_Lu_hτ,_xD^β_Rv_hτ)_L^2(Ω)
      +(_xD^β_Ru_hτ,_xD^β_Lv_hτ)_L^2(Ω)]dt.


Let

    ℒ_0(t)={   t_1-t/τ_1, t ∈ I_1 
        0,     t∈Ĩ_n∖ I_1
        ., ℒ̃_0(t)=1/Γ(1-α)∫_t_0^tdℒ_0(s)/(t-s)^α, ℒ̂_0(t)=1/Γ(1-α)∫_t_0^t_1dℒ_0(s)/(t-s)^α,


    ℒ_k(t)={   t_k+1-t/τ_k+1, t ∈ I_k+1
       t-t_k-1/τ_k, t ∈ I_k 
        0,        t∈Ĩ_n∖ (I_k ∪ I_k+1)
        ., ℒ̂_k(t)=1/Γ(1-α)∫_t_k-1^t_k+1dℒ_k(s)/(t-s)^α, 
      k=1,⋯,n-1

and

    ℒ_n(t)={   t-t_n-1/τ_n, t ∈ I_n 
        0,        t∈Ĩ_n∖ I_n
        ., ℒ̃_n(t)=1/Γ(1-α)∫_t_n-1^tdℒ_n(s)/(t-s)^α.

Note that

    𝒱^*_n={ϕ_l(x)× 1, l=1,⋯,M-1},

where ϕ_l(x) is the shape function at x_l∈Ω_h. Using

    u_hτ(x,t)=u_h^0(x)ℒ_0(t)+∑_k=1^n-1u_h^k(x)ℒ_k(t)+u_h^n(x)ℒ_n(t),

we have

    (_0^CD^α_tu_hτ,ϕ_l× 1)_Q_1
      =(u_h^0,ϕ_l)_L^2(Ω)(ℒ̃_0,1)_L^2(I_1)
      +(u_h^1,ϕ_l)_L^2(Ω)(ℒ̃_1,1)_L^2(I_1), 
       (_0^CD^α_tu_hτ,ϕ_l× 1)_Q_n
      =(u_h^0,ϕ_l)_L^2(Ω)(ℒ̂_0,1)_L^2(I_n)+ 
           ∑_k=1^n-1(u_h^k,ϕ_l)_L^2(Ω)(ℒ̂_k,1)_L^2(I_n)
      +(u_h^n,ϕ_l)_L^2(Ω)(ℒ̃_n,1)_L^2(I_n), n>1, 
       ∫_t_n-1^t_n1×(_xD^β_Lu_hτ,_xD^β_Rϕ_l)_L^2(Ω)dt
      =(_xD^β_Lu_h^0,_xD^β_Rϕ_l)_L^2(Ω)
      (ℒ_0,1)_L^2(I_n)+
           ∑_k=1^n-1(_xD^β_Lu_h^k,_xD^β_Rϕ_l)_L^2(Ω)
      (ℒ_k,1)_L^2(I_n) +
      (_xD^β_Lu_h^n,_xD^β_Rϕ_l)_L^2(Ω)
      (ℒ_n,1)_L^2(I_n), 
       ∫_t_n-1^t_n1×
      (_xD^β_Ru_hτ,_xD^β_Lϕ_l)_L^2(Ω)dt
      =(_xD^β_Ru_h^0,_xD^β_Lϕ_l)_L^2(Ω)
      (ℒ_0,1)_L^2(I_n)+ 
           ∑_k=1^n-1(_xD^β_Ru_h^k,_xD^β_Lϕ_l)_L^2(Ω)
      (ℒ_k,1)_L^2(I_n)+
      (_xD^β_Ru_h^n,_xD^β_Lϕ_l)_L^2(Ω)
      (ℒ_n,1)_L^2(I_n).


Substituting (<ref>)-(<ref>) into (<ref>), yields

    C^n_hτU^n_hτ = G^n_hτ,

where the coefficient matrix

    C^n_hτ = M_h+Γ(3-α)/2τ^α_n A_h^β,

the right-hand side vector

    G^n_hτ=Γ(3-α)τ_n^α-1F^n_hτ+
      [M_h-Γ(3-α)/2τ^α_n A_h^β]U^n-1_hτ - ∑_k=1^n-1τ_n^α-1×
         (t_n-t_k-1)^2-α-(t_n-1-t_k-1)^2-α -(t_n-t_k)^2-α+(t_n-1-t_k)^2-α/τ_k
      M_h (U^k_hτ - U^k-1_hτ),

the mass matrix

    M_h = h/6(
       [ 4 1      ; 1 4 1    ;   ⋱ ⋱ ⋱  ;     1 4 1;       1 4 ])_(M-1)×(M-1),

the stiffness matrix A_h^β = (a^h_i,j)_(M-1)×(M-1) with its entries

    {    a^h_i,i = h^1-2β(2^4-2β-8)/2cos(βπ)Γ(4-2β),                                  i=1,⋯,M-1
        a^h_j,j+1=a^h_j+1,j=h^1-2β(3^3-2β-2^5-2β+7)/2cos(βπ)Γ(4-2β),                j=1,⋯,M-2
        a^h_k,k+l= a^h_k+l,k = h^1-2β/2cos(βπ)Γ(4-2β)[(l+2)^3-2β
             -4(l+1)^3-2β+6l^3-2β-4(l-1)^3-2β+(l-2)^3-2β], k=1,⋯,M-l-1
        .,

the vector 

    F^n_hτ=(f^n_1,f^n_2,⋯,f^n_M-1)^T, f^n_l=(f,ϕ_l× 1)_Q_n, l=1,⋯,M-1

and the fully FE approximations

    U^k_hτ=(u^k_1,u^k_2,⋯,u^k_M-1)^T, u^0_j=ψ_0,I(x_j), 
      u^k_j=u_h^k(x_j), k=1,⋯,n, j=1,⋯,M-1.



  (<ref>) is reduced via dividing both sides of (<ref>) by the factor
  τ^1-α_n / Γ(3-α), in case of the severe loss in convergence of the fully discrete FE scheme.


Next, a number of 
characterizations are established regarding A_h^β just defined by (<ref>).


  The stiffness matrix A_h^β is symmetric and satisfies

  
  * a^h_i,i > 0 for i=1,⋯,M-1;
  
  * a^h_i,j < 0 for i j, i,j=1,⋯,M-1;
  
  * ∑_j=1^M-1 a^h_i,j>0 for i=1,⋯,M-1;
  
  * The following relation holds for the particular case when h≤1/7

    ∑_j=1^M-1 a^h_i,j≥{    -h^1-2β(4-2^3-2β)/2cos(βπ)Γ(4-2β), i=1,M-1 
        -2^2βh(2β-1)/cos(βπ)Γ(2-2β),   i=2,⋯,M-2
        .;

  
  * A_h^β is an M-matrix.




  The symmetric property of A_h^β is an obvious fact by (<ref>).
  Since β∈(1/2,1), then 4-2β<3 and cos(βπ)<0, which give immediately a^h_ii > 0, i=1,⋯,M-1.
  This proves the first part of the theorem.
  The second part is an immediate consequence of the facts that on the interval β∈(1/2,1),
  f(β)=3^3-2β-2^5-2β+7 is a strictly increasing function, and the bivariate function

    f_β(l)=(l+2)^3-2β-4(l+1)^3-2β+6l^3-2β-4(l-1)^3-2β+(l-2)^3-2β > 0, 2≤ l≤ M-2.

  In fact, it is evident that

    1.5^2β > 1.5 > 3^3/2^5ln3/ln2⇒ f'(β)=-2ln3 · 3^3-2β + 2ln2 · 2^5-2β > 0,

  and

    f_β(l)=h^2β-3[g(x_l+2)-4g(x_l+1)+6g(x_l)-4g(x_l-1)+g(x_l-2)] > 0

  using Taylor's expansion with

    (l/l+1)^2+2β-(l/l-1)^2+2β > -30l/2β+1⇒
      l^-1-2β+2β+1/30[(l+1)^-2-2β-(l-1)^-2-2β] > 0,

  where g(x)=(x-a)^3-2β and x_l=lh+a.

  To prove the third part, use

    (_xD^β_Lϕ_i,_xD^β_Rϕ_j)_L^2(Ω) = -(_xD_L^2β-1ϕ_i,dϕ_j/dx)_L^2(Ω)

  and

    ϕ̃:=∑_j=1^M-1ϕ_j=1-ϕ_0-ϕ_M

  to obtain the relation

    ∑_j=1^M-1 a^h_i,j=-(_xD_L^2β-1ϕ̃,1)_Ω_i-(_xD_L^2β-1ϕ̃,1)_Ω_i+1
      +(_xD_L^2β-1ϕ_i,1)_Ω_1-(_xD_L^2β-1ϕ_i,1)_Ω_M_h/2hcos(βπ),

  where

    _xD^2β-1_Lϕ̃={   (x-a)^2-2β/hΓ(3-2β),                                a < x < x_1 
       (x-a)^2-2β-(x-x_1)^2-2β/hΓ(3-2β),                   x_1 < x < x_M_h-1
       (x-a)^2-2β-(x-x_1)^2-2β-(x-x_M_h-1)^2-2β/hΓ(3-2β), x_M_h-1 < x < x_M_h.

  and

    _xD_L^2β-1ϕ_i(x)={    0,                                           x < x_i-1
       (x-x_i-1)^2-2β/hΓ(3-2β),                               x_i-1 < x < x_i 
       (x-x_i-1)^2-2β-2(x-x_i)^2-2β/hΓ(3-2β),                 x_i < x < x_i+1
       (x-x_i-1)^2-2β-2(x-x_i)^2-2β+(x-x_i+1)^2-2β/hΓ(3-2β), x > x_i+1..

  Assume that Ω=(0,1) without loss of generality, one can easily derive

    ∑_j=1^M-1 a^h_i,j=-(4-2^3-2β)h^3-2β-1+3(1-h)^3-2β
      -3(1-2h)^3-2β+(1-3h)^3-2β/2cos(βπ)h^2Γ(4-2β), i=1,M-1,
       ∑_j=1^M-1 a^h_i,j=-3(ih)^3-2β-3[(i-1)h]^3-2β
      +[(i-2)h]^3-2β-[(i+1)h]^3-2β/2cos(βπ)h^2Γ(4-2β)
        - 3(1-ih)^3-2β-[1-(i-1)h]^3-2β-3[1-(i+1)h]^3-2β
      +[1-(i+2)h]^3-2β/2cos(βπ)h^2Γ(4-2β), i=2,⋯,M-2

  and deduce ∑_j=1^M-1 a^h_i,j>0 by Taylor's formula and β∈(1/2,1).

  Another step to do in the proof is the result 4, which follows from

    h≤1/7⇒7-4(1-ξ)^-1-2β > -1/β h⇒-1+3(1-h)^3-2β-3(1-2h)^3-2β+(1-3h)^3-2β>0

  for all ξ∈(0,2h),

    (i-1/i)^3+2β > β+1/21(β+1)+30i⇒3(ih)^3-2β-3(ih-h)^3-2β+(ih-2h)^3-2β-
      (ih+h)^3-2β/h^3(3-2β)(2-2β)(2β-1)(ih)^-2β>1

  and

    3(1-ih)^3-2β-[1-(i-1)h]^3-2β-3[1-(i+1)h]^3-2β
      +[1-(i+2)h]^3-2β/h^3(3-2β)(2-2β)(2β-1)(1-ih)^-2β>1

  for i=2,⋯,M-2, together with the inequality (ih)^-2β+(1-ih)^-2β≥ 2^1+2β.

  Finally, according to properties 1 and 2, the result 5 will be proved by showing that
  (A_h^β)^-1 is nonnegative, which can be easily proved by contradiction with property 3.


Observe from (<ref>)-(<ref>) that M_h and A_h^β
are both symmetric Toeplitz matrices independent of any time terms.
The under-mentioned corollaries are natural consequences of Theorem <ref>.


  The coefficient matrix C^n_hτ is a symmetric Toeplitz matrix. Furthermore,
  it will be independent of time level n if the temporal discretization is also uniform.



  The coefficient matrix C^n_hτ is an M-matrix, if and only if

    τ_n^α/h^2β > -2cos(βπ)Γ(4-2β)/3Γ(3-α)(3^3-2β-2^5-2β+7).




  This result will follow from Theorem <ref>, if we can show that

    h/6+Γ(3-α)/2τ^α_n
      h^1-2β(3^3-2β-2^5-2β+7)/2cos(βπ)Γ(4-2β)<0,

  which is an immediate application of the condition (<ref>).




 §.§ Numerical experiments and the saturation error order



  Consider (<ref>)-(<ref>) with Ω=(0,1), T=1, ψ_0(x)=0 and

    f(x,t) =Γ(3-α)/Γ(3-2α)t^2-2αx^2(1-x)^2 +
      t^2-α/cos(βπ)[x^2-2β+(1-x)^2-2β/Γ(3-2β)-
                              6x^3-2β+6(1-x)^3-2β/Γ(4-2β) +
      12x^4-2β+12(1-x)^4-2β/Γ(5-2β)].



The exact solution is u(x,t)=t^2-αx^2(1-x)^2. In the case of uniform temporal and spatial meshes,
Tables <ref> and <ref>
present errors e_0 := u(·,1)-u_hτ(·,1)_L^2(Ω) and convergence rates.




















From Tables <ref> and <ref>, we can obtain that
the fully FE solution u_hτ achieves the saturation error order 𝒪(τ^2+h^2) under ·_0 norm.

Fig. <ref> illustrates the comparisons of exact solutions and numerical solutions
of α=0.2, 0.4 and β=0.6, 0.8 with t=1 and h=τ=1/32.





§ CONDITION NUMBER ESTIMATION


This section is devoted to deriving the condition number estimation on the coefficient matrix of (<ref>)
in uniform temporal and and spatial discretizations.


  For the linear system (<ref>), we have

    κ(C^n_hτ)= 1+𝒪(τ^α h^-2β).




  Let C_α = Γ(3-α)/2, we divide our proof in three steps.
  First, it is trivially true that C^n_hτ
  is spectrally equivalent to the matrix I+C_ατ^α M_h^-1/2A_h^β M_h^-1/2, i.e.

    κ(C^n_hτ)≃κ(I+C_ατ^α M_h^-1/2A_h^β M_h^-1/2).


  The next thing to do in the proof is to verify

    λ_min(M_h^-1/2A_h^β M_h^-1/2)≳ 1, λ_max(M_h^-1/2A_h^β M_h^-1/2)≲ h^-2β,

  which is equivalent to

    (v⃗_h,v⃗_h)≲ (M_h^-1/2A_h^β M_h^-1/2v⃗_h,v⃗_h)
      ≲ h^-2β(v⃗_h,v⃗_h), ∀v⃗_h∈ℝ^M-1.

  Set u⃗_h = M_h^-1/2v⃗_h := (u^h_1,⋯,u^h_M-1)^T, rewrite (<ref>) as
  (M_hu⃗_h,u⃗_h)≲ (A_h^βu⃗_h,u⃗_h)≲ h^-2β(M_hu⃗_h,u⃗_h).
  It is sufficient to verify that
  (M_hu⃗_h,u⃗_h)≃ h(u⃗_h,u⃗_h).
  It follows by (<ref>) and the Cauchy-Schwarz inequality that

    h/3(u⃗_h, u⃗_h)≤(M_h u⃗_h, u⃗_h)=h[2/3∑_l=1^M-1(u^h_l)^2
      +1/3∑_l=1^M-2u^h_lu^h_l+1]≤ h(u⃗_h, u⃗_h).

  Thus (<ref>) will follow if we can show that
  h(u⃗_h,u⃗_h)≲ (A_h^βu⃗_h,u⃗_h)≲ h^1-2β(u⃗_h,u⃗_h).
  We start by showing the second inequality. Utilizing Theorem <ref> and the Cauchy-Schwarz inequality, we arrive at

    (A_h^βu⃗_h, u⃗_h)    ≤   ∑_i=1^M-1 a^h_i,i (u^h_i)^2 -
     1/2∑_i=1^M-1∑_j i a^h_i,j [(u^h_i)^2 + (u^h_j)^2] 
       =   ∑_i=1^M-1 a^h_i,i (u^h_i)^2 - ∑_i=1^M-1∑_j=i+1^M-1 a^h_i,j (u^h_i)^2
      - ∑_i=1^M-1∑_j=i+1^M-1 a^h_i,j (u^h_j)^2 
       =   ∑_i=1^M-1 (u^h_i)^2 [a^h_i,i - ∑_j=i+1^M-1 a^h_i,j - ∑_j=1^i-1 a^h_i,j]
       ≤    2 a^h_1,1 (u⃗_h, u⃗_h) = (2^4-2β-8)/cos(βπ)Γ(4-2β) h^1-2β (u⃗_h, u⃗_h),

  which proves the second inequality. To prove the left inequality, set u_h:= Φ_hu⃗_h, we rewrite it as

    1/cos(βπ)(_xD^β_Lu_h,_xD^β_Ru_h)_L^2(Ω)=(A_h^βu⃗_h,u⃗_h)
      ≳ h(u⃗_h,u⃗_h)≃ (M_hu⃗_h,u⃗_h)= (u_h, u_h)_L^2(Ω)

  which can be deduced by (<ref>)-(<ref>), where Φ_h = (ϕ_1,⋯,ϕ_M-1).

  Finally, we have to show that

    κ(I+C_ατ^α M_h^-1/2A_h^β M_h^-1/2)   =   λ_max(I+C_ατ^α M_h^-1/2A_h^β
      M_h^-1/2)/λ_min(I+C_ατ^α M_h^-1/2A_h^β M_h^-1/2)
       =   1+C_ατ^αλ_max(M_h^-1/2A_h^β
      M_h^-1/2)/1+C_ατ^αλ_min(M_h^-1/2A_h^β M_h^-1/2)
       ≤    1+C_α(2^4-2β-8)/cos(βπ)Γ(4-2β)τ^α h^-2β.

  This completes the proof based on the spectral equivalence relation (<ref>).



  The estimation (<ref>) is compatible with the correlative result 1+𝒪(τ h^-2) of
  integer order parabolic differential equations.


An important particular case of Theorem <ref> is singled out in the following corollary.


  Let τ be proportional to h^μ with μα≥2β. Then

    κ(C^n_hτ)=𝒪(1).



In what follows, we examine the correctness of (<ref>) concerning Example <ref> with typical
α and β for three specific cases: τ = h, τ = h^2 and τ is fixed (doesn't change along with h).
In under-mentioned tables, λ_min and λ_max respectively indicate the smallest and largest eigenvalues,
κ represents the condition number and ratio is the quotient of the condition number in fine grid divided by that in coarse grid.







It is observed from Tables <ref>-<ref> that numerical results are in good agreement with our theoretical estimation.



§ AMG'S CONVERGENCE ANALYSIS AND AN ADAPTIVE AMG METHOD


Within the section, involving FFT to perform Toeplitz matrix-vector multiplications,
we introduce the so-called Ruge-Stüben or classical AMG method <cit.>
with low algorithmic complexity, fulfill its theoretical investigation,
and then propose an adaptive AMG method through Corollary <ref>.

 The classical AMG method for the linear system (<ref>).
  

    
  Step 1 Perform the Setup phase to the coefficient matrix C^n_hτ.
    

      
  1.1 Set the strength-of-connection tolerance θ;

      
  1.2 Build the ingredients required by a hierarchy of levels, coarsest to finest, including the grid transfer operator P.

    

    
  Step 2 Invoke the classical V(ϱ_1,ϱ_2)-cycle to solve (<ref>) until convergence.
    Below is the description of two-grid V(ϱ_1,ϱ_2)-cycle.
    

      
  2.1 Do ϱ_1 pre-smoothing steps on (<ref>);

      
  2.2 Compute and restrict the residual: r^c=P^T(G^n_hτ - C^n_hτU^n_hτ);

      
  2.3 Solve the residual equation on coarse level: (P^T C^n_hτ P) e^c = r^c;

      
  2.4 Interpolation and correction: U^n_hτ=U^n_hτ+Pe^c;

      
  2.5 Do ϱ_2 post-smoothing steps on (<ref>).

    

  



  In pre- and post-smoothing processes, damped-Jacobi iterative methods are favorable choices,
  which can maintain the low computational cost 𝒪(Mlog M) calculated by FFT.


For theoretical investigations, we rewrite (<ref>) and
the grid transfer operator P in block form regarding a given C/F splitting

    C^n_hτU^n_hτ =
    (
    [ A_FF A_FC; A_CF A_CC ])
    (
    [ u_F; u_C ])
    =(
    [ f_F; f_C ]) = G^n_hτ, 
    P=(
    [ I_FC; I_CC ]),

and introduce the following inner products

    (u_F,v_F)_0,F = (D_FFu_F,v_F), (u,v)_1 =(C^n_hτu,v), (u,v)_2 =(D^-1_hτC^n_hτu,C^n_hτv)

with their associated norms ·_0,F=√((·,·)_0,F) and ·_i=√((·,·)_i) (i=1,2),
where I_CC is the identity operator, D_FF= diag(A_FF) and D_hτ= diag(C^n_hτ).

For simplicity, we here denote C^n_hτ=(c_ij)_(M-1)×(M-1),
and only consider the two-grid V(0,1)-cycle, whose iteration matrix has the form

    M_h,H=S[I-P(P^T C^n_hτ P)^-1P^TC^n_hτ],

where S is a relaxation operator usually chosen as damped-Jacobi or Gauss-Seidel iterative method.

Combining Corollary <ref> and the two-level convergence theory in the work <cit.>, leads to the following lemmas and theorem.


  Under the condition (<ref>), for all e_h∈ℝ^M-1, damped-Jacobi and Gauss-Seidel relaxations satisfy the smoothing property

    Se_h^2_1≤e_h^2_1-σ_1e_h^2_2

  with σ_1 independent of e_h and step sizes h and τ_n.



  On the strength of Theorem A.3.1 and A.3.2 in <cit.>,
  we produce that damped-Jacobi relaxation with parameter 0<ω<2/η satisfies (<ref>) with

    σ_1=ω(2-ωη),

  and Gauss-Seidel relaxation satisfies (<ref>) with

    σ_1=1/(1+γ_-)(1+γ_+),

  both independent of e_h, where

    η≥ρ(D^-1_hτC^n_hτ), γ_-=max_i{1/w_ic_ii∑_j<iw_j|c_ij|}, γ_+=max_i{1/w_ic_ii∑_j>iw_j|c_ij|},

  and w=(w_i) is an arbitrary positive vector with C^n_hτw being also positive. 

  By exploiting (<ref>)-(<ref>), the assumption (<ref>) and Corollary <ref>,
  we conclude that C^n_hτ is strictly diagonally dominant. Recall that C^n_hτw is a positive vector,
  yield γ_-<1, γ_+<1 and

    ρ(D^-1_hτC^n_hτ)≤|D^-1_hτC^n_hτ|_w=max_i{1/w_i∑_jw_j|c_ij|/c_ii}<2,

  which implicitly mean that η, γ_- and γ_+ can be chosen to be independent of h and τ_n, and complete the proof.



  The inequality (<ref>) implies that there exists ϵ>0 such that ρ(D^-1_hτC^n_hτ)=2-3ϵ.
  Then η=2-2ϵ>ρ(D^-1_hτC^n_hτ) and hence the upper bound of parameter ω: 2/η=1/(1-ϵ)>1,
  which suggests that Jacobi relaxation with ω=1 is available in such a case.



  For all symmetric M-matrices, σ_1≤1/η<1 holds for damped-Jacobi relaxation,
  while σ_1∈(1/4,1) for Gauss-Seidel relaxation.



  Under the condition (<ref>) and a given C/F splitting,
  for all e_h=(e_F^T,e_C^T)^T∈ℝ^M-1,
  the direct interpolation I_FC satisfies

    e_F-I_FCe_C^2_0,F≤σ_2e_h^2_1

  with σ_2 independent of e_h, h and τ_n.



  According to Theorem A.4.3 in <cit.>, I_FC satisfies (<ref>) with
  σ_2 of the form regarding a given C/F splitting

    σ_2 ≥max_i∈ F{∑_j∈ N_ic_ij/∑_j∈ C_ic_ij}

  independent of e_h, where N_i={j i:c_ij0}, C_i is the subset of N_i
  whose values will be used to interpolate at F-point i. As a result of (<ref>)
  and the fact that c_ij (j∈ N_i) are all negative, the following relation holds: σ_2 > 1.

  Notice here that the classical Ruge-Stüben based coarsening strategy generates
  at least one of points i-1 and i+1 to be C-points and strongly influence i
  — viz. it retains i-1∈ C_i or i+1∈ C_i. Therefore, it can be seen that

    ∑_j∈ N_ic_ij/∑_j∈ C_ic_ij<-c_ii/c_ii-1=-c_ii/c_ii+1
      =-16cos(βπ)Γ(4-2β)+6Γ(3-α)τ^α_nh^-2β(2^4-2β-
      8)/4cos(βπ)Γ(4-2β)+6Γ(3-α)τ^α_nh^-2β(3^3-2β-2^5-2β+7),

  indicating that σ_2 is independent of h and τ_n by plugging (<ref>), and thus prove the theorem.



  Let any C/F splitting be given. Under the condition (<ref>),
  there exist positive constants σ_1 and σ_2 independent of h and τ_n and satisfying σ_2>1>σ_1,
  such that a uniform two-grid convergence is achieved as follows

    M_h,H_1≤√(1-σ_1/σ_2).




  The proof of this result is straightforward and is based on Theorem A.4.1 and A.4.2 in <cit.>, Lemma <ref> and <ref>.


Now observe from Theorem <ref> that, despite the independence of h and τ_n,
σ_2 relies ruinously on θ in Step 1.1 of Algorithm <ref>
due to the fact that C^n_hτ is nearly dense leading to a quite complicated adjacency graph.
In addition, it is found that σ_2(θ) may be much larger than 1 as θ approaches zero,
with that comes a sharp pullback in convergence rate.
Hence, an appropriate θ is a critical component of Algorithm <ref> to handle fractional diffusion equations.

We now turn to reveal a reference formula on θ. Note the heuristic that the distribution of ratios of off-diagonal elements
relative to the maximum absolute off-diagonal element (namely the minor diagonal element for C^n_hτ) plays a major role
in the choice of θ. Since C^n_hτ is a symmetric Toeplitz matrix from Corollary <ref>,
its first row involving all off-diagonal elements of C^n_hτ is deserved to be the representative row.
Taking β=0.8 as an example, Fig. <ref> shows the distribution of the ratios c_1j/c_12 (j≥2),
which reminds us of the attenuation in off-diagonal elements, states

    c_13/c_12≈ 0.160426,  
     c_1j/c_12 < c_14/c_12≈ 0.034394, j=5,6,⋯,M-1,

and suggests that c_1j (j≥4) should be viewed as weak couplings (wouldn't be used for interpolation)
because they are less than 5% of c_12. Besides, for a better complexity and higher efficiency,
only the nearest neighbors are potentially used to limit the interpolation matrix on each grid level
to at most 3 coefficients per row, although c_13 reaches around 16% of c_12.



It thus appears that the strength-of-connection tolerance θ should be of the form

    θ = c_13/c_12+ϵ_0,

where ϵ_0 is some small number, which can be chosen to be 10^-5 in one-dimensional realistic problems.



As is known, Algorithm <ref> is much more expensive for well-conditioned problems than basic iterative techniques,
such as conjugate gradient (CG) or (plain) Jacobi iterative method. For the purpose of solving (<ref>)
in an optimal way, an adaptive AMG method is proposed below by combining Algorithm <ref>, the reference formula (<ref>)
and the condition number estimation (<ref>)
in Corollary <ref> as the clear distinction to adaptively pick an appropriate solver.

 An adaptive AMG method S_ad for the linear system (<ref>).
  

    
  Step 1 If the condition (<ref>) is unsatisfied, then goto Step 2,
    else set S_ad as the CG or Jacobi iterative method;

    
  Step 2 Set S_ad as the classical AMG method described in Algorithm <ref>,
    with θ chosen via the reference formula (<ref>).

  




§ PERFORMANCE EVALUATION


Let us illustrate the effectiveness of Algorithms <ref> and <ref>. Numerical experiments are performed in a 64 bit Fedora 18 platform,
double precision arithmetic on Intel Xeon (W5590) with 24.0 GB RAM, 3.33 GHz,
with an -O2 optimization parameter. In the following tables, dashed entries (-) indicate the solutions
either diverge or fail to converge after 1000 iterations,
Its is the number of iterations until the stopping criterion 10^-12 is reached,
T_c represents the CPU time including both Setup and Solve phases with second as its unit,
C_g and C_o respectively denote grid and operator complexities, which are defined as
sums of the number of degrees of freedom and nonzero elements on all grid levels divided by those of the finest grid level,
and used as measures for memory requirements, aritmetic operations and the execution time in Setup and Solve phases.


  Comparisons of the classical AMG over CG and Jacobi iterative methods
  for the case when (<ref>) is satisfied with two different fractional orders.




As expected, the results in Table <ref> show that Jacobi, CG and AMG methods
are robust with respect to the mesh size and fractional order,
which indicates indirectly the correctness of (<ref>).
In addition, CG method runs 3.28 and 3.03 times faster than Jacobi and AMG methods for M=4096 and α=β=0.6, respectively.


  Comparisons between the classical AMG method and CG method for the case when (<ref>) is unsatisfied.




As shown in Table <ref>, AMG method converges robustly regarding to the mesh size and may be weakly dependent of β,
while the number of iterations of CG method is quite unstable, and sometimes CG method even break down.
Furthermore AMG method runs 12.24 times faster than CG method for M=4096 and β=0.8.



Table <ref> shows the results of τ=h.
Despite the advantage in computational cost and robustness over CG method,
AMG method is nearly independent of α and β in this circumstance.
Meanwhile, by an investigation in terms of number of iterations in Tables <ref> and <ref>,
CG method converges faster because of the improvement in condition number
from 𝒪(h^-2β) to 𝒪(h^α-2β).


  Comparisons of S_ad over the classical AMG and CG methods when the i-th time step size τ_i is chosen to be

    τ_i = {    h^2,    i=1,⋯,K_1 
        1/32, i=K_1+1,⋯,K_1+K_2
        ..





We can observe from Table <ref> that S_ad and AMG methods are fairly robust as to the mesh size,
roughly 10 and 6 on the average. Yet the average number of iterations of CG method varies from 85 to 142.
Moreover S_ad has a considerable advantage over others in CPU time,
runs 1.72 and 6.09 times faster than AMG and CG methods for M=2048 and K_2=100.


  Analyze effects of the strength-of-connection tolerance θ on the performance of the classical AMG method.






It is seen from Tables <ref> and <ref> that there is a unique threshold θ_0 independent of h
which guarantees the robustness of the classical AMG method, and makes number of iterations of the classical AMG monotonically decreasing
when θ<θ_0, or even the classical AMG possibly diverge when θ is small enough,
e.g., θ_0=0.16043 and θ_0=0.00685 for cases β=0.8 and β=0.99.
By direct calculations, we have c_13/c_12≈ 0.160426 and c_13/c_12≈ 0.006846.
Utilizing the relation (<ref>) and ϵ_0=10^-5,
the corresponding values of θ are respectively larger than those of θ_0.
This confirms the reasonability of the reference formula (<ref>).



§ CONCLUSION


In this paper, we propose the variational formulation for a class of time-space Caputo-Riesz fractional diffusion equations,
prove that the resulting matrix is a symmetric Toeplitz matrix, an M-matrix by appending a very weak constraint
and its condition number is bounded by 1+𝒪(τ^α h^-2β),
introduce the classical AMG method and prove rigorously that its convergence rate is independent of time and space step sizes,
provide explicitly a reference formula of the strength-of-connection tolerance to guarantee the robustness and predictable behavior of AMG method in all cases,
and develop an adaptive AMG method via our condition number estimation to decrease the computation cost.
Numerical results are all in conformity with the theoretical results, and verify the reasonability of the reference formula
and the considerable advantage of the proposed adaptive AMG algorithm over other traditional iterative methods, e.g. Jacobi,
CG and the classical AMG methods.



§ ACKNOWLEDGMENTS


This work is under auspices of National Natural Science Foundation of China (11571293, 11601460, 11601462)
and the General Project of Hunan Provincial Education Department of China (16C1540, 17C1527).




§ REFERENCES


99

l-002
F. Liu, P. Zhuang, V. Anh, I. Turner, K. Burrage,
Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation,
Appl. Math. Comput. 191 (1) (2007) 12-20.

d-001
Z. Q. Ding, A. G. Xiao, M. Li,
Weighted finite difference methods for a class of space fractional partial differential equations with variable coefficients,
J. Comput. Appl. Math. 233 (8) (2010) 1905-1914.

g-001
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.

y-007
S. P. Yang, A. G. Xiao, X. Y. Pan,
Dependence analysis of the solutions on the parameters of fractional delay differential equations,
Adv. Appl. Math. Mech. 3 (5) (2011) 586-597.

c-002
X. N. Cao, J. L. Fu, H. Huang,
Numerical method for the time fractional Fokker-Planck equation,
Adv. Appl. Math. Mech. 4 (6) (2012) 848-863.

w-002
H. Wang, T. S. Basu,
A fast finite difference method for two-dimensional space-fractional diffusion equations,
SIAM J. Sci. Comput. 34 (5) (2012) A2444-A2458.

c-001
M. H. Chen, W. H. Deng, Y. J. Wu,
Superlinearly convergent algorithms for the two-dimensional space-time Caputo-Riesz fractional diffusion equation,
Appl. Numer. Math. 70 (2013) 22-41.

w-003
D. L. Wang, A. G. Xiao, H. L. Liu,
Dissipativity and stability analysis for fractional functional differential equations,
Fract. Calc. Appl. Anal. 18 (6) (2015) 1399-1422.

w-001
D. L. Wang, A. G. Xiao, W. Yang,
Maximum-norm error analysis of a difference scheme for the space fractional CNLS,
Appl. Math. Comput. 257 (2015) 241-251.

y-004
W. Yang, D. L. Wang, L. Yang,
A stable numerical method for space fractional Landau-Lifshitz equations,
Appl. Math. Lett. 61 (2016) 149-155.

e-001
V. J. Ervin, J. P. Roop,
Variational formulation for the stationary fractional advection dispersion equation,
Numer. Methods Partial Differ. Equ. 22 (3) (2006) 558-576.

z-002
H. Zhang, F. Liu, V. Anh,
Galerkin finite element approximation of symmetric space-fractional partial differential equations,
Appl. Math. Comput. 217 (6) (2010) 2534-2545.

b-002
K. Burrage, N. Hale, D. Kay,
An efficient implicit FEM scheme for fractional-in-space reaction-diffusion equations,
SIAM J. Sci. Comput. 34 (4) (2012) A2145-A2172.

b-001
W. P. Bu, Y. F. Tang, J. Y. Yang,
Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations,
J. Comput. Phys. 276 (2014) 26-38.

m-001
K. Mustapha, B. Abdallah, K. M. Furati,
A discontinuous Petrov-Galerkin method for time-fractional diffusion equations,
SIAM J. Numer. Anal. 52 (5) (2014) 2512-2529.

f-001
L. B. Feng, P. Zhuang, F. Liu, I. Turner, Y. T. Gu,
Finite element method for space-time fractional diffusion equation,
Numer. Algor. 72 (3) (2016) 749-767.

b-004
W. P. Bu, A. G. Xiao, W. Zeng,
Finite difference/finite element methods for distributed-order time fractional diffusion equations,
J. Sci. Comput. 72 (1) (2017) 422-441.

l-003
F. Liu, P. Zhuang, I. Turner, K. Burrage, V. Anh,
A new fractional finite volume method for solving the fractional diffusion equation,
Appl. Math. Model. 38 (15-16) (2014) 3871-3878.

l-001
Y. M. Lin, C. J. Xu,
Finite difference/spectral approximations for the time-fractional diffusion equation,
J. Comput. Phys. 225 (2) (2007) 1533-1552.

z-001
M. Zayernouri, W. R. Cao, Z. Q. Zhang, G. E. Karniadakis,
Spectral and discontinuous spectral element methods for fractional delay equations,
SIAM J. Sci. Comput. 36 (6) (2014) B904-B929.

y-005
Y. Yang, Y. P. Chen, Y. Q. Huang,
Spectral-collocation method for fractional Fredholm integro-differential equations,
J. Korean Math. Soc. 51 (1) (2014) 203-224.

y-006
Y. Yang, Y. P. Chen, Y. Q. Huang,
Convergence analysis of the Jacobi spectral-collocation method for fractional integro-differential equations.
Acta Math. Sci. 34B (3) (2014) 673-690.

y-001
Y. Yang,
Jacobi spectral Galerkin methods for fractional integro-differential equations,
Calcolo 52 (4) (2015) 519-542.

y-002
Y. Yang,
Jacobi spectral Galerkin methods for Volterra integral equations with weakly singular kernel,
Bull. Korean Math. Soc. 53 (1) (2016) 247-262.

y-003
S. P. Yang, A. G. Xiao,
An efficient numerical method for fractional differential equations with two Caputo derivatives,
J. Comput. Math. 34 (2) (2016) 113-134.

l-004
S. L. Lei, H. W. Sun,
A circulant preconditioner for fractional diffusion equations,
J. Comput. Phys. 242 (2013) 715-725.

m-002
T. Moroney, Q. Q. Yang,
A banded preconditioner for the two-sided, nonlinear space-fractional diffusion equation,
Comput. Math. Appl. 66 (5) (2013) 659-667.

j-001
J. H. Jia, H. Wang,
Fast finite difference methods for space-fractional diffusion equations with fractional derivative boundary conditions,
J. Comput. Phys. 293 (2015) 359-369.

g-002
X. M. Gu, T. Z. Huang, X. L. Zhao, H. B. Li, L. Li,
Strang-type preconditioners for solving fractional diffusion equations by boundary value methods,
J. Comput. Appl. Math. 277 (2015) 73-86.

d-002
M. Donatelli, M. Mazza, S. Serra-Capizzano,
Spectral analysis and structure preserving preconditioners for fractional diffusion equations,
J. Comput. Phys. 307 (2016) 262-279.

p-001
H. K. Pang, H. W. Sun,
Multigrid method for fractional diffusion equations,
J. Comput. Phys. 231 (2012) 693-703.

b-003
W. P. Bu, X. T. Liu, Y. F. Tang, J. Y. Yang,
Finite element multigrid method for multi-term time fractional advection diffusion equations,
Int. J. Model. Simul. Sci. Comput. 6 (2015) 1540001.

j-002
Y. J. Jiang, X. J. Xu,
Multigrid methods for space fractional partial differential equations,
J. Comput. Phys. 302 (2015) 374-392.

c-003
L. Chen, R. H. Nochetto, E. Otárola, A. J. Salgado,
Multilevel methods for nonuniformly elliptic operators and fractional diffusion,
Math. Comput. 85 (302) (2016) 2583-2607.

z-003
X. Zhao, X. Z. Hu, W. Cai, G. E. Karniadakis,
Adaptive finite element method for fractional differential equations using hierarchical matrices,
arXiv: 1603.01358v2.

x-001
J. Xu,
Iterative methods by space decomposition and subspace correction,
SIAM Rev. 34 (4) (1992) 581-613.

r-001
J. W. Ruge, K. Stüben,
Algebraic multigrid, in Multigrid Methods,
Front. Appl. Math. 3 (1987) 73-130.

s-001
K. Stüben,
An introduction to algebraic multigrid,
in Multigrid, U. Trottenberg, C. W. Oosterlee and A. Schüller, eds., Academic Press, Singapore, 2001.








 