Superconvergence of a New Nonconforming Mixed Finite Element Scheme for Elliptic Problem

Abstract

A new nonconforming mixed finite element scheme for the second-order elliptic problem is proposed based on a new mixed variational form. It has the lowest degrees of freedom on rectangular meshes. The superclose property is proven by employing integral identity technique. Then global superconvergence result is derived through interpolation postprocessing operators. At last, some numerical experiments are carried out to verify the theoretical analysis.

1. Introduction

Mixed finite element method (MFEM) is an important branch of FEMs and has been used widely in numerical computation of practical problems. A lot of studies on this aspect have been devoted to the second-order elliptic problems [1–5], in which two approximation spaces of MFEM should satisfy the famous B-B condition [6]. However, since the variable 𝑢 and flux 𝜓=∇𝑢 belong to 𝐿²(Ω) and 𝐻(div,Ω), respectively, it is not easy to construct a stable MFE space pair. In order to circumvent or ameliorate this deficiency, many approaches have been proposed, such as the least squares FEM [7], stabilization FEM [8], and 𝐻¹-Galerkin FEM [9]. Recently, [10] presented a new MFEM, in which B-B condition is automatically satisfied when 𝑀_ℎ and 𝐻_ℎ meet a relation of inclusion; that is, ∇𝑀_ℎ⊂𝐻_ℎ, where 𝑀_ℎ and 𝐻_ℎ are finite element approximation spaces of 𝑢 and flux 𝜓, respectively. This advantage makes the construction of stable MFE space pair extremely simple and convenient. A family of triangular and rectangular conforming MFE space pairs with lower degrees of freedom is constructed in [10], in which the total degrees of freedom of first-order and second-order MFE schemes are about 5𝑁_𝑝 and 16𝑁_𝑝, respectively; herein 𝑁_𝑝 denotes the number of nodal points in subdivision. Later, [11] derived the similar results with [10] for conforming MFEM and gave a numerical example.

In this paper, motivated by the idea of [10, 11], we first construct a new nonconforming MFEM (NMFEM). The original variable 𝑢 is approximated by the constrained 𝑄₁^rot element space [12] and the flux 𝜓 by piecewise constant vectors space, respectively. Note that the total degrees of freedom of the NMFEM are only about 3𝑁_𝑝, the lowest on rectangular meshes. We prove that it satisfies B-B condition. Then by the use of integral identity technique, we derive the superclose property for 𝑢 in energy norm and flux 𝜓 in 𝐿² norm. Furthermore, the global superconvergence result with order O(ℎ²) is obtained through interpolation postprocessing operators. Finally, some numerical results are provided to verify the theoretical analysis. It is observed that, compared with FEM using 4-node quadrilateral (FEM-Q4) and MFEM (𝑢 and flux 𝜓 are approximated by piecewise constants and the Raviart-Thomas element, resp.), NMFEM behaves well and has higher rates than MFEM, and it has almost the same rate as FEM-Q4 for 𝑢 in 𝐿² norm, but a higher rate than FEM-Q4 for 𝑢 in energy norm. Moreover, NMFEM is effective and accurate for the diffusion problem studied in [13].

The remainder of this paper is organized as follows. In Section 2, we introduce NMFEM and derive the superclose property and superconvergence results. In Section 3, we carry out some numerical experiments to show the performance of NMFEM.

We will use standard notations for the Sobolev spaces 𝐻^𝑚(Ω) with norm ‖⋅‖_𝑚 and seminorm |⋅|_𝑚, 𝐻^𝑚(𝐾) with norm ‖⋅‖_𝑚,𝐾 and seminorm |⋅|_𝑚,𝐾, where 𝑚>0 is an integer. Besides, let ‖⋅‖₀ and ‖⋅‖_0,𝐾 be the 𝐿²(Ω) norm and 𝐿²(𝐾) norm, respectively. Throughout the paper, 𝐶 denotes a positive constant independent of the mesh parameter ℎ and may be different at each appearance.

2. Superconvergence Analysis for NMFEM

Consider the following elliptic problem: −Δ𝑢=𝑓,inΩ,𝑢=0,on𝜕Ω, (1) where Ω⊂𝑅² is a bounded convex polygon domain, 𝑓∈𝐿²(Ω).

Let 𝜓=∇𝑢, and then problem (1) is equivalent to the following equations: 𝜓−∇𝑢=0,inΩ,−div𝜓=𝑓,inΩ,𝑢=0,on𝜕Ω. (2) We adopt a new mixed variational form in [10] of problem (2). Find (𝜓,𝑢)∈𝐻×𝑀 such that 𝑎𝑏(𝑣)𝜓,𝜑+𝑏𝜑,𝑢=0,∀𝜑∈𝐻,𝜓,𝑣=𝐺,∀𝑣∈𝑀, (3) where 𝐻=(𝐿²(Ω))², 𝑀=𝐻¹₀(Ω), 𝑎(𝜓,𝜑)=∫_Ω𝜓⋅𝜑𝑑𝑥𝑑𝑦, 𝑏(𝜑,𝑣)=−∫_Ω𝜑⋅∇𝑣𝑑𝑥𝑑𝑦, 𝐺(𝑣)=−∫_Ω𝑓𝑣𝑑𝑥𝑑𝑦.

Obviously, 𝑎(⋅,⋅) and 𝑏(⋅,⋅) are continuous bilinear functionals, 𝐺(⋅) is a continuous linear functional, and for all 𝜑∈𝐻,𝑎(𝜑,𝜑)=∫_Ω𝜑⋅𝜑𝑑𝑥𝑑𝑦=‖𝜑‖²₀. Moreover, for 𝑣∈𝑀, we have ∇𝑣∈𝐻. So there exists a constant 𝛽>0 such that sup_𝜑∈𝐻𝑏𝜑,𝑣𝜑₀≥𝑏(−∇𝑣,𝑣)‖−∇𝑣‖₀=|𝑣|₁≥𝛽‖𝑣‖₁, (4) that is, the B-B condition is satisfied, and therefore (3) has a unique solution (𝜓,𝑢).

Let 𝑄((𝜓,𝑢),(𝜑,𝑣))=𝑎(𝜓,𝜑)+𝑏(𝜑,𝑢)−𝑏(𝜓,𝑣), and then (3) can be written as 𝑄𝜓,𝑢,𝜑,𝑣=−𝐺(𝑣),∀𝜑∈𝐻,𝑣∈𝑀. (5)

Let 𝑇_ℎ be a rectangular partition of the domain Ω. For a given element 𝐾∈𝑇_ℎ, its four vertices are denoted by 𝑎₁=(𝑥_𝐾−ℎ_𝑥,𝑦_𝐾−ℎ_𝑦), 𝑎₂=(𝑥_𝐾+ℎ_𝑥,𝑦_𝐾−ℎ_𝑦), 𝑎₃=(𝑥_𝐾+ℎ_𝑥,𝑦_𝐾+ℎ_𝑦), 𝑎₄=(𝑥_𝐾−ℎ_𝑥,𝑦_𝐾+ℎ_𝑦), and four edges by 𝑙_𝑖=𝑎_𝑖𝑎_𝑖+1(𝑖=1,2,3,4mod4), ℎ_𝐾=max_{𝐾∈𝑇ℎ}{ℎ_𝑥,ℎ_𝑦},ℎ=max_{𝐾∈𝑇ℎ}ℎ_𝐾. Let 𝐾=[−1,1]×[−1,1] be the reference element with nodes 𝑎₁=(−1,−1), 𝑎₂=(1,−1), 𝑎₃=(1,1), 𝑎₄=(−1,1) and edges 𝑙_𝑖̂=𝑎_𝑖𝑎_𝑖+1(𝑖=1,2,3,4mod4).

Define the affine mapping 𝐹∶𝐾→𝐾 by 𝑥=𝑥_𝐾+ℎ_𝑥𝜉,𝑦=𝑥_𝐾+ℎ_𝑦𝜂. (6) Let 𝑃₁(𝐾) be the space of polynomials with degrees ≤1 defined on 𝐾, and then the constrained 𝑄₁^rot element space CR^ℎ₀ is defined by [12]: CR^ℎ₀=𝑣;𝑣|_𝐾=𝑣̂∘𝐹⁻¹,𝑣̂∈𝑃₁𝐾,_𝑙[𝑣]𝑑𝑠=0,𝑙⊂𝜕𝐾,∀𝐾∈𝑇_ℎ, (7) where [𝑣] denotes the jump value of 𝑣 across the boundary 𝑙, and [𝑣]=𝑣 if 𝑙⊂𝜕Ω.

Let 𝑁^𝑉_𝑖 denote the number of interior nodes. It has been proven in [12] that dim(CR^ℎ₀)=𝑁^𝑉_𝑖 and 𝑃₁(𝐾)=span{𝜙_𝑖,𝑖=1∼4}, where 𝜙_𝑖 are defined associated with nodes 𝑎_𝑖(𝑖=1,2,3,4) of 𝐾 as 𝜙₁=141−𝜉−𝜂,𝜙₂=14,𝜙1+𝜉−𝜂₃=141+𝜉+𝜂,𝜙₄=141−𝜉+𝜂. (8)

We choose the following FE spaces 𝑀_ℎ and 𝐻_ℎ to approximate 𝑀 and 𝐻, respectively: 𝑀_ℎ=CR^ℎ₀,𝐻_ℎ=𝑞=𝑞₁,𝑞₂;𝑞_𝐾∈𝑄₀(𝐾)²,∀𝐾∈𝑇_ℎ, (9) where 𝑄₀(𝐾) is the space of constants on 𝐾. Obviously, the total degree of freedoms of the nonconforming MFE space pair 𝐻_ℎ×𝑀_ℎ is 3𝑁_𝑝, and ‖⋅‖_ℎ=(∑_{𝐾∈𝑇ℎ}|⋅|²_1,𝐾)^1/2 is a norm over 𝑀_ℎ.

Then the MFE approximation of problem (3) is to find (𝜓_ℎ,𝑢_ℎ)∈𝐻_ℎ×𝑀_ℎ such that 𝑎𝜓_ℎ,𝜑_ℎ+𝑏_ℎ𝜑_ℎ,𝑢_ℎ=0,∀𝜑_ℎ∈𝐻_ℎ,𝑏_ℎ𝜓_ℎ,𝑣_ℎ=𝐺𝑣_ℎ,∀𝑣_ℎ∈𝑀_ℎ, (10) where 𝑏_ℎ(𝜑_ℎ,𝑣_ℎ)=−∑_{𝐾∈𝑇ℎ}∫_𝐾𝜑_ℎ⋅∇𝑣_ℎ𝑑𝑥𝑑𝑦.

For 𝑣_ℎ∈𝑀_ℎ, from the affine mapping 𝐹 and the definition of 𝑀_ℎ, we can get 𝑣_ℎ|_𝐾∈𝑃₁(𝐾), and then ∇𝑣_ℎ|_𝐾∈(𝑄₀(𝐾))², that is, ∇𝑣_ℎ∈𝐻_ℎ. Thus sup_{𝜑ℎ∈𝐻ℎ}𝑏_ℎ𝜑_ℎ,𝑣_ℎ𝜑_ℎ0≥𝑏_ℎ−∇𝑣_ℎ,𝑣_ℎ−∇𝑣_ℎ0=∑_{𝐾∈𝑇ℎ}∫_𝐾∇𝑣_ℎ⋅∇𝑣_ℎ𝑑𝑥𝑑𝑦−∇𝑣_ℎ0=𝑣_ℎ²_ℎ𝑣_ℎℎ=𝑣_ℎℎ, (11) that is, the discrete B-B condition holds, and (10) has a unique solution (𝜓_ℎ,𝑢_ℎ).

For all (𝑞_ℎ,𝜆_ℎ)∈𝐻_ℎ×𝑀_ℎ, let 𝑄_ℎ((𝑞_ℎ,𝜆_ℎ),(𝜑_ℎ,𝑣_ℎ))=𝑎(𝑞_ℎ,𝜑_ℎ)+𝑏_ℎ(𝜑_ℎ,𝜆_ℎ)−𝑏_ℎ(𝑞_ℎ,𝑣_ℎ), and then we can prove the following two important lemmas.

Lemma 1. For all (𝑞_ℎ,𝜆_ℎ)∈𝐻_ℎ×𝑀_ℎ, the following inequality holds: sup_{(𝜑ℎ,𝑣ℎ)∈𝐻ℎ×𝑀ℎ}𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝜑_ℎ,𝑣_ℎ𝜑_ℎ0+𝑣_ℎℎ≥𝐶𝑞_ℎ0+𝜆_ℎℎ. (12)

Proof. In order to prove (12), we will construct a pair (𝜑_ℎ,𝑣_ℎ)∈𝐻_ℎ×𝑀_ℎ satisfying 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝜑_ℎ,𝑣_ℎ≥𝐶𝑞_ℎ0+𝜆_ℎℎ𝜑_ℎ0+𝑣_ℎℎ. (13) Obviously, 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝑞_ℎ,𝜆_ℎ=𝑞_ℎ²₀. (14) On the other hand, for a given arbitrary but fixed 𝜆_ℎ∈𝑀_ℎ, by the discrete B-B condition, there exists a 𝑤_ℎ∈𝐻_ℎ such that _{𝐾∈𝑇ℎ𝐾}𝑤_ℎ⋅∇𝜆_ℎ𝑑𝑥𝑑𝑦≥𝑤_ℎ0𝜆_ℎℎ,𝑤_ℎ0=𝜆_ℎℎ. (15) So, for 0<𝑟<1, there holds 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,−𝑟𝑤_ℎ𝑞,0=−𝑟𝑎_ℎ,𝑤_ℎ+𝑟_{𝐾∈𝑇ℎ𝐾}𝑤_ℎ⋅∇𝜆_ℎ𝑞𝑑𝑥𝑑𝑦≥−𝑟_ℎ0𝜆_ℎℎ+𝑟𝜆_ℎ²_ℎ≥−𝑟𝑞_ℎ²₀+34𝑟𝜆_ℎ²_ℎ. (16) As a result, setting (𝜑_ℎ,𝑣_ℎ)=(𝑞_ℎ−𝑟𝑤_ℎ,𝜆_ℎ), we have 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝑞_ℎ−𝑟𝑤_ℎ,𝜆_ℎ=𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝑞_ℎ,𝜆_ℎ+𝑄_ℎ𝑞_ℎ,𝜆_ℎ,−𝑟𝑤_ℎ≥(𝑞,01−𝑟)_ℎ²₀+34𝑟𝜆_ℎ²_ℎ, (17) which implies that 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝑞_ℎ−𝑟𝑤_ℎ,𝜆_ℎ≥𝐶₁𝑞_ℎ²₀+𝜆_ℎ²_ℎ≥𝐶₁2𝑞_ℎ0+𝜆_ℎℎ², (18) where 𝐶₁=min{1−𝑟,(3/4)𝑟}.
Note that 𝑞_ℎ−𝑟𝑤_ℎ0+𝜆_ℎℎ≤𝑞_ℎ0+𝑟𝑤_ℎ0+𝜆_ℎℎ=𝑞_ℎ0+(1+𝑟)𝜆_ℎℎ≤(1+𝑟)𝑞_ℎ0+𝜆_ℎℎ, (19) and we have 𝑄_ℎ𝑞_ℎ,𝜆_ℎ,𝑞_ℎ−𝑟𝑤_ℎ,𝜆_ℎ≥𝐶₁2()𝑞1+𝑟^ℎ0+𝜆_ℎℎ𝑞_ℎ−𝑟𝑤_ℎ0+𝜆_ℎℎ, (20) which follows the desired result (12).
Let Π_ℎ and 𝐼_ℎ denote the associated interpolation operators of 𝑀_ℎ and conforming bilinear element space, respectively. Let 𝜋_ℎ be the 𝑄₁^rot element interpolation operator [14]; that is, for all 𝑣∈𝐻¹(𝐾), 𝜋_ℎ𝑣∈𝑄₁^rot satisfying ∫_𝑙𝜋_ℎ𝑣𝑑𝑠=∫_𝑙𝑣𝑑𝑠, where 𝑄₁^rot={𝑣_ℎ;𝑣_ℎ|_𝐾∈span{1,𝑥,𝑦,𝑥²−𝑦²},∫_𝑙[𝑣_ℎ]𝑑𝑠=0,𝑙⊂𝜕𝐾}. It has been shown in [12] that, for all 𝑣∈𝐻²(𝐾), Π_ℎ𝑣=𝜋_ℎ𝐼_ℎ𝑣∈𝑀_ℎ because 𝑥𝑦 is a bubble function for 𝜋_ℎ.

Lemma 2. Assume that 𝑢∈𝐻³(Ω), 𝜓∈(𝐻²(Ω))², we have 𝑏_ℎ𝜓−𝐽_ℎ𝜓,𝑣_ℎ=0,∀𝑣_ℎ∈𝑀_ℎ,_{𝐾∈𝑇ℎ𝐾}∇𝑢−Π_ℎ𝑢⋅𝜑_ℎ𝑑𝑥𝑑𝑦≤𝐶ℎ²‖𝑢‖₃𝜑_ℎ0,∀𝜑_ℎ∈𝐻_ℎ,_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ𝑑𝑠≤𝐶ℎ²𝜓₂𝑣_ℎℎ,∀𝑣_ℎ∈𝑀_ℎ, (21)(22)(23) where 𝐽_ℎ∶(𝐿²(Ω))²→𝐻_ℎ is a interpolation operator satisfying ∫_𝐾(𝜓−𝐽_ℎ𝜓)𝑑𝑥𝑑𝑦=0.

Proof. Since (23) has been proven by one of the authors in [15], we only need to prove (21) and (22).
In fact, because 𝜕𝑣_ℎ/𝜕𝑥 and 𝜕𝑣_ℎ/𝜕𝑦 are constants on each 𝐾∈𝑇_ℎ, we have 𝑏_ℎ𝜓−𝐽_ℎ𝜓,𝑣_ℎ=−_{𝐾∈𝑇ℎ𝐾}𝜓−𝐽_ℎ𝜓⋅∇𝑣_ℎ𝑑𝑥𝑑𝑦=0, (24) which is (21).
On the other hand, note that _𝐾∇𝑢−Π_ℎ𝑢⋅𝜑_ℎ=𝑑𝑥𝑑𝑦_𝐾∇𝑢−𝐼_ℎ𝑢⋅𝜑_ℎ𝑑𝑥𝑑𝑦+_𝐾∇𝐼_ℎ𝑢−Π_ℎ𝑢⋅𝜑_ℎ𝑑𝑥𝑑𝑦. (25) Since 𝜑_ℎ=(𝜑_1ℎ,𝜑_2ℎ) is a constant vector on 𝐾 and Π_ℎ𝑣=𝜋_ℎ𝐼_ℎ𝑣, it follows from integration by parts that _𝐾∇𝐼_ℎ𝑢−Π_ℎ𝑢⋅𝜑_ℎ𝑑𝑥𝑑𝑦=𝜑_{ℎ𝐾𝜕𝐾}𝐼_ℎ𝑢−𝜋_ℎ𝐼_ℎ𝑢𝑑𝑠=𝜑_{ℎ𝐾𝜕𝐾}𝐼_ℎ𝑢−𝐼_ℎ𝑢𝑑𝑠=0. (26) Furthermore, let 𝐹(𝑦)=(1/2)((𝑦−𝑦_𝐾)²−ℎ²_𝑦). Note that (𝑢−𝐼_ℎ𝑢)(𝑥_𝐾±ℎ_𝑥,𝑦_𝐾±ℎ_𝑦)=0 and 𝐹(𝑦_𝐾±ℎ_𝑦)=0. Employing integral identity technique [16], we have _𝐾𝜕𝑢−𝐼_ℎ𝑢𝜑𝜕𝑥^1ℎ𝑑𝑥𝑑𝑦=𝜑_1ℎ𝐾𝐾𝐹^′′𝑦𝜕𝑢−𝐼_ℎ𝑢𝜕𝑥𝑑𝑥𝑑𝑦=𝜑_1ℎ𝐾𝐹^′𝑦_𝐾+ℎ_𝑦𝑙3𝜕𝑢−𝐼_ℎ𝑢𝜕𝑥𝑑𝑥−𝐹^′𝑦_𝐾−ℎ_𝑦𝑙1𝜕𝑢−𝐼_ℎ𝑢−𝜕𝑥𝑑𝑥_𝐾𝐹^′𝑦𝜕²𝑢−𝐼_ℎ𝑢𝜕𝑥𝜕𝑦𝑑𝑥𝑑𝑦=−𝜑_1ℎ𝐾𝐹𝑦_𝐾+ℎ_𝑦𝑙3𝜕²𝑢−𝐼_ℎ𝑢𝑦𝜕𝑥𝜕𝑦𝑑𝑥−𝐹_𝐾−ℎ_𝑦𝑙1𝜕²𝑢−𝐼_ℎ𝑢−𝜕𝑥𝜕𝑦𝑑𝑥_𝐾𝐹𝑦𝜕³𝑢−𝐼_ℎ𝑢𝜕𝑥𝜕𝑦²𝑑𝑥𝑑𝑦=𝜑_1ℎ𝐾𝐾𝐹𝑦𝜕³𝑢−𝐼_ℎ𝑢𝜕𝑥𝜕𝑦²𝑑𝑥𝑑𝑦≤𝐶ℎ²‖𝑢‖_3,𝐾𝜑_1ℎ0,𝐾. (27) Similarly, _𝐾𝜕𝑢−𝐼_ℎ𝑢𝜑𝜕𝑦^2ℎ𝑑𝑥𝑑𝑦≤𝐶ℎ²‖𝑢‖_3,𝐾𝜑_2ℎ0,𝐾. (28) Therefore, from (25)–(28), we have _{𝐾∈𝑇ℎ𝐾}∇𝑢−Π_ℎ𝑢⋅𝜑_ℎ=𝑑𝑥𝑑𝑦_{𝐾∈𝑇ℎ𝐾}∇𝑢−𝐼_ℎ𝑢⋅𝜑_ℎ𝑑𝑥𝑑𝑦=_{𝐾∈𝑇ℎ𝐾}𝜕𝑢−𝐼_ℎ𝑢𝜑𝜕𝑥^1ℎ𝑑𝑥𝑑𝑦+_𝐾𝜕𝑢−𝐼_ℎ𝑢𝜑𝜕𝑦^2ℎ𝑑𝑥𝑑𝑦≤𝐶ℎ²‖𝑢‖₃𝜑_ℎ0. (29) The proof is completed.

Now we start to state the following superclose property.

Theorem 3. Assume that (𝜓,𝑢) and (𝜓_ℎ,𝑢_ℎ) are the solutions of (3) and (10), respectively; 𝑢∈𝐻³(Ω), 𝜓∈(𝐻²(Ω))², there holds 𝜓_ℎ−𝐽_ℎ𝜓₀+𝑢_ℎ−Π_ℎ𝑢_ℎ≤𝐶ℎ²‖𝑢‖₃+𝜓₂. (30)

Proof. For (𝜑_ℎ,𝑣_ℎ)∈𝐻_ℎ×𝑀_ℎ, from (2) and (10), we have 𝑎𝜓,𝜑_ℎ+𝑏_ℎ𝜑_ℎ,𝑢=𝑎𝜓_ℎ,𝜑_ℎ+𝑏_ℎ𝜑_ℎ,𝑢_ℎ𝑏=0,_ℎ𝜓,𝑣_ℎ−𝑏_ℎ𝜓_ℎ,𝑣_ℎ=𝑏_ℎ𝜓,𝑣_ℎ−𝐺𝑣_ℎ=−_{𝐾∈𝑇ℎ𝐾}𝜓⋅∇𝑣_ℎ𝑑𝑥𝑑𝑦+_{𝐾∈𝑇ℎ𝐾}𝑓𝑣_ℎ𝑑𝑥𝑑𝑦=−_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ𝑑𝑥𝑑𝑦+_{𝐾∈𝑇ℎ𝐾}div𝜓+𝑓𝑣_ℎ𝑑𝑥𝑑𝑦=−_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ𝑑𝑥𝑑𝑦. (31) Applying (31) yields 𝑄_ℎ𝜓−𝜓_ℎ,𝑢−𝑢_ℎ,𝜑_ℎ,𝑣_ℎ=_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ𝑑𝑥𝑑𝑦,∀𝜑_ℎ,𝑣_ℎ∈𝐻_ℎ×𝑀_ℎ. (32) Using (12) in Lemma 1 and (32), we can obtain 𝜓_ℎ−𝐽_ℎ𝜓₀+𝑢_ℎ−Π_ℎ𝑢_ℎ≤𝐶sup_{(𝜑ℎ,𝑣ℎ)∈𝐻ℎ×𝑀ℎ}𝑄_ℎ𝜓_ℎ−𝐽_ℎ𝜓,𝑢_ℎ−Π_ℎ𝑢,𝜑_ℎ,𝑣_ℎ𝜑_ℎ0+𝑣_ℎℎ≤𝐶sup_{𝜑ℎ,𝑣ℎ∈𝐻ℎ×𝑀ℎ}𝑄_ℎ𝜓−𝐽_ℎ𝜓,𝑢−Π_ℎ𝑢,𝜑_ℎ,𝑣_ℎ−_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ×𝜑𝑑𝑥𝑑𝑦_ℎ0+𝑣_ℎℎ⁻¹=𝐶sup_{(𝜑ℎ,𝑣ℎ)∈𝐻ℎ×𝑀ℎ}𝑎𝜓−𝐽_ℎ𝜓,𝜑_ℎ+𝑏_ℎ𝜑_ℎ,𝑢−Π_ℎ𝑢−𝑏_ℎ𝜓−𝐽_ℎ𝜓,𝑣_ℎ−_{𝐾∈𝑇ℎ𝜕𝐾}𝜓⋅𝑛𝑣_ℎ×𝜑𝑑𝑥𝑑𝑦_ℎ0+𝑣_ℎℎ⁻¹. (33) Hence the desired result follows from the interpolation theorem and Lemma 2.
In order to derive global superconvergence for 𝑢 and flux 𝜓, we introduce the following postprocessing operators 𝐼_2ℎ∶𝑢∈(𝐻²(Ω)∩𝐻¹₀(Ω))↦𝐼_2ℎ𝑢∈𝑄₂(𝜅) and 𝐽_2ℎ∶𝜓∈(𝐻¹(Ω))²↦𝐽_2ℎ𝜓∈(𝑄₁(𝜅))² as 𝐼_2ℎ𝑢(𝑎_𝑖)=𝑢(𝑎_𝑖), 𝑖=1,…,9, and ∫_𝐾𝑗(𝜓−𝐽_2ℎ𝜓)=0, 𝑗=1,2,3,4, where 𝑢(𝑎_𝑖) are the value of 𝑢 on the nodes 𝑎_𝑖 and 𝑎_𝑖 are nodes of 𝑇_ℎ on macroelement 𝜅, while 𝜅∈𝑇_2ℎ consists of the four small elements 𝐾_𝑗 in 𝑇_ℎ(𝑗=1,2,3,4), and 𝑄₁(𝜅) and 𝑄₂(𝜅) are bilinear and biquadratic piecewise polynomials spaces, respectively. It can be checked that the following properties hold: 𝐼_2ℎΠ_ℎ𝑢=𝐼_2ℎ𝑢,𝐼_2ℎ𝑢−𝑢_ℎ≤𝐶ℎ²|𝑢|₃,𝐼_2ℎ𝑣_ℎℎ≤𝐶𝑣_ℎℎ,∀𝑣_ℎ∈𝑀_ℎ,𝐽_2ℎ𝐽_ℎ𝜓=𝐽_2ℎ𝜓,𝐽_2ℎ𝜓−𝜓₀≤𝐶ℎ²𝜓₂,𝐽_2ℎ𝜑_ℎ0≤𝐶𝜑_ℎ0,∀𝜑_ℎ∈𝐻_ℎ. (34) Then we can have the following superconvergence result.

Theorem 4. Under the assumptions in Theorem 3, there holds 𝑢−𝐼_2ℎ𝑢_ℎℎ+𝜓−𝐽_2ℎ𝜓_ℎ0≤𝐶ℎ²‖𝑢‖₃+𝜓₂. (35)

Proof. It follows from (30) in Theorem 3, (34), and the triangle inequality that 𝑢−𝐼_2ℎ𝑢_ℎℎ+𝜓−𝐽_2ℎ𝜓_ℎ0≤𝑢−𝐼_2ℎΠ_ℎ𝑢_ℎ+𝐼_2ℎΠ_ℎ𝑢−𝐼_2ℎ𝑢_ℎℎ+𝜓−𝐽_2ℎ𝐽_ℎ𝜓₀+𝐽_2ℎ𝐽_ℎ𝜓−𝐽_2ℎ𝜓_ℎ0≤𝑢−𝐼_2ℎΠ_ℎ𝑢_ℎ+𝐼_2ℎΠ_ℎ𝑢−𝑢_ℎℎ+𝜓−𝐽_2ℎ𝐽_ℎ𝜓₀+‖𝐽_2ℎ𝐽_ℎ𝜓−𝜓_ℎ0≤𝑢−𝐼_2ℎ𝑢_ℎ+𝐶Π_ℎ𝑢−𝑢_ℎℎ+𝜓−𝐽_2ℎ𝜓₀+𝐶𝐽_ℎ𝜓−𝜓_ℎ0≤𝐶ℎ²‖𝑢‖₃+𝜓₂, (36) which is the desired result.

3. Numerical Experiments

In this section, some numerical examples and comparison with other methods are presented to confirm theoretical analysis and good performance of NMFEM.

We consider the problem (1) with Ω=[0,1]² and the exact solution is 𝑢=𝑥𝑦(1−𝑥)(1−𝑦), and then the flux field 𝜓=(𝜓₁,𝜓₂) can be expressed as 𝜓₁=𝑦(1−2𝑥)(1−𝑦),𝜓₂=𝑥(1−𝑥)(1−2𝑦). We divide the domain Ω into a family of quasiuniform rectangles with number of 𝑁. The figures of exact solution (𝑢,𝜓) of problem (1) and finite element approximation (𝑢_ℎ,𝜓_ℎ) with 𝑁=64² are plotted in Figures 1, 2, and 3, respectively.

(a)

(b)

Figure 1 

The exact flux 𝜓1 (a) Numerical approximation 𝜓1ℎ (b).

(a)

(b)

Figure 2 

The exact flux 𝜓2 (a) and numerical approximation 𝜓2ℎ (b).

(a)

(b)

Figure 3 

The exact displacement 𝑢 (a) and numerical approximation 𝑢ℎ (b).

In Tables 1 and 2, we present the superclose and superconvergence results of the original variable 𝑢 in energy norm and flux 𝜓 in 𝐿² norm with 𝑁=8²;16²;32²;64²;128²;256², respectively. It is clearly that ‖𝜓_ℎ−𝐽_ℎ𝜓‖₀, ‖𝑢_ℎ−Π_ℎ𝑢‖_ℎ, ‖𝜓−𝐽_2ℎ𝜓_ℎ‖₀, and ‖𝑢−𝐼_2ℎ𝑢_ℎ‖_ℎ are converged at order 2 with respect to ℎ, which coincide with our theoretical analysis in Theorems 3 and 4. In order to describe the results more intuitively, we plot the errors in the logarithm scales in Figure 4.

Table 1 

Superclose results of 𝜓 and 𝑢.

Table 2 

Superconvergence results of 𝜓 and 𝑢.

(a)

(b)

Figure 4 

The superclose results (a) and the superconvergence results (b).

Moreover, we compare the results of NMFEM with those of FEM using 4-node quadrilateral (FEM-Q4) and MFEM (𝑢 and flux 𝜓 are approximated by piecewise constants and the Raviart-Thomas element, resp.).

The convergence rates of errors of 𝑢 in 𝐿² and energy norm are shown in Figures 5 and 6. The comparison of the flux 𝜓 in 𝐿² norm of NMFEM with MFEM is also given in Figure 7. It is observed that (a) the convergence rates of 𝑢 and flux 𝜓 in 𝐿² norm of NMFEM are better than those of MFEM; (b) NMFEM has almost the same rate as FEM-Q4 for 𝑢 in 𝐿² norm, but a higher rate than FEM-Q4 for 𝑢 in energy norm.

Figure 5 

Convergence rates of 𝑢 in 𝐿2 norm.

Figure 6 

Convergence rates of 𝑢 in energy norm.

Figure 7 

Convergence rates of the flux 𝜓 in 𝐿2 norm.

Furthermore, a numerical experiment is carried out to demonstrate the effectiveness and accuracy of NMFEM for the following diffusion problem: −div(𝐷∇𝑢)=𝑓,inΩ,𝑢=𝑢_𝜕,on𝜕Ω, (37) where 𝑓 is the source term and 𝑓∈𝐿²(Ω), 𝑢_𝜕 is the boundary data, the permeability 𝐷 is a symmetric tensor-valued function such that (a) 𝐷 is piecewise Lipschitz-continuous on Ω and (b) the set of the eigenvalues of 𝐷 is included in [𝜆_min,𝜆_max] (with 𝜆_min>0) for all (𝑥,𝑦)∈Ω.

As [13], we consider the problem (37) with Ω=[0,1]², 𝐷=1001if𝑥≤0.5,𝐷=100000.01if𝑥>0.5, (38) and the analytical solution ()()𝑢=cos𝜋𝑥sin𝜋𝑦if𝑥≤0.5,𝑢=0.01cos𝜋𝑥sin𝜋𝑦if𝑥>0.5. (39) Let 𝜓=𝐷∇𝑢, we can get the superclose and superconvergence results of 𝑢 in energy norm and flux 𝜓 in 𝐿² norm. The errors are listed in Tables 3 and 4 and plotted in the logarithm scales in Figure 8, respectively.

Table 3 

Superclose results of 𝜓 and 𝑢.

Table 4 

Superconvergence results of 𝜓 and 𝑢.

(a)

(b)

Figure 8 

The superclose results (a) and the superconvergence results (b).

Acknowledgments

The authors would like to express their sincere thanks to the anonymous referee for his many helpful suggestions, which contribute significantly to the improvement of the paper. The research is supported by the NSF of China (no. 10971203; no. 11271340), Research Fund for the Doctoral Program of Higher Education of China (no. 20094101110006), and Foundation of He’nan Educational Committee (no. 13B110144).