C^{0} discontinuous Galerkin finite element methods for second order linear elliptic partial differential equations in non-divergence formThis work was partial supported by the NSF through grants DMS-1016173 and DMS-1318486 (Feng), and DMS-1417980 (Neilan), and the Alfred Sloan Foundation (Neilan).

# C0 discontinuous Galerkin finite element methods for second order linear elliptic partial differential equations in non-divergence form††thanks: This work was partial supported by the NSF through grants DMS-1016173 and DMS-1318486 (Feng), and DMS-1417980 (Neilan), and the Alfred Sloan Foundation (Neilan).

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

This paper is concerned with finite element approximations of strong solutions of second-order linear elliptic partial differential equations (PDEs) in non-divergence form with continuous coefficients. A nonstandard (primal) finite element method, which uses finite-dimensional subspaces consisting globally continuous piecewise polynomial functions, is proposed and analyzed. The main novelty of the finite element method is to introduce an interior penalty term, which penalizes the jump of the flux across the interior element edges/faces, to augment a nonsymmetric piecewise defined and PDE–induced bilinear form. Existence, uniqueness and error estimate in a discrete energy norm are proved for the proposed finite element method. This is achieved by establishing a discrete Calderon–Zygmund–type estimate and mimicking strong solution PDE techniques at the discrete level. Numerical experiments are provided to test the performance of proposed finite element method and to validate the convergence theory.

## 1 Introduction

In this paper we consider finite element approximations of the following linear elliptic PDE in non-divergence form:

 (1.1a) Lu:=−A:D2u =f in Ω, (1.1b) u =0 on ∂Ω.

Here, is an open bounded domain with boundary , is given, and is a positive definite matrix on , but not necessarily differentiable. Problems such as (1.1) arise in fully nonlinear elliptic Hamilton-Jacobi-Bellman equations, a fundamental problem in the field of stochastic optimal control [12, 17]. In addition, elliptic PDEs in non-divergence form appear in the linearization and numerical methods of fully nonlinear second order PDEs [6, 11, 20].

Since is not smooth, the PDE (1.1a) cannot be written in divergence form, and therefore notions of weak solutions defined by variational principles are not applicable. Instead, the existence and uniqueness of solutions are generally sought in the classical or strong sense. In the former case, Schauder theory states the existence of a unique solution to (1.1) provided the coefficient matrix and source function are Hölder continuous, and if the boundary satisfies . In the latter case, the Calderon-Zygmund theory states the existence and uniqueness of satisfying (1.1) almost everywhere provided , and . In addition, the existence of a strong solution to (1.1) in two-dimensions and on convex domains is proved in [18, 3, 2].

Due to their non-divergence structure, designing convergent numerical methods, in particular, Galerkin-type methods, for problem (1.1) has been proven to be difficult. Very few such results are known in the literature. Nevertheless, even problem (1.1) does not naturally fit within the standard Galerkin framework, several finite element methods have been recently proposed. In [19] the authors considered mixed finite element methods using Lagrange finite element spaces for problem (1.1). An analogous discontinuous Galerkin (DG) method was proposed in [9]. The convergence analysis of these methods for non-smooth remains open. A least-squares-type discontinuous Galerkin method for problem (1.1) with coefficients satisfying the Cordes condition was proposed and analyzed in [23]. Here, the authors established optimal order estimates in with respect to a -type norm.

The primary goal of this paper is to develop a structurally simple and computationally easy finite element method for problem (1.1). Our method is a primal method using Lagrange finite element spaces. The method is well defined for all polynomials degree greater than one and can be easily implemented on current finite element software. Moreover, our finite element method resembles interior penalty discontinuous Galerkin (DG) methods in its formulation and its bilinear form, which contains an interior penalty term penalizing the jumps of the fluxes across the element edges/faces. Hence, it is a DG finite element method. In addition, we prove that the proposed method is stable and converges with optimal order in a discrete -type norm on quasi-uniform meshes provided that the polynomial degree of the finite element space is greater than or equal to two.

While the formulation and implementation of the finite element method is relatively simple, the convergence analysis is quite involved, and it requires several nonstandard arguments and techniques. The overall strategy in the convergence analysis is to mimic, at the discrete level, the stability analysis of strong solutions of PDEs in non-divergence form (see [14, Section 9.5]). Namely, we exploit the fact that locally, the finite element discretization is a perturbation of a discrete elliptic operator in divergence form with constant coefficients; see Lemma 3.1. The first step of the stability argument is to establish a discrete Calderon-Zygmund-type estimate for the Lagrange finite element discretization of the elliptic operator in (1.1) with constant coefficients, which is equivalent to a global inf-sup condition for the discrete operator. The second step is to prove a local version of the global estimate and inf-sup condition. With these results in hand, local stability estimate for the proposed DG discretization of (1.1) can be easily obtained. We then glue these local stability estimates to obtain a global Gärding-type inequality. Finally, to circumvent the lack of a (discrete) maximum principle which is often used in the PDE analysis, we use a nonstandard duality argument to obtain a global inf-sup condition for the proposed DG discretization for problem (1.1). See Figure 1 for an outline of the convergence proof. Since the method is linear and consistent, the stability estimate naturally leads to the well-posedness of the method and the energy norm error estimate.

The organization of the paper is as follows. In Section 2 the notation is set, and some preliminary results are given. Discrete stability properties, including a discrete Calderon-Zygmund-type estimate, of finite element discretizations of PDEs with constant coefficients are established. In Section 3, we present the motivation and the formulation of our discontinuous finite element method for problem (1.1). Mimicking the PDE analysis from [14] at the discrete level, we prove a discrete stability estimate for the discretization operator. In addition, we derive an optimal order error estimate in a discrete -norm. Finally, in Section 4, we give several numerical experiments which test the performance of the proposed DG finite element method and validate the convergence theory.

## 2 Notation and preliminary results

### 2.1 Mesh and space notation

Let be an bounded open domain. We shall use to denote a generic subdomain of and denotes its boundary. denotes the standard Sobolev spaces for and , and to denote the subspace of consisting functions whose traces vanish up to order on . denotes the standard inner product on and . To avoid the proliferation of constants, we shall use the notation to represent the relation for some constant independent of mesh size .

Let be a quasi-uniform, simplical, and conforming triangulation of the domain . Denote by the set of interior edges in , the set of boundary edges in , and , the set of all edges in . We define the jump and average of a vector function on an interior edge as follows:

 [[w]]∣∣e =w+⋅n+∣∣e+w−⋅n−∣∣e, =12(w+⋅n+∣∣e−w−⋅n−∣∣e),

where and is the outward unit normal of .

For a normed linear space , we denote by its dual and the pairing between and . The Lagrange finite element space with respect to the triangulation is given by

 (2.1) Vh:={vh∈H10(Ω): vh|T∈Pk(T) ∀T∈Th},

where denotes the set of polynomials with total degree not exceeding on . We also define the piecewise Sobolev space with respect to the mesh

 Ws,p(Th):=∏T∈ThWs,p(T), W(p)h:=W2,p(Th)∩W1,p0(Ω), Lph(Th):=∏T∈ThLp(T), Ws,ph(D):=Ws,p(Th)∣∣D,Lph(D):=Lp(Th)∣∣D.

For a given subdomain , we also define and as the subspaces that vanish outside of by

 Vh(D): ={v∈Vh;v|Ω∖D=0},W(p)h(D):={v∈W(p)h; v|Ω∖D=0}.

We note that is non-trivial for .

Associated with , we define a semi-norm on for

 (2.2) ∥v∥W2,ph(D)

Here, denotes the piecewise Hessian matrix of , i.e., for all .

Let be the projection defined by

 (2.3) (Qhw,vh)=(w,vh)∀w∈L2(Ω),vh∈Vh.

It is well known that [8] satisfies for any

 (2.4) ∥Qhw∥Wm,p(Ω)≲∥w∥Wm,p(Ω)m=0,1;1

For any domain and any , we also introduce the following mesh-dependent semi-norm

 (2.5)

By (2.3), it is easy to see that is a norm on . Moreover by (2.4)

 (2.6) ∥wh∥Lp(Ω) =supv∈Lp′(Ω)(wh,v)∥v∥Lp′(Ω)=supv∈Lp′(Ω)(wh,Qhv)∥v∥Lp′(Ω) ≲supv∈Lp′(Ω)(wh,Qhv)∥Qhv∥Lp′(Ω)≤∥wh∥Lph(Ω)∀wh∈Vh.

### 2.2 Some basic properties of W(p)h functions

In this subsection we cite or prove some basic properties of the broken Sobolev functions in , and in particular, for piecewise polynomial functions. These results, which have independent interest in themselves, will be used repeatedly in the later sections. We begin with citing a familiar trace inequality followed by proving an inverse inequality.

###### Lemma 2.1 ([4])

For any , there holds

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

for any . Therefore by scaling, there holds

 (2.8) ∑e∈EIhhe∥v∥pLp(e∩¯D)≲{∥v∥pLp(D)∀v∈Vh(D),∥v∥pLp(D)+hp∥∇v∥pLp(D)∀v∈W(p)h(D).
###### Lemma 2.2

For any , and , there holds

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

where

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

By (2.2), (2.7) and inverse estimates [7, 4], we have

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

The next lemma states a very simple fact about the discrete norm on .

###### Lemma 2.3

For any , there holds

 (2.11) ∥φ∥W2,ph(Ω)≤∥φ∥W2,p(Ω)∀φ∈W2,p(Ω).

Next, we state some super-approximation results of the nodal interpolant with respect to the discrete semi-norm. The derivation of the following results is standard [21], but for completeness we give the proof in Appendix A

###### Lemma 2.4

Denote by the nodal interpolant onto . Let with for . Then for each with , there holds

 (2.12) hm∥ηvh−Ih(ηvh)∥Wm,p(D) ≲hd∥vh∥Lp(Dh)for m=0,1, (2.13) ∥ηvh−Ih(ηvh)∥W2,ph(D) ≲1d2∥vh∥W1,p(Dh),

Moreover, if , there holds

 (2.14) ∥ηvh−Ih(ηvh)∥W2,ph(D)≲hd3∥vh∥W2,p(Dh).

Here, satisfy the conditions in Lemma 2.2.

To conclude this subsection, we state and prove a discrete Sobolev interpolation estimate.

###### Lemma 2.5

There holds for all ,

 ∥∇w∥2Lp(Ω)≲∥w∥Lp(Ω)∥w∥W2,ph(Ω)∀w∈W(p)h.
{proof}

Writing and integrating by parts, we find

 ∥∇w∥pLp(Ω)=−∫Ω(|∇w|p−2Δw+(p−2)|∇w|p−4(D2hw∇w)⋅∇w)wdx +∑e∈EIh∫e[[|∇w|p−2∇w]]wds (2.15)

To bound the first term in (2.15) we apply Hölder’s inequality to obtain

 (2.16) ∫Ω|∇w|p−2|D2hw||w|dx ≤∥∥|∇w|p−2∥∥Lpp−2(Ω)∥D2w∥Lp(Ω)∥w∥Lp(Ω) =∥∇w∥p−2Lp(Ω)∥D2hw∥Lp(Ω)∥w∥Lp(Ω).

Likewise, by Lemma 2.1 we have

 (2.17) ∑e∈EIh∫e[[|∇w|p−2∇w]]wds ≤∑e∈EIh(h1pe∥∥[[∇w]]∥∥Lp(e))p−2(h1−ppe∥[[∇w]]∥Lp(e))(h1pe∥w∥Lp(e)) ≲∥∇w∥p−2Lp(Ω)∥w∥Lp(Ω)(∑e∈EIhh1−pe∥∥[[∇w]]∥∥pLp(e))1p ≲∥∇w∥p−2Lp(Ω)∥w∥Lp(Ω)∥w∥W2,ph(Ω).

Combining (2.15)–(2.17) we obtain the desired result. The proof is complete.

### 2.3 Stability estimates for auxiliary PDEs with constant coefficients

In this subsection, we consider a special case of (1.1a) when the coefficient matrix is a constant matrix, . We introduce the finite element approximation (or projection) of on and extend to the broken Sobolev space . We then establish some stability results for the operator . These stability results will play an important role in our convergence analysis of the proposed DG finite element method in Section 3.

Let be a positive definite matrix and set

 (2.18) L\rm 0w

The operator induces the following bilinear form:

 (2.19) a0(w,v):=⟨L\rm 0w,v⟩=∫ΩA\rm 0∇w⋅∇vdx∀w,v∈H10(Ω),

and the Lax-Milgram Theorem (cf. [10]) implies that exists and is bounded. Moreover, if , the Calderon-Zygmund theory (cf. [14, Chapter 9]) infers that exists and there holds

 (2.20) ∥L−1\rm 0φ∥W2,p(Ω)≲∥φ∥Lp(Ω)∀φ∈Lp(Ω).

Equivalently,

 (2.21) ∥w∥W2,p(Ω)≲∥∥L\rm 0% w∥∥Lp(Ω)∀w∈W2,p∩W1,p0(Ω).

The bilinear form naturally leads to a finite element approximation (or projection) of on , that is, we define the operator by

 (2.22) (L\rm 0,hwh,vh):=a0(wh,vh)∀vh,wh∈Vh.
###### Remark 2.1

When , the identity matrix, is exactly the finite element the discrete Laplacian that is, . By finite element theory [4], we know that is one-to-one and onto, and therefore exists.

Recall the following DG integration by parts formula:

 (2.23) ∫Ωτ⋅∇hvdx +∫e{{τ}}⋅[[v]]ds)+∑e∈EBh∫e(τ⋅ne)vds,

which holds for any piecewise scalar–valued function and vector–valued function . Here, is defined piecewise, i.e., for all . For any , using (2.23) with , we obtain

 (2.24)

We note that the above new form of is not well defined on . However, it is well defined on with . Hence, we can easily extend the domain of the operator to broken Sobolev space . Precisely, (abusing the notation) we define to be the operator induced by the bilinear form on , namely,

 (2.25) ⟨L\rm 0,hw,v⟩:=a0(w,v)∀w∈W(p)h, v∈W(p′)h.

A key ingredient in the convergence analysis of our finite element methods for PDEs in non–divergence form is to establish global and local discrete Calderon–Zygmund-type estimates similar to (2.21) for . These results are presented in the following two lemmas.

###### Lemma 2.6

There exists such that for all there holds

 (2.26) ∥wh∥W2,ph(Ω)≲∥L\rm 0% ,hwh∥Lp(Ω)∀wh∈Vh.
{proof}

First note that (2.26) is equivalent to

 (2.27) ∥∥L−1\rm 0,hφh∥∥W2,ph(Ω)≲∥φh∥Lp(Ω)∀φh∈Vh.

For any fixed , let and . Therefore, and , respectively, are the solutions of the following two problems:

 (2.28) a0(w,v)=(φh,v)∀v∈H10(Ω),a0(wh,vh)=(φh,vh)∀vh∈Vh,

and thus, is the elliptic projection of .

By (2.21) we have

 (2.29) ∥w∥W2,p(Ω)≲∥φh∥Lp(Ω).

Using well–known finite element estimate results [4, Theorem 8.5.3], finite element interpolation theory, and (2.29) we obtain that there exists such that for all

 (2.30) ∥w−wh∥W1,p(Ω)≲∥w−Ihw∥W1,p(Ω)≲h∥w∥W2,p(Ω)≲h∥φh∥Lp(Ω).

It follows from the triangle inequality, an inverse inequality (see Lemma 2.2), the stability of , (2.29) and (2.30) that

 ∥w−wh∥W2,ph(Ω) ≲∥w−Ihw∥W2,ph(Ω)+∥Ihw−wh∥W2,ph(Ω) ≲∥w∥W2,p(Ω)+h−1∥Ihw−wh∥W1,p(Ω) ≲∥φh∥Lp(Ω)+h−1∥Ihw−w∥W1,p(Ω)+h−1∥w−wh∥W1,p(Ω) ≲∥φh∥Lp(Ω).

Thus,

 ∥wh∥W2,ph(Ω)≲∥w−wh∥W2,ph(Ω)+∥w∥W2,p(Ω)≲∥φh∥Lp(Ω),

which yields (2.27), and hence, (2.26).

###### Lemma 2.7

For and , define

 (2.31) BR(x\rm 0):={x∈Ω: |x−x\rm 0|

Let with . Then there holds

 (2.32) ∥wh∥W2,ph(BR(x\rm 0))≲∥∥L\rm 0,hwh∥Lph(BR′(x\rm 0))∀wh∈Vh(BR(x\rm 0)),
{proof}

To ease notation, set and . Recalling (2.6), we have by Lemma 2.6,

 ∥wh∥W2,ph(BR) =∥wh∥W2,ph(Ω)≲∥L\rm 0% ,hwh∥Lp(Ω)≲∥L\rm 0,hwh∥Lph(Ω)=supvh∈Vha0(wh,vh)∥vh∥Lp′(Ω).

Set , so that . Denote by the indicator function of . Since on , we have

 a0(wh,vh)=a0(wh,χBR′′vh)=a0(wh,Ih(χBR′′vh))∀vh∈Vh.

Moreover, we have and

 ∥Ih(χBR′′vh)∥Lp′(BR′)=∥Ih(χBR′′vh)∥Lp′(Ω)≲∥χBR′′vh∥Lp′(Ω)≲∥vh∥Lp′(Ω).

Consequently,

 ∥wh∥W2,ph(BR) ≲supvh∈Vha0(wh,Ih(χBR′′vh))∥Ih(χBR′′vh)∥Lp′(BR′)≤supvh∈Vh(BR′)a0(wh,vh)∥vh∥Lp′(BR′) =supvh∈Vh(BR′)(L% \rm 0,hwh,vh)∥vh∥Lp′(BR′)=∥L\rm 0,hwh∥Lph(BR′).

## 3 C0 DG finite element methods and convergence analysis

### 3.1 The PDE problem

To make the presentation clear, we state the precise assumptions on the non-divergence form PDE problem (1.1). Let be a positive definite matrix-valued function with

 (3.1) λ|ξ|2A(x)ξ⋅ξ≤Λ|ξ|2∀ξ∈Rn,x∈¯¯¯¯Ω

and constants Under the above assumption, is known to be uniformly elliptic, hence, strong solutions (i.e., solutions) of problem (1.1) must satisfy the Aleksandrov maximum principle [14, 10, 16].

By the theory for the second order non-divergence form uniformly elliptic PDEs [14, Chapter 9], we know that if , for any with , there exists a unique strong solution to (1.1) satisfying

 (3.2) ∥u∥W2,p(Ω)≲∥f∥Lp(Ω).

Moreover, when and , it is also known that [15, 13, 3, 18, 2] the above conclusion holds if is a convex domain.

For the remainder of the paper, we shall always assume that is positive definite satisfying (3.1), and problem (1.1) has a unique strong solution which satisfies the Calderon-Zygmund estimate (3.2).

### 3.2 Formulation of C0 DG finite element methods

The formulation of our DG finite element method for non-divergence form PDEs is relatively simple, which is inspired by the finite element method for divergence form PDEs and relied only on an unorthodox integration by parts.

To motivate its derivation, we first look at how one would construct standard finite element methods for problem (1.1) when the coefficient matrix belongs to . In this case, since the divergence of (taken row-wise) is well defined, we can rewrite the PDE (1.1a) in divergence form as follows:

 (3.3) −∇⋅(A∇u)+(∇⋅A)⋅∇u=f.

Hence, the original non-divergence form PDE is converted into a “diffusion-convection equation” with the “diffusion coefficient” and the “convection coefficient” .

A standard finite element method for problem (3.3) is readily defined as seeking such that

 (3.4) ∫Ω(A∇uh)⋅∇vhdx+∫Ω(∇⋅A)⋅∇uhvhdx=∫Ωfvhdx∀vh∈Vh.

Now come back to the case where only belongs to . In our setting, the formulation (3.4) is not viable any more because does not exist as a function. To circumvent this issue, we apply the DG integration by parts formula (2.23) to the first term on the left-hand side of (3.4) with and in (3.4) is understood piecewise, we get

 (3.5) −∫Ω(A:D2huh)vhdx+∑e∈EIh∫e[[A∇uh]]vhds=∫Ωfvhdx∀vh∈Vh.

Here we have used the fact that and .

No derivative is taken on in (3.5), so each of the terms is well defined on . This indeed yields the DG formulation of this paper.

{definition}

The discontinuous Galerkin (DG) finite element method is defined by seeking such that

 (3.6) ah(uh,vh)=(f,vh)∀vh∈Vh,

where

 (3.7) ah(wh,vh) :=−∫Ω(A:D2hwh)vhdx+∑e∈EIh∫e[[A∇wh]]vhds, (3.8) (f,vh) :=∫Ωfvhdx∀vh∈Vh.

A few remarks are given below about the proposed DG finite element method.

###### Remark 3.1

(a) The above method is also defined for