# Finite Element Multigrid Method for the Boundary Value Problem of Fractional Advection Dispersion Equation

Authors: Zhiqiang Zhou; Hongying Wu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The stationary fractional advection dispersion equation is discretized by linear finite element scheme, and a full V-cycle multigrid method (FV-MGM) is proposed to solve the resulting system. Some useful properties of the approximation and smoothing operators are proved. Using these properties we derive the convergence results in bothL2 norm and energy norm for FV-MGM. Numerical examples are given to demonstrate the convergence rate and efficiency of the method.


---

## Body

## 1. Introduction

We investigate the finite element full V-cycle multigrid method (FV-MGM) to the boundary value problem of linear stationary fractional advection dispersion equation (FADE).Problem 1.
GivenΩ=(0,1),f:Ω¯→ℝ, find u:Ω¯→ℝ such that
(1)-12D(Dx-β0+Dx1-β)Du+c(x)u=f(x),inΩ,u=0,on∂Ω,
where 0<β<1, c(x)≥0, D represents a single spatial derivative, and Dx-β0 and Dx1-β represent Riemann-Liouville left and right fractional integral operators, respectively [1–3].FADE has been used in modeling physical phenomena exhibiting anomalous diffusion, that is, diffusion not accurately modeled by the usual advection dispersion equation [4–7]. For example, solutes moving through aquifers do not generally follow a Fickian, second-order, and governing equation because of large deviations from the stochastic process of Brownian motion [8–10]. Many scholars developed numerical methods, including finite difference method [11], finite element method [12–14], spectral method [15] and moving collocation method [16] to solve FADEs. Most of them used Gauss elimination method or conjugate gradient norm residual method to solve the resulting system, so the computational complexity is 𝒪(N3) or 𝒪(Nlog2N). To date only a few of them considered MGM numerical methods. For example, Pang and Sun [11] developed an MGM to solve the linear system with Toeplitz-like structure, but the fractional derivatives are defined in the Grünwald-Letnikov form and the discretized system is obtained by difference scheme. It motivates us to design a fast MGM algorithm to deal with FE equations of FADE.In this paper, we follow the ideas in [17, 18] to develop a FV-MGM for solving the resulting system of Problem 1 discretized by linear finite element method. By selecting appropriate iteration operator and iteration numbers, we prove that FV-MGM has the same convergence rate as classic FEM and the computational cost increases linearly with respect to the increasing of unknown variables.The remaining of this paper is organized as follows. The next section recalls the variational formulation of the above FADE and corresponding convergence results. In Section3 we describe the multigrid algorithm, estimate the spectral radius of FE equations, show the properties of approximation and soothing operators, and prove the convergence theorems for FV-MGM. Numerical examples demonstrating convergence rate and computational work are presented in Section 4. Concluding remarks are given in the final section.


## 2. Variational Formulation and Convergence Results

Letα:=(2-β)/2, so that 1/2<α≤1. Define the associate bilinear form B:H0α(Ω)×H0α(Ω)→ℝ as
(2)B(u,v):=12〈D0x-βDu,Dv〉+12〈Dx1-βDu,Dv〉+(c(x)u,v),
and linear functional  F:H0α(Ω)→ℝ  as
(3)F(v):=〈f,v〉,
where (·,·) denotes the inner product on L2(Ω) and 〈·,·〉 the duality pairing of H-μ(Ω) and H0μ(Ω),μ≥0.Thus, the Galerkin variational solution of Problem1 may be defined as follows.Definition 2 (variational solution).
A functionu∈H0α(Ω) is a variational solution of Problem 1 provided that
(4)B(u,v)=F(v),∀v∈H0α(Ω).Ervin and Roop ([17], 2006) proved that the bilinear form B(·,·) satisfies the coercivity and continuity; that is, there exist two positive constants C1  and  C2 such that
(5)B(u,u)≥C1∥u∥Hα(Ω)2,∀u∈H0α(Ω),(6)B(u,v)≤C2∥u∥Hα(Ω)∥v∥Hα(Ω),∀u,v∈H0α(Ω).
Hence, they claimed that there exists a unique solution to (4). Note that from (5) and (6) we have norm equivalence of ∥·∥Hα(Ω)  and energy norm  ∥·∥E=B(·,·)1/2, that is,
(7)C1∥v∥Hα(Ω)≤∥v∥E≤C2∥v∥Hα(Ω),∀v∈H0α(Ω).LetSh denote a partition of Ω such that Ω¯={⋃‍K:K∈Sh}. Assume that there exists a positive constant ρ such that ρhk≤hK≤hk, where hk=maxK∈ShhK. Let Pm(K) denote the space of polynomials of degree less than or equal to m on K∈Sh. Associated with Sh, define the finite-dimensional subspace Xh⊂H0α(Ω) as
(8)Xh:={v∈H0α(Ω)∩C0(Ω¯):v∣K∈Pm-1(K),∀K∈Sh}.
Let uh be the solution to the finite-dimensional variational problem:
(9)B(uh,vh)=F(vh),∀vh∈Xh.Note that the existence and uniqueness of solution to (9) follow from the fact that Xh is a subset of the space H0α(Ω)∩H0m(Ω). Ervin and Roop have obtained convergence estimate in the energy norm ∥·∥E. At the same time, under the assumption concerning the regularity of the solution to the adjoint problem of Problem 1 [17, Assumption 4.1], convergence estimate in the L2 norm was proved.Theorem 3 (see [17, Corollary 4.3]).
Letu∈H0α(Ω)⋂‍Hr(Ω)(α≤r≤m) solve (4) and uh solve (9). Then there exists a constant  C3>0 such that  e=u-uh satisfies
(10)∥e∥Hα(Ω)≤C3hkr-α∥u∥Hr(Ω).Theorem 4 (see [17, Theorem 4.4]).
Letu∈H0α(Ω)⋂‍Hr(Ω)(α≤r≤m) solve (4) and uh solve (9), where m-1 is the degree of Galerkin finite element model. Then, if the regularity of the solution to the adjoint problem is satisfied, there exists a constant C4 (or C4′) such that the error e=u-uh satisfies
(11)∥e∥L2(Ω)≤C4hkα∥u∥E,α≠34,∥e∥L2(Ω)≤C4′hkα-ε∥u∥E,α=34,0<∀ε<12.From Theorems3 and 4 and the equivalence of  ∥·∥Hα(Ω) and ∥·∥E, we have the error estimates in the L2 norm.Theorem 5.
Under the same assumptions in Theorem4, one has
(12)∥e∥L2(Ω)≤C5hkr∥u∥Hr(Ω),α≠34,∥e∥L2(Ω)≤C5hkr-ε∥u∥Hr(Ω),α=34,0<∀ε<12,
where C5=C2C3C4 (or C5=C2C3C4′).


## 3. Multigrid Method for FADE

For the simpleness, we only discuss the case of linear finite element, and it is necessary to restrict the mesh partition. Let{𝒯k} denote a sequence of partitions of Ω¯ and hk be the mesh size of 𝒯k. From now on, we assume that {𝒯k} is a quasiuniform family, that is,
(13)ρhk≤xi-xi-1≤hk,i=1,2,…,Nk,
with positive constant ρ. At the same time, let 𝒯k be obtained from 𝒯k-1 via a regular subdivision, that is, 𝒯k-1⊂𝒯k for k≥2. Let Vk denote C0 piecewise linear functions with respect to 𝒯k that vanish on ∂Ω.


### 3.1. Mesh-Dependent Norms

Definition 6.
The mesh-dependent inner product(·,·)k on  Vk  is defined by
(14)(v,w)k=hk∑i=1Nk-1v(xi)w(xi).Definition 7.
The operatorAk:Vk→Vk is defined by
(15)(Akv,w)k=B(v,w),∀v,w∈Vk.In terms of the operatorAk, the discretized equation of (9) can be written as
(16)Akuk=fk,
where fk∈Vk satisfies
(17)(fk,vk)=(f,v),∀v∈Vk.Remark 8.
(i) From the Riesz representation theorem, we point out that the operatorAk is defined uniquely. (ii) Since B(v,w) is symmetric and satisfies the proposition of coercivity, Ak is symmetric positive definite woth respect to (·,·)k.The mesh-dependent norms are defined by(18)∥|v|∥s,k=(Aksv,v)k.
Observe that the energy norm ∥·∥E coincides with ∥|·|∥1,k norm on Vk. Similarly,  ∥|·|∥0,k  is the norm associated with the mesh-dependent inner product (14). The following lemma shows that ∥|·|∥0,k is equivalent to the L2 norm.Lemma 9.
The norm∥·∥L2(Ω) is equivalent to mesh-depentent norm ∥|·|∥0,k, that is,
(19)ρ∥|v|∥0,k≤∥v∥L2(Ω)≤∥|v|∥0,k,∀v∈Vk.
Particularly,  ∥|·|∥0,k=∥·∥L2(Ω)  for the uniform mesh.Proof.
The proof is analogous to Lemma6.2.7 in reference [18].In order to estimate the spectral radius,Λ(Ak), of Ak, we need the following norm interpolation lemma (see [19, Theorem 1.3.7]).Lemma 10 (interpolation theorem).
Assume thatΩ  is an open set of  ℝd with a Lipchitz continuous boundary. Let s1<s2 be two real numbers, and  μ=(1-θ)s1+θs2,0≤θ≤1. Then there exists a constant C6>0 such that
(20)∥u∥Hμ(Ω)≤C6∥u∥Hs1(Ω)1-θ∥u∥Hs2(Ω)θ.Lemma 11 (spectral radius theorem).
We have the estimation for spectral radiusΛ(Ak):
(21)Λ(Ak)≤C7hk-2α,
where C7 is a positive constant independent of k.Proof.
Letϕ∈Vk. Applying the continuity (6) of  B(·,·), norm interpolation Lemma 10, and inverse estimation, we have
(22)B(ϕ,ϕ)≤C2∥ϕ∥Hα(Ω)2≤C2(C6∥ϕ∥H0(Ω)1-α·∥ϕ∥H1(Ω)α)2≤C2(C6∥ϕ∥H0(Ω)1-α·C8hk-α∥ϕ∥H0(Ω)α)2=Chk-2α∥ϕ∥H0(Ω)2=Chk-2α∥ϕ∥L2(Ω)2,
where C8>0 is the constant of inverse estimate and C=C2C62C82.
Letλ be an eigenvalue of Ak with corresponding eigenvector ϕ∈Vk. From Lemma 9, we have
(23)B(ϕ,ϕ)=(Akϕ,ϕ)k=λ(ϕ,ϕ)k=λ∥|ϕ|∥0,k2=λ∥ϕ∥L2(Ω)2.
Combining (22) and (23), the following inequality holds:
(24)λ=B(ϕ,ϕ)∥ϕ∥L2(Ω)2≤Chk-2α∥ϕ∥L2(Ω)2∥ϕ∥L2(Ω)2=C7hk-2α,
where C7=C2C62C82.


### 3.2. The V-Cycle Multigrid Algorithm

Firstly, we give the intergrid transfer operators which play a very important role in convergence analysis.Definition 12 (intergrid transfer operator).
The coarse-to-fine intergrid transfer operatorIk-1k:Vk-1→Vk is taken to be the natural injection, that is,
(25)Ik-1kv=v,∀v∈Vk-1.
The fine-to-coarse operator Ikk-1:Vk→Vk-1 is defined by
(26)(Ikk-1w,v)k-1=(w,Ik-1kv)k=(w,v)k,∀v∈Vk-1,w∈Vk.Now, we describe thekth level of  V-cycle MGM (V-MGM) and full V-cycle MGM. Let m be the smoothing number, nk  the iteration number of kth level V-MGM, and Λk be a parameter dependent on k. Denote MG(k,z0,g) the approximate solution of the Akz=g obtained from the kth level iteration with initial guess z0. The discussion of parameters Λk and nk is left to the next subsection and Section 4.Algorithm 13 (thekth level V-cycle multigrid algorithm).
BEGIN:
Fork=1:  MG(1,z0,g)=A1-1g.
Fork>1, MG(k,z0,g) is obtained recursively in three steps.
Presmoothing Step. For1≤l≤m, let zl=zl-1+(1/Λk)(g-Akzl-1).
Error Correction.
g
¯
=
I
k
k
-
1
(
g
-
A
k
z
m
),
q
1
=
M
G
(
k
-
1,0
,
g
¯
),
z
m
+
1
=
z
m
+
I
k
-
1
k
q
1.
Postsmoothing step. Form+1≤l≤2m+1, let zl=zl-1+(1/Λk)(g-Akzl-1).
The output of thekth level iteration is  MG(k,z0,g):=z2m+1.
END.Algorithm 14 (the full V-cycle multigrid algorithm).
BEGIN:
Fork=1, u^1=A1-1f1.
Fork≥2,
u
0
k
=
I
k
-
1
k
u
^
k
-
1,
u
l
k
=
M
G
(
k
,
u
l
-
1
k
,
f
k
)
,
1
≤
l
≤
n
k,
u
^
k
=
u
n
k
k.
END.


### 3.3. Approximation and Smoothing Properties

In this subsection, we prove some properties of projection operatorPk-1 and smoothing operator Rk, which are the key ingredients for V-MGM and FV-MGM algorithm.Definition 15.
LetPk:V→Vk be the orthogonal projection with respect to B(·,·), that is,
(27)B(v,w)=B(Pkv,w),∀w∈Vk.Definition 16.
Define the relaxation operator:(28)Rk:=I-1ΛkAk.Throughout this paper, we assume thatΛk denotes some upper bound for the spectral radius of Ak satisfying  Λk≤C7hk-2α.Lemma 17.
Letv∈Vk, g¯=Ikk-1Akv,  and q∈Vk-1 satisfy the equation Ak-1q=g¯. Then q=Pk-1v.Proof.
Forw∈Vk-1,
(29)B(q,w)=(Ak-1q,w)k-1=(g¯,w)k-1=(Ikk-1Akv,w)k-1=(Akv,w)k=B(v,w).
Therefore, from the definition of Pk-1 we have q=Pk-1v.Remark 18.
Letg=Akv. From the Riesz representation theorem there exists a unique f∈H-α(Ω) such that
(30)(g,w)k=〈f,w〉,∀w∈Vk,
which means that v is a finite element solution of variational problem (4) on Vk. Observing that
(31)(g¯,w)k-1=(Ikk-1g,w)k-1=(g,w)k=〈f,w〉,∀w∈Vk-1,
we claim that q=Pk-1v is also a finite element solution of variational problem (4) on Vk-1. Therefore, noting  Vk-1⊂Vk (since 𝒯k-1⊂𝒯k) and applying Theorem 4, we have
(32)∥v-Pk-1v∥L2(Ω)≤C4hkα∥v-Pk-1v∥E,α≠34,∥v-Pk-1v∥L2(Ω)≤C4′hkα-ε∥v-Pk-1v∥E,α=34.Lemma 19.
There exists a positive constantC9  such that
(33)∥|(I-Pk-1)v|∥0,k≤C9hkα∥|(I-Pk-1)v|∥1,k,∀v∈Vk,α≠34,∥|(I-Pk-1)v|∥0,k≤C9hkα-ε∥|(I-Pk-1)v|∥1,k,∀v∈Vk,α=34,0<ε<12.Proof.
We only prove the caseα≠3/4, and the proof for the other case is analogous. Applying Remark 18 and Lemma 9 we have
(34)∥|(I-Pk-1)v|∥0,k≤1ρ∥(I-Pk-1)v∥L2(Ω)≤C4ρhkα∥(v-Pk-1v)∥E=C9hkα∥|(I-Pk-1)v|∥1,k,
where C9=C4/ρ.Lemma 20.
There exists a positive constantC9 such that
(35)∥|(I-Pk-1)v|∥1,k≤C9hkα∥|v|∥2,k,∀v∈Vk,α≠34,(36)∥|(I-Pk-1)v|∥1,k≤C9hkα-ε∥|v|∥2,k,∀v∈Vk,α=34,0<ε<12.Proof.
Applying projection operatorPk-1, generalized Cauchy-Schwarz inequality (see [18, Lemma 6.2.10]), and Lemma 19(37)∥|(I-Pk-1)v|∥1,k2=∥(I-Pk-1)v∥E2=B(v-Pk-1v,v-Pk-1v)=B(v-Pk-1v,v)≤∥|(I-Pk-1)v|∥0,k∥|v|∥2,k≤C9hkα∥|(I-Pk-1)v|∥1,k∥|v|∥2,k=C9hkα∥|(I-Pk-1)v|∥1,k∥|v|∥2,k
if α≠3/4. Therefore, dividing through by ∥|(I-Pk-1)v|∥1,k yields (35).  The case α=3/4 is analogous.Lemma 21.
Letm be the number of smoothing steps. Then
(38)B((I-Rk)Rk2mv,v)≤12mB((I-Rk2m)v,v).Proof.
See Lemma6.6.7 in [18].Lemma 22.
Letm be the number of smoothing steps. Then there exists a positive constant C10 such that
(39)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)≤C10mB((I-Rk2m)v,v).Proof.
Applying Lemmas20, 11, and 21, we have
(40)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)=∥(I-Pk-1)Rkmv∥E2=∥|(I-Pk-1)Rkmv|∥1,k2≤C92hk2α∥|Rkmv|∥2,k2=C92hk2α(Ak2Rkmv,Rkmv)k=C92hk2αB(AkRkmv,Rkmv)=C92hk2αΛkB((I-Rk)Rkmv,Rkmv)≤C92C7B((I-Rk)Rkmv,Rkmv)=C92C7B((I-Rk)Rk2mv,v)≤C92C7(2m)B((I-Rk2m)v,v)=C10mB((I-Rk2m)v,v)
if α≠3/4, where C10=C92C7/2. Note that (39) also holds for α=3/4 since ε>0.


### 3.4. Convergence Analysis

In this subsection, we prove convergence theorems for V-cycle FE multigrid method.Definition 23.
The error operator of V-cycle multigridEk:Vk→Vk is defined recursively by
(41)E1=0,Ek=Rkm[I-(I-Ek-1)Pk-1]Rkm,k>1.Proposition 24.
The operatorEk has the following propositions (see [18, Lemma 6.6.2 and Proposition 6.6.9]):(i)
if  z,g∈Vk  satisfy  Akz=g, then(42)Ek(z-z0)=z-MG(k,z0,g),(ii)
Ek  is symmetric positive semidefinite with respect toB(·,·)  for  k≥1, that is,(43)B(Ekv,w)=B(v,Ekw),∀v,w∈Vk,B(Ekv,v)≥0.Theorem 25.
Letm be the number of smoothing steps. Then
(44)B(Ekv,v)≤C10m+C10B(v,v),∀v∈Vk.Proof.
The proof is by induction. Fork=1, (44) is trivially true since Ek=0. Assume that (44) holds for k-1. Let γ=C10/(m+C10); therefore 1-γ=γm/C10. Applying induction hypothesis, Proposition 24, and Lemma 22, we have
(45)B(Ekv,v)=B(Rkmv,Rkmv)-B(Pk-1Rkmv,Pk-1Rkmv)+B(Ek-1Pk-1Rkmv,Pk-1Rkmv)≤B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB(Pk-1Rkmv,Pk-1Rkmv)=(1-γ)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB(Pk-1Rkmv,Pk-1Rkmv)≤(1-γ)C10mB((I-Rk2m)v,v)+γB(Rkmv,Rkmv)=γB((I-Rk2m)v,v)+γB(Rkmv,Rkmv)=γB(v,v),
which ends the proof.Corollary 26.
Letm be the number of smoothing steps. Then one has the estimation for spectral radius Λ(Ek):
(46)Λ(Ek)≤C10m+C10.Theorem 27 (convergence of thekth level iteration for V-cycle).
Letm be the number of smoothing steps. Then
(47)∥z-MG(k,z0,g)∥E≤C10m+C10∥z-z0∥E.
Hence, the kth level iteration for any m is a contraction with the contraction number independent of k.Proof.
From (42) in Proposition 24, it suffices to show
(48)∥Ekv∥E≤C10m+C10∥v∥E,
which can be proved by Corollary 26.Theorem 28 (full multigrid convergence).
Letu∈H0α(Ω)⋂‍Hr(Ω)(α≤r≤2) solve (4). If the iteration numbers in full multigrid algorithm nk  are  large enough, there exists a positive constant C11 such that
(49)∥uk-u^k∥E≤C11hkr-α∥u∥Hr(Ω).Proof.
For the simpleness we assume thatnk≥n≥1. Define e^k=uk-u^k. In particular, e^1=0. Let γ=C10/(m+C10) (C10 is the constant in Theorem 27). We have
(50)∥e^k∥E≤γn∥uk-u^k-1∥E≤γn(∥uk-u∥E+∥u-uk-1∥E+∥uk-1-u^k-1∥E)≤γn[C2C3hkr-α∥u∥Hr(Ω)+C2C3hk-1r-α∥u∥Hr(Ω)+∥e^k-1∥E]=γn[3C2C3hkr-α∥u∥Hr(Ω)+∥e^k-1∥E].
In the following analysis, we assume that  2γn<1 (therefore n>logγ1/2). By iterating the above inequality and from continuity (6) and Theorem 3, we have
(51)∥e^k∥E≤3C2C3γnhkr-α∥u∥Hr(Ω)+3C2C3γ2nhk-1r-α∥u∥Hr(Ω)+⋯+3C2C3γknh1r-α∥u∥Hr(Ω)≤3C2C3γn1-2γnhkr-α∥u∥Hr(Ω)=C11hkr-α∥u∥Hr(Ω),
where C11=3C2C3γn/(1-2γn).The following Corollary is a natural conclusion of Theorems5 and 28.Corollary 29.
Under the assumptions in Theorem28, the convergence rate of multigrid solution u^k in L2 norm is r;  that is, there exists a constant C12>0 such that
(52)∥u^k-u∥L2(Ω)≤C12hkr∥u∥Hr(Ω),α≠34,∥u^k-u∥L2(Ω)≤C12hkr-ε∥u∥Hr(Ω),α=34,0<∀ε<12.We end this section with a proposition that the work involved in the full multigrid algorithm is𝒪(dk), where dk=Nk-2=dimVk (see [18, Proposition 6.7.4]).


## 3.1. Mesh-Dependent Norms

Definition 6.
The mesh-dependent inner product(·,·)k on  Vk  is defined by
(14)(v,w)k=hk∑i=1Nk-1v(xi)w(xi).Definition 7.
The operatorAk:Vk→Vk is defined by
(15)(Akv,w)k=B(v,w),∀v,w∈Vk.In terms of the operatorAk, the discretized equation of (9) can be written as
(16)Akuk=fk,
where fk∈Vk satisfies
(17)(fk,vk)=(f,v),∀v∈Vk.Remark 8.
(i) From the Riesz representation theorem, we point out that the operatorAk is defined uniquely. (ii) Since B(v,w) is symmetric and satisfies the proposition of coercivity, Ak is symmetric positive definite woth respect to (·,·)k.The mesh-dependent norms are defined by(18)∥|v|∥s,k=(Aksv,v)k.
Observe that the energy norm ∥·∥E coincides with ∥|·|∥1,k norm on Vk. Similarly,  ∥|·|∥0,k  is the norm associated with the mesh-dependent inner product (14). The following lemma shows that ∥|·|∥0,k is equivalent to the L2 norm.Lemma 9.
The norm∥·∥L2(Ω) is equivalent to mesh-depentent norm ∥|·|∥0,k, that is,
(19)ρ∥|v|∥0,k≤∥v∥L2(Ω)≤∥|v|∥0,k,∀v∈Vk.
Particularly,  ∥|·|∥0,k=∥·∥L2(Ω)  for the uniform mesh.Proof.
The proof is analogous to Lemma6.2.7 in reference [18].In order to estimate the spectral radius,Λ(Ak), of Ak, we need the following norm interpolation lemma (see [19, Theorem 1.3.7]).Lemma 10 (interpolation theorem).
Assume thatΩ  is an open set of  ℝd with a Lipchitz continuous boundary. Let s1<s2 be two real numbers, and  μ=(1-θ)s1+θs2,0≤θ≤1. Then there exists a constant C6>0 such that
(20)∥u∥Hμ(Ω)≤C6∥u∥Hs1(Ω)1-θ∥u∥Hs2(Ω)θ.Lemma 11 (spectral radius theorem).
We have the estimation for spectral radiusΛ(Ak):
(21)Λ(Ak)≤C7hk-2α,
where C7 is a positive constant independent of k.Proof.
Letϕ∈Vk. Applying the continuity (6) of  B(·,·), norm interpolation Lemma 10, and inverse estimation, we have
(22)B(ϕ,ϕ)≤C2∥ϕ∥Hα(Ω)2≤C2(C6∥ϕ∥H0(Ω)1-α·∥ϕ∥H1(Ω)α)2≤C2(C6∥ϕ∥H0(Ω)1-α·C8hk-α∥ϕ∥H0(Ω)α)2=Chk-2α∥ϕ∥H0(Ω)2=Chk-2α∥ϕ∥L2(Ω)2,
where C8>0 is the constant of inverse estimate and C=C2C62C82.
Letλ be an eigenvalue of Ak with corresponding eigenvector ϕ∈Vk. From Lemma 9, we have
(23)B(ϕ,ϕ)=(Akϕ,ϕ)k=λ(ϕ,ϕ)k=λ∥|ϕ|∥0,k2=λ∥ϕ∥L2(Ω)2.
Combining (22) and (23), the following inequality holds:
(24)λ=B(ϕ,ϕ)∥ϕ∥L2(Ω)2≤Chk-2α∥ϕ∥L2(Ω)2∥ϕ∥L2(Ω)2=C7hk-2α,
where C7=C2C62C82.


## 3.2. The V-Cycle Multigrid Algorithm

Firstly, we give the intergrid transfer operators which play a very important role in convergence analysis.Definition 12 (intergrid transfer operator).
The coarse-to-fine intergrid transfer operatorIk-1k:Vk-1→Vk is taken to be the natural injection, that is,
(25)Ik-1kv=v,∀v∈Vk-1.
The fine-to-coarse operator Ikk-1:Vk→Vk-1 is defined by
(26)(Ikk-1w,v)k-1=(w,Ik-1kv)k=(w,v)k,∀v∈Vk-1,w∈Vk.Now, we describe thekth level of  V-cycle MGM (V-MGM) and full V-cycle MGM. Let m be the smoothing number, nk  the iteration number of kth level V-MGM, and Λk be a parameter dependent on k. Denote MG(k,z0,g) the approximate solution of the Akz=g obtained from the kth level iteration with initial guess z0. The discussion of parameters Λk and nk is left to the next subsection and Section 4.Algorithm 13 (thekth level V-cycle multigrid algorithm).
BEGIN:
Fork=1:  MG(1,z0,g)=A1-1g.
Fork>1, MG(k,z0,g) is obtained recursively in three steps.
Presmoothing Step. For1≤l≤m, let zl=zl-1+(1/Λk)(g-Akzl-1).
Error Correction.
g
¯
=
I
k
k
-
1
(
g
-
A
k
z
m
),
q
1
=
M
G
(
k
-
1,0
,
g
¯
),
z
m
+
1
=
z
m
+
I
k
-
1
k
q
1.
Postsmoothing step. Form+1≤l≤2m+1, let zl=zl-1+(1/Λk)(g-Akzl-1).
The output of thekth level iteration is  MG(k,z0,g):=z2m+1.
END.Algorithm 14 (the full V-cycle multigrid algorithm).
BEGIN:
Fork=1, u^1=A1-1f1.
Fork≥2,
u
0
k
=
I
k
-
1
k
u
^
k
-
1,
u
l
k
=
M
G
(
k
,
u
l
-
1
k
,
f
k
)
,
1
≤
l
≤
n
k,
u
^
k
=
u
n
k
k.
END.


## 3.3. Approximation and Smoothing Properties

In this subsection, we prove some properties of projection operatorPk-1 and smoothing operator Rk, which are the key ingredients for V-MGM and FV-MGM algorithm.Definition 15.
LetPk:V→Vk be the orthogonal projection with respect to B(·,·), that is,
(27)B(v,w)=B(Pkv,w),∀w∈Vk.Definition 16.
Define the relaxation operator:(28)Rk:=I-1ΛkAk.Throughout this paper, we assume thatΛk denotes some upper bound for the spectral radius of Ak satisfying  Λk≤C7hk-2α.Lemma 17.
Letv∈Vk, g¯=Ikk-1Akv,  and q∈Vk-1 satisfy the equation Ak-1q=g¯. Then q=Pk-1v.Proof.
Forw∈Vk-1,
(29)B(q,w)=(Ak-1q,w)k-1=(g¯,w)k-1=(Ikk-1Akv,w)k-1=(Akv,w)k=B(v,w).
Therefore, from the definition of Pk-1 we have q=Pk-1v.Remark 18.
Letg=Akv. From the Riesz representation theorem there exists a unique f∈H-α(Ω) such that
(30)(g,w)k=〈f,w〉,∀w∈Vk,
which means that v is a finite element solution of variational problem (4) on Vk. Observing that
(31)(g¯,w)k-1=(Ikk-1g,w)k-1=(g,w)k=〈f,w〉,∀w∈Vk-1,
we claim that q=Pk-1v is also a finite element solution of variational problem (4) on Vk-1. Therefore, noting  Vk-1⊂Vk (since 𝒯k-1⊂𝒯k) and applying Theorem 4, we have
(32)∥v-Pk-1v∥L2(Ω)≤C4hkα∥v-Pk-1v∥E,α≠34,∥v-Pk-1v∥L2(Ω)≤C4′hkα-ε∥v-Pk-1v∥E,α=34.Lemma 19.
There exists a positive constantC9  such that
(33)∥|(I-Pk-1)v|∥0,k≤C9hkα∥|(I-Pk-1)v|∥1,k,∀v∈Vk,α≠34,∥|(I-Pk-1)v|∥0,k≤C9hkα-ε∥|(I-Pk-1)v|∥1,k,∀v∈Vk,α=34,0<ε<12.Proof.
We only prove the caseα≠3/4, and the proof for the other case is analogous. Applying Remark 18 and Lemma 9 we have
(34)∥|(I-Pk-1)v|∥0,k≤1ρ∥(I-Pk-1)v∥L2(Ω)≤C4ρhkα∥(v-Pk-1v)∥E=C9hkα∥|(I-Pk-1)v|∥1,k,
where C9=C4/ρ.Lemma 20.
There exists a positive constantC9 such that
(35)∥|(I-Pk-1)v|∥1,k≤C9hkα∥|v|∥2,k,∀v∈Vk,α≠34,(36)∥|(I-Pk-1)v|∥1,k≤C9hkα-ε∥|v|∥2,k,∀v∈Vk,α=34,0<ε<12.Proof.
Applying projection operatorPk-1, generalized Cauchy-Schwarz inequality (see [18, Lemma 6.2.10]), and Lemma 19(37)∥|(I-Pk-1)v|∥1,k2=∥(I-Pk-1)v∥E2=B(v-Pk-1v,v-Pk-1v)=B(v-Pk-1v,v)≤∥|(I-Pk-1)v|∥0,k∥|v|∥2,k≤C9hkα∥|(I-Pk-1)v|∥1,k∥|v|∥2,k=C9hkα∥|(I-Pk-1)v|∥1,k∥|v|∥2,k
if α≠3/4. Therefore, dividing through by ∥|(I-Pk-1)v|∥1,k yields (35).  The case α=3/4 is analogous.Lemma 21.
Letm be the number of smoothing steps. Then
(38)B((I-Rk)Rk2mv,v)≤12mB((I-Rk2m)v,v).Proof.
See Lemma6.6.7 in [18].Lemma 22.
Letm be the number of smoothing steps. Then there exists a positive constant C10 such that
(39)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)≤C10mB((I-Rk2m)v,v).Proof.
Applying Lemmas20, 11, and 21, we have
(40)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)=∥(I-Pk-1)Rkmv∥E2=∥|(I-Pk-1)Rkmv|∥1,k2≤C92hk2α∥|Rkmv|∥2,k2=C92hk2α(Ak2Rkmv,Rkmv)k=C92hk2αB(AkRkmv,Rkmv)=C92hk2αΛkB((I-Rk)Rkmv,Rkmv)≤C92C7B((I-Rk)Rkmv,Rkmv)=C92C7B((I-Rk)Rk2mv,v)≤C92C7(2m)B((I-Rk2m)v,v)=C10mB((I-Rk2m)v,v)
if α≠3/4, where C10=C92C7/2. Note that (39) also holds for α=3/4 since ε>0.


## 3.4. Convergence Analysis

In this subsection, we prove convergence theorems for V-cycle FE multigrid method.Definition 23.
The error operator of V-cycle multigridEk:Vk→Vk is defined recursively by
(41)E1=0,Ek=Rkm[I-(I-Ek-1)Pk-1]Rkm,k>1.Proposition 24.
The operatorEk has the following propositions (see [18, Lemma 6.6.2 and Proposition 6.6.9]):(i)
if  z,g∈Vk  satisfy  Akz=g, then(42)Ek(z-z0)=z-MG(k,z0,g),(ii)
Ek  is symmetric positive semidefinite with respect toB(·,·)  for  k≥1, that is,(43)B(Ekv,w)=B(v,Ekw),∀v,w∈Vk,B(Ekv,v)≥0.Theorem 25.
Letm be the number of smoothing steps. Then
(44)B(Ekv,v)≤C10m+C10B(v,v),∀v∈Vk.Proof.
The proof is by induction. Fork=1, (44) is trivially true since Ek=0. Assume that (44) holds for k-1. Let γ=C10/(m+C10); therefore 1-γ=γm/C10. Applying induction hypothesis, Proposition 24, and Lemma 22, we have
(45)B(Ekv,v)=B(Rkmv,Rkmv)-B(Pk-1Rkmv,Pk-1Rkmv)+B(Ek-1Pk-1Rkmv,Pk-1Rkmv)≤B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB(Pk-1Rkmv,Pk-1Rkmv)=(1-γ)B((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB((I-Pk-1)Rkmv,(I-Pk-1)Rkmv)+γB(Pk-1Rkmv,Pk-1Rkmv)≤(1-γ)C10mB((I-Rk2m)v,v)+γB(Rkmv,Rkmv)=γB((I-Rk2m)v,v)+γB(Rkmv,Rkmv)=γB(v,v),
which ends the proof.Corollary 26.
Letm be the number of smoothing steps. Then one has the estimation for spectral radius Λ(Ek):
(46)Λ(Ek)≤C10m+C10.Theorem 27 (convergence of thekth level iteration for V-cycle).
Letm be the number of smoothing steps. Then
(47)∥z-MG(k,z0,g)∥E≤C10m+C10∥z-z0∥E.
Hence, the kth level iteration for any m is a contraction with the contraction number independent of k.Proof.
From (42) in Proposition 24, it suffices to show
(48)∥Ekv∥E≤C10m+C10∥v∥E,
which can be proved by Corollary 26.Theorem 28 (full multigrid convergence).
Letu∈H0α(Ω)⋂‍Hr(Ω)(α≤r≤2) solve (4). If the iteration numbers in full multigrid algorithm nk  are  large enough, there exists a positive constant C11 such that
(49)∥uk-u^k∥E≤C11hkr-α∥u∥Hr(Ω).Proof.
For the simpleness we assume thatnk≥n≥1. Define e^k=uk-u^k. In particular, e^1=0. Let γ=C10/(m+C10) (C10 is the constant in Theorem 27). We have
(50)∥e^k∥E≤γn∥uk-u^k-1∥E≤γn(∥uk-u∥E+∥u-uk-1∥E+∥uk-1-u^k-1∥E)≤γn[C2C3hkr-α∥u∥Hr(Ω)+C2C3hk-1r-α∥u∥Hr(Ω)+∥e^k-1∥E]=γn[3C2C3hkr-α∥u∥Hr(Ω)+∥e^k-1∥E].
In the following analysis, we assume that  2γn<1 (therefore n>logγ1/2). By iterating the above inequality and from continuity (6) and Theorem 3, we have
(51)∥e^k∥E≤3C2C3γnhkr-α∥u∥Hr(Ω)+3C2C3γ2nhk-1r-α∥u∥Hr(Ω)+⋯+3C2C3γknh1r-α∥u∥Hr(Ω)≤3C2C3γn1-2γnhkr-α∥u∥Hr(Ω)=C11hkr-α∥u∥Hr(Ω),
where C11=3C2C3γn/(1-2γn).The following Corollary is a natural conclusion of Theorems5 and 28.Corollary 29.
Under the assumptions in Theorem28, the convergence rate of multigrid solution u^k in L2 norm is r;  that is, there exists a constant C12>0 such that
(52)∥u^k-u∥L2(Ω)≤C12hkr∥u∥Hr(Ω),α≠34,∥u^k-u∥L2(Ω)≤C12hkr-ε∥u∥Hr(Ω),α=34,0<∀ε<12.We end this section with a proposition that the work involved in the full multigrid algorithm is𝒪(dk), where dk=Nk-2=dimVk (see [18, Proposition 6.7.4]).


## 4. Numerical Examples

In this section we employ the FV-MGM listed as Algorithm14 in Section 3 to solve FADEs, which demonstrate the convergence rate and involved work in Theorem 28 and its Corollary 29. The parameter m is taken 5, Λk=C7hk-2α, and nk is controlled by the stopping criterion:
(53)∥Aku^nkk-g∥L2(Ω)≤𝒪(hkα).
We note that the positive constants C7 can be estimated by testing examples on the 1th and 2th level meshes, that is,
(54)C7=12(h12αΛ(A1)+h22αΛ(A2)).
Here we calculate C7=0.05, and, from experimental experience, let the right side of stopping criterion be 0.085hkα. All numerical experiments are run in MATLAB 7.0.0.19920(R14) on a PC with the configuration: Intel(R)Core(TM)i3-2100 CPU @ 3.2 GHz and 1.88 GB RAM.Example 30.
Letc(x)=1. It can be verified that u(x)=x2-x3 is the exact solution to the boundary value problem:
(55)-12D(D0x-β+Dx1-β)Du+u=f,u(0)=0,u(1)=0,
where
(56)f(x)=x2-x3+12f1(x)+12f2(x),f1(x)=-2xβΓ(1+β)+6x1+βΓ(2+β),f2(x)=-β(1-x)β-1Γ(1+β)+4(β+1)(1-x)βΓ(2+β)-6(β+2)(1-x)β+1Γ(3+β).Asu∈H02(Ω), Corollary 29 predicts a rate of convergence of 2 in L2 norm. Table 1 includes numerical results over a uniform partition of [0,1], which support the predicted rates of convergence for different values of β.Table 1
Experimental error∥u-u^k∥L2(Ω) for different β.
h
k
β
=
0.3
β
=
0.6
β
=
0.9
Error
Conv. rate
Error
Conv. rate
Error
Conv. rate
1/4
9.3587
e
-
3
8.3460
e
-
3
7.6044
e
-
3
1/8
2.1751
e
-
3
2.1388
1.8444
e
-
3
2.1779
1.6348
e
-
3
2.2177
1/16
5.0679
e
-
4
2.1016
4.1362
e
-
4
2.1568
3.6349
e
-
4
2.1692
1/32
1.1967
e
-
4
2.0823
9.4937
e
-
5
2.1233
8.3872
e
-
5
2.1157
1/64
2.9496
e
-
5
2.0205
2.2629
e
-
5
2.0688
2.0343
e
-
5
2.0436
1/128
7.0818
e
-
6
2.0583
5.4727
e
-
6
2.0478
5.0852
e
-
6
2.0002As comparisons, we also carry out the Gauss elimination (GE) and conjugate gradient normal residual (CGNR) method with the same stopping criterion of the FV-MGM, that is,(57)∥Aku^nkk-g∥L2(Ω)≤10-n,
to solve the corresponding system. For escaping “out of memory,” we define that the data type of Ak is short float in our program. Table 2 includes CPU time (without including stiffness matrix calculating time) and iteration numbers for each of the numerical methods with n=4. The rate of the increasing CPU time is defined by(58)Rcpu=log[CPUtimeonfinermesh/CPUtimeoncoarsermesh]log[dimoffinermesh/dimofcoarsermesh].From the table, we see that the numbers of iterations by our FV-MGM are independent of hk. In contrast, the CGNR method (with the initial value z0=0) needs more iterations when hk decrease. Similarly, the CPU time needed for GE and CGNR method increases much faster than that of the FV-MGM.Table 2
Comparisons for solving Example30 with  β=0.7  by the GE method, the CGNR method, and the FV-MGM, respectively.
h
k
FVMGM
GE
CGNR
FVMGM
Error
CPU(s)
R
cpu
Iter
CPU(s)
R
cpu
Iter =nk
CPU(s)
R
cpu
1/128
2.0270
e
-
5
0.0313
29
0.0000
4
0.0000
1/256
1.3811
e
-
5
0.0625
0.995
46
0.0156
5
0.0156
1/512
1.1507
e
-
5
0.1719
1.451
72
0.0313
0.999
5
0.0469
1.579
1/1024
1.1101
e
-
5
0.5156
1.580
112
0.2187
2.796
3
0.1250
1.417
1/2048
1.1002
e
-
5
2.5000
2.274
174
1.5625
2.833
1
0.2688
1.103
1/4096
1.0849
e
-
5
15.484
2.627
270
9.2031
2.554
1
0.5719
1.088


## 5. Concluding Remarks

In general the discretized system of FADE has a full and ill-conditioned coefficient matrix, so FV-MGM is a high-efficient algorithm for solving these equations. Theorems and examples in this paper show that the convergence rate of FV-MGM is the same as classic FEM under the stopping criterion (53), and the computational work is only 𝒪(dimVk) while the stopping criterion is taken (57). Different from integer-order equations, all the convergence analysis is on fractional Sobolev spaces Hα(Ω). The nonsymmetric form of FADE will be studied in the future.



---

*Source: 385463-2013-06-04.xml*