A new nonconforming mixed finite element scheme for second order eigenvalue problem

Applied Mathematics and Computation 273 (2016) 842–855

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A new nonconforming mixed finite element scheme for second order eigenvalue problem Dongyang Shi∗, Lele Wang, Xin Liao School of Mathematics and Statistics, Zhengzhou University, Zhengzhou 450001, China

a r t i c l e

i n f o

Keywords: Eigenvalue problem Nonconforming mixed element Superconvergence and extrapolation

a b s t r a c t A new nonconforming mixed finite element method (MFEM for short) is established for the Laplace eigenvalue problem. Firstly, the optimal order error estimates for both the eigenvalue and eigenpair (the original variable u and the auxiliary variable p = ∇ u) are deduced, the lower bound of eigenvalue is estimated simultaneously. Then, by use of the special property of the nonconforming EQ1rot element (the consistency error is of order O(h2 ) in broken H1 -norm, which is one order higher than its interpolation error), the techniques of integral identity and interpolation postprocessing, we derive the superclose and superconvergence results of order O(h2 ) for u in broken H1 -norm and p in L2 -norm. Furthermore, with the help of asymptotic expansions, the extrapolation solution of order O(h3 ) for eigenvalue is obtained. Finally, some numerical results are presented to validate our theoretical analysis. © 2015 Elsevier Inc. All rights reserved.

1. Introduction In this paper, we discuss the following eigenvalue problem



−u = λu, in , u = 0, on ∂ ,

(1)

where  ⊂ R2 is a bounded convex polygonal domain with Lipschitz continuous boundary ∂ . Eigenvalue problems play a very important role in mathematical physics and engineering technology which have won extensive attention by scholars. For the second order elliptic eigenvalue problem, the standard FEMs were studied in [1–15]. Where [1–5] estimated the lower and upper bounds for the eigenvalues; [6–8] derived asymptotic expansions and extrapolations for the eigenvalues; Dari et al. [9] and the authors in [10–12] proposed posteriori error estimates and adaptive FEM, respectively; Yang and Bi [13] and Xie [14] discussed the two-grid and multigrid methods, respectively; Lin and Xie [15] introduced a multilevel correction scheme. The standard MFEMs were researched in [16–22]. In which Mercier et al. [16] analyzed the error estimates for eigenvalue and eigenvectors by mixed and hybrid FEMs; Gardini and Lin and Xie [17,18] obtained the superconvergence results of order O(h2 ) for the eigenvector u in L2 -norm; Lin and Xie [19] investigated the error estimates for eigenvalue and eigenvectors and derived the asymptotic expansions and extrapolation for eigenvalue; Durán et al. and Jia et al. [20,21] discussed a posteriori error estimates; Chen et al. [22] gave the two-grid method for MFEMs. However, most of the above MFEMs for elliptic eigenvalue problem focused on the conforming elements without consideration on the superconvergence and extrapolation for nonconforming elements. ∗

Corresponding author. Tel.: +86 13838117268. E-mail addresses: [email protected], [email protected] (D. Shi), [email protected] (L. Wang), [email protected] (X. Liao).

http://dx.doi.org/10.1016/j.amc.2015.10.064 0096-3003/© 2015 Elsevier Inc. All rights reserved.

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

843

Recently, a new mixed variational form for second elliptic problem was constructed in [23,24]. Compared with the standard MFE scheme, the new mixed formulation avoids the divergence space H (div) which makes the theoretical analysis simpler; and the LBB condition is automatically satisfied when the gradient of approximation space for the original variable is included in the approximation space for the auxiliary variable, which leads to easier choice for the mixed approximation spaces. Subsequently, this method was further applied to different equations (see [25–29]). For problem (1), Weng et al. [26] provided the convergence for eigenvalue and eigenpair in two-grid new MFE scheme via triangular elements, but it did not cover the superconvergence and extrapolation; Shi et al. [27] studied the new formulation with conforming MFEM, and obtained the error estimates for the eigenpair and extrapolation solution for the eigenvalue without analysis of the superconvergence and numerical experiments. The main purpose of this article is to apply the new MFE scheme to problem (1) with nonconforming elements. The outline of this paper is organized as follows: In Sections 2 and 3, the optimal error estimates for eigenvalue and eigenpair and the approximation from below for eigenvalue are obtained. In Section 4, the superclose and superconvergence results of order O(h2 )  in L2 -norm based on the special property are derived for both the original variable u in broken H1 -norm and auxiliary variable p of the nonconforming EQ1rot element (when u ∈ H3 (), the consistency error is of order O(h2 ) which is one order higher than that of the interpolation error O(h)) and the techniques of integral identity and interpolation postprocessing. In Section 5, the extrapolation solution with accuracy O(h3 ) for eigenvalue is received with the help of the asymptotic expansions. In Section 6, some numerical results are provided to illustrate the validity of the theoretical analysis and effectiveness of the proposed method. 2. New MFE scheme and convergence analysis Let  be a polygon domain with edges parallel to the coordinate axes, Th be a rectangular subdivision of  which need not to satisfy the regular condition [30]. For all K ∈ Th , K = [xK − hxK , xK + hxK ] × [yK − hyK , yK + hyK ], we denote the barycenter of element K by (xK , yK ), the four vertices and four sides are zi , li = zi zi+1 (mod4)(i = 1, 2, 3, 4), respectively. hK = max{hxK , hyK }, h = maxK∈Th hK . The associated MFE spaces Mh (nonconforming EQ1rot element space [31–34]) and Hh (the lowest order Raviart–Thomas element space [30]) are defined by



 Mh = {vh : vh |K ∈ span{1, x, y, x2 , y2 }, ∀ K ∈ Th , F [vh ]ds = 0, F ⊂ ∂ K }, h |K ∈ span{1, x} × span{1, y}, ∀ K ∈ Th }, h = (q1h , q2h ) : q Hh = Hh1 × Hh2 = {q

where [vh ] stands for the jump of vh across the boundary F and [vh ] = vh if F ⊂ ∂ .  1 Then we denote the norms on Mh and Hh as  · h = ( K∈T | · |21,K ) 2 and  · 0 , respectively. h  ∈ Hh ,  ∈ H → h p Similar to [25], the corresponding interpolation operators are defined as Ih : v ∈ M → Ih v ∈ Mh and h : p Ih |K = IK , h |K = K satisfying



li

(v − IK v)ds = 0,



K

(v − IK v)dxdy = 0,



li

( p − K p) · n i ds = 0, ∀ K ∈ Th ,

 i is the unit outer normal vector on li (i = 1, 2, 3, 4). here M = H01 (), H = (L2 ())2 , and n  , u) ∈ R × H × M,  u 0 = 1, such that  = ∇ u, then the new mixed variational formulation for (1) is to find (λ, p Let p



( p, q) = (∇ u, q), ∀ q ∈ H, ( p, ∇v) = λ(u, v), ∀ v ∈ M,

(2)

 where (∗, ) =  ∗ · dxdy. There exist eigenvalues λi (i ∈ N) for Eq. (2) as follows [35]

0 < λ1 ≤ λ2 ≤ · · · ≤ λi ≤ · · · , lim λi = ∞ i→∞

and associated eigenfunctions

( p 1 , u1 ), ( p 2 , u2 ), · · · , ( p i , ui ), · · · ,  i = ∇ ui . where (ui , u j ) = δi j (δ ij denotes the Kronecker symbol) and p  h , uh ) ∈ R × Hh × Mh ,  uh 0 = 1, such that The new nonconforming mixed discrete scheme for (2) is: find (λh , p



∀ qh ∈ Hh , ( p h , qh ) = (∇ uh , qh )h , ( p h , ∇vh )h = λh (uh , vh ), ∀ vh ∈ Mh ,

(3)

( p, q) = (∇ u, q), ∀ q ∈ H, ( p, ∇v) = ( f, v), ∀ v ∈ M,

(4)

  where (∗, )h = K K ∗ · dxdy.  , u) ∈ H × M, such that Now we introduce the boundary value problem corresponding to (2): for all f ∈ G, find ( p



where G = L2 (), and the regularity of solutions shows that

 u 2 +  p 1 ≤ C  f 0 .

(5)

844

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

Throughout this article, C denotes a positive constant whose value may be different in different places but remains independent of the mesh parameter h.  h , uh ) ∈ Hh × Mh , satisfying The discrete scheme for (4) is: find ( p

 ( p h , qh ) = (∇ uh , qh )h , ∀ qh ∈ Hh , ( p h , ∇vh )h = ( f, vh ), ∀ vh ∈ Mh .

(6)

The existence and uniqueness of the solution for problem (4) and (6) can be found in [23]. On the one hand, similar to [16] we define the bounded linear operator

T : G → M, T f = u, , S : G → H, S f = p

Th : G → Mh , Th f = uh , h. Sh : G → Hh , Sh f = p

Then (4) and (6) can be rewritten as

 (S f, q) − (∇ T f, q) = 0, ∀ q ∈ H, (S f, ∇v) = ( f, v), ∀ v ∈ M.  (Sh f, qh ) − (∇ Th f, qh )h = 0, ∀ qh ∈ Hh , (Sh f, ∇vh )h = ( f, vh ), ∀ vh ∈ Mh .

(7)

(8)

The corresponding operator equations for (2) and (3) are

λTu = u, , S(λu) = p

λh Th uh = uh , h. Sh (λh uh ) = p

(9) (10)

) ∈ Mh × Hh and On the other hand, like [16], we define the projection operators Ah and Bh satisfy (Ah u, Bh p

 (Bh p, qh ) − (∇ Ah u, qh )h = 0, ∀ qh ∈ Hh , (Bh p, ∇vh )h = λ(u, vh ), ∀ vh ∈ Mh .

(11)

Now we provide the following lemmas which will play an important role in the error estimates. Lemma 2.1. [25, 36] For all vh ∈ Mh , qh ∈ Hh , we have

(∇(u − Ih u), ∇vh )h = 0,

(12)

(∇(u − Ih u), qh )h = 0,

(13)

 vh 0 ≤ C vh h ,

(14)

( p − h p, qh ) ≤ Ch2 | p |2  qh 0 ,

(15)

 K

∂K

 Ch | u |2 vh h , u ∈ H 2 (), ∂u vh ds ≤ ∂n Ch2 | u |3 vh h , u ∈ H 3 ().

(16)

 ∈ (H 1 ())2 , then we arrive at Lemma 2.2. Assume that u ∈ H 2 (), p

 S − Sh 0 ≤ Ch, T − Th h ≤ Ch,

(17)

 T − Th 0 ≤ Ch2 ,

(18)

where 

(S−Sh ) f 0 (T −Th ) f h S − Sh 0 = sup f ∈G , T − Th h = sup f ∈G .  f 0  f 0

Proof. For the boundary value problem (4) and (6), similar to the proof in [23], we get

u − uh h +  p − p h 0 ≤ Ch(  u 2 +  p 1 ) ≤ Ch  f 0 . Then there hold

T f − Th f h = u − uh h ≤ Ch  f 0 ,  S f − Sh f 0 = p − p h 0 ≤ Ch  f 0 , which lead to (17).

(19)

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

845

 , ϕ) ∈ H × M be the solution of the following auxiliary problem: for g ∈ G, In order to prove (18), let (ψ



(ψ , ∇v) = (g, v), ∀ v ∈ M, (ψ , q) = (∇ϕ , q), ∀ q ∈ H,

  +ϕ  ≤C g . then there holds  ψ 1 2 0 Define the projection operator Rh : H → Hh satisfying

(ψ − Rh ψ , qh ) = 0,  ψ − Rh ψ 0 ≤ Ch  ψ 1 , ∀ qh ∈ Hh . Then we have

(g, u − uh ) = (ψ , ∇(u − uh ))h +

 K

=

∂K

(20)

ψ · uh · n ds

(ψ − Rh ψ , ∇(u − Ih u))h + (ψ − Rh ψ , ∇(Ih u − uh ))h + (Rh ψ , p − p h ) + 

=

(ψ − Rh ψ , ∇(u − Ih u))h − (ψ − Rh ψ , p − p h ) + (∇ϕ , p − p h ) +

=

(ψ − Rh ψ , ∇(u − Ih u))h − (ψ − Rh ψ , p − p h ) + (∇(ϕ − Ih ϕ), p − p h )h

K

+



∂K

K

 ds +  · (Ih ϕ − ϕ) · n p

 K

∂K

∂K

 K

∂K

ψ · uh · n ds

ψ · uh · n ds

ψ · (uh − u) · n ds.

By use of the interpolation theory, (5), (16 a), (19) and (20), we derive

(g, u − uh ) ≤ Ch2  ψ 1  u 2 +Ch2  ψ 1  p 1 +Ch2  ϕ 2  p 1 +Ch2  ψ 1  u 2 ≤ Ch2  f 0  g 0 . Therefore

 u − uh 0 = sup g∈G

(u − uh , g) ≤ Ch2  f 0 ,  g 0

which leads to

 T f − Th f 0 ≤ Ch2  f 0 . The proof is completed.  Lemma 2.3. For f ∈ G, there holds

((T − Th ) f, f ) = −(Sh f − S f, Sh f − S f ) − 2(∇(T f − Th f ), Sh f − S f )h + 2

 K

∂K

 ds. S f · Th f · n

(21)

Proof. Let vh = Th f, according to the second formulas of (7) and (8), we have

(S f − Sh f, ∇ Th f )h =

 K

∂K

 ds. S f · Th f · n

(22)

Due to (7) and (22), we deduce that

(Sh f, ∇(Th f − T f ))h = (S f, ∇ Th f )h − (S f − Sh f, ∇ Th f )h − (Sh f − S f, ∇ T f ) − (S f, ∇ T f ) = (S f, ∇(Th f − T f ))h + (S f − Sh f, ∇ T f ) − = ( f, Th f − T f ) + (S f − Sh f, ∇ T f ).

 K

∂K

 ds S f · Th f · n

Consequently,

( f, T f − Th f ) = (S f − Sh f, ∇ T f ) − (Sh f, ∇(Th f − T f ))h  A1 − A2 . Now we start to estimate A2 . Using (7) and (8), it is easy to see that

A2 = (Sh f, ∇ Th f )h − (Sh f, ∇ T f ) = (Sh f, Sh f − S f ) = (Sh f − S f, Sh f − S f ) + (∇ T f, Sh f − S f ). Substituting A2 into (23), yields

( f, T f − Th f ) = (S f − Sh f, ∇ T f ) − (Sh f − S f, Sh f − S f ) − (∇ T f, Sh f − S f ) = 2(S f − Sh f, ∇(T f − Th f ))h + 2(S f − Sh f, ∇ Th f )h − (Sh f − S f, Sh f − S f ) = 2(S f − Sh f, ∇(T f − Th f ))h + 2

 K

∂K

 ds − (Sh f − S f, Sh f − S f ). S f · Th f · n

(23)

846

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

The proof is completed.  Then we give the following error estimates. ) and (λh , uh , p  h ) be the solutions of (2) and (3), respectively, there hold Theorem 2.1. Let (λ, u, p

 u − uh 0 ≤ Ch2 ,

(24)

|λ − λh | ≤ Ch2 ,

(25)

 p − p h 0 ≤ Ch.

(26)

Proof. Firstly, we let uh =

∞

j=1

(uh , u j )u j , Euh =



λ j =λ (uh , u j )u j . Obviously,

Euh

(λ,

Euh 0 ,

Eu

Euh ∇ Eu  ) is the solution of (2). For h 0

simplicity, let u = Eu h , from (9) we have h 0

(Tuh − Th uh , u) = (T − 1/λh )

∞ 

((uh , u j )u j , u) =

j=1

∞ 

(1/λ j − 1/λh )((uh , u j )u j , u) = (1/λ − 1/λh )(uh , u)

(27)

j=1

and

 uh − u 20 = 2 − 2(uh , u) = 2 − 2  Euh 0 = 2(1 −



1−  uh − Euh

20 ) ≤ uh − Euh 20 (1+  uh − Euh 20 ),

(28)

1

here the inequality (1 − a) 2 ≥ 1 − 12 a(1 + a) (for all a ∈ R and |a| < 1) is used in the last step of (28). Let d(λ) = minλ j =λ | λ1 − λ1 |, i.e., | λ1 − λ1 |≥ d(λ)(λ j = λ). According to [16], we have j

h

 Euh − uh 20 = 



h

(uh , u j )u j 20 =

λ j =λ

≤

∞  j=1

(

j



(uh , u j )2 ≤

λ j =λ

1

λh

−

1

λj

)2 (uh , u j )2 /d2 (λ) ≤



(

λ j =λ

1

λh

−

1

λj

)2 (uh , u j )2 /d2 (λ)

 Tuh − Th uh 20 . d2 (λ)

Therefore,

 u − uh 0 ≤ C  Tuh − Th uh 0 = C  (T − Th )u + (T − Th )(uh − u) 0 ≤ C  (T − Th )u 0 +C  T − Th 0  uh − u 0 , which together with (18), yields

 u − uh 0 ≤ C  (T − Th )u 0 ≤ Ch2 , this is the desired result (24). Secondly, using (27), we obtain

λh − λ =

λh λ λλ λλ (Tu − T u , u) = h ((T − Th )u, u) + h ((T − Th )(uh − u), u). (u, uh ) h h h (u, uh ) (u, uh )

Choosing f = λu in (21) and substituting it into (29), we find that

λh − λ =

λh λ

(u, uh ) +

[−(Sh u − Su, Sh u − Su) − 2(∇(Tu − Th u), Sh u − Su) + 2

 K

∂K

 ds] Su · Th u · n

λh λ ((T − Th )(uh − u), u). (u, uh )

Then, the result (25) follows from the Schwarz inequality, (17) and (18). Finally, according to (10), there holds

 p − p h 0 =  S(λu) − Sh (λh uh ) 0 = S(λu) − Sh (λu) + Sh (λu) − Sh (λh uh ) 0 ≤ C  S − Sh 0 +C  λu − λh uh 0 ≤ C (  S − Sh 0 +  T − Th 0 +|λh − λ|). By use of Lemma 2.2 and (25) we get (26), which ends the proof. 

(29)

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

847

3. The approximation of eigenvalue from below  h , uh ) be the solutions of (2) and (3), respectively, for all vh ∈ Mh , we have the expansion as follows , u), (λh , p Lemma 3.1. Let (λ, p

λ − λh = p − p h 20 +2( p, p h ) − 2( p h , ∇vh )h − λh  uh − vh 20 +λh (  vh 20 −  u 20 ).

(30)

Proof. Note that

 p − p h 20 = p 20 +  p h 20 −2( p, p h ) = λ + λh − 2( p, p h ).

(31)

 h , ∇vh )h to both sides of (31), it follows from (3) that Adding −2( p

 h , ∇vh )h +  p −p h −2( p

20 = λ + λh − 2( p, p h ) − 2λh (uh , vh ) , p  h ) + λh  uh − vh 20 −λh (  vh 20 −  u 20 ) − 2λh . = λ + λh − 2( p

Thus

λ − λh = p − p h 20 +2( p, p h ) − λh  uh − vh 20 +λh (  vh 20 −  u 20 ) − 2( p h , ∇vh )h , which ends the proof.  Theorem 3.1. If u ∈ H 2 ()



H01 (), we have the lower approximation below

λh ≤ λ.

(32)

Proof. On the one hand, we define the piecewise constant projection operator P0 as

P0 v|K =

1 |K |



K

v dxdy, ∀ K ∈ Th ,

which leads to

v − P0 v0 ≤ Ch|v|1 , ∀ v ∈ H 1 ().

(33)

On the other hand, let vh = Ih u in (30), then from (13) and (33), the terms on the right hand of (30) can be estimated as

 h )h = 0, , p  h ) − 2( p  h , ∇ Ih u)h = 2(∇ u − ∇ Ih u, p 2( p

 uh − Ih u 20 ≤ uh − u 20 +  u − Ih u 20 ≤ Ch4 + Ch4  u 22 , | Ih u 20 −  u 20 |=|





(Ih u − u)(Ih u + u)dxdy |=|





(Ih u − u)[(Ih u + u) − P0 (Ih u + u)]dxdy |≤ Ch3  u 22 Ih u + uh .

With the help of (26) and (30), we get the desired result (32).  4. Superclose and superconvergence analysis  h ) be the solutions of (2) and (3), respectively. If u ∈ H 3 () ∩ H01 (), p  ∈ (H 2 ())2 , we have  ), (λh , uh , p Theorem 4.1. Let (λ, u, p the following superclose conclusions

uh − Ih uh ≤ Ch2 ,

(34)

 p h − h p 0 ≤ Ch2 .

(35)

Proof. The error equations are derived from (1) and (3) as follows



( p − p h , qh ) = (∇(u − uh ), qh )h ,   ( p − p h , ∇vh )h = K ∂ K ∂∂ un vh ds + (λu − λh uh , vh ).

 − h p  ) + (h p −p h)  δ + χ . −p h = ( p Let u − uh = (u − Ih u) + (Ih u − uh )  w + σ , p h = −∇σ in (36a) and vh = σ in (36b), then adding them, we have Taking q

(∇(w + σ ), ∇σ )h =

 K

∂K

∂u σ ds + (λ(u − uh ) + (λ − λh )uh , σ ). ∂n

By use of (12), (14), (16 b), (24) and (25), we arrive at

σ 2h ≤ Ch2 | u |3 σ h + Ch2  σ 0 +Ch2  σ 0  uh 0 ≤ Ch2 σ h , which leads to (34). h = χ in (36a), we get On the other hand, choosing q

(δ + χ , χ) = (∇(w + σ ), χ )h = (∇σ , χ )h .

(36)

848

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

˜ Fig. 1. Big element K. Table 1 Numerical approximation λh on regular meshes. m×n

4×4

8×8

16 × 16

32 × 32

λh

18.8353

19.4931

19.6762

19.7233



Applying (15), there holds

 χ 20 = (∇σ , χ)h − (δ, χ ) ≤ Ch2  χ 0 +Ch2 | p |2  χ 0 ≤ Ch2  χ 0 , which ends the proof.  For the purpose of obtaining the superconvergence results, by applying the method introduced in [6], we combine the ad ˜ i.e., K˜ = 4i=1 Ki (see Fig. 1). The corresponding subdivision is noted by jacent four elements K1 , K2 , K3 , K4 into a big element K, T2h . Constructing the interpolation postprocessing operators I2h and 2h on K˜ as in [6,37,38]:



 I2h u |K˜ ∈ P2 (K˜ ), i = 1, 2, 3, 4, ∀ K˜ ∈ T2h , Li (I2h u − u)ds = 0,  ( I u − u ) dxdy = 0, ( I u K1 K3 2h K2 K4 2h − u)dxdy = 0,

 2h p |K˜ ∈ Q1,1 (K˜ ) × Q1,1 (K˜ ), ∀ K˜ ∈ T2h ,   1 1 1 2 2 2 i = 3, 4, 7, 8, i = 1, 2, 5, 6, li (2h p − p )ds = 0, li (2h p − p )ds = 0, where P2 (K˜ ) denotes the space of polynomial on K˜ with degree less than or equal to 2, Q1, 1 defines the bilinear element space on ˜ Li (i = 1, 2, 3, 4) are the four lines of K. ˜ K, From [6,25] we know that I2h and 2h satisfy the following properties

I2h Ih u = I2h u,

I2h u − uh ≤ Ch2 | u |3 , I2h vh h ≤ C vh h , ∀ vh ∈ Mh ,

(37)

and

2h h p = 2h p,

 2h p − p 0 ≤ Ch2 | p |2 ,  2h qh 0 ≤ C  qh 0 , ∀ qh ∈ Hh .

(38)

Theorem 4.2. Under the assumptions of Theorem 4.1, we can get the superconvergence results as follows

u − I2h uh h ≤ Ch2 ,

(39)

 p − 2h p h 0 ≤ Ch2 .

(40)

Proof. By triangle inequality, Lemma 4.1 and (37), there holds

u − I2h uh h = u − I2h Ih u + I2h Ih u − I2h uh h ≤ u − I2h uh + C Ih u − uh h ≤ Ch2 , similarly, we can derive (40). The proof is completed.  Remark 4.1. When the standard MFE scheme is employed (see [17–19]), we can only obtain the convergence and superconvergence results for u in L2 -norm instead of broken H1 -norm. Remark 4.2. From the analysis of Theorems 3.1–4.1, we can see that Lemma 2.3 is crucial to the superclose and superconvergence results. It can be checked that (12)–(16b) are also hold true when Mh is replaced by the Q1rot element space [39] and Hh remains unchanged on the square meshes. Furthermore, if Mh and Hh are selected as rectangular constrained Q1rot element space [40] and the piecewise constant space respectively, it has been proved in [29] that (12) and (13) can be estimated as (∇(u − Ih u), ∇vh )h = h 0 , respectively, the conclusions in Theorems 3.1–4.1 are also valid in h )h = O(h2 ) | u |3  q O(h2 ) | u |3 vh h and (∇(u − Ih u), q this case. However, for the famous rectangular Wilson element, (16b) cannot be obtained because the consistency error is only of order O(h) and cannot be improved any more (see [41]). On the other hand, if Mh is chosen as the quasi-Wilson element space [42–44], though the consistency error is of order O(h2 ), the choice of appropriate space Hh that matches Mh is difficult. This further indicates that how to select the proper approximation spaces to obtain the superclose and superconvergence results is not an easy work.

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

849

Table 2 Numerical approximation λh on anisotropic meshes. m×n

4 × 40

8 × 80

16 × 160

32 × 320

λh

19.2604

19.6133

19.7073

19.7312



Table 3 Numerical results of λ on regular meshes. m×n

|λ − λh |

order

|λ − λ˜ h |

order

4×4 8×8 16 × 16 32 × 32

0.9038 0.2460 0.06293 0.01582

– 1.8768 1.9673 1.9917

0.2933 0.02684 0.001879 0.0001210

– 3.4501 3.8363 3.9567

Table 4 Numerical results of λ on anisotropic meshes. m×n

|λ − λh |

order

|λ − λ˜ h |

order

4 × 40 8 × 80 16 × 160 32 × 320

0.4787 0.1258 0.03188 0.007997

– 1.9275 1.9808 1.9951

0.09346 0.008213 0.0005612 3.58917e−05

– 3.5083 3.8711 3.9670

Table 5 Numerical results of u on regular meshes. m×n

 u − uh 0

order

u − uh h

order

uh − Ih uh

order

u − I2h uh h

order

4×4 8×8 16 × 16 32 × 32

0.05071 0.01281 3.2104e−3 8.0304e−4

– 1.9849 1.9967 1.9992

1.0006 0.5026 0.2516 0.1258

– 0.9933 0.9978 0.9994

0.05519 0.007368 9.3801e−04 1.1780e−04

– 2.9051 2.9736 2.9932

0.8589 0.2249 0.05688 0.01426

– 1.9330 1.9836 1.9959

Table 6 Numerical results of u on anisotropic meshes. m×n

 u − uh 0

order

u − uh h

order

uh − Ih uh

order

u − I2h uh h

order

4 × 40 8 × 80 16 × 160 32 × 320

0.005957 0.001362 3.2664e−4 8.0663e−5

– 2.1286 2.0602 2.0177

0.7113 0.3572 0.1788 0.08946

– 0.9936 0.9978 0.9994

0.04070 0.005269 6.6519e−04 8.3359e−05

– 2.9495 2.9858 2.9963

0.1100 0.02431 0.005808 0.001433

– 2.1782 2.0657 2.0184

Table 7  on regular meshes. Numerical results of p m×n

 p − p h 0

order

 p h − h p 0

order

 p − 2h p h 0

order

4×4 8×8 16 × 16 32 × 32

0.9993 0.5025 0.2516 0.1258

– 0.9915 0.9977 0.9994

0.2264 0.05708 0.01427 0.003568

– 1.9882 1.9995 2.0000

0.6878 0.1713 0.04282 0.01070

– 2.0047 2.0007 2.0001

5. Asymptotic expansions and extrapolation  h ) be the solutions of (2) and (3), respectively, there holds  ), (λh , uh , p Lemma 5.1. Let (λ, u, p

λh − λ = λh (u − Ih u, u¯ h ) + ( p¯ h , p − h p) + (∇ u¯ h , h p − p)h +

 K

u

 p

where u¯ h = (u,uh ) , p¯ h = ( p¯ 1h , p¯ 2h ) = (u,uh ) . h h Proof. The definition of u¯ h gives that (u, u¯ h ) = 1, therefore

λh = λh (u, u¯ h ) = λh (Ah u, u¯ h ) + λh (u − Ah u, u¯ h )  I1 + I2 .

∂K

∂u · u¯ ds, ∂n h

(41)

850

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855 Table 8  on anisotropic meshes. Numerical results of p m×n

 p − p h 0

order

 p h − h p 0

order

 p − 2h p h 0

order

4 × 40 8 × 80 16 × 160 32 × 320

0.7102 0.3571 0.1788 0.0894

– 0.9916 0.9977 0.9994

0.1583 0.04020 0.01008 0.002522

– 1.9772 1.9954 1.9989

0.4861 0.1211 0.03028 0.007570

– 2.0039 2.0006 2.0001

2

2

5 1.5

1.5

1

1

5 0.5

0.5

0 1

0 1 0.5 0 0

0.2

0.4

0.6

0.8

1 0.5 0 0

0.2

0.4

0.6

0.8

1

Fig. 2. The exact solution u (left) and approximate solution uh (right) on regular meshes.

10

10

5

5

0

0

−5

−5

−10 1

−10 1 0.5 0 0

0.2 02

0.4

0.6

0.8

1 0.5 0 0

0 2 0.2

0.4

0.6

0.8

1

Fig. 3. The exact solution p1 (left) and approximate solution p1h (right) on regular meshes.

From (3) and (11), we have

I1 =

1 1 1 λh (A u, u ) = ( p , ∇ Ah u)h = (∇ uh , Bh p)h = λ(u, uh ) = λ(u, u¯ h ) (u, uh ) h h (u, uh ) h (u, uh ) (u, uh )

(42)

and

λh (u − Ih u, u¯ h ) + λh (Ih u − Ah u, u¯ h ) λh (u − Ih u, u¯ h ) + ( p¯ h , ∇(Ih u − u))h + ( p¯ h , ∇(u − Ah u))h = λh (u − Ih u, u¯ h ) + ( p¯ h , ∇ u) − ( p¯ h , ∇ Ah u)h )  − Bh p = λh (u − Ih u, u¯ h ) + ( p¯ h , p  ) + (∇ u¯ h , h p  − Bh p )h  − h p = λh (u − Ih u, u¯ h ) + ( p¯ h , p  − h p  )h − (∇ u¯ h , Bh p  ) + (∇ u¯ h , h p −p  )h + (∇ u¯ h , p  )h = λh (u − Ih u, u¯ h ) + ( p¯ h , p   ∂u  − h p  ) + (∇ u¯ h , h p −p  )h + · u¯ h ds. = λh (u − Ih u, u¯ h ) + ( p¯ h , p ∂K ∂ n K

I2 = =

(43)

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

10

10

5

5

0

0

−5

−5

−10 1

−10 1 0.5 0.2 0 2

0 0

0.6

0.4

0.8

851

1 0.5 0 0

0.2 0 2

0.4

0.6

0.8

1

Fig. 4. The exact solution p2 (left) and approximate solution p2h (right) on regular meshes.

2

2

1.5

1.5

1

1

0.5

0.5

0 1

0 1 0.5 0.2 02

0 0

0.4

0.6

0.8

1 0.5 0 0

0 2 0.2

0.4

0.6

0.8

1

Fig. 5. The exact solution u (left) and approximate solution uh (right) on anisotropic meshes.

The desired result comes from the combination of (42) and (43).  Lemma 5.2. [6] When the mesh is uniform subdivision, i.e., hxK ≡ h1 , hyK ≡ h2 , there hold

(u − Ih u, u¯ h ) = O(h4 )  u 4 | u¯ h |h ,  K

∂K

(44)

  h2  h2  ∂u · u¯ h ds = 1 uxxy u¯ hy dxdy + 2 uxyy u¯ hx dxdy + O(h3 ) | u |4 u¯ h h , ∂n 3 3 K K K K

( p − h p, p¯ h ) = −

  h21  h2   |3  p¯ h p1xx p¯ 1h dxdy − 2 p2yy p¯ 2h dxdy + O(h3 ) | p 3 3 K K K

( p − h p, ∇ u¯ h )h = −

0 ,

(45)

(46)

K

  h21  h2   |3 u¯ h h . p1xx u¯ hx dxdy − 2 p2yy u¯ hy dxdy + O(h3 ) | p 3 3 K K K

(47)

K

Lemma 5.3. Under the assumptions of Lemma 5.2, we have

λh − λ =

  h21  h2  uxxy uy dxdy + 2 uxyy ux dxdy + O(h3 ). 3 3 K K K

(48)

K

Proof. [45] has proved that

 u¯ h − uh 0 ≤ Ch2 ,  p¯ h − p h 0 ≤ Ch2 .

(49)

852

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

10

10

5

5

0

0

−5

−5

−10 1

−10 1 0.5 0 0

0.2

0.4

0.6

0.8

1 0.5 0 0

0.2

0.4

0.6

0.8

1

Fig. 6. The exact solution p1 (left) and approximate solution p1h (right) on anisotropic meshes.

10

10

5

5

0

0

−5

−5

−10 1

−10 1 0.5 0 0

0.2

0.4

0.6

0.8

1 0.5 0 0

0.2

0.4

0.6

0.8

Fig. 7. The exact solution p2 (left) and approximate solution p2h (right) on anisotropic meshes.

From Lemmas 5.1–5.2 and (49), we obtain

    h21  h2  h2  h2  p1xx p¯ 1h dxdy − 2 p2yy p¯ 2h dxdy + 1 p1xx u¯ hx dxdy + 2 p2yy u¯ hy dxdy 3 3 3 3 K K K K K K K K   h22  h21  uxxy u¯ hy dxdy + uxyy u¯ hx dxdy + O(h3 ) + 3 3 K K K K   h2  h2  = f1 (h) − 1 p1xx [( p¯ 1h − p1h ) + ( p1h − p1 )]dxdy − 2 p2yy [( p¯ 2h − p2h ) + ( p2h − p2 )]dxdy 3 3 K K K K   h2  h2  + 1 p1xx [(u¯ hx − uhx ) + (uhx − ux )]dxdy + 2 p2yy [(u¯ hy − uhy ) + (uhy − uy )]dxdy 3 3 K K K K   h22  h21  uxxy [(u¯ hy − uhy ) + (uhy − uy )]dxdy + uxyy [(u¯ hx − uhx ) + (uhx − ux )]dxdy + O(h3 ) + 3 3 K K

λh − λ = −

K

K

= f1 (h) + O(h3 ), where

    h21  h2  h2  h2  p1xx p1 dxdy − 2 p2yy p2 dxdy + 1 p1xx ux dxdy + 2 p2yy uy dxdy 3 3 3 3 K K K K K K K K     h22  h21  h22  h21  uxxy uy dxdy + uxyy ux dxdy = uxxy uy dxdy + uxyy ux dxdy. + 3 3 3 3 K K K K

f1 (h) = −

K

K

K

K

1

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

853

0

0

10

10

2

h

4

h λ

2

h

4

h λ

−1

−1

1

10

1

10

λ

λ

2

2

−2

10

error

error

−2

10

−3

−3

10

−4

10

10

−4

10

−5

10

−5

2

3

10

4

10

10

5

10

10

2

10

number of elements

3

10

4

5

10

10

number of elements

Fig. 8. Error reduction results for λ on regular meshes(left) and anisotropic meshes (right). 0

10 0

h

10

h

h2

h2

h3 u

−1

h3 u

−1

10

1

10

1

u

u

2

2

u

3

error

error

u4

−3

−3

10

−4

10

10

−4

10

−5

10

3

10

u4

−2

10

u

−2

−5

2

3

10

4

10

10

5

10

10

2

10

number of elements

3

10

4

10

5

10

number of elements

Fig. 9. Error reduction results for u on regular meshes(left) and anisotropic meshes (right).

The proof is completed.  In order to construct the extrapolation solution, with a similar analysis as in [6], we divide each element K ∈ Th into four  h/2 ) ∈ R × Mh/2 × Hh/2 . uniform parts, and calculate eigenvalue pair (λh/2 , uh/2 , p  ∈ (H 3 ())2 , there holds Theorem 5.1. Assume that u ∈ H 4 () ∩ H01 (), p

λ˜ h − λ = O(h3 ),

(50)

˜ = 4λh/2 −λh . where λ h 3 Proof. From Lemma 5.3 we can get 4 f1 (h/2) = f1 (h), therefore

λ˜ h − λ =

4λh/2 −λh −3λ 3

=

4(λh/2 −λ)−(λh −λ) 3

=

4 f1 (h/2)− f1 (h) 3

+ O(h3 ) = O(h3 ).

The proof is completed.  Remark 5.1. The postprocessing method here is an idea to improve the convergence rates of numerical solutions. But how to couple it with other numerical methods such as adaptive and multigrid methods (see [12,15]) and how to extend the results of this paper to triangular meshes (including conforming and nonconforming elements) are the topics of our future study. Remark 5.2. The conclusions in this paper are also valid for anisotropic meshes. In next section, we will present some numerical results to validate the effectiveness of the new MFEM under anisotropic meshes.

854

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855 0

10

0

10

h

h

2

2

h p

h p

p

p

1

1

2

2

p

3

p

3

−1

10

−1

error

error

10

−2

−2

10

10

−3

10

−3

2

10

3

4

10

10

5

10

10

2

10

number of elements

3

10

4

10

5

10

number of elements

Fig. 10. Error reduction results for  p on regular meshes(left) and anisotropic meshes (right).

6. Numerical results Consider the problem (1) with  = [0, 1] × [0, 1]. The exact solutions are

u = 2 sin π x sin π y,

 = 2π ( cos π x sin π y, sin π x cos π y), p

λ = 2π 2 = 19.7392.

We divide the domain  into m × n rectangles (m and n are positive integers), and present the numerical results under regular meshes (n = m) and anisotropic meshes (n = 10m) respectively. From Tables 1 and 2 we find the numerical approximation is lower bound of the exact eigenvalue. Then, the convergence, superclose, superconvergence and extrapolation results are listed in Tables 3, 4,5, 6, 7, and 8. It is −p  h 0 are convergent at rate of O(h), |λ − λh |,  u − uh 0 , u − I2h uh h ,  p  h − h p  0 and  clearly that u − uh h and  p ˜ |  − 2h p  h 0 are convergent at rate of O(h2 ), which coincide with our theoretical analysis. Furthermore, uh − Ih uh and |λ − λ p h are convergent at rate of O(h3 ) and O(h4 ), respectively, which are one order higher than the theoretical results.  h on regular meshes and anisotropic meshes, respectively  and finite element solutions uh , p We plot the exact solutions u, p (see Figs. 2, 3, 4, 5, 6, and 7). At the same time, we describe the error reduction results in Figs. 8, 9, and 10, where λ1 and λ2 denote |λ − λh | and |λ − λ˜ h |, respectively; u1 , u2 , u3 and u4 represent  u − uh 0 , u − uh h , uh − Ih uh and u − I2h uh h respectively; −p  h 0 ,  p  h − h p  0 and  p  − 2h p  h 0 , respectively. p1 , p2 and p3 denote  p Acknowledgments The research was supported by the National Natural Science Foundation of China (grant no. 11271340,11101381). The authors would like to thank the referee for his valuable suggestions on our manuscript. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.amc.2015.10.059. References [1] Y.D. Yang, Z.M. Zhang, F.B. Lin, Eigenvalue approximation from below using nonconforming finite elements, Sci. China Math. 53 (2010) 137–150. [2] Q. Lin, H. H. Xie, J. C. Xu, 2011, Lower bounds of the discretization for piecewise polynomials, arXiv preprint arXiv:1106.4395. [3] F.S. Luo, Q. Lin, H.H. Xie, Computing the lower and upper bounds of Laplace eigenvalue problem: by combining conforming and nonconforming finite element methods, Sci. China Math. 55 (2012) 1069–1082. [4] H. Jun, Y.Q. Huang, Q. Shen, The lower/upper bound property of approximate eigenvalues by nonconforming finite element methods for elliptic operators, J. Sci. Comput. 58 (2014a) 574–591. [5] H. Jun, Y.Q. Huang, Q. Lin, Lower bounds for eigenvalues of elliptic operators: by nonconforming finite element methods, J. Sci. Comput. 61 (2014b) 196–221. [6] Q. Lin, J.F. Lin, Finite Element Methods: Accuracy and Improvement, Science Press, Beijing, 2006. [7] Q. Lin, H.T. Huang, Z.C. Li, New expansions of numerical eigenvalues for −u = λρ u by nonconforming elements, Math. Comput. 77 (2008) 2061–2084. [8] D.Y. Shi, Y.C. Peng, S.C. Chen, Error estimates for rotated Qrot 1 element approximation of the eigenvalue problem on anisotropic meshes, Appl. Math. Lett. 22 (2009) 952–959. [9] E.A. Dari, R.G. Durán, C. Padra, A posteriori error estimates for nonconforming approximation of eigenvalue problems, Appl. Numer. Math. 62 (2012) 580– 591. [10] Z.M. Chen, S.B. Dai, On the efficiency of adaptive finite element methods for elliptic problems with discontinuous coefficients, SIAM J. Sci. Comput. 24 (2002) 443–462. [11] S. Giani, I.G. Graham, A convergent adaptive method for elliptic eigenvalue problems, SIAM J. Numer. Anal. 47 (2009) 1067–1091.

D. Shi et al. / Applied Mathematics and Computation 273 (2016) 842–855

855

[12] H. Li, Y.D. Yang, The adaptive finite element method based on multi-scale discretizations for eigenvalue problems, Comput. Math. Appl. 65 (2013) 1086–1102. [13] Y.D. Yang, H. Bi, Two-grid finite element discretization schemes based on shifted-inverse power method for elliptic eigenvalue problems, SIAM J. Numer. Anal. 49 (2011) 1602–1624. [14] H.H. Xie, A multigrid method for eigenvalue problem, J. Comput. Phys. 274 (2014) 550–561. [15] Q. Lin, H.H. Xie, A multi-level correction scheme for eigenvalue problems, Math. Comput. 84 (2015) 71–88. [16] B. Mercier, J. Osborn, J. Rappaz, P.A. Raviart, Eigenvalue approximation by mixed and hybrid methods, Math. Comp. 36 (1981) 427–453. [17] F. Gardini, Mixed approximation of eigenvalue problems: a superconvergence result, ESAIM: Math. Model. Numer. Anal. 43 (2009) 853–865. [18] Q. Lin, H.H. Xie, A superconvergence result for mixed finite element approximations of the eigenvalue problem, ESAIM: Math. Model. Numer. Anal. 46 (2012) 797–812. [19] Q. Lin, H.H. Xie, Asymptotic error expansion and Richardson extrapolation of eigenvalue approximations for second order elliptic problems by the mixed finite element method, Appl. Numer. Math. 59 (2009) 1884–1893. [20] R.G. Durán, L. Gastaldi, C. Padra, A posteriori error estimators for mixed approximations of eigenvalue problems, Math. Models Methods Appl. Sci. 9 (1999) 1165–1178. [21] S.H. Jia, H.T. Chen, H.H. Xie, A posteriori error estimator for eigenvalue problems by mixed finite element method, Sci. China Math. 56 (2013) 887–900. [22] H.T. Chen, S.H. Jia, H.H. Xie, Postprocessing and higher order convergence for the mixed finite element approximations of the eigenvalue problem, Appl. Numer. Math. 61 (2011) 615–629. [23] S.C. Chen, H.R. Chen, New mixed element schemes for second order elliptic problem, Math. Numer. Sin. 32 (2010) 213–218. [24] F. Shi, J.P. Yu, K.T. Li, A new mixed finite element scheme for elliptic equations, Chinese J. Eng. Math. 28 (2011) 231–237. [25] D.Y. Shi, Y.D. Zhang, High accuracy analysis of a new nonconforming mixed finite element scheme for Sobolev equations, Appl. Math. Comput. 218 (2011) 3176–3186. [26] Z.F. Weng, X.L. Feng, S.Y. Zhai, Investigations on two kinds of two-grid mixed finite element methods for the elliptic eigenvalue problem, Comput. Math. Appl. 64 (2012) 2635–2646. [27] D.Y. Shi, J.H. Mo, M.H. Li, New mixed finite element scheme for second order eigenvalue problem, Math. Appl. 26 (2013) 506–514. [28] D.Y. Shi, Q.L. Tang, Y.D. Zhang, A new characteristic nonconforming mixed finite element scheme for convection-dominated diffusion problem, J. Appl. Math. (2013) Art. ID 951692. [29] L.F. Pei, D.Y. Shi, Superconvergence of a new nonconforming mixed finite element scheme for elliptic problem, Math. Probl. Eng. Art. ID 829820 (2013). [30] P.G. Ciarlet, The Finite Element Method for Elliptic Problem, North Holland, Amsterdam, 1978. [31] Q. Lin, H.H. Xie, F.S. Luo, Y. Li, Y.D. Yang, Stokes eigenvalue approximations from below with nonconforming mixed finite element methods, Math. Practice Theory 40 (2010) 157–168. [32] D.Y. Shi, C. Xu, J.H. Chen, Anisotropic nonconforming EQrot 1 quadrilateral finite element approximation to second order elliptic problems, J. Sci. Comput. 56 (2013) 637–653. [33] J.Z. Wu, X.Z. Xing, D.Y. Shi, Composite penalty method of a low order anisotropic nonconforming finite element for the Stokes problem, Appl. Math. J. Chinese Univ. Ser. B. 28 (2013) 49–56. [34] D.Y. Shi, C. Xu, EQrot 1 nonconforming finite element approximation to Signorini problem, Sci. China Math. 56 (2013) 1301–1311. [35] F. Chatelin, Spectral Approximation of Linear Operators, Academic Press, New York, 1983. [36] Q. Lin, L. Tobiska, A.H. Zhou, Superconvergence and extrapolation of non-conforming low order finite elements applied to the Poisson equation, IMA J. Numer. Anal. 25 (2005) 160–181. [37] Y.M. Zhao, F.L. Wang, D.Y. Shi, EQrot 1 nonconforming finite element method for nonlinear dual phase lagging heat conduction equations, Acta Math. Appl. Sin. Engl. Ser. 29 (2013) 201–214. [38] D.Y. Shi, Q.L. Tang, Nonconforming H1 -Galerkin mixed finite element method for strongly damped wave equations, Numer. Funct. Anal. Optim. 34 (2013) 1348–1369. [39] R. Rannacher, S. Turek, Simple nonconforming quadrilateral Stokes element, Numer. Methods Partial Differ. Equ. 8 (1992) 97–111. [40] J. Hu, H.Y. Man, Z.C. Shi, Constrained nonconforming Qrot 1 element for stokes flow and planar elasticity, J. Comput. Math. 27 (2005) 311–324. [41] Z.C. Shi, A remark on the optimal order of convergence of Wilson nonconforming element, Math. Numer. Sinica 8 (1986) 159–163. [42] S.C. Chen, D.Y. Shi, Accuracy analysis for quasi-Wilson element, Acta Math. Sci. Ser. B Engl. Ed. 20 (2000) 44–48. [43] D.Y. Shi, D. Zhang, Approximation of nonconforming quasi-Wilson element for sine-Gordon equations, J. Comput. Math. 31 (2013) 271–282. [44] D.Y. Shi, L.F. Pei, Nonconforming quadrilateral finite element method for a class of nonlinear sine-Gordon equations, Appl. Math. Comput. 219 (2013) 9447– 9460. [45] S.H. Jia, H.H. Xie, X.B. Yin, S.Q. Gao, Approximation and eigenvalue extrapolation of biharmonic eigenvalue problem by nonconforming finite element methods, Numer. Methods Partial Differ. Equ. 24 (2008) 435–448.