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

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


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).
FADE has been used in modeling physical phenomena exhibiting anomalous diffusion, that is, diffusion not accurately modeled by the usual advection dispersion equation [4][5][6][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][9][10].Many scholars developed numerical methods, including finite difference method [11], finite element method [12][13][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 O( 3 ) or O( log 2 ).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.

Definition 2 (variational solution). A function 𝑢 ∈ 𝐻
Ervin and Roop ( [17], 2006) proved that the bilinear form (⋅, ⋅) satisfies the coercivity and continuity; that is, there exist two positive constants  1 and  2 such that 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, 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  ℎ ⊂   0 (Ω) as Let  ℎ be the solution to the finite-dimensional variational problem: Note that the existence and uniqueness of solution to (9) follow from the fact that  ℎ is a subset of the space   0 (Ω) ∩   0 (Ω).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  2 norm was proved.

Theorem 5. Under the same assumptions in Theorem 4, one has
where

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 {T  } denote a sequence of partitions of Ω and ℎ  be the mesh size of T  .From now on, we assume that {T  } is a quasiuniform family, that is, with positive constant .At the same time, let T  be obtained from T −1 via a regular subdivision, that is, Let   denote  0 piecewise linear functions with respect to T  that vanish on Ω.

Mesh-Dependent Norms
Definition 6.The mesh-dependent inner product (⋅, ⋅)  on   is defined by Definition 7. The operator   :   →   is defined by In terms of the operator   , the discretized equation of ( 9) can be written as where   ∈   satisfies Remark 8. (i) From the Riesz representation theorem, we point out that the operator   is defined uniquely.(ii) Since (V, ) is symmetric and satisfies the proposition of coercivity,   is symmetric positive definite woth respect to (⋅, ⋅)  .
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 [ Lemma 11 (spectral radius theorem).We have the estimation for spectral radius Λ(  ): where  7 is a positive constant independent of . Proof.

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, The fine-to-coarse operator  −1  :   →  −1 is defined by 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 (,  0 , ) the approximate solution of the    =  obtained from the th level iteration with initial guess  0 .The discussion of parameters Λ  and   is left to the next subsection and Section 4.
Algorithm 13 (the th level V-cycle multigrid algorithm).

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, Definition 16.Define the relaxation operator: Throughout this paper, we assume that Λ  denotes some upper bound for the spectral radius of   satisfying Λ  ≤  7 ℎ −2  .
Remark 18.Let  =   V. From the Riesz representation theorem there exists a unique  ∈  − (Ω) such that which means that V is a finite element solution of variational problem (4) on   .Observing that we claim that  =  −1 V is also a finite element solution of variational problem (4) on  −1 .Therefore, noting  −1 ⊂   (since T −1 ⊂ T  ) and applying Theorem 4, we have (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 where  9 =  4 /.
Lemma 21.Let  be the number of smoothing steps.Then Proof.See Lemma 6.6.7 in [18].

Convergence Analysis.
In this subsection, we prove convergence theorems for V-cycle FE multigrid method.
Corollary 26.Let  be the number of smoothing steps.Then one has the estimation for spectral radius Λ(  ): to solve the corresponding system.For escaping "out of memory, " we define that the data type of   is short float in our program.
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 = 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.

Concluding Remarks
In general the discretized system of FADE has a full and ill-conditioned coefficient matrix, so FV-MGM is a highefficient 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 O(dim   ) while the stopping criterion is taken (57).Different from integerorder equations, all the convergence analysis is on fractional Sobolev spaces   (Ω).The nonsymmetric form of FADE will be studied in the future.

Table 2 :
Comparisons for solving Example 30 with  = 0.7 by the GE method, the CGNR method, and the FV-MGM, respectively.uniformpartition of [0, 1], which support the predicted rates of convergence for different values of .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 − , a 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 [CPU time on finer mesh/CPU time on coarser mesh] log [dim of finer mesh/ dim of coarser mesh] .