# An Efficient Algorithm with Stabilized Finite Element Method for the Stokes Eigenvalue Problem

Authors: Zhifeng Weng; Yaoxiong Cai
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper provides a two-space stabilized mixed finite element scheme for the Stokes eigenvalue problem based on local Gauss integration. The two-space strategy contains solving one Stokes eigenvalue problem using theP1-P1 finite element pair and then solving an additional Stokes problem using the P2-P2 finite element pair. The postprocessing technique which increases the order of mixed finite element space by using the same mesh can accelerate the convergence rate of the eigenpair approximations. Moreover, our method can save a large amount of computational time and the corresponding convergence analysis is given. Finally, numerical results are presented to confirm the theoretical analysis.


---

## Body

## 1. Introduction

The Stokes eigenvalue problem is one of the most important eigenvalue problems and plays an important role in the stability analysis of nonlinear partial differential equations [1]. The eigenvalue problems are used in many application areas: structural mechanics and fluid mechanics. Thus, development of the efficient numerical methods for studying the eigenvalue problems has practical meanings and has been noticed by many researchers. At the time of writing, numerous works are devoted to these problems (see [2–10] and the references cited therein).Many effective postprocessing strategies that improve the convergence rate for the approximations of the eigenvalue problems by the finite element methods have been well developed. The two-grid method is one of these efficient postprocessing methods. The basic idea of two-grid scheme is first introduced by Xu [11, 12] for the nonsymmetric and nonlinear elliptic problems. Hence, it can be seen as a postprocessing technique and can take less CPU time compared to the one grid methods. To the best of our knowledge, some details of the two-grid scheme can be found in the works of Xu and Zhou [13], Chien and Jeng [14, 15], Chen et al. [7, 16], Hu and Cheng [17], Yang et al. [18, 19], Huang et al. [8], and Weng et al. [20, 21]. The two-space method is actually the iterative Galerkin method, which was first used for solving integral equation eigenvalue problems by Sloan [22] and differential equation eigenvalue problems by Lin and Xie [23]. Particularly, Racheva and Andreev [24] have proposed a postprocessing method for the 2m-order self-adjoint eigenvalue problems by two-grid method or the two-space method. A similar method has been given for the Stokes eigenvalue problem [7, 25], elliptic eigenvalue problem [16], and the biharmonic eigenvalue problem [26] by mixed finite element methods.In fact, two-space method [27–29] can be cast in the framework of Xu’s work regarding the two-grid method. However, the two-space method is different from the two-grid method. This two-space method consists in solving the original Stokes eigenvalue problem in the k-order mixed finite element space and one additional Stokes source problem in an augmented mixed finite element space by a k+1-order mixed finite element space on the same mesh. Besides, the two-space method only needs one mesh size while the two-grid method needs two mesh sizes, a coarse mesh, and a fine mesh. In fact, the two-space method can avoid the discussion on the relation of the coarse and fine meshes. For this reason, in the present paper we establish a two-space discretization scheme for the Stokes eigenvalue problem.Recently, more attention has been paid to the lowest equal order finite element pairs for simulating the incompressible flow. The lowest equal order finite element pairs offer some computational advances; for example, they are simple and have practical uniform data structure and adequate accuracy, because they show an identical degree distribution for both the velocity and pressure. Moreover, they are of practical importance in scientific computation owing to their very convenient computational cost. However, the lowest equal order mixed finite element pairs do not satisfy the inf-sup condition. Numerical tests show that the violation of the inf-sup condition often brings about unphysical pressure oscillations. In order to avoid the instability problem, the stabilized finite element methods are applied to the incompressible flow. Therefore, a lot of work focuses on stabilization (see [30–37]) of the lowest equal order pairs. Particularly, based on the work of Bochev et al. [30], Li et al. [31, 32] used the projection of the pressure onto the piecewise constant space to add the stabilized term for P1-P1 element and Zheng et al. [35] used the projection of the pressure-gradient onto the piecewise constant space to add the stabilized term for P2-P2 element.Influenced by the work mentioned above, the paper focuses on the method, which combines two-space discretization scheme with a stabilized finite element method based on local Gauss integration technique for the Stokes eigenvalue problem. The paper is organized as follows. In Section2, we introduce the studied Stokes eigenvalue problem and the notations and some well-known results used throughout this paper. Some stabilized finite element strategies based on two local Gauss integrations are recalled in Section 3. In Section 4, a two-space stabilized finite element algorithm is constructed and its error estimates are discussed. In Section 5, numerical experiments are reported for illustrating the theoretical results and the high efficiency of the proposed method. Finally, we will conclude our presentation in Section 6 with a few comments and also possible future research topics.


## 2. Preliminaries

In this paper, we consider the following Stokes eigenvalue problem:(1)-Δu+∇p=λuinΩ,div⁡u=0inΩ,u=0onΓ,where Ω⊂R2 is a bounded and convex domain with a Lipschitz-continuous boundary Γ, p(x) represents the pressure, u(x) is the velocity vector, and λ∈R is the eigenvalue.We shall introduce the following Hilbert spaces:(2)V=H01Ω2,Y=L2Ω2,W=L02Ω=q∈L2Ω,∫Ωqdx=0. The spaces [L2(Ω)]m, m=1,2, are equipped with the L2-scalar product (·,·) and L2-norm ·0. The norm and seminorm in [Hk(Ω)]2 are denoted by ·k and ·k, respectively. The space V is equipped with the norm ∇·0 or its equivalent norm ·1 due to Poincaré inequality. Spaces consisting of vector-valued functions are denoted in boldface. Furthermore, the norm in the space dual to V is given by(3)u-1=supv∈V,v1=1u,v.Therefore, we define the following bilinear forms a(·,·), d(·,·), and b(·,·) on V×V, V×W, and V×V, respectively, by (4)au,v=∇u,∇v,∀u,v∈V,dv,q=div⁡v,q,∀v∈V,∀q∈W,bu,v=u,v,∀u,v∈V and a generalized bilinear form B((·,·),(·,·)) on (V×W)×(V×W); that is,(5)Bu,p,v,q=au,v-dv,p+du,q,∀u,p,v,q∈V×W.With the above notations, the variational formulation of problem (1) reads as follows: Find (u,p;λ)∈(V×W)×R with u0=1, such that(6)Bu,p,v,q=λbu,v,∀v,q∈V×W.From [1], we know that eigenvalue problem (6) has an eigenvalue sequence {λj}(7)0<λ1≤λ2≤λ3≤⋯ and corresponding eigenvectors (8)u1,u2,u3,…, with the orthogonal property b(ui,uj)=δij.Let(9)Mλi=u∈V,uis  an  eigenvector  of6corresponding  toλi.Moreover, the bilinear formd(·,·) satisfies the inf-sup condition for all q∈W(10)supv∈Vdv,qv1≥βq0,where β>0 is a constant depending only on Ω. Therefore, the generalized bilinear form B satisfies the continuity property and coercive condition(11)Bu,p,v,q≤Cu1+p0×v1+q0,supv,q∈V,WBu,p,v,qv1+q0≥β1u1+p0,where C and β1 are the positive constants depending only on Ω. Throughout the paper we use c or C to denote a generic positive constant whose value may change from place to place, which remains independent of the mesh parameter.


## 3. A Stabilized Mixed Finite Element Method

From now on,h is a real positive parameter tending to 0. The finite element subspaces V¯h×W¯h,Vh×Wh of V×W are characterized by Th, a partitioning of Ω into triangles T with the mesh size h, assumed to be uniformly regular in the usual sense [38]. Then we define them as follows:(12)V¯h=vh=v1,v2∈C0Ω2∩V:viT∈P1T,∀T∈Th,i=1,2,W¯h=w∈C0∩W:wT∈P1T,∀T∈Th,Vh=vh=v1,v2∈C0Ω2∩V:viT∈P2T,∀T∈Th,i=1,2,Wh=w∈C0∩W:wT∈P2T,∀T∈Th,Mh1=vh=v1,v2∈C0Ω2∩V∣viT∈P1T⊕BT,∀T∈Th,i=1,2,Mh2=vh=v1,v2∈C0Ω2∩V∣viT∈P2T⊕BT,∀T∈Th,i=1,2, where Pk(T) represents the set of all polynomials on T of degree less than k∈N and B(T) denotes the space of bubble functions. The bubble functions are defined as follows: (13)BT=vh∈CT∣vh∈Span⁡λ0λ1λ2∀T∈Th, where λi are area coordinates on T, i=0,1,2. The area coordinate is also known as a triangle barycentre coordinate, where the three components (λ0,λ1,λ2) are of the ratio between the area of the three triangles and the area of the mother triangle.It is known that this choice of the approximate spacesMh1×W¯h or Mh2×Wh satisfies the inf-sup condition in [38], but this choice of the approximate spaces V¯h×W¯h or Vh×Wh does not satisfy the inf-sup condition [30, 32, 35]. As a consequence, we give a stabilized finite element approximation based on local Gauss integration technique (see [32, 35]). The idea is as follows.LetΠ:L2(Ω)→R0 be the standard L2-projection:(14)p,qΠp,q,∀p∈W,q∈R0, where R0={q∈W:qT∈P0(T),∀T∈Th}.The projection operatorΠ has the following properties:(15)Πp0≤cp0,∀p∈W,p-Πp0≤chp1,∀p∈H1Ω. The P1-P1 stabilized bilinear terms are used by(16)G¯p¯h,q=p¯h-Πp¯h,q-Πq,p¯h,q∈W¯h,and the P2-P2 stabilization term is given by(17)Gp,q=∇p-Π∇p,∇q-Π∇q,∀p,q∈Wh.The stabilized term which is defined by local Gaussian quadrature can be rewritten as(18)Gp,q=∑T∈Th∫T,2∇p·∇qdx-∫T,1∇p·∇qdx,∀p,q∈Wh, where ∫T,i‍g(x,y)dx denotes a Gaussian quadrature over T which is exact for polynomials of degree i,i=1,2. In particular, when i=1, the trial function ∇p∈Wh is projected to the piecewise constant space. Besides, the stabilized term G¯(p¯h,q) can be rewritten as (19)G¯p¯h,q=∑T∈Th∫T,2p¯h·qdxdy-∫T,1p¯h·qdxdy,∀p¯h,q∈W¯h, where the trial function p¯h∈W¯h must be projected to R0 when i=1 for any q∈W¯h. Indeed, Becker and Hansbo have found [33] that the stabilized methods of [30, 32] are identical from a numerical point of view for these low-order approximations.By adding the stabilization term into the generalized bilinear formB((·,·),(·,·)), we define (20)Bhuh,ph,v,q=Buh,ph,v,q-Gph,q,B¯hu¯h,p¯h;v,q=Bu¯h,p¯h;v,q-G¯p¯h,q.Then the corresponding discrete variational formulation for the Stokes eigenvalue problem reads as follows: find(uh,ph;λh)∈(Vh×Wh)×R with uh0=1, such that(21)Bhuh,ph,v,q=λhruh,v,∀v,q∈Vh×Whand find (u¯h,p¯h;λ¯h)∈(V¯h×W¯h)×R with u-h0=1, such that(22)B¯hu¯h,p¯h;v,q=λ¯hru¯h,v,∀v,q∈V¯h×W¯h.Remark 1.
For thePkb-Pk (Mhk×Wh,k=1,2) pair which satisfy inf-sup condition, there are points of difference between them. The Pk-Pk stabilized method in this article only adds the stabilized term with respect to the pressure space. However, the Pkb-Pk method has the implicit stabilized term in the velocity space.Withv=uh,q=ph and thanks to the positive definiteness of a(uh,uh), we deduce that the discrete eigenvalues λjh are positive. Let the eigenvalue of (21) be ordered as follows: (23)0<λ1h≤λ2h≤λ3h≤⋯≤λNh, and let us consider the corresponding eigenfunctions (24)u1h,p1h,u2h,p2h,u3h,p3h,…,uNh,pNh, where r(uih,ujh)=δij,1≤i,j≤N,N denotes the dimension of the finite element space.Similarly, letMh(λih) be the eigenspace associated with λih; that is, (25)Mhλih=uh∈Vh,uhis  an  eigenfunction  of21corresponding  toλih.For (22) with P1-P1 pairs, it can be given similarly. The corresponding nature of the eigenvalues is omitted for the sake of simplicity.The next theorem shows the continuity property and the weak coercivity property of the bilinear formBh((uh,ph),(v,q)) for the finite element space Vh×Wh in [35] and B¯h((u¯h,p¯h);(v,q)) for the finite element space V¯h×W¯h in [30, 32].Theorem 2.
For all(uh,ph),(v,q)∈Vh×Wh, there exist positive constants C and β2, independent of h, such that(26)Bhuh,ph,v,q≤Cuh1+ph0×v1+q0,supv,q∈Vh,WhBhuh,ph,v,qv1+q0≥β2uh1+ph0. Moreover, for all (u¯h,p¯h),(v,q)∈(V¯h×W¯h), there exist positive constants C1 and β3, independent of h, such that(27)B¯hu¯h,p¯h;v,q≤C1u¯h1+p¯h0×v1+q0,supv,q∈V¯h,W¯hB¯hu¯h,p¯h;v,qv1+q0≥β3u¯h1+p¯h0.The next theorem contains the convergence result of eigenfunctions and eigenvalues for the Stokes eigenvalue problem in [8, 20].Theorem 3.
With(u,p,λ) belonging to (H3(Ω)2∩V)×(H2(Ω)∩W)×R and being the exact solution of (6), one deduces that there exists a discrete eigenpair (uh,ph;λh) of (21) which satisfies the following error estimates:(28)u-uh0+hu-uh1+p-ph0≤ch3,λ-λh≤ch4.Furthermore, if the exact solution(u,p;λ)∈(H2(Ω)2∩V)×(H1(Ω)∩W)×R, then (u¯h,p¯h;λ¯h)∈V¯h×W¯h×R of problem (22) satisfies(29)u-u¯h0+hu-u¯h1+p-p¯h0≤ch2,λ-λ¯h≤ch2.


## 4. Two-Space Stabilized Finite Element Scheme and Error Estimates

In this section, we shall present a two-space stabilized finite element algorithm to reduce the computational cost. The two-space stabilized finite element approximation consists of three steps.Step 1.
On the mesh sizeh, solve the following Stokes eigenvalue problem by P1-P1 pair and find (u¯h,p¯h;λ¯h)∈(V¯h×W¯h)×R and u¯h0=1, such that, for all (v,q)∈V¯h×W¯h,(30)B¯hu¯h,p¯h;v,q=λ¯hu¯h,v.Step 2.
On the same mesh sizeh, solve the following Stokes problem by P2-P2 pair and find (uh,ph)∈Vh×Wh such that for all (vh,qh)∈Vh×Wh(31)Bhuh,ph;vh,qh=λ¯hu¯h,vh.Step 3.
Compute the eigenvalue by the Rayleigh quotient(32)λh=Bhuh,ph;uh,phuh,uh,where uh∈Vh∖{0}.Next, we will study the convergence of the two-space stabilized finite element solution. To do this, we define the Galerkin projection operator(Rh,Qh):(V,W)→(Vh,Wh) by(33)BhRhv,q,Qhv,q,vh,qh=Bv,q,vh,qh,∀vh,qh∈Vh,Wh.By Theorem 2, (Rh,Qh) is well defined and the following approximation properties are fulfilled in [20].Lemma 4.
For all(u,p)∈(H3(Ω)2∩V,H2(Ω)∩W), one has(34)u-Rhu,p1+p-Qhu,p0≤ch2u3+p2.The following identity that relates the errors in the eigenvalue and eigenvector can be found in [3].Lemma 5.
Let(u,p;λ) be an eigenpair of (6); for any s∈V∖{0} and w∈W, one has(35)Bs,w,s,wrs,s-λ=Bs-u,w-p,s-u,w-prs,s-λrs-u,s-urs,s.The next theorem provides the error estimates for our two-space scheme.Theorem 6.
Let(uh,ph;λh) be the ith discrete eigenpair. Then the ith eigenpair (u,p;λ) of the Stokes operator is such that(36)u-uh1+p-ph0≤ch2u3+p2,(37)λ-λh≤ch4u3+p22.Proof.
Denoted by(θh,ρh)=(Rh(u,p)-uh,Qh(u,p)-ph), subtracting (6) from (31), we derive from (33)(38)Bhθh,ρh;v,q=λu-u¯h,v+λ-λ¯hu¯h,v.Let (v,q)=(θh,ρh) in (38), by using Theorem 2, Sobolev embedding theorem, and Theorems 3, we obtain (39)ν∇θh0+ρh0≤β2-1λu-u¯h-1+λ-λ¯hu¯h-1≤ch2u2+p1.
Combining the triangle inequality with Lemma4, we deduce(40)u-uh1+p-ph0≤θh1+u-Rhu,p1+ρh0+p-Qhu,p0≤ch2u3+p2,and finally we obtain (36).
Moreover, using (32) and Lemma 5, we have (41)λh-λ=Bhuh,ph;uh,phuh,uh-λ=Buh-u,ph-p;uh-u,ph-p-Gph,phuh,uh-λuh-u,uh-uuh,uh.
Taking the norm and using (15) and (36), we conclude the proof that is (42)λ-λh≤cu-uh12+h2ph-Πph02+cp-ph02≤cu-uh12+h2ph-p+p-Πp+Πp-Πph02+cp-ph02≤cu-uh12+h2p-Πp02+cp-ph02≤ch4u3+p22.Remark 7.
From Theorem3, for the usual P2-P2 stabilized finite element solution (uh,ph;λh) which involves solving a Stokes eigenvalue problem with mesh size h, we have the following error estimates:(43)u-uh0+hu-uh1+p-ph0≤ch3,λ-λh≤ch4.Furthermore, if we use the two-space stabilized finite element method, then we get the convergence rate of the same order as the usual stabilized finite element method from Theorem6. However, our method is more efficient than the P2-P2 stabilized finite element scheme in the same mesh because our method for solving Stokes eigenvalue problem is to compute an initial approximation based on a lower number of nodes, which takes less CPU time.Remark 8.
For the two-space algorithm withP2b-P2(Mh2×Wh) pair which satisfies inf-sup condition, we obtain the same result. The procedure of the two-space method with P2b-P2 pair could be described in the following manner: Firstly, we solve the following Stokes eigenvalue problem by P1b-P1 (Mh1×Wh) pair on the mesh size h; then, we should solve the new Stokes problem by P2b-P2 pair on the same mesh size h. Finally we can compute the eigenvalue by the Rayleigh quotient.


## 5. Numerical Experiments

In this section we present numerical results to check the theoretical analysis contained in Theorem6. Our goal is to confirm the theoretical results of the new two-space stabilized finite element method for the two-dimensional Stokes eigenvalue approximated by the equal order finite element pairs based on local Gauss integration.In our numerical experiments,Ω is the unit square domain [0,1]×[0,1] in R2. The domain Ω is uniformly divided by the triangulations of mesh size h. Here, we just consider the first eigenvalue of the Stokes eigenvalue problem for the sake of simplicity. Following [4], we employ the approximation λ1=52.3446911 as the reference solution for the first eigenvalue. Note that in these computations we set ν=1.When solving the Stokes problem with a mesh sizeh, we need the solutions λ¯h and u¯h generated by a lower finite element pair P1-P1. To do this we interpolate the solutions λ¯h and u¯h onto the grid with the same mesh size h, but increasing the order of the mixed finite element space. In conclusion, the solution of the two-space method is obtained by one simple eigenvalue problem by a lower finite element pair and one time interpolation by a higher finite element pair P2-P2.Our goal in this test is to validate the merit of the two-space method as compared with theP2-P2 stabilized method and the two-space method with P2b-P2 pair. The eigenvalue approximation λh, the eigenvalue error, the convergence rates, and the CPU time for the stabilized mixed finite element methods for different values of h are tabulated in Tables 1, 2, and 3. From Tables 1, 2, and 3, we can see that the three methods work well and keep the convergence rates just as predicted by the theoretical analysis, but our two-space method can take less CPU time. For the two-space method with P2b-P2 pair, the two-space method with P2-P2 pair approximates the velocity variable with a lower number of nodes, so our method can save a lot of time.Table 1
Relative error and convergence rate forP2-P2 pair.
1
/
h
λh
λ-λh/λ
Rate
CPU time
8
52.4269
1.570
E
-
3
0.156
16
52.3505
1.111
E
-
4
3.821
0.672
32
52.3451
7.345
E
-
6
3.919
3.469
64
52.3447
5.245
E
-
7
3.808
27.001Table 2
Relative error and convergence rate of two-space method withP2-P2 pair.
1
/
h
λh
λ-λh/λ
Rate
CPU time
8
52.4594
2.191
E
-
3
0.14
16
52.3529
1.570
E
-
4
3.803
0.609
32
52.3452
1.029
E
-
5
3.930
2.766
64
52.3447
6.536
E
-
7
3.977
13.875Table 3
Relative error and convergence rate of two-space method withP2b-P2 pair.
1
/
h
λh
λ-λh/λ
Rate
CPU time
8
52.518
3.311
E
-
3
0.218
16
52.3594
2.817
E
-
4
3.555
0.875
32
52.3458
2.033
E
-
5
3.793
3.984
64
52.3448
1.338
E
-
6
3.926
22.916Next numerical test is about the second, third, and fourth eigenvaluesλ2,3,4h. The reference values are computed over a fine mesh h=1/64 and the results are λ2,3,4=92.1245411,92.1245843,128.209971. Then, in Figure 1, we exhibit the O(h4) convergence rate as has been predicted in Theorem 6 with the two-space method.Figure 1
The convergence rate of the eigenvalue forλ2,3,4 on the unit square with the two-space method.Moreover, we give two plots of numerical solutions of two kinds of two-space schemes at the mesh1/h=48 in Figure 2 for the details. Figure 2 shows the stability of two schemes.Figure 2
Plots of the velocity and pressure ath=1/48: numerical solution of two-space method with P2-P2(a–c) and numerical solution of two-space method with P2b-P2(d–f) for u1h,u2h,ph.
(a)
(b)
(c)
(d)
(e)
(f)


## 6. Conclusions

In this paper, we presented the two-space algorithm for the Stokes eigenvalue problem discretized by stabilized mixed finite element scheme, based on local Gauss integration technique. The main feature of our method is to combine two equal order stabilized methods, then use the first-order mixed finite element space to solve the original Stokes eigenvalue problem, and solve the Stokes source problem in the second-order mixed finite element space on the same mesh. Moreover, the related error estimates have been derived. Finally, numerical tests show that the two-space stabilized mixed finite element method is numerically efficient for solving the Stokes eigenvalue problem. The two-space algorithm can achieve the same accuracy as the stabilized finite element solution as theP2-P2 stabilized method by taking less CPU time. Obviously, this method can be extended to the case of three dimensions easily. And there are some open questions including the possible extension of the method to other linear and nonlinear eigenvalue problems.



---

*Source: 6362505-2017-12-31.xml*