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 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.

## 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.

## 1Introduction

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:

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 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] which seeks solutions in the Hölder space for when and . The second one is the (strong) solution theory [13] 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]. 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] 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] 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]). 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]).

The primary goal of this paper is to develop convergent interior penalty discontinuous Galerkin (IP-DG) methods for approximating the strong solution of problem under the assumption that and its solution satisfies the Calderon-Zygmund estimate [13]. 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 . 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 –; 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 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.

## 2Preliminary Results

### 2.1Notation

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

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

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

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

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 :

where and denote the piecewise gradient and Hessian of .

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

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

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

### 2.2Properties 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].

Next, we prove an inverse inequality between the -norm and the -norm.

To show , we use , , and standard inverse estimates [3] to obtain

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

Using the relation and the definition of , we find that

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

The proof is complete.

The following lemma shows that the broken Sobolev norms are controlled by their corresponding Sobolev norms.

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.

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 ,

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

Likewise, we find

The proof is complete.

Next we establish a discrete Sobolev interpolation estimate for DG functions.

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

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

Applying estimates , we conclude that

Likewise, by and an inverse estimate,

Combining – 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.

## 3DG 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 ? can be regarded as a DG discrete Calderon-Zygmund estimate for .

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

From this we gather the standard PDE weak form:

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

and therefore

Define by

where the IP-DG bilinear form is defined by

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:

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

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.1DG 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

then the quasi-optimal error estimate

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

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

For , and , we define the following weighted norms

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 with , then

for all . Similar to pointwise estimates of finite element approximations (e.g., [25]), the ingredients to prove 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 ?). Using these estimates and applying similar arguments in [25] then yield the estimate

for certain values of . Integrating this expression with respect to and applying Fubini’s theorem (cf. Lemma ?) 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

and then derive estimates for . We note that the standard duality argument to derive estimates yields

which is of little benefit. Rather, our strategy is to modify the arguments given in [6] and estimate in terms of (cf. Lemma ?) and then apply Fubini’s theorem. We note that it is due to this term that the factor appears in Theorem ?.

(i) Let and extend to by zero. Denote by the ball of radius and center . Then by a change of variables and interchanging integrals, we find

This proves .

(ii) To prove we again extend to by zero and make a change of variables to obtain

where is a dilation of . Therefore,

For , there holds

Combining this identity with yields the inequality

If , then we find by a direct calculation that

and therefore by ,

Next, by trace inequalities given in Lemma ?, we have

Noting that

we obtain

Likewise we have

Combining – yields

Finally applying the identities – to yields the desired result –. The proof is complete.

Step 1: Set-up. By the triangle inequality, an inverse estimate, and Hölder’s inequality we obtain

Therefore by standard approximation theory, and since on , we have

Replacing and by and , respectively, yields

Next, define by

and let be the regularized Green’s function satisfying

Setting to be the DG approximation of , i.e., , and , we have by Galerkin orthogonality and the continuity of the bilinear form,

Consequently, by , we have

Thus, the proof will be completed once it is shown that . This result is derived in the following steps.

Step 2: Dyadic decomposition of . To estimate