- Email: [email protected]
A modified weak Galerkin finite element method
Journal of Computational and Applied Mathematics 271 (2014) 319–327
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A modified weak Galerkin finite element method X. Wang a , N.S. Malluwawadu a , F. Gao b , T.C. McMillan a,∗ a
Department of Mathematics and Statistics, University of Arkansas at Little Rock, Little Rock, AR 72204, United States
b
School of Mathematics, Shandong University, Jinan, Shandong 250100, China
article
abstract
info
Article history: Received 28 September 2013 Received in revised form 18 March 2014 Keywords: Weak Galerkin Finite element
In this paper we introduce a new discrete weak gradient operator and a new weak Galerkin (WG) finite element method for second order Poisson equations based on this new operator. This newly defined discrete weak gradient operator allows us to use a single stabilizer which is similar to the one used in the discontinuous Galerkin (DG) methods without having to worry about choosing a sufficiently large parameter. In addition, we will establish the optimal convergence rates and validate the results with numerical examples. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In this paper we describe a finite element method that is a combination of processes used in the weak Galerkin (WG) [1] and discontinuous Galerkin (DG) [2] methods for second order Poisson equations. First, we introduce a new discrete weak gradient operator. This newly defined discrete weak gradient operator allows us to draw from the strengths of both WG and DG. The bilinear form defined here, a(u, v), contains only one stabilizer term. In contrast to DG, once this method is implemented there is no need to experiment when applying it to different problems in order to find an α large enough to obtain convergence. On the other hand, the freedom to introduce a positive number as a parameter makes it possible to further optimize this method. The convergence rates of the present method are optimal. We consider as our model problem the boundary value problem
−1u = f u=g
in Ω
(1)
on ∂ Ω ,
where Ω ⊆ Rd (d = 2, 3) is a polygonal region. We seek an approximate solution to this problem by applying a finite element method. Since the main ideas are the same, we focus on the case d = 2. We begin with some notational definitions. Let Th be a shape regular triangulation of Ω with elements T and edges e. Let hT be the diameter of T and h = maxT ∈Th hT . Given two elements T1 and T2 , with a common edge e, we denote by {v} the function defined on e that is the average (vo |T1 + vo |T2 )/2. We denote by [[v]] the function on e that is the jump experienced by v across the edge e: [[v]] = vo |T1 n1 +vo |T2 n2 , where n1 and n2 are the outward unit normal vectors on edge e for elements T1 and T2 , respectively. Our weak formulation will use the following vector spaces of functions on Ω : Vh = {v : v|T ∈ Pk (T ) for T ∈ Th }, Vh0 = {v ∈ Vh : v|e = 0 for e ∈ ∂ Ω },
and
Wk = {q : q|T ∈ [Pk (T )] for T ∈ Th }. 2
∗
Corresponding author. Tel.: +1 5015698102. E-mail addresses: [email protected] (X. Wang), [email protected] (N.S. Malluwawadu), [email protected], [email protected] (F. Gao), [email protected] (T.C. McMillan). http://dx.doi.org/10.1016/j.cam.2014.04.014 0377-0427/© 2014 Elsevier B.V. All rights reserved.
320
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
The notation e ∈ ∂ Ω means that e is an edge inthe boundary of Ω . Define (v, w)D = D vw dA and (v, w)∂ D = ∂ D vw ds. The projection operators Qo and Qb , defined piecewise on the interior and boundary of each of the elements of Th , are the L2 -projection onto Pk (T ), and L2 -projection onto Pk (e), respectively. Let Rh be the projection operator onto [Pk−1 (T )]2 whose components are each the L2 -projection onto Pk−1 (T ). We now define the new discrete gradient used in the modified WG method. Definition 1.1. Given a partition Th of Ω and a piecewise smooth function u on Ω , for all T ∈ Th , the discrete gradient of u on T is the unique element ∇w u|T in [Pk−1 (T )]2 such that
(∇w u, q)T = −(u, ∇ · q)T + ({u}, q · n)∂ T for all q ∈ [Pk−1 (T )]2 , where n denotes the outward facing unit normal perpendicular to an e ∈ ∂ T . Note that, when u is continuous on Ω , {u} = u on ∂ T . Thus
(∇w u, q)T = −(u, ∇ · q)T + (u, q · n)∂ T = (∇ u, q)T (2) 2 for all q ∈ [Pk−1 (T )] . Note also that the discrete gradient defined here is different from the one defined in [1]. The weak formulation for our boundary value problem can now be given. Find uh ∈ Vh such that uh |e = Qb g for e ∈ ∂ Ω , and a(uh , v) =
(∇w uh , ∇w v)T + h−1 ([[uh ]], [[v]])e = (f , v) for all v ∈ Vh0 .
(3)
e
T
Theorem 1.2. There is a unique uh ∈ Vh0 satisfying a(uh , v) = (f , v) for all v ∈ Vh0 . Proof. Let f = g = 0. Since Vh0 is finite dimensional, it suffices to show that a(uh , v) = 0 for all v ∈ Vh0 implies that uh = 0 on Ω . Under this assumption, let v = uh . Then a(uh , uh ) =
(∇w uh , ∇w uh )T +
h−1 ([[uh ]], [[uh ]])e = 0.
e
T
Since all the terms in this sum are nonnegative, they must be zero. Since ([[uh ]], [[uh ]])e = 0, and hence [[uh ]] = 0, on each edge e, uh is continuous on Ω . It follows from (∇w uh , ∇w uh )T = 0 on each T ∈ Th that ∇ uh = ∇w uh = 0 on each T . Therefore uh is constant on each T . Since [[uh ]] = 0, uh has the same constant value on each T . Now, since g = 0 on ∂ Ω , uh = 0 on any boundary element T . Therefore, uh = 0 on Ω . The |||v||| is defined by |||v|||2 = a(v, v). The positivity property of ||| · ||| follows from the above theorem. That |||v||| has the additional properties of a norm is easily verified. 2. Properties of the discrete gradient In this section, we establish properties of the discrete gradient and the interaction with the projection operators that will be used in the error analysis. Lemma 2.1. For all u ∈ H 1 (Ω ), ∇w u = Rh ∇ u. Proof. We establish the result on each element T ∈ Th . Let q ∈ [Pk−1 (T )]2 . It follows from (2) that (∇w u, q)T = (∇ u, q)T = (Rh ∇ u, q)T . Therefore (∇w u − Rh ∇ u, q)T = 0 for all q ∈ [Pk−1 (T )]2 . In particular, letting q = ∇w u − Rh ∇ u ∈ [Pk−1 (T )]2 , yields (∇w u − Rh ∇ u, ∇w u − Rh ∇ u)T = 0. Therefore, ∇w u = Rh ∇ u on T for all T ∈ Th . Lemma 2.2. For all u ∈ H k+1 (Ω ),
12
∥∇w u − ∇ Qo u∥
2 T
12
=
∥Rh ∇ u − ∇ Qo u∥
2 T
≤ chk |u|k+1 .
T
T
Proof. The first equality follows from Lemma 2.1.
12
∥Rh ∇ u − ∇ Qo u∥
2 T
12
≤
T
∥ R h ∇ u − ∇ u∥
T
≤ chk |u|k+1 . Lemma 2.3. For all u ∈ H k+1 (Ω ),
21
e
∥Qo u − u∥
2 e
1
≤ chk+ 2 |u|k+1 .
2 T
12
+
T
∥∇ u − ∇ Qo u∥
2 T
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
321
Proof. It follows from the trace inequality that
21
∥Qo u − u∥
12
≤c (h−1 ∥Qo u − u∥2T + h∥∇(Qo u − u)∥2T )
2 e
e
T
≤ ch
k+ 21
|u|k+1 .
Lemma 2.4. For all w ∈ Wk−1 and for all v such that v|T ∈ H 1 (T ),
(∇v, w)T = (∇w v, w)T + ([[v]], {w})e . T
e
T
Proof.
(∇w v, w)T = − (v, ∇ · w)T + ({v}, w · n)∂ T , T
T
and
T
(∇v, w)T = − (v, ∇ · w)T + (v, w · n)∂ T . T
T
T
Subtracting the first of these equations from the second yields
(∇v − ∇w v, w)T =
([[v]], {w})e . e
T
Lemma 2.5. If u is the solution to (1), then
(∇w (u − uh ), ∇w v)T +
([[v]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e = 0, e
T
for all v ∈ Vh0 . Proof. Test Eq. (1) with v ∈ Vh0 :
−(1u, v) = (f , v). Therefore,
(f , v) =
[−(1u, v)T ] T
(∇ u · n, v)∂ T (∇ u, ∇v)T − = T
T
(∇ u · n, v)∂ T (Rh ∇ u, ∇v)T − = T
T
= (Rh ∇ u, ∇w v)T + ([[v]], {Rh ∇ u})e − (∇ u · n, v)∂ T by Lemma 2.4. e
T
T
Reindexing the last sum over the edges of the triangulation and using Lemma 2.1 yields
(∇w v, ∇w u)T + ([[v]], {Rh ∇ u − ∇ u})e = (f , v). T
(4)
e
From the weak formulation (3), we have
(h−1 [[uh ]], [[v]])e = (f , v). (∇w uh , ∇w v)T + T
e
Subtracting (5) from (4) completes the proof:
(∇w (u − uh ), ∇w v)T + ([[v]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e = 0. T
e
3. Error analysis Now, let E = u − uh . Theorem 3.1. There is a constant c ∈ R such that |||E ||| ≤ chk |u|k+1 .
(5)
322
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
Proof. We begin with
|||E |||2 = a(E , E ) = (∇w E , ∇w E )T + ([[E ]], h−1 [[E ]])e . e
T
Letting E˜ = Qo u − uh ,
(∇w E , ∇w E )T = (Rh ∇ u − ∇w uh , ∇w E )T (by Lemma 2.1) T
T
= (∇w Qo u − ∇w uh , ∇w E )T + (Rh ∇ u − ∇w Qo u, ∇w E )T T
T
(Rh ∇ u − ∇w Qo u, ∇w E )T , = (∇w E˜, ∇w E )T + T
T
and letting v = E˜ in Lemma 2.5, we get
([[E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e . (∇w E˜, ∇w E )T = − e
T
Our estimation for |||E ||| follows:
|||E |||2 =
(Rh ∇ u − ∇w Qo u, ∇w E )T ([[E ]], h−1 [[E ]])e − ([[E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e + e
=
e
T
([[E ]], −{Rh ∇ u − ∇ u} + h−1 ([[E ]] + [[uh ]]))e e
+
(Rh ∇ u − ∇w Qo u, ∇w E )T ([[E − E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e + e
=
T
([[E ]], −{Rh ∇ u − ∇ u} + h−1 [[u]])e + ([[u − Qo u]], {Rh ∇ u − ∇ u} + h−1 [[u − uh ]])e e
e
(Rh ∇ u − ∇w Qo u, ∇w E )T +
(6)
T
=
([[E ]], −{Rh ∇ u − ∇ u} + h−1 [[u − Qo u]])e + ([[u − Qo u]], {Rh ∇ u − ∇ u})e e
e
(Rh ∇ u − ∇w Qo u, ∇w E )T +
(7)
T
21
≤
∥[[E ]]∥
12
2 e
e
+h
e
12
∥[[u − Qo u]]∥
e
∥[[E ]]∥
2 e
∥{Rh ∇ u − ∇ u}∥
∥[[u − Qo u]]∥ 12
2 e
2 e
e
12
12
2 e
e
+
∥{Rh ∇ u − ∇ u}∥
12
−1
2 e
+
e
∥Rh ∇ u − ∇w Qo u∥
12
2 T
T
∥∇w E ∥
. (8)
2 T
T
Eq. (6) follows from the fact that [[u]] = [[u − uh ]] + [[uh ]] = [[E ]] + [[uh ]], and Eq. (7) follows from [[u]] = 0. We will deal with each of the four terms in (8) separately. For the first term we have
21
∥[[E ]]∥
12
2 e
e
∥{Rh ∇ u − ∇ u}∥
2 e
12
=
e
h
−1
∥[[E ]]∥
21
2 e
e
h∥{Rh ∇ u − ∇ u}∥
2 e
e
12
1 2
≤ |||E |||h c
h
−1
2 T
∥Rh ∇ u − ∇ u∥ +
T
h∥∇(Rh ∇ u − ∇ u)∥
T
≤ |||E |||chk |u|k+1 .
(9)
For the second term, we have
12
h
−1
e
∥[[E ]]∥
2 e
12 e
∥[[u − Qo u]]∥
2 e
21
=h
−1
e
h
−1
2 T
∥[[E ]]∥
2 e
12
h∥[[u − Qo u]]∥
e
≤ |||E |||chk |u|k+1 (by Lemma 2.3).
2 e
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
323
For the third term,
21
12
∥[[u − Qo u]]∥2e
e
∥{Rh ∇ u − ∇ u}∥2e
e
12
≤ ch
k+ 12
|u|k+1
∥{Rh ∇ u − ∇ u}∥
(by Lemma 2.3)
2 e
e
12
≤ ch
k+ 12
|u|k+1
h
−1
2 T
∥Rh ∇ u − ∇ u∥ +
T
h∥∇(Rh ∇ u − ∇ u)∥
2 T
T
≤ ch2k |u|2k+1 .
(10)
For the fourth term,
12
12
∥∇w E ∥
2 T
∥Rh ∇ u − ∇w Qo u∥
2 T
=
T
T
12
∥∇w E ∥
2 T
12
T
∥Rh ∇ u − ∇ Qo u∥
2 T
T
≤ |||E |||chk |u|k+1 (by Lemma 2.2).
(11)
Lines (9) and (10) result from applications of the trace inequality. Now combining our results for the four terms in (8) yields the estimate
|||E |||2 ≤ |||E |||chk |u|k+1 + ch2k |u|2k+1 . Now, using the inequality 2ab ≤ a2 + b2 , yields
|||E ||| ≤ chk |u|k+1 . To obtain an L2 estimate of the error E = u − uh consider
− 1φ = E
on Ω ,
(12)
where φ ∈ H02 (Ω ). Assume, for some c > 0,
∥φ∥2 ≤ c ∥E ∥. Now, testing (12) with E , we have
∥E ∥2 = (−1φ, E ) = (∇φ, ∇ E )T − (∇φ · n, E )∂ T T
T
= (∇φ, ∇ E )T − (∇φ, [[E ]])e e
T
= (Rh ∇φ, ∇ E )T − (∇φ, [[E ]])e + (∇φ − Rh ∇φ, ∇ E )T e
T
T
=
(Rh ∇φ, ∇w E )T + ([[E ]], {Rh ∇φ})e − (∇φ, [[E ]])e + (∇φ − Rh ∇φ, ∇ E )T
=
(∇w φ, ∇w E )T + ([[E ]], {Rh ∇φ − ∇φ})e
e
T
T
e
T
+
e
(∇φ − Rh ∇φ, Rh ∇ E )T + (∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T T
T
= (∇w Qo φ, ∇w E )T + (∇w (φ − Qo φ), ∇w E )T T
T
+ ([[E ]], {Rh ∇φ − ∇φ})e + (∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T . e
We will bound each term in (13).
T
(13)
324
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
Applying Lemma 2.5 to the first term in (13),
−1 (∇w Qo φ, ∇w (u − uh ))T = ([[Qo φ − φ]], {Rh ∇ u − ∇ u} + h [[uh − u]])e T e 1 12 2 2 −2 2 2 ∥{Rh ∇ u − ∇ u}∥e + h ∥[[uh − u]]∥e ≤ ∥[[Qo φ − φ]]∥e e
e
12
≤c
h
−1
∥Qo φ − φ∥ + h∥∇(Qo φ − φ)∥ 2 T
2 T
T
12 1 · + h−1 |||E |||2 2 h−1 ∥Rh ∇ u − ∇ u∥2T + h∥∇(Rh ∇ u − ∇ u)∥2T T 1
≤ ch3/2 ∥φ∥2 hk− 2 |u|k+1 ≤ chk+1 |u|k+1 ∥E ∥. For the second term, we have
21 12 2 2 ∥∇w (φ − Qo φ)∥T ∥∇w E ∥T (∇w (φ − Qo φ), ∇w E )T ≤ T T T ≤ ch∥φ∥2 |||E ||| ≤ chk+1 |u|k+1 ∥E ∥. For the third term, we have
21 12 −1 2 2 h ∥[[E ]]∥e h ∥{Rh ∇φ − ∇φ}∥e ([[E ]], {Rh ∇φ − ∇φ})e ≤ e e e 1 2 −1 2 2 ≤ |||E ||| h h ∥Rh ∇φ − ∇φ∥T + h∥∇(Rh ∇φ − ∇φ)∥T T
≤ |||E |||ch∥φ∥2 ≤ chk+1 |u|k+1 ∥E ∥. For the fourth term, we have
12
(∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T ≤ T
∥∇φ − Rh ∇φ∥
2 T
T
12
∥∇ E − Rh ∇ E ∥
2 T
T
≤ ch∥φ∥2 hk |E |k+1 ≤ ch∥E ∥hk |u|k+1 ≤ chk+1 |u|k+1 ∥E ∥, st using the fact that the (k + 1) derivative of uh = 0. Combining these results, we have
∥E ∥2 ≤ chk+1 ∥E ∥|u|k+1 . This establishes Theorem 3.2. There is a constant c ∈ R such that ∥E ∥ ≤ chk+1 |u|k+1 . 4. Numerical experiments In this section, we give three numerical examples using scheme (3) constructed in Section 1 to verify the error estimate in Theorem 3.2. We construct triangular mesh as follows. First we partition the square domain Ω = [0, 1] × [0, 1] into N × N subsquares uniformly to obtain the square mesh. Next we divide each square element into two triangles by the diagonal line with a positive slope, completing the construction of a triangular mesh. Let h = 1/N (N = 2, 4, 8, 16, 32, 64) be mesh sizes for different triangular meshes. All of the examples given below will use these triangulations of Ω , and will apply the method to find a solution uh = {uo , ub } where u0 |T ∈ P1 (T ), and ub |e ∈ P1 (e).
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
325
Table 1 Discrete L2 -norm error and convergence rate for Example 4.1. N
∥u − uh ∥
Order
2 4 8 16 32 64
3.1239E−02 1.1339E−02 3.3681E−03 9.1461E−04 2.3842E−04 6.0878E−05
1.4621 1.7513 1.8807 1.9396 1.9695
Table 2 Discrete L2 -norm error and convergence rate Example 4.2. N
∥u − uh ∥
Order
2 4 8 16 32 64
8.0959E−03 1.8468E−03 4.0990E−04 8.8522E−05 1.9811E−05 4.6291E−06
2.1321 2.1717 2.2112 2.1598 2.0975
Table 3 Discrete L2 -norm error and convergence rate Example 4.3. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.6041E−02 3.7817E−03 8.2975E−04 1.6945E−04 3.4859E−05 7.4966E−06
2.0846 2.1883 2.2918 2.2813 2.2172
4.1. Homogeneous boundary cases We consider the following elliptic problem
−1u = f , u = 0,
in Ω
(14)
on ∂ Ω .
In all three numerical examples, we choose source term f (x) according to the corresponding analytical solution of each example. Example 4.1. The analytical solution to (14) is u = sin(π x) sin(π y). The finite element method (3), with the different mesh sizes h = 1/N (N = 2, 4, 8, 16, 32), is applied, and the corresponding discrete L2 -norms errors and convergence rates are listed in Table 1. Example 4.2. The analytical solution to (14) is u = x(1 − x)y(1 − y). Numerical error results and convergence rate are listed in Table 2 and are based on the same mesh sizes as that of the first example. Example 4.3. The analytical solution to (14) is u = x(1 − x)y(1 − y) exp(x − y). Numerical error results and convergence rate are listed in Table 3 based on the same mesh sizes as that of the first example. All the three numerical examples given above are in good agreement with the theoretical analysis in Section 3, which validates that the WG finite element method (3) is stable and convergent with second order convergence rate in discrete L2 norm.
326
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327 Table 4 Discrete L2 -norm error and convergence rate for Example 4.4. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.6701E−01 5.4979E−02 1.5817E−02 4.3953E−03 1.1840E−03 3.0968E−04
1.6030 1.7974 1.8475 1.8922 1.9349
Table 5 Discrete L2 -norm error and convergence rate Example 4.5. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.2026E−02 4.0292E−03 1.1195E−03 3.0289E−04 8.0659E−05 2.1013E−05
1.5776 1.8476 1.8860 1.9089 1.9405
Table 6 Discrete L2 -norm error and convergence rate Example 4.6. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.4564E−02 5.2053E−03 1.4748E−03 3.9698E−04 1.0504E−04 2.7270E−05
1.4844 1.8195 1.8933 1.9181 1.9456
4.2. Nonhomogeneous boundary cases In this subsection, we consider examples of the same elliptic problem (1) with a nonhomogeneous boundary condition. Example 4.4. The analytical solution to (1) is u = sin(π x) sin(π y) + x. The finite element method (3), with the different mesh sizes h = 1/N (N = 2, 4, 8, 16, 32) is applied, and the corresponding discrete L2 -norms errors and convergence rates are listed in Table 4. Example 4.5. The analytical solution to (1) is u = x(1 − x)y(1 − y) + x. Numerical error results and convergence rate are listed in Table 5 based on the same mesh sizes as that of the first example. Example 4.6. The analytical solution to (1) is u = x(1 − x)y(1 − y) exp(x − y) + x. Numerical error results and convergence rate are listed in Table 6 based on the same mesh sizes as that of the first example. All the three numerical examples given above are in good agreement with the theoretical analysis in Section 3, which validates that the WG finite element method (3) is stable and convergent with second order convergence rate in discrete L2 norm for nonhomogeneous boundary case.
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
327
Acknowledgment The third author’s research was partially supported by the Natural Science Foundation of Shandong Province of China grant ZR2013AM023. References [1] J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, arXiv:1104.2897v1. [2] D.N. Arnold, F. Brezzi, B. Cockburn, D. Marini, Discontinuous Galerkin methods for elliptic problems, in: B. Cockburn, G. Karniadakis, C.W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, in: Lecture Notes in Computational Science and Engineering, vol. 11, SpringerVerlag, New York, Heidelberg, Berlin, 2000, pp. 89–101.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A modified weak Galerkin finite element method X. Wang a , N.S. Malluwawadu a , F. Gao b , T.C. McMillan a,∗ a
Department of Mathematics and Statistics, University of Arkansas at Little Rock, Little Rock, AR 72204, United States
b
School of Mathematics, Shandong University, Jinan, Shandong 250100, China
article
abstract
info
Article history: Received 28 September 2013 Received in revised form 18 March 2014 Keywords: Weak Galerkin Finite element
In this paper we introduce a new discrete weak gradient operator and a new weak Galerkin (WG) finite element method for second order Poisson equations based on this new operator. This newly defined discrete weak gradient operator allows us to use a single stabilizer which is similar to the one used in the discontinuous Galerkin (DG) methods without having to worry about choosing a sufficiently large parameter. In addition, we will establish the optimal convergence rates and validate the results with numerical examples. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In this paper we describe a finite element method that is a combination of processes used in the weak Galerkin (WG) [1] and discontinuous Galerkin (DG) [2] methods for second order Poisson equations. First, we introduce a new discrete weak gradient operator. This newly defined discrete weak gradient operator allows us to draw from the strengths of both WG and DG. The bilinear form defined here, a(u, v), contains only one stabilizer term. In contrast to DG, once this method is implemented there is no need to experiment when applying it to different problems in order to find an α large enough to obtain convergence. On the other hand, the freedom to introduce a positive number as a parameter makes it possible to further optimize this method. The convergence rates of the present method are optimal. We consider as our model problem the boundary value problem
−1u = f u=g
in Ω
(1)
on ∂ Ω ,
where Ω ⊆ Rd (d = 2, 3) is a polygonal region. We seek an approximate solution to this problem by applying a finite element method. Since the main ideas are the same, we focus on the case d = 2. We begin with some notational definitions. Let Th be a shape regular triangulation of Ω with elements T and edges e. Let hT be the diameter of T and h = maxT ∈Th hT . Given two elements T1 and T2 , with a common edge e, we denote by {v} the function defined on e that is the average (vo |T1 + vo |T2 )/2. We denote by [[v]] the function on e that is the jump experienced by v across the edge e: [[v]] = vo |T1 n1 +vo |T2 n2 , where n1 and n2 are the outward unit normal vectors on edge e for elements T1 and T2 , respectively. Our weak formulation will use the following vector spaces of functions on Ω : Vh = {v : v|T ∈ Pk (T ) for T ∈ Th }, Vh0 = {v ∈ Vh : v|e = 0 for e ∈ ∂ Ω },
and
Wk = {q : q|T ∈ [Pk (T )] for T ∈ Th }. 2
∗
Corresponding author. Tel.: +1 5015698102. E-mail addresses: [email protected] (X. Wang), [email protected] (N.S. Malluwawadu), [email protected], [email protected] (F. Gao), [email protected] (T.C. McMillan). http://dx.doi.org/10.1016/j.cam.2014.04.014 0377-0427/© 2014 Elsevier B.V. All rights reserved.
320
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
The notation e ∈ ∂ Ω means that e is an edge inthe boundary of Ω . Define (v, w)D = D vw dA and (v, w)∂ D = ∂ D vw ds. The projection operators Qo and Qb , defined piecewise on the interior and boundary of each of the elements of Th , are the L2 -projection onto Pk (T ), and L2 -projection onto Pk (e), respectively. Let Rh be the projection operator onto [Pk−1 (T )]2 whose components are each the L2 -projection onto Pk−1 (T ). We now define the new discrete gradient used in the modified WG method. Definition 1.1. Given a partition Th of Ω and a piecewise smooth function u on Ω , for all T ∈ Th , the discrete gradient of u on T is the unique element ∇w u|T in [Pk−1 (T )]2 such that
(∇w u, q)T = −(u, ∇ · q)T + ({u}, q · n)∂ T for all q ∈ [Pk−1 (T )]2 , where n denotes the outward facing unit normal perpendicular to an e ∈ ∂ T . Note that, when u is continuous on Ω , {u} = u on ∂ T . Thus
(∇w u, q)T = −(u, ∇ · q)T + (u, q · n)∂ T = (∇ u, q)T (2) 2 for all q ∈ [Pk−1 (T )] . Note also that the discrete gradient defined here is different from the one defined in [1]. The weak formulation for our boundary value problem can now be given. Find uh ∈ Vh such that uh |e = Qb g for e ∈ ∂ Ω , and a(uh , v) =
(∇w uh , ∇w v)T + h−1 ([[uh ]], [[v]])e = (f , v) for all v ∈ Vh0 .
(3)
e
T
Theorem 1.2. There is a unique uh ∈ Vh0 satisfying a(uh , v) = (f , v) for all v ∈ Vh0 . Proof. Let f = g = 0. Since Vh0 is finite dimensional, it suffices to show that a(uh , v) = 0 for all v ∈ Vh0 implies that uh = 0 on Ω . Under this assumption, let v = uh . Then a(uh , uh ) =
(∇w uh , ∇w uh )T +
h−1 ([[uh ]], [[uh ]])e = 0.
e
T
Since all the terms in this sum are nonnegative, they must be zero. Since ([[uh ]], [[uh ]])e = 0, and hence [[uh ]] = 0, on each edge e, uh is continuous on Ω . It follows from (∇w uh , ∇w uh )T = 0 on each T ∈ Th that ∇ uh = ∇w uh = 0 on each T . Therefore uh is constant on each T . Since [[uh ]] = 0, uh has the same constant value on each T . Now, since g = 0 on ∂ Ω , uh = 0 on any boundary element T . Therefore, uh = 0 on Ω . The |||v||| is defined by |||v|||2 = a(v, v). The positivity property of ||| · ||| follows from the above theorem. That |||v||| has the additional properties of a norm is easily verified. 2. Properties of the discrete gradient In this section, we establish properties of the discrete gradient and the interaction with the projection operators that will be used in the error analysis. Lemma 2.1. For all u ∈ H 1 (Ω ), ∇w u = Rh ∇ u. Proof. We establish the result on each element T ∈ Th . Let q ∈ [Pk−1 (T )]2 . It follows from (2) that (∇w u, q)T = (∇ u, q)T = (Rh ∇ u, q)T . Therefore (∇w u − Rh ∇ u, q)T = 0 for all q ∈ [Pk−1 (T )]2 . In particular, letting q = ∇w u − Rh ∇ u ∈ [Pk−1 (T )]2 , yields (∇w u − Rh ∇ u, ∇w u − Rh ∇ u)T = 0. Therefore, ∇w u = Rh ∇ u on T for all T ∈ Th . Lemma 2.2. For all u ∈ H k+1 (Ω ),
12
∥∇w u − ∇ Qo u∥
2 T
12
=
∥Rh ∇ u − ∇ Qo u∥
2 T
≤ chk |u|k+1 .
T
T
Proof. The first equality follows from Lemma 2.1.
12
∥Rh ∇ u − ∇ Qo u∥
2 T
12
≤
T
∥ R h ∇ u − ∇ u∥
T
≤ chk |u|k+1 . Lemma 2.3. For all u ∈ H k+1 (Ω ),
21
e
∥Qo u − u∥
2 e
1
≤ chk+ 2 |u|k+1 .
2 T
12
+
T
∥∇ u − ∇ Qo u∥
2 T
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
321
Proof. It follows from the trace inequality that
21
∥Qo u − u∥
12
≤c (h−1 ∥Qo u − u∥2T + h∥∇(Qo u − u)∥2T )
2 e
e
T
≤ ch
k+ 21
|u|k+1 .
Lemma 2.4. For all w ∈ Wk−1 and for all v such that v|T ∈ H 1 (T ),
(∇v, w)T = (∇w v, w)T + ([[v]], {w})e . T
e
T
Proof.
(∇w v, w)T = − (v, ∇ · w)T + ({v}, w · n)∂ T , T
T
and
T
(∇v, w)T = − (v, ∇ · w)T + (v, w · n)∂ T . T
T
T
Subtracting the first of these equations from the second yields
(∇v − ∇w v, w)T =
([[v]], {w})e . e
T
Lemma 2.5. If u is the solution to (1), then
(∇w (u − uh ), ∇w v)T +
([[v]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e = 0, e
T
for all v ∈ Vh0 . Proof. Test Eq. (1) with v ∈ Vh0 :
−(1u, v) = (f , v). Therefore,
(f , v) =
[−(1u, v)T ] T
(∇ u · n, v)∂ T (∇ u, ∇v)T − = T
T
(∇ u · n, v)∂ T (Rh ∇ u, ∇v)T − = T
T
= (Rh ∇ u, ∇w v)T + ([[v]], {Rh ∇ u})e − (∇ u · n, v)∂ T by Lemma 2.4. e
T
T
Reindexing the last sum over the edges of the triangulation and using Lemma 2.1 yields
(∇w v, ∇w u)T + ([[v]], {Rh ∇ u − ∇ u})e = (f , v). T
(4)
e
From the weak formulation (3), we have
(h−1 [[uh ]], [[v]])e = (f , v). (∇w uh , ∇w v)T + T
e
Subtracting (5) from (4) completes the proof:
(∇w (u − uh ), ∇w v)T + ([[v]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e = 0. T
e
3. Error analysis Now, let E = u − uh . Theorem 3.1. There is a constant c ∈ R such that |||E ||| ≤ chk |u|k+1 .
(5)
322
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
Proof. We begin with
|||E |||2 = a(E , E ) = (∇w E , ∇w E )T + ([[E ]], h−1 [[E ]])e . e
T
Letting E˜ = Qo u − uh ,
(∇w E , ∇w E )T = (Rh ∇ u − ∇w uh , ∇w E )T (by Lemma 2.1) T
T
= (∇w Qo u − ∇w uh , ∇w E )T + (Rh ∇ u − ∇w Qo u, ∇w E )T T
T
(Rh ∇ u − ∇w Qo u, ∇w E )T , = (∇w E˜, ∇w E )T + T
T
and letting v = E˜ in Lemma 2.5, we get
([[E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e . (∇w E˜, ∇w E )T = − e
T
Our estimation for |||E ||| follows:
|||E |||2 =
(Rh ∇ u − ∇w Qo u, ∇w E )T ([[E ]], h−1 [[E ]])e − ([[E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e + e
=
e
T
([[E ]], −{Rh ∇ u − ∇ u} + h−1 ([[E ]] + [[uh ]]))e e
+
(Rh ∇ u − ∇w Qo u, ∇w E )T ([[E − E˜]], {Rh ∇ u − ∇ u} − h−1 [[uh ]])e + e
=
T
([[E ]], −{Rh ∇ u − ∇ u} + h−1 [[u]])e + ([[u − Qo u]], {Rh ∇ u − ∇ u} + h−1 [[u − uh ]])e e
e
(Rh ∇ u − ∇w Qo u, ∇w E )T +
(6)
T
=
([[E ]], −{Rh ∇ u − ∇ u} + h−1 [[u − Qo u]])e + ([[u − Qo u]], {Rh ∇ u − ∇ u})e e
e
(Rh ∇ u − ∇w Qo u, ∇w E )T +
(7)
T
21
≤
∥[[E ]]∥
12
2 e
e
+h
e
12
∥[[u − Qo u]]∥
e
∥[[E ]]∥
2 e
∥{Rh ∇ u − ∇ u}∥
∥[[u − Qo u]]∥ 12
2 e
2 e
e
12
12
2 e
e
+
∥{Rh ∇ u − ∇ u}∥
12
−1
2 e
+
e
∥Rh ∇ u − ∇w Qo u∥
12
2 T
T
∥∇w E ∥
. (8)
2 T
T
Eq. (6) follows from the fact that [[u]] = [[u − uh ]] + [[uh ]] = [[E ]] + [[uh ]], and Eq. (7) follows from [[u]] = 0. We will deal with each of the four terms in (8) separately. For the first term we have
21
∥[[E ]]∥
12
2 e
e
∥{Rh ∇ u − ∇ u}∥
2 e
12
=
e
h
−1
∥[[E ]]∥
21
2 e
e
h∥{Rh ∇ u − ∇ u}∥
2 e
e
12
1 2
≤ |||E |||h c
h
−1
2 T
∥Rh ∇ u − ∇ u∥ +
T
h∥∇(Rh ∇ u − ∇ u)∥
T
≤ |||E |||chk |u|k+1 .
(9)
For the second term, we have
12
h
−1
e
∥[[E ]]∥
2 e
12 e
∥[[u − Qo u]]∥
2 e
21
=h
−1
e
h
−1
2 T
∥[[E ]]∥
2 e
12
h∥[[u − Qo u]]∥
e
≤ |||E |||chk |u|k+1 (by Lemma 2.3).
2 e
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
323
For the third term,
21
12
∥[[u − Qo u]]∥2e
e
∥{Rh ∇ u − ∇ u}∥2e
e
12
≤ ch
k+ 12
|u|k+1
∥{Rh ∇ u − ∇ u}∥
(by Lemma 2.3)
2 e
e
12
≤ ch
k+ 12
|u|k+1
h
−1
2 T
∥Rh ∇ u − ∇ u∥ +
T
h∥∇(Rh ∇ u − ∇ u)∥
2 T
T
≤ ch2k |u|2k+1 .
(10)
For the fourth term,
12
12
∥∇w E ∥
2 T
∥Rh ∇ u − ∇w Qo u∥
2 T
=
T
T
12
∥∇w E ∥
2 T
12
T
∥Rh ∇ u − ∇ Qo u∥
2 T
T
≤ |||E |||chk |u|k+1 (by Lemma 2.2).
(11)
Lines (9) and (10) result from applications of the trace inequality. Now combining our results for the four terms in (8) yields the estimate
|||E |||2 ≤ |||E |||chk |u|k+1 + ch2k |u|2k+1 . Now, using the inequality 2ab ≤ a2 + b2 , yields
|||E ||| ≤ chk |u|k+1 . To obtain an L2 estimate of the error E = u − uh consider
− 1φ = E
on Ω ,
(12)
where φ ∈ H02 (Ω ). Assume, for some c > 0,
∥φ∥2 ≤ c ∥E ∥. Now, testing (12) with E , we have
∥E ∥2 = (−1φ, E ) = (∇φ, ∇ E )T − (∇φ · n, E )∂ T T
T
= (∇φ, ∇ E )T − (∇φ, [[E ]])e e
T
= (Rh ∇φ, ∇ E )T − (∇φ, [[E ]])e + (∇φ − Rh ∇φ, ∇ E )T e
T
T
=
(Rh ∇φ, ∇w E )T + ([[E ]], {Rh ∇φ})e − (∇φ, [[E ]])e + (∇φ − Rh ∇φ, ∇ E )T
=
(∇w φ, ∇w E )T + ([[E ]], {Rh ∇φ − ∇φ})e
e
T
T
e
T
+
e
(∇φ − Rh ∇φ, Rh ∇ E )T + (∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T T
T
= (∇w Qo φ, ∇w E )T + (∇w (φ − Qo φ), ∇w E )T T
T
+ ([[E ]], {Rh ∇φ − ∇φ})e + (∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T . e
We will bound each term in (13).
T
(13)
324
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
Applying Lemma 2.5 to the first term in (13),
−1 (∇w Qo φ, ∇w (u − uh ))T = ([[Qo φ − φ]], {Rh ∇ u − ∇ u} + h [[uh − u]])e T e 1 12 2 2 −2 2 2 ∥{Rh ∇ u − ∇ u}∥e + h ∥[[uh − u]]∥e ≤ ∥[[Qo φ − φ]]∥e e
e
12
≤c
h
−1
∥Qo φ − φ∥ + h∥∇(Qo φ − φ)∥ 2 T
2 T
T
12 1 · + h−1 |||E |||2 2 h−1 ∥Rh ∇ u − ∇ u∥2T + h∥∇(Rh ∇ u − ∇ u)∥2T T 1
≤ ch3/2 ∥φ∥2 hk− 2 |u|k+1 ≤ chk+1 |u|k+1 ∥E ∥. For the second term, we have
21 12 2 2 ∥∇w (φ − Qo φ)∥T ∥∇w E ∥T (∇w (φ − Qo φ), ∇w E )T ≤ T T T ≤ ch∥φ∥2 |||E ||| ≤ chk+1 |u|k+1 ∥E ∥. For the third term, we have
21 12 −1 2 2 h ∥[[E ]]∥e h ∥{Rh ∇φ − ∇φ}∥e ([[E ]], {Rh ∇φ − ∇φ})e ≤ e e e 1 2 −1 2 2 ≤ |||E ||| h h ∥Rh ∇φ − ∇φ∥T + h∥∇(Rh ∇φ − ∇φ)∥T T
≤ |||E |||ch∥φ∥2 ≤ chk+1 |u|k+1 ∥E ∥. For the fourth term, we have
12
(∇φ − Rh ∇φ, ∇ E − Rh ∇ E )T ≤ T
∥∇φ − Rh ∇φ∥
2 T
T
12
∥∇ E − Rh ∇ E ∥
2 T
T
≤ ch∥φ∥2 hk |E |k+1 ≤ ch∥E ∥hk |u|k+1 ≤ chk+1 |u|k+1 ∥E ∥, st using the fact that the (k + 1) derivative of uh = 0. Combining these results, we have
∥E ∥2 ≤ chk+1 ∥E ∥|u|k+1 . This establishes Theorem 3.2. There is a constant c ∈ R such that ∥E ∥ ≤ chk+1 |u|k+1 . 4. Numerical experiments In this section, we give three numerical examples using scheme (3) constructed in Section 1 to verify the error estimate in Theorem 3.2. We construct triangular mesh as follows. First we partition the square domain Ω = [0, 1] × [0, 1] into N × N subsquares uniformly to obtain the square mesh. Next we divide each square element into two triangles by the diagonal line with a positive slope, completing the construction of a triangular mesh. Let h = 1/N (N = 2, 4, 8, 16, 32, 64) be mesh sizes for different triangular meshes. All of the examples given below will use these triangulations of Ω , and will apply the method to find a solution uh = {uo , ub } where u0 |T ∈ P1 (T ), and ub |e ∈ P1 (e).
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
325
Table 1 Discrete L2 -norm error and convergence rate for Example 4.1. N
∥u − uh ∥
Order
2 4 8 16 32 64
3.1239E−02 1.1339E−02 3.3681E−03 9.1461E−04 2.3842E−04 6.0878E−05
1.4621 1.7513 1.8807 1.9396 1.9695
Table 2 Discrete L2 -norm error and convergence rate Example 4.2. N
∥u − uh ∥
Order
2 4 8 16 32 64
8.0959E−03 1.8468E−03 4.0990E−04 8.8522E−05 1.9811E−05 4.6291E−06
2.1321 2.1717 2.2112 2.1598 2.0975
Table 3 Discrete L2 -norm error and convergence rate Example 4.3. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.6041E−02 3.7817E−03 8.2975E−04 1.6945E−04 3.4859E−05 7.4966E−06
2.0846 2.1883 2.2918 2.2813 2.2172
4.1. Homogeneous boundary cases We consider the following elliptic problem
−1u = f , u = 0,
in Ω
(14)
on ∂ Ω .
In all three numerical examples, we choose source term f (x) according to the corresponding analytical solution of each example. Example 4.1. The analytical solution to (14) is u = sin(π x) sin(π y). The finite element method (3), with the different mesh sizes h = 1/N (N = 2, 4, 8, 16, 32), is applied, and the corresponding discrete L2 -norms errors and convergence rates are listed in Table 1. Example 4.2. The analytical solution to (14) is u = x(1 − x)y(1 − y). Numerical error results and convergence rate are listed in Table 2 and are based on the same mesh sizes as that of the first example. Example 4.3. The analytical solution to (14) is u = x(1 − x)y(1 − y) exp(x − y). Numerical error results and convergence rate are listed in Table 3 based on the same mesh sizes as that of the first example. All the three numerical examples given above are in good agreement with the theoretical analysis in Section 3, which validates that the WG finite element method (3) is stable and convergent with second order convergence rate in discrete L2 norm.
326
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327 Table 4 Discrete L2 -norm error and convergence rate for Example 4.4. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.6701E−01 5.4979E−02 1.5817E−02 4.3953E−03 1.1840E−03 3.0968E−04
1.6030 1.7974 1.8475 1.8922 1.9349
Table 5 Discrete L2 -norm error and convergence rate Example 4.5. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.2026E−02 4.0292E−03 1.1195E−03 3.0289E−04 8.0659E−05 2.1013E−05
1.5776 1.8476 1.8860 1.9089 1.9405
Table 6 Discrete L2 -norm error and convergence rate Example 4.6. N
∥u − uh ∥
Order
2 4 8 16 32 64
1.4564E−02 5.2053E−03 1.4748E−03 3.9698E−04 1.0504E−04 2.7270E−05
1.4844 1.8195 1.8933 1.9181 1.9456
4.2. Nonhomogeneous boundary cases In this subsection, we consider examples of the same elliptic problem (1) with a nonhomogeneous boundary condition. Example 4.4. The analytical solution to (1) is u = sin(π x) sin(π y) + x. The finite element method (3), with the different mesh sizes h = 1/N (N = 2, 4, 8, 16, 32) is applied, and the corresponding discrete L2 -norms errors and convergence rates are listed in Table 4. Example 4.5. The analytical solution to (1) is u = x(1 − x)y(1 − y) + x. Numerical error results and convergence rate are listed in Table 5 based on the same mesh sizes as that of the first example. Example 4.6. The analytical solution to (1) is u = x(1 − x)y(1 − y) exp(x − y) + x. Numerical error results and convergence rate are listed in Table 6 based on the same mesh sizes as that of the first example. All the three numerical examples given above are in good agreement with the theoretical analysis in Section 3, which validates that the WG finite element method (3) is stable and convergent with second order convergence rate in discrete L2 norm for nonhomogeneous boundary case.
X. Wang et al. / Journal of Computational and Applied Mathematics 271 (2014) 319–327
327
Acknowledgment The third author’s research was partially supported by the Natural Science Foundation of Shandong Province of China grant ZR2013AM023. References [1] J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, arXiv:1104.2897v1. [2] D.N. Arnold, F. Brezzi, B. Cockburn, D. Marini, Discontinuous Galerkin methods for elliptic problems, in: B. Cockburn, G. Karniadakis, C.W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, in: Lecture Notes in Computational Science and Engineering, vol. 11, SpringerVerlag, New York, Heidelberg, Berlin, 2000, pp. 89–101.