2021 UNC Greensboro PDE Conference, Electronic Journal of Differential Equations, Conference 26 (2022), pp. 123–138. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu PENALTY PARAMETER AND DUAL-WIND DISCONTINUOUS GALERKIN APPROXIMATION METHODS FOR ELLIPTIC SECOND ORDER PDES THOMAS LEWIS, AARON RAPP, YI ZHANG Abstract. This article analyzes the effect of the penalty parameter used in symmetric dual-wind discontinuous Galerkin (DWDG) methods for approx- imating second order elliptic partial differential equations (PDE). DWDG methods follow from the DG differential calculus framework that defines dis- crete differential operators used to replace the continuous differential operators when discretizing a PDE. We establish the convergence of the DWDG approxi- mation to a continuous Galerkin approximation as the penalty parameter tends towards infinity. We also test the influence of the regularity of the solution for elliptic second-order PDEs with regards to the relationship between the penalty parameter and the error for the DWDG approximation. Numerical experiments are provided to validate the theoretical results and to investigate the relationship between the penalty parameter and the L2-error. 1. Introduction Let Ω ⊂ Rd be a bounded convex polygonal domain, f ∈ L2(Ω), and g ∈ H1/2(∂Ω). We consider the second order elliptic partial differential problem: find u ∈ H2(Ω) such that −∆u = f in Ω, (1.1a) u = g on ∂Ω, (1.1b) where ∆u = ∑d i=1 ∂2 ∂x2 i u. Dual-Wind Discontinuous Galerkin (DWDG) methods have been applied to the second order elliptic PDE (1.1) as well as its Neumann boundary condition counterpart in [14, 10, 11]. The focus of these papers was to establish a priori results for the proposed DWDG methods initially analyzed in a discrete H1-space inspired by the weak form of (1.1): find u ∈ H1 g (Ω) such that (∇u,∇v)Ω = (f, v)Ω ∀v ∈ H1 0 (Ω), (1.2) where (v, w)Ω = ∫ Ω vw dx. In this paper we will further investigate properties of the DWDG method and the effects of adding a jump stabilization term since one of 2020 Mathematics Subject Classification. 65N30, 65N99. Key words and phrases. Discontinuous Galerkin methods; DWDG methods; penalty parameter; Poisson problem. ©2022 This work is licensed under a CC BY 4.0 license. Published August 25, 2022. 123 124 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 the key features of DWDG is the fact that it is naturally stable without penalizing jumps. There are many continuous Galerkin (CG) methods and discontinuous Galerkin (DG) methods that accurately approximate (1.2) with both weak and strong en- forcement of the boundary information. For DG methods, a jump penalization sta- bility term scaled by a penalty parameter γ is introduced in the discrete variational formulation to ensure coercivity; such terms are not needed for CG approximations. It was shown in [5, 6, 12] that a continuous Galerkin method is the limit of a DG interior penalty method as the penalty parameter tends towards infinity. Since DWDG methods are more similar to the more general class of flux-based DG meth- ods such as local discontinuous Galerkin (LDG) methods, a goal of this paper will be to show that DWDG methods will also approach the CG approximation as its limit as γ →∞. We also note that DWDG methods allow for natural enforcement of the Dirichlet boundary condition instead of the typical weak or strong approach. The penalty parameter γ is an artificial parameter that is traditionally intro- duced into the discrete formulation of a problem to guarantee the stability of the method in a discontinuous space. For interior-penalty methods, the penalty term naturally occurs in the derivation of a discrete Friedrich’s inequality (see [2, 4]). For LDG methods, the parameter does not arise as naturally, and it can be eliminated under certain assumptions in the formulation of the problem [9, 13]. For DWDG methods, the inclusion of an LDG-like penalty term was introduced to ensure sta- bility of the methods for non quasi-uniform meshes. Investigation of the stability and convergence of the dual-wind derivatives and the penalty term for second order elliptic PDEs found that the penalty term was not necessary for a shape-regular mesh. If the mesh was quasi-uniform and each triangle did not have more than one edge on the boundary, then the penalty term could be removed by setting γe = 0 for all edges of a mesh e, or it could be negative as long as γe > −C∗ for some constant C∗ > 0 (see [13, 14, 10]). The possibility of a non-positive penalty parameter indicates that DWDG meth- ods naturally weight jumps in the discontinuous approximation space without adding jump terms. The discrete derivatives ∇±h,g themselves control the jumps of the approximation eliminating the need for γe 6= 0 for all edges of the mesh e. This allows us to interpret γ as an artificial penalty parameter for the DWDG methods that can be controlled or eliminated. We refer the reader to [16] which contains a wide range of numerical experiments that first explored the relationship between the penalty parameter γ and the L2- error for the DWDG approximation. The numerical experiments in [16] showed a strictly increasing relationship between γ and the L2-error given that the solution was smooth and the approximation was found on a fine mesh. There was no sig- nificant impact from the minimal angle of the mesh, or whether the problem had homogeneous or non-homogeneous boundary conditions. The solutions for the nu- merical experiments in [16] were either in C∞(Ω) or in the test function space V rh . A second goal of this paper will be to extend the initial tests in [16] by experimen- tally studying the effect of the regularity of the solution to (1.2) on the relationship between γ and the L2-error. The rest of the paper is organized as follows. Section 2 introduces notation and briefly discusses the discrete DG interior calculus that will be used to formulate the symmetric dual-wind DG methods in Section 3. In Section 4, we will present EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 125 analysis that shows the DWDG methods will converge to a CG method as γ → ∞. Numerical experiments verifying the DWDG limit result will be presented in Section 5. Further investigations of the relationship between γ and the L2-error of the DWDG approximation with regards to the regularity of the solution of (1.2) will also be presented in Section 5. Throughout this article, C will denote a generic positive constant independent of the mesh size and γ. As such, the constant C can take different values at different occurrences. 2. Notation and DG differential calculus 2.1. Notation. We will use the standard space and function notation from [1] and [3] in this paper. Let Wm,p(Ω) denote the set of all Lp(Ω) functions whose distributional derivatives up to order m are in Lp(Ω) and Wm,p 0 (Ω) denote the set of Wm,p(Ω) functions whose traces vanish up to order m − 1 on ∂Ω. In the special case where p = 2, we denote this as Hm := Wm,2 and Hm 0 := Wm,2 0 . Bold-face format will be used for the corresponding vector-valued Sobolev spaces Wm,p(Ω) := [Wm,p(Ω)]d,Hm(Ω) := [Hm(Ω)]d, etc. We choose to introduce the discrete derivatives for a general dimension d ≥ 1 that will be used to formulate the numerical methods and establish our analytic results in Section 4. However, the numerical experiments presented in Section 5 will be done with d = 2. Let Th denote a shape-regular simplicial triangulation of Ω [3, 8]. Let EIh be the set of interior (d − 1)-dimensional simplices for the triangulation and EBh the set of boundary (d − 1)-dimensional simplices so that Eh := EIh ∪ EBh . We will denote the diameter of the simplex K ∈ Th as hK , and we set h := maxK∈Th hK . Let dK represent the diameter of the inscribed circle of a triangle K. A mesh Th is shape-regular if there is a constant C > 0 such that max K∈Th hK dK ≤ C. A mesh Th is quasi-uniform if there exist a constant ρ > 0 such that ρ < hK/hK′ < ρ−1 ∀K,K ′ ∈ Th. We define the piecewise vector spaces with respect to the triangulation Wm,p(Th) := ∏ K∈Th Wm,p(K), Wm,p(Th) := ∏ K∈Th Wm,p(K), etc. with the special cases Vh := W 1,1(Th) ∩ C0(Th) and Vh := [Vh]d. Let (v, w)Th := ∑ K∈Th ∫ K vw dx be the piecewise L2 inner product with respect to the triangulation, and let 〈v, w〉Sh := ∑ e∈Sh ∫ e vw ds be the piecewise L2 inner product over any subset Sh ⊆ Eh. To further simplify notation used later in this paper we will use 〈cev, w〉Sh := ∑ e∈Sh ∫ e cevw ds 126 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 where ce is a constant that depends on the edge e ∈ Sh ⊆ Eh. We will also use the following notation for the L2 norm over the triangularization in this paper: ‖v‖2L2(Th) := (v, v)Th and ‖v‖2L2(Sh) := 〈v, v〉Sh for Sh ⊆ Eh. For an integer r ≥ 1, we define the DG polynomial space V rh to be V rh := ∏ K∈Th Pr(K), where Pr(K) denotes the space of polynomials with domain K and degree not exceeding r. We will denote the DG polynomial space with zero-trace over ∂Ω as V rh,0 = {vh ∈ V rh : vh = 0 on ∂Ω}. As with the vector-valued Sobolev spaces, we will define the vector-valued polynomial spaces using bold-face format; Vr h := [V rh ]d and Vr h,0 := [V rh,0]d. For other spaces throughout this paper, we will use bold-face notation to indicate the vector-valued space. We note the inclusions V rh ⊂ Vh and Vr h ⊂ Vh. Let K+,K− ∈ Th, and let e = ∂K− ∩ ∂K+ ∈ EIh. Without loss of generality we will assume the global labeling number of K+ is larger than that of K−. Letting v± := v|K± , we define the jumps and averages across the (d−1)-dimensional simplex e as [[v]]e := v+ − v−, {v}e := 1 2 ( v+ + v− ) ∀v ∈ Vh. If K+ contains an edge in EBh , then for the boundary (d− 1)-dimensional simplex e = ∂K+ ∩ ∂Ω, we define [[v]]e := v+ and {v}e := v+. Set (n (1) e , n (2) e , . . . , n (d) e )t = ne := nK+ |e = −nK− |e to be the unit normal on e ∈ Eh. For i ∈ {1, 2, . . . , d} and v ∈ Vh, we define the labeling-dependent upwind trace operator, Q+ i (v), and the labeling dependent downwind trace operator, Q−i (v), on e ∈ EIh in the direction xi as Q±i (v) := {v} ± 1 2 sgn(n(i) e )[[v]], where sgn(n(i) e ) =  1 if n (i) e > 0, −1 if n (i) e < 0, 0 if n (i) e = 0, and Qi(v) := 1 2 ( Q−i (v) + Q+ i (v) ) . The operators Q−i (v) and Q+ i (v) can be re- garded as the “backward” and “forward” limit of v in the xi direction on e ∈ EIh, , respectively. On a boundary simplex e ∈ Ebh, we define Q±i (v) = Qi(v) = v. 2.2. Discontinuous Galerkin calculus. We use the trace operators Q−i , Q+ i , and Qi to define our DG discrete partials derivative operators ∂−h,xi , ∂ + h,xi , ∂h,xi : Vh → V hr that will be used to formulate the DWDG method for approximating (1.2). Definition 2.1. Let v ∈ Vh, g ∈ L1(∂Ω), i ∈ {1, 2, . . . , d}, let ∂xi be the usual (weak) partial derivative operator in the direction xi, and let n(i) be the piece- wise constant vector-valued function satisfying n(i)|e = [ne]i. The discrete partial derivative operators ∂±h,xi , ∂h,xi : Vh → V rh are defined by( ∂±h,xiv, ϕh ) Th := 〈 Q±i (v)n(i), [[ϕh]] 〉 Eh − ( v, ∂xiϕh ) Th ∀ϕh ∈ V rh ,( ∂h,xiv, ϕh ) Th := 〈 Qi(v)n(i), [[ϕh]] 〉 Eh − ( v, ∂xiϕh ) Th ∀ϕh ∈ V rh . EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 127 The discrete partial derivatives with given boundary data ∂±,gh,xi , ∂ g h,xi : Vh → V rh are defined as( ∂±,gh,xi v, ϕh ) Th := ( ∂±h,xiv, ϕh)Th + 〈(g − v)n(i), ϕh〉EBh ∀ϕh ∈ V rh ,( ∂ g h,xiv, ϕh ) Th := ( ∂h,xiv, ϕh)Th + 〈(g − v)n(i), ϕh〉EBh ∀ϕh ∈ V rh . Remark 2.2. The definition of ∂±,gh,xi and ∂ g h,xi is equivalent to setting Q±i (v) = g and Qi(v) = g on EBh implying that the trace data g is naturally incorporated into the discrete partial derivative operator. Definition 2.3. The discrete gradient operators ∇±h ,∇h,∇ ± h,g,∇h,g : Vh → Vr h are defined as[ ∇±h v ] i := ∂±h,xiv, [ ∇hv ] i := ∂h,xiv, [ ∇±h,gv ] i := ∂±,gh,xi v, [ ∇h,gv ] i := ∂ g h,xiv. The discrete divergence operators div±h ,divh : Vh → V rh are defined as div±h v := d∑ i=1 ∂±h,xi [v]i, divhv = 1 2 ( div+ h v + div−h v ) . We end this section by listing some properties of the DG discrete derivatives that will be used in the formulation of our proposed methods. From Definition 2.3, we can relate ∇±h,g and ∇±h,0 for all v, w ∈ Vh by( ∇±h,gv,ϕh)Th = ( ∇±h,0v,ϕh)Th + 〈g,ϕh · n〉EBh ∀ϕh ∈ Vr h, (2.1)( ∇±h,gv −∇ ± h,gw,ϕh ) Th = ( ∇±h,0(v − w),ϕh ) Th ∀ϕh ∈ Vr h. (2.2) The following discrete analog of integration by parts holds (i = 1, 2, . . . , d), see [11], (∂±h,xivh, ϕh)Th = −(vh, ∂ ∓ h,xi ϕh)Th + 〈vh, ϕhn(i)〉EBh ∀vh, ϕh ∈ V rh (2.3) yielding (∇±h vh,ϕh)Th = −(vh,div∓h ϕh)Th + 〈vh,ϕh · n〉EBh ∀vh ∈ V rh , ϕh ∈ Vr h, (2.4) (∇±h,0vh,ϕh)Th = −(vh,div∓h ϕh)Th ∀vh ∈ V rh , ϕh ∈ Vr h. (2.5) An important property of the DG discrete derivatives is that they are the L2 projections of the derivative of a function v if v ∈ H1(Ω). Lemma 2.4 ([11]). For any v ∈ H1(Ω), both ∂±h,xiv and ∂h,xiv are the L2 projec- tions of ∂xiv onto V rh ; that is, (∂±h,xiv, wh)Th = (∂h,xiv, wh)Th = (∂xiv, wh)Th for all wh ∈ V rh . Moreover, if v ∈ H1(Ω) satisfies v|∂Ω = g for some g ∈ L2(Ω), then both ∂±,gh,xi v and ∂ g h,xiv are the L2 projection of ∂xiv onto V rh . From the above lemma, it is possible to derive another useful identity for the discrete differential operators. Lemma 2.5 ([14]). Let v ∈ H2(Ω) with v = g on the boundary ∂Ω. Then it holds − (∆v, ϕh)Th = (∇h,gv,∇h,0ϕ)Th + 〈 { ∇h,gv −∇v } · n, [[ϕh]]〉Eh ∀ϕh ∈ V rh . (2.6) 128 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 3. Dual-wind discontinuous Galerkin approximation method for Poisson’s equation Recall that the strong form of Poisson’s equation (1.1) is: find a function u such that −∆u = f in Ω, u = g on ∂Ω. Now we define the discrete Laplacian operator ∆h,g : V rh −→ V rh by ∆h,gvh := div+ h ∇ − h,gvh + div−h ∇ + h,gvh 2 ∀vh ∈ V rh , (3.1) which involves an up-wind gradient operator and a down-wind gradient operator utilized in a symmetric way. The DWDG method for (1.1) seeks a function uγh ∈ V rh such that −∆h,gu γ h + jh,g(u γ h) = Phf, (3.2) where Phf ∈ V rh is the L2 projection of f onto V rh defined by (Phf, vh)Th = (f, vh)Th for all vh ∈ V rh and jh,g : Vh −→ V rh is a jump penalization/stabilization operator defined by (jh,g(v), wh)Th := γ〈h−1 e [[v]], [[wh]]〉Eh − γ〈h−1 e g, wh〉EBh ∀wh ∈ V rh . (3.3) Traditionally for DWDG methods, we can let γ = γe be a piecewise constant for each e ∈ Eh. However, for this paper we will assume that the “penalty” parameter is constant across all edges, including boundary edges. We say “penalty” because we can set γ ≤ 0 and still have an accurate approximation method [10, 14, 15, 16]. Using (2.5) and (2.1), a direct calculation shows that problem (3.2) is equivalent to finding uγh ∈ V rh such that Bh,γ(uγh, vh) = Fh,γ(vh) ∀vh ∈ V rh , (3.4) where Bh,γ(vh, wh) := 1 2 (( ∇+ h,0vh,∇ + h,0wh ) Th + ( ∇−h,0vh,∇ − h,0wh ) Th ) + γ〈h−1 e [[vh]], [[wh]]〉Eh ∀vh, wh ∈ V rh , (3.5) Fh,γ(v) := (f, v)Th + γ〈h−1 e g, v〉EBh − 〈g,∇h,0v · n〉EBh ∀v ∈ Vh. (3.6) We define the associated discrete energy norms by ‖vh‖2h := 1 2 ( ‖∇+ h,0vh‖ 2 L2(Ω) + ‖∇−h,0vh‖ 2 L2(Ω) ) , (3.7) |||vh|||2h,γ := Bh,γ(vh, vh) = ‖vh‖2h + γ ∑ e∈Eh ∥∥h−1/2 e [[vh]] ∥∥2 L2(e) . (3.8) Note that ‖|| · |||h,γ is a norm on V rh whenever γ ≥ 0 and that the bilinear form Bh,γ(·, ·) is stable with respect to ||| · |||h,γ [14, 15]. 4. DWDG converges to CG In this section, we prove analytically that the DWDG approximation (3.4) will converge to a continuous Galerkin (CG) approximation as γ →∞. For this section, EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 129 we will assume that γ > 0. Note that this only requires the mesh to be shape- regular. We will also assume that g = Πhg on ∂Ω where Πhg is the standard nodal interpolation of g into V rh . We define the spaces V ch = {v ∈ C(Ω) : v|K ∈ Pr(K), ∀K ∈ Th}, V ch,0 = {vh ∈ V ch : vh = 0 on ∂Ω}. The following lemma is adapted from [12, Lemma 3.1], and it shows an important equivalence of norms that will be used for our main analytical result. Note that we choose γ = 1 for (3.5) to define the energy norm in our analysis, and to state the following lemma based on the energy norm ||| · |||h,1. Lemma 4.1. There is a positive constant C independent of h and γ such that |||wh|||2V rh /V ch ≤ C ∑ e∈Eh ‖h−1/2 e [[wh]]‖2L2(e) ∀wh ∈ V rh /V ch , where ||| · |||V rh /V ch denotes the quotient norm |||wh|||V rh /V ch = inf v∈V ch |||wh + v|||h,1. Remark 4.2. The argument in the proof of [12, Lemma 3.1] is for an energy norm derived from a flux-based DG method. However, the proof does not invoke any special properties of the energy norm, and it can be generalized to any flux-based DG energy norm defined over the finite dimensional DG space V rh . Thus, we omit the proof. Recalling that we have assumed that g = Πhg on ∂Ω, the CG approximation corresponds to finding uch ∈ V ch with uch|∂Ω = Πhg such that B(uch, vh) = (f, vh) ∀vh ∈ V ch,0, (4.1) where B(vh, wh) = (∇vh,∇wh)Th . To prove our main result, we first show uγh satisfies a Galerkin orthogonality condition. Lemma 4.3. Let uγh be the solution to the DWDG method (3.4) and uch be the solution to the CG method (4.1). Then Bh,1(uch − u γ h, vh) = 0 ∀vh ∈ V ch,0. Proof. Let uγh be the solution (3.4), uch be the solution to (4.1), and vh ∈ V rh,0 ⊂ V rh . By (3.5), (3.4), the fact that [[vh]] = 0 for all e ∈ EIh, and the fact that vh = 0 on ∂Ω, we have Bh,1(uch − u γ h, vh) = Bh,1(uch, vh)−Bh,1(uγh, vh) = Bh,1(uch, vh)−Bh,γ(uγh, vh) + 〈γ − 1 he [[uγh]], [[vh]]〉Eh = Bh,1(uch, vh)− Fh,γ(vh) = 1 2 (( ∇+ h,0u c h,∇+ h,0vh ) Th + ( ∇−h,0u c h,∇−h,0vh ) Th ) + 〈h−1 e [[uch]], [[vh]]〉Eh − Fh,γ(vh) = 1 2 ( ∇+ h,0u c h,∇+ h,0vh ) Th + 1 2 ( ∇−h,0u c h,∇−h,0vh ) Th − Fh,γ(vh). (4.2) 130 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 By (2.1) and Lemma 2.4, equation (4.2) becomes Bh,1(uch − u γ h, vh) = 1 2 ( ∇+ h,gu c h,∇+ h,0vh ) Th + 1 2 ( ∇−h,gu c h,∇−h,0vh ) Th − 〈g,∇h,0vh · n〉EBh − Fh,γ(vh) = 1 2 ( ∇uch,∇+ h,0vh ) Th + 1 2 ( ∇uch,∇−h,0vh ) Th − 〈g,∇h,0vh · n〉EBh − Fh,γ(vh) (4.3) = ( ∇uch,∇h,0vh ) Th − 〈g,∇h,0vh · n〉EBh − Fh,γ(vh) = ( ∇uch,∇vh ) Th − 〈g,∇h,0vh · n〉EBh − Fh,γ(vh) = ( f, vh ) Th − 〈g,∇h,0vh · n〉EBh − Fh,γ(vh) Lastly, by applying (3.6) to (4.3), we have Bh,1(uch − u γ h, vh) = ( f, vh ) Th − 〈g,∇h,0vh · n〉EBh − (f, vh)Th − γ〈h−1 e g, vh〉EBh + 〈g,∇h,0vh · n〉EBh = −γ〈h−1 e g, vh〉EBh = 0. � Now we prove the main result of our paper. In [16], it was shown through nu- merical experiments that the L2-error was following a check-mark trend or strictly increasing as the the choice of γ increased. Based on the DWDG bilinear form (3.5), any jumps in the DWDG approximation will introduce more energy into in discrete system for large values of γ. Intuitively, this would encourage the DWDG approximation to minimize any jumps, making it resemble a continuous Galerkin approximation for large enough values of γ. In the next theorem, we will show that the difference of the DWDG approximation (3.4) and the CG approximation (4.1) in a discrete energy norm will be bounded, and the bound will explicitly depend on the choice of γ and h. Theorem 4.4. If u ∈ Hr+1(Ω) for r ≥ 1 is the solution to (1.1), uγh is the solution to the DWDG method (3.4), and uch is the solution to the CG method (4.1), then |||uch − u γ h|||h,1 ≤ Chr γ ‖u‖Hr+1(Ω). Proof. Let u ∈ Hr+1(Ω) be the solution to (1.1), let uγh be the solution to (3.4) and uch be the solution to (4.1). By Lemma 4.3 and the boundedness of ||| · |||h,1, we have |||uch − u γ h||| 2 h,1 = Bh,1(uch − u γ h, u c h − u γ h) = Bh,1(uch − u γ h, u c h − u γ h + v) ≤ |||uch − u γ h|||h,1|||u c h − u γ h + v|||h,1 (4.4) for any v ∈ V ch,0. Dividing both sides of (4.4) by |||uch − u γ h|||h,1, we obtain |||uch − u γ h|||h,1 ≤ |||u c h − u γ h + v|||h,1 (4.5) EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 131 for all v ∈ V ch,0. Since this holds for all v ∈ V ch,0, we have by Lemma 4.1 that |||uch − u γ h|||h,1 ≤ inf v∈V ch,0 |||uch − u γ h + v|||h,1 = |||uch − u γ h|||V rh /V ch,0 ≤ C ( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 . (4.6) By (3.5) we have γ ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) ≤ Bh,γ(uch − u γ h, u c h − u γ h) = Bh,γ(uch − u, uch − u γ h) +Bh,γ(u− uγh, u c h − u γ h). (4.7) Now we consider the first term on the right-hand side of (4.7). Recall our assump- tion Πhg = g, which implies uch − u ∈ H1 0 (Ω). By the Cauchy-Schwarz inequality and the fact that [[uch − u]]e = 0 for all e ∈ Eh, we have Bh,γ(uch − u, uch − u γ h) ≤ 1 2 ‖∇+ h,0(uch − u)‖L2(Th)‖∇+ h,0(uch − u γ h)‖L2(Th) + 1 2 ‖∇−h,0(uch − u)‖L2(Th)‖∇−h,0(uch − u γ h)‖L2(Th) + γ ∑ e∈Eh 〈h−1 e [[uch − u]], [[uch − u γ h]]〉L2(e) = 1 2 ‖∇+ h,0(uch − u)‖L2(Th)‖∇+ h,0(uch − u γ h)‖L2(Th) + 1 2 ‖∇−h,0(uch − u)‖L2(Th)‖∇−h,0(uch − u γ h)‖L2(Th). By Lemma 2.4, (3.8) with our assumption that γ > 0, and (4.6), we have Bh,γ(uch − u, uch − u γ h) = 1 2 ‖∇+ h,0(uch − u)‖L2(Th)‖∇+ h,0(uch − u γ h)‖L2(Th) + 1 2 ‖∇−h,0(uch − u)‖L2(Th)‖∇−h,0(uch − u γ h)‖L2(Th) = 1 2 ‖∇+ h,0(uch − u)‖L2(Th) ( ‖∇+ h,0(uch − u γ h)‖L2(Th) + ‖∇−h,0(uch − u γ h)‖L2(Th) ) ≤ 1 2 ‖uch − u‖h (2‖uch − u γ h‖h) ≤ |||uch − u|||h,1|||uch − u γ h|||h,1 ≤ C|||uch − u|||h,1 ( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 . (4.8) 132 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 We will next consider the second term on the right-hand side of (4.7). By (3.4) and (3.5), we have Bh,γ(u− uγh, u c h − u γ h) = Bh,γ(u, uch − u γ h)−Bh,γ(uγh, u c h − u γ h) = Bh,γ(u, uch − u γ h)− Fh,γ(uch − u γ h) = 1 2 (( ∇+ h,0u,∇ + h,0(uch − u γ h) ) Th + ( ∇−h,0u,∇ − h,0(uch − u γ h) ) Th ) + γ〈h−1 e u, (uch − u γ h)〉EBh − Fh,γ(uch − u γ h). (4.9) Next, by applying (3.6), (2.1), and Lemma 2.5 to (4.9), we have Bh,γ(u− uγh, u c h − u γ h) = 1 2 ( ∇+ h,0u,∇ + h,0(uch − u γ h) ) Th + 1 2 ( ∇−h,0u,∇ − h,0(uch − u γ h) ) Th + γ〈h−1 e u, uch − u γ h〉EBh − (f, uch − u γ h)Th − γ〈h−1 e g, uch − u γ h〉EBh + 〈g,∇h,0(uch − u γ h) · n〉EBh = 1 2 ( ∇+ h,gu,∇ + h,0(uch − u γ h) ) Th + 1 2 ( ∇−h,gu,∇ − h,0(uch − u γ h) ) Th − (f, uch − u γ h)Th = ( ∇h,gu,∇h,0(uch − u γ h) ) Th − (f, uch − u γ h)Th = ( ∇h,gu,∇h,0(uch − u γ h) ) Th − (−∆u, uch − u γ h)Th = ( ∇h,gu,∇h,0(uch − u γ h) ) Th − (∇h,gu,∇h,0(uch − u γ h))Th + 〈 { ∇h,gu−∇u } · n, [[uch − u γ h]]〉Eh = 〈 { ∇h,gu−∇u } · n, [[uch − u γ h]]〉Eh . Applying the Cauchy-Schwarz inequality, the trace theorem with scaling, and Lemma 2.4, we have Bh,γ(u− uγh, u c h − u γ h) = 〈 { ∇h,gu−∇u } · n, [[uch − u γ h]]〉Eh ≤ ( ∑ e∈Eh ‖h1/2 e { ∇h,gu−∇u } ‖2L2(e) )1/2( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 ≤ Chr‖u‖Hr+1(Ω) ( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 . (4.10) Combining (4.7), (4.8), and (4.10), we have γ ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) ≤ C ( |||uch − u|||h,1 + Chr‖u‖Hr+1(Ω) )( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 . Dividing both sides above by (∑ e∈Eh ‖h −1 e [[uch − u γ h]]‖2L2(e) )1/2 and γ gives us( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 ≤ C γ ( |||uch − u|||h,1 + Chr‖u‖Hr+1(Ω) ) . (4.11) EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 133 Lastly, recalling that uch − u ∈ H1 0 (Ω), from (2.1) and Lemma 2.4, we have |||uch − u|||2h,1 = 1 2 ‖∇+ h,0(uch − u)‖2L2(Th) + 1 2 ‖∇−h,0(uch − u)‖2L2(Th) + ∑ e∈Eh ‖h−1/2 e [[uch − u]]‖2L2(e) = 1 2 ‖∇+ h,0(uch − u)‖2L2(Th) + 1 2 ‖∇−h,0(uch − u)‖2L2(Th) ≤ 1 2 ‖∇uch −∇u‖2L2(Th) + 1 2 ‖∇uch −∇u‖2L2(Th) = ‖∇uch −∇u‖2L2(Ω) ≤ Ch2r‖u‖2Hr+1(Ω) . (4.12) By combining (4.6), (4.11), and (4.12), we obtain the result |||uch − u γ h|||h,1 ≤ ( ∑ e∈Eh ‖h−1/2 e [[uch − u γ h]]‖2L2(e) )1/2 ≤ C γ ( |||uch − u|||h,1 + Chr‖u‖Hr+1(Ω) ) ≤ C γ ( Chr‖u‖Hr+1(Ω) + Chr‖u‖Hr+1(Ω) ) = Chr γ ‖u‖Hr+1(Ω). � Since uch, u γ h ∈ V rh and ||| · |||h,1 is a norm on V rh , an immediate consequence of Theorem 4.4 is the following. Corollary 4.5. If u ∈ Hr+1(Ω) for r ≥ 1 is the solution to (1.1), uγh is the solution to the DWDG method (3.4), and uch is the solution to the CG method (4.1), then lim γ→∞ |||uγh − u c h|||h,1 = 0. 5. Numerical experiments In this section, we will verify the results of Theorem 4.4. We will also investigate the effect of the regularity of solutions to (1.2) on the relationship between the penalty parameter γ, the energy error, and the L2-error of the DWDG approxima- tion. 5.1. DWDG Converges to CG. To validate Theorem 4.4, we consider the ho- mogeneous problem −∆u = 2π2 sin(πx) sin(πy) in Ω = [0, 1]2, u = 0 on ∂Ω. (5.1) The solution to this equation is u = sin(πx) sin(πy). The approximations and error calculations were done on a fixed mesh with h = 1 32 for linear basis polynomials and h = 1 16 for quadratic and cubic basis polynomials. Based on the proof of Theorem 4.4 and [16, Lemma IV.4], we measured the error in • the energy norm |||uch − u γ h|||h,1, • the H1 semi-norm ‖∇(uch − u γ h)‖L2(Th), 134 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 • the jump error (∑ e∈Eh ‖h −1/2 e [[uch − u γ h]] ‖2L2(e) )1/2 . The error and rates using linear, quadratic, and cubic basis polynomials for the DWDG approximation with varying values of the penalty parameter γ can be found in Table 1. In all three error measurements, we are able to achieve a rate of 1 with respect to γ as γ → ∞. Note that the matrix corresponding to the bilinear form Bh,γ(·, ·) will require special treatment for larger values of γ due to the spectral radius increasing instep with the value of γ. Table 1. Energy error, H1 semi-norm error, jump error and their respective rates for (5.1). Poly γ |||uch − u γ h|||h,1 Rate ‖∇(uch − u γ h)‖L2(Th) Rate Jump Error Rate r = 1 1 1.3514e-02 — 1.1934e-02 — 9.4287e-03 — 10 7.2250e-03 0.2719 6.3040e-03 0.2772 4.9789e-03 0.2773 102 1.2874e-03 0.7491 1.1171e-03 0.7515 8.8320e-04 0.7511 103 1.3971e-04 0.9645 1.2114e-04 0.9648 9.5790e-05 0.9647 104 1.4094e-05 0.9962 1.2237e-05 0.9956 9.6608e-06 0.9963 105 1.4440e-06 0.9895 1.4043e-06 0.9402 9.6691e-07 0.9996 r = 2 1 1.0215e-03 — 9.9343e-04 — 5.4330e-04 — 10 6.9298e-04 0.1685 7.2228e-04 0.1384 3.4854e-04 0.1928 102 1.7932e-04 0.5871 2.0400e-04 0.5491 8.4431e-05 0.6157 103 2.1642e-05 0.9183 2.5125e-05 0.9095 1.0034e-05 0.9250 104 2.2109e-06 0.9907 2.5726e-06 0.9897 1.0232e-06 0.9915 105 2.2157e-07 0.9991 2.5771e-07 0.9992 1.0253e-07 0.9991 r = 3 1 1.0741e-05 — 1.4247e-05 — 4.9554e-06 — 10 7.7927e-06 0.1393 1.0155e-05 0.1470 3.5359e-06 0.1466 102 2.1274e-06 0.5638 2.7326e-06 0.5701 9.4950e-07 0.5710 103 2.5835e-07 0.9156 3.3106e-07 0.9167 1.1482e-07 0.9175 104 2.6405e-08 0.9905 3.3830e-08 0.9906 1.1730e-08 0.9907 105 2.6645e-09 0.9961 3.4807e-09 0.9876 1.1755e-09 0.9991 5.2. Relationship between γ and the L2-error. For our experiments, we will test the DWDG methods for approximating (1.2). We will find the L2-error on a fixed mesh with h = 1 16 . We will run the DWDG approximation with two sets of γ’s to get an idea of the behavior of the L2-error near γ = 0 as well as for large values of γ. To determine the behavior of the L2-error around γ = 0 we chose γ ∈ {−2,−1.8,−1.6,−1.4,−1.2, · · · , 10}. Since the constant C∗ has a complex deriva- tion, the exact value is generally unknown. Going only as low as γ = −2 should guarantee that the method is stable, and it gives us more information about the ef- fect of γ decreasing towards zero and becoming negative. To determine the behavior of the L2-error as γ →∞ we chose γ ∈ {10, 101.5, 102, 102.5, 103, 103.5, 104, 104.5, 105} to reflect the results found in Table 1. We measure the error of the DWDG approx- imation in the the L2 norm ‖u−uγh‖L2(Th). Since the choice of starting mesh, based on minimal angle, did not have a noticeable effect on the L2-error of the DWDG approximation [16], we have chosen to preform our experiments on a uniform criss- cross mesh. LDG methods are known to require γ > 0 for cris-cross meshes [13]. Thus, the choice represents a potentially challenging meshing strategy for DWDG methods without penalization. The numerical experiment was run on the problem −∆uα = f in Ω = [−1, 1]2, uα|∂Ω = g on ∂Ω, (5.2) EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 135 where α > 0 is a constant. The functions f and g are chosen so that the solution for (5.2) is uα(x, y) = { cos(π2 y) if x < 0, cos(π2 y) + xα if x ≥ 0. This solution belongs to Hα+ 1 2 (Ω) but does not belong to Hα+1/2+ε(Ω) for all ε > 0 [7]. To test the impact of the regularity of the solution on the choice of γ we will choose α ∈ {1.5, 2.5, 3.5, 4.5}. We have also plotted the CG approximation L2-error in each plot of Figure 1. Note that the CG L2-error was larger than the DWDG L2-errors in Figure 2, even though the CG L2-error was not plotted. The CG L2- error values are given instead of plotted for each combination of α and the smaller γ in Figure 2, as the relatively large difference in the L2-error for the DWDG and CG method prevented both errors to be plotted while visually retaining the shape of the DWDG’s L2-error curve. For the linear DWDG methods, we see a check-mark relationship between γ and the L2-error of the approximation when γ is close to zero (see Figure 2). When the solution u ∈ H2(Ω), we see an upward trend in the L2-error after γ = 5 creating a wide check-mark relationship. As the regularity of the solution increased, we see the bottom of the check-mark relationship move to the left, indicating that if the regularity was high enough, we would see a strictly increasing relationship between γ and the DWDG L2-error. This is a similar observation for the linear DWDG approximation in [16], where refining the mesh also moved the bottom of the check- mark relationship to the left. For the larger values of γ (see Figure 1), we see a strictly increasing relationship that tapers off as the DWDG L2-error approaches the CG L2-error. From this, we conjecture that there is a strictly increasing relationship between the L2-error and γ for the linear DWDG approximation for values of γ that are large enough or if the regularity of the solution is high enough. For the DWDG methods with quadratic and cubic basis polynomials, we see a strictly increasing relationship between γ and the L2-error of the DWDG approxi- mation when looking at the smaller range of γ (see Figure 2). For the larger range of γ (see Figure 1) a strictly increasing relationship was persistent for quadratic basis functions. For cubic basis functions, we see a strictly increasing trend until the largest values of γ when u ∈ H3(Ω) or u ∈ H4(Ω). Overall, we observe that as the regularity of the solution to (1.1) increased, we could see a strictly increasing trend between the L2-error and γ for γ ≥ 0. It is of interest to note that as γ → ∞, the spectral radius of the DWDG matrix from the bilinear form Bh,γ(·, ·) will scale with the choice of γ. Though we choose large values of γ to show the convergence of the DWDG approximation to the CG approximation in our numerical results, the DWDG approximation always had lower errors than the CG approximation. Furthermore, the best L2-errors were found near γ = 0. The choice of γ = 0 simplifies the approximation, and it appears to yield a more accurate approximation. Acknowledgements. This work was supported by the National Science Founda- tion under Grant DMS-2111059. Y. Zhang was supported by the National Science Foundation under Grant DMS-2111004. References [1] R. A. Adams, J. JF. Fournier; Sobolev spaces, Elsevier, 2003. 136 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 1 2 3 4 5 0.01185 0.0119 0.01195 0.012 0.01205 0.0121 0.01215 P1; Solution in H 2 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 8.614 8.6145 8.615 8.6155 8.616 8.6165 10 -3 P2; Solution in H 2 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 6.3672 6.3674 6.3676 6.3678 6.368 6.3682 6.3684 10 -3 P3; Solution in H 2 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 5 6 7 8 9 10 10 -4 P1; Solution in H 3 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 1.44 1.46 1.48 1.5 10 -5 P2; Solution in H 3 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 5.625 5.63 5.635 5.64 5.645 5.65 5.655 5.66 10 -6 P3; Solution in H 3 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 1 1.1 1.2 1.3 10 -3 P1; Solution in H 4 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 5.5 6 6.5 7 7.5 8 8.5 10 -6 P2; Solution in H 4 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 3 4 5 6 7 10 -8 P3; Solution in H 4 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 1.4 1.6 1.8 2 10 -3 P1; Solution in H 5 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 1.1 1.2 1.3 1.4 1.5 1.6 1.7 10 -5 P2; Solution in H 5 ( ) DWDG L2 Err CG L2 Error 1 2 3 4 5 0.6 0.8 1 1.2 1.4 1.6 10 -7 P3; Solution in H 5 ( ) DWDG L2 Err CG L2 Error Figure 1. DWDG L2-error vs large γ plots for linear (left col- umn), quadratic (middle column) and cubic (right column) basis polynomials on a uniform mesh with h = 1/16. [2] D. Arnold, F. Brezzi, B. Cockburn, L. Marini; Unified analysis of discontinuous Galerkin methods for elliptic problems, SINUM 39 (2002), no. 5, 1749–1779. [3] S. Brenner, R. Scott; The mathematical theory of finite element methods, vol. 15, Springer Science & Business Media, 2007. [4] S. C. Brenner; Poincaré–Friedrichs inequalities for piecewise H1 functions, SIAM J Numer Anal 41 (2003), no. 1, 306–19. [5] E. Burman, A. Quarteroni, B. Stamm; Interior penalty continuous and discontinuous finite element approximations of hyperbolic equations, J Sci Comput 43 (2010), no. 3, 293–312. [6] A. Cangiani, J. Chapman, E. Georgoulis, M. Jensen; On local super-penalization of interior penalty discontinuous galerkin methods, International Journal of Numerical Analysis and Modeling 11 (2013), 478–495. [7] Paul Castillo, Bernardo Cockburn, Ilaria Perugia, Dominik Schötzau; An a priori error anal- ysis of the local discontinuous galerkin method for elliptic problems, SIAM Journal on Nu- merical Analysis 38 (2000), no. 5, 1676–1706. [8] P. Ciarlet; The finite element method for elliptic problems, Classics in applied mathematics, no. 40, SIAM, Philadelphia, PA, 2002. EJDE-2018/CONF/26 DWDG PENALTY PARAMETER STUDY 137 -2 0 2 4 6 8 10 0.011888 0.01189 0.011892 0.011894 0.011896 P1; Solution in H 2 ( ) DWDG L2 Err CG L2 Err = 1.215e-02 -2 0 2 4 6 8 10 8.61411 8.61412 8.61413 8.61414 8.61415 8.61416 8.61417 8.61418 10 -3 P2; Solution in H 2 ( ) DWDG L2 Err CG L2 Err = 8.616e-03 -2 0 2 4 6 8 10 6.3672736 6.3672737 6.3672738 6.3672739 6.367274 6.3672741 6.3672742 10 -3 P3; Solution in H 2 ( ) DWDG L2 Err CG L2 Err = 6.368e-03 -2 0 2 4 6 8 10 2 3 4 5 6 10 -4 P1; Solution in H 3 ( ) DWDG L2 Err CG L2 Err = 9.347e-04 -2 0 2 4 6 8 10 1.425 1.43 1.435 1.44 1.445 1.45 10 -5 P2; Solution in H 3 ( ) DWDG L2 Err CG L2 Err = 1.506e-05 -2 0 2 4 6 8 10 5.6278 5.6279 5.628 5.6281 10 -6 P3; Solution in H 3 ( ) DWDG L2 Err CG L2 Err = 5.655e-06 -2 0 2 4 6 8 10 3 4 5 6 7 8 9 10 -4 P1; Solution in H 4 ( ) DWDG L2 Err CG L2 Err = 1.388e-03 -2 0 2 4 6 8 10 4 4.5 5 5.5 6 10 -6 P2; Solution in H 4 ( ) DWDG L2 Err CG L2 Err = 8.367e-06 -2 0 2 4 6 8 10 2.95 3 3.05 3.1 3.15 3.2 10 -8 P3; Solution in H 4 ( ) DWDG L2 Err CG L2 Err = 6.737e-08 -2 0 2 4 6 8 10 4 6 8 10 12 14 10 -4 P1; Solution in H 5 ( ) DWDG L2 Err CG L2 Err = 2.027e-03 -2 0 2 4 6 8 10 8 9 10 11 10 -6 P2; Solution in H 5 ( ) DWDG L2 Err CG L2 Err = 1.669e-05 -2 0 2 4 6 8 10 5.6 5.8 6 6.2 6.4 6.6 10 -8 P3; Solution in H 5 ( ) DWDG L2 Err CG L2 Err = 1.484e-07 Figure 2. DWDG L2-error vs small γ plots for linear (left col- umn), quadratic (middle column) and cubic (right column) basis polynomials on a uniform mesh with h = 1/16. [9] B. Cockburn, B. Dong; An analysis of the minimal dissipation local discontinuous Galerkin method for convection-diffusion problems, J Sci Comput 32 (2007), no. 2, 233–262. [10] W. Feng, T. L. Lewis, S. M. Wise; Discontinuous Galerkin derivative operators with applica- tions to second-order elliptic problems and stability, Math Meth Appl Sci 38 (2015), no. 18, 5160–5182. [11] X. Feng, T. Lewis, M. Neilan; Discontinuous Galerkin finite element differential calculus and applications to numerical solutions of linear and nonlinear partial differential equations, J Comput Appl Math 299 (2016), 68–91. [12] M. Larson, A. Niklasson; Conservation properties for the continuous and discontinuous Galerkin methods, Chalmers Finite Element Center Preprint 8, 2000. [13] T. Lewis; Distributional derivatives and stability of discontinuous Galerkin finite element approximation methods, Electron J Differ Eq (2016), 59–76. [14] T. Lewis M. Neilan; Convergence analysis of a symmetric dual-wind discontinuous Galerkin method: Convergence analysis of DWDG, J Sci Comput 59 (2014), no. 3, 602–625. [15] T. Lewis, A. Rapp, Y. Zhang; Convergence analysis of symmetric dual-wind discontinuous galerkin approximation methods for the obstacle problem, J Math Anal Appl (2020), 123840. 138 T. LEWIS, A. RAPP, Y. ZHANG EJDE-2022/CONF/26 [16] A. Rapp; Symmetric dual-wind discontinuous Galerkin methods for elliptic variational in- equalities, Ph.D. thesis, The University of North Carolina at Greensboro, 2020. Thomas Lewis Department of Mathematics and Statistics, The University of North Carolina at Greens- boro, Greensboro, NC 27402, USA Email address: tllewis3@uncg.edu Aaron Rapp Department of Mathematical Sciences, The University of the Virgin Islands, Charlotte Amalie West, St. Thomas, 00820, United States Virgin Islands Email address: aaron.rapp@uvi.edu Yi Zhang Department of Mathematics and Statistics, The University of North Carolina at Greens- boro, Greensboro, NC 27402, USA Email address: y zhang7@uncg.edu 1. Introduction 2. Notation and DG differential calculus 2.1. Notation 2.2. Discontinuous Galerkin calculus 3. Dual-wind discontinuous Galerkin approximation method for Poisson's equation 4. DWDG converges to CG 5. Numerical experiments 5.1. DWDG Converges to CG 5.2. Relationship between and the L2-error Acknowledgements References