Analysis of mixed discontinuous Galerkin formulations for quasilinear elliptic problems

TL;DR

This paper develops a unified framework for analyzing discontinuous Galerkin methods applied to nonlinear quasilinear elliptic problems, establishing stability, existence, uniqueness, and optimal error estimates.

Contribution

It extends well-known DG schemes to nonlinear problems in a canonical way and proves their stability and optimal error bounds.

Findings

Stability of DG solutions for nonlinear problems established

Existence and uniqueness proven for monotone, Lipschitz cases

Optimal error estimates in energy and L2 norms obtained

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 $h$-optimality of the error estimate in the energy norm as well as in the $L_2$ norm.

Equations

$$-\nabla \cdot \mathbf{a}(x,u,\nabla u)$$

$$u$$

$$|a_{i}(x,\bm{\zeta})|\leq c_{1}\Big{(}1+\sum_{i=0}^{2}|\zeta_{i}|\Big{)}+|\phi_{i}(x)|,\qquad\forall\bm{\zeta}\in\mathbb{R}^{3},\forall x\in\Omega,$$

$$0<\lambda \sum _{k=1}^2{\psi }_k^2\le \sum _{\begin{array}{c}j=1\\ k=0\end{array}}^2\frac{\partial a_j(\cdot ,\mathbf{\zeta })}{\partial {\zeta }_k}{\psi }_j{\psi }_k\le \Lambda \sum _{k=1}^2{\psi }_k^2.$$

$$0<\lambda \sum _{k=1}^2{\psi }_k^2\le \sum _{\begin{array}{c}j=1\\ k=1\end{array}}^2\frac{\partial a_j(\cdot ,\mathbf{\zeta })}{\partial {\zeta }_k}{\psi }_j{\psi }_k\le \Lambda \sum _{k=1}^2{\psi }_k^2.$$

$${\big{(}\bm{a}(x,\xi^{\prime},\bm{\eta}^{\prime})-\bm{a}(x,\xi,\bm{\eta})\big{)}}\cdot{(\bm{\eta}^{\prime}-\bm{\eta})}\geq C_{\text{sm}}|\bm{\eta}-\bm{\eta}^{\prime}|^{2}$$

$$\bm{a}(x,\xi^{\prime},\bm{\eta}^{\prime})-\bm{a}(x,\xi,\bm{\eta})\leq C_{\text{lc}}\big{[}|\xi-\xi^{\prime}|^{2}+|\bm{\eta}-\bm{\eta}^{\prime}|^{2}\big{]}^{1/2},$$

$$\mathbf{a}(x,u,\nabla u)=\mathbf{a}(x,u)\nabla u$$

$$\mathbf{a}(x,u,\nabla u)=\mathbf{a}(x,|\nabla u|)\nabla u$$

$$\mathbf{a}(x,u,\nabla u)=\frac{\nabla u}{(1+|\nabla u|{)}^{1/2}}$$

$$\mathbf{a}(v,\mathbf{z})-\mathbf{a}(u,\mathbf{w})$$

$$={\tilde{\mathbf{a}}}_u(u,\mathbf{w})(v-u)+{\tilde{\mathbf{a}}}_{\mathbf{z}}(u,\mathbf{w})(\mathbf{z}-\mathbf{w})$$

$$R_{{\mathbf{a}}_i}={\tilde{\mathbf{a}}}_{uu}(u-v{)}^2+(\mathbf{w}-\mathbf{z}{)}^t{\tilde{\mathbf{a}}}_{\mathbf{z}\mathbf{z}}(u,\mathbf{w})(\mathbf{w}-\mathbf{z})+2{\tilde{\mathbf{a}}}_{u\mathbf{z}}(u,\mathbf{w})(\mathbf{z}-\mathbf{w})(u-v).$$

$${\tilde{\mathbf{a}}}_u(u,\mathbf{w})$$

$${\tilde{\mathbf{a}}}_{uu}(u,\mathbf{w})$$

$${\tilde{\mathbf{a}}}_{\mathbf{z}\mathbf{z}}(u,\mathbf{w})$$

$$l^{-1}\le \frac{h_{\kappa }}{h_{{\kappa }^{\prime }}}<l,$$

$$h\le C{h}_{\kappa },\forall \kappa \in {\mathcal{T}}_h.$$

$$\{\{w\}\}=\frac{1}{2}(w_{\kappa ,e}+w_{{\kappa }^{\prime },e}),[[w]]=w_{\kappa ,e}{\nu }_{\kappa }+w_{{\kappa }^{\prime },e}{\nu }_{{\kappa }^{\prime }}$$

$$\{\{\mathbf{\tau }\}\}=\frac{1}{2}({\mathbf{\tau }}_{\kappa ,e}+{\mathbf{\tau }}_{{\kappa }^{\prime },e}),[[\mathbf{\tau }]]={\mathbf{\tau }}_{\kappa ,e}\cdot {\nu }_{\kappa }+{\mathbf{\tau }}_{{\kappa }^{\prime },e}\cdot {\nu }_{{\kappa }^{\prime }}.$$

$$[[w]]=w_{\kappa ,e}{\nu }_{\kappa },\{\{\mathbf{\tau }\}\}={\mathbf{\tau }}_{\kappa ,e},$$

$$W_r^s(\Omega ,{\mathcal{T}}_h)=\{v\in L_r:v{|}_{\kappa }\in W_r^s,\forall \kappa \in {\mathcal{T}}_h\},$$

$$\|v\|_{W_{r}^{s}(\Omega,\mathcal{T}_{h})}=\Big{(}\sum_{\kappa\in\mathcal{T}_{h}}\|v\|^{r}_{W_{r}^{s}(\kappa)}\Big{)}^{1/r},\qquad|v|_{W_{r}^{s}(\Omega,\mathcal{T}_{h})}=\Big{(}\sum_{\kappa\in\mathcal{T}_{h}}|v|^{r}_{W_{r}^{s}(\kappa)}\Big{)}^{1/r},$$

$$\Vert v{\Vert }_{W_{\infty }^s(\Omega ,{\mathcal{T}}_h)}=\underset{\kappa \in {\mathcal{T}}_h}{\max }\Vert v{\Vert }_{W_{\infty }^s(\kappa )},|v{|}_{W_{\infty }^s(\Omega ,{\mathcal{T}}_h)}=\underset{\kappa \in {\mathcal{T}}_h}{\max }|v{|}_{W_{\infty }^s(\kappa )},$$

$$V_{h,q}$$

$${\Sigma }_{h,p}$$

$${\int }_{\Omega }r(\phi )\cdot \mathbf{\tau }\mathrm{d}x=-\sum _{e\in {\mathcal{E}}_h}{\int }_e\phi \cdot \{\{\mathbf{\tau }\}\}\mathrm{d}s,{\int }_{\Omega }l(\phi )\cdot \mathbf{\tau }\mathrm{d}x=-\sum _{e\in {\mathcal{E}}_{h,I}}{\int }_e\phi [[\mathbf{\tau }]]\mathrm{d}s,$$

$${\int }_{\Omega }{r}^e(\phi )\cdot \mathbf{\tau }\mathrm{d}x$$

$${\int }_{\Omega }{l}^e(\phi )\cdot \mathbf{\tau }\mathrm{d}x$$

$$\Vert r(\phi ){\Vert }_{L_2(\Omega )}^2=\Vert \sum _{e\in {\mathcal{E}}_h}{r}^e(\phi ){\Vert }_{L_2(\Omega )}^2\le N_l\sum _{e\in {\mathcal{E}}_h}\Vert r^e(\phi ){\Vert }_{L_2(\Omega )}^2.$$

$$\Vert l(\phi ){\Vert }_{L_2(\Omega )}^2=\Vert \sum _{e\in {\mathcal{E}}_h}{l}^e(\phi ){\Vert }_{L_2(\Omega )}^2\le N_l\sum _{e\in {\mathcal{E}}_h}\Vert l^e(\phi ){\Vert }_{L_2(\Omega )}^2\le 4N_l\sum _{e\in {\mathcal{E}}_h}\Vert r^e(\phi {\nu }_{\kappa }){\Vert }_{L_2(\Omega )}^2.$$

$$|||v||{|}_h^2:=|v{|}_{1,\Omega }^2+\sum _{e\in {\mathcal{E}}_h}||r^e([[v_h]])|{|}_{L_2(\Omega )}^2,\forall v\in V(h)$$

$$|v{|}_{\ast ,h}^2:=\sum _{e\in {\mathcal{E}}_h}||r^e([[v]])|{|}_{L_2(\Omega )}^2,\forall v\in V(h).$$

Full text

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

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 $L^{2}$-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 $L_{2}$ 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 ${L_{2}}$ norms.

2 Quasi-linear elliptic problems

We consider the following quasilinear elliptic problem

[TABLE]

where $f\in{L_{2}}(\Omega)$ and $u_{D}\in H^{1/2}(\partial\Omega)$. For simplicity we will set $u_{D}\equiv 0$ in later analysis. Also $\Omega$ is a bounded and simply connected domain in $\mathbb{R}^{2}$. For brevity one might use the notation $\bm{a}(\cdot,\bm{\zeta})=\bm{a}(\cdot,u,\nabla u)$ where $\bm{\zeta}\in\mathbb{R}^{3}$, such that $\zeta_{0}\equiv u$ and $\zeta_{i}\equiv u_{x_{i}},i=1,2$. Moreover, by $\bm{a}_{u}$ and \bm{a}_{\bm{z}}:=\big{[}\frac{\partial\bm{a}_{i}(\cdot,\bm{\zeta})}{\partial\zeta_{j}}\big{]}_{i,j=1,2}, we denote the derivative of $\bm{a}(x,u,\nabla u)$ 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 $\bm{a}(\cdot,\bm{\zeta})=(a_{1}(\cdot,\bm{\zeta}),a_{2}(\cdot,\bm{\zeta}))$ satisfies some conditions:

1. A(i)

The functions $a_{i}(x,\bm{\zeta}),i=1,2$ are continuous in $\Omega\times\mathbb{R}^{3}$ and satisfy the following growth condition;

[TABLE]

where $c_{1}>0$ and $\phi_{i}\in{L_{2}}(\Omega)$, for $i=1,2$. 2. A(ii)

There exist two constants $0<\lambda\leq\Lambda<\infty$ such that for all $\bm{\zeta},\bm{\psi}\in\mathbb{R}^{3}$

[TABLE]

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

[TABLE]

For simplicity in the later analysis we assume that $\bm{a}_{\bm{z}}$ is symmetric. Then we can interpret (5) as imposition of lower and upper bound on the eigenvalues of $\bm{a}_{\bm{z}}$; $\lambda$ and $\Lambda$, respectively. 3. A(iii)

The functions $a_{i}(x,\bm{\zeta}),i=1,2$ are in $\mathcal{C}_{b}^{2}(\Omega\times\mathbb{R}^{3})$, i.e., they are twice continuously differentiable functions with all the derivatives through second order being bounded.

Note that 3 is required to guarantee the meaningfulness of the of the definition of the discrete problem. On the other hand the assumptions A(ii) and A(iii) 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 A(ii) 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 1**.**

We say a diffusion operator $\bm{a}(x,u,\nabla u)$ satisfies strong monotonicity property, if there exists a constant $C_{\text{sm}}>0$ such that

[TABLE]

for all $\xi,\xi^{\prime}\in\mathbb{R}$ and $\bm{\eta},\bm{\eta}^{\prime}\in\mathbb{R}^{2}$. Moreover we call it a globally Lipschitz operator if there exists a constant $C_{\text{lc}}<\infty$ such that

[TABLE]

for all $\xi,\xi^{\prime}\in\mathbb{R}$ and $\bm{\eta},\bm{\eta}^{\prime}\in\mathbb{R}^{2}$.

One can show that if the left and right inequalities in (4) are satisfied, the diffusion operator $\bm{a}$ is strongly monotone and globally Lipschitz continuous in the sense of Definition 1. 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. (a)

Newtonian flow model

[TABLE] 2. (b)

Non-Newtonian flow model

[TABLE]

with an specific case of mean curvature flow

[TABLE]

We remark that for the case (b), following [11, 28], we typically assume the monotonicity of the diffusion operator (see Definition 1), 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 A(ii); 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 $x$; henceforth we write $\bm{a}(v,\bm{z})$ instead of $\bm{a}(x,v,\bm{z})$, while we still allow explicit dependence on $x$.

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 $u,v\in\mathbb{R}$ and $\bm{z},\bm{w}\in\mathbb{R}^{2}$ one has

[TABLE]

where $R_{\bm{a}}(u-v,\bm{w}-\bm{z})=(R_{\bm{a}_{1}}(u-v,\bm{w}-\bm{z}),R_{\bm{a}_{2}}(u-v,\bm{w}-\bm{z}))$ is defined as

[TABLE]

for $i=1,2$. Moreover we define $\tilde{\bm{a}}_{u},\tilde{\bm{a}}_{\bm{z}},\tilde{\bm{a}}_{uu},\tilde{\bm{a}}_{u\bm{z}},\tilde{\bm{a}}_{\bm{z}\bm{z}}$ as

[TABLE]

where $v(t)=u+t(v-u)$ and $\bm{z}(t)=\bm{w}+t(\bm{z}-\bm{w})$.

Seeking clarity, sometimes we denote $\tilde{\bm{a}}_{u},\tilde{\bm{a}}_{\bm{z}},\tilde{\bm{a}}_{uu},\tilde{\bm{a}}_{u\bm{z}},\tilde{\bm{a}}_{\bm{z}\bm{z}}$ by all four arguments, e.g. $\tilde{\bm{a}}_{u}(u,\bm{w},v,\bm{z})$ instead of $\tilde{\bm{a}}_{u}(u,\bm{w})$.

3.1 Triangulation and finite element space

Here we are going to consider the boundary of the domain $\Omega$ 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 $\Omega$ as $\mathcal{T}_{h}=\{\kappa\}$ composed of (non-overlapping) triangular or rectangular elements (with possible hanging nodes) and $h_{\kappa}$ is the diameter of each $\kappa\in\mathcal{T}_{h}$. Also we define $h:=\max_{\kappa\in\mathcal{T}_{h}}h_{\kappa}$ and $\nu_{\kappa}$ is the outward normal to $\partial\kappa$. In the following we assume that $\mathcal{T}_{h}$ is of bounded variation, that is, there exists a constant $l>1$ such that

[TABLE]

where $\kappa,\kappa^{\prime}\in\mathcal{T}_{h}$ share an edge. This bounded variation property means that there is an upper bound for the number of neighboring elements of each $\kappa\in\mathcal{T}_{h}$, denoted by $N_{l}$. In case that $\mathcal{T}_{h}$ has no hanging nodes $N_{l}=3$ and $N_{l}=4$ for triangular and rectangular elements, respectively. We also need the following quasi-uniformity property of the mesh in the ${L_{2}}$ error analysis in section 8.2, that is

[TABLE]

We denote the skeleton of the triangulation $\mathcal{T}_{h}$, i.e. the set of all edges of $\kappa\in\mathcal{T}_{h}$, by $\mathcal{E}_{h}$. Also we denote the set of boundary and interior edges of $\mathcal{T}_{h}$ by $\mathcal{E}_{h,\partial}$ and $\mathcal{E}_{h,I}$, respectively, and the length of edge $e$ by $h_{e}$. Following [3], let us fix some definitions for the jumps and average of the discontinuous functions on the skeleton $\mathcal{E}_{h}$. Let us set the trace values as $w_{\kappa,e}=w_{\kappa}|_{e}$. For any interior edge $e\in\mathcal{E}_{h,I}$, where $e$ is the common edge of $\kappa,\kappa^{\prime}\in\mathcal{T}_{h}$, and for all $w\in\prod_{\kappa\in\mathcal{T}_{h}}{L_{2}}(\partial\kappa)$, we define

[TABLE]

and similarly for all $\bm{\tau}\in\prod_{\kappa\in\mathcal{T}_{h}}[{L_{2}}(\partial\kappa)]^{2}$

[TABLE]

For any boundary edge $e\in\mathcal{E}_{h,\partial}$ we define

[TABLE]

for all $w\in\prod_{\kappa\in\mathcal{T}_{h}}{L_{2}}(\partial\kappa)$ and $\bm{\tau}\in\prod_{\kappa\in\mathcal{T}_{h}}[{L_{2}}(\partial\kappa)]^{2}$.

Moreover, let us consider the following broken Sobolev space on the triangulation $\mathcal{T}_{h}$; for $1\leq r<\infty$

[TABLE]

with the corresponding norm and seminorm

[TABLE]

and for the case $r=\infty$ the associated norm and seminorm are defined as

[TABLE]

when $\|\cdot\|_{W_{r}^{s}(\kappa)}$ and $|\cdot|_{W_{r}^{s}(\kappa)}$ are the standard Sobolev norms on $\kappa$. Also we denote $W_{2}^{s}$ as $H^{s}$ by tradition.

Now let define the following finite dimensional approximation spaces

[TABLE]

with $q\geq 1$ and $p=q$ or $p=q-1$ (To satisfy the inclusion property $\nabla V_{h,q}\subset\Sigma_{h,p}$ as [3].) Here by $\mathcal{P}^{q}(\kappa)$ we denote the space of the polynomials of total degree $q$ on $\mathbb{R}^{2}$ and restricted to $\kappa$.

Moreover, let us define two lifting operator $r:[{L_{2}}(\mathcal{E}_{h})]^{2}\to\Sigma_{h,p}$ and $l:{L_{2}}(\mathcal{E}_{h,I})\to\Sigma_{h,p}$ as

[TABLE]

for all $\bm{\tau}\in\Sigma_{h,p}$. 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 $r^{e}:[{L_{2}}(e)]^{d}\to\Sigma_{h,p}$ and $l^{e}:{L_{2}}(e)\to\Sigma_{h,p}$, such that

[TABLE]

for all edges $e\in\mathcal{E}_{h}$. Also by noting that $r(\varphi)=\sum_{e\in\mathcal{E}_{h}}r^{e}(\varphi)$ and, by applying Cauchy-Schwarz inequality and taking the norm over $\mathcal{T}_{h}$, one has

[TABLE]

Furthermore, one might note that for any $e\in\mathcal{E}_{h,I}$ and $e\subset\partial\kappa,\kappa\in\mathcal{T}_{h}$ it holds $l^{e}(\varphi)=2r^{e}(\varphi\nu_{\kappa})$ and consequently

[TABLE]

using Cauchy-Schwarz inequality.

Also let us introduce the space $V(h):=V_{h,q}+H^{2}(\Omega)\bigcap H^{1}_{0}(\Omega)$ and the corresponding energy norm $|\!|\!|\cdot|\!|\!|_{h}:V(h)\to\mathbb{R}$ as

[TABLE]

as well as the following seminorm $|{\cdot}|_{{*},h}:V(h)\to\mathbb{R}$ as

[TABLE]

From the structure of (29), using (26) and (27) reads, there exists $C_{s}<\infty$ such that

[TABLE]

for all $v\in V(h)$. Moreover, we have the following estimate from [10, Lemma 2]

Lemma 2**.**

There exist two positive constants $C_{r},C_{R}>0$, such that for all $e\in\mathcal{E}_{h}$ we have

[TABLE]

where the constants are $h$ independent and only depend on the minimum angle of the triangles and the polynomial degree $q$.

Proof.

The original proof in [10] was presented for $w_{h}\in V_{h,q}$, while in [3] the extended version to $V(h)$ was given, exploiting the Sobolev embedding to deduce $V(h)\setminus V_{h,q}\subset H^{2}\subset\mathcal{C}(\Omega)$ in $\mathbb{R}^{2}$. Then (31) trivially holds for all $w\in V(h)\setminus V_{h,q}$. For the details on the explicit value of $C_{r}$ and $C_{r}$ we refer to [10] and [40]. ∎

In order to see the relation between $|\!|\!|\cdot|\!|\!|_{h}$ and $\|\cdot\|_{{L_{2}}(\Omega)}$, let us remark the following relation from [2, Lemma 2.1], for all $v\in H^{1}(\Omega,\mathcal{T}_{h})$

[TABLE]

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

Now, restricted to $v\in V(h)$ and by using (31), (32) reduces to the following Poincaré-type inequality

[TABLE]

with some $C_{\text{en}}<\infty$, and by the definition of the energy norm (28) one can write

[TABLE]

Also we will need the following inverse inequality

Lemma 3**.**

Let consider $v_{h}\in V_{h,q}$, then for $r\geq 2$ there exists a constant $C_{\text{inv}}>0$ such that

[TABLE]

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

[23, Lemma 2.1]** For $\phi\in H^{s}(\kappa)$, there exists a positive constant $C_{A}$, depending on $s$ and $q$, but independent of $\phi$ and $h_{\kappa}$, and a sequence $\phi_{h_{\kappa}}\in\mathcal{P}^{q}(\kappa)$, $q=1,2,\cdots$ such that

1. (i)

for any $s\geq l+\frac{1}{2}$

[TABLE] 2. (ii)

for any $0\leq l\leq s-1+\frac{2}{r}$

[TABLE]

where $\mu=\min(q+1,s)$.

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

Using Lemma 4 and the properties of the energy norm (28), we present the following approximation result

Lemma 5**.**

For $\forall\phi\in V(h)$ there exists a constant $C_{A}^{\prime}>0$ independent of $h$ and $\phi$, and a mapping $\pi_{h}:V(h)\to V_{h,q}$ such that

[TABLE]

where $\mu=\min(q+1,s)$.

Proof.

The proof exploits Lemmas 2 and 4 and the definition of $V(h)$ which provides continuity of $w\in V(h)\setminus V_{h,q}$. 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 $\pi_{h}$. Hence we leave it undefined with merely assuming that it satisfies the approximation property described in Lemmas 4 and 5. We refer to [11] and [23] for explicit examples of such projections. Let us remark that the Galerkin $[{L_{2}}]^{2}$ 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 $(\bm{\sigma},\bm{\theta},u)$:

[TABLE]

Our goal is to approximate the exact solution $(\bm{\sigma},\bm{\theta},u)$ by discrete functions $(\bm{\sigma}_{h},\bm{\theta}_{h},u_{h})$ in the finite element space $\Sigma_{h,p}\times\Sigma_{h,p}\times V_{h,q}$. The weak formulation can be written as

[TABLE]

Here $\nabla_{h}$ and $\nabla_{h}\cdot$ are the element-wise version of the gradient and divergence operator, respectively.

This flux formulation is complete but the definition of the numerical fluxes $\hat{u}$ and $\hat{\bm{\sigma}}$ which depends on $(\bm{\sigma}_{h},\bm{\theta}_{h},u_{h})$, 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 $u_{h}$ as its only unknown.

The only requirement we impose here is that the $\hat{u}$ should be independent of $\bm{\theta}_{h}$ and $\bm{\sigma}_{h}$, i.e., $\hat{u}=\hat{u}(u_{h})$, 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 $\bm{\sigma}_{h}$ and $\bm{\theta}_{h}$ in terms of $u_{h}$. In the first step, by using (43), one can solve for $\bm{\sigma}_{h}$ as the Galerkin $[{L_{2}}(\Omega)]^{2}$ projection,

[TABLE]

where $\mathcal{G}_{h}:[{L_{2}}(\Omega)]^{2}\to\Sigma_{h,p}$ has the following property; for all $\bm{\xi}\in[{L_{2}}(\Omega)]^{2}$

[TABLE]

Moreover, let us remark the following identity; for any $v\in\prod_{\kappa\in\mathcal{T}_{h}}{L_{2}}(\partial\kappa)$ and $\bm{\xi}\in\prod_{\kappa\in\mathcal{T}_{h}}[{L_{2}}(\partial\kappa)]^{2}$ the following holds

[TABLE]

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

[TABLE]

Now let us require the numerical fluxes $\hat{u}$ and $\hat{\bm{\sigma}}$ to be conservative, which leads to

[TABLE]

on any $e\in\mathcal{E}_{h,I}$. Also in accordance with the Dirichlet boundary condition (42), we set $\hat{u}=0$ on $e\in\mathcal{E}_{h,\partial}$. Using (49), one can simplify the weak formulation as

[TABLE]

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

[TABLE]

Now $(\bm{\theta}_{h},\bm{\sigma}_{h})$ can be solved locally in terms of $u_{h}$ by inserting (52) and (53) in (51). Then one can obtain the primal formulation as the following; find $u_{h}\in V_{h,q}$ such that

[TABLE]

where $\displaystyle\mathcal{F}(v_{h})=\int_{\Omega}fv_{h}\,\mathrm{d}x$, and

[TABLE]

Note that in the rest of the paper we might abuse the notation $\bm{\theta}_{h}$ defined in (52) in the operator form as $\bm{\theta}_{h}(\cdot):=\nabla_{h}(\cdot)+r(\llbracket\cdot\rrbracket)+l(\{\!\!\{\cdot\}\!\!\}-\hat{u}(\cdot))$.

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 $\hat{u}$ and $\hat{\bm{\sigma}}$:

1. (i)

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

[TABLE]

Hence, $\bm{\theta}_{h}=\nabla_{h}u_{h}+r(\llbracket u_{h}\rrbracket)$. This leads to the following quasi-linear formulation

[TABLE] 2. (ii)

BR2. The BR2 formulation, inherited from the original definition [8], is defined as

[TABLE]

with some $\eta_{e}>0$ as the stabilization parameter. Here, $\bm{\theta}_{h}$ is the same as in BR1, and the primal formulation reads

[TABLE] 3. (iii)

SIPG. Similar to [3], we choose the fluxes as

[TABLE]

with some penalty parameter $\mu_{e}>0$. Similarly to BR2, the primal formulation is given as

[TABLE] 4. (iv)

LDG. The LDG formulation, inherited from the original version in [16], can be obtained by setting

[TABLE]

where $\bm{\beta}\in[{L_{2}}(\mathcal{E}_{h,I})]^{2}$ is some mesh-dependent parameter and constant on each edge. Also like SIPG, we have the penalty parameter $\mu_{e}>0$. Note that in LDG formulation $\bm{\theta}_{h}=\nabla_{h}u_{h}+r(\llbracket u_{h}\rrbracket)+l({\bm{\beta}}\cdot{\llbracket u_{h}\rrbracket})$. Finally, the primal formulation reads

[TABLE]

Remark 1**.**

Seeking comparison of the current formulations with the already proposed versions in the literature, let us record a few observations. Firstly, it is worth pointing out that those fluxes which depend only on $u_{h}$ and $\bm{\sigma}_{h}$ (BR1, LDG) do not require any formal modification compared to the linear case. Furthermore, our LDG formulation is similar to [11] and the discretization proposed in [22] for a simpler form of nonlinearity in the diffusion.

For the SIPG formulation, looking at (58) and exploiting Taylor’s expansion (3), we can write

[TABLE]

Ignoring the boundary terms, the primal formulation proposed in [23] (or [27] in simpler case) reads as

[TABLE]

which shows two differences in the second term: in [23] there is $\bm{a}_{\bm{z}}$ instead of the average $\tilde{\bm{a}}_{\bm{z}}$, which is a consequence of the direct consideration of the primal form. Also in our formulation we have an additional Galerkin projection, which seems essential to later proofs in the paper. Also one might notice that while $\bm{a}_{\bm{z}}$ needs to be computed explicitly in [23], in our formulation $\tilde{\bm{a}}_{\bm{z}}$ only appears in the analysis and there is no need to compute and implement this term.

Furthermore, comparing the current primal BR2 formulation with the version proposed in [6] for the case of nonlinear Helmholtz equation, and by following similar arguments to those shown for SIPG, one observes that the only difference of two formulations is in the application of the additional Galerkin projection in our formulation.

5 Consistency and adjoint consistency

It is clear that due to the discrete nature of the lifting operator and the Galerkin projection, the primal formulation (54) is inconsistent with the exact solution (for all four presented methods). This is also the case for the formulations presented in [11] and [35]. Moreover, the scheme is inconsistent with the smooth solution of its corresponding adjoint problem constructed by the linearization in the neighborhood of the exact solution.

In order to investigate these consistency errors, we consider a smooth solution to (1). In sections 5.1 and 5.2 we prove that, despite of the primal and adjoint inconsistency of the formulations, and in case of sufficiently regular solutions, one obtains asymptotic consistency in the mesh refining limit, for primal as well as for the adjoint problem.

In the following of this paper, due to the required regularity in the adjoint analysis, we assume that the exact solution of the problem (1), denoted by $u$, belongs to $W_{\infty}^{2}(\Omega)$.

5.1 Consistency

Due to the assumption 3, the fact that $f\in{L_{2}}(\Omega)$ and $-\nabla\cdot\bm{a}(\cdot,u,\nabla u)=f$ on $\Omega$ for the exact solution $u$, we know that $\bm{a}(u,\nabla u)\in H(\mathrm{div};\Omega)$. Hence, $\llbracket\bm{a}(u,\nabla u)\rrbracket=0$ on each $e\in\mathcal{E}_{h,I}$. Also note that we have $\bm{\theta}_{h}(u)=\nabla u$ for all schemes. Let us first consider the schemes BR1, BR2, and SIPG. For these schemes, and for any $v\in V(h)$, one can rewrite the left hand side of (55) as

[TABLE]

Applying the divergence theorem on the first two terms on the right hand side, and using the continuity of $\bm{a}(u,\nabla u)\cdot\nu_{\kappa}$ on the interfaces yield

[TABLE]

Noting (1), the consistency error for the primal formulation is the same for BR1, BR2 and SIPG; that is

[TABLE]

Similarly, for LDG scheme, following the same lines as [11] we obtain

[TABLE]

In general, even for very regular solutions, or in case of linear diffusion $\bm{a}(u,\nabla u)=\nabla u$, this consistency error is not zero and in fact is equal to the Galerkin projection error into the polynomial space $\Sigma_{h,p}$. However, provided that the diffusion $\bm{a}$ is regular enough, the formulations are asymptotically consistent; i.e., $\|\mathcal{E}_{p}(u,\cdot)\|_{V^{\prime}(h)}\to 0$ as $h\to 0$. Note that $\|\cdot\|_{V^{\prime}(h)}$ is the dual norm on space $V(h)$ defined as

[TABLE]

where $A:V(h)\to V^{\prime}(h)$ is a linear continuous operator on $V(h)$.

In order to prove this asymptotic consistency, one might find an upper bound for $\mathcal{E}_{p}(u,v)$ which vanishes as $h$ goes to zero. Here we present a generalized form of [11, Lemma 5.2] for different discretizations, which provides us with such an estimate:

Lemma 6**.**

Assume $\bm{a}(u,\nabla u)\in H^{s_{*}}(\Omega,\mathcal{T}_{h})$ with some non-negative integer $s_{*}$. Then, there exists $C_{\text{con}}>0$, independent of $h$ but dependent on $q$ and $s_{*}$, such that for the schemes presented in section 4 the following holds

[TABLE]

for all $w\in V(h)$, where $\mu_{*}=\min(s_{*},q+1)$.

Proof.

The proof is provided in [11] for the LDG scheme, which in case $\bm{\beta}=0$ gives the desired result for BR1, BR2, and SIPG. ∎

5.2 Adjoint consistency

Adjoint consistency is an important property in obtaining the optimal ${L_{2}}$ convergence rate for the solution, and super-convergence of the target functionals (cf. [26] or [36]). In order to apply Aubin-Nitsche duality argument, let us consider the following auxiliary dual problem

[TABLE]

where $u$ is the exact solution of (1). From assumption A(iii) and $u\in W_{\infty}^{2}(\Omega)$, and provided that $u_{h}\in{L_{2}}(\Omega)$ (cf. section 6) one can check that

[TABLE]

Using (66) combined with the convexity of $\Omega$, it is a classical result that the unique solution of the adjoint problem, $\psi\in H^{2}\cap H^{1}_{0}$ satisfies the following elliptic regularity [25, Theorem 9.1.22]

[TABLE]

Moreover, from (66) and the structure of (64), one can conclude that $\bm{a}_{\bm{z}}(u,\nabla u)\nabla\psi\in H(\mathrm{div},\Omega)$. Therefore one gets $\llbracket\bm{a}_{\bm{z}}(u,\nabla u)\nabla\psi\rrbracket=0$, on any interior edge $e\in\mathcal{E}_{h,I}$.

In order to do the linearization, we take the Fréchet derivative of $\mathcal{B}(u,v)$ around the exact solution $u$. For the BR1 formulation one might get

[TABLE]

Similarly, for BR2 and SIPG one has

[TABLE]

where the term $F$ is the Fréchet derivative of the penalty term; there holds

[TABLE]

for SIPG and BR2, respectively. For the LDG formulation we have

[TABLE]

Now, let us consider $\mathcal{B}^{\prime}[u](w,\psi)$ for the smooth exact solution of the dual problem, $\psi\in H^{2}(\Omega)$; one might write for BR1, BR2 and SIPG,

[TABLE]

for all $v\in V(h)$, with the following definition of the consistency error for the dual problem (64)

[TABLE]

Similarly, for the LDG formulation one can show

[TABLE]

Presence of a non-zero consistency error for the dual problem reveals the dual inconsistency of the proposed scheme, but as we are going to show that one can obtain the asymptotic dual consistency of the scheme as proposed in [30]; the quasilinear form $\mathcal{B}$ is called asymptotic dual consistent with the target functional $\mathcal{J}:\mathbb{R}\to\mathbb{R}$ if the following holds

[TABLE]

where $u$ and $\psi$ are the exact solutions of the primal and dual problems, (1) and (64) respectively. Here $\mathcal{J}^{\prime}[u]$ is the Fréchet derivative of $\mathcal{J}$, which in our analysis and according to (64) is,

[TABLE]

As we will prove later in section 8.2, such an asymptotic adjoint consistency leads to the optimal convergence rate in ${L_{2}}$ norm. This property has been already investigated e.g., in [33] by numerical tests. Let us present the following lemma for the adjoint consistency error

Lemma 7**.**

Assume that the exact solution of the adjoint problem (64) has elliptic regularity, i.e., (67) holds. Then, there exists $\tilde{C}_{\text{con}}>0$, independent of $h$ but dependent on $q$, such that

[TABLE]

for all $v\in V(h)$.

Proof.

Using assumption A(iii) and Lemmas 2 and 4, for BR2 and SIPG we have

[TABLE]

For LDG there exists additional term which can be bounded similarly as

[TABLE]

Using the boundedness of parameter $\bm{\beta}$ and the fact that

[TABLE]

completes the proof of the lemma with sufficiently large $\tilde{C}_{\text{con}}$. ∎

Using the result of the Lemma 7, smoothness of $\bm{a}_{z}$ and elliptic regularity of $\psi$ yields

[TABLE]

which shows the asymptotic adjoint consistency of the proposed formulations.

6 Stability

In this section we prove the stability of the approximated solution of (54). Let us remark that in this section we do not assume neither the monotonicity nor the Lipschitz continuity of the diffusion operator. The only additional assumption we need is mentioned as the following remark:

Remark 2**.**

Based on the examples in section 1, and our main interest to interpret these problems as nonlinear diffusion phenomenon, it looks natural to assume the following

[TABLE]

which tells that there is no diffusive behavior in the absence of the gradient.

Using (76) and the Taylor’s expansion (3), as well as assumption A(iii) in form of (5), one may write

[TABLE]

We will exploit this result in the rest of this section.

It is well-known that the BR1 method is only weakly-stable (cf. [3]), so we only present the stability result for BR2, SIPG and LDG as the following lemma:

Lemma 8**.**

Let us assume $u_{h}\in V_{h,q}$ is the solution of the primal formulation (54). Then for LDG, BR2 and SIPG methods the following stability result holds

[TABLE]

Here, $C_{\text{F}}$ is the continuity constant of linear form $\mathcal{F}(\cdot)$, that is

[TABLE]

while $C_{\text{co}}>0$ is the following coercivity constant of the operator $\mathcal{B}$ as

[TABLE]

Proof.

For different formulations, we insert the test function $v_{h}=u_{h}$ and investigate if (80) actually holds. Using the fact that $\mathcal{B}(u_{h},u_{h})=\mathcal{F}(u_{h})$, the rest of the proof is complete by the following line

[TABLE]

So the only missing step is the coercivity proof, for which we check $\mathcal{B}(u_{h},u_{h})$ for different formulations:

1. (i)

LDG: From (59) and (52) one has

[TABLE]

From the property (76) we have

[TABLE]

using Young’s inequality with $\delta\in(0,1)$. Now using (30) and the boundedness of $\bm{\beta}$, one has

[TABLE]

for some constant $C$. Then (31) gives the following

[TABLE]

which indicates the condition $\mu_{e}>0$ as the coercivity requirement of LDG formulation. This coincides with the result for the Poisson problems (cf. [3]) and the version of LDG discussed in [31]. 2. (ii)

BR2, SIPG: We combine the proof for these two methods by writing them as a variation of BR1 method. This treatment will be used once more in section 7.1. Using (56), (57) and (58), one may write the following decomposition of primal formulation; for all $v,w\in V(h)$

[TABLE]

where $\mathcal{B}_{\mathrm{BR1}}$ is the primal formulation of BR1 scheme and the next two terms are defined as

[TABLE]

Using (77) one can easily show the positive semi-definiteness of $\mathcal{B}_{\text{BR1}}(u_{h},u_{h})$ as (56). For $T^{(1)}(u_{h},u_{h})$ using the Taylor’s expansion (3) we have

[TABLE]

and for the remaining term $T^{(2)}_{{\mathrm{BR2}},{\mathrm{SIPG}}}$ one can obtain

[TABLE]

Here we have used (31) and (77) in the first and second equation, respectively. Henceforth, we adopt the following notation for a constant $C_{\mathrm{T}}>0$

[TABLE]

Hence, using (85), (77) and (89), one might write

[TABLE]

Applying Young inequality with $0<\delta<1$ and (26), we arrive at

[TABLE]

This leads to the following criterion $N_{l}(\Lambda+\lambda(\dfrac{1}{\delta}-1))\leq C_{\mathrm{T}}$, and since $\delta$ can be arbitrary chosen close to $1$, the margin for stability is $C_{\mathrm{T}}>N_{l}\Lambda$ which can be written separately for BR2 and SIPG as

[TABLE]

This shows that by choosing sufficiently large penalty parameter, as (90) suggests, one can ensure the stability of the solution of SIPG and BR2 formulations.

This completes the proof of coercivity as well the proof of the lemma. ∎

Remark 3**.**

The result of Lemma 8 shows that the stability estimate is available in the degenerate case, i.e. where $\lambda=0$, for SIPG method while for BR2 the uniform lower bound $\lambda$ should be uniformly larger than zero. In the special case when $\bm{a}(u,\nabla u)=a(u)\nabla u$, by the same lines of argument as in [43], one can check that the stability criteria for SIPG remains unchanged while the BR2 one reduces to $\eta_{e}\geq N_{l}$, which is consistent with the result of Poisson problem (cf. [3]).

7 Existence and uniqueness of the discrete solution

In the following sections of this paper, we only consider the BR2 and SIPG formulations. The reason for doing this is the fact that our discussions on the LDG method will follow basically the arguments presented in [11]. Moreover, the BR1 method is not so common due to non-compact stencil and instability [3].

We first prove the strong monotonicity and Lipschitz continuity of the nonlinear operator corresponding to the primal formulation in sections 7.1 and 7.2, respectively. Then we prove the uniqueness of the discrete solution and present a Strang type error estimate.

7.1 Strong monotonicity

Let us start with the following lemma

Lemma 9**.**

If the diffusion operator $\bm{a}$ satisfies the strong monotonicity property as (6), then there exists $C_{\text{SM}}>0$ such that

[TABLE]

for all $w,v\in V(h)$.

Before stating the proof, let us remark the following estimate which later will be used in the analysis:

Lemma 10**.**

For all $w\in V(h)$, the following holds

[TABLE]

for $\eta\geq N_{l}+\frac{1}{2}$.

Proof.

Using Young’s inequality with $0<\delta<1$, (52) and (26), one can write

[TABLE]

Hence, by setting $\delta=1/2$ we arrive at

[TABLE]

Using the definition of energy norm (28), one can add $\pm 1/2|{w}|_{{*},h}^{2}$ to the right hand side and the result is obtained. ∎

Now we are ready to present the proof of Lemma 9:

of Lemma 9.

Using the decomposition (85), and setting $\xi_{h}=v_{h}-w_{h}$ one gets

[TABLE]

Strong monotonicity (6) and (56) readily imply

[TABLE]

For the term $F_{1}$, using (7), (86) and a Young inequality with $\delta>0$, one has

[TABLE]

Using the fact that $-|\nabla_{h}\xi_{h}|^{2}\geq-2|\bm{\theta}_{h}(\xi_{h})|^{2}-2|r(\llbracket\xi_{h}\rrbracket)|^{2}$ and (26), we arrive at

[TABLE]

For the remaining term $F_{2}$, using the definitions (87) and (88), as well (31) and (6), one has

[TABLE]

for SIPG and BR2 methods, respectively. Similar to (89), we define $C_{\mathrm{T}}^{\prime}$ such that, for both methods

[TABLE]

Hence, combining all terms one can write

[TABLE]

Let us set $\eta$ to be a constant larger than $N_{l}+1/2$; then using Lemma 10 and (33) give

[TABLE]

Now consider arbitrary $0<C_{\text{SM}}<\frac{C_{\text{sm}}}{2}$ and set $\delta$ large enough such that the following holds

[TABLE]

The only remaining free parameter is the stabilization parameter $C_{\mathrm{T}}^{\prime}$ and one can choose it sufficiently large, such that

[TABLE]

Finally we arrive at $\mathcal{B}(v_{h},\xi_{h})-\mathcal{B}(w_{h},\xi_{h})\geq C_{\text{SM}}|\!|\!|\xi_{h}|\!|\!|_{h}^{2}$, for both BR2 and SIPG and the proof completes. ∎

7.2 Lipschitz continuity

For the Lipschitz continuity we have the following lemma

Lemma 11**.**

If the diffusion operator $\bm{a}$ satisfies the Lipschitz continuity property as (6), there exists $C_{\text{LC}}<\infty$ independent of the mesh size such that

[TABLE]

for all $z,v,w\in V(h)$.

Proof.

Using the decomposition (85) and similar to Lemma 9 one gets

[TABLE]

where

[TABLE]

It is straightforward to show the Lipschitz continuity of $\mathcal{B}_{\mathrm{BR1}}(v,w)$ using (7). For $F_{1}$, noting the fact that one may write

[TABLE]

and readily

[TABLE]

For BR2 scheme one may note, for any $e\in\mathcal{E}_{h}$

[TABLE]

On the other hand, using the fact the $r^{e}$ vanishes outside two neighbor elements, for the corresponding $F_{2}$ term we have,

[TABLE]

Finally, for SIPG method one can easily write

[TABLE]

Combining the results for $\mathcal{B}_{\mathrm{BR1}}$, $F_{1}$ and $F_{2}$ concludes the proof. ∎

Now we present the existence, uniqueness and stability result for the approximated solution $u_{h}$ as well as a Strang type error estimate. Though we discussed the stability of the solution in section 6, note that the result of Lemma 12 requires stronger condition on the diffusion operator (strong monotonicity and global Lipschitz continuity) than the more general result already discussed in section 6.

Lemma 12**.**

There exists a unique $u_{h}\in V_{h,q}$ solution of (54), which satisfies

[TABLE]

Moreover, the following error estimate holds

[TABLE]

where $u$ is the exact solution of (1).

Proof.

Using the strong monotonicity and Lipschitz continuity (Lemmas 9 and 11) the unique solvability of (54) can be proved by a well-known result as [32, Theorem 3.2.23] or [45, Theorem 35.4], also see [28]. For the stability proof we refer to [11, Theorem 4.5] and here we only present the proof of Strang type error estimate due to its application in the rest of our analysis.

Consider the error $e=u-u_{h}$ and decompose it as $e=\eta+\xi$, where $\eta=u-v_{h}$ and $\xi=v_{h}-u_{h}$ where $v_{h}\in V_{h,q}$. From Lemma 9 one has

[TABLE]

While the first group of terms is the consistency error, applying Lemma 11 gives

[TABLE]

Applying a triangle inequality completes the proof. ∎

8 A priori error estimates

In this section, we provide the error estimate of the BR2 and SIPG methods illustrated in section 4. In section 8.1, we prove the optimal error estimate in the energy norm $|\!|\!|\cdot|\!|\!|_{h}$ using the Strang type error estimate in Lemma 12 and the asymptotic consistency result in section 5.1. In section 8.2, we prove the optimal error estimate in ${L_{2}}$ norm exploiting the result on asymptotic adjoint consistency already provided in section 5.2.

8.1 Energy norm error estimate

Combining the result of Lemmas 12 and 6, gives the following corollary for the error estimate in the energy norm:

Corollary 8.1**.**

Let $u_{h}$ and $u$ be the solution of (54) and the exact solution of (1), respectively. Also assume that $u\in H^{s}(\Omega,\mathcal{T}_{h})$ and $\bm{a}(u,\nabla u)\in H^{s_{*}}(\Omega,\mathcal{T}_{h})$. Then the following holds

[TABLE]

with some $C_{\text{err}}>0$, $\mu=\min(s,q+1)$ and $\mu_{*}=\min(s_{*},q+1)$.

Proof.

Using Lemmas 12 and 6 one can easily write

[TABLE]

Choosing $v_{h}=\pi_{h}u$ and applying the approximation result in Lemma 5 completes the proof with choosing a sufficiently large $C_{\text{err}}$. ∎

Using the error decomposition $e=\eta+\xi$ (as in Lemma 12), we also are interested in obtaining an estimate for $\xi=\pi_{h}u-u_{h}$. By setting $v_{h}=\pi_{h}u$ in (101) and similar to the proof of Corollary 8.1 one might get

[TABLE]

which yields, for some $\tilde{C}_{\text{err}}>0$

[TABLE]

Moreover, using the definition of $\bm{\sigma}_{h}$ and $\bm{\theta}_{h}$ as (53) and (52), one can prove the corresponding lemma for their error estimate

Lemma 13** (Theorem 5.5 in [11]).**

Under the same assumptions as those of Corollary 8.1, there exists $\tilde{C}_{\text{err}}>0$ independent of the mesh size, such that

[TABLE]

and, not necessarily with the same $\tilde{C}_{\text{err}}$,

[TABLE]

where $\mu=\min(s,q+1)$ and $\mu_{*}=\min(s_{*},q+1)$.

The proof follows the same lines as [11] and we skip it here.

8.2 $L_{2}$ norm error estimate

In this section, we present the error estimate of the solution $u_{h}$ in the ${L_{2}}$ norm. Recalling the adjoint problem (64), and by inserting $w=u-u_{h}$ as the (infinite dimensional) test function in (5.2) and (73), to obtain

[TABLE]

where $\mathcal{E}_{p}$ and $\mathcal{E}_{d}$ are defined in (60) and (71) as the consistency errors of the primal and the dual problem, respectively, while $\mathcal{N}(u,u_{h},\psi)$ is the second order linearization error. Here, $\psi_{h}=\pi_{h}\psi\in V_{h,q}$. By the Taylor’s formula in section 3.2 we have (note that $\bm{\theta}_{h}(e)=\nabla_{h}e+r(\llbracket e\rrbracket)$ as already defined)

[TABLE]

First let us handle the error estimate of the last three terms on the right hand side of (8.2), $\mathcal{N}(u,u_{h},\psi)$, $\mathcal{E}_{p}(u,\psi_{h})$ and $\mathcal{E}_{d}(u,u_{h},\psi)$, in the following steps:

1. (i)

Using (109) and very similar arguments as [23, Lemma 3.10] we have

[TABLE]

From the embedding theorem, the last term is bounded since we know $H^{2}(\Omega)\subset W_{4}^{1}(\kappa)$, i.e., $\|\psi\|_{W_{4}^{1}(\Omega,\mathcal{T}_{h})}\leq C\|\psi\|_{H^{2}(\Omega,\mathcal{T}_{h})}$ by a uniform constant $C$. Using (34) gives, since $u-u_{h}\in V(h)$

[TABLE]

Now, the only term to handle (to obtain an additional order of $h$) is $\|u-u_{h}\|_{W_{4}^{1}(\Omega,\mathcal{T}_{h})}$ and $\|r(\llbracket e\rrbracket)\|_{{L_{4}}(\Omega)}$. By writing $u-u_{h}=u-\pi_{h}u+\pi_{h}u-u_{h}=\eta+\xi$, one gets

[TABLE]

In case of $q\geq 2$ and $u\in H^{5/2}(\Omega)$ (note that $q+1>5/2$), application of Lemma 4 gives

[TABLE]

On the other hand using the inverse inequality (35) one has

[TABLE]

Then, application to (105) for term $\|\xi\|_{W_{2}^{1}(\kappa)}$ with $\mu=\min\{s,q+1\}$ and $\mu_{*}=\min\{s_{*},q+1\}$, and quasi-uniformity condition (14) yield

[TABLE]

Hence, $\|\xi\|_{W_{4}^{1}(\Omega,\mathcal{T}_{h})}\leq Ch$ for $u\in H^{5/2}(\Omega,\mathcal{T}_{h})$ and $q\geq 2$, provided that $\bm{a}(u,\nabla u)\in H^{3/2}(\Omega,\mathcal{T}_{h})$. Similarly, since $r(\llbracket e\rrbracket)\in\Sigma_{h,p}$ one can employ (35) to write

[TABLE]

Using (26), (14), Corollary 8.1, and with similar arguments as ((i)) yield

[TABLE]

Combining all terms we have

[TABLE] 2. (ii)

Using Lemma 6 and noticing that $\mathcal{E}_{p}(u,\psi)=0$ (see (60) for smooth $\psi$ as well as the boundary condition of (64)), one has the following upper bound for $\mathcal{E}_{p}(u,\psi_{h})$

[TABLE]

The approximation result of Lemma 4 and $H^{2}$-regularity of $\psi$ give

[TABLE] 3. (iii)

For the adjoint consistency error $\mathcal{E}_{d}(u,e,\psi)$, from Lemma 7 one has

[TABLE]

Combining steps (i)-(iii), with the energy error estimate in Corollary 8.1, one can write

[TABLE]

The remaining term in (8.2) can be bounded using Lipschitz continuity (7) and Lemma 5

[TABLE]

Application to the elliptic regularity of $\psi$ and Corollary 8.1 gives

[TABLE]

which can be summarized as the following lemma for the optimal ${L_{2}}$ error estimate

Lemma 14**.**

Assume that $u_{h}$ and $u$ are the solutions of (54) and the exact solution of (1), respectively. Also assume that $u\in H^{s}(\Omega,\mathcal{T}_{h})$ and $\bm{a}(u,\nabla u)\in H^{s_{*}}(\Omega,\mathcal{T}_{h})$ with $s\geq\frac{5}{2}$ and $s_{*}\geq\frac{3}{2}$. Then, there exists $C_{\text{err}}^{\prime}>0$ such that the following holds

[TABLE]

where $\mu=\min(s,q+1)$ and $\mu_{*}=\min(s_{*},q+1)$ and $q\geq 2$.

9 Conclusion

In this work, we have analyzed different DG formulations of a quasilinear elliptic problem by introducing appropriate numerical flux functions inspired by their original version in linear problems. We showed that in spite of the fact that all of these formulations are inconsistent, they have the asymptotic consistency property for both primal and dual problem. Moreover, we also proved the stability of the solution in ${L_{2}}$ norm under mild assumptions on the problem.

Furthermore, for BR2 and SIPG discretizations, we proved the existence and uniqueness of the discrete solution in case of monotone and globally Lipschitz diffusion operator. Afterwards, under regularity assumptions for the exact solution, we proved the optimal convergence rate in energy norm as well as ${L_{2}}$ norm.