Numerical Investigations on Several Stabilized Finite Element Methods for the Stokes Eigenvalue Problem

Abstract

Several stabilized finite element methods for the Stokes eigenvalue problem based on the lowest equal-order finite element pair are numerically investigated. They are penalty, regular, multiscale enrichment, and local Gauss integration method. Comparisons between them are carried out, which show that the local Gauss integration method has good stability, efficiency, and accuracy properties, and it is a favorite method among these methods for the Stokes eigenvalue problem.

1. Introduction

It is well known that numerical approximation of eigenvalue problems plays an important role in the analysis of the stability of nonlinear PDEs. Meanwhile, they are wildly used in many application areas: structural mechanics and fluid mechanics (see [1]). Thus, development of an efficient and effective computational method for investigating these problems has practical significance and has drawn the attention of many peoples. At the time of writing, numerous works are devoted to these problems (see [2–9], and the references cited therein). Yin et al. [10] derived a general procedure to produce an asymptotic expansion for eigenvalues of the Stokes problem on rectangular mesh. Chen and Lin [11] proposed the stream function-vorticity-pressure method for the eigenvalue problem. Rate of convergence estimates were derived for the approximation of eigenvalues and eigenvectors by Mercier et al. [12]. Xu and Zhou [13] proposed a two-grid discretization scheme for solving eigenvalue problems.

Mixed finite element methods are a natural choice for solving fluid mechanics equations because these equations naturally appear in mixed form in terms of velocity and pressure [14, 15]. In the analysis and practice of employing mixed finite element methods in solving the Stokes equations, the inf-sup condition has played an important role because it ensures a stability and accuracy of the underlying numerical schemes. A pair of finite element spaces that are used to approximate the velocity and the pressure unknowns are said to be stable if they satisfy the inf-sup condition. Intuitively speaking, the inf-sup condition is something that enforces a certain correlation between two finite element spaces so that they both have the required properties when employed for the Stokes equations. However, due to computational convenience and efficiency in practice, some mixed finite element pairs which do not satisfy the inf-sup condition are also popular. Thus, much attention has been paid to the study of the stabilized method for the Stokes problem.

Recent studies have focused on stabilization techniques, which include penalty [16–18], regular [19], multiscale enrichment [20], and local Gauss integration method [21, 22]. There exist a lot of theoretical results for the stabilized mixed finite element methods for the Stokes equations, and the comparisons between them are also given (see [17, 19, 20, 22–24], and the references cited therein). In this paper, we mainly focus on the Stokes eigenvalue problem solved by these stabilized finite element methods based on the lowest equal-order finite element space pair. Moreover, we present the comparisons between these methods for the considered problem.

The remainder of this paper is organized as follows. In the Section 2, we introduce the studied Stokes eigenvalue problem, the notations, and some well-known results used throughout this paper. Then several stabilized mixed finite element methods are reviewed, and their key stabilization techniques are recalled in Section 3. Comparisons between these stabilized methods are performed numerically in Section 4. Finally, we end with some short conclusions in Section 5.

2. Preliminaries

Let Ω be a bounded, convex, and open subset of ℝ² with a Lipschitz continuous boundary 𝜕Ω. We consider the Stokes eigenvalue problem as follows. Find (𝑢,𝑝;𝜆), such that−𝜈Δ𝑢+∇𝑝=𝜆𝑢inΩ,div𝑢=0inΩ,𝑢=0on𝜕Ω, (2.1) where 𝑢=(𝑢₁(𝑥),𝑢₂(𝑥)) represents the velocity vector, 𝑝=𝑝(𝑥) the pressure, 𝜈>0 the viscosity, and 𝜆∈ℝ the eigenvalue.

We will introduce the following Hilbert spaces: 𝑋=𝐻¹₀(Ω)²,𝑌=𝐿²(Ω)²,𝑀=𝐿²₀(Ω)=𝑞∈𝐿²(Ω)∶_Ω𝑞d𝑥=0. (2.2) The spaces 𝐿²(Ω)^𝑚, 𝑚=1,2, are equipped with the 𝐿²-scalar product (⋅,⋅) and 𝐿²-norm ‖⋅‖_𝐿² or ‖⋅‖₀. The space 𝑋 is endowed with the usual scalar product (∇𝑢,∇𝑣) and the norm ‖∇𝑢‖₀. Standard definitions are used for the Sobolev spaces 𝑊^𝑚,𝑝(Ω), with the norm ‖⋅‖_𝑚,𝑝, 𝑚,𝑝≥0. We will write 𝐻^𝑚(Ω) for 𝑊^𝑚,2(Ω) and ‖⋅‖_𝑚 for ‖⋅‖_𝑚,2.

We define the continuous bilinear forms 𝑎(⋅,⋅) and 𝑑(⋅,⋅) on 𝑋×𝑋 and 𝑋×𝑀, respectively, by 𝑎()()𝑑=𝑢,𝑣=𝜈∇𝑢,∇𝑣,∀𝑢,𝑣∈𝑋,𝑣,𝑞𝑞,div𝑣,∀𝑣∈𝑋,∀𝑞∈𝑀, (2.3) and a generalized bilinear form 𝐵((⋅,⋅);(⋅,⋅)) on (𝑋×𝑀)×(𝑋×𝑀) by 𝐵𝑢,𝑝;𝑣,𝑞=𝑎(𝑢,𝑣)−𝑑𝑣,𝑝−𝑑𝑢,𝑞,∀𝑢,𝑝,𝑣,𝑞∈𝑋×𝑀. (2.4)

With the above notations, the variational formulation of problem (2.1) reads as follows. Find (𝑢,𝑝;𝜆)∈(𝑋×𝑀)×ℝ with ‖𝑢‖₀=1, such that for all (𝑣,𝑞)∈𝑋×𝑀,𝐵𝑢,𝑝;𝑣,𝑞=𝜆(𝑢,𝑣). (2.5)

Moreover, the bilinear form 𝑑(⋅,⋅) satisfies the inf-sup condition [25] for all 𝑞∈𝑀:sup_𝑣∈𝑋𝑑𝑣,𝑞‖∇𝑣‖₀≥𝛽₁𝑞₀, (2.6) where 𝛽₁ is a positive constant depending only on Ω. Therefore, the generalized bilinear form 𝐵((⋅,⋅);(⋅,⋅)) satisfies the continuity property and inf-sup condition: 𝐵𝑢,𝑝;𝑣,𝑞≤𝑐𝜈‖∇𝑢‖₀+𝑝₀‖∇𝑣‖₀+𝑞₀,∀𝑢,𝑝,𝑣,𝑞∈𝑋×𝑀,sup_{𝑣,𝑞∈𝑋×𝑀}𝐵𝑢,𝑝;𝑣,𝑞‖∇𝑣‖₀+𝑞₀≥𝛽₂𝜈‖∇𝑢‖₀+𝑝₀,∀𝑢,𝑝∈𝑋×𝑀, (2.7) where 𝑐 and 𝛽₂ are the positive constants depending only on Ω.

3. Stabilized Mixed Finite Element Methods

For ℎ>0, we introduce the finite-dimensional subspaces 𝑋_ℎ×𝑀_ℎ⊂𝑋×𝑀, which are characterized by 𝐾_ℎ, a partitioning of Ω into triangles 𝐾 with the mesh size ℎ, assumed to be uniformly regular in the usual sense.

Then we define 𝑋_ℎ=𝑢∈𝐶⁰Ω²∩𝑋∶𝑢|_𝐾∈𝑃₁(𝐾)²,∀𝐾∈𝐾_ℎ,𝑀_ℎ=𝑞∈𝐶⁰Ω∩𝑀∶𝑞|_𝐾∈𝑃₁(𝐾),∀𝐾∈𝐾_ℎ, (3.1) where 𝑃₁(𝐾) represents the space of linear functions on 𝐾.

Remark 3.1 (Nonconforming finite element space). Denote the boundary edge by Γ_𝑗=𝜕Ω∩𝜕𝐾_𝑗 and the interior boundary by Γ_𝑗𝑘=Γ_𝑘𝑗=𝜕𝐾_𝑗∩𝜕𝐾_𝑘. Set the centers of Γ_𝑗 and Γ_𝑗𝑘 by 𝜁_𝑗 and 𝜁_𝑗𝑘, respectively. The nonconforming finite element space NC_ℎ for the velocity will be taken to be NC_ℎ=𝑣∶𝑣_𝑗=𝑣∣_𝐾𝑗∈𝑃₁𝐾_𝑗²,𝑣_𝑗𝜁_𝑗𝑘=𝑣_𝑘𝜁_𝑘𝑗,𝑣𝜁_𝑗=0,𝐾_𝑗∈𝐾_ℎ,∀𝑗,𝑘, (3.2) where 𝑃₁(𝐾_𝑗) is the set of all polynomials on 𝐾_𝑗 of degree less than 1. Note that the nonconforming finite element space NC_ℎ is not a subspace of 𝑋. In this nonconforming case, the pair of finite element spaces is NC_ℎ×𝑀_ℎ [26]; that is, the conforming space is still used for pressure.

Moreover, the discrete mixed finite element formulation for the Stokes eigenvalue problem reads: find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(𝑋_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈𝑋_ℎ×𝑀_ℎ,𝐵𝑢_ℎ,𝑝_ℎ;𝑣,𝑞=𝜆_ℎ𝑢_ℎ,𝑣. (3.3)

Next, we denote by 𝑈 the array of the velocity and by 𝑃 the array of the pressure. It is easy to see that (3.3) can be written in matrix form:𝐵𝐴𝐵_𝑇0𝑈𝑃=𝜆_ℎ𝑈𝑃𝐸000, (3.4) where the matrices 𝐴, 𝐵, and 𝐸 are deduced in the usual manner, using the bases for 𝑋_ℎ and 𝑀_ℎ, from the bilinear forms 𝑎(⋅,⋅), 𝑑(⋅,⋅), and (⋅,⋅), respectively, and 𝐵^𝑇 is the transpose of matrix 𝐵.

Remark 3.2. In the nonconforming case, we define the bilinear forms 𝑎^′(𝑢,𝑣)=𝜈_{𝐾∈𝐾ℎ}(∇𝑢,∇𝑣)_𝐾,𝑑^′𝑣,𝑝=_{𝐾∈𝐾ℎ}div𝑣,𝑞_𝐾, (3.5) and a generalized bilinear form 𝐵^′((⋅,⋅);(⋅,⋅)) on (𝑋×𝑀)×(𝑋×𝑀) by 𝐵^′𝑢,𝑝;𝑣,𝑞=𝑎^′(𝑢,𝑣)−𝑑^′𝑣,𝑝−𝑑^′𝑢,𝑞,∀𝑢,𝑝,𝑣,𝑞∈𝑋×𝑀. (3.6)

Accordingly, the discrete nonconforming formulation for the Stokes eigenvalue problem is: find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(NC_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈NC_ℎ×𝑀_ℎ, 𝐵^′𝑢_ℎ,𝑝_ℎ;𝑣,𝑞=𝜆_ℎ𝑢_ℎ,𝑣. (3.7)

Note that the lowest equal-order pair does not satisfy the discrete inf-sup condition: sup_{𝑣ℎ∈𝑋ℎ}𝑑𝑣_ℎ,𝑞_ℎ∇𝑣_ℎ0≥𝛽₃𝑞_ℎ0,orsup_{𝑣ℎ∈NCℎ}𝑑^′𝑣_ℎ,𝑞_ℎ∇𝑣_ℎ0,ℎ≥𝛽₃𝑞_ℎ0,∀𝑞_ℎ∈𝑀_ℎ, (3.8) where the constant 𝛽₃>0 is independent of ℎ, and ‖∇𝑣_ℎ‖_0,ℎ is the discrete energy seminorm in the nonconforming case. In order to fulfill this condition, several ways have been used to stabilize the lowest equal-order finite element space pair.

Method 1 (Penalty method). Find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(𝑋_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈𝑋_ℎ×𝑀_ℎ, 𝐵𝑢_ℎ,𝑝_ℎ;+𝜀𝑣,𝑞𝜈𝑝_ℎ,𝑞=𝜆_ℎ𝑢_ℎ,𝑣, (3.9) where 𝜀>0 is a penalty parameter. The performance of this method obviously depends on the choice of the penalty parameter 𝜀. Then the matrix form of (3.9) can be expressed as 𝐵𝐴𝐵_𝑇𝜀𝜈𝐷𝑈𝑃=𝜆_ℎ𝑈𝑃𝐸000, (3.10) where the matrix 𝐷 is deduced in the usual manner, using the base for 𝑀_ℎ, from (𝑝_ℎ,𝑞). Because the coefficient matrix of (3.10) is usually large and sparse, it is not easy to compute the numerical solution using a direct method. In general, one can use the Uzawa algorithm. By some simple calculation, we get 𝐸^𝑃_𝑖+1=𝐼−𝜏𝐵^𝑇𝐴−𝜆_ℎ𝐸⁻¹𝜀𝐵−𝜈𝐷𝐸^𝑃_𝑖, (3.11) where 𝐸^𝑃_𝑖=𝑃−𝑃_𝑖 and 𝜏 is a positive constant.

Method 2 (Regular method). Find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(𝑋_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈𝑋_ℎ×𝑀_ℎ, 𝐵𝑢_ℎ,𝑝_ℎ;𝑣,𝑞−𝛿_{𝐾∈𝐾ℎ}∇𝑝_ℎ−𝜆_ℎ𝑢_ℎ,∇𝑞_𝐾=𝜆_ℎ𝑢_ℎ,𝑣, (3.12) where 𝛿=ℎ²/(𝛼𝜈) is a stabilization parameter and 𝛼>0. The regular method uses a simple way to stabilize the mixed finite element approximation without a loss of accuracy. In fact, it can be treated the regular method as a special Douglas-Wang's scheme [19]. Then the matrix form of this stabilized version can be expressed as 𝐵𝐴𝐵_𝑇−𝛿𝐷₁𝑈𝑃=𝜆_ℎ𝑈𝑃𝐸0𝐹0, (3.13) where additional blocks 𝐷₁ and 𝐹 correspond to the following respective terms: _{𝐾∈𝐾ℎ}∇𝑝_ℎ,∇𝑞_𝐾,−𝛿_{𝐾∈𝐾ℎ}𝑢_ℎ,∇𝑞_𝐾. (3.14) As (3.11),we also have 𝐸^𝑃_𝑖+1=𝐼−𝜏𝐵^𝑇−𝜆_ℎ𝐹𝐴−𝜆_ℎ𝐸⁻¹𝐵+𝛿𝐷₁𝐸^𝑃_𝑖. (3.15)

Method 3 (Multiscale enrichment method). Find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(𝑋_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈𝑋_ℎ×𝑀_ℎ, 𝐵𝑢_ℎ,𝑝_ℎ;𝑣,𝑞−𝛿_{1𝐾∈𝐾ℎ}∇𝑝_ℎ−𝜆_ℎ𝑢_ℎ,∇𝑞_𝐾+𝛿_{2𝑒=𝜕𝐾𝑗∩𝜕𝐾𝑘}𝜈𝜕_𝑛𝑢_ℎ,𝜈𝜕_𝑛𝑣_𝑒=𝜆_ℎ𝑢_ℎ,𝑣, (3.16) where 𝛿₁=ℎ²/(𝛼₁𝜈), 𝛿₂=ℎ/(𝛼₂𝜈)are the stabilization parameters, 𝛼₁,𝛼₂>0, 𝑛 is the normal outward vector, 𝜕_𝑛 is normal derivative operator, and [𝑣] denotes the jump of 𝑣 across 𝑒. This stabilized method includes the usual Galerkin least squares stabilized terms on each finite element and positive jump terms at interelement boundaries. Moreover, a direct algebraic manipulation leads to the matrix form 𝐵𝐴𝐵_𝑇+𝛿₂𝐷₂−𝛿₁𝐷₁𝑈𝑃=𝜆_ℎ𝑈𝑃𝐸0𝐹0, (3.17) where the matrix 𝐷₂ is deduced in the usual manner, using the bases for 𝑋_ℎ, from the term ∑_{𝑒=𝜕𝐾𝑗∩𝜕𝐾𝑘}([𝜈𝜕_𝑛𝑢_ℎ],[𝜈𝜕_𝑛𝑣])_𝑒. As (3.11), we have 𝐸^𝑃_𝑖+1=𝐼−𝜏𝐵^𝑇+𝛿₂𝐷₂−𝜆_ℎ𝐹𝐴−𝜆_ℎ𝐸⁻¹𝐵+𝛿𝐷₁𝐸^𝑃_𝑖. (3.18)

Method 4 (Local Gauss integration method). Find (𝑢_ℎ,𝑝_ℎ;𝜆_ℎ)∈(𝑋_ℎ×𝑀_ℎ)×ℝ with ‖𝑢_ℎ‖₀=1, such that for all (𝑣,𝑞)∈𝑋_ℎ×𝑀_ℎ, 𝐵𝑢_ℎ,𝑝_ℎ;𝑣,𝑞−𝐺𝑝_ℎ,𝑞=𝜆_ℎ𝑢_ℎ,𝑣, (3.19) where 𝐺(𝑝_ℎ,𝑞) is defined by 𝐺𝑝_ℎ,𝑞=_{𝐾∈𝐾ℎ𝐾,2}𝑝_ℎ𝑞d𝑥−_𝐾,1𝑝_ℎ𝑞d𝑥,∀𝑝_ℎ,𝑞∈𝑀_ℎ, (3.20) where ∫_𝐾,𝑖𝑔(𝑥)d𝑥 indicates a local Gauss integral over 𝐾 that is exact for polynomials of degree 𝑖, 𝑖=1,2. In particular, the trial function 𝑝_ℎ∈𝑀_ℎ must be projected to piecewise constant space when 𝑖=1. This stabilization technique is free of stabilization parameters and does not require any calculation of high-order derivatives, a specification of any mesh-dependent parameter or edge-based data structures. Then the corresponding matrix form of this stabilized method is 𝐵𝐴𝐵_𝑇𝑈−𝐺𝑃=𝜆_ℎ𝑈𝑃𝐸000, (3.21) where the matrix 𝐺 is deduced in the usual manner, using the bases for 𝑀_ℎ, from the term 𝐺(𝑝_ℎ,𝑞). As (3.11), we get 𝐸^𝑃_𝑖+1=𝐼−𝜏𝐵^𝑇𝐴−𝜆_ℎ𝐸⁻¹𝐵+𝐺𝐸^𝑃_𝑖. (3.22)

Remark 3.3. It is well known that if 𝜏 is well chosen, then 𝑈_𝑖 and 𝑃_𝑖 converge, respectively, to 𝑈 and 𝑃 with a rate of convergence based on (3.11), (3.15), (3.18), and (3.22). From these equations, we can find that Method 1 converges faster than Method 4. Methods 2 and 3, whose coefficient matrices are not symmetric, may cost more time to converge.

Remark 3.4. By using the regularity assumptions and well-established techniques for eigenvalue approximation [1, 8, 10, 12], the theoretical convergence rates should be of order 𝑂(ℎ²) and 𝑂(ℎ) for the velocity in the 𝐿²- and 𝐻¹-norms, respectively, of order 𝑂(ℎ) for the pressure in the 𝐿²-norm, and of order 𝑂(ℎ²) for the eigenvalue by using all these stabilized methods.

4. Numerical Experiments

In this section we numerically compare the performance of the various stabilized mixed finite element methods discussed in the previous section. In all experiments, the algorithms are implemented using public domain finite element software [27] with some of our additional codes.

Let the computation be carried out in the region Ω={(𝑥,𝑦)∣0<𝑥,𝑦<1}. We consider the Stokes eigenvalue problem in the case of the viscidity 𝜈=1, and it will be numerically solved by the stabilized mixed methods on uniform mesh (see Figure 1). Here, we just consider the first eigenvalue of the Stokes eigenvalue problem for the sake of simplicity. The exact solution of this problem is unknown. Thus, we take the numerical solution by the standard Galerkin method (𝑃₂-𝑃₁ element) computed on a very fine mesh (6742 grid points) as the “exact” solution for the purpose of comparison. Here, we take 𝜆=52.3447 as the first exact eigenvalue.

Figure 1 

Uniform finite element partitioning of the unit square.

As we know, the stabilized term of the regular and multiscale enrichment methods must be controlled by carefully designed stabilization parameters (i.e., 𝛿,𝛿₁,𝛿₂), whose optimal values are often unknown. Hence, in Figures 2 and 3 we show the effect on the error of varying 𝛿, 𝛿₁, and 𝛿₂ on a fixed mesh ℎ=1/64 for the regular and multiscale enrichment approximations, respectively. An interesting thing can be observed that these two methods have an analogous convergence pattern with respect to the parameter 𝛿 and 𝛿₁. Because they involve a similar stabilization term with respect to these two parameters, we note that the errors can become large with some values of the stabilization parameters.

Figure 2 

Effect of varying 𝛿 for the regular method.

(a)

(b)

Figure 3 

Effect of varying 𝛿1 (a) and 𝛿2 (b) for the multiscale enrichment method.

Results gotten from the penalty, regular, multiscale enrichment, local Gauss integration, and nonconforming local Gauss integration methods are presented in Tables 1–5, respectively. Here, we choose 𝜀=1.0𝑒−5, 𝛼=8, 𝛼₁=8, and 𝛼₂=12. Because they can deal with the considered problem well. From Tables 1, 2, 3, 4, and 5, we can see that these methods work well and keep the convergence rates just like the theoretical analysis except the multiscale enrichment method. Meanwhile, it can be seen that the penalty method requires the least CPU-time, which validates the analysis in Remark 3.3. As expected, we have an interesting observation that the error of the nonconforming local Gauss integration method is better than the conforming version, which is not surprising since the degrees of freedom of the nonconforming method are nearly three times than that of conforming one on uniform mesh (see Figure 1). Hence, it is natural that the nonconforming local Gauss integration method is more accurate and costs more CPU-time.

Table 1 

Results get from the penalty method with 𝜀=1.0𝑒−5.

Table 2 

Results get from the regular method with 𝛼=8.

Table 3 

Results get from the multiscale enrichment method with 𝛼₁=8 and 𝛼₂=12.

Table 4 

Results get from the local Gauss integration method.

Table 5 

Results get from the local Gauss integration method with the nonconforming element.

Besides, to show the stability and efficiency of these methods for the considered problem, we present the velocity streamlines and the pressure contours with ℎ=1/64 in Figures 5, 6, 7, 8, and 9. Meanwhile, we present the results by the standard Galerkin method (𝑃₂-𝑃₁ element) computed on a very fine mesh (6742 grid points) for the purpose of comparison in Figure 4. From Figures 5(a)–9(a), five resolved vortices are captured, which is consistent with that in Figure 4(a). For the pressure, the penalty method is divergent from Figure 5(b), although it costs the least time. From Figures 6(b)-9(b), the nonconforming local Gauss integration method shows the best numerical stability.

(a)

(b)

Figure 4 

Velocity streamlines (a) and pressure level lines (b) for the standard Galerkin method (𝑃2-𝑃1 element).

(a)

(b)

Figure 5 

Velocity streamlines (a) and pressure level lines (b) for the penalty method with 𝜀=1.0𝑒−5.

(a)

(b)

Figure 6 

Velocity streamlines (a) and pressure level lines (b) for the regular method with 𝛼=8.

(a)

(b)

Figure 7 

Velocity streamlines (a) and pressure level lines (b) for the multiscale enrichment method with 𝛼1=8 and 𝛼2=12.

(a)

(b)

Figure 8 

Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method.

(a)

(b)

Figure 9 

Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method with the nonconforming element.

5. Conclusions

In this paper we have presented several stabilized mixed finite element methods in solving the Stokes eigenvalue problem based on the lowest equal-order finite element space pair. By being compared numerically, we get some conclusions as follows.(i)The stability and efficiency of the penalty method depend on the penalty parameter. The smaller this parameter, the more stable the method. However, if this parameter is too small, the condition number of the system matrix arising from this method will become too large to solve. (ii)The performance of the regular and multiscale enrichment method heavily depends on the choice of the stabilization parameters, which is a difficult task in reality. Meanwhile, a poor choice of these stabilization parameters can also lead to serious deterioration in the convergence rates. (iii)The local integration method is free of stabilization parameters and shows numerically the best performance among the methods considered for the given problem. (iv)From Tables 4 and 5, we have an interesting observation that the value of 𝜆_ℎ by conforming method becomes small to converge to the exact solution and the one by nonconforming method becomes large to converge to the exact solution. We hold that it can obtain more accurate numerical solution by Lagrange interpolation between conforming and nonconforming results based on the same degrees of freedom of these two methods. It may get superconvergence result on this triangular mesh.

Acknowledgments

The authors would like to thank the editor and reviewers for their valuable comments and suggestions which helped us improve the results of this paper. This work is in parts supported by the NSF of China (no. 10901131, no. 10971166), the National High Technology Research and Development Program of China (863 Program, no. 2009AA01A135), the China Postdoctoral Science Foundation (no. 200801448, no. 20070421155), and the Natural Science Foundation of Xinjiang Province (no. 2010211B04).