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

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 both 𝐿² norm and energy norm for FV-MGM. Numerical examples are given to demonstrate the convergence rate and efficiency of the method.

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),𝑓∶Ω→ℝ, find 𝑢∶Ω→ℝ such that −12𝐷₀𝐷_𝑥^−𝛽+_𝑥𝐷₁^−𝛽𝐷𝑢+𝑐(𝑥)𝑢=𝑓(𝑥),inΩ,𝑢=0,on𝜕Ω, (1) where 0<𝛽<1, 𝑐(𝑥)≥0, 𝐷 represents a single spatial derivative, and ₀𝐷_𝑥^−𝛽 and _𝑥𝐷₁^−𝛽 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 𝒪(𝑁³) or 𝒪(𝑁log²𝑁). 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 Section 3 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 𝐵∶𝐻^𝛼₀(Ω)×𝐻^𝛼₀(Ω)→ℝ as 𝐵()1𝑢,𝑣∶=2⁰𝐷_𝑥^−𝛽+1𝐷𝑢,𝐷𝑣2^𝑥𝐷₁^−𝛽𝐷𝑢,𝐷𝑣+(𝑐(𝑥)𝑢,𝑣), (2) and linear functional 𝐹∶𝐻^𝛼₀(Ω)→ℝ as 𝐹(𝑣)∶=𝑓,𝑣, (3) where (⋅,⋅) denotes the inner product on 𝐿²(Ω) and ⟨⋅,⋅⟩ the duality pairing of 𝐻^−𝜇(Ω) and 𝐻^𝜇₀(Ω),𝜇≥0.

Thus, the Galerkin variational solution of Problem 1 may be defined as follows.

Definition 2 (variational solution). A function 𝑢∈𝐻^𝛼₀(Ω) is a variational solution of Problem 1 provided that 𝐵(𝑢,𝑣)=𝐹(𝑣),∀𝑣∈𝐻^𝛼₀(Ω). (4)

Ervin and Roop ([17], 2006) proved that the bilinear form 𝐵(⋅,⋅) satisfies the coercivity and continuity; that is, there exist two positive constants 𝐶₁ and 𝐶₂ such that 𝐵(𝑢,𝑢)≥𝐶₁‖𝑢‖²_𝐻^𝛼_(Ω),∀𝑢∈𝐻^𝛼₀(Ω),𝐵(𝑢,𝑣)≤𝐶₂‖𝑢‖_𝐻^𝛼_(Ω)‖𝑣‖_𝐻^𝛼_(Ω),∀𝑢,𝑣∈𝐻^𝛼₀(Ω). (5)(6) Hence, they claimed that there exists a unique solution to (4). Note that from (5) and (6) we have norm equivalence of ‖⋅‖_𝐻^𝛼_(Ω) and energy norm ‖⋅‖_𝐸=𝐵(⋅,⋅)^1/2, that is, 𝐶₁‖𝑣‖_𝐻^𝛼_(Ω)≤‖𝑣‖_𝐸≤𝐶₂‖𝑣‖_𝐻^𝛼_(Ω),∀𝑣∈𝐻^𝛼₀(Ω). (7)

Let 𝑆_ℎ denote a partition of Ω such that Ω={⋃𝐾∶𝐾∈𝑆_ℎ}. Assume that there exists a positive constant 𝜌 such that 𝜌ℎ_𝑘≤ℎ_𝐾≤ℎ_𝑘, where ℎ_𝑘=max_{𝐾∈𝑆ℎ}ℎ_𝐾. Let 𝑃_𝑚(𝐾) denote the space of polynomials of degree less than or equal to 𝑚 on 𝐾∈𝑆_ℎ. Associated with 𝑆_ℎ, define the finite-dimensional subspace 𝑋_ℎ⊂𝐻^𝛼₀(Ω) as 𝑋_ℎ∶=𝑣∈𝐻^𝛼₀(Ω)∩𝐶⁰Ω∶𝑣∣_𝐾∈𝑃_𝑚−1(𝐾),∀𝐾∈𝑆_ℎ. (8) Let 𝑢_ℎ be the solution to the finite-dimensional variational problem: 𝐵𝑢_ℎ,𝑣_ℎ=𝐹𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ. (9)

Note that the existence and uniqueness of solution to (9) follow from the fact that 𝑋_ℎ is a subset of the space 𝐻^𝛼₀(Ω)∩𝐻^𝑚₀(Ω). Ervin and Roop have obtained convergence estimate in the energy norm ‖⋅‖_𝐸. 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 𝐿² norm was proved.

Theorem 3 (see [17, Corollary 4.3]). Let 𝑢∈𝐻^𝛼₀(Ω)⋂𝐻^𝑟(Ω)(𝛼≤𝑟≤𝑚) solve (4) and 𝑢_ℎ solve (9). Then there exists a constant 𝐶₃>0 such that 𝑒=𝑢−𝑢_ℎ satisfies ‖𝑒‖_𝐻^𝛼_(Ω)≤𝐶₃ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω). (10)

Theorem 4 (see [17, Theorem 4.4]). Let 𝑢∈𝐻^𝛼₀(Ω)⋂𝐻^𝑟(Ω)(𝛼≤𝑟≤𝑚) solve (4) and 𝑢_ℎ solve (9), where 𝑚−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 𝐶₄ (or 𝐶^′₄) such that the error 𝑒=𝑢−𝑢_ℎ satisfies ‖𝑒‖_𝐿²_(Ω)≤𝐶₄ℎ^𝛼_𝑘‖𝑢‖_𝐸3,𝛼≠4,‖𝑒‖_𝐿²_(Ω)≤𝐶^′₄ℎ_𝑘^𝛼−𝜀‖𝑢‖_𝐸3,𝛼=41,0<∀𝜀<2. (11)

From Theorems 3 and 4 and the equivalence of ‖⋅‖_𝐻^𝛼_(Ω) and ‖⋅‖_𝐸, we have the error estimates in the 𝐿² norm.

Theorem 5. Under the same assumptions in Theorem 4, one has ‖𝑒‖_𝐿²_(Ω)≤𝐶₅ℎ^𝑟_𝑘‖𝑢‖_𝐻^𝑟_(Ω)3,𝛼≠4,‖𝑒‖_𝐿²_(Ω)≤𝐶₅ℎ_𝑘^𝑟−𝜀‖𝑢‖_𝐻^𝑟_(Ω)3,𝛼=41,0<∀𝜀<2, (12) where 𝐶₅=𝐶₂𝐶₃𝐶₄ (or 𝐶₅=𝐶₂𝐶₃𝐶^′₄).

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 {𝒯_𝑘} denote a sequence of partitions of Ω and ℎ_𝑘 be the mesh size of 𝒯_𝑘. From now on, we assume that {𝒯_𝑘} is a quasiuniform family, that is, 𝜌ℎ_𝑘≤𝑥_𝑖−𝑥_𝑖−1≤ℎ_𝑘,𝑖=1,2,…,𝑁_𝑘, (13) with positive constant 𝜌. At the same time, let 𝒯_𝑘 be obtained from 𝒯_𝑘−1 via a regular subdivision, that is, 𝒯_𝑘−1⊂𝒯_𝑘 for 𝑘≥2. Let 𝑉_𝑘 denote 𝐶⁰ piecewise linear functions with respect to 𝒯_𝑘 that vanish on 𝜕Ω.

3.1. Mesh-Dependent Norms

Definition 6. The mesh-dependent inner product (⋅,⋅)_𝑘 on 𝑉_𝑘 is defined by (𝑣,𝑤)_𝑘=ℎ_𝑘^𝑁𝑘−1_𝑖=1𝑣𝑥_𝑖𝑤𝑥_𝑖. (14)

Definition 7. The operator 𝐴_𝑘∶𝑉_𝑘→𝑉_𝑘 is defined by 𝐴_𝑘𝑣,𝑤_𝑘=𝐵(𝑣,𝑤),∀𝑣,𝑤∈𝑉_𝑘. (15)

In terms of the operator 𝐴_𝑘, the discretized equation of (9) can be written as 𝐴_𝑘𝑢_𝑘=𝑓_𝑘, (16) where 𝑓_𝑘∈𝑉_𝑘 satisfies 𝑓_𝑘,𝑣_𝑘=𝑓,𝑣,∀𝑣∈𝑉_𝑘. (17)

Remark 8. (i) From the Riesz representation theorem, we point out that the operator 𝐴_𝑘 is defined uniquely. (ii) Since 𝐵(𝑣,𝑤) is symmetric and satisfies the proposition of coercivity, 𝐴_𝑘 is symmetric positive definite woth respect to (⋅,⋅)_𝑘.

The mesh-dependent norms are defined by ‖|𝑣|‖_𝑠,𝑘=𝐴^𝑠_𝑘𝑣,𝑣_𝑘. (18) Observe that the energy norm ‖⋅‖_𝐸 coincides with ‖|⋅|‖_1,𝑘 norm on 𝑉_𝑘. Similarly, ‖|⋅|‖_0,𝑘 is the norm associated with the mesh-dependent inner product (14). The following lemma shows that ‖|⋅|‖_0,𝑘 is equivalent to the 𝐿² norm.

Lemma 9. The norm ‖⋅‖_𝐿²_(Ω) is equivalent to mesh-depentent norm ‖|⋅|‖_0,𝑘, that is, 𝜌‖|𝑣|‖_0,𝑘≤‖𝑣‖_𝐿²_(Ω)≤‖|𝑣|‖_0,𝑘,∀𝑣∈𝑉_𝑘. (19) Particularly, ‖|⋅|‖_0,𝑘=‖⋅‖_𝐿²_(Ω) for the uniform mesh.

Proof. The proof is analogous to Lemma 6.2.7 in reference [18].

In order to estimate the spectral radius, Λ(𝐴_𝑘), of 𝐴_𝑘, we need the following norm interpolation lemma (see [19, Theorem 1.3.7]).

Lemma 10 (interpolation theorem). Assume that Ω is an open set of  ℝ^𝑑 with a Lipchitz continuous boundary. Let 𝑠₁<𝑠₂ be two real numbers, and 𝜇=(1−𝜃)𝑠₁+𝜃𝑠₂,0≤𝜃≤1. Then there exists a constant 𝐶₆>0 such that ‖𝑢‖_𝐻^𝜇_(Ω)≤𝐶₆‖𝑢‖_𝐻^1−𝜃𝑠_1(Ω)‖𝑢‖^𝜃_𝐻^𝑠_2(Ω). (20)

Lemma 11 (spectral radius theorem). We have the estimation for spectral radius Λ(𝐴_𝑘): Λ𝐴_𝑘≤𝐶₇ℎ_𝑘^−2𝛼, (21) where 𝐶₇ is a positive constant independent of 𝑘.

Proof. Let 𝜙∈𝑉_𝑘. Applying the continuity (6) of 𝐵(⋅,⋅), norm interpolation Lemma 10, and inverse estimation, we have 𝐵𝜙,𝜙≤𝐶₂𝜙²_𝐻𝛼(Ω)≤𝐶₂𝐶₆𝜙_𝐻^1−𝛼_0(Ω)⋅𝜙^𝛼_𝐻1(Ω)²≤𝐶₂𝐶₆𝜙_𝐻^1−𝛼_0(Ω)⋅𝐶₈ℎ_𝑘^−𝛼𝜙^𝛼_𝐻0(Ω)²=𝐶ℎ_𝑘^−2𝛼𝜙²_𝐻0(Ω)=𝐶ℎ_𝑘^−2𝛼𝜙²_𝐿2(Ω), (22) where 𝐶₈>0 is the constant of inverse estimate and 𝐶=𝐶₂𝐶²₆𝐶²₈.
Let 𝜆 be an eigenvalue of 𝐴_𝑘 with corresponding eigenvector 𝜙∈𝑉_𝑘. From Lemma 9, we have 𝐵𝜙,𝜙=𝐴_𝑘𝜙,𝜙_𝑘=𝜆𝜙,𝜙_𝑘=𝜆𝜙²_0,𝑘=𝜆𝜙²_𝐿2(Ω). (23) Combining (22) and (23), the following inequality holds: 𝐵𝜆=𝜙,𝜙𝜙²_𝐿2(Ω)≤𝐶ℎ_𝑘^−2𝛼𝜙²_𝐿2(Ω)𝜙²_𝐿2(Ω)=𝐶₇ℎ_𝑘^−2𝛼, (24) where 𝐶₇=𝐶₂𝐶²₆𝐶²₈.

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 operator 𝐼^𝑘_𝑘−1∶𝑉_𝑘−1→𝑉_𝑘 is taken to be the natural injection, that is, 𝐼^𝑘_𝑘−1𝑣=𝑣,∀𝑣∈𝑉_𝑘−1. (25) The fine-to-coarse operator 𝐼_𝑘^𝑘−1∶𝑉_𝑘→𝑉_𝑘−1 is defined by 𝐼_𝑘^𝑘−1𝑤,𝑣_𝑘−1=𝑤,𝐼^𝑘_𝑘−1𝑣_𝑘=(𝑤,𝑣)_𝑘,∀𝑣∈𝑉_𝑘−1,𝑤∈𝑉_𝑘. (26)

Now, we describe the 𝑘th level of V-cycle MGM (V-MGM) and full V-cycle MGM. Let 𝑚 be the smoothing number, 𝑛_𝑘 the iteration number of 𝑘th level V-MGM, and Λ_𝑘 be a parameter dependent on 𝑘. Denote 𝑀𝐺(𝑘,𝑧₀,𝑔) the approximate solution of the 𝐴_𝑘𝑧=𝑔 obtained from the 𝑘th level iteration with initial guess 𝑧₀. The discussion of parameters Λ_𝑘 and 𝑛_𝑘 is left to the next subsection and Section 4.

Algorithm 13 (the 𝑘th level V-cycle multigrid algorithm).  BEGIN: For 𝑘=1: 𝑀𝐺(1,𝑧₀,𝑔)=𝐴₁⁻¹𝑔. For 𝑘>1, 𝑀𝐺(𝑘,𝑧₀,𝑔) is obtained recursively in three steps. Presmoothing Step. For 1≤𝑙≤𝑚, let 𝑧_𝑙=𝑧_𝑙−1+(1/Λ_𝑘)(𝑔−𝐴_𝑘𝑧_𝑙−1). Error Correction.  𝑔=𝐼_𝑘^𝑘−1(𝑔−𝐴_𝑘𝑧_𝑚),  𝑞₁=𝑀𝐺(𝑘−1,0,𝑔),  𝑧_𝑚+1=𝑧_𝑚+𝐼^𝑘_𝑘−1𝑞₁.  Postsmoothing step. For 𝑚+1≤𝑙≤2𝑚+1, let 𝑧_𝑙=𝑧_𝑙−1+(1/Λ_𝑘)(𝑔−𝐴_𝑘𝑧_𝑙−1). The output of the 𝑘th level iteration is 𝑀𝐺(𝑘,𝑧₀,𝑔)∶=𝑧_2𝑚+1. END.

Algorithm 14 (the full V-cycle multigrid algorithm).  BEGIN: For 𝑘=1, 𝑢₁=𝐴₁⁻¹𝑓₁. For 𝑘≥2,  𝑢^𝑘₀=𝐼^𝑘_𝑘−1𝑢_𝑘−1,  𝑢^𝑘_𝑙=𝑀𝐺(𝑘,𝑢^𝑘_𝑙−1,𝑓_𝑘),1≤𝑙≤𝑛_𝑘,  𝑢_𝑘=𝑢^𝑘_𝑛𝑘. END.

3.3. Approximation and Smoothing Properties

In this subsection, we prove some properties of projection operator 𝑃_𝑘−1 and smoothing operator 𝑅_𝑘, which are the key ingredients for V-MGM and FV-MGM algorithm.

Definition 15. Let 𝑃_𝑘∶𝑉→𝑉_𝑘 be the orthogonal projection with respect to 𝐵(⋅,⋅), that is, 𝐵(𝑣,𝑤)=𝐵𝑃_𝑘𝑣,𝑤,∀𝑤∈𝑉_𝑘. (27)

Definition 16. Define the relaxation operator: 𝑅_𝑘1∶=𝐼−Λ_𝑘𝐴_𝑘. (28)

Throughout this paper, we assume that Λ_𝑘 denotes some upper bound for the spectral radius of 𝐴_𝑘 satisfying Λ_𝑘≤𝐶₇ℎ_𝑘^−2𝛼.

Lemma 17. Let 𝑣∈𝑉_𝑘, 𝑔=𝐼_𝑘^𝑘−1𝐴_𝑘𝑣, and 𝑞∈𝑉_𝑘−1 satisfy the equation 𝐴_𝑘−1𝑞=𝑔. Then 𝑞=𝑃_𝑘−1𝑣.

Proof. For 𝑤∈𝑉_𝑘−1, 𝐵𝑞,𝑤=𝐴_𝑘−1𝑞,𝑤_𝑘−1=𝑔,𝑤_𝑘−1=𝐼_𝑘^𝑘−1𝐴_𝑘𝑣,𝑤_𝑘−1=𝐴_𝑘𝑣,𝑤_𝑘=𝐵(𝑣,𝑤). (29) Therefore, from the definition of 𝑃_𝑘−1 we have 𝑞=𝑃_𝑘−1𝑣.

Remark 18. Let 𝑔=𝐴_𝑘𝑣. From the Riesz representation theorem there exists a unique 𝑓∈𝐻^−𝛼(Ω) such that 𝑔,𝑤_𝑘=𝑓,𝑤,∀𝑤∈𝑉_𝑘, (30) which means that 𝑣 is a finite element solution of variational problem (4) on 𝑉_𝑘. Observing that 𝑔,𝑤_𝑘−1=𝐼_𝑘^𝑘−1𝑔,𝑤_𝑘−1=𝑔,𝑤_𝑘=𝑓,𝑤,∀𝑤∈𝑉_𝑘−1, (31) we claim that 𝑞=𝑃_𝑘−1𝑣 is also a finite element solution of variational problem (4) on 𝑉_𝑘−1. Therefore, noting 𝑉_𝑘−1⊂𝑉_𝑘 (since 𝒯_𝑘−1⊂𝒯_𝑘) and applying Theorem 4, we have 𝑣−𝑃_𝑘−1𝑣_𝐿2(Ω)≤𝐶₄ℎ^𝛼_𝑘𝑣−𝑃_𝑘−1𝑣_𝐸3,𝛼≠4,𝑣−𝑃_𝑘−1𝑣_𝐿2(Ω)≤𝐶^′₄ℎ_𝑘^𝛼−𝜀𝑣−𝑃_𝑘−1𝑣_𝐸3,𝛼=4. (32)

Lemma 19. There exists a positive constant 𝐶₉ such that 𝐼−𝑃_𝑘−1𝑣_0,𝑘≤𝐶₉ℎ^𝛼_𝑘𝐼−𝑃_𝑘−1𝑣_1,𝑘,∀𝑣∈𝑉_𝑘3,𝛼≠4,𝐼−𝑃_𝑘−1𝑣_0,𝑘≤𝐶₉ℎ_𝑘^𝛼−𝜀𝐼−𝑃_k−1𝑣_1,𝑘,∀𝑣∈𝑉_𝑘3,𝛼=41,0<𝜀<2. (33)

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 𝐼−𝑃_𝑘−1𝑣_0,𝑘≤1𝜌𝐼−𝑃_𝑘−1𝑣_𝐿2(Ω)≤𝐶₄𝜌ℎ^𝛼_𝑘𝑣−𝑃_𝑘−1𝑣_𝐸=𝐶₉ℎ^𝛼_𝑘𝐼−𝑃_𝑘−1𝑣_1,𝑘, (34) where 𝐶₉=𝐶₄/𝜌.

Lemma 20. There exists a positive constant 𝐶₉ such that 𝐼−𝑃_𝑘−1𝑣_1,𝑘≤𝐶₉ℎ^𝛼_𝑘‖|𝑣|‖_2,𝑘,∀𝑣∈𝑉_𝑘3,𝛼≠4,𝐼−𝑃_𝑘−1𝑣_1,𝑘≤𝐶₉ℎ_𝑘^𝛼−𝜀‖|𝑣|‖_2,𝑘,∀𝑣∈𝑉_𝑘3,𝛼=41,0<𝜀<2. (35)(36)

Proof. Applying projection operator 𝑃_𝑘−1, generalized Cauchy-Schwarz inequality (see [18, Lemma 6.2.10]), and Lemma 19(𝐼−𝑃_𝑘−1)𝑣²_1,𝑘=(𝐼−𝑃_𝑘−1)𝑣²_𝐸=𝐵𝑣−𝑃_𝑘−1𝑣,𝑣−𝑃_𝑘−1𝑣=𝐵𝑣−𝑃_𝑘−1≤𝑣,𝑣𝐼−𝑃_𝑘−1𝑣_0,𝑘‖|𝑣|‖_2,𝑘≤𝐶₉ℎ^𝛼_𝑘𝐼−𝑃_𝑘−1𝑣_1,𝑘‖|𝑣|‖_2,𝑘=𝐶₉ℎ^𝛼_𝑘𝐼−𝑃_𝑘−1𝑣_1,𝑘‖|𝑣|‖_2,𝑘 (37) if 𝛼≠3/4. Therefore, dividing through by ‖|(𝐼−𝑃_𝑘−1)𝑣|‖_1,𝑘 yields (35). The case 𝛼=3/4 is analogous.

Lemma 21. Let 𝑚 be the number of smoothing steps. Then 𝐵𝐼−𝑅_𝑘𝑅_𝑘^2𝑚≤1𝑣,𝑣𝐵2𝑚𝐼−𝑅_𝑘^2𝑚𝑣,𝑣. (38)

Proof. See Lemma 6.6.7 in [18].

Lemma 22. Let 𝑚 be the number of smoothing steps. Then there exists a positive constant 𝐶₁₀ such that 𝐵𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣≤𝐶₁₀𝑚𝐵𝐼−𝑅_𝑘^2𝑚𝑣,𝑣. (39)

Proof. Applying Lemmas 20, 11, and 21, we have 𝐵𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣=(𝐼−𝑃_𝑘−1)𝑅^𝑚_𝑘𝑣²_𝐸=𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣²_1,𝑘≤𝐶²₉ℎ_𝑘^2𝛼𝑅^𝑚_𝑘𝑣²_2,𝑘=𝐶²₉ℎ_𝑘^2𝛼𝐴²_𝑘𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣_𝑘=𝐶²₉ℎ_𝑘^2𝛼𝐵𝐴_𝑘𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣=𝐶²₉ℎ_𝑘^2𝛼Λ_𝑘𝐵𝐼−𝑅_𝑘𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣≤𝐶²₉𝐶₇𝐵𝐼−𝑅_𝑘𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣=𝐶²₉𝐶₇𝐵𝐼−𝑅_𝑘𝑅_𝑘^2𝑚≤𝐶𝑣,𝑣₂₉𝐶₇(2𝑚)𝐵𝐼−𝑅_𝑘^2𝑚=𝐶𝑣,𝑣₁₀𝑚𝐵𝐼−𝑅_𝑘^2𝑚𝑣,𝑣 (40) if 𝛼≠3/4, where 𝐶₁₀=𝐶²₉𝐶₇/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 multigrid 𝐸_𝑘∶𝑉_𝑘→𝑉_𝑘 is defined recursively by 𝐸₁𝐸=0,_𝑘=𝑅^𝑚_𝑘𝐼−𝐼−𝐸_𝑘−1𝑃_𝑘−1𝑅^𝑚_𝑘,𝑘>1. (41)

Proposition 24. The operator 𝐸_𝑘 has the following propositions (see [18, Lemma 6.6.2 and Proposition 6.6.9]):(i)if  𝑧,𝑔∈𝑉_𝑘  satisfy  𝐴_𝑘𝑧=𝑔, then𝐸_𝑘𝑧−𝑧₀=𝑧−𝑀𝐺𝑘,𝑧₀,𝑔, (42)(ii)  𝐸_𝑘  is symmetric positive semidefinite with respect to 𝐵(⋅,⋅)  for  𝑘≥1, that is,𝐵𝐸_𝑘𝑣,𝑤=𝐵𝑣,𝐸_𝑘𝑤,∀𝑣,𝑤∈𝑉_𝑘,𝐵𝐸_𝑘𝑣,𝑣≥0. (43)

Theorem 25. Let 𝑚 be the number of smoothing steps. Then 𝐵𝐸_𝑘≤𝐶𝑣,𝑣¹⁰𝑚+𝐶₁₀𝐵(𝑣,𝑣),∀𝑣∈𝑉_𝑘. (44)

Proof. The proof is by induction. For 𝑘=1, (44) is trivially true since 𝐸_𝑘=0. Assume that (44) holds for 𝑘−1. Let 𝛾=𝐶₁₀/(𝑚+𝐶₁₀); therefore 1−𝛾=𝛾𝑚/𝐶₁₀. Applying induction hypothesis, Proposition 24, and Lemma 22, we have 𝐵𝐸_𝑘𝑣,𝑣=𝐵𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣−𝐵𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝑃_𝑘−1𝑅^𝑚_𝑘𝑣+𝐵𝐸_𝑘−1𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝑃_𝑘−1𝑅^𝑚_𝑘𝑣≤𝐵𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣+𝛾𝐵𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝑃_𝑘−1𝑅^𝑚_𝑘𝑣=1−𝛾𝐵𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣+𝛾𝐵𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝐼−𝑃_𝑘−1𝑅^𝑚_𝑘𝑣+𝛾𝐵𝑃_𝑘−1𝑅^𝑚_𝑘𝑣,𝑃_𝑘−1𝑅^𝑚_𝑘𝑣≤1−𝛾𝐶₁₀𝑚𝐵𝐼−𝑅_𝑘^2𝑚𝑅𝑣,𝑣+𝛾𝐵^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣=𝛾𝐵𝐼−𝑅_𝑘^2𝑚𝑣,𝑣+𝛾𝐵𝑅^𝑚_𝑘𝑣,𝑅^𝑚_𝑘𝑣=𝛾𝐵(𝑣,𝑣), (45) which ends the proof.

Corollary 26. Let 𝑚 be the number of smoothing steps. Then one has the estimation for spectral radius Λ(𝐸_𝑘): Λ𝐸_𝑘≤𝐶₁₀𝑚+𝐶₁₀. (46)

Theorem 27 (convergence of the 𝑘th level iteration for V-cycle). Let 𝑚 be the number of smoothing steps. Then 𝑧−𝑀𝐺𝑘,𝑧₀,𝑔_𝐸≤𝐶₁₀𝑚+𝐶₁₀𝑧−𝑧_0𝐸. (47) Hence, the 𝑘th level iteration for any 𝑚 is a contraction with the contraction number independent of 𝑘.

Proof. From (42) in Proposition 24, it suffices to show 𝐸_𝑘𝑣_𝐸≤𝐶₁₀𝑚+𝐶₁₀‖𝑣‖_𝐸, (48) which can be proved by Corollary 26.

Theorem 28 (full multigrid convergence). Let 𝑢∈𝐻^𝛼₀(Ω)⋂𝐻^𝑟(Ω)(𝛼≤𝑟≤2) solve (4). If the iteration numbers in full multigrid algorithm 𝑛_𝑘 are large enough, there exists a positive constant 𝐶₁₁ such that 𝑢_𝑘−𝑢_𝑘𝐸≤𝐶₁₁ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω). (49)

Proof. For the simpleness we assume that 𝑛_𝑘≥𝑛≥1. Define 𝑒_𝑘̂=𝑢_𝑘−𝑢_𝑘. In particular, 𝑒₁̂=0. Let 𝛾=𝐶₁₀/(𝑚+𝐶₁₀) (𝐶₁₀ is the constant in Theorem 27). We have 𝑒̂_𝑘𝐸≤𝛾^𝑛𝑢_𝑘−𝑢_𝑘−1𝐸≤𝛾^𝑛𝑢_𝑘−𝑢_𝐸+𝑢−𝑢_𝑘−1𝐸+𝑢_𝑘−1−𝑢_𝑘−1𝐸≤𝛾^𝑛𝐶₂𝐶₃ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω)+𝐶₂𝐶₃ℎ^𝑟−𝛼_𝑘−1‖𝑢‖_𝐻^𝑟_(Ω)+𝑒̂_𝑘−1𝐸=𝛾^𝑛3𝐶₂𝐶₃ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω)+𝑒̂_𝑘−1𝐸. (50) In the following analysis, we assume that 2𝛾^𝑛<1 (therefore 𝑛>log_𝛾^1/2). By iterating the above inequality and from continuity (6) and Theorem 3, we have 𝑒̂_𝑘𝐸≤3𝐶₂𝐶₃𝛾^𝑛ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω)+3𝐶₂𝐶₃𝛾^2𝑛ℎ^𝑟−𝛼_𝑘−1‖𝑢‖_𝐻^𝑟_(Ω)+⋯+3𝐶₂𝐶₃𝛾^𝑘𝑛ℎ₁^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω)≤3𝐶₂𝐶₃𝛾^𝑛1−2𝛾^𝑛ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω)=𝐶₁₁ℎ_𝑘^𝑟−𝛼‖𝑢‖_𝐻^𝑟_(Ω), (51) where 𝐶₁₁=3𝐶₂𝐶₃𝛾^𝑛/(1−2𝛾^𝑛).

The following Corollary is a natural conclusion of Theorems 5 and 28.

Corollary 29. Under the assumptions in Theorem 28, the convergence rate of multigrid solution 𝑢_𝑘 in 𝐿² norm is 𝑟; that is, there exists a constant 𝐶₁₂>0 such that 𝑢_𝑘−𝑢_𝐿2(Ω)≤𝐶₁₂ℎ^𝑟_𝑘‖𝑢‖_𝐻^𝑟_(Ω)3,𝛼≠4,𝑢_𝑘−𝑢_𝐿2(Ω)≤𝐶₁₂ℎ_𝑘^𝑟−𝜀‖𝑢‖_𝐻^𝑟_(Ω)3,𝛼=41,0<∀𝜀<2. (52)

We end this section with a proposition that the work involved in the full multigrid algorithm is 𝒪(𝑑_𝑘), where 𝑑_𝑘=𝑁_𝑘−2=dim𝑉_𝑘 (see [18, Proposition 6.7.4]).

4. Numerical Examples

In this section we employ the FV-MGM listed as Algorithm 14 in Section 3 to solve FADEs, which demonstrate the convergence rate and involved work in Theorem 28 and its Corollary 29. The parameter 𝑚 is taken 5, Λ_𝑘=𝐶₇ℎ_𝑘^−2𝛼, and 𝑛_𝑘 is controlled by the stopping criterion: 𝐴_𝑘𝑢^𝑘_𝑛𝑘−𝑔_𝐿2(Ω)≤𝒪ℎ^𝛼_𝑘. (53) We note that the positive constants 𝐶₇ can be estimated by testing examples on the 1th and 2th level meshes, that is, 𝐶₇=12ℎ₁^2𝛼Λ𝐴₁+ℎ₂^2𝛼Λ𝐴₂. (54) Here we calculate 𝐶₇=0.05, and, from experimental experience, let the right side of stopping criterion be 0.085ℎ^𝛼_𝑘. 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. Let 𝑐(𝑥)=1. It can be verified that 𝑢(𝑥)=𝑥²−𝑥³ is the exact solution to the boundary value problem: −12𝐷₀𝐷_𝑥^−𝛽+_𝑥𝐷₁^−𝛽𝑢(0)(1)𝐷𝑢+𝑢=𝑓,=0,𝑢=0, (55) where 𝑓(𝑥)=𝑥²−𝑥³+12𝑓₁(𝑥)+12𝑓₂(𝑥),𝑓₁(𝑥)=−2𝑥^𝛽Γ+1+𝛽6𝑥^1+𝛽Γ,𝑓2+𝛽₂(𝑥)𝛽()=−1−𝑥^𝛽−1Γ+4()1+𝛽𝛽+11−𝑥^𝛽Γ−6()2+𝛽𝛽+21−𝑥^𝛽+1Γ.3+𝛽 (56)

As 𝑢∈𝐻²₀(Ω), Corollary 29 predicts a rate of convergence of 2 in 𝐿² 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 ‖𝑢−𝑢_𝑘‖_𝐿²_(Ω) for different 𝛽.

As 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, 𝐴_𝑘𝑢^𝑘_𝑛𝑘−𝑔_𝐿2(Ω)≤10^−𝑛, (57) to solve the corresponding system. For escaping “out of memory,” we define that the data type of 𝐴_𝑘 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 𝑛=4. The rate of the increasing CPU time is defined by𝑅_cpu=log[CPUtimeonfinermesh/CPUtimeoncoarsermesh][].logdimoffinermesh/dimofcoarsermesh (58)From the table, we see that the numbers of iterations by our FV-MGM are independent of ℎ_𝑘. In contrast, the CGNR method (with the initial value 𝑧₀=0) needs more iterations when ℎ_𝑘 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 Example 30 with 𝛽=0.7 by the GE method, the CGNR method, and the FV-MGM, respectively.

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 𝒪(dim𝑉_𝑘) while the stopping criterion is taken (57). Different from integer-order equations, all the convergence analysis is on fractional Sobolev spaces 𝐻^𝛼(Ω). The nonsymmetric form of FADE will be studied in the future.