Analysis of mixed discontinuous Galerkin formulations for quasilinear elliptic problems

# Analysis of mixed discontinuous Galerkin formulations for quasilinear elliptic problems

Mohammad Zakerzadeh22footnotemark: 2    Georg May22footnotemark: 2 The research of the authors was supported by the Deutsche Forschungsgemeinschaft (German Research Association) through grant GSC 111.
###### Abstract

In this manuscript we present an approach to analyze the discontinuous Galerkin solution for general quasilinear elliptic problems. This approach is sufficiently general to extend most of the well-known discretization schemes, including BR1, BR2, SIPG and LDG, to nonlinear cases in a canonical way, and to establish the stability of their solution. Furthermore, in case of monotone and globally Lipschitz problems, we prove the existence and uniqueness of the approximated solution and the -optimality of the error estimate in the energy norm as well as in the norm.

d
11footnotetext: Aachen Institute for Advanced Study in Computational Engineering Science, RWTH Aachen, 52062 Aachen, Germany ({zakerzadeh, may}@aices.rwth-aachen.de).

iscontinuous Galerkin methods, nonlinear elliptic problems, optimal estimates

{AMS}

65N12, 65N15, 65N30

## 1 Introduction

A great amount of research has been devoted to the analysis of discontinuous Galerkin (DG) schemes for (non)linear elliptic and parabolic problems. This goes back to work of Babuska [4], Wheeler [41], Arnold [2] and Dupont et al. [29] for interior penalty (IP) methods. In later work, e.g. [9, 37], a non-symmetric IP scheme was presented and analyzed. Other classes of discretization techniques have also been introduced, e.g., the first and second method of Bassi–Rebay [7, 8], or Shu and Cockburn’s LDG method [16]. For a comprehensive literature study of these problems we refer to [15] and [3].

In their seminal paper [3], Arnold et al. provided a unified framework for obtaining different classes of DG methods for the linear Poisson problem. For nonlinear problems the picture is more complicated, and obtaining different classes of methods usually follows very different approaches; while IP methods often have been introduced and analyzed in the primal form [17, 23, 24, 34, 28], for LDG methods the usual approach is to start from a mixed formulation, e.g., [11, 22, 12] or [42]. Originating from this, different techniques are usually employed in the analysis of each family. The LDG methods [22, 11, 12, 42] have been formulated by using three unknowns in the mixed form as proposed in [14], which yields an inconsistent primal formulation as analyzed, e.g., in [11] (and [35] for the simpler linear case). In contrast, IP methods like [24, 23, 28, 34, 17] usually enjoy a consistent primal form.

On the other hand, the nonlinear versions of the Bassi-Rebay methods are often obtained by ad-hoc extension, mimicking the linear counter part, e.g. [6], and to the best knowledge of the authors, there does not exist rigorous analysis of the second Bassi–Rebay method for nonlinear problems, although they are widely used for nonlinear problems like the Navier–Stokes equations.

This manuscript aims to extend the discussion presented in [11, 24, 22, 23, 28, 34] in two ways: Firstly, we try to treat different types of discretization in a canonical way, starting from a mixed formulation as [3]. We will show that our formulation leads to a slightly modified version of previously proposed nonlinear formulations for symmetric IP (SIPG) [23], and Bassi–Rebay [6], while we recover the LDG formulation of [11]. Here we are interested in these methods since, among the methods investigated in [3], these are the only ones that are stable and consistent for primal and adjoint solutions, and we might hope for extending these properties to the corresponding nonlinear discretization. Also we prove that the approximate solution of the discontinuous Galerkin problem is well-posed; i.e, unique and stable, provided that the diffusion operator is strongly monotone and globally Lipschitz continuous. This is comparable to analysis of [11] for LDG and somewhat similar to [17] for IP; however our results cover different formulations. It is plausible that the same kind of analysis is applicable to cases which do not satisfy either strong monotonicity or global Lipschitz continuity as in [34, 24, 23]; nevertheless, we do not address this extension here and restrict our analysis to these two assumptions.

Secondly, we are going to provide an optimal error estimate for SIPG and Bassi–Rebay methods in terms of mesh size, both in energy norm and in the -norm. Such estimates have previously been derived for LDG methods, in [11] for monotone and globally Lipschtitz continuous problem, and in [22] for cases which neither of these two assumptions hold. For the SIPG method, error estimates have been presented in [23, 24] and for incomplete penalty (IIPG) in [34], for non-monotone and not globally Lipschitz operators. Also we have the results of [17] for IIPG in case of monotone and globally Lipschitz operators. Let us remark that since the formulations presented in [34, 17] are adjoint inconsistent, they derived the error estimate only in the energy norm and did not deal with the corresponding adjoint problem.

Besides the canonical approach of discretization we present here, the rigorous analysis of the second method of Bassi–Rebay is novel in the literature, as well as the different approach we adopted in the error estimate for an asymptotically adjoint consistent formulation. Furthermore, we present explicit conditions of stability for all formulations. To the best knowledge of the authors this has not been done before for a nonlinear version of Bassi–Rebay methods and is similar to the explicit bounds for the SIPG method in [39, 18] and [31]. It is worth mentioning that for the SIPG method, unlike the cited literature, our stability analysis remains valid for degenerate diffusion.

The structure of this paper is as follows: In §2 we provide a short review on the quasi-linear elliptic problems and the properties of the diffusion operator. In §3 we review some approximation results as well as the properties of triangulation of the computational domain. Section 4 deals with our canonical approach for obtaining different DG formulation in the primal form, and later in §5, we analyze the consistency and adjoint consistency of the presented formulations. The sections 6 and 7 are devoted to the stability and uniqueness of the DG approximate solution, respectively. Note that before reaching section 7 we do not require the monotonicity of the operator and the stability result is valid for more general problems. Section 8 presents the optimal error convergence estimate in both energy and norms.

## 2 Quasi-linear elliptic problems

We consider the following quasilinear elliptic problem

 (1) −∇⋅a(x,u,∇u) =f, in Ω, (2) u =uD, on ∂Ω,

where and . For simplicity we will set in later analysis. Also is a bounded and simply connected domain in . For brevity one might use the notation where , such that and . Moreover, by and , we denote the derivative of with respect to its second and third arguments.

In the general theory of nonlinear elliptic problems (see [32, 45]), it is usually assumed that the function satisfies some conditions:

1. The functions are continuous in and satisfy the following growth condition;

 (3) |ai(x,ζ)|≤c1(1+2∑i=0|ζi|)+|ϕi(x)|,∀ζ∈R3,∀x∈Ω,

where and , for .

2. There exist two constants such that for all

 (4) 0<λ2∑k=1ψ2k≤2∑j=1k=0∂aj(⋅,ζ)∂ζkψjψk≤Λ2∑k=1ψ2k.

We also might consider a relaxed version of (4) as the following

 (5) 0<λ2∑k=1ψ2k≤2∑j=1k=1∂aj(⋅,ζ)∂ζkψjψk≤Λ2∑k=1ψ2k.

For simplicity in the later analysis we assume that is symmetric. Then we can interpret (5) as imposition of lower and upper bound on the eigenvalues of ; and , respectively.

3. The functions are in , i.e., they are twice continuously differentiable functions with all the derivatives through second order being bounded.

Note that 1 is required to guarantee the meaningfulness of the of the definition of the discrete problem. On the other hand the assumptions 2 and 3 need to be valid only for the arguments provided in sections 7 and 8 for uniqueness of the discrete solution and error estimate, respectively. Furthermore, the relaxed version of assumption 2 in (5) is required in the stability analysis in section 6.

Let us note the following important definition for nonlinear elliptic PDEs, the so-called strong monotonicity and global Lipschitz continuity property of the diffusion operator (as [19], [44, p. 476]): {definition} We say a diffusion operator satisfies strong monotonicity property, if there exists a constant such that

 (6)

for all and . Moreover we call it a globally Lipschitz operator if there exists a constant such that

 (7) a(x,ξ′,η′)−a(x,ξ,η)≤C%lc[|ξ−ξ′|2+|η−η′|2]1/2,

for all and .

One can show that if the left and right inequalities in (4) are satisfied, the diffusion operator is strongly monotone and globally Lipschitz continuous in the sense of Definition 2. We skip the proof here and refer to [5, Lemma 2.1] or [17, Lemma 4].

As examples of problems settled in the category of (1) we have

1. Newtonian flow model

 (8) a(x,u,∇u)=a(x,u)∇u
2. Non-Newtonian flow model

 (9) a(x,u,∇u)=a(x,|∇u|)∇u

with an specific case of mean curvature flow

 (10) a(x,u,∇u)=∇u(1+|∇u|)1/2

We remark that for the case (b), following [11, 28], we typically assume the monotonicity of the diffusion operator (see Definition 2), while the diffusion operator in case (a) is not usually a monotone operator (only under very restrictive assumptions, see [1]). Since the presentation of the primal form of DG formulations in section 4 and stability result in section 6 are not affected by the monotonicity property, we do not exclude this case now. But in obtaining a priori error estimate in section 8, we will only consider monotone cases and one can refer to [34], [24] or [23] for a priori estimate for non-monotone cases. One may check that the mean curvature flow example satisfies the assumption 2; hence it is strongly monotone and globally Lipschitz continuous, see [23, section 5].

## 3 Preliminaries

First let us set a notation convention and suppress the dependence of the diffusion operator on the space variable ; henceforth we write instead of , while we still allow explicit dependence on .

Now, as a tool that we will use later in sections 6 and 8, let us consider the following integral form of Taylor’s formula; for and one has

 a(v,z)−a(u,w) =au(u,w)(v−u)+az(u,w)(z−w)+Ra(u−v,w−z) (11) =~au(u,w)(v−u)+~az(u,w)(z−w)

where is defined as

 (12) Rai=~auu(u−v)2+(w−z)t~azz(u,w)(w−z)+2~auz(u,w)(z−w)(u−v).

for . Moreover we define as

 ~au(u,w) =∫10au(v(t),z(t))dt, ~az(u,w)=∫10az(v(t),z(t))dt ~auu(u,w) =∫10(1−t)au(v(t),z(t))dt, ~auz(u,w)=∫10(1−t)au(v(t),z(t))dt, ~azz(u,w) =∫10(1−t)au(v(t),z(t))dt

where and .

Seeking clarity, sometimes we denote by all four arguments, e.g. instead of .

### 3.1 Triangulation and finite element space

Here we are going to consider the boundary of the domain sufficiently smooth (in order to apply the duality argument, see section 5.2), e.g. a convex polygon. Then we consider a shape-regular triangulation on as composed of (non-overlapping) triangular or rectangular elements (with possible hanging nodes) and is the diameter of each . Also we define and is the outward normal to . In the following we assume that is of bounded variation, that is, there exists a constant such that

 (13) l−1≤hκhκ′

where share an edge. This bounded variation property means that there is an upper bound for the number of neighboring elements of each , denoted by . In case that has no hanging nodes and for triangular and rectangular elements, respectively. We also need the following quasi-uniformity property of the mesh in the error analysis in section 8.2, that is

 (14) h≤Chκ,∀κ∈Th.

We denote the skeleton of the triangulation , i.e. the set of all edges of , by . Also we denote the set of boundary and interior edges of by and , respectively, and the length of edge by . Following [3], let us fix some definitions for the jumps and average of the discontinuous functions on the skeleton . Let us set the trace values as . For any interior edge , where is the common edge of , and for all , we define

 (15) {{w}}=12(wκ,e+wκ′,e),\llbracketw\rrbracket=wκ,eνκ+wκ′,eνκ′

and similarly for all

 (16) {{τ}}=12(τκ,e+τκ′,e),\llbracketτ\rrbracket=τκ,e⋅νκ+τκ′,e⋅νκ′.

For any boundary edge we define

 (17) \llbracketw\rrbracket=wκ,eνκ,{{τ}}=τκ,e,

for all and .

Moreover, let us consider the following broken Sobolev space on the triangulation ; for

 (18) Wsr(Ω,Th)={v∈Lr:v|κ∈Wsr,∀κ∈Th},

with the corresponding norm and seminorm

 (19) ∥v∥Wsr(Ω,Th)=(∑κ∈Th∥v∥rWsr(κ))1/r,|v|Wsr(Ω,Th)=(∑κ∈Th|v|rWsr(κ))1/r,

and for the case the associated norm and seminorm are defined as

 (20) ∥v∥Ws∞(Ω,Th)=maxκ∈Th∥v∥Ws∞(κ),|v|Ws∞(Ω,Th)=maxκ∈Th |v|Ws∞(κ),

when and are the standard Sobolev norms on . Also we denote as by tradition.

Now let define the following finite dimensional approximation spaces

 (21) Vh,q :={uh∈L2(Ω):uh|κ∈Pq(κ),∀κ∈Th}, (22) Σh,p :={θh∈[L2(Ω)]2:θh|κ∈[Pp(κ)]2,∀κ∈Th},

with and or (To satisfy the inclusion property as [3].) Here by we denote the space of the polynomials of total degree on and restricted to .

Moreover, let us define two lifting operator and as

 (23) ∫Ωr(φ)⋅τdx=−∑e∈Eh∫eφ⋅{{τ}}ds,∫Ωl(φ)⋅τdx=−∑e∈Eh,I∫eφ\llbracketτ\rrbracketds,

for all . Using the Riesz representation theorem one can prove the existence and uniqueness of the lifting operators introduced by (23) (see e.g., [11, Lemma 3.3]). Also we define an edge-wise version of right and left lifting operators as and , such that

 (24) ∫Ωre(φ)⋅τdx =−∫eφ⋅{{τ}}ds, ∀τ∈Σh,p,∀e∈Eh, (25) ∫Ωle(φ)⋅τdx =−∫eφ\llbracketτ\rrbracketds, ∀τ∈Σh,p,∀e∈Eh,I

for all edges . Also by noting that and, by applying Cauchy-Schwarz inequality and taking the norm over , one has

 (26) ∥r(φ)∥2L2(Ω)=∥∑e∈Ehre(φ)∥2L2(Ω)≤Nl∑e∈Eh∥re(φ)∥2L2(Ω).

Furthermore, one might note that for any and it holds and consequently

 (27) ∥l(φ)∥2L2(Ω)=∥∑e∈Ehle(φ)∥2L2(Ω)≤Nl∑e∈Eh∥le(φ)∥2L2(Ω)≤4Nl∑e∈Eh∥re(φνκ)∥2L2(Ω).

using Cauchy-Schwarz inequality.

Also let us introduce the space and the corresponding energy norm as

 (28) |||v|||2h:=|v|21,Ω+∑e∈Eh||re(\llbracketvh\rrbracket)||2L2(Ω),∀v∈V(h)

as well as the following seminorm as

 (29) |v|2∗,h:=∑e∈Eh||re(\llbracketv\rrbracket)||2L2(Ω),∀v∈V(h).

From the structure of (29), using (26) and (27) reads, there exists such that

 (30) ∥r(v)∥2L2(Ω)+∥l(v)∥2L2(Ω)≤Cs|v|2∗,h

for all . Moreover, we have the following estimate from [10, Lemma 2] {lemma} There exist two positive constants , such that for all we have

 (31) Crh−1/2e∥\llbracketw\rrbracket∥L2(e)≤∥re(\llbracketw\rrbracket)∥L2(Ω)≤CRh−1/2e∥\llbracketw\rrbracket∥L2(e),∀w∈V(h),

where the constants are independent and only depend on the minimum angle of the triangles and the polynomial degree . {proof} The original proof in [10] was presented for , while in [3] the extended version to was given, exploiting the Sobolev embedding to deduce in . Then (31) trivially holds for all . For the details on the explicit value of and we refer to [10] and [40].

In order to see the relation between and , let us remark the following relation from [2, Lemma 2.1], for all

 (32) ∥v∥2L2(Ω)≤C(Ω,Th)[∥∇v∥2L2(Ω)+∑e∈Ehh−1e∥\llbracketv\rrbracket∥2L2(e)],

which holds when is convex. For a similar result on non-convex domains we refer to [11, 21, 38].

Now, restricted to and by using (31), (32) reduces to the following Poincaré-type inequality

 (33) ∥vh∥2L2(Ω)≤Cen|||vh|||2h,

with some , and by the definition of the energy norm (28) one can write

 (34) ∥vh∥2H1(Ω,Th)≤(Cen+1)|||vh|||2h.

Also we will need the following inverse inequality {lemma} Let consider , then for there exists a constant such that

 (35) ∥vh∥Lr(κ)≤Cinvh2/r−1κ∥vh∥L2(κ).

The proof of this lemma can be found, e.g., in [13, p. 140] and we skip it.

### 3.2 Approximation properties

First, let state some approximation properties in the next two lemmas {lemma} [23, Lemma 2.1] For , there exists a positive constant , depending on and , but independent of and , and a sequence , such that

1. for any

 (36) ∥ϕ−ϕhκ∥Hl(e)≤CAhκμ−l−1/2∥ϕ∥Hs(κ)
2. for any

 (37) ∥ϕ−ϕhκ∥Wlr(κ)≤CAhκμ−l−1+2/r∥ϕ∥Hs(κ)

where .

For the proof of this lemma we refer to [23, Lemma 2.1] and the references cited therein.

Using Lemma 3.2 and the properties of the energy norm (28), we present the following approximation result {lemma} For there exists a constant independent of and , and a mapping such that

 (38) |||ϕ−πhϕ|||h≤C′A(∑κ∈Thh2(μ−1)κ|ϕ|2Hs(κ))1/2,

where . {proof} The proof exploits Lemmas 3.1 and 3.2 and the definition of which provides continuity of . Then the proof follows the same lines as [3, section 4.3.].

Note that in our analysis in this work, we do not need to specify the explicit type of projection . Hence we leave it undefined with merely assuming that it satisfies the approximation property described in Lemmas 3.2 and 3.2. We refer to [11] and [23] for explicit examples of such projections. Let us remark that the Galerkin projection that we introduce in (47) also satisfies these approximation properties.

## 4 Discontinuous Galerkin formulation

Here we follow the formulation presented in [11] (also see [14] and [35]). Let us start with writing the nonlinear elliptic problem (1) as a system of first order nonlinear PDEs in terms of new variables :

 (39) −∇⋅σ =f, in Ω, (40) σ =a(u,θ), in Ω, (41) θ =∇u, in Ω, (42) u =0, on ∂Ω.

Our goal is to approximate the exact solution by discrete functions in the finite element space . The weak formulation can be written as

 (43) ∫Ωa(uh,θh)⋅ζhdx=∫Ωσh⋅ζhdx, ∀ζh∈Σh,p, (44) ∫Ωθh⋅τhdx+∫Ωuh(∇h⋅τh)dx=∑κ∈Th∫∂κ^uτh⋅νds, ∀τh∈Σh,p, (45) ∫Ωσh⋅∇hvhdx−∑κ∈Th∫∂κvh(^σ⋅ν)ds=∫Ωfvhdx, ∀vh∈Vh,q.

Here and are the element-wise version of the gradient and divergence operator, respectively.

This flux formulation is complete but the definition of the numerical fluxes and which depends on , and needs to be designed carefully. We postpone the explicit definition of these numerical fluxes till the next section, where we present the equivalent primal formulation; i.e., a formulation which has as its only unknown.

The only requirement we impose here is that the should be independent of and , i.e., , which provides the availability of the primal formulation. All the formulations we are going to present have this local property. For examples of other forms of non-local formulation and their analysis we refer to [42], [22] and [20].

### 4.1 Primal formulation

In order to obtain the primal formulation we need to solve the unknowns and in terms of . In the first step, by using (43), one can solve for as the Galerkin projection,

 (46) σh=Gh(a(uh,θh)),

where has the following property; for all

 (47) ∑κ∈Th∫κξ⋅τdx=∑κ∈Th∫κGh(ξ)⋅τdx,∀τ∈Σh,p.

Moreover, let us remark the following identity; for any and the following holds

 (48) ∑κ∈Th∫∂κvξ⋅νds=∑e∈Eh,I∫e{{v}}\llbracketξ\rrbracketds+∑e∈Eh∫e\llbracketv\rrbracket⋅{{ξ}}ds.

Applying (48) in (44) and (45), one can write

 ∫Ωθh⋅τhdx−∫Ω∇huh⋅τhdx+∑e∈Eh,I∫e{{uh−^u}}\llbracketτh\rrbracketds+∑e∈Eh∫e\llbracketuh−^u\rrbracket⋅{{τh}}ds=0, ∫Ωσh⋅∇hvhdx−∑e∈Eh∫e{{^σ}}⋅\llbracketvh\rrbracketds−∑e∈Eh,I∫e\llbracket^σ\rrbracket{{vh}}ds=∫Ωfvhdx.

Now let us require the numerical fluxes and to be conservative, which leads to

 (49) \llbracket^u\rrbracket=0,{{^u}}=^u,\llbracket^σ\rrbracket=0,{{^σ}}=^σ,

on any . Also in accordance with the Dirichlet boundary condition (42), we set on . Using (49), one can simplify the weak formulation as

 (50) ∫Ωθh⋅τhdx−∫Ω∇huh⋅τhdx+∑e∈Eh,I∫e({{uh}}−^u)\llbracketτh\rrbracketds+∑e∈Eh∫e\llbracketuh\rrbracket⋅{{τh}}ds=0, (51) ∫Ωσh⋅∇hvhdx−∑e∈Eh∫e^σ⋅\llbracketvh\rrbracketds=∫Ωfvhdx.

Using the definition of lifting operator (23) and (46), we arrive at

 (52) θh =∇huh+r(\llbracketuh\rrbracket)+l({{uh}}−^u), (53) σh =Gh(a(uh,∇huh+r(\llbracketuh\rrbracket)+l({{uh}}−^u))).

Now can be solved locally in terms of by inserting (52) and (53) in (51). Then one can obtain the primal formulation as the following; find such that

 (54) B(uh,vh)=F(vh),∀vh∈Vh,q,

where , and

 (55) B(uh,vh)=∫Ωa(uh,∇huh+r(\llbracketuh\rrbracket)+l({{uh}}−^u))⋅∇hvhdx−∑e∈Eh∫e^σ⋅\llbracketvh\rrbracketds.

Note that in the rest of the paper we might abuse the notation defined in (52) in the operator form as .

The only undefined part in the primal formulation (55) is the explicit definition of the numerical fluxes. We are going to consider four different formulations by adopting different and :

1. BR1. The BR1 formulation is defined as in [7],

 ^u={{{uh}}on Eh,I0on Eh,∂,^σ={{σh}}on Eh.

Hence, . This leads to the following quasi-linear formulation

 (56) B(uh,vh)= ∫Ωa(uh,θh(uh))⋅θh(vh)dx.
2. BR2. The BR2 formulation, inherited from the original definition [8], is defined as

 ^u={{{uh}}on Eh,I0on Eh,∂,^σ={{Gh(a(uh,∇uh)+ηea(uh,re(\llbracketuh\rrbracket)))}}on Eh

with some as the stabilization parameter. Here, is the same as in BR1, and the primal formulation reads

 (57) B(uh,vh) =∫Ωa(uh,θh)⋅∇hvh+a(uh,∇huh)⋅r(\llbracketvh\rrbracket)dx +∑e∈Ehηe∫Ωa(uh,re(\llbracketuh\rrbracket))⋅re(\llbracketvh\rrbracket)dx.
3. SIPG. Similar to [3], we choose the fluxes as

 ^u={{{uh}}on Eh,I0on Eh,∂,^σ={{Gh(a(uh,∇uh)}}−μehe\llbracketuh\rrbracketon Eh

with some penalty parameter . Similarly to BR2, the primal formulation is given as

 (58) B(uh,vh) =∫Ωa(uh,θh)⋅∇hvh+a(uh,∇huh)⋅r(\llbracketvh\rrbracket)dx+∑e∈Ehμehe∫e\llbracketuh\rrbracket⋅\llbracketvh\rrbracketds.
4. LDG. The LDG formulation, inherited from the original version in [16], can be obtained by setting

where is some mesh-dependent parameter and constant on each edge. Also like SIPG, we have the penalty parameter . Note that in LDG formulation