Interior Penalty Discontinuous Galerkin Methods for Second Order Linear Non-Divergence Form Elliptic PDEsThe work of the first and third authors was partial supported by the NSF through grant DMS-1318486 and the work of the second author was partially supported by the NSF grant DMS-1417980 and the Alfred Sloan Foundation.

# Interior Penalty Discontinuous Galerkin Methods for Second Order Linear Non-Divergence Form Elliptic PDEs††thanks: The work of the first and third authors was partial supported by the NSF through grant DMS-1318486 and the work of the second author was partially supported by the NSF grant DMS-1417980 and the Alfred Sloan Foundation.

Xiaobing Feng Department of Mathematics, University of Tennessee, Knoxville, TN 37996 (xfeng@math.utk.edu).    Michael Neilan Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260 (neilan@pitt.edu).    Stefan Schnake Department of Mathematics, University of Tennessee, Knoxville, TN 37996 (schnake@math.utk.edu).
###### Abstract

This paper develops interior penalty discontinuous Galerkin (IP-DG) methods to approximate strong solutions of second order linear elliptic partial differential equations (PDEs) in non-divergence form with continuous coefficients. The proposed IP-DG methods are closely related to the IP-DG methods for advection-diffusion equations, and they are easy to implement on existing standard IP-DG software platforms. It is proved that the proposed IP-DG methods have unique solutions and converge with optimal rate to the strong solution in a discrete -norm. The crux of the analysis is to establish a DG discrete counterpart of the Calderon-Zygmund estimate and to adapt a freezing coefficient technique used for the PDE analysis at the discrete level. As a byproduct of our analysis, we also establish broken -norm error estimates for IP-DG approximations of constant coefficient elliptic PDEs. Numerical experiments are provided to gauge the performance of the proposed IP-DG methods and to validate the theoretical convergence results.

{AMS}

65N30, 65N12, 35J25

## 1 Introduction

In this paper we develop interior penalty discontinuous Galerkin methods for approximating the strong solution to the following second order linear elliptic PDE in non-divergence form:

 (1a) Lu(x):=−A(x):D2u(x) =f(x) x∈Ω, (1b) u(x) =0 x∈∂Ω,

where is an open, bounded domain with boundary , with , and is positive definite in .

Non-divergence form elliptic PDEs arrive naturally from many applications such as stochastic optimal control and game theory [12]; they are also encountered in the linearization of fully nonlinear PDEs such as Monge-Ampère-type equations [5]. If is differentiable, then it is easy to check that equation (1a) can be rewritten as a diffusion-convection equation with as the diffusion coefficient and () as the convection coefficient. However, if , then this formulation is not possible since () does not exist as a function, but rather only as a measure. We recall that in the literature there are three well-established PDE theories for non-divergence form elliptic PDEs depending on the smoothness of and (as well as the smoothness of the boundary which we assume sufficiently smooth here to ease the presentation). The first theory is the classical solution (or Schauder’s) theory [13, Chapter 6] which seeks solutions in the Hölder space for when and . The second one is the (strong) solution theory [13, Chapter 9] which seeks solutions in the Sobolev space for that satisfy the PDE almost everywhere in under the assumptions that and . If the coefficient matrix satisfies the Cordes condition and if the domain is convex, then the notion of strong solutions may be extended to the case [26, 16]. The third one is the viscosity (weak) solution theory [8] which seeks solutions in (or even in , the space of bounded functions in ) that satisfy the PDE in the viscosity sense under the assumptions that and . We note that all the three solution concepts and PDE theories are non-variational; this is a main difference between divergence form PDEs [13, Chapter 8] and non-divergence form PDEs. We also note that in the case of the viscosity solution theory, the uniqueness and regularity of solutions had been the main focus of the study for second order non-divergence form PDEs (see [8, 26, 17, 23] and the references therein).

In contrast to the advances of the PDE analysis, almost no progress on numerical methods and numerical analysis was achieved until very recently for second order elliptic PDEs in non-divergence form with non-differentiable coefficient matrix (cf. [9, 11, 15, 18, 20, 26, 28]). The main difficulty is caused by the non-divergence structure of the PDEs which prevents any straightforward application of Galerkin-type numerical methodologies such as finite element methods, discontinuous Galerkin methods and spectral methods. Moreover, the non-variational nature of the strong and viscosity solution concepts make convergence analysis and/or error estimates of any convergent numerical methods very delicate and difficult. Nonstandard numerical techniques are often required to do the job (cf. [11, 18, 20, 26, 28]).

The primary goal of this paper is to develop convergent interior penalty discontinuous Galerkin (IP-DG) methods for approximating the strong solution of problem (1) under the assumption that and its solution satisfies the Calderon-Zygmund estimate [13, Chapter 9]. This paper can be viewed as a DG counterpart of [11] where non-standard convergent finite element methods were proposed and analyzed for approximating the strong solution of (1). The reason for undertaking such an extension is twofold. First, we intend to take advantage of some of the features of DG methods such as simplicity and ease of computation, and flexibility of mesh and ease for adaptivity, to design better numerical methods for problem (1a)–(1b); in the meantime, we develop new numerical analysis techniques and machineries for non-divergence form PDEs. Second, since the jumps of the normal derivatives across element edges and mesh-dependent bilinear forms are already used in the finite element methods of [11], it is natural to use totally discontinuous piecewise polynomial approximation spaces to explore the full potential of the interior penalty technique and idea. It should be noted that although the generalization of the finite element formulations of [11] to the DG case is not very hard, the DG convergence analysis is more involved because of the extra difficulty caused by the non-conformity of the DG finite element space. The crux of the analysis of this paper is to establish a discrete Calderon-Zygmund estimate for the proposed IP-DG methods and to adapt a freezing coefficient technique used in the PDE analysis at the discrete level. Moreover, in order to prove the desired discrete Calderon-Zygmund estimate, we need to establish the stability and error estimates for the IP-DG approximations of constant coefficient elliptic PDEs. Such estimates seem to be new and have independent interest.

The remainder of this paper is organized as follows. In Section 2 we provide the notation and a collection of preliminary estimates, in particular, some properties of DG functions are either cited or proved. Section 3 analyzes the IP-DG methods for constant coefficient elliptic PDEs. The stability and error estimates are derived, which in turn lead to global stability estimates. The latter estimate can be regarded as a discrete Calderon-Zygmund estimate for the proposed IP-DG methods. To the best of our knowledge, such an estimate is new and of independent interest. Section 4 is devoted to the formulation, stability and convergence analysis and error estimate for the proposed IP-DG methods. Here, using the continuity of the coefficient matrix , the freezing coefficient technique is adapted to establish local stability (or a left-side inf-sup condition) for the IP-DG discrete operators, which together with a covering argument leads to a global Gärding-type inequality for the formal adjoint operators of the IP-DG operators. Next, using a duality argument we obtain a global left-side inf-sup condition for the IP-DG discrete operators. Finally, by employing another duality argument, we derive the desired discrete Calderon-Zygmund estimate for the DG discrete operators. Once this stability estimate is shown, well-posedness and convergence of the IP-DG methods follow easily. In Section 5 we present a number of numerical experiments to verify our theoretical results and to gauge the performance of the proposed IP-DG methods, even for the case which is not covered by our convergence theory.

## 2 Preliminary Results

### 2.1 Notation

Let be an open and bounded domain in . For a subdomain of with boundary , let and for and denote the standard Lebesgue and Sobolev spaces respectively, and be the closure of in . Let be the inner product on and . To improve the readability of the paper, we adopt the convention that stands for for some which does not depend on any discretization parameters.

Let be a shape-regular and conforming triangulation of . Let and denote respectively the sets of all interior and boundary edges/faces of , and set . We introduce the broken Sobolev spaces

 Ws,p(Th) :=∏T∈ThWs,p(T), Lp(Th):=W0,p(Th), Ws,ph(D) :=Ws,p(Th)∣∣D, Lph(D):=Lp(Th)∣∣D.

For any interior edge/face we define the jump and average of a scalar or vector valued function as

 [v]∣∣e:=v+−v−,{v}∣∣e:=12(v++v−),

where . On a boundary edge/face with , we set . For any we use to denote the unit outward normal vector pointing in the direction of the element with the smaller global index. For we set to be the outward normal to restricted to . The standard DG finite element space is defined as

 Vh :={vh∈W2,p(Th);vh∣∣T∈Pk(T)∀T∈Th},

where denotes the set of all polynomials of degree less than or equal to on . We also introduce for any

 Vh(D):={v∈Vh;v∣∣Ω∖¯¯¯¯D≡0}.

Note that is nontrivial provided that there exists an inscribed ball with radius such that . We also adopt the convention .

For each , let be constant on . We define the following mesh-dependent norms on and :

 (2) ∥v∥W2,ph(D) :=∥D2hv∥Lp(D)+(∑e∈EIhh1−pe∥∥|[∇v]|∥∥pLp(e∩¯D))1p +(∑e∈Ehγpeh1−2pe∥[v]∥pLp(e∩¯D))1p, (3) ∥v∥W1,ph(D) :=∥∇hv∥Lp(D)+(∑e∈Ehγpeh1−pe∥[v]∥pLp(e∩¯D))1p +(∑e∈Ehhe∥{∇v}∥pLp(e∩¯D))1p,

where and denote the piecewise gradient and Hessian of .

In addition, we define the discrete -norm and -norm as follows:

 (4) ∥q∥W−2,ph(D) :=sup0≠vh∈Vh(q,vh)D∥vh∥W2,p′h(D), (5) ∥q∥W−1,ph(D) :=sup0≠v∈W1,p′h(D)(q,v)D∥v∥W1,p′h(D),

where . Finally, for any domain and any , we introduce the following mesh-dependent semi-norm

 (6)

It can be proved that (cf. [11])

 (7) ∥wh∥Lp(Ω)≲∥wh∥Lph(Ω)∀wh∈Vh.

### 2.2 Properties of the DG space Vh

In this subsection we collect some technical lemmas that cover the basic properties of functions in the DG space . These facts will be used many times in the later sections. We first state the standard trace inequalities for broken Sobolev functions, a proof of this lemma can be found in [3]. {lemma} For any and , there holds

 (8) ∥v∥pLp(∂T)≲(hp−1T∥∇v∥pLp(T)+h−1T∥v∥pLp(T))∀v∈W1,p(T).

Therefore by scaling we have

 (9) ∑e∈EIhhe∥v∥pLp(e∩¯¯¯¯D)≲{∥v∥pLp(D)∀v∈Vh(D),∥v∥pLp(D)+hp∥∇v∥pLp(D)∀v∈W2,ph(D).

Next, we prove an inverse inequality between the -norm and the -norm. {lemma} For any , , there holds for

 (10) ∥vh∥W2,ph(D) ≲h−1∥vh∥W1,ph(Dh),

where

 (11) Dh={x∈Ω;dist(x,D)≤h}.
{proof}

To show (10), we use (2), (8), and standard inverse estimates [3] to obtain

 ∥vh∥W2,ph(D) ≲∥D2hvh∥Lp(D)+(∑e∈Ehγpeh1−2pe∥[v]∥pLp(e∩¯D))1p +∑T∈ThT⊂Dh(h1−pT(hp−1T∥D2vh∥pLp(T)+h−1T∥∇vh∥pLp(T)))1p ≲∥D2hvh∥Lp(D)+(∑e∈Ehγpeh−peh1−pe∥[v]∥pLp(e∩¯D))1p +∑T∈ThT⊂Dh(∥D2vh∥pLp(T)+h−pT∥∇vh∥pLp(T))1p ≲h−1∥vh∥W1,ph(Dh)+h−1(∑e∈EIhγpeh1−pe∥[v]∥pLp(e∩¯D))1p ≲h−1∥vh∥W1,ph(Dh).

We also prove an inverse inequality between the -norm and the -norm.

{lemma}

Let . For any and subdomain we have

 (12) ∥vh∥Lp(D)≲h−1∥vh∥W−1,ph(D).
{proof}

Using the relation (7) and the definition of , we find that

 ∥vh∥Lp(D)≤∥vh∥Lp(Ω)≲∥vh∥Lph(Ω)=sup0≠wh∈Vh(vh,wh)D∥wh∥Lp′(Ω)∀vh∈Vh(D).

Therefore, by the standard inverse estimate and noting that , we obtain

 ∥vh∥Lp(D) ≲h−1sup0≠wh∈Vh(vh,wh)D∥wh∥W1,p′h(D)≤h−1sup0≠w∈W1,p′h(D)(vh,w)D∥w∥W1,p′h(D)=∥vh∥W−1,ph(D).

The proof is complete.

The following lemma shows that the broken Sobolev norms are controlled by their corresponding Sobolev norms. {lemma} For any there holds the following inequality:

 ∥φ∥W2,ph(Ω) ≤∥φ∥W2,p(Ω)∀φ∈W2,p(Ω)∩W1,p0(Ω).
{proof}

Since the inequality holds for all and it can be extended to all by a density argument.

The next lemma establishes a Poincaré-Friedrichs’ inequality for DG functions.

{lemma}

Let such that and . Then for any there holds the following inequalities:

 (13) ∥vh∥Lp(D) ≲diam(D)∥vh∥W1,ph(D), (14) ∥vh∥W1,ph(D) ≲diam(D)∥vh∥W2,ph(D).
{proof}

Let denote the generalized Hsiegh–Clough–Tochner space [10], and let be the reconstruction operator constructed in [14]. The arguments given in [14] show that, for ,

 (15) Ehvh∈H20(Dh), ∥vh−Ehvh∥Lp(Ω)≲h∥vh∥W1,ph(D), ∥vh−Ehvh∥Wm,ph(Ω)≲hs−m∥vh∥Ws,ph(D),1≤m≤s≤2,

where is the same as in Lemma 2.2. Therefore, by the triangle inequality, the Poincarè-Friedrichs inequality, and the assumption ,

 ∥vh∥Lp(D) ≤∥Ehvh∥Lp(D)+∥vh−Ehvh∥Lp(Ω) ≲diam(D)∥Ehvh∥W1,p(Dh)+h∥vh∥W1,ph(D)≲diam(D)∥vh∥W1,ph(D).

Likewise, we find

 ∥vh∥W1,ph(D) ≤∥Ehvh∥W1,p(Dh)+∥vh−Ehvh∥W1,ph(Ω), ≲diam(D)∥Ehvh∥W2,p(Dh)+h∥vh∥W2,ph(Ω)≲diam(D)∥vh∥W2,ph(D).

The proof is complete.

Next we establish a discrete Sobolev interpolation estimate for DG functions. {lemma} Let . For all we have

 (16) ∥vh∥2W1,ph(Ω)≲∥vh∥Lp(Ω)∥vh∥W2,ph(Ω).
{proof}

Let be the enriching operator in the proof of Lemma 2.2. By the triangle inequality and scaling we find

 (17) ∥vh∥2W1,ph(Ω)≲∥vh−Ehvh∥2W1,ph(Ω)+∥Ehvh∥2W1,p(Ω).

Since we can apply the Gagliardo–Nirenberg estimate [4] to get

 ∥Ehvh∥2W1,p(Ω)≲∥Ehvh∥W2,p(Ω)∥Evh∥Lp(Ω).

Applying estimates (15), we conclude that

 (18) ∥Ehvh∥2W1,p(Ω)≲∥vh∥W2,ph(Ω)∥vh∥Lp(Ω).

Likewise, by (15) and an inverse estimate,

 (19) ∥vh−Ehvh∥2W1,ph(Ω)≲h2∥vh∥2W2,ph(Ω)≲∥vh∥Lp(Ω)∥vh∥W2,ph(Ω).

Combining (17)–(19) completes the proof.

Next we prove some local super approximation estimates for the DG nodal interpolation in various discrete norms. The derivation of the lemma is standard (cf. [19]); for completeness we provide the proof in the Appendix. {lemma} Let denote the nodal interpolation operator, and with for . Then for any and we have

 (20) ∥ηvh−Ih(ηvh)∥Lp(D) ≲hd∥vh∥Lp(Dh), (21) h∥∇h(ηvh−Ih(ηvh))∥Lp(D) ≲hd∥vh∥Lp(Dh), (22) h2∥D2h(ηvh−Ih(ηvh))∥Lp(D) ≲hd∥vh∥Lp(Dh), (23) ∥ηvh−Ih(ηvh)∥W2,ph(D) ≲1d2(∥vh∥Lp(Dh)+∥∇hvh∥Lp(Dh)),

where is the same as in Lemma 2.2. Moreover, there holds

 (24) ∥ηvh−Ih(ηvh)∥W2,ph(D)≲hd3∥vh∥W2,ph(Dh)

if the polynomial degree .

## 3 DG discrete W1,p and Calderon-Zygmund estimates for PDEs with constant coefficients

In this section we consider the constant coefficient case, that is, on . We define three interior-penalty discontinuous Galerkin discretizations to the PDE operator and extend their domains to the broken Sobolev space . Our goal in this subsection is to prove global stability estimates for which will be crucially used in the next section. The final global stability estimate given in Theorem 3.2 can be regarded as a DG discrete Calderon-Zygmund estimate for .

Let be a constant, positive-definite matrix in and define

 (25) L0w:=−A0:D2w=−∇⋅(A0∇w).

From this we gather the standard PDE weak form:

 (26) a0(w,v):=∫ΩA0∇w⋅∇vdx∀w,v∈H10(Ω).

The Lax-Milgram theorem [3] yields the existence and boundedness of . Moreover if we have from Calderon-Zygmund theory [13] that exists and

 ∥L−10φ∥W2,p(Ω)≲∥φ∥Lp(Ω)∀φ∈Lp(Ω),

and therefore

 ∥w∥W2,p(Ω)≲∥L0w∥Lp(Ω)∀w∈W2,p(Ω)∩W1,p0(Ω).

Define by

 (27) (Lε0,hwh ,vh):=aε0,h(wh,vh)∀vh,wh∈Vh,

where the IP-DG bilinear form is defined by

 (28) a0,h(wh,vh):= ∫ΩA0∇hwh⋅∇hvhdx−∑e∈Eh∫e{A0∇wh⋅νe}[vh]dS −ε∑e∈Eh∫e{A0∇vh⋅νe}[wh]dS+∑e∈Eh∫eγehe[wh][vh]dS,

and is a penalization parameter. The parameter choices give respectively the SIP-DG, IIP- DG, and NIP-DG formulations. For the sake of clarity and readability we shall assume for the rest of the paper that may be either , , or unless otherwise stated.

We recall the following well-known DG integration by parts formula:

 (29) ∫Ωτ⋅∇hvdx=−∫Ω (∇h⋅τ)vdx+∑e∈EIh∫e[τ⋅νe]{v}dS+∑e∈Eh∫e{τ⋅νe}[v]dS,

which holds for any piecewise scalar-valued function and vector-valued function . Applying (29) to the first term on the right-hand side of (28) yields

 (30) aε0,h(wh,vh)= −∫Ω(A0:D2hwh)vhdx+∑e∈EIh∫e[A0∇wh⋅νe]{vh}dS −ε∑e∈Eh∫e{A0∇vh⋅νe}[wh]dS+∑e∈Eh∫eγehe[wh][vh]dS

for any . By Hölder’s inequality, it is easy to check that the above new form of is also well-defined on with . As a result, this new form enables us to extend the domain of to and .

### 3.1 DG discrete W1,p error estimates

From the standard IP-DG theory [22], there exists depending only on the shape regularity of the mesh and on such that is invertible on provided ; in the non–symmetric case , can be any positive number. Moreover, if and satisfy

 (31) aε0,h(w−wh,vh)=0∀vh∈Vh,

then the quasi-optimal error estimate

 (32) ∥w−wh∥W1,2h(Ω)≲infvh∈Vh∥w−vh∥W1,2h(Ω)

is satisfied. The goal of this subsection is to generalize this result to general exponent for the SIP-DG method. In particular, we have

{theorem}

Suppose and satisfy (31) with . Then there holds

 (33) ∥w−wh∥W1,ph(Ω)≲h|logh|t∥w∥W2,p(Ω),

where if and if .

To prove Theorem 3.1 we introduce some notation given in [6] (also see [25]). For given , we define the weight function as

 (34) σz(x)=h|x−z|+h.

For , and , we define the following weighted norms

 ∥v∥Lp(D),z,s =(∫D∣∣σsz(x)v(x)∣∣pdx)1/p, ∥v∥W1,p(D),z,s =∥v∥Lp(D),z,s+∥∇hv∥Lp(D),z,s, (35) ∥v∥W1,ph(D),z,s =∥v∥W1,p(D),z,s+(∑e∈Ehh1−pe∥∥σsz[v]∥∥pLp(e∩¯¯¯¯D))1/p +(∑e∈Ehhe∥∥σsz{∇hv}∥∥pLp(e∩¯¯¯¯D))1/p.

The weighted norms in the case are defined analogously.

The derivation of error estimates of DG approximations is based on the work [6], where localized pointwise estimates of DG approximations are obtained. There it was shown that if and satisfy (31) with , then

 (36) |∇(w−wh)(z)|≲infvh∈Vh∥w−vh∥W1,∞h(Ω),z,s0≤s

for all . Similar to pointwise estimates of finite element approximations (e.g., [25, 3]), the ingredients to prove (36) include duality arguments and DG approximation estimates of regularized Green functions in a weighted (discrete) -norm. These results are rather technical and involve dyadic decompositions of , local DG error estimates, and Green function estimates.

Here, we follow a similar argument to derive estimates; the main difference being that we derive DG approximation estimates of regularized Green functions in a weighted (discrete) -norm with (cf. Lemma 3.1). Using these estimates and applying similar arguments in [25, 6] then yield the estimate

 |∇(w−wh)(z)|p≲h−ninfvh∈Vh∥w−vh∥pW1,ph(Ω),z,s

for certain values of . Integrating this expression with respect to and applying Fubini’s theorem (cf. Lemma 3.1) then yields estimates of the piecewise gradient error.

Unfortunately, the strategy just described does not immediately give us estimates for the terms appearing in the -norm. To bypass this difficulty, we first use the trace inequality

 ∑e∈Ehh1−pe∥[w−wh]∥pLp(e)≲∥∇h(w−wh)∥pLp(Ω)+h−p∥w−