Robust residual-based a posteriori error estimators for mixed finite element methods for fourth order elliptic singularly perturbed problems

Robust residual-based a posteriori error estimators for mixed finite element methods for fourth order elliptic singularly perturbed problems

Shaohong Du School of Mathematics and Statistics, Chongqing Jiaotong University, Chongqing 400074, China. E-mail: duzheyan.student@sina.com.    Runchang Lin Department of Mathematics and Physics, Texas A&M International University (TAMIU), Laredo, Texas 78041, USA. E-mail: rlin@tamiu.edu. This research is partially supported by the US National Science Foundation under grant DMS-1217268, the National Natural Science Foundation of China under grant 11428103, and a University Research Development Award of TAMIU.    Zhimin Zhang Beijing Computational Science Research Center, Beijing 100193, China. E-mail: zmzhang@csrc.ac.cn. Department of Mathematics, Wayne State University, Detroit, MI 48202, USA. E-mail: ag7761@wayne.edu. This research is partially supported by the US National Science Foundation through grant DMS-1419040 and by the National Natural Science Foundation of China under grants 11471031, 91430216, and U1530401.
Abstract

We consider mixed finite element approximation of a singularly perturbed fourth-order elliptic problem with two different boundary conditions, and present a new measure of the error, whose components are balanced with respect to the perturbation parameter. Robust residual-based a posteriori estimators for the new measure are obtained, which are achieved via a novel analytical technique based on an approximation result. Numerical examples are presented to validate our theory.

Key words. fourth order elliptic singularly perturbed problems, mixed finite element methods, a new measure of the error, robust residual-based a posteriori error estimators

AMS subject classifications. 65N15, 65N30, 65J15

1 Introduction

Let be a bounded polygonal or polyhedral domain with Lipschitz boundary in . Consider the following fourth-order singularly perturbed elliptic equation

 ε2△2u−△u=f   in  Ω (1)

with boundary conditions

 u=△u=0   on  Γ (2)

or

 u=∂u∂n=0   on  Γ, (3)

where , is the standard Laplace operator, and denotes the outer normal derivative on . In two dimensional cases, the boundary value problems (1)-(2) and (1)-(3) arise in the context of linear elasticity of thin bucking plate with representing the displacement of the plate. The dimensionless positive parameter , assumed to be small (i.e., ), is defined by

 ε=t3E12(1−ν2)l2T,

where, is the thickness of the plate, is the Young modulus of the elastic material, is the Poisson ratio, is the characteristic diameter of the plate, and is the absolute value of the density of the isotropic stretching force applied at the end of the plate [20]. In three dimensions, problems (1)-(2) and (1)-(3) can be a gross simplification of the stationary Cahn-Hilliard equations with being the length of the transition region of phase separation.

Conforming, nonconforming, and mixed finite element methods for fourth order problem have been extensively studied [2, 3, 7, 8, 15, 19, 28, 29, 31, 32, 34, 35, 36]. However, its a posteriori error estimation is a much less explored topic. Even for the Kirchhoff plate bending problem, the finite element a posteriori error analysis is still in its infancy. In 2007, Beiro et al. [5] developed an estimator for the Morley element approximation using the standard technique for nonconforming element. Later, Hu et al. [25] improved the methods of [5, 37, 38] by dropping two edge jump terms in both the energy norm of the error and the estimator, and by dropping the normal component in the estimators of [5, 37]. Therefore, a naive extension of the estimators in [5, 37, 38] to the current problem may probably not be robust in the parameter .

Designing robust a posteriori estimators is challenging, especially for singularly perturbed problems, since constants occurring in estimators usually depend on the small perturbation parameter . This motivates us to think about the question: What method and norm are suitable for the singularly perturbed fourth-order elliptic problem? In the literature, is a widely used measure for the primal weak formulation. We recall a priori estimates in [23] for boundary condition (3) and convex domain :

 ∥u∥Hs(Ω)≤Cε32−s∥f∥L2(Ω),   for s=2,3. (4)

Hereafter, we use to denote a generic constant independent of with different value at different occurrence. This leads to

 ε∥△u∥L2(Ω)≤Cε∥u∥H2(Ω)≤Cε1/2∥f∥L2(Ω).

Multiply both sides of (1) by , and then integrate over . Using integration by parts and boundary condition (2) or (3), we have from the Poincaré inequality that

 ε2(△u,△u)+(∇u,∇u)=(f,u)≤∥f∥L2(Ω)∥u∥L2(Ω)≤C∥f∥L2(Ω)∥∇u∥L2(Ω).

As a consequence, . This suggests that the two components of are unbalanced with respect to if .

Furthermore, if we set , then problem (1) is written as

 −ε2△ψ+ψ=f.

Note that has boundary layer, but usually does not have one. Thus,

 (ε2∥△u∥2L2(Ω)+∥∇u∥2L2(Ω))1/2=(ε2∥ψ∥2L2(Ω)+|u|2H1(Ω))1/2

approaches to as , which fails to describe the layer of .

An observation of the two decoupled equations and suggests that the two measures and can portray the layer of and the first and second derivatives of . From (4), we have

 ε∥∇ψ∥L2(Ω)≤Cε∥u∥H3(Ω)≤Cε−1/2∥f∥L2(Ω).

Notice that

 ∥ψ∥L2(Ω)≤C∥u∥H2(Ω)≤Cε−1/2∥f∥L2(Ω).

If , then and are balanced with respect to for the boundary condition (3). These inspire us to think about the mixed finite element method for the problem (1) and the two aforementioned measures.

However, the mixed finite element method for the problem (1) is a much less explored topic, since there exist some special problems such as the fourth order problem, where attempts at using the results of Brezzi and Babuška were not entirely successful since not all of the stability conditions were satisfied, cf. [3] and the reference therein. To overcome this difficulty, Falk et al. developed abstract results from which optimal error estimates for these (biharmonic equation) and other problems could be derived ([19, 8]). However, it is not easy to extend the results of [19] to the problem (1), because of the existence of an extra term and the singular perturbation parameter . Recently, for a fourth order reaction diffusion equation, the error estimates of its mixed finite element method was derived in [15]. We refer to [10, 22] about the a posteriori estimation of Ciarlet-Raviart methods for the biharmonic equation.

In this work, our goal is to develop robust residual-type a posteriori estimators for a mixed finite element method for the problem (1) in the two aforementioned measures. The main difficulty lies in the fact that the boundary condition (3) does not include any information on the immediate variable . In order to overcome this difficulty, we develop a novel technique to analyze residual-based a posteriori error estimator. The key idea is to replace a function (such that ) without boundary restriction by a function with boundary restriction, which catches at least “ times” of in the -weighted energy norm (see Lemma 3.3 below). Combining this novel design with standard tools, we develop uniformly robust residual-type a posteriori estimators with respect to the singularly perturbed parameter in the two aforementioned measures. We refer to the reference [26] on balanced norm for mixed formulation for singularly perturbed reaction-diffusion problems.

The rest of this paper is organized as follows: In Section 2, we introduce mixed weak formulations and some notations, and prove an equivalent relation between the primal weak solution and the weak solution determined by its mixed formulation. Some preliminary results are provided in Section 3. Residual-type a posteriori estimators are developed and proven to be reliable in Section 4. An efficient lower bound is proved in Section 5. In Section 6, numerical tests are provided to support our theory.

2 The mixed weak formulations

Setting , and employing the boundary condition (2), we attain the Ciarlet-Raviart mixed problem :

 ⎧⎪⎨⎪⎩−ε2△ψ+ψ=fin  Ω  −△u=ψin Ω u=ψ=0on Γ. (5)

Similarly, using the boundary condition (3), we arrive at the Ciarlet-Raviart mixed formulation :

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩−ε2△ψ+ψ=fin  Ω  −△u=ψin Ω u=∂u∂n=0% on Γ. (6)

For any bounded open subset of with Lipschitz boundary , let and be the standard Lebesgue and Sobolev spaces equipped with standard norms and , (see [1] for details). Note that . We denote the semi-norm in . Similarly, denote and the inner products on and , respectively. We shall omit the symbol in the notations above if .

The weak formulation of problem reads: find such that

 {(ε2∇ψ,∇φ)+(ψ,φ)=(f,φ)∀ φ∈H10(Ω),(∇u,∇v)=(ψ,v)∀ v∈H10(Ω). (7)

The weak formulation of the problem reads: find such that

 {(ε2∇ψ,∇φ)+(ψ,φ)=(f,φ)∀ φ∈H10(Ω),(∇u,∇v)=(ψ,v)∀ v∈H1(Ω). (8)

Note that, by the Lax-Milgram lemma, both systems (7) and (8) have a unique solution. In fact, by regularity theory for elliptic problems [21], if is convex and , then and . Thus (7) has solution, which is unique since its homogeneous system has only one solution satisfying . Similar conclusion can be drawn for the system (8).

It is well known that the primal weak formulation of (1)-(2) is: find such that

 (ε2△~u,△v)+(∇~u,∇v)=(f,v),   ∀v∈H2(Ω)∩H10(Ω), (9)

and that the one of (1)-(3) is: find such that

 (ε2△~u,△v)+(∇~u,∇v)=(f,v),   ∀v∈H20(Ω). (10)

The classical results of PDEs imply that (9) and (10) have unique solutions (see [21]). A natural question is whether the determined by (7) (or (8)) is the solution of (9) (or (10)). In [39], for biharmonic equation on a reentrant corners polygon, a counterexample is shown. The following theorems answer this question.

{theorem}

The solution of (9) and the determined by (7) are identical if and only if . {proof} The necessity is trivial. If the solution of (7) is such that , then we have from the second equation

 (−△u,w)=(ψ,w),   ∀w∈H10(Ω).

Notice that is dense in . It follows that

 (−△u,w)=(ψ,w),   ∀w∈L2(Ω).

Integration by parts yields

 (−△u,△φ)=(ψ,△φ)=−(∇ψ,∇φ),   ∀φ∈H2(Ω)∩H10(Ω),

which implies

 (ε2△u,△φ)=(ε2∇ψ,∇φ),   ∀φ∈H2(Ω)∩H10(Ω). (11)

We obtain from (11) and the second equation of (7) that

 (ε2△u,△φ)+(∇u,∇φ)=(ε2∇ψ,∇φ)+(ψ,φ),   ∀φ∈H2(Ω)∩H10(Ω).

In terms of (9), we proved that is the solution of (9).

{theorem}

The solution of (10) and the determined by (8) are identical if and only if . {proof} The necessity is trivial. From the second equation of (8), integration by parts, and variational principle, we know that the Neumann boundary condition on is automatically satisfied. Following the proof of Theorem 2, we know that if the solution of (8) is in , then is the solution of (10).

Let be a shape regular partition of into triangles (tetrahedra for ) or parallelograms (parallelepiped for ) satisfying the angle condition [13], i.e., there exists a constant such that

 (12)

where . Let be the space of polynomials of total degree at most if is a simplex, or the space of polynomials with degree at most for each variable if is a parallelogram/parallelepiped. Define the finite element spaces and by

 Vh:={vh∈C(¯¯¯¯Ω):vh|K∈Pk(K), ∀K∈Th}

and

 V0h:={vh∈Vh:vh|Γ=0},

respectively.

We introduce the mixed finite element method for problem : find such that

 (13)

For problem , the mixed problem reads: find such that

 (14)

By standard arguments, problem (13) possesses a unique solution provided there exist functions and satisfying

 (15)

then is the trivial solution to the system. In fact, taking in the first equation of (15), one gets . Setting in the second equation of (15), one obtains . Similarly, it is verified that problem (14) has also a unique solution.

We define a measure of the error between the exact solution and the numerical solution by

 ∥(ψ−ψh,u−uh)∥2:=∥ψ−ψh∥2E+|u−uh|21,

where

 ∥ψ−ψh∥E=(ε2|ψ−ψh|21+∥ψ−ψh∥2)1/2

is the standard energy norm of the numerical error . In this paper, we aim at robust a posterior error estimators for the numerical errors , , and .

We next introduce some notations that will be used later. We denote the set of interior sides (if ) or faces (if ) in , the set of sides or faces of , and the union of all elements in sharing at least one point with . For a side or face in , which is the set of element sides or faces in , let be the diameter of , and be the union of all elements in sharing . For a function in the “broken Sobolev space” , we define as the jump of across an interior side or face , where and are the two neighboring elements such that .

Throughout of this paper, we denote by a constant depending only on , and denote by constants depending on the mesh shape regularity and . In what follows we use the notation to represent with a generic constant independent of mesh size. In addition, abbreviates .

3 Preliminary results

For problem , and are decoupled. However, does not obtain any information directly from boundary conditions. It will be difficult to develop residual-based a posteriori error estimates if the residual on the boundary is not clear. To overcome this difficulty, we shall develop a novel analytical technique (see Section 4), which is based on the following approximation result. {lemma} Let which satisfies (operator to be understood in weak sense). Then it holds

 infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2=1.
{proof}

Consider the functional

 J(w)=ε2|v−w|21+||v−w||2,  ∀ w∈H10(Ω).

Minimization of such functional in immediately leads to the following variational problem: Find such that

 ∫Ω(ε2∇(v−~v)⋅∇ϕ+(v−~v)ϕ)dx=0,  ∀ ϕ∈H10(Ω). (16)

Integrating by parts, we arrive at

 ∫Ω(−ε2△(v−~v)+v−~v)ϕdx=0,  ∀ ϕ∈H10(Ω),

which implies

 −ε2△~v+~v=−ε2△v+v=0.

So is the solution to the following problem:

 {−ε2△~v+~v=0%in Ω,~v=0on Γ. (17)

Since the problem (17) has only trivial solution, we have

 infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2=infw∈H10(ω)J(w)ε2|v|21+||v||2 = ε2|v−~v|21+||v−~v||2ε2|v|21+||v||2=ε2|v|21+||v||2ε2|v|21+||v||2=1,

which completes the proof.

{lemma}

If there holds the following relation

 infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2=1,

then satisfies . {proof} Since , the condition

 1=infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2

can be satisfied when . On the other hand, is the solution to the variational problem: Find such that

 ∫Ω(ε2∇(v−~v)⋅∇ϕ+(v−~v)ϕ)dx=0,  ∀ ϕ∈H10(Ω).

By integrating by parts, we obtain

 ∫Ω(−ε2△(v−~v)+(v−~v))ϕdx=0,  ∀ ϕ∈H10(Ω),

which leads to , this means

 −ε2△v+v=−ε2△~v+~v=0.
{lemma}

Let such that . Then there exists such that

 infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2≤γ.
{proof}

We only prove the case using proof by contradiction, since is obvious. Assume that there does not exist such that

 infw∈H10(ω)(ε2|v−w|21+||v−w||2)≤γ(ε2|v|21+||v||2),

which means

 infw∈H10(ω)(ε2|v−w|21+||v−w||2)>γ(ε2|v|21+||v||2)

for all . From the proof of Lemma 3, there exists such that

 ε2|v|21+||v||2≥infw∈H10(Ω)(ε2|v−w|21+||v−w||2)=ε2|v−~v|21+||v−~v||2>γ(ε2|v|21+||v||2).

In particular, for any , let , it holds that

 ε2|v|21+||v||2≥ε2|v−~v|21+||v−~v||2>n−1n(ε2|v|21+||v||2).

Since is the upper bound for with respect to the positive integer number , and is its supremum. Therefore,

 ε2|v|21+||v||2≤ε2|v−~v|21+||v−~v||2≤ε2|v|21+||v||2,

which means

 ε2|v−~v|21+||v−~v||2=ε2|v|21+||v||2.

 infw∈H10(Ω)ε2|v−w|21+||v−w||2ε2|v|21+||v||2=ε2|v−~v|21+||v−~v||2ε2|v|21+||v||2=1.

From Lemma 3, we have , this leads to a contradiction. We complete the proof.

Remark 3.1

Note that we can not prove that is dense in by recursion by using Lemma 3, because of .

Denote by the quasi-interpolation operator of Clément (cf. [13, 14, 33]). {lemma} For all , define and the weighted factors by

 αT:=min{hTε−1,1}  and  αE:=ε−1/2min{hTε−1,1},

respectively. Then the following local error estimates hold for :

 ∥v−Ihv∥T≲αT∥v∥E,~ωT (18)

and

 ∥v−Ihv∥E≲αE∥v∥E,~ωT. (19)
{proof}

Following the line of the proof of Lemma 3.2 in [33], we obtain the desired estimates (18) and (19).

For , denote and the two bubble functions defined in [33], and a continuation operator introduced in [33] by

 PE:L∞(E)→L∞(ωE),

which maps polynomials onto piecewise polynomials of the same degree. {lemma} The following estimates hold for all (the set of polynomials of degree at most ) and

 ∥v∥2T≲(v,ψTv)T, (20)
 ∥vψT∥T≤∥v∥T, (21)
 ∥vψT∥E,T≲α−1T∥v∥T. (22)

Furthermore, for , set . Then there hold the following estimates for all and .

 ∥σ∥2E≲(σ,ψE,θEPEσ)E, (23)
 ∥ψE,θEPEσ∥ωE≲ε1/2min{hEε−1,1}1/2∥σ∥E, (24)
 ∥ψE,θEPEσ∥E,ωE≲ε1/2min{hEε−1,1}−1/2∥σ∥E. (25)
{proof}

Following the line of the proof of Lemma 3.3 in [33], we attain (20)-(25).

4 A reliable upper bound

For all , define and the elementwise indicators of and , respectively, by

 ηψ,T:={α2T∥f+ε2△ψh−ψh∥2T+12∑E∈ET∩E0hα2E∥∥[ε2∂ψh∂n]∥∥2E}1/2

and

 ηu,T:={h2T∥△uh+ψh∥2T+12∑E∈ET∩E0hhE∥∥[∂uh∂n]∥∥2E}1/2.
{theorem}

Let and be the solutions to (7) and (13), respectively. Then there exist positive constants , and , independent of the mesh-size function and , such that

 ∥ψ−ψh∥E≤C1{∑T∈Thη2ψ,T}1/2, (26)
 |u−uh|1≤C2{∑T∈Thη2ψ,T+η2u,T}1/2, (27)
 ∥(ψ−ψh,u−uh)∥≤C3{∑T∈Thη2ψ,T+η2u,T}1/2. (28)
{proof}

From the definition of the measure , (28) follows from (26) and (27). We need to prove (26) and (27). We have from the first equations of (7) and (13) that

 (ε2∇(ψ−ψh),∇φh)+(ψ−ψh,φh)=0,   ∀ φh∈V0h. (29)

For any , let be the Clemént interpolation of in , i.e., . Applying integration by parts and (29), we get

 (ε2∇(ψ−ψh),∇φ)+(ψ−ψh,φ) = (ε2∇(ψ−ψh),∇(φ−φh))+(ψ−ψh,φ−φh) = ∑T∈Th∫T(−ε2△ψ+ε2△ψh+ψ−ψh)(φ−φh)+∫∂Tε2∂(ψ−ψh)∂n(φ−φh) = ∑T∈Th∫T(f+ε2△ψh−ψh)(φ−φh)−∫∂Tε2∂ψh∂n(φ−φh) ≤ ∑T∈Th{∥f+ε2△ψh−ψh∥T∥φ−φh∥T+12∑E∈ET∩E0h∥∥[ε2∂ψh∂n]∥∥E∥φ−φh∥E}. (30)

Notice that

 ∥ψ−ψh∥E=(ε2∇(ψ−ψh),∇(ψ−ψh))+(ψ−ψh,ψ−ψh)∥ψ−ψh∥E≤sup0≠φ∈H10(Ω)(ε2∇(ψ−ψh),∇φ)+(ψ−ψh,φ)∥φ∥E. (31)

The first estimate (26) follows from a combination of (31), (30), and (18)-(19).

We next prove (27). From the second equation of (7) and (13), we get

 (∇(u−uh),∇vh)=(ψ−ψh,vh),   ∀vh∈V0h. (32)

Similarly, we have, for any and ,

 (∇(u−uh),∇v)=(∇(u−uh),∇(v−vh))+(∇(u−uh),∇vh) = ∑T∈Th∫T(−△u+△uh)(v−vh)+∫∂T∂(u−uh)∂n(v−vh)+(ψ−ψh,vh) = ∑T∈Th∫T(−△u+△uh−ψ+ψh)(v−vh)+∫∂T∂(u−uh)∂n(v−vh)+(ψ−ψh,v) = ∑T∈Th∫T(△uh+ψh)(v−vh)−∫∂T∂uh∂n(v−vh)+(ψ−ψh,v) ≤ ∑T∈Th{∥△u