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

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


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.


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

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

Since is not smooth, the PDE 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 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 almost everywhere provided , and . In addition, the existence of a strong solution to in two-dimensions and on convex domains is proved in [18].

Due to their non-divergence structure, designing convergent numerical methods, in particular, Galerkin-type methods, for problem has been proven to be difficult. Very few such results are known in the literature. Nevertheless, even problem 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 . 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 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 . 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.

= [rectangle, draw, fill=white!20, text width=13em, text centered, rounded corners, minimum height=4em] = [rectangle, draw, fill=white!20, text width=17em, text centered, rounded corners, minimum height=4em] = [draw, -latex’]

Outline of the convergence proof.
Outline of the convergence proof.

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]). 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 ?. 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 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 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 . See Figure ? 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 . 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.

2Notation and preliminary results

2.1Mesh 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:

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

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

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

We note that is non-trivial for .

Associated with , we define a semi-norm on for

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

Let be the projection defined by

It is well known that [8] satisfies for any

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

By , it is easy to see that is a norm on . Moreover by

2.2Some basic properties of 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.

By , and inverse estimates [7], we have

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

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

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

Writing and integrating by parts, we find

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

Likewise, by Lemma ? we have

Combining – we obtain the desired result. The proof is complete.

2.3Stability estimates for auxiliary PDEs with constant coefficients

In this subsection, we consider a special case of 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

The operator induces the following bilinear form:

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


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

Recall the following DG integration by parts formula:

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

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,

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 for . These results are presented in the following two lemmas.

First note that is equivalent to

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

and thus, is the elliptic projection of .

By we have

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

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


which yields , and hence, .

To ease notation, set and . Recalling , we have by Lemma ?,

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

Moreover, we have and


3 DG finite element methods and convergence analysis

3.1The PDE problem

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

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

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

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

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

3.2Formulation of 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 when the coefficient matrix belongs to . In this case, since the divergence of (taken row-wise) is well defined, we can rewrite the PDE in divergence form as follows:

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 is readily defined as seeking such that

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

Here we have used the fact that and .

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

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

3.3Stability analysis and well-posedness theorem

As in Section 2.3, using the bilinear form we can define the finite element approximation (or projection) of on , that is, we define by

Trivially, can be rewritten as: Find such that

Similar to the argument for , we can extend the domain of to the broken Sobolev space , that is, (abusing the notation) we define by

The main objective of this subsection is to establish a stability estimate for the operator on the finite element space . From this result, the existence, uniqueness and error estimate for will naturally follow. The stability proof relies on several technical estimates which we derive below. Essentially, the underlying strategy, known as a perturbation argument in the PDE literature, is to treat the operator locally as a perturbation of a stable operator with constant coefficients. The following lemma quantifies this statement.

Since is continuous on , it is uniformly continuous. Therefore for every there exists such that if satisfy , there holds . Consequently for any

with .

Set and consider , and . Since , it follows from , , , , and that

The desired inequality now follows from the definition of .

For to be determined below, let as in Lemma ?. Let and set . Then by Lemmas ? and ? with and , we have for any

For sufficiently small (depending only on ), we can kick back the first term on the right-hand side. This completes the proof.

Set . By the definition of , , and , we have for any