Simplified weak Galerkin and new finite difference schemes for the Stokes equation

Journal of Computational and Applied Mathematics 361 (2019) 176–206

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Simplified weak Galerkin and new finite difference schemes for the Stokes equation ∗

Yujie Liu a , ,1 , Junping Wang b ,2 a b

School of Data and Computer Science, Sun Yat-sen University, Guangzhou, 510275, China Division of Mathematical Sciences, National Science Foundation, Alexandria, VA 22314, United States

article

info

a b s t r a c t

Article history: Received 5 August 2018 Received in revised form 8 April 2019 Keywords: Finite difference methods Superconvergence Simplified weak Galerkin Stokes equations

This article presents a simplified formulation for the weak Galerkin finite element method for the Stokes equations without using the degrees of freedom associated with the unknowns in the interior of each element as formulated in the original weak Galerkin algorithm. The simplified formulation preserves the important mass conservation property locally on each element and allows the use of general polygonal partitions. A particular application of the simplified weak Galerkin on rectangular partitions yields a new class of 5- and 7-point finite difference schemes for the Stokes equations. An explicit formula is presented for the computation of the element stiffness matrices on arbitrary polygonal elements. Error estimates of optimal order are established for the simplified weak Galerkin finite element method in the H 1 and L2 norms. Furthermore, a superconvergence of order O(h3/2 ) is established on rectangular partitions for the velocity approximation in the H 1 norm at cell centers, and a similar superconvergence is derived for the pressure approximation in the L2 norm at cell centers. Some numerical results are reported to confirm the convergence and superconvergence theory. © 2019 Elsevier B.V. All rights reserved.

1. Introduction This paper is concerned with the development of a simplified formulation for the weak Galerkin finite element method for the Stokes equations developed in [1]. For simplicity, consider the following 2D Stokes problem which seeks a velocity field u and a pressure unknown p satisfying

{ −∆u + ∇ p divu u

=f =0 =0

in Ω , in Ω , on Γ ,

(1)

(1) (2) t where Ω is a bounded polygonal domain with boundary Γ = ∫ ∂ Ω , f = (f , f ) represents external source, and the pressure function p is assumed to have mean value zero; i.e., Ω pdΩ = 0.

∗ Corresponding author. E-mail addresses: [email protected] (Y. Liu), [email protected] (J. Wang). 1 The research of Liu was partially supported by Guangdong Provincial Natural Science Foundation (No. 2017A030310285), Shandong Provincial natural Science Foundation (No. ZR2016AB15) and Youthful Teacher Foster Plan Of Sun Yat-Sen University (No. 171gpy118). 2 The research of Wang was supported in part by the NSF IR/D program, while working at National Science Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. https://doi.org/10.1016/j.cam.2019.04.024 0377-0427/© 2019 Elsevier B.V. All rights reserved.

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

177

The weak Galerkin finite element method (WG-FEM) is a natural extension of the standard conforming finite element method where differential operators are approximated as discrete distributions or discrete weak derivatives. WG-FEM has the flexibility of dealing with discontinuous elements while sharing the same weak formulation with the classical conforming finite element methods. Following its first development in [2,3] for second order elliptic equations, the method has gained a lot attention and popularity and has been applied to various PDEs including the Stokes equations, biharmonic equation, elasticity equation, and Maxwell’s equations, see [4–6] and the references therein for more details. In this paper, we propose a simplified formulation for the weak Galerkin finite element method for the Stokes equations. In the simplified weak Galerkin (SWG), the degrees of freedom involves only those on the element boundary; i.e., the unknowns associated with the interior of each element in the original weak Galerkin are not used in SWG. This simplified formulation still preserves the important mass conservation property locally on each element and allows the use of general polygonal partitions. The stability and convergence is established by following the general framework developed in [1], but with a non-trivial modification due to the absence of the interior unknowns employed in the original weak Galerkin. A comprehensive formula for the computation of the element stiffness matrices is presented for a better understanding and dissemination of the method. As a particular case study, the global stiffness matrix for the simplified weak Galerkin is assembled on uniform rectangular partitions, yielding a new class of 5- and 7-point finite difference schemes for the Stokes equations. Furthermore, a superconvergence is established in the H 1 norm for the velocity approximation at cell centers, as well as for the pressure approximation in the L2 norm at cell centers on (nonuniform) rectangular partitions. A section is devoted to the description of the new class of 5- and 7-point finite difference schemes in order to facilitate the practical use of this scheme. Due to the connection between the finite difference scheme and the simplified weak Galerkin finite element method, the convergence theory developed for SWG can be easily extended to the new finite difference schemes. The local conservation property is also inherited by the finite difference scheme. The SWG method indeed provides a way to generalize the new finite difference scheme from regular Cartesian grids to irregular polygonal grids. There are existing finite difference schemes in literature for solving the Stokes equations [7–9]. The MAC scheme [7] is one of the finite difference methods that has enjoyed a lot attention and popularity in computational fluid dynamics. Nicolaides [10] showed that the MAC method can be seen as a special case of the co-volume formulation. Han and Wu [11] derived the MAC scheme from a newly proposed mixed finite element method, in which the two components of the velocity and the pressure are defined on different meshes (an approach that resembles the MAC idea). Recently, Kanschat [12] showed that the MAC scheme is algebraically equivalent to the divergence-conforming discontinuous Galerkin method with a modified interior penalty formulation. In [13], Minev demonstrated that the MAC scheme can be interpreted within the framework of the local discontinuous Galerkin (LDG) methods. The 5- and 7-point finite difference schemes of this paper are based on different stencils from the MAC scheme, and therefore represent a new class of numerical schemes which are stable and have error estimates of optimal order. Finite difference schemes can also be obtained from the finite volume element methods (FVEM) [14–16]. The FVEM also known as finite volume methods [17,18], covolume methods [19], combined finite volume-finite element methods [20,21] or box methods [22,23], have been extensively developed for the Stokes equations [11,24–26]. In these approaches, the numerical methods were largely devised by using a volume integral formulation of the problem with a finite partitioning set of volumes, with admissible functions restricted to a finite element space. A complementary dual mesh is often constructed by connecting the barycenters of the triangles in the FE mesh, with the midpoints of the associated edges [16,20], or the vertices of the triangles [17,27] for 2D problems. Compared with the FVEM, our scheme is based on polygonal partitions which do not require any dual mesh. Our finite element functions are defined on the boundary of each polygonal element, while the FVEMs are often based on a classical conforming finite element space with test functions defined on the dual mesh. SWG preserves the important mass conservation property, while FVEMs may or may not retain the mass conservation property, depending on the selection of the pressure space. The paper is organized as follows: In Section 2, we present a simplified version of the weak Galerkin finite element method on arbitrary polygonal partitions. In Section 3, we present a formula with comprehensive derivation for the computation of the element stiffness matrices in SWG. In Section 4, we apply the SWG scheme to uniform rectangular partitions and derive a class of 5- and 7-point finite difference schemes. In particular, Section 4.2 is devoted to a presentation of the class of 5-point finite difference scheme for the model problem (1). In Section 5, we provide a stability analysis for the simplified weak Galerkin. In Section 6, we derive some error estimates for the SWG approximations in various Sobolev norms. In Section 7, a superconvergence result is developed for the numerical velocity and pressure arising from rectangular partitions. Finally, in Section 8, we present a few numerical results to demonstrate the efficiency and accuracy of SWG. 2. A simplified weak Galerkin Assume that the domain is of polygonal type and is partitioned into unstructured polygons Th = {T }. Let T ∈ Th be a polygonal element of N sides (e.g., a hexahedron shown in Fig. 1). Denote by {ei }Ni=1 the edge set of T . Denote by Mi the midpoint of the edge ei , and ni the outward normal direction of ei ; see Fig. 1.

178

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 1. An illustrative hexagonal element.

Given any piecewise constant function ub on the boundary of T ∈ Th ; i.e., ub |ei = ub,i , we define the weak gradient of ub on T by

∇w ub :=

N 1 ∑

|T |

ub,i |ei |ni ,

(2)

i=1

where |ei | is the length of the edge ei , |T | is the area of the element T . For convenience, we denote by W (T ) the space of piecewise constant functions on ∂ T . The global finite element space W (Th ) is constructed by patching together all the local elements W (T ) through single values on interior edges. The subspace of W (Th ) consisting of functions with vanishing boundary value is denoted as W0 (Th ). We use the conventional notation of Pj (T ) for the space of polynomials of degree j ≥ 0 on T . For each ub ∈ W (T ), we associate it with a linear extension in T , denoted by s(ub ) ∈ P1 (T ), satisfying N ∑

(s(ub )(Mi ) − ub,i )v (Mi )|ei | = 0 ∀ v ∈ P1 (T ).

(3)

i=1

It is easy to see that s(ub ) is well defined by (3), and its computation is simple and local. In fact, s(ub ) can be viewed as an extension of ub from ∂ T to T through a least-squares fitting. On each element T ∈ Th , we introduce the following stabilizer: ST (ub , vb ) := h−1

N ∑

(s(ub )(Mi ) − ub,i )(s(vb )(Mi ) − vb,i )|ei |

i=1

(4)

= h−1 ⟨Qb s(ub ) − ub , Qb s(vb ) − vb ⟩∂ T for ub , vb ∈ W (T ), where Qb is the L2 projection operator onto W (T ) and ⟨·, ·⟩∂ T stands for the usual L2 inner product in L2 (∂ T ). The weak gradient and the local stabilizer defined in (2) and (4) can be extended naturally to vector-valued function, for example in the two dimensional space we would have

[

] ∇ w ub ∇w ub := , ∇w v b

(5)

ST (ub , wb ) := ST (ub , wb ) + ST (vb , zb )

(6)

wb ∈ [W (T )]2 and wb = ∈ [W (T )]2 . vb zb For ub ∈ [W (T )]2 , the weak divergence of ub is defined as a constant on T satisfying the following equation:

[ ]

for ub =

[

ub

]

(∇w · ub , 1)T := ⟨ub · n, 1⟩∂ T .

(7)

The weak divergence is therefore given by

∇w · ub |T =

N 1 ∑

|T |

i=1

ub |ei ·ni |ei |.

(8)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

179

The Simplified Weak Galerkin (SWG) method for the Stokes equations (1) seeks ub ∈ [W0 (Th )]2 , ph ∈ P0 (Th ) such that

{ κ S(ub , vb ) + (∇w ub , ∇w vb ) − (∇w · vb , ph ) = (f , s(vb )), (9) (∇w · ub , w ) = 0, ∑ for all vb ∈ [W0 (Th )]2 and w ∈ P0 (Th ), where S(ub , vb ) := T ST (ub , vb ) and κ > 0 is any prescribed stabilization parameter. Remark 2.1. The solution ub of the numerical algorithm (9) is identical to the second component of the solution {u0 , ub } obtained from the weak Galerkin finite element method introduced in [1] when the lowest order of WG elements are employed. A proof can be given by a straightforward elimination of the variable u0 from the corresponding numerical scheme with details omitted. Such an elimination is not known to be possible for WG elements of higher orders. 3. On the computation of element stiffness matrices The simplified weak Galerkin finite element method (9) can be reformulated as follows:

⎧ ⎛ ⎡ ⎤ ⎞ ⎪ ∑ ∑ ⎪ ⎪ w ⎪ ⎝κ ST (ub , wb ) + (∇w ub , ∇w wb )T − (∇w · ⎣ b ⎦ , ph )T ⎠ = (f (1) , s(wb ))T , ⎪ ⎪ ⎪ ⎪ T ⎪ T 0 ⎪ ⎪ ⎞ ⎡ ⎤ ⎛ ⎪ ⎪ ⎨ ∑ ∑ ⎝κ ST (vb , zb ) + (∇w vb , ∇w zb )T − (∇w · ⎣ 0 ⎦ , ph )T ⎠ = (f (2) , s(zb ))T , ⎪ ⎪ T T ⎪ zb ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ ⎪ ∑ ⎪ u ⎪ ⎪ (∇w · ⎣ b ⎦ , χ )T = 0, ⎪ ⎪ ⎩ T vb

(10)

for all wb , zb ∈ W0 (Th ) and χ ∈ P0 (Th ). From (10), the element stiffness matrix on T ∈ Th consists of three sub-matrices corresponding to the following forms: ST (ub , wb ), (∇w ub , ∇w wb )T , and (∇w ·

[ ] ub

vb

, χ )T .

The goal of this section is to derive the matrix for each of the three forms. Denote by Xub , Xub , and P the vector representation of ub , vb , and ph given by ub,1 ⎢ ub,2 ⎥

⎡

⎤ vb,1 ⎢ vb,2 ⎥ ⎥ =⎢ ⎣ .. ⎦ , and P = [ph ]. . v b ,N ⎡

⎤

⎥ Xu b = ⎢ ⎣ .. ⎦ , Xvb .

ub,N

Here ub,i and vb,i , i = 1 . . . N, represent the values of ub and vb at the midpoint Mi of the edge ei , and P is the value of ph at the center of T . The following theorem gives a computational formula for the element stiffness matrix and the element load vector. Theorem 1. as follows:

The element stiffness matrix and the element load vector for the SWG scheme are given in a block matrix form

[ κA + B 0 Q1t

0

Q1 Q2 0

κA + B Q2t

][

Xub Xvb P

]

F1

[ ] ∼ = F2 , 0

where the block components in (11) are computed explicitly by: A = {ai,j }N ×N = h−1 (E − EM(M t EM)−1 M t E), B = {bi,j }N ×N , bi,j = ni · nj

|ei ||ej | , |T |

D = {di,j }3×N = (M t EM)−1 M t E , (1)

F1 = {fi }N ×1 , fi

(1)

∫

f (1) (x, y)(d1,i + d2,i (x − xT ) + d3,i (y − yT ))dT ,

= T

(2)

F2 = {fi }N ×1 , fi

(2)

∫

f (2) (x, y)(d1,i + d2,i (x − xT ) + d3,i (y − yT ))dT ,

= T

(11)

180

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Q1 = {qi }N ×1 , qi = |ei | Q2 = {qi }N ×1 , qi = |ei |

[ ]

1 · ni , 0

[ ]

0 · ni . 1

Here Q1t and Q2t stand for the transpose of Q1 and Q2 , respectively. The matrices M and E are given by 1 ⎢1

⎡

M=⎢ ⎣ ..

.

1

x1 − xT x2 − xT

.. . xN − xT

y1 − yT y2 − yT ⎥

⎤

⎡

|e1 |

⎢ , E=⎢ .. ⎥ ⎦ ⎣ . yN − yT N ×3

⎤ |e2 |

..

⎥ ⎥ ⎦

.

|eN |

,

N ×N

where MT = (xT , yT ) is any point on the plane (e.g., the center of T as a specific case), (xi , yi ) is the midpoint of ei , |ei | is the length of edge ei , ni is the outward normal vector on ei , and |T | is the area of the element T . The rest of this section is devoted to a derivation of the formula (11). 3.1. Element stiffness matrix for ST (ub , wb ) For any ub ∈ W (T ), the linear extension of ub , s(ub ), can be represented as follows:

s(ub ) = α0 + α1 (x − xT ) + α2 (y − yT ). By (3), we have N ∑

(s(ub )(Mi ) − ub,i )φ (Mi )|ei | = 0 ∀ φ ∈ P1 (T ).

i=1

It follows that N ∑ ) ( α0 + α1 (xi − xT ) + α2 (yi − yT ) − ub,i φ (Mi )|ei | = 0 ∀ φ ∈ P1 (T ), i=1

where (xi , yi ) is the midpoint of ei . In particular, we may choose φ = 1, x − xT , y − yT to obtain

⎧ N ∑( ) ⎪ ⎪ ⎪ α0 + α1 (xi − xT ) + α2 (yi − yT ) − ub,i |ei | = 0, ⎪ ⎪ ⎪ ⎪ i=1 ⎪ ⎪ N ⎨∑ ( ) α0 + α1 (xi − xT ) + α2 (yi − yT ) − ub,i (xi − xT )|ei | = 0, ⎪ ⎪ i=1 ⎪ ⎪ ⎪ N ⎪ ∑ ( ) ⎪ ⎪ ⎪ α0 + α1 (xi − xT ) + α2 (yi − yT ) − ub,i (yi − yT )|ei | = 0, ⎩ i=1

which gives

⎧ N N ∑ ∑ ⎪ ⎪ ⎪ ub,i |ei |, (α 0 + α1 (xi − xT ) + α2 (yi − yT )) |ei | = ⎪ ⎪ ⎪ ⎪ i=1 i=1 ⎪ ⎪ N N ⎨∑ ∑ ub,i (xi − xT )|ei |, (α0 + α1 (xi − xT ) + α2 (yi − yT )) (xi − xT )|ei | = ⎪ ⎪ i=1 i=1 ⎪ ⎪ ⎪ N N ⎪ ∑ ∑ ⎪ ⎪ ⎪ ub,i (yi − yT )|ei |. (α0 + α1 (xi − xT ) + α2 (yi − yT )) (yi − yT )|ei | = ⎩ i=1

i=1

Then (12) can be reformulated as

⎡ ⎤ ub , 1 [ ] α0 ⎢ ub , 2 ⎥ ⎥ M t EM α1 = M t E ⎢ ⎣ .. ⎦ , . α2 ub,N

(12)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

181

which leads to

⎡ ⎤ ub,1 [ ] α0 ⎢ ub,2 ⎥ ⎥ α1 = (M t EM)−1 M t E ⎢ ⎣ .. ⎦ . . α2

(13)

ub,N

Since s(ub ) = α0 + α1 (x − xT ) + α2 (y − yT ), we have

s(ub )(M1 ) ⎢ s(ub )(M2 ) ⎥

⎡

1 ⎢1

⎤

.. .

⎢ ⎣

x1 − xT x2 − xT

⎡

⎥ = ⎢. ⎦ ⎣. .

s(ub )(MN )

1

⎡

.. . xN − xT ⎤

y1 − yT [ ] y2 − yT ⎥ α0

⎤

.. ⎥ ⎦ α1 . α2 yN − yT

[ ] α0 = M α1 α2

ub,1

⎢ ub,2 ⎥ ⎥ = M(M t EM)−1 M t E ⎢ ⎣ .. ⎦ . . ub,N Let wb be the basis function of W (T ) corresponding to the edge ej of T ; i.e.,

wb =

{

1, on ej , 0, on other edges,

that is wb,j = 1, wb,{1,...,N }/j = 0, then the coefficient (α0 , α1 , α2 )t for s(wb ) is given by

⎡w

b,1

⎡0⎤

⎤

⎢ .. ⎥ [ ] ⎢ .. ⎥ [ ] ⎢.⎥ ⎢ . ⎥ d1,j α0 ⎢ ⎥ ⎥ ⎢ α1 = (M t EM)−1 M t E ⎢ wb,j ⎥ = (M t EM)−1 M t E ⎢1⎥ ≜ d2,j . ⎢.⎥ ⎢ . ⎥ d3,j α2 ⎣.⎦ ⎣ . ⎦ . . wb,N

0

Thus, we have ST (ub , wb ) =

N ∑

h−1 (s(ub )(Mi ) − ub,i )(s(wb )(Mi ) − wb,i )|ei |

i=1

= h− 1

N ∑

(ub,i − s(ub )(Mi ))wb,i |ei | (14)

i=1

⎛

⎡

ub,1

⎤⎞

⎜ ⎢ ub,2 ⎥⎟ t −1 t ⎢ . ⎥⎟ · wb,j |ej | = h−1 ⎜ (I − M(M EM) M E) N ⎝ ⎣ . ⎦⎠ . ub,N

j

where IN is the N × N identity matrix. The identity (14) gives rise to the element stiffness matrix for the term ST (ub , wb ). The result can be summarized as follows. Theorem 2. For the polygonal element T depicted in Fig. 1, the element stiffness matrix corresponding to the bilinear form ST (ub , wb ) is given by A = {ai,j }Ni,j=1 := h−1 (E − EM(M t EM)−1 M t E),

(15)

where 1 ⎢1

⎡

M=⎢ ⎣ ..

.

1

x1 − xT x2 − xT

.. . xN − xT

y1 − yT y2 − yT ⎥

⎤

.. ⎥ ⎦ . yN − yT N ×3

⎡ |e1 | ⎢ , E=⎢ ⎣

⎤ |e2 |

..

.

,

⎥ ⎥ ⎦ |eN |

N ×N

MT = (xT , yT ) is the center of T or any point on the plane, (xi , yi ) is the midpoint of ei , |ei | is the length of edge ei , ni is the outward normal vector, |T | is the area of the element T.

182

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

3.2. Element stiffness matrix for (∇w ub , ∇w wb )T From the weak gradient formulation (2), we have (∇w ub , ∇w wb )T = (

=(

=

N 1 ∑

|T |

ub,i ni |ei |,

i=1

N ∑ |ei ||ej | i=1

=

i=1

N 1 ∑

|T |

ub,i ni |ei |,

N ∑

|T |

N 1 ∑

|T | 1

|T |

wb,j nj |ej |)T

(16)

j=1

wb,j nj |ej |)T

(ni · nj )ub,i wb,j ,

bi,j ub,i wb,j .

i=1

Thus, the corresponding element stiffness matrix is given by B := {bi,j }Ni,j=1 , bi,j = (ni · nj )

|ei ||ej | . |T |

(17)

3.3. Element stiffness matrix for (∇w · ub , χ )T Note that χ ∈ P0 (T ) is a constant. Thus, we have (∇ w · u b , χ ) T =

∫

∇w · ub χ dT = χ

T N

=

∫ ∂T

ub · nds

(18)

∑ [u ] b,i · ni |ei |χ vb,i i=1

= Q1t Xub χ + Q2t Xvb χ,

(19)

which gives the contribution of Q1 and Q2 in the element stiffness matrix. 3.4. Element load vector For a computation of the element load vector corresponding to (f (1) , s(wb ))T , let wb be the basis function of W (T ) corresponding to edge ej of T ; i.e.,

wb =

{

1, on ej , 0, on other edges.

The linear extension of wb is given by

s(wb ) = d1,j + d2,j (x − xT ) + d3,j (y − yT ), where

⎡w

b ,1

d1,j d2,j d3,j

[

]

⎤

⎡0⎤

⎢ .. ⎥ ⎢ .. ⎥ ⎢ . ⎥ ⎢.⎥ ⎢ ⎥ t −1 t ⎢ ⎥ = (M EM) M E ⎢ wb,j ⎥ = (M EM) M E ⎢1⎥ . ⎢ . ⎥ ⎢.⎥ ⎣ . ⎦ ⎣.⎦ . . wb,N 0 t

−1

t

It follows that (f (1) , s(wb ))T =

∫

f (1) s(wb )dT T

∫

(20) f (1) (x, y)(d1,j + d2,j (x − xT ) + d3,j (y − yT ))dT ,

= T

from which an element load vector can be easily computed. In practical computation, the integral in (20) should be approximated numerically by some quadrature rules.

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

183

4. SWG on Cartesian grids Consider the Stokes problem (1) defined on the unit square domain Ω = (0, 1) × (0, 1). A uniform mesh of the domain can be given by the Cartesian product of two one-dimensional grids: xi+µ = (i + µ − 0.5)h,

i = 0, 1, . . . , n,

yj+µ = (j + µ − 0.5)h,

j = 0 , 1 , . . . , n,

µ = 0, 1/2, µ = 0, 1/2,

where n is a positive integer and h = 1/n is the meshsize. Denote by Tij := [xi− 1 , xi+ 1 ] × [yj− 1 , yj+ 1 ], 2

2

2

2

i, j = 1, . . . , n

the square element centered at (xi , yj ) for i, j = 1, . . . , n. The collection of all such elements forms a uniform square partition Th of the domain. The collection of all the element edges is denoted as Eh . The goal of this section is to assemble the global stiffness matrix and the load vector for the SWG scheme associated with uniform rectangular partitions Th . The resulting matrix problem can be seen as a finite difference scheme for the Stokes equations. In particular, this will result in a new 5-point finite difference scheme when the stabilization parameter has the value of κ = 4. Let T ∈ Th be a rectangular element depicted in Fig. 2 with center MT = (xT , yT ). From (13), the linear extension of ub ∈ W (T ) (i.e., s(ub )) can be verified to be (see Lemma 6.1 in [6] for details):

s(ub ) = α0 + α1 (x − xT ) + α2 (y − yT ),

(21)

where

⎧ |e1 |(ub,1 + ub,2 ) + |e3 |(ub,3 + ub,4 ) ⎪ ⎪ α = , ⎪ ⎨ 0 2|e1 | + 2|e3 | α1 = (ub,2 − ub,1 )/|e3 |, ⎪ ⎪ ⎪ ⎩ α2 = (ub,4 − ub,3 )/|e1 |.

(22)

From (21) we have (s(vb ) − vb )(M1 ) = (s(vb ) − vb )(M2 ) |e3 | = (vb3 + vb4 − vb1 − vb2 ), 2(|e1 | + |e3 |)

(23)

(s(vb ) − vb )(M3 ) = (s(vb ) − vb )(M4 ) |e1 | =− (vb3 + vb4 − vb1 − vb2 ). 2(|e1 | + |e3 |)

(24)

Hence,

|e1 |(s(vb ) − vb )(M1 ) = −|e3 |(s(vb ) − vb )(M3 ).

(25)

A by-product of (22) is the following estimate:

∥∇ s(v )∥T ≤ C ∥∇w v∥T .

(26)

In addition, the linear function s(vb ) − vb satisfies (s(vb ) − vb )|e1 (y) = (s(vb ) − vb )|e2 (y), (s(vb ) − vb )|e3 (x) = (s(vb ) − vb )|e4 (x).

(27)

4.1. A 7-point finite difference scheme Two sets of grid points can be generated by using the 1-d grid points {xi } and {yj }. The first consists of the mid-points of all edges in Eh (i.e., the dotted points colored in red in Fig. 3), which are used to approximate the velocity field. The second set of grid points consists of all cell centers (i.e., the dotted points colored in blue in Fig. 3), which are used to approximate the pressure unknown. Let ukℓ = (ukℓ , vkℓ )t be the velocity approximation at the red-dotted grid points (xk , yℓ ) and pij be the approximate pressure at the blue-dotted grid points (xi , yj ). It can be seen that a red-dotted grid point (xk , yℓ ) is located on the boundary if either k or ℓ takes the value of 21 or n + 12 . The following is the main result of this section.

184

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 2. An illustrative square element.

Fig. 3. Grid points for the Stokes equations: (a) for the velocity (left, red dots), and (b) for the pressure (right, blue dots) . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Theorem 3. The simplified weak Galerkin scheme (9) is algebraically equivalent to the following finite difference scheme on uniform square partitions:

⎧ c1 ui+ 3 ,j + c2 ui+ 1 ,j + c3 ui− 1 ,j + c4 [ui+1,j− 1 + ui+1,j+ 1 + ui,j− 1 + ui,j+ 1 ] ⎪ ⎪ 2 2 2 2 2 ⎪ ]2 [ 2 ⎪ ⎪ ⎪ h2 p − p ⎪ i + 1 , j i , j ⎪ = f i + 1 ,j , +h ⎪ ⎪ 0 2 2 ⎪ ⎨ c1 ui,j+ 3 + c2 ui,j+ 1 + c3 ui,j− 1 + c4 [ui− 1 ,j+1 + ui+ 1 ,j+1 + ui− 1 ,j + ui+ 1 ,j ] 2 2 2 2 2 ⎪ [ 2 ]2 2 ⎪ ⎪ ⎪ h 0 ⎪ ⎪ = f i,j+ 1 +h ⎪ ⎪ pi,j+1 − pi,j 2 2 ⎪ ⎪ ⎪ ⎩ h[u 1 − u 1 ] + h[v 1 − v 1 ] = 0. i+ 2 ,j

i− 2 ,j

κ

i,j+ 2

(28)

i,j− 2

κ

κ

− 1), c2 = + 2, c4 = − , κ > 0 is the stabilization parameter. For κ = 4, we have c1 = c3 = 0, 4 2 4 c2 = 4 and c4 = −1 so that the scheme (28) gives rise to the five-point finite difference scheme (43).

where c1 = c3 = (

The rest of this subsection is devoted to a derivation of the finite difference scheme (28). To this end, let wb be the basis function of W (T ) corresponding to edge e1 of T (see Fig. 2) so that

wb =

{

1, on e1 , 0, on ei , i = 2, 3, 4.

On square elements where |e1 | = |e3 |, from (24) it is not hard to see that (wb − s(wb ))|M1 = (wb − s(wb ))|M2 =

1 4

,

1 (wb − s(wb ))|M3 = (wb − s(wb ))|M4 = − . 4

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

185

Thus, we have

⟨ub − s(ub ), wb − s(wb )⟩∂ T = ⟨ub , wb − s(wb )⟩∂ T 4 ∑ = |ei |ub,i (wb − s(wb ))|Mi

(29)

i=1

h

=

( (∇w ub , ∇w wb )T =

4 1 ∑

|T | ( [ i=1

4

(ub,1 + ub,2 − ub,3 − ub,4 ),

ub,i |ei |ni ,

]

1 ub,2 − ub,1 h ub,4 − ub,3

=

4 1 ∑

|T | ,

) wb,i |ei |ni

(30)

i=1

[

1 0−1 h 0−0

]) T

( )2 1

= −|T |

h − ub,2 ,

(f (1) , s(wb ))T =

= ub,1 ( ∫ f (1)

T 2

=

h

[

6

3 4

1 4

(ub,2 − ub,1 )

) 1 − (x − xT ) dT

(31)

h

1 f (1) |M1 +f (1) |MT − f (1) |M2 4

]

+ O(h4 ).

Analogously, the same techniques can be applied to the computation of ST (vb , wb ), (∇w vb , ∇w wb )T , and (f (2) , s(wb ))T . Next, let vb be the vector basis function of [W (T )]2 associated with edge e1 , that is

vb =

[ ] wb 0

, or vb =

[

0

]

wb

, with wb =

{

1, on e1 , 0, on other edges.

Then we have (∇w · vb , ph )T = ph |e1 |

[ ]

[ ]

1 0 · n1 , or = ph |e1 | · n1 0 1

(32)

= −hph , or = 0. (∇ w · u b , χ ) T =

N ∑

N

[ ] ub,i

[ ]

∑ 0 1 · ni |ei | · ni |ei | + vb,i 1 0

(33)

i=1

i=1

= h(ub,2 − ub,1 ) + h(vb,4 − vb,3 ). Here χ is the characteristic function of the element T . Eqs. (29), (30), (31), (32) and (33) comprise the discrete scheme corresponding to the basis function wb on edge e1 of the element T :

⎧κ ⎪ (ub,1 + ub,2 − ub,3 − ub,4 ) + ub,1 − ub,2 − hph = (f (1) , s(wb ))T ⎪ ⎪ ⎨4 κ (vb,1 + vb,2 − vb,3 − vb,4 ) + vb,1 − vb,2 = (f (2) , s(wb ))T ⎪ ⎪4 ⎪ ⎩ h(ub,2 − ub,1 ) + h(vb,4 − vb,3 ) = 0

(34)

The two local equations (34) corresponding to the elements Ti,j and Ti+1,j that share ei+ 1 ,j as a common edge (see 2 Fig. 4(a)) are given by

κ 4 ∫

(ui+ 1 ,j + ui+ 3 ,j − ui+1,j− 1 − ui+1,j+ 1 ) + ui+ 1 ,j − ui+ 3 ,j + hpi+1,j 2

2

2

2

2

2

(35)

f (1) s(wb )dx

= Ti+1,j

and

κ 4 ∫

(ui+ 1 ,j + ui− 1 ,j − ui,j− 1 − ui,j+ 1 ) + ui+ 1 ,j − ui− 1 ,j − hpi,j 2

2

f (1) s(wb )dx.

= Ti,j

2

2

2

2

(36)

186

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 4. Stencils for the finite difference scheme (28).

Summing up Eqs. (35) and (36) yields the following global linear equation corresponding to the degree of freedom ui+ 1 ,j : 2

c1 ui+ 3 ,j + c2 ui+ 1 ,j + c3 ui− 1 ,j + c4 [ui+1,j− 1 + ui+1,j+ 1 + ui,j− 1 + ui,j+ 1 ] 2

2

2

2

2

2

2

+ h(pi+1,j − pi,j )

∫

(37)

f (1) s(wb )dT ,

= Ti,j

⋃

Ti+1,j

κ

κ

κ

− 1, c2 = + 2, c4 = − . 4 2 4 The right-hand side of (37) can be approximated by using proper numerical integrations (the Simpson rule in the

where c1 = c3 =

x-direction and the midpoint rule in the y- direction) as follows:

∫

f (1) s(wb )dT Ti,j

⋃

∫ f

=

Ti+1,j

(1)

=

h

+ = = =

[f (1) s(wb )|M

6

h

h2

2

6 h2 24 h2 2

∫

2

Ti,j 2

s(wb,(i+ 1 ,j) )|Ti,j dT +

6

i− 1 ,j 2

[f (1) s(wb )|M

Ti+1,j

2

+4f (1) s(wb )|Mi,j +f (1) s(wb )|M

i + 1 ,j 2

1

3

i− 2 ,j

i + 1 ,j 2

]

+4f (1) s(wb )|Mi+1,j +f (1) s(wb )|M

(1) [− f (1)1 + fi,(1) f 1 ]+ j +

4

f (1) s(wb,(i+ 1 ,j) )|Ti+1,j dT

4

i+ 2 ,j

h

2

3

i+ 3 ,j 2

1

(1) (1) [ f (1)1 + fi+ f 3 ] + O(h4 ) 1,j −

6 4

i+ 2 ,j

4

(1) (1) (1) 4 [−f (1)1 + 4fi,(1) j + 6f 1 + 4fi+1,j − f 3 ] + O (h ) i− 2 ,j

f

(1) i+ 12 ,j

i+ 2 ,j

+ O(h4 ),

where we have used the fact that

s(wb,(i+ 1 ,j) )|Ti,j = 2

s(wb,(i+ 1 ,j) )|Ti+1,j = 2

1 4 1 4

1

+ (x − xTi,j ), h 1

− (x − xTi+1,j ). h

] + O(h4 )

i+ 2 ,j

i+ 2 ,j

(38)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

187

It follows that (37) can be rewritten as c1 ui+ 3 ,j + c2 ui+ 1 ,j + c3 ui− 1 ,j + c4 [ui+1,j− 1 + ui+1,j+ 1 + ui,j− 1 + ui,j+ 1 ] 2

2

2

2

2

2

2

+ h(pi+1,j − pi,j ) =

h2 2

f

(1)

(39)

+ O(h ). 4

i+ 12 ,j

Analogously, for the degree of freedom ui,j+ 1 , we may obtain 2

c1 ui,j+ 3 + c2 ui,j+ 1 + c3 ui,j− 1 + c4 [ui− 1 ,j+1 + ui+ 1 ,j+1 + ui− 1 ,j + ui+ 1 ,j ] 2

2

h2

=

f

(1)

2

2

2

2

(40)

+ O(h ).

i,j+ 12

2

2

4

The stencil for the unknown u, or more precisely for ui+ 1 ,j (respectively ui,j+ 1 ) is the seven dotted-points shown in 2 2 Fig. 4(a) (respectively Fig. 4(b)) with weights A=(

κ 4

− 1,

κ 2

+ 2,

κ 4

κ κ κ κ − 1, − , − , − , − ) 4

4

4

(41)

4

for κ > 0. The same calculation can be carried out for the second component v of the velocity variable; details are omitted. 4.2. A 5-point finite difference scheme For the particular value of κ = 4, the weight A in (41) for the 7-point stencil becomes to be A = (0, 4, 0, −1, −1, −1, −1)

(42)

so that (28) is reduced to a 5-point finite difference scheme described as follows: 5-Point Finite Difference Algorithm 4.1. Find ukℓ and pij such that (1) the discrete homogeneous Dirichlet boundary condition of ukℓ = 0 is satisfied at all the red dotted grid points (xk , yℓ ) on the domain boundary, and (2) the following set of finite difference equations are satisfied:

⎧ ⎤ ⎡ pi+1,j − pi,j ⎪ 4ui+ 1 ,j − ui+1,j− 1 − ui+1,j+ 1 − ui,j− 1 − ui,j+ 1 ⎪ ⎪ 2 2 2 2 2 ⎪ ⎦ = 1f 1 , h ⎪ +⎣ ⎪ ⎪ h2 2 i+ 2 ,j ⎪ ⎪ 0 ⎪ ⎨ ⎡ ⎤ − u − u − u 4u 1 1 1 1 1 − u 0 1 i− 2 ,j+1 i+ 2 ,j+1 i− 2 ,j i+ 2 ,j i,j+ 2 ⎪ ⎪ + ⎣ pi,j+1 − pi,j ⎦ = f i,j+ 1 , ⎪ ⎪ 2 ⎪ h2 2 ⎪ ⎪ ⎪ h ⎪ ⎪ ⎩ ui+ 1 ,j − ui− 1 ,j + vi,j+ 1 − vi,j− 1 = 0, 2

2

2

(43)

2

where

⎡ f k,ℓ =

(1)

⎤

⎣fk,ℓ ⎦ (2)

fk,ℓ

⎡ :=

⎤

(1) ⎣f (xk , yℓ )⎦

f (2) (xk , yℓ )

.

By using the component notation for the velocity u = (u, v )t , we may rewrite the finite difference equations (43) as follows:

⎧ ⎪ ⎪ ⎪ ⎪(4ui+ 12 ,j − ui+1,j− 12 − ui+1,j+ 12 − ui,j− 12 − ui,j+ 12 ) + h(pi+1,j − pi,j ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (4vi+ 1 ,j − vi+1,j− 1 − vi+1,j+ 1 − vi,j− 1 − vi,j+ 1 ) ⎪ ⎪ 2 2 2 2 2 ⎨ (4ui,j+ 1 − ui− 1 ,j+1 − ui+ 1 ,j+1 − ui− 1 ,j − ui+ 1 ,j ) ⎪ 2 2 2 2 2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (4vi,j+ 1 − vi− 1 ,j+1 − vi+ 1 ,j+1 − vi− 1 ,j − vi+ 1 ,j ) + h(pi,j+1 − pi,j ) ⎪ ⎪ 2 2 2 2 2 ⎪ ⎪ ⎪ ⎩ u −u +v −v i+ 12 ,j

i− 12 ,j

i,j+ 21

i,j− 21

= = = =

h2 2 h2 2 h2 2 h2 2

f f f f

(1) i+ 12 ,j (2) i+ 12 ,j (1) i,j+ 12 (2) i,j+ 12

, , (44)

, ,

= 0.

Fig. 5 shows the stencil of (43) for the velocity field as a five-point finite difference scheme with weights A˜ = (4, −1, −1, −1, −1); the left figure is for the scheme centered at (xi+ 1 ,j , yj ) and the right one is for (xi , yj+ 1 ). 2

2

188

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 5. Stencils for the velocity field in the finite difference scheme (43).

5. Stability For simplicity, we introduce two bilinear forms: a(u, v) := κ S(u, v) + (∇w u, ∇w v),

(45)

b(v, w ) := (∇w · v, w ),

(46)

where u, v ∈ [W0 (Th )] and w ∈ P0 (Th ). The SWG scheme (9) is a typical saddle-point problem which can be analyzed by using the well known theory developed by Babu˘ska [28] and Brezzi [29]. The core of the theory is to verify two properties: 2

(i) boundedness and a certain coercivity for the bilinear form a(·, ·), (ii) boundedness and inf–sup condition for the bilinear form b(·, ·), which are to be established in the rest of this section. The space [W0 (Th )]2 is a normed linear space equipped with the following triple-bar norm

|||w|||2 = (∇w w, ∇w w) + κ S(w, w) (47) [ ] w for w = ∈ [W0 (Th )]2 . To show that |||·||| is indeed a norm in [W0 (Th )]2 , it suffices to verify the length property for z

|||·|||. To this end, assume |||w||| = 0 for some w ∈ [W0 (Th )]2 . It follows that ∇w w = 0 ,

(w|ei −s(w)(Mi ))|ei | = 0

for each edge ei ∈ Eh . From (w|ei −s(w)(Mi ))|ei | = 0, we have

∇ s(w) =

N 1 ∑

|T |

s(w)(Mi )|ei | ⊗ ni =

i=1

N 1 ∑

|T |

w|ei |ei | ⊗ ni = ∇w w = 0,

i=1

where x ⊗ y = {xi yj }2×2 for any two vectors x = (x1 , x2 ) and y = (y1 , y2 ). Therefore, s(w) has constant value on each element T , and so does w on the boundary of the element. Moreover, since w = 0 on ∂ Ω , we then have w ≡ 0 on Ω as Ω is assumed to be connected. It is clear that a(w, w) = |||w|||2 for any w ∈ [W0 (Th )]2 . It follows from the definition of |||·||| and the usual Cauchy–Schwarz inequality that the following boundedness and coercivity hold true for the bilinear form a(·, ·). Lemma 4.

For any u, w ∈ [W0 (Th )]2 , we have

|a(u, w)| ≤ |||u||||||w|||, a(w, w) ≧ |||w|||2 .

(48) (49)

For the bilinear form b(·, ·), we have the following result on the inf–sup condition. Lemma 5. sup

There exists a positive constant β independent of h such that

v∈W0 (Th )

b(v, w )

|||v|||

for all w ∈ P0 (Th ).

≥ β∥w∥,

(50)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

189

Proof. For any given w ∈ P0 (Th ) ⊂ L20 (Ω ), it is well known [30–33] that there exists a vector-valued function v˜ ∈ [H01 (Ω )]2 such that (∇ · v˜ , w ) ≥ C1 ∥w∥2 ,

∥˜v∥1 ≤ C2 ∥w∥.

(51)

where Ci > 0 are two constants depending only on the domain Ω . Set v = Qb v˜ ∈ [W0 (Th )] where Qb is the L projection operator from [H01 (Ω )]2 to [W0 (Th )]2 . We define also Q0 and Qh the L2 projection operators onto [P1 (Th )]2 and [P0 (Th )]2×2 respectively. Then, for all q ∈ [P0 (Th )]2×2 , we have 2

2

(∇w Qb v˜ , q)T = ⟨Qb v˜ , q · n⟩∂ T

(52)

= ⟨˜v, q · n⟩∂ T = (∇ v˜ , q)T = (Qh (∇ v˜ ), q)T , which implies ∇w Qb v˜ = Qh (∇ v˜ ). Thus, the first part of the norm |||·||| verifies:

∑

∥∇w v∥2T =

T ∈ Th

∑

∥∇w Qb v˜ ∥2T =

T ∈ Th

∑

∥Qh ∇ v˜ ∥2T ≤ C ∥˜v∥21 .

(53)

T ∈Th

Next, we have

∥Qb (s(v) − Qb v)∥2∂ T = ⟨Qb (s(Qb v˜ )) − Qb v˜ , Qb (s(Qb v˜ )) − Qb v˜ ⟩∂ T = ⟨Qb (s(Qb v˜ )) − Qb v˜ , Qb (Q0 v˜ ) − Qb v˜ ⟩∂ T ≤ ∥Qb (s(Qb v˜ )) − Qb v˜ ∥∂ T ∥Qb (Q0 v˜ ) − Qb v˜ ∥∂ T , which gives

∥Qb (s(v) − Qb v)∥2∂ T ≤ ≤ ≤ ≤

∥Qb (Q0 v˜ ) − Qb v˜ ∥2∂ T ∥Q0 v˜ − v˜ ∥2∂ T C1 (h−1 ∥Q0 v˜ − v˜ ∥2T + h∥∇ (Q0 v˜ − v˜ )∥2T ) C2 h∥˜v∥21 ,

thus

κ S(v, v) = κ h−1

∑

∥Qb (s(v) − Qb v)∥2∂ T ≤ C3 ∥˜v∥1 .

(54)

T ∈ Th

Combining the two equations (53) and (54) we arrive at

|||v||| ≤ C0 ∥˜v∥1 ,

(55)

for some constant C0 . On the other side, from (51) we have b(v, w ) = (∇w · v, w ) = (∇w · Qb v˜ , w )

= = = ≥

(56)

⟨Qb v˜ · n, w⟩ ⟨˜v · n, w⟩ (∇ · v˜ , w ) C1 ∥w∥2 .

Using the above equation, (55), and (51), we obtain C1 ∥w∥2 |b(v, w)| ≥ ≥ β∥w∥, |||v||| C0 ∥˜v∥1 for a positive constant β . This completes the proof of the lemma.

(57)

It follows from Lemmas 4 and 5 that the following solvability holds true for the simplified weak Galerkin algorithm (9). Readers are referred to [1–3] for a detailed discussion on the original weak Galerkin finite element method. Theorem 6. The numerical scheme (9) has one and only one solution for any positive stabilization parameter κ > 0. Due to the connection between the finite difference scheme (28) and the SWG finite element method, we have the following result for the solvability of the finite difference method. Theorem 7. For any given stabilization parameter κ > 0, the finite difference scheme (28) has one and only one solution. The same conclusion holds true for the five-point finite difference scheme (43).

190

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

6. Error estimates

[ ] Let u and p be the exact solution of the model problem (1), and denote by ub =

ub

vb

∈ [W0 (Th )]2 and ph ∈ P0 (Th ) the

numerical approximation arising from the SWG scheme (9). Let Qb u and Qh p be the L2 projection of u and p in the spaces [W0 (Th )]2 and P0 (Th ), respectively. By the error function we mean the difference between the L2 projection and the SWG approximations: e = Qb u − ub , η = Qh p − ph .

(58)

We now derive two equations for which the error functions e and η must satisfy. The resulting equations are called the error equations, which play an important role in the convergence analysis of the SWG scheme. Lemma 8. Let (u; p) ∈ [H 1 (Ω )]2 × L2 (Ω ) be sufficiently smooth and satisfy the following equation

− ∆u + ∇ p = f in Ω .

(59)

The following equation holds true

∑

(∇w Qb u, ∇w v)T −

∑

T

(∇w · v, Qh p)T

(60)

T

= (f , s(v)) +

∑ ∂u ⟨ − Qh (∇ u) · n, s(v) − v⟩∂ T ∂n T

∑ − ⟨(s(v) − v) · n, p − Qh p⟩∂ T , T

for all v ∈ [W0 (Th )]2 , where Qh is the L2 projection operator onto [P0 (Th )]2 . Proof. From Eq. (52), we have

∑ =

T ∑

(∇w Qb u, ∇w v)T (Qh (∇ u), ∇w v)T

(61)

T

=

∑

⟨Qh (∇ u) · n, v⟩∂ T

T

=

∑ ⟨Qh (∇ u) · n, v⟩∂ T − ⟨Qh (∇ u) · n, s(v)⟩∂ T + ⟨Qh (∇ u) · n, s(v)⟩∂ T

=

T ∑ ⟨Qh (∇ u) · n, v − s(v)⟩∂ T + (∇ u, ∇ s(v))T T

∑ ∂u = ⟨Qh (∇ u) · n, v − s(v)⟩∂ T + (−∇ · (∇ u), s(v))T + ⟨ , s(v)⟩∂ T ∂n T ∑ ∂u = (−∆u, s(v)) + ⟨ − Qh (∇ u) · n, s(v) − v⟩∂ T , ∂n T

for all v ∈ [W0 (Th )] , where v|∂ Ω = 0 has been used. For (∇w · v, Qh p)T , we have 2

∑

(∇w · v, Qh p)T

T

= = = = =

∑

−(s(v), ∇ Qh p)T + ⟨v · n, Qh p⟩∂ T

T ∑

(∇ s(v), Qh p)T + ⟨(v − s(v)) · n, Qh p⟩∂ T

T ∑

(∇ s(v), p)T + ⟨(v − s(v)) · n, Qh p⟩∂ T

T ∑ T ∑

−(s(v), ∇ p)T + ⟨s(v) · n, p⟩∂ T + ⟨(v − s(v)) · n, Qh p⟩∂ T −(∇ p, s(v))T + ⟨(s(v) − v) · n, p − Qh p⟩∂ T

T

= − (∇ p, s(v)) +

∑ T

⟨(s(v) − v) · n, p − Qh p⟩∂ T

(62)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

191

for all v ∈ [W0 (Th )]2 , where we have ∑used the definition for the weak divergence operator (∇w ·) and the fact that ∇ Qh p = 0 in the second line, the identity of T ⟨v · n, p⟩∂ T = 0 in the sixth line. We now test (59) with s(v), v ∈ [W0 (Th )]2 , to obtain (−∆u, s(v)) + (∇ p, s(v)) = (f , s(v)).

(63)

Substituting Eqs. (61) and (62) into (63) yields

∑

∑

T

T

(∇w Qb u, ∇w v)T −

= (f , s(v)) +

(∇w · v, Qh p)T

∑ ∑ ∂u − Qh (∇ u) · n, s(v) − v⟩∂ T − ⟨(s(v) − v) · n, p − Qh p⟩∂ T , ⟨ ∂n T

T

which completes the proof of the lemma. The following is a result on the error equation for the SWG scheme (9). Lemma 9. Let e and η be the error functions for the numerical solution arising from the SWG scheme (9), as given in (58). Then, we have a(e, v) − b(v, η) = ϕu,p (v),

(64)

b(e, w ) = 0,

(65)

for all v ∈ [W0 (Th )] and w ∈ P0 (Th ), where 2

ϕu,p (v) = κ S(Qb u, v) +

∑ ∂u − Qh (∇ u) · n, s(v) − v⟩∂ T ⟨ ∂n T

∑ − ⟨(s(v) − v) · n, p − Qh p⟩∂ T .

(66)

T

Proof. By adding κ S(Qb u, v) to both sides of (60) we obtain a(Qb (u), v) − b(v, Qh p)

= κ S(Qb u, v) +

∑

∑

T

T

(∇w Qb u, ∇w v)T −

= (f , s(v)) + κ S(Qb u, v) +

(∇w · v, Qh p)T

∑ ∂u ⟨ − Qh (∇ u) · n, s(v) − v⟩∂ T ∂n

(67)

T

−

∑ ⟨(s(v) − v) · n, p − Qh p⟩∂ T . T

Now subtracting the first equation of (9) from (67) yields a(e, v) − b(v, η) = ϕu,p (v)

∀v ∈ [W0 (Th )]2 .

This completes the derivation of the error equation (64). As to the second one (65), we have b(e, w ) = (∇w · e, w ) = (∇w · (Qb u − ub ), w )

= = = =

(∇w · Qb u, w ) − (∇w · ub , w ) (Qh ∇ · u, w ) − 0 (Qh (0), w ) − 0 0,

for all w ∈ P0 (Th ). This completes the proof of the lemma. 6.1. Error estimates in H 1 The goal here is to establish an optimal-order error estimate for the numerical solution (ub ; ph ) in a discrete H 1 -norm for the velocity and the L2 norm for the pressure. The result can be stated as follows. Theorem 10. Let (u; p) ∈ [H01 (Ω ) ∩ H 2 (Ω )]2 × (L20 (Ω ) ∩ H 1 (Ω )) be the solution of (1), and (ub ; ph ) ∈ [W0 (Th )]2 × P0 (Th ) be the solution of the SWG scheme (9), respectively. Then, the following error estimate holds true

|||Qb u − ub ||| + ∥Qh p − ph ∥ ≤ Ch(∥u∥2 + ∥p∥1 ).

(68)

192

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Proof. By letting v = e in (64) and then using (65) with w = η, we have

|||e|||2 = ϕu,p (e) = κ S(Qb u, e) +

∑ ∂u ⟨ − Qh (∇ u) · n, s(e) − e⟩∂ T ∂n T

(69)

∑ − ⟨(s(e) − e) · n, p − Qh p⟩∂ T . T

From the Cauchy–Schwarz inequality and (54), for the first term in (69), we have

⏐ ⏐ ⏐ −1 ∑ ⏐ ⏐ |S(Qb u, e)| = ⏐h ⟨Qb s(Qb u) − Qb u, Qb s(e) − Qb e⟩∂ T ⏐⏐ T ⏐ ⏐ ⏐ −1 ∑ ⏐ ⏐ = ⏐h ⟨Qb (Q0 u) − Qb u, Qb s(e) − Qb e⟩∂ T ⏐⏐ T ⏐ ⏐ ⏐ −1 ∑ ⏐ = ⏐⏐h ⟨Qb (Q0 u − u), Qb s(e) − Qb e⟩∂ T ⏐⏐

(70)

T

)1 1 ≤ h ∥Q0 u − u∥ + ∥∇ (Q0 u − u)∥ 2 S(e, e) 2 ≤ Ch∥u∥2 |||e|||. (

−2

Next, for any constant tensor φ, we have from the definition of the weak gradient that (∇w · e, φ)T = ⟨e, φ · n⟩∂ T

= ⟨e − s(e), φ · n⟩∂ T + ⟨s(e), φ · n⟩∂ T = ⟨e − Qb s(e), φ · n⟩∂ T + (∇ s(e), φ)T , which gives

∥∇ s(e)∥2T ≤ ∥∇w e∥2T + Ch−1 ∥e − Qb s(e)∥2∂ T . It follows that

∥e − s(e)∥∂ T ≤ ∥e − Qb s(e)∥∂ T + ∥(I − Qb )s(e)∥∂ T ≤ ∥e − Qb s(e)∥∂ T + Ch1/2 ∥∇ s(e)∥T ≤ C ∥e − Qb s(e)∥∂ T + Ch

1/2

(71)

∥∇w e∥T .

From the Cauchy–Schwarz inequality and the estimate (71) we obtain

⏐ ⏐ ⏐ ∂u ⏐ ⏐⟨ ⏐ ⏐ ∂ n − Q0 (∇ u) · n, e − s(e)⟩∂ T ⏐ ∂u ≤∥ − Q0 (∇ u) · n∥∂ T ∥e − s(e)∥∂ T ∂n ( ) ∂u ≤ C∥ − Q0 (∇ u) · n∥∂ T ∥e − Qb s(e)∥∂ T + h1/2 ∥∇w e∥T . ∂n Now summing over all the elements yields an estimate for the second term in (69):

⏐ ∑ ⏐⏐ ∂ u ⏐ ⏐ ⏐⟨ ⏐ ∂ n − Q0 (∇ u) · n, e − s(e)⟩∂T ⏐ T ∑ ∂u ( ) ≤ ∥ − Q0 (∇ u) · n∥∂ T ∥e − Qb s(e)∥∂ T + Ch1/2 ∥∇w e∥T ∂n T ( )1 ( )1 ≤C ∥∇ u − Q0 (∇ u)∥20 + h2 ∥∇ 2 u∥0 2 ∥∇w e∥2 + κ S(e, e) 2 ≤Ch∥u∥2 |||e|||.

(72)

Similarly, we have

⏐ ⏐ ⏐∑ ⏐ ⏐ ⏐ ⏐ ⟨(s(e) − e) · n, p − Qh p⟩∂ T ⏐ ≤ Ch∥p∥1 |||e|||. ⏐ ⏐

(73)

T

Combining the estimates (70), (72), and (73) with (69) yields the first part of the error estimate (68):

|||Qb u − ub ||| ≤ Ch(∥u∥2 + ∥p∥1 ).

(74)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

193

To estimate ∥η∥, we use Eq. (64) to obtain b(v, η) = a(e, v) − ϕu,p (v). Using the equation above, (70)–(73), and (48), we arrive at

|b(v, η)| ≤ Ch(∥u∥2 + ∥p∥1 )|||v|||. Combining the above estimate with the inf–sup condition (50) gives

∥η∥ ≤ Ch(∥u∥2 + ∥p∥1 ), which, together with (74), yields the desired estimate (68). 6.2. Error estimates in L2 To derive an L2 -error estimate for the velocity approximation, we consider the problem of seeking (ψ; ξ ) such that

{ −∆ψ + ∇ξ = g divψ = 0 ψ =0

in Ω in Ω on Γ = ∂ Ω .

(75)

Assume that the problem (75) has the [H 2 (Ω )]2 × H 1 (Ω )-regularity in the sense that the solution (ψ; ξ ) ∈ [H 2 (Ω )]2 × H 1 (Ω ) and the following a priori estimate holds true:

∥ψ∥2 + ∥ξ ∥1 ≤ C ∥g ∥

(76)

for any g ∈ [L2 (Ω )]2 . Theorem 11. Let (u; p) ∈ [H01 (Ω ) ∩ H 2 (Ω )]2 × (L20 (Ω ) ∩ H 1 (Ω )) be the solution of (1), and (ub ; ph ) ∈ [W0 (Th )]2 × P0 (Th ) be the solution of the SWG scheme (9), respectively. Then, the following error estimate holds true:

∥s(Qb u) − s(ub )∥ ≤ Ch2 (∥u∥2 + ∥p∥1 ).

(77)

Proof. Let (ψ; ξ ) be the solution of (75) with g = s(e) = s(Qh u − ub ). From (60), we have for any v ∈ [W0 (Th )]2

∑

∑

T

T

(∇w Qb ψ, ∇w v)T −

=(s(e), s(v)) +

(∇w · v, Qh ξ )T

∑ ∂ψ ⟨ − Qh (∇ψ) · n, s(v) − v⟩∂ T ∂n

(78)

T

−

∑ ⟨(s(v) − v) · n, ξ − Qh ξ ⟩∂ T . T

In particular, by letting v = e we obtain

∑

(∇w Qb ψ, ∇w e)T −

T

∑

(∇w · e, Qh ξ )T

T

= (s(e), s(e)) +

∑ ∂ψ ⟨ − Qh (∇ψ) · n, s(e) − e⟩∂ T ∂n T

∑ − ⟨(s(e) − e) · n, ξ − Qh ξ ⟩∂ T . T

By adding and subtracting κ S(Qb ψ, e) in the above equation we arrive at

∥s(e)∥2 = a(Qb ψ, e) − b(e, Qh ξ ) ∑ ∂ψ − κ S(Qb ψ, e) − ⟨ − Qh (∇ψ) · n, s(e) − e⟩∂ T ∂n T ∑ + ⟨(s(e) − e) · n, ξ − Qh ξ ⟩∂ T

(79)

T

= a(Qb ψ, e) − b(e, Qh ξ ) − ϕψ,ξ (e). From (65) and the second equation of (75) we have b(e, Qh ξ ) = 0,

b(Qb ψ, η) = 0.

(80)

194

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Combining the above two equations gives

∥s(e)∥2 = a(e, Qb ψ) − b(Qb ψ, η) − ϕψ,ξ (e).

(81)

Using (64) and the equation above, we obtain

∥s(e)∥2 = ϕu,p (Qb ψ) − ϕψ,ξ (e).

(82)

The rest of the proof will establish some estimates for the right-hand side of (82). The following are two useful estimates in the forthcoming analysis. For any v ∈ [W (T )]2 , we have

∥v − Qb s(v)∥∂ T ≤ ∥v − Qb φ∥∂ T ∥s(v)∥∂ T ≤ C ∥v∥∂ T .

∀φ ∈ [P1 (T )]2 ,

(83) (84)

The inequality (83) is a direct consequence of the definition of s(v) as the projection of v into the space of linear functions with respect to the norm ∥ · ∥∂ T . The inequality (84) can be derived by turning ∥s(v)∥∂ T into a discrete norm using the midpoints Mi . Each of the terms in ϕu,p (Qb ψ ); namely,

ϕu,p (Qb ψ) = κ S(u, Qb ψ) +

∑ ∂u ⟨ − Qh (∇ u) · n, s(Qb ψ) − Qb ψ⟩∂ T ∂n T

∑ − ⟨(s(Qb ψ) − Qb ψ) · n, p − Qh p⟩∂ T , T

can be handled as follows:

• For the stability term κ S(u, Qb ψ), one has from the orthogonality property of the extension operator s(·) ⏐ ⏐ ∑ ⏐ ⏐ |S(u, Qb ψ)| = ⏐⏐h−1 ⟨Qb s(u) − Qb u, Qb s(Qb ψ) − Qb ψ⟩∂ T ⏐⏐ T ⏐ ⏐ ⏐ ⏐ −1 ∑ ⏐ = ⏐h ⟨Qb (Q0 u) − Qb u, Qb s(Qb ψ) − Qb ψ⟩∂ T ⏐⏐ T

≤ h−1

∑

∥Q0 u − u∥∂ T ∥Qb s(Qb ψ) − Qb ψ∥∂ T .

T

From (83) with v = Qb ψ and φ = Q0 ψ we have

|S(u, Qb ψ)| ≤ h−1

∑

∥Q0 u − u∥∂ T ∥Qb s(Qb ψ) − Qb ψ∥∂ T

T

≤ h−1

∑

∥Q0 u − u∥∂ T ∥Qb (Q0 ψ) − Qb ψ∥∂ T

T

) 21 (

( ≤ h−1

∑

) 12 ∑

∥Q0 u − u∥2∂ T

T

∥Q0 ψ − ψ∥2∂ T

T

≤ Ch2 ∥u∥2 ∥ψ∥2 . • For the second term in ϕu,p (Qb ψ), since ψ = 0 on ∂ Ω , we have ∑ ∂u ∑ ∂u ⟨ − Qh (∇ u) · n, Qb ψ − ψ⟩∂ T = ⟨ , Qb ψ − ψ⟩∂ T = 0, ∂n ∂n T

T

which leads to

⏐ ⏐∑ ⏐ ⏐ ∂u ⏐ ⟨ − Qh (∇ u) · n, s(Qb ψ) − Qb ψ⟩∂ T ⏐⏐ ⏐ ∂n T ⏐∑ ⏐ ⏐ ⏐ ∂u ⏐ =⏐ ⟨ − Qh (∇ u) · n, s(Qb ψ) − ψ⟩∂ T ⏐⏐ ∂n T ⏐ ⏐∑ ⏐ ⏐ ∂u = ⏐⏐ ⟨ − Qh (∇ u) · n, s(Qb ψ) − s(Qb (Q0 ψ)) + Q0 ψ − ψ⟩∂ T ⏐⏐, ∂n T

where we have used the fact that

s(Qb (Q0 ψ )) = Q0 ψ.

(85)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

195

It follows from (84) that

⏐∑ ⏐ ⏐ ⏐ ∂u ⏐ ⏐ − Q ( ∇ u) · n , s (Q ψ ) − Q ψ⟩ ⟨ h b b ∂T ⏐ ⏐ ∂n T

) 12

( ∑

≤

2

∥∇ u · n − Qh (∇ u) · n∥∂ T

T

) 12

( ∑

·

(86)

,

∥Q0 ψ − ψ∥2∂ T + ∥s(Qb (ψ − Q0 ψ))∥2∂ T

T

) 12 (

( ∑

≤C

∥∇ u · n − Qh (∇ u) · n∥2∂ T

) 12 ∑

T

∥Q0 ψ − ψ∥2∂ T

T

≤Ch2 ∥u∥2 ∥ψ∥2 . • As to the last term in ϕu,p (Qb ψ), we have ⏐ ⏐∑ ⏐ ⏐ ⏐ ⏐ ⟨ ( s (Q ψ ) − Q ψ ) · n , p − Q p ⟩ b b h ∂ T ⏐ ⏐ T ⏐ ⏐∑ ⏐ ⏐ =⏐⏐ ⟨(s(Qb ψ) − ψ) · n, p − Qh p⟩∂ T ⏐⏐ T ⏐∑ ⏐ ⏐ ⏐ ⏐ =⏐ ⟨(s(Qb ψ) − s(Qb Q0 ψ)) · n + (Q0 ψ − ψ) · n, p − Qh p⟩∂ T ⏐⏐

(87)

T

) 12 (

( ≤C

∑

2

h∥p − Qh p∥∂ T

) 12 ∑

T

h

−1

∥Q0 ψ − ψ∥∂ T 2

T

≤Ch2 ∥p∥1 ∥ψ∥2 . The estimates (85), (86) and (87), together with the regularity (76), collectively yield

|ϕu,p (Qb ψ)| ≤ Ch2 (∥u∥2 + ∥p∥1 )∥ψ∥2 ≤ Ch2 (∥u∥2 + ∥p∥1 )∥s(e)∥.

(88)

Next, using the same argument as in the proof of Theorem 10 and the H 2 -regularity assumption (76) we arrive at

|ϕψ,ξ (e)| ≤ Ch(∥ψ∥2 + ∥ξ ∥1 )|||e||| ≤ Ch∥s(e)∥|||e|||. It then follows from the estimate (68) that

|ϕψ,ξ (e)| ≤ Ch2 (∥u∥2 + ∥p∥1 )∥s(e)∥.

(89)

Finally, substituting (88) and (89) into (82) gives

∥s(e)∥2 ≤ Ch2 (∥u∥2 + ∥p∥1 )∥s(e)∥,

(90)

which gives rise to the desired error estimate (77). 7. Superconvergence Next we shall derive a superconvergence for the approximate velocity and pressure when the finite element partition consists of only rectangular elements. The result is based on the framework developed in [6] for the second order elliptic equation, with a particular attention paid to the pressure term. Theorem 12. Let (u; p) ∈ [H01 (Ω ) ∩ H 3 (Ω )]2 × (L20 (Ω ) ∩ H 2 (Ω )) be the solution of (1), and (ub ; ph ) ∈ [W0 (Th )]2 × P0 (Th ) be the solution of the SWG scheme (9) on rectangular meshes, respectively. Then, the following superconvergence holds true:

∥∇w (Qh u − ub )∥0 + ∥Qh p − ph ∥0 ≤ Ch3/2 (∥u∥3 + ∥p∥2 ).

(91)

Proof. The proof ∑ relies on a refined treatment of the linear functional ϕu,p (·) in the error equation (64). To this end, 2 consider the term T ⟨(s(v) − v) · n, p − Qh p⟩∂ T on the right-hand side of (66), where v ∈ [W0 (Th )] is an arbitrary test function. Recall that, from (23) and (24), s(v) − v has the same value at the center of e1 and e2 (respectively, e3 and e4 ).

196

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

As shown in Fig. 2, we have n1 = −n2 (respectively, n3 = −n4 for the two horizontal edges). As Qh p is constant-valued on T , we thus have

⟨(s(v) − v) · n, Qh p⟩∂ T = 0. Furthermore, with v = (v1 , v2 )t , we obtain

⟨(s(v) − v) · n, p − Qh p⟩∂ T = ⟨(s(v) − v) · n, p⟩∂ T = ⟨s(v1 ) − v1 , p⟩e2 − ⟨s(v1 ) − v1 , p⟩e1 + ⟨s(v2 ) − v2 , p⟩e4 − ⟨s(v2 ) − v2 , p⟩e3 = (χ1 , ∂x p)T + (χ2 , ∂y p)T ,

(92)

where χ1 is the constant extension of s(v1 ) − v1 along the x-direction; i.e.,

χ1 (x, y) := (s(v1 ) − v1 )|e1 (y) = (s(v1 ) − v1 )|e2 (y),

(x, y) ∈ T

and χ2 is defined analogously as the constant extension of s(v2 ) − v2 along the y-direction:

χ2 (x, y) := (s(v2 ) − v2 )|e3 (x) = (s(v2 ) − v2 )|e4 (x),

(x, y) ∈ T .

Eqs. (27) show that the above extensions are well-defined. The two terms on the right-hand side of (92) can be handled as follows. As χ1 is linear in y and constant in x, we have (χ1 , ∂x p)T = χ1 (M1 )

∫

∂x pdT + R11 (T ),

(93)

T

where Mi is the center of the edge ei and R11 (T ) =

∫ T

2 E1 (y)∂y χ1 ∂xy pdT with the kernel satisfying |E1 (y)| ≤ Ch2 . It is easy

to see that ∂y χ1 = ∂y s(v1 ). Thus, the error term R11 (T ) has the following estimate:

|R11 (T )| ≤ Ch2 ∥∇ s(v1 )∥T ∥∇ 2 p∥T .

(94)

For the first term on the right-hand of (93), we apply the trapezoidal rule in the x-direction to obtain

χ1 (M1 )

∫

∂x pdT = T

h 2

χ1 (M1 )

∫

h

∂x pdy + χ1 (M2 ) 2

e1

∫

∂x pdy + R12 (T ),

(95)

e2

where we have used χ1 (M1 ) = χ1 (M2 ) and R12 (T ) = χ1 (M1 )

∫

2 E2 (x)∂xx pdT T

with the kernel function satisfying |E2 (x)| ≤ Ch. The remainder term R12 (T ) can then be bounded as follows:

|R12 (T )| ≤ Ch2 |χ1 (M1 )| ∥∇ 2 p∥T .

(96)

Substituting (95) into (93) yields the following: (χ1 , ∂x p)T =

h 2

χ1 (M1 )

∫

h

∂x pdy + χ1 (M2 )

∫

2

e1

∂x pdy + R11 (T ) + R12 (T ).

(97)

e2

By introducing the following function

ρb(1)

⎧ ∫ 1 ⎪ ⎪ ∂x pdy ⎪ ⎪ ⎪ ⎨ |e1 | ∫e1 1 = ∂x pdy ⎪ ⎪ |e2 | e2 ⎪ ⎪ ⎪ ⎩

0

on e1 , on e2 , on e3 and e4 ,

we may rewrite (97) as follows: (χ1 , ∂x p)T =

=

h 2 h 2

⟨χ1 , ρb(1) ⟩∂ T + R11 (T ) + R12 (T ) (98)

⟨s(v1 ) − v1 , ρb(1) ⟩∂ T + R11 (T ) + R12 (T ).

Analogously, for the second term on the right-hand side of (92), there are remainder terms R21 (T ) and R22 (T ) such that (χ2 , ∂y p)T =

h 2

⟨s(v2 ) − v2 , ρb(2) ⟩∂ T + R21 (T ) + R22 (T ),

(99)

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

197

where

ρb(2) =

⎧ ∫ ⎪ 1 ⎪ ⎪ ∂y pdx ⎪ ⎪ ⎪ ⎨ |e3 | ∫e3 1

on e3 ,

∂y pdx

⎪ ⎪ |e4 | ⎪ ⎪ ⎪ ⎪ ⎩

on e4 ,

e4

on e1 and e2 ,

0

By setting ρb = (ρ

(1) b

(2) t b ) ,

,ρ

we have from combining (92) with (98) and (99) that

⟨(s(v) − v) · n, p − Qh p⟩∂ T = where R(T ) =

∑2

i,j=1

∥∇ s(v)∥T +

|R(T )| ≤ Ch

⟨s(v) − v, ρb ⟩∂ T + R(T ),

(100)

Rij (T ) is the remainder satisfying

( 2

h 2

4 ∑

) |(s(v) − v)(Mi )| ∥∇ 2 p∥T

(101)

i=1

( ) ≤ Ch ∥∇w v∥T + ST (v, v)1/2 ∥∇ 2 p∥T . 2

Here we have used the estimate (26) and the expression (4) in the second inequality. Summing over all the elements in Th gives

∑

⟨(s(v) − v) · n, p − Qh p⟩∂ T =

T ∈Th

h ∑ 2

⟨s(v) − v, ρb ⟩∂ T + R,

(102)

T ∈Th

where

⏐ ⏐ ⏐ ⏐∑ ⏐ ⏐ ( ) R(T )⏐⏐ ≤ Ch2 ∥∇w v∥ + S(v, v)1/2 ∥∇ 2 p∥. |R| = ⏐⏐ ⏐ ⏐T ∈Th

(103)

As to the first two terms on the right-hand side of (66), it has been shown in [6] that there exist another function

σb = (σb(1) , σb(2) )t and a remainder R5 such that κ S(Qb u, v) +

∑ ∂u ⟨ − Qh (∇ u) · n, s(v) − v⟩∂ T ∂n T

=

h ∑ 2

(104)

⟨v − s(v), σb ⟩∂ T + R5 ,

T ∈Th

where

( ) |R5 | ≤ Ch2 ∥∇w v∥ + S(v, v)1/2 ∥∇ 3 u∥.

(105)

Combining (66) with (104) and (102) gives

ϕu,p (v) = = =

h ∑ 2

⟨v − s(v), σb + ρb ⟩∂ T + R5 − R

T ∈Th

h2 2 h2 2κ

(106)

S(σb + ρb , v) + R5 − R [a(σb + ρb , v) − (∇w (σb + ρb ), ∇w v)] + R5 − R, 2

where we have used (45) in the last equation. By letting e˜ = e − 2hκ (σb +ρb ), the error equations (64)–(65) can be rewritten as a(e˜ , v) − b(v, η) = R5 − R − b(e˜ , w ) = −

h

2

2κ

h2 2κ

(∇w (σb + ρb ), ∇w v),

b(σb + ρb , w ),

(107) (108)

198

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

where the right-hand sides have the following bound:

⏐ ⏐ 2 ⏐ ⏐ ⏐R5 − R − h (∇w (σb + ρb ), ∇w v)⏐ ≤ Ch2 (∥u∥3 + ∥p∥2 )a(v, v)1/2 , ⏐ ⏐ 2κ ⏐ ⏐ 2 ⏐ ⏐h ⏐ b(σb + ρb , w)⏐ ≤ Ch2 (∥u∥3 + ∥p∥2 )∥w∥. ⏐ ⏐ 2κ

(109)

Eqs. (107)–(108) with O(h2 ) load functions on the right-hand sides lead immediately to the superconvergence estimate (91) by selecting a particular test function v. The estimate for η is given by the inf–sup condition. The drop of half order 2 is due to the fact that the small perturbation term 2hκ (σb + ρb ) applied to the error function e may not vanish on the boundary. Details on how the technique works can be found in [6]. 8. Numerical experiments In this section, we numerically verify the theoretical error estimates developed in the previous sections for the SWG scheme (9) and the finite difference schemes (28) and (43) when applied to the Stokes problem (1). The following metrics are used to measure the magnitude of the error:

• For the SWG scheme: L2 -norm for the velocity w = (u, v )T :

∥wb − w∥20 =

∑

h2 |wb,e − w(xe , ye )|2 ,

e∈Eh

H -norm for the velocity w = (u, v )T : 1

∥wb − w∥21 =

∑

⏐2

h2 ⏐∇w wb,T − ∇w(xc , yc )⏐ ,

⏐

T ∈Th 2

L -norm for the pressure p :

∥ph − p∥20 =

n n ∑ ∑

h2 |pi,j − p(xc , yc )|2 ,

i=1 j=1

where (xe , ye ) is the midpoints of the element edge e and (xc , yc ) is the barycenter of element T . • For the finite difference schemes: L2 -norm for the velocity w = u, v :

∥wb − w∥20 =

n+1 n ∑ ∑

h2 |wi− 1 ,j − w (xi− 1 , yj )|2 + 2

n n+1 ∑ ∑

2

h2 |wi,j− 1 − w (xi , yj− 1 )|2 , 2

2

i=1 j=1

i=1 j=1

H 1 -norm for the velocity w = u, v :

⏐2 ⏐w 1 − w 1 ⏐ ⏐ i+ 2 ,j ∂w i − 2 ,j ⏐ − (xi , yj )⏐⏐ ∥wb − w∥ = h ⏐ h ∂x i=1 j=1 ⏐ ⏐2 n n ∑ ∑ ⏐ wi,j+ 1 − wi,j− 1 ⏐ ∂w 2 2 + h2 ⏐⏐ − (xi , yj )⏐⏐ , h ∂y 2 1

n n ∑ ∑

2

i=1 j=1

2

L -norm for the pressure p:

∥ph − p∥20 =

n n ∑ ∑

h2 |pi,j − p(xi , yj )|2 ,

i=1 j=1

8.1. Test Case 1 The computational domain in the first test case is given by Ω = (0, π ) × (0, π ), which is partitioned into n × n small squares of equal-size. The exact solution of the model problem is given by

) ( ⎧ sin2 (x) cos(y) sin(y) ⎪ ⎨u = , − cos(x) sin(x) sin2 (y) ⎪ ⎩ p = cos(x) cos(y). The right-hand side function and the Dirichlet boundary data are chosen to match the exact solution. Note that this example has been considered in [5].

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

199

Fig. 6. Numerical results of test case 1 obtained with the FD scheme (28), κ = 4 using uniform rectangular meshes with h = 1/32. Table 1 Error and convergence performance of the finite difference scheme (28) for the Stokes equations on uniform meshes of test case 1. r refers to the order of convergence in O(hr ).

κ = 0.4 n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

2.78e−01 1.01e−01 3.07e−02 8.23e−03

– 1.47 1.71 1.90

4.34e−01 1.84e−01 5.95e−02 1.63e−02

– 1.23 1.63 1.87

2.52e−01 8.87e−02 2.66e−02 7.06e−03

– 1.51 1.74 1.91

2.85e−01 1.21e−01 3.95e−02 1.08e−02

– 1.23 1.61 1.86

5.72e−01 2.96e−01 1.03e−01 2.92e−02

– 0.95 1.52 1.82

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

1.35e−01 4.26e−02 1.17e−02 3.02e−03

– 1.66 1.86 1.96

2.38e−01 8.36e−02 2.38e−02 6.18e−03

– 1.51 1.82 1.94

1.36e−01 4.24e−02 1.15e−02 2.95e−03

– 1.69 1.88 1.96

1.62e−01 5.83e−02 1.67e−02 4.37e−03

– 1.47 1.80 1.94

3.85e−01 1.50e−01 4.45e−02 1.18e−02

– 1.36 1.76 1.91

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

2.35e−02 6.26e−03 1.60e−03 4.01e−04

– 1.91 1.97 1.99

5.90e−02 1.64e−02 4.25e−03 1.08e−03

– 1.85 1.95 1.98

5.69e−02 1.53e−02 3.89e−03 9.78e−04

– 1.90 1.97 1.99

6.61e−02 1.92e−02 5.01e−03 1.27e−03

– 1.78 1.94 1.98

1.48e−01 4.29e−02 1.13e−02 2.88e−03

– 1.79 1.92 1.97

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

3.29e−02 8.31e−03 2.08e−03 5.21e−04

– 1.98 2.00 2.00

5.11e−02 1.28e−02 3.20e−03 8.00e−04

– 2.00 2.00 2.00

3.92e−02 9.83e−03 2.46e−03 6.15e−04

– 1.99 2.00 2.00

5.85e−02 1.47e−02 3.68e−03 9.21e−04

– 1.99 2.00 2.00

8.32e−02 2.15e−02 5.43e−03 1.36e−03

– 1.95 1.99 2.00

κ = 1.0

κ=4

κ = 40

8.1.1. Numerical results obtained with the FD schemes The numerical approximation pressure field obtained with the FD scheme (28) is shown in Fig. 6(a), while the velocity vector and magnitude fields are plotted in Fig. 6(b). Table 1 illustrates the numerical performance of the finite difference scheme (28) for different values of κ = 0.4, 1.0, 4.0, 10.0. We note that, when κ = 4, the FD scheme (28) becomes (43). The numerical velocity is clearly converging at the optimal order of r = 2 in the L2 norm. For the numerical pressure, the error at the cell center is decreasing at the order of r = 2, and the same conclusion can be obtained for the numerical velocity in a discrete H 1 norm defined by using the cell centers. This implies that the pressure and the numerical velocity are of superconvergent to the exact solution at the center of cells. The numerical results outperform the theoretical prediction of r = 1.5 shown in Theorem 12. It is noted that the convergence rate may decrease for coarse meshes when the value of κ is small. But the optimal rate of convergence seems to be regained for fine meshes. 8.1.2. Numerical results of the SWG schemes on polygonal meshes The numerical pressure obtained from the SWG scheme (9) for test case 1 on hexagonal mesh is shown in Fig. 7(a), while the velocity magnitude is plotted in Fig. 7(b).

200

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 7. Numerical results of test case 1 obtained with the SWG scheme (9) using hexagonal meshes with h = 1/32, κ = 4. Table 2 Error and convergence performance of the SWG scheme (9) with κ = 4, for the Stokes equation on polygonal meshes with hexagonal elements of test case 1. r refers to the order of convergence in O(hr ). SWG for Test case 1 on Hexagonal mesh n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32 64

3.23e−01 1.15e−01 3.35e−02 9.33e−03 2.45e−03 6.17e−04

– 1.49 1.78 1.84 1.93 1.99

2.70e−01 1.19e−01 3.91e−02 1.18e−02 3.48e−03 1.06e−03

– 1.19 1.60 1.73 1.76 1.71

3.89e−01 1.79e−01 6.63e−02 1.86e−02 4.88e−03 1.25e−03

– 1.12 1.43 1.83 1.93 1.96

Table 2 illustrates the numerical performance of the SWG scheme (28). The numerical velocity is clearly converging at the optimal order of r = 2 in the L2 norm. For the numerical pressure, the error at the cell center is decreasing at the order of r = 1.9. For the numerical velocity, the error is decreasing at the order of r = 1.7 in the discrete H 1 norm defined by using the cell centers. The results indicate that the pressure and the velocity are superconvergent to the exact solution at the center of cells in L2 and H 1 norms, respectively. Moreover, numerical results concerning triangular, quadrilateral, and octagonal partitions have been included and remarked in the Appendix. 8.2. Test Case 2 In this test case, the analytical solution is given by

( ) ⎧ −256x2 (x − 1)2 y(y − 1)(2y − 1) ⎪ ⎨u = , 2 2 256y (y − 1) x(x − 1)(2x − 1) ⎪ ⎩ p = 150(x − 0.5)(y − 0.5), and the domain is the unit square Ω = (0, 1) × (0, 1). The right-hand side function and the Dirichlet boundary data are chosen to match the exact solution. Note that this example has been considered in [34]. 8.2.1. Numerical results obtained from the FD schemes The numerical approximation pressure field obtained from the FD scheme (28) is shown in Fig. 8(a), while the velocity vector and magnitude fields are plotted in Fig. 8(b). Table 3 illustrates the numerical performance of the finite difference scheme (28) for Test Case 2 with varying values of κ = 0.4, 1.0, 4.0, 10.0 included. It can be seen that the numerical velocity is converging at the optimal order of r = 2 in the usual L2 norm. For the pressure approximation, the error at the cell centers seems to be decreasing at the rate of O(h1.7 ) or better. It appears that the numerical velocity is also converging at the order of r = 1.7 in the discrete H 1 norm defined on the cell centers. The results clearly indicate a superconvergence of the finite difference scheme for the velocity in H 1 and the pressure in L2 on cell centers. Once again, the numerical results outperform the theoretical prediction of r = 1.5 shown in Theorem 12.

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

201

Fig. 8. Numerical results of test case 2 obtained with the SWG scheme using uniform rectangular meshes with h = 1/32. Table 3 Error and convergence performance of the finite difference scheme (28) for the Stokes equations on uniform meshes of test case 2. r refers to the order of convergence in O(hr ).

κ = 0.4 n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

4.21e−01 1.70e−01 5.63e−02 1.57e−02

– 1.31 1.60 1.84

1.83e+00 9.60e−01 3.63e−01 1.15e−01

– 0.93 1.40 1.66

5.19e−01 2.15e−01 6.90e−02 1.89e−02

– 1.27 1.64 1.87

2.51e+00 1.20e+00 4.21e−01 1.28e−01

– 1.07 1.51 1.72

4.47e+00 2.49e+00 9.67e−01 3.14e−01

– 0.84 1.36 1.62

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

2.39e−01 8.57e−02 2.49e−02 6.57e−03

– 1.48 1.78 1.92

1.26e+00 5.31e−01 1.75e−01 5.22e−02

– 1.25 1.60 1.74

2.79e−01 9.82e−02 2.80e−02 7.31e−03

– 1.51 1.81 1.94

1.49e+00 5.87e−01 1.87e−01 5.48e−02

– 1.35 1.65 1.77

3.20e+00 1.37e+00 4.66e−01 1.42e−01

– 1.22 1.56 1.71

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

1.03e−01 2.90e−02 7.55e−03 1.91e−03

– 1.82 1.94 1.98

6.26e−01 1.97e−01 5.73e−02 1.60e−02

– 1.67 1.78 1.84

7.17e−02 2.07e−02 5.43e−03 1.38e−03

– 1.79 1.93 1.98

4.78e−01 1.60e−01 4.88e−02 1.41e−02

– 1.57 1.72 1.79

1.39e+00 4.68e−01 1.43e−01 4.14e−02

– 1.58 1.71 1.79

n

∥uh − u∥0

r =

∥uh − u∥1

r =

∥vh − v∥0

r =

∥vh − v∥1

r =

∥ph − p∥0

r =

8 16 32 64

5.47e−02 1.37e−02 3.42e−03 8.54e−04

– 2.00 2.00 2.00

3.11e−01 7.97e−02 2.02e−02 5.11e−03

– 1.96 1.98 1.98

4.14e−02 1.04e−02 2.60e−03 6.51e−04

– 1.99 2.00 2.00

2.41e−01 6.24e−02 1.59e−02 4.05e−03

– 1.95 1.97 1.98

2.37e−01 6.87e−02 1.92e−02 5.25e−03

– 1.79 1.84 1.87

κ = 1.0

κ=4

κ = 40

8.2.2. Numerical results of the SWG schemes on polygonal meshes The numerical pressure obtained from the SWG scheme (9) on hexagonal meshes for test case 2 is shown in Fig. 9(a), while the velocity magnitude function is plotted in Fig. 9(b). Table 4 illustrates the numerical performance of the SWG scheme (9). The numerical velocity is clearly converging at the optimal order of r = 2 in the L2 norm. For the numerical pressure, the error at the cell center is decreasing at the order of r ≈ 1.7. For the velocity approximation, the convergence is seen to be at the order of r ≈ 1.6 in a discrete H 1 norm defined on the set of cell centers. The results suggest that the pressure and the velocity are superconvergent to the exact solution at the cell centers. More numerical results on triangular, quadrilateral, and octagonal meshes are reported in the Appendix. 8.3. Test Case 3: Lid driven Cavity The lid-driven cavity problem is a benchmark test case for Stokes flow, which has been tested in many existing literatures, including [5,9,35,36]. In this test case, a uniform mesh with meshsize h = 1/32 is employed in our new finite difference scheme. The source term is f = 0, and the Dirichlet boundary condition is given as u = (1, 0)T for

202

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. 9. Numerical results of test case 2 obtained with the SWG scheme (9) using hexagonal meshes with h = 1/32, κ = 4.

Fig. 10. Numerical results of lid-driven problem obtained with the FD scheme (28) using uniform rectangular meshes with h = 1/32.

Table 4 Error and convergence performance of the SWG scheme (9) with κ = 4, for the Stokes equation on polygonal meshes with hexagonal elements of test case 2. r refers to the order of convergence in O(hr ). SWG for Test case 2 on Hexagonal mesh n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32 64

1.18e+00 5.31e−01 1.79e−01 5.22e−02 1.40e−02 3.58e−03

– 1.15 1.56 1.78 1.90 1.97

4.85e+00 2.61e+00 1.08e+00 3.73e−01 1.22e−01 4.01e−02

– 0.89 1.27 1.53 1.61 1.61

5.50e+00 4.27e+00 2.08e+00 6.05e−01 1.78e−01 5.41e−02

– 0.37 1.04 1.78 1.77 1.72

y = 1, x ∈ (0, 1) and u = (0, 0)T on the rest of the boundary. The exact solution of lid-driven problem is not known, but the solution is known to have singularity at the corner points (0, 1) and (1, 1). The velocity field and streamlines obtained with the FD scheme (28) are plotted in Figs. 10(a) and 10(b), respectively. The pressure contour plot is shown in Fig. 11. The shape of streamlines is similar to the result given in [5,9,35,36], and the results look quite reasonable. The numerical results obtained with the SWG scheme (9) on hexagonal mesh are plotted in Fig. 12. The numerical results are similar to the results obtained by the FD scheme (28), and the results look quite reasonable.

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

203

Fig. 11. Numerical results — Pressure contours of lid-driven problem obtained with the FD scheme (28) using uniform rectangular meshes with h = 1/32.

Fig. 12. Numerical results of lid cavity test case 3 obtained with the SWG scheme (9) using hexagonal meshes with h = 1/32, κ = 4.

Appendix. Numerical results of the SWG scheme on polygonal meshes

A.1. Test case 1 The numerical pressures obtained from the SWG scheme (9) on polygonal meshes for the test case 1 are shown in Figs. A.13(a)–A.13(d), while the velocity magnitude function is plotted in Figs. A.13(e)–A.13(h). Table A.5 illustrates the numerical performance of the SWG scheme (9). The numerical velocity is converging at the optimal order of r = 2 in the L2 norm for each finite element partitions. For the pressure approximations, the error at the cell centers is decreasing at the order of r ≈ 1.3 for triangular meshes and r ≈ 1.7 for other type of meshes. The optimal

204

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

Fig. A.13. Numerical results of test case 1 obtained from the SWG scheme (9), numerical pressure: (a), (b), (c), (d), velocity magnitude: (e), (f), (g), (h), using polygonal meshes with h = 1/16.

Fig. A.14. Numerical results of test case 2 obtained from the SWG scheme (9), pressure field: (a), (b), (c), (d), velocity magnitude: (e), (f), (g), (h), using polygonal meshes with h = 1/16.

order of r = 1 can be observed for the numerical velocity in the discrete H 1 norm for triangular meshes and r ≈ 1.5 for other type of meshes. This implies that the pressure and the velocity are superconvergent to the exact solution at the cell centers for rectangular and hexagonal meshes.

A.2. Test case 2 The numerical pressure obtained from the SWG scheme (9) on various polygonal partitions for test case 2 is shown in Figs. A.14(a)–A.14(d), while the velocity magnitude is plotted in Figs. A.14(e)–A.14(h). Table A.6 illustrates the numerical performance of the SWG scheme (9). The numerical velocity is converging at the optimal order of r = 2 in the L2 norm for each tested polygonal partition. For the pressure approximation, the error at the cell centers is decreasing at the order of r ≈ 1.5 for triangular meshes and r = 1.7 for all other tested partitions. The optimal order of r = 1 can be observed for the numerical velocity in the discrete H 1 norm for all the polygonal partitions, while a superconvergence of r ≈ 1.6 can be observed for some polygonal meshes. In particular, the results suggest that the pressure and the velocity approximations are superconvergent to the exact solution at the cell centers for rectangular and hexagonal meshes.

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

205

Table A.5 Error and convergence performance of the SWG scheme (9) for the Stokes equation on polygonal meshes of test case 1. r refers to the order of convergence in O(hr ). Triangle

Square

n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32

7.39e−01 1.31e−01 3.75e−02 1.02e−02 2.61e−03

0.00 2.50 1.80 1.88 1.96

1.04e+00 3.41e−01 1.74e−01 8.88e−02 4.47e−02

0.00 1.61 0.97 0.97 0.99

4.00e−01 2.39e−01 1.03e−01 3.88e−02 1.59e−02

0.00 0.75 1.21 1.41 1.29

6.44e−01 1.40e−01 3.94e−02 1.03e−02 2.63e−03

0.00 2.20 1.83 1.93 1.98

2.18e−01 1.51e−01 5.28e−02 1.49e−02 3.87e−03

0.00 0.53 1.51 1.83 1.94

3.30e−01 2.04e−01 7.78e−02 2.29e−02 6.06e−03

0.00 0.69 1.39 1.77 1.92

Hexagonal

Polygonal

n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32

3.23e−01 1.15e−01 3.35e−02 9.33e−03 2.45e−03

– 1.49 1.78 1.84 1.93

2.70e−01 1.19e−01 3.91e−02 1.18e−02 3.48e−03

– 1.19 1.60 1.73 1.76

3.89e−01 1.79e−01 6.63e−02 1.86e−02 4.88e−03

– 1.12 1.43 1.83 1.93

7.11e−01 1.63e−01 4.65e−02 1.26e−02 3.22e−03

– 2.12 1.81 1.89 1.97

1.94e−01 1.33e−01 5.22e−02 1.74e−02 6.56e−03

– 0.55 1.34 1.59 1.41

4.41e−01 2.31e−01 7.93e−02 2.37e−02 6.56e−03

– 0.93 1.54 1.74 1.85

Table A.6 Error and convergence performance of the SWG scheme (9) for the Stokes equation on polygonal meshes of test case 2. r refers to the order of convergence in O(hr ). Triangle

Square

n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32

1.58e+00 5.56e−01 2.09e−01 6.28e−02 1.68e−02

– 1.51 1.41 1.73 1.90

8.27e+00 5.14e+00 3.14e+00 1.72e+00 8.87e−01

– 0.68 0.71 0.87 0.95

7.05e+00 5.22e+00 2.48e+00 9.24e−01 3.25e−01

– 0.43 1.07 1.43 1.51

1.46e+00 5.24e−01 1.85e−01 5.24e−02 1.37e−02

– 1.48 1.50 1.82 1.94

2.52e+00 2.81e+00 1.14e+00 3.68e−01 1.08e−01

–

−0.16 1.30 1.64 1.76

6.71e+00 4.77e+00 2.01e+00 6.73e−01 2.05e−01

– 0.49 1.25 1.58 1.72

Hexagonal

Polygonal

n

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

∥e(u,v) ∥0

r =

∥e(u,v) ∥1

r =

∥ep ∥0

r =

2 4 8 16 32

1.18e+00 5.31e−01 1.79e−01 5.22e−02 1.40e−02

– 1.15 1.56 1.78 1.90

4.85e+00 2.61e+00 1.08e+00 3.73e−01 1.22e−01

– 0.89 1.27 1.53 1.61

5.50e+00 4.27e+00 2.08e+00 6.05e−01 1.78e−01

– 0.37 1.04 1.78 1.77

1.85e+00 8.50e−01 2.92e−01 8.54e−02 2.27e−02

– 1.13 1.54 1.77 1.91

3.09e+00 2.84e+00 1.32e+00 5.31e−01 2.27e−01

– 0.12 1.11 1.31 1.23

7.63e+00 5.68e+00 2.25e+00 7.85e−01 2.47e−01

– 0.43 1.33 1.52 1.67

References [1] J. Wang, X. Ye, A weak Galerkin finite element method for the Stokes equations, Adv. Comput. Math. 42 (1) (2016) 155–174, http: //dx.doi.org/10.1007/s10444-015-9415-2, arXiv:13022707v1. [2] J. Wang, X. Ye, A weak Galerkin mixed finite element method for second-order ellliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115, available at arXiv:11042897vl. [3] J. Wang, X. Ye, A weak Galerkin mixed finite element method for second-order elliptic problems, Math. Comp. 83 (2014) 2101–2126. [4] C. Wang, J. Wang, A primal–dual weak Galerkin finite element method for second elliptic equations in non-divergence form, Math. Comp. (2017) http://dx.doi.org/10.1090/mcom/3220. [5] R. Wang, X. Wang, Q. Zhai, R. Zhang, A weak Galerkin finite element scheme for solving the stationary Stokes equations, J. Comput. Appl. Math. 302 (2016) 171–185. [6] D. Li, C. Wang, J. Wang, Superconvergence of the gradient approximation for weak Galerkin finite element methods on nonuniform rectangular partitions, https://arxiv.org/pdf/1804.03998v2.pdf. [7] J.E. Welch, F.H. Harlow, J.P. Shannon, et al., The MAC method. A computing technique for solving viscous, incompressible, transient fluid-flow problems involving free surfaces, 1966. [8] J.C. Strikwerda, Finite difference methods for the Stokes and Navier–Stokes equations, SIAM J. Sci. Stat. Comput. 5 (1) (1984) 56–68. [9] M. Li, T. Tang, B. Fornberg, A compact fourth-order finite difference scheme for the steady incompressible Navier–Stokes equations, Internat. J. Numer. Methods Fluids 20 (10) (1995) 1137–1151. [10] R. Nicolaides, Flow discretization by complementary volume techniques, in: Proc. 9th AIAA CFD Meeting, AIAA paper, Buffalo, NY, 1989, pp. 89–1978. [11] H. Han, X. Wu, A new mixed finite element formulation and the MAC method for the Stokes equations, SIAM J. Numer. Anal. 35 (2) (1998) 560–571. [12] G. Kanschat, Divergence-free discontinuous Galerkin schemes for the Stokes equations and the MAC scheme, Int. J. Numer. Methods Fluids 56 (7) (2008) 941–950. [13] P.D. Minev, Remarks on the links between low-order DG methods and some finite-difference schemes for the Stokes problem, Internat. J. Numer. Methods Fluids 58 (3) (2008) 307–317. [14] R.E. Ewing, T. Lin, Y. Lin, On the accuracy of the finite volume element method based on piecewise linear polynomials, SIAM J. Numer. Anal. 39 (2002) 1865–1888. [15] A. Quarteroni, R. Ruiz-Baier, Analysis of a finite volume element method for the Stokes problem, Numer. Math. 118 (2011) 737–764. [16] Z. Cai, On the finite volume element method, Numer. Math. 58 (1991) 713–735. [17] X. Ye, A discontinuous finite volume method for the Stokes problems, SIAM J. Numer. Anal. 44 (2006) 183–198.

206

Y. Liu and J. Wang / Journal of Computational and Applied Mathematics 361 (2019) 176–206

[18] J. Li, Z. Chen, A new stabilized finite volume method for the stationary Stokes equations, Adv. Comput. Math. 30 (2009) 141–152. [19] S.C. Chou, Analysis and convergence of a covolume method for the generalized Stokes problem, Math. Comp. 66 (1997) 85–104. [20] M. Feistauer, J. Felcman, M. Lukà˘ecovà-Medvidóvà, G. Warnecke, Error estimates for a combined finite volume-finite element method for nonlinear convection–diffusion problems, SIAM J. Numer. Anal. 36 (1999) 1528–1548. [21] T. Gallouët, R. Herbin, J.C. Latché, A convergent finite element-finite volume scheme for the compressible Stokes problem. Part I: the isothermal case, Math. Comp. 78 (2009) 1333–1352. [22] R.E. Bank, D.J. Rose, Some error estimates for the box method, SIAM J. Numer. Anal. 24 (1987) 777–787. [23] L. El Alaoui, An adaptive finite volume box scheme for solving a class of nonlinear parabolic equations, Appl. Math. Lett. 22 (2009) 291–296. [24] V. Girault, A combined finite element and marker and cell method for solving Navier–Stokes equations, Numer. Math. 26 (1976) 39–59. [25] S. Nicaise, K. Djadel, Convergence analysis of a finite volume method for the Stokes system using non-conforming arguments, IMA J. Numer. Anal. 25 (2005) 523–548. [26] X. Ye, On the relationship between finite volume and finite element methods applied to the Stokes equations, Numer. Methods Partial Differential Equations 5 (2001) 440–453. [27] M. Cui, X. Ye, Superconvergence of finite volume methods for the Stokes equations, Numer. Methods Partial Differential Equations 25 (2009) 1212–1230. [28] I. Babu˘ska, The finie element method with lagrange multipliers, Numer. Math. 20 (1973) 173–192. [29] F. Brezzi, On the existence, uniqueness, and approximation of saddle point problems arising from lagrange multipliers, RAIRO 8 (1974) 129–151. [30] S. Brenner, R. Scott, Mathematical Theory of Finite Element Methods, Springer, 2007. [31] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [32] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, 1986. [33] M.D. Gunzburger, Finite Element Methods for Viscous Incompressible Flows: A Guide To Theory, Practice, and Algorithms, Academic, San Diego, 1989. [34] J. Li, S. Sun, The superconvergence phenomenon and proof of the MAC scheme for the Stokes equations on non-uniform rectangular meshes, J. Sci. Comput. 65 (1) (2015) 341–362. [35] J. Liu, Penalty-factor-free discontinuous Galerkin methods for 2-dim Stokes problems, SIAM J. Numer. Anal. 49 (5) (2011) 2165–2181. [36] J. Wang, Y. Wang, X. Ye, A robust numerical method for Stokes equations based on divergence-free h (div) finite element methods, SIAM J. Sci. Comput. 31 (4) (2009) 2784–2802.