# A Non-gradient DG method for second-order Elliptic Equations in the Non-
divergence Form

Weifeng Qiu  Department of Mathematics, City University of Hong Kong, 83 Tat
Chee Avenue, Hong Kong, China<EMAIL_ADDRESS>,  Jin Ren  Department
of Mathematics and Statistics, Old Dominion University, Norfolk, VA 23529,
USA<EMAIL_ADDRESS>,  Ke Shi  Department of Mathematics and Statistics, Old
Dominion University, Norfolk, VA 23529, USA<EMAIL_ADDRESS>and  Yuesheng Xu
Department of Mathematics and Statistics, Old Dominion University, Norfolk, VA
23529, USA<EMAIL_ADDRESS>
###### Abstract.

$L^{1}$ based optimization is widely used in image denoising, machine learning
and related applications. One of the main features of such approach is that it
naturally provide a sparse structure in the numerical solutions. In this
paper, we study an $L^{1}$ based mixed DG method for second-order elliptic
equations in the non-divergence form. The elliptic PDE in nondivergence form
arises in the linearization of fully nonlinear PDEs. Due to the nature of the
equations, classical finite element methods based on variational forms can not
be employed directly. In this work, we propose a new optimization scheme
coupling the classical DG framework with recently developed $L^{1}$
optimization technique. Convergence analysis in both energy norm and
$L^{\infty}$ norm are obtained under weak regularity assumption. Such $L^{1}$
models are nondifferentiable and therefore invalidate traditional gradient
methods. Therefore all existing gradient based solvers are no longer feasible
under this setting. To overcome this difficulty, we characterize solutions of
$L^{1}$ optimization as fixed-points of proximity equations and utilize matrix
splitting technique to obtain a class of fixed-point proximity algorithms with
convergence analysis. Various numerical examples are displayed to illustrate
the numerical solution has sparse structure with careful choice of the bases
of the finite dimensional spaces. Numerical examples in both smooth and
nonsmooth settings are provided to validate the theoretical results.

###### Key words and phrases:

discontinuous Galerkin, non-gradient optimization, non-divergence form

###### 2000 Mathematics Subject Classification:

65N30, 65L12

W. Qiu’s research is partially supported by the Research Grants Council of the
Hong Kong Special Administrative Region, China. (Project Nos. CityU 11302219,
CityU 11300621). K. Shi is partially supported by National Science Foundation
(Award Number: 2012235). Y. Xu is supported in part by the US National Science
Foundation under grants DMS-1912958 and DMS-220838, and by the US National
Institutes of Health under grant R21CA263876. As a convention the names of the
authors are alphabetically ordered. All authors contributed equally in this
article.

##  1\. Introduction

One of the main features of using $L^{1}$ norm based optimizations is that the
approximation naturally has a sparse representation under specified
transformations or bases. This feature is particularly important in large
scale computations. For instance, in multi-physics applications, one needs to
use the solution from one model as inputs/parameters for the subsequent
physical model. A typical challenge lays in the transition of computation
between models with different scales. Due to the curse of dimensionality, the
computational complexity grows exponantially as the number of independent
inputs/paramters grow linearly. In other words, too many independent inputs
will lead to extremely large system to solve. Therefore, it is desirable to
develop sparsity promoting algorithms such that the input, the solution from
the foremal model can be represented by a small number of basis functions.

In the area of numerical PDEs, numerical schemes based on $L^{1}$ stablization
are less explored. The main challenge is that the $L^{1}$ norm based scheme
will result in minimizing an non-differentiable functional. In [16], the
author studied a finite element method for solving first order PDEs under the
$L^{1}$ norm. The author overcome the difficulty of minimizing the non-
differentiable functional by regularizing the functional with a sequence of
differentiable minimization problems.

In this paper, we propose a numerical scheme with two main features: (1) using
a sparsity promoting functional in the minimization scheme and adopting fixed-
point proximity algorithm [22, 23, 20, 32] with exact convergence analysis of
non-smooth optimization; (2) A sparsity promoting base of the finite element
space [5, 6, 24, 25]. The first condition is met by using $L^{1}$
stabilization while the second is met by using DG finite elements. In [7] the
authors systematically studied the construction of multiscale bases under the
DG setting. It is worth to mention that the current study is an ongoing effort
to develop new numerical schemes such that the numerical solutions with sparse
structure in the finite dimensional space. As a step-stone, we consider the
following elliptic equations in non-divergence form:

$\displaystyle A:D^{2}u$ $\displaystyle=f\quad\text{in}\quad\Omega,$ (1.1)
$\displaystyle u$ $\displaystyle=0,\quad\text{on}\quad\partial\Omega.$ (1.2)

Here the domain $\Omega\in\mathbb{R}^{d}$ $(d=2,3)$ is an open and bounded
polytope and the coefficient matrix
$A=[A_{ij}]\in[L^{\infty}(\Omega)]^{d\times d}$ is a symmetric positive
definite matrix with eigenvalues bounded between $(0,\Lambda)$ with
$\Lambda>0$. $D^{2}u=[\frac{\partial^{2}u}{\partial x_{i}\partial x_{j}}]$
denotes the Hessian matrix of $u$ and
$A:D^{2}u:=\sum_{ij}A_{ij}\frac{\partial^{2}u}{\partial x_{i}\partial x_{j}}$.
Note that $A$ is not necessarily differentiable. Indeed if
$A\in[C^{1}(\Omega)]^{d\times d}$ we may recast the above system into a
divergence form [30]. Hence, the standard steady-state diffusion problem is a
special case of (1.1).

The PDEs in non-divergence form like (1.1) arises in the linearization of
fully nonlinear PDEs such as stochastic control problems [13], Monge-Ampére
equations [27], and the fully nonlinear Hamilton-Jacobi-Bellman equations
[17]. One of the main challenges around this problem is the fact that the
coefficient matrix $A(x)$ in general is not smooth in many practical settings.
In this case, the standard weak variational formulation of the PDE is no
longer available. Indeed, analytical results on the PDE level such as the
existence and uniqueness of solutions for (2.1) are often characterized in the
classical or strong sense.

Due to the possible non-smooth cofficient matrix, developing standard
numerical methods for problem (1.1) with variational form is very difficult.
In [19] the authors study a mixed finite element method based on Galerkin
framework by introducing a so-called “finite element Hessian” and recently in
[9] they extend the idea with discontinuous Galerkin elements. The analysis of
the DG method is provided in [26]. In [11, 12] the authors consider two DG
type methods based on variational form. The key ingredient for the error
estimates is on the discrete Calderon-Zygmond estimates which assumes a full
$H^{2,p}$ regularity for the PDE operator:

$\displaystyle\|w\|_{2,p,\Omega}\leq C_{p}\|A:D^{2}w\|_{p,\Omega},\quad\forall
w\in W_{0}^{1,p}(\Omega)\cap W^{2,p}(\Omega).$ (1.3)

Here $1\leq p<\infty$. This regularity assumption holds if the coefficient
matrix $A\in[C^{0}(\overline{\Omega})]^{d\times d}$ and $\partial\Omega$ is
$C^{1,1}$. In contrast, several methods are studied [34, 14] for equations
with non-smooth coefficient matrix. they considered DG method for
discontinuous $A$ with Cordés condition and convex polyhedral domain $\Omega$.
A primal-dual weak Galerkin method is introduced by Wang and Wang [35].
Roughly speaking, Cordés condition requires that the possibly discontinuous
coefficient matrix $A$ is not too far away from identity matrix multiplied by
a constant. Consequently the difference between the nondivergence form and a
diffusive form is bounded by the $|D^{2}u|$ pointwisely. Recently, a
$L^{p}$-weak Galerkin method for the problem is developed in [4] based on
primal dual weak Galerkin framework. They recast the problem into a Min-Max
optimization problem. All the above mentioned works require a full $W^{2,p}$
regularity assumption (1.3) for the PDE operator in order to establish the
stability of the methods (methods in [34, 14, 35] are valid only for $p=2$).

In a recent work by one of the authors Qiu [30], they develop a first order
system least squares method for (1.1). It is worth to mention that for the
stability of the method, their method only requires $W^{1,2}$ regularity of
the PDE operator:

$\displaystyle\|w\|_{1,2,\Omega}\leq C\|A:D^{2}w\|_{2,\Omega},$ (1.4)

for any $w\in\\{v\in W_{0}^{1,2}(\Omega):A:D^{2}v\in L^{2}(\Omega)\\}$. All
the above mentioned work are formulated and implemented as Galerkin type
schemes. Nochetto and Zhang [29] studied a two-scale method, which is based on
the integro-differential approach and focuses on $L^{\infty}$-error estimates.

In this paper we advocate a different approach to discretize the equation
(1.1) under the minimization framework. Namely, we employ standard DG element
for the primary unknown $u$ and its gradient $\nabla u$. We augment the
minimization problem with $L^{1}$-type penalty terms. Comparing with existing
finite element methods for solving (2.1) based on variational formulations,
our approach has a few distinct features. First, the discrete problem does not
need to use integration by parts which fits the original PDE naturally. No
artificial differential operators are needed. Consequently, the error
estimates is relatively neat and simple. Secondly, we are able to relax the
regularity requirement of the PDE operator from $W^{2,p}$ (1.3) to $W^{1,2}$
(1.4). We only need the exact solution of the PDE $u\in W^{2+\delta,2}$ to
achieve convergence result in the energy norm comparable with above mentioned
existing work. Furthermore, we establish pointwise error estimate based on the
Aleksandrov-Bakelman-Pucci estimates [15], when the domain
$\Omega\in\mathbb{R}^{d}(d=2,3)$ is an open and bounded polytope and the
coefficient matrix $A=A(x)\in[L^{\infty}(\Omega)]^{d\times d}$ is a symmetric
positive definite matrix with eigenvalues bounded between $(0,\Lambda)$ with
$\Lambda>0$ (we need to assume $A$ is piecewise polynomial if $d=3$). Finally,
conventional methods such as least squares or Galerkin type methods based on
variational forms are feasible for PDEs with smooth solutions. Nevertheless,
when the PDE operator is lack of strong regularity such as (1.1) and/or the
desired energy functional is not differentiable ($L^{1}$ for instance), the
weak formulation of a PDE is not application under such settings.

Unlike smooth optimization, the proposed $L^{1}$ optimization cannot be solved
via gradient-based algorithms due to its nondifferentiability. To overcome the
difficulty brought by nonsmoothness, here we adopt the framework of Fixed-
Point Proximity Algorithms (FP2A) originated from image processing and data
science [22, 23, 20, 32] to solve the proposed $L^{1}$ optimization problem.
FP2A first characterizes solutions of $L^{1}$ optimization to be fixed-points
of proximity equations of nonsmooth functions, then a matrix splitting
technique provides a class of fixed-point proximity algorithms. The framework
of firmly nonexpansive operators are then applied to prove the convergence of
FP2A. It is also proven in [20] that FP2A covers popular algorithms such as
first-order primal-dual algorithms and alternating direction method of
multipliers.

The rest of the paper is organized as follows. In Section 2, we present the
details of the mixed DG method for the problem (1.1) with $L^{1}$
stabilization and the main error estimates results. In Section 3 we present
the detailed proofs of the error estimstes. In Section 4, we review some
preliminary approximation results and present the main error estimates of the
method. Given the fact that the scheme is based on $L^{1}$ stabilization, the
implementation of far from trivial, in Section 5, we will present the details
on the iterative algorithm for the scheme and present the convergence analysis
of the algorithm. Several numerical examples will be performed in Section 5
with concluding remarks in Section 6.

##  2\. Mixed formulation for non-divergence form

In this section we present the mixed DG method for the problem (1.1) as a
minimization problem. We begin by specializing (1.1) as a first order system:

$\displaystyle\nabla{\boldsymbol{q}}-u$
$\displaystyle=0\quad\text{in}\quad\Omega,$ (2.1a) $\displaystyle
A:\nabla{\boldsymbol{q}}$ $\displaystyle=f\quad\text{in}\quad\Omega,$ (2.1b)
$\displaystyle u$ $\displaystyle=0,\quad\text{on}\quad\partial\Omega.$ (2.1c)

We define the finite element spaces for the numerical scheme for solving
(2.1). We adopt the notation and norms for the spaces as in [1]. We consider a
family of conforming triangulations ${\mathcal{T}}_{h}$ made of shape-regular
simplexes (triangles in 2D, tetrahedra in 3D). Here conforming means that for
any two element $K,K^{\prime}\in{\mathcal{T}}_{h}$,
$\overline{K}\cap\overline{K^{\prime}}$ is either $\emptyset$, a common
vertex, a common edge or a common face of both elements. Let
${\mathcal{E}}_{h}^{I}$ denote the set of all interior faces of
${\mathcal{T}}_{h}$, and ${\mathcal{E}}_{h}^{B}$ denotes the set of all
boundary faces. We define
${\mathcal{E}}_{h}:={\mathcal{E}}_{h}^{I}\cup{\mathcal{E}}_{h}^{B}$. We use
$h_{K}$ to denote the diameter of the element $K$, and $h_{F}$ the diameter of
the edge $F\in{\mathcal{E}}_{h}^{B}$. The mesh size of ${\mathcal{T}}_{h}$ is
defined as $h:=\max_{K\in{\mathcal{T}}_{h}}h_{K}$. We denote by
${\boldsymbol{n}}_{K}$ the unit outward normal vector on $\partial K$. We will
drop the sub-index $K$ to denote a generic outer normal vector
${\boldsymbol{n}}$. We also introduce the average and jump operators. Let
$F=\partial K\cap\partial K^{\prime}$ be an interior face shared by $K$ and
$K^{\prime}$. $\phi$ is a generic piecewise smooth function (scalar- or
vector-valued). We define the average of $\phi$ on $F$ as

$\\{\\{\phi\\}\\}:=\frac{1}{2}(\phi+\phi^{\prime})$

where $\phi$ and $\phi^{\prime}$ denote the trace of $\phi$ from the interior
of $K$ and $K^{\prime}$ respectively. Furthermore, let $w$ be a piecewise
smooth function and ${\boldsymbol{v}}$ a piecewise smooth vector-valued field.
Analogously, we define the following jumps on $F$:

$\displaystyle\llbracket
w\rrbracket:=w{\boldsymbol{n}}_{K}+w^{\prime}{\boldsymbol{n}}_{K^{\prime}},$
$\displaystyle\llbracket{\boldsymbol{v}}\cdot{\boldsymbol{n}}\rrbracket:={\boldsymbol{v}}\cdot{\boldsymbol{n}}_{K}+{\boldsymbol{v}}^{\prime}\cdot{\boldsymbol{n}}_{K^{\prime}}.$

On a boundary face $F=\partial K\cap\partial\Omega$, we set accordingly

$\\{\\{\phi\\}\\}:=\phi,\quad\llbracket
w\rrbracket:=w{\boldsymbol{n}},\quad\text{and}\quad\llbracket{\boldsymbol{v}}\rrbracket_{N}:={\boldsymbol{v}}\cdot{\boldsymbol{n}}.$

Let $P_{k}(D)$ denote the space of polynomials of degree no more than $k$ over
the domain $D$. We define the following spaces for $k\geq 2$:

$\displaystyle W_{h}$ $\displaystyle=\\{w\in L^{2}(\Omega)|w_{K}\in
P_{k}(K),\forall K\in{\mathcal{T}}_{h}\\},$
$\displaystyle{\boldsymbol{V}}_{h}$
$\displaystyle=\\{{\boldsymbol{v}}\in[L^{2}(\Omega)]^{d}|{\boldsymbol{v}}_{K}\in[P_{k-1}(K)]^{d},\forall
K\in{\mathcal{T}}_{h},\\},$ $\displaystyle W^{c}_{h}$ $\displaystyle=W_{h}\cap
H^{1}_{0}(\Omega),\quad{\boldsymbol{V}}^{c}_{h}={\boldsymbol{V}}_{h}\cap
H(\text{div};\Omega).$

Finally, we use the standard notation for Sobolev norms: $\|\cdot\|_{k,p,D}$
and $|\cdot|_{k,p}$ are used for the standard $W^{k,p}(D)$ norms and semi-
norms, respectively. In addition, for discrete functions we define:

$\|w\|_{k,p,{\mathcal{T}}_{h}}:=\left(\sum_{K\in{\mathcal{T}}_{h}}\|w\|^{p}_{k,p,K}\right)^{\frac{1}{p}},$
$\|w\|_{k,p,{\mathcal{E}}_{h}}:=\left(\sum_{F\in{\mathcal{E}}_{h}}\|w\|^{p}_{k,p,F}\right)^{\frac{1}{p}},$
$\|w\|_{k,p,{\mathcal{E}}^{I}_{h}}:=\left(\sum_{F\in{\mathcal{E}}^{I}_{h}}\|w\|^{p}_{k,p,F}\right)^{\frac{1}{p}},$

We drop the index $k$ when $k=0$ and we drop both indices $k,p$ for the
standard $L^{2}$ norm. i.e. $\|w\|_{p,D}:=\|w\|_{0,p,D}$ and
$\|w\|_{D}:=\|w\|_{0,2,D}$.

We next describe our numerical scheme. For $(w,v)\in
W_{h}\times{\boldsymbol{V}}_{h}$, we introduce the penalty term

$\displaystyle{\mathcal{S}}({\boldsymbol{v}},w):=(h^{-1}\|\llbracket{\boldsymbol{v}}\cdot{\boldsymbol{n}}\rrbracket\|_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
w\rrbracket\|_{1,{\mathcal{E}}_{h}})$ (2.2)

and the energy functional

$\mathcal{J}_{h}(w,{\boldsymbol{v}}):=\|A:\nabla{\boldsymbol{v}}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{v}}-\nabla
w\|^{2}_{{\mathcal{T}}_{h}}+\tau{\mathcal{S}}({\boldsymbol{v}},w),$

where the $h$-independent parameter $\tau$ can be chosen as any positive real
number bigger or equal to $1$. We seek an approximation
$(u_{h},{\boldsymbol{q}}_{h})\in W_{h}\times{\boldsymbol{V}}_{h}$ by solving
the minimization problem

$\displaystyle(u_{h},{\boldsymbol{q}}_{h})$
$\displaystyle=\operatorname*{arg\,min}_{(w,v)\in
W_{h}\times{\boldsymbol{V}}_{h}}\mathcal{J}_{h}(w,{\boldsymbol{v}}).$ (2.3)

The objective function of minimization problem (2.3) is convex and thus, it
has a (global) minimizer. However, since it involves a non-smooth penalty
term, standard gradient type algorithm cannot be applied to this problem.
Fixed-point iterative algorithms will be developed in section 4 for solving
this minimization problem.

######  Remark 2.1.

For non-homogeneous boundary condition case, we should modify the penalty term
(2.2) as:

${\mathcal{S}}({\boldsymbol{v}},w):=\|\llbracket{\boldsymbol{v}}\cdot{\boldsymbol{n}}\rrbracket\|_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
w\rrbracket\|_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|w-g\|_{1,{\mathcal{E}}_{h}^{B}}.$

######  Remark 2.2.

Note that in the energy functional $\mathcal{J}_{h}$ we can replace the
$L^{2}$ norm with any $L^{p}$ norm with $1\leq p<\infty$; the error analysis
should follow almost the same lines as $L^{2}$ since we do not rely on
variational formulations. Nevertheless, we use the $L^{2}$ norm here for the
sake of simplicity in the implementation. As part of our future work, we will
develop efficient algorithms to treat general $L^{p}$ norms in the functional.

######  Remark 2.3.

We believe that the framework can be formulated with the primary variable $u$
only. It should not have substantial difficulty in the analysis of the method.
Nevertheless, for the sake of simplicity in the implementation step, we recast
the problem into the first order system (2.1) and approximate $u$ and
${\boldsymbol{q}}$ together.

In this paper, we provide two types of error estimates of the above scheme. We
first provide the error estimate with the standard discrete $H^{1}$ norms. In
this case we assume that the elliptic equation has the $H^{1}$ regularity: for
any $w\in\\{v\in H_{0}^{1}(\Omega):A:D^{2}v\in L^{2}(\Omega)\\}$,

$\|w\|_{1,2,\Omega}\leq C\|A:D^{2}w\|_{\Omega}.$ (2.4)

It is worth to mention that the above regularity assumption is weaker than the
well-known Cordés condition for the full $H^{2}$ regularity [34]. Namely,
Cordés condition requires $\Omega$ must be convex polyhedral domain and
$A\in[L^{\infty}(\Omega)]^{d\times d}$ is symmetric positive definite and
there exists $\epsilon\in(0,1]$ such that

$\frac{\sum_{i,j}a^{2}_{ij}}{({\bf
Tr}A)^{2}}\leq\frac{1}{n-1+\epsilon}\quad\text{a.e. in}\quad\Omega.$

In contrast, the relaxed $H^{1}$ regularity (2.4) holds in more general
settings. For instance, if $\Omega\in\mathbb{R}^{3}$ is a polyhedral domain
(possibly nonconvex), we can establish (2.4) with similar assumption as Condés
condition as follows:

######  Proposition 2.1.

Let $\Omega\in\mathbb{R}^{3}$ be a polyhedral domain. If
$A\in[L^{\infty}(\Omega)]^{3\times 3}$ is symmetric positive definite and
there exists $\mu>0$ such that

$\|A-\mu\text{Id}\|_{\infty}\leq\mathcal{C}_{\Omega},$

where $\mathcal{C}_{\Omega}$ is a positive constant only depends on the domain
$\Omega$. Then it holds that:

$\|w\|_{1,2,\Omega}\leq C\|A:D^{2}w\|_{\Omega},$

for all $w\in H^{1}_{0}(\Omega)$.

###### Proof.

By triangle inequality, for any $\mu>0$ we have

$\displaystyle\|\mu\Delta w\|_{p,\Omega}$
$\displaystyle\leq\|A:D^{2}w-\mu\Delta
w\|_{p,\Omega}+\|A:D^{2}w\|_{p,\Omega}\leq\|A-\mu\text{Id}\|_{\infty}\|w\|_{2,p,\Omega}+\|A:D^{2}w\|_{p,\Omega}.$

If $\Omega\in\mathbb{R}^{3}$ be a polyhedral domain, in [8] it shows that for
any $p\in[\frac{6}{5},\frac{4}{3})$, there exists $C_{p}>0$ such that

$\|w\|_{2,p,\Omega}\leq C_{p}\|\Delta w\|_{p,\Omega}.$

Inserting this inequality in the previous step we have

$C_{p}^{-1}\mu\|w\|_{2,p,\Omega}\leq\|A-\mu\text{Id}\|_{\infty}\|w\|_{2,p,\Omega}+\|A:D^{2}w\|_{p,\Omega}.$

If there exists a $\mu>0$ such that
$\|A-\mu\text{Id}\|_{\infty}<\frac{1}{2}C_{p}^{-1}\mu\|w\|_{2,p,\Omega}$, we
can have a bound:

$\|w\|_{2,p,\Omega}\leq\frac{2C_{p}}{\mu}\|A:D^{2}w\|_{p,\Omega}.$

Finally, by the Sobolev embedding, if $p=\frac{6}{5}$ with $d=3$ we have

$\|w\|_{1,2,\Omega}\leq
C\|w\|_{2,p,\Omega}\leq\mathcal{C}_{\Omega}\|A:D^{2}w\|_{p,\Omega}.$

Here $\mathcal{C}_{\Omega}$ only depends on the domain $\Omega$. This
completes the proof. ∎

Next we present the first error estimate of the scheme (2.3) which is
summarized as follows:

######  Theorem 2.2.

Let $(u,{\boldsymbol{q}})$ be the exact solution of problem (2.1) that
satisfies the regularity assumption (2.4) and $u\in H^{2+\delta}(\Omega)$ with
$\delta>0$. If $(u_{h},{\boldsymbol{q}}_{h})$ is a minimizer of (2.3), then

$\|u-u_{h}\|_{1,2,{\mathcal{T}}_{h}}+\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\|_{{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$
(2.5)

In addition, if the stronger regularity assumption (1.3) holds for the non-
divergence operator, then we have the optimal convergence with respect to the
$H^{2}$ norm:

$\|u-u_{h}\|_{2,2,{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$
(2.6)

Here $\mathcal{C}_{f}:=\max(\tau^{-2}\|f\|_{\Omega}^{2},1)$.

######  Remark 2.4.

Here we only assume the exact solution has the strong $H^{2+\delta}$
regularity. The linear non-divergence operator is only assumed to have a weak
$H^{1}$ regularity (2.4) for the energy norm error estimate. In addition if we
have full $H^{2}$ regularity for the nondivergence operator, we then obtain
the optimal convergence under the discrete $H^{2}$ norm.

Moreover, we can establish the error estimate with the well-known Aleksandrov-
Bakelman-Pucci estimates (ABP) [15, Chapter 9, Theorem 9.1]: If $w\in
C^{0}(\bar{\Omega})\cap W^{2,d}_{loc}(\Omega)\cap H_{0}^{1}(\Omega)$, then it
holds:

$\|w\|_{\infty,\Omega}\leq C\|A:D^{2}w\|_{d,\Omega}.$ (2.7)

The above estimate holds if the domain $\Omega\in\mathbb{R}^{d}$ $(d=2,3)$ is
an open and bounded polytope, and the coefficient matrix
$A=A(x)\in[L^{\infty}(\Omega)]^{d\times d}$ is a symmetric positive definite
matrix with eigenvalues bounded between $(0,\Lambda)$ with $\Lambda>0$.

######  Theorem 2.3.

Let $(u,{\boldsymbol{q}})$ be the exact solution of the problem (2.1) and
$(u_{h},{\boldsymbol{q}}_{h})$ be a minimizer of (2.3). If $u\in
H^{2+\delta}(\Omega)$ with $\delta>0$, then the following statements hold:

(i) In the 2D case with any polynomial degree $k\geq 2$ we have:

$\displaystyle\|u-u_{h}\|_{\infty,\Omega}$ $\displaystyle\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$

(ii) In the 3D case with $k=2,3$, $\delta>\frac{1}{2}$ and $A$ is piecewise
polynomial with respect to ${\mathcal{T}}_{h}$, we have:

$\|u-u_{h}\|_{\infty,\Omega}\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}-\frac{1}{2}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$

######  Remark 2.5.

In the 3D case, we make two additional restrictions, one on the the polynomial
degree $k=2,3$ and the other on the coefficient matrix $A$. This is due to the
use of the discrete inverse inequality in the analysis. Nevertheless, for the
2D case, this requirement is not needed.

##  3\. Proofs of the error estimates

This section is devoted to providing the proof of the error estimates stated
in the last section. We begin by gathering several auxiliary results to be
used in the proof.

######  Lemma 3.1.

If $(u_{h},{\boldsymbol{q}}_{h})$ is a minimizer of $\mathcal{J}_{h}$ with
$\tau>0$, then

$\tau{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})\leq\|f\|^{2}_{\Omega}$

and

$\tau(h^{-2}\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-4}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}})\leq
2\|f\|^{2}_{\Omega}\mathcal{S}(u_{h},{\boldsymbol{q}}_{h}).$

###### Proof.

Since $(u_{h},{\boldsymbol{q}}_{h})$ is a minimizer of $\mathcal{J}_{h}$, we
observe that

$\tau{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})\leq\mathcal{J}_{h}(u_{h},{\boldsymbol{q}}_{h})\leq\mathcal{J}_{h}(0,{\bf{0}})=\|f\|^{2}_{\Omega}.$

The second inequality is a direct consequence of the first assertion. ∎

For discrete functions, we frequently use the following inverse and trace
inequalities [2]: For $w_{h}\in W(D)$ with $D\in\mathbb{R}^{d}$, we have that

$\displaystyle\|w_{h}\|_{l,p,D}$ $\displaystyle\leq
Ch_{K}^{m-l+\frac{d}{p}-\frac{d}{q}}\|w_{h}\|_{m,q,D},\quad\text{for all
$K\in{\mathcal{T}}_{h}$,}\quad 0\leq m\leq l,1\leq p,q\leq\infty,$ (3.1)
$\displaystyle\|w_{h}\|_{\partial K}$ $\displaystyle\leq
Ch^{-\frac{1}{2}}_{K}\|w_{h}\|_{K},\quad\text{for all
$K\in{\mathcal{T}}_{h}$.}$ (3.2) The global version of the above inequalities
are as follows: for $w_{h}\in W_{h}$ it holds
$\displaystyle\|w_{h}\|_{l,p,{\mathcal{E}}_{h}}$ $\displaystyle\leq
Ch^{m-l+\min\\{0,\frac{d}{p}-\frac{d}{q}\\}}\|w_{h}\|_{m,q,{\mathcal{E}}_{h}},\quad\text{for}\quad
0\leq m\leq l,1\leq p,q\leq\infty,$ (3.3)
$\displaystyle\|w_{h}\|_{{\mathcal{E}}_{h}}$ $\displaystyle\leq
Ch^{-\frac{1}{2}}\|w_{h}\|_{{\mathcal{T}}_{h}},$ (3.4)
$\displaystyle\|\llbracket w_{h}\rrbracket\|_{{\mathcal{E}}_{h}}$
$\displaystyle\leq Ch^{-\frac{1}{2}}\|w_{h}\|_{{\mathcal{T}}_{h}}.$ (3.5)

In the error estimates, the essential step is to bridging the discontinuous
numerical solution $(u_{h},{\boldsymbol{q}}_{h})\in
W_{h}\times{\boldsymbol{V}}_{h}$ with the exact solution
$(u,{\boldsymbol{q}})\in H^{1}(\Omega)\times H(\text{div};\Omega)$ so that we
can apply the regularity of the PDE (2.4) on the error. To this end, we use
several lifting operators on discrete functions as follows.

The first lifting operator is to construct a $u^{c}_{h}\in W^{c}_{h}$ such
that their difference is controlled by the jumps of $u_{h}$. Namely, for an
arbitrary $w_{h}\in W_{h}$, we can use the averaging operator introduced in
[18] to define $w^{c}_{h}$ such that $w^{c}_{h}\in W^{c}_{h}$ and it satisfies
the property that [18][Theorem 2.1]:

$\|w_{h}-w^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}+h^{2}\|\nabla(w_{h}-w^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}\leq
Ch\|\llbracket w_{h}\rrbracket\|^{2}_{{\mathcal{E}}_{h}}.$

With the inverse inequality (3.3) for discrete functions, we have the
following result:

$\|w_{h}-w^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}+h^{2}\|\nabla(w_{h}-w^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}\leq
C\|\llbracket w_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}.$ (3.6)

Likewise, for vector-valued functions
${\boldsymbol{v}}_{h}\in{\boldsymbol{V}}_{h}$, we can construct a
$H(\text{div})$ conforming function
${\boldsymbol{v}}_{h}^{c}\in{\boldsymbol{V}}_{h}\cap H(\text{div};\Omega)$
such that the following approximation property is true:

$\|{\boldsymbol{v}}_{h}-{\boldsymbol{v}}^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}+h^{2}\|\nabla({\boldsymbol{v}}_{h}-{\boldsymbol{v}}^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}\leq
C\|\llbracket{\boldsymbol{v}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}.$
(3.7)

Notice that $u^{c}_{h}\in W^{c}_{h}$ is only in $H^{1}(\Omega)$ which is not
sufficient for the regularity condition $w\in H^{1}_{0}(\Omega),D^{2}w\in
L^{2}(\Omega)$. Hence, we need further lift the function to a function in
$H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$. We will need two different lifting
operators for the $H^{1}$ and $L^{\infty}$ error estimates respectively.

The first one is due to [3]. They established the following result:

######  Lemma 3.2.

For any function $w^{c}_{h}\in W^{c}_{h}$, there exists
$\widetilde{w}^{c}_{h}\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ such that:

$\|D^{l}(w^{c}_{h}-\widetilde{w}^{c}_{h})\|^{2}_{\Omega}\leq
Ch^{3-2l}\|\llbracket\nabla
w^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}^{I}_{h}},\quad
l=0,1,2.$

It is worth to mention that here the construction of $\widetilde{w}^{c}_{h}$
is done in the context of conforming virtual element spaces. This leads to the
fact that $\widetilde{w}^{c}_{h}$ in general is not a piecewise polynomial in
$W_{h}$. Nevertheless, our analysis for the first error estimate does not
require the lifted function in the discrete space $W_{h}$.

The next auxiliary result is on the relation between
$(w_{h},{\boldsymbol{v}}_{h})$ and its conforming counterpart
$(w^{c}_{h},{\boldsymbol{v}}^{c}_{h})$.

######  Lemma 3.3.

If $(w^{c}_{h},{\boldsymbol{v}}^{c}_{h})\in
W^{c}_{h}\times{\boldsymbol{V}}^{c}_{h}$ is the lifting of
$(w_{h},{\boldsymbol{v}}_{h})\in W_{h}\times{\boldsymbol{V}}_{h}$, then

$\displaystyle\|\llbracket\nabla
w^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}_{h}^{I}}$
$\displaystyle\leq Ch^{-1}\|{\boldsymbol{v}}^{c}_{h}-\nabla
w^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}},$
$\displaystyle\|{\boldsymbol{v}}^{c}_{h}-\nabla
w^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}$ $\displaystyle\leq
C(\|{\boldsymbol{v}}_{h}-\nabla
w_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\llbracket{\boldsymbol{v}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
w_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}).$

###### Proof.

The first assertion is due to the fact that
${\boldsymbol{v}}_{h}^{c}\in{\boldsymbol{V}}^{c}_{h}$ and the trace inequality
(3.5):

$\|\llbracket\nabla
w^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}_{h}^{I}}=\|\llbracket({\boldsymbol{v}}^{c}_{h}-\nabla
w^{c}_{h})\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}_{h}^{I}}\leq
Ch^{-1}\|({\boldsymbol{v}}^{c}_{h}-\nabla
w^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}.$

For the second assertion we begin with triangle inequality:

$\displaystyle\|{\boldsymbol{v}}^{c}_{h}-\nabla
w^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}$ $\displaystyle\leq
C(\|{\boldsymbol{v}}_{h}-\nabla
w_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|{\boldsymbol{v}}^{c}_{h}-{\boldsymbol{v}}_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\nabla(w^{c}_{h}-w_{h})\|^{2}_{{\mathcal{T}}_{h}})$
$\displaystyle\leq C(\|{\boldsymbol{v}}_{h}-\nabla
w_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\llbracket{\boldsymbol{v}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
w_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}),$

the last step is due to the approximation property of the lifting operators
(3.6), (3.7). ∎

We are ready to present the proof of Theorem 2.2:

###### Proof.

We begin by proving

$\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq
C\max(\tau^{-2}\|f\|_{\Omega}^{2},1)\mathcal{J}_{h}(u_{h},{\boldsymbol{q}}_{h}).$
(3.8)

To this end, we let $(u^{c}_{h},{\boldsymbol{q}}^{c}_{h})\in
W^{c}_{h}\times{\boldsymbol{V}}^{c}_{h}$ be the lifting of
$(u_{h},{\boldsymbol{q}}_{h})$. In addition, we let
$\widetilde{u}^{c}_{h}\in\widetilde{W}^{c}_{k+2}$ be the $C^{1}$-lifting of
$u^{c}_{h}$, and set

$T_{1}:=C\|u-\widetilde{u}^{c}_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}},\
T_{2}:=C\|\widetilde{u}^{c}_{h}-u^{c}_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}},\
T_{3}:=C\|u^{c}_{h}-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}.$

We use the triangle inequality to obtain

$\displaystyle\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq T_{1}+T_{2}+T_{3}.$

By the approximation property of the lifting operator (3.6), we have that

$T_{3}\leq Ch^{-2}\|\llbracket u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}.$

Similarly, by Lemma 3.2, Lemma 3.3 and inverse inequality (3.3), we can bound
$T_{2}$ as:

$\displaystyle T_{2}$ $\displaystyle\leq
Ch\|\llbracket\nabla{u}^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}^{I}_{h}}\leq
C\|{\boldsymbol{q}}^{c}_{h}-\nabla u^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}$
$\displaystyle\leq C(\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}})$

For $T_{1}$, we recall the regularity assumption (2.4) to have that

$T_{1}\leq
C\|A:D^{2}\widetilde{u}^{c}_{h}-A:D^{2}u\|^{2}_{\Omega}=C\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega}.$
(3.9)

Introducing the notation

$T_{11}:=\|A:(\nabla{\boldsymbol{q}}_{h}-\nabla{\boldsymbol{q}}^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}},\
T_{12}:=\|A:(\nabla{\boldsymbol{q}}^{c}_{h}-D^{2}u^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}},$

and

$T_{13}:=\|A:(D^{2}u^{c}_{h}-D^{2}\widetilde{u}^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}),$

by the triangle inequality, we further split the right hand side of equation
(3.9) as

$\displaystyle T_{1}\leq
C(\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+T_{11}+T_{12}+T_{13}).$

For $T_{11}$, by the lifting property (3.7) we have:

$T_{11}\leq
C\|\nabla({\boldsymbol{q}}_{h}-{\boldsymbol{q}}^{c}_{h})\|^{2}_{{\mathcal{T}}_{h}}\leq
h^{-2}\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}.$

For $T_{12}$, with inverse inequality (3.3), we have:

$\displaystyle T_{12}$ $\displaystyle\leq
C\|\nabla{\boldsymbol{q}}^{c}_{h}-D^{2}u^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}\leq
Ch^{-2}\|{\boldsymbol{q}}^{c}_{h}-\nabla u^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}$
then with Lemma 3.3, we further have: $\displaystyle T_{12}$
$\displaystyle\leq C(h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-4}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}})$

For $T_{13}$, by virtue of Lemma 3.2 we have:

$T_{13}\leq Ch^{-1}\|\llbracket\nabla
u^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}_{h}^{I}}.$

We can bound the right side by Lemma 3.3 to have:

$T_{13}\leq Ch^{-2}(\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}})$

Now if we combine the above estimates for $T_{1},T_{2},T_{3}$ we have:

$\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq
C(\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-4}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}).$

By virtue of Lemma 3.1, we have:

$\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq
C(\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+2\tau^{-1}\|f\|^{2}_{\Omega}{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})).$

We notice that

$\displaystyle\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+2\tau^{-1}\|f\|^{2}_{\Omega}{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})$
$\displaystyle\leq$
$\displaystyle\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\max(\dfrac{2\|f\|_{\Omega}^{2}}{\tau^{2}},1)\tau{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})$
$\displaystyle\leq$
$\displaystyle\max(\dfrac{2\|f\|_{\Omega}^{2}}{\tau^{2}},1)(\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\tau{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h})).$

Then we have:

$\displaystyle\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq$ $\displaystyle
C\mathcal{C}_{f}(\|A:\nabla{\boldsymbol{q}}_{h}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\tau{\mathcal{S}}(u_{h},{\boldsymbol{q}}_{h}))$
$\displaystyle=$ $\displaystyle
C\mathcal{C}_{f}{\mathcal{J}}_{h}(u_{h},{\boldsymbol{q}}_{h}).$

This completes the proof of (3.8). Finally, we recall that
$(u_{h},{\boldsymbol{q}}_{h})$ is the minimizer of the functional
${\mathcal{J}}_{h}$ (2.3). This implies that:

$\displaystyle\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}{\mathcal{J}}_{h}(u_{h},{\boldsymbol{q}}_{h})\leq
C\mathcal{C}_{f}{\mathcal{J}}_{h}(w_{h},{\boldsymbol{v}}_{h})$

for all $(w_{h},{\boldsymbol{v}}_{h})\in W_{h}\times{\boldsymbol{V}}_{h}$. Let
$\Pi_{W}u$ be the Lagrange interpolant of $u$ over $W^{c}_{h}$ and
$\boldsymbol{\Pi}_{V}{\boldsymbol{q}}$ be the BDM projection of
${\boldsymbol{q}}$ in ${\boldsymbol{V}}_{h}^{c}$. By taking
$(w_{h},{\boldsymbol{v}}_{h})=(\Pi_{W}u,\boldsymbol{\Pi}_{V}{\boldsymbol{q}})$
we have that ${\mathcal{S}}(\Pi_{W}u,\boldsymbol{\Pi}_{V}{\boldsymbol{q}})=0$
since they are in the comforming subspaces
$W^{c}_{h}\times{\boldsymbol{V}}^{c}_{h}$. The above estimate implies:

$\displaystyle\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}{\mathcal{J}}_{h}(u_{h},{\boldsymbol{q}}_{h})\leq
C\mathcal{C}_{f}{\mathcal{J}}_{h}(\Pi_{W}u,\boldsymbol{\Pi}_{V}{\boldsymbol{q}})$
$\displaystyle=$ $\displaystyle
C\mathcal{C}_{f}(\|A:\nabla\boldsymbol{\Pi}_{V}{\boldsymbol{q}}-f\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|\boldsymbol{\Pi}_{V}{\boldsymbol{q}}-\nabla\Pi_{W}u\|^{2}_{{\mathcal{T}}_{h}}+\tau{\mathcal{S}}(\Pi_{W}u,\boldsymbol{\Pi}_{V}{\boldsymbol{q}})$
$\displaystyle=$ $\displaystyle
C\mathcal{C}_{f}(\|A:\nabla(\boldsymbol{\Pi}_{V}{\boldsymbol{q}}-{\boldsymbol{q}})\|^{2}_{{\mathcal{T}}_{h}}+h^{-2}\|(\boldsymbol{\Pi}_{V}{\boldsymbol{q}}-{\boldsymbol{q}})-\nabla(\Pi_{W}u-u)\|^{2}_{{\mathcal{T}}_{h}})$
$\displaystyle\leq$ $\displaystyle
C\mathcal{C}_{f}h^{2\min\\{k-1,\delta\\}}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega}.$

The last step we used the standard approximation properties of the
interpolations $(\Pi_{W}u,\boldsymbol{\Pi}_{V}{\boldsymbol{q}})$. This
completes the proof of the error estimates for
$\|u-u_{h}\|_{1,2,{\mathcal{T}}_{h}}$. For the error in ${\boldsymbol{q}}$,
notice that (3.8) implies that:

$\|u-u_{h}\|^{2}_{1,2,{\mathcal{T}}_{h}}+h^{-2}\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}\mathcal{J}_{h}(u_{h},{\boldsymbol{q}}_{h}).$

With the same argument as the above, we have:

$h^{-1}\|{\boldsymbol{q}}_{h}-\nabla u_{h}\|_{{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$

By triangle inequality, we have:

$\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\|_{{\mathcal{T}}_{h}}\leq\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|_{{\mathcal{T}}_{h}}+\|\nabla u-\nabla u_{h}\|_{{\mathcal{T}}_{h}}\leq
C\mathcal{C}_{f}h^{\min\\{k-1,\delta\\}}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$

This completes the proof for (2.5).

Finally, if we assume the full $H^{2}$ regularity (1.3), the above proofs can
be seamlessly modified to get the optimal convergence rate for the discrete
$H^{2}$ norm. Namely, we can use the same splitting:

$\displaystyle\|u-u_{h}\|^{2}_{2,2,{\mathcal{T}}_{h}}$ $\displaystyle\leq
C(\|u-\widetilde{u}^{c}_{h}\|^{2}_{2,2,{\mathcal{T}}_{h}}+\|\widetilde{u}^{c}_{h}-u^{c}_{h}\|^{2}_{2,2,{\mathcal{T}}_{h}}+\|u^{c}_{h}-u_{h}\|^{2}_{2,2,{\mathcal{T}}_{h}})$
$\displaystyle\leq\widetilde{T}_{1}+Ch^{-2}T_{2}+Ch^{-2}T_{3}.$

Here we applied the standard dicsrete inverse inequality (3.1). This extra
$h^{-2}$ power is not going to degenerate the order of convergence since in
the above estimate $T_{1}$ is the dominating error while $T_{2},T_{3}$ are
$h^{2}$ higher order terms.

For $\widetilde{T}_{1}$, by virtue of the full $H^{2}$ regularity (1.3), we
can bound this term as:

$\widetilde{T}_{1}\leq
C\|A:D^{2}\widetilde{u}^{c}_{h}-A:D^{2}u\|^{2}_{\Omega}=C\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega}$

and the rest estimate is exactly the same as $T_{1}$. Therefore, the final
convergence rate remains the same as the $H^{1}$ norm without losing any
order. This completes the proof of (2.6).

∎

For the second error estimates Theorem 2.3, the main difference from the above
proof lays in the fact that we relax the regularity assumption to (2.7):

$\|w\|_{\infty,\Omega}\leq C\|A:D^{2}w\|_{d,\Omega}.$

For $d=3$, note that the RHS of the above ABP estimates become a stronger
$L^{3}$ norm of the source term. This forces us to use a discrete inverse
inequality in the analysis. Consequently, we need to use another lifting
operator such that the resulting function remains to be piecewise polynomials.
Namely, for any function $w_{h}\in W^{c}_{h}$, let $\widetilde{W}^{c}_{h}$ to
be the classical $C^{1}$-conforming Hsieh-Clough-Tocher macroelement [10] of
degree $k+2$. In [28] they prove that for $k\geq 2$ with $d=2$ and $k=2,3$
with $d=3$, there exists a $\widehat{w}^{c}_{h}\in\widetilde{W}^{c}_{h}$ such
that the same approximation property as in Lemma 3.2 holds:

$\|D^{l}(w^{c}_{h}-\widetilde{w}^{c}_{h})\|^{2}_{\Omega}\leq
Ch^{3-2l}\|\llbracket\nabla
w^{c}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{{\mathcal{E}}^{I}_{h}},\quad
l=0,1,2.$

It is worth to mention that the use of virtual element lifting in Lemma 3.2 is
necessary since the above Hsieh-Clough-Tocher macroelement lifting is only
defined for $k=2,3$ in 3D.

Next we define the $H^{2}$ conforming projection of $u$ in two steps:

######  Definition 3.4.

Let $u\in H^{2+\delta}(\Omega)\cap H^{1}_{0}(\Omega)$. We use $\Pi_{h}u\in
W^{c}_{h}$ to denote the standard Lagrange interpolant of $u$.
$\widetilde{\Pi}_{h}^{c}u$ is the $H^{2}$ projection of $\Pi^{c}_{h}u$ by
using the Hsieh-Clough-Tocher macroelement lifting.

We have the following approxmation property of $\widetilde{\Pi}_{h}^{c}u$:

######  Lemma 3.5.

Let $u\in H^{2+\delta}(\Omega)\cap H^{1}_{0}(\Omega)$ and
$\widetilde{\Pi}_{h}^{c}u\in\widetilde{W}^{c}_{h}\cap H^{1}_{0}(\Omega)$ be
the $H^{2}$ interpolant of $u$ defined as above. We have:

$\|u-\widetilde{\Pi}_{h}^{c}u\|_{l,2,\Omega}\leq
Ch^{\min\\{k+1,2+\delta\\}-l}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega},\quad
l=0,1,2.$

###### Proof.

By the triangle inequality, standard interpolant property and Lemma 3.2 we
have:

$\displaystyle\|u-\widetilde{\Pi}_{h}^{c}u\|_{l,2,\Omega}$
$\displaystyle\leq\|u-\Pi_{h}^{c}u\|_{l,2,\Omega}+\|\Pi_{h}^{c}u-\widetilde{\Pi}_{h}^{c}u\|_{l,2,\Omega}$
$\displaystyle\leq
Ch^{\min\\{k+1,2+\delta\\}-l}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}+Ch^{\frac{3}{2}-l}\|\llbracket\nabla\Pi^{c}_{h}u\cdot{\boldsymbol{n}}\rrbracket\|_{{\mathcal{E}}^{I}_{h}}.$

The second term on the right hand side can be written as:

$\displaystyle
Ch^{\frac{3}{2}-l}\|\llbracket\nabla\Pi^{c}_{h}u\cdot{\boldsymbol{n}}\rrbracket\|_{{\mathcal{E}}^{I}_{h}}$
$\displaystyle=Ch^{\frac{3}{2}-l}\|\llbracket\nabla(u-\Pi^{c}_{h}u)\cdot{\boldsymbol{n}}\rrbracket\|_{{\mathcal{E}}^{I}_{h}}$
$\displaystyle\leq
Ch^{\min\\{k+1,2+\delta\\}-l}\|u\|_{\min\\{k+1,2+\delta\\},2,\Omega}.$

The last step is due to the trace inequality and the approximation property of
Lagrange interpolant. This completes the proof. ∎

###### Proof.

of Theorem 2.3 The ABP regularity (2.7) now reads as:

$\|w\|_{\infty,\Omega}\leq C\|A:D^{2}w\|_{\Omega},\quad\forall w\in
C^{0}(\overline{\Omega})\cap W_{loc}^{2,2}(\Omega)\cap H_{0}^{1}(\Omega).$
(3.10)

Similar as in the above proof, we split the error into three terms:

$\displaystyle\|u-u_{h}\|^{2}_{\infty,\Omega}$ $\displaystyle\leq
C(\|u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega}+\|\widetilde{u}^{c}_{h}-u^{c}_{h}\|^{2}_{\infty,\Omega}+\|u^{c}_{h}-u_{h}\|^{2}_{\infty,\Omega})$
(3.11) $\displaystyle:=T_{1}+T_{2}+T_{3}.$

For $T_{2},T_{3}$, we can use the same arguments as in the above proof to
obtain:

$\displaystyle T_{2}$ $\displaystyle\leq
Ch^{-d}\|\widetilde{u}^{c}_{h}-u^{c}_{h}\|^{2}_{{\mathcal{T}}_{h}}\leq
Ch^{-2}(\|{\boldsymbol{q}}_{h}-\nabla
u_{h}\|^{2}_{{\mathcal{T}}_{h}}+\|\llbracket{\boldsymbol{q}}_{h}\cdot{\boldsymbol{n}}\rrbracket\|^{2}_{1,{\mathcal{E}}^{I}_{h}}+h^{-2}\|\llbracket
u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}),$ $\displaystyle T_{3}$
$\displaystyle\leq Ch^{-d}\|u^{c}_{h}-u_{h}\|^{2}_{{\mathcal{T}}_{h}}\leq
Ch^{-4}\|\llbracket u_{h}\rrbracket\|^{2}_{1,{\mathcal{E}}_{h}}.$

For $T_{1}$, we further split it by inserting the projection of $u$ defined in
Lemma 3.5:

$\displaystyle T_{1}$ $\displaystyle\leq
C(\|u-\widetilde{\Pi}^{c}_{h}u\|^{2}_{\infty,\Omega}+\|\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega})$
$\displaystyle\leq
C(\|u-\widetilde{\Pi}^{c}_{h}u\|^{2}_{2,2,\Omega}+\|\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega}).$

We first consider the case $d=2$. Since
$\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\in\widetilde{W}^{c}_{h}\cap
H^{1}_{0}(\Omega)\subset C^{0}(\bar{\Omega})\cap W^{2,2}_{loc}(\Omega)$, we
can use the ABP condition (3.10) to obtain:

$\displaystyle\|\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega}\leq$
$\displaystyle
C\|A:D^{2}(\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h})\|^{2}_{\Omega}$
$\displaystyle\leq$ $\displaystyle
C\|A:D^{2}(u-\widetilde{\Pi}^{c}_{h}u)\|^{2}_{\Omega}+C\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega}.$

So overall we have:

$\displaystyle T_{1}$ $\displaystyle\leq
C(\|u-\widetilde{\Pi}^{c}_{h}u\|^{2}_{2,2,\Omega}+\|A:D^{2}(u-\widetilde{\Pi}^{c}_{h}u)\|^{2}_{\Omega}+\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega})$
$\displaystyle\leq
Ch^{2\min\\{k-1,\delta\\}}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega}+C\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega}.$

At this point, the second term on the right hand side in the above estimate is
the same as in the estimates of $T_{1}$ in the proof of Theorem 2.1.
Therefore, we can proceed in the same way from here to have the final
estimates as:

$\displaystyle\|u-u_{h}\|^{2}_{\infty,\Omega}\leq$ $\displaystyle
Ch^{2\min\\{k-1,\delta\\}}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega}+C\mathcal{C}_{f}\mathcal{J}_{h}(u_{h},{\boldsymbol{q}}_{h})$
$\displaystyle\leq$ $\displaystyle
C\mathcal{C}_{f}h^{2\min\\{k-1,\delta\\}}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega},$

if the exact solution is smooth enough $u\in H^{2+\delta}(\Omega)$ for some
$\delta>0$.

Finally we consider the case $d=3$. The burden for the sacrifice in the
convergence rate lays in the ABP estimates. Namely, now the ABP estimate (2.7)
becomes:

$\|w\|_{\infty,\Omega}\leq C\|A:D^{2}w\|_{3,\Omega},\quad\forall w\in
C^{0}(\overline{\Omega})\cap W_{\text{loc}}^{2,3}(\Omega)\cap
H_{0}^{1}(\Omega).$ (3.12)

The $L^{3}$ norm on the right side is stronger than the one we have in the
minimizer. We bypass this issue with standard discrete inverse inequality
which degenerates the order of convergence in the final stage. We lay out the
details of the proof below.

Begin with the spitting as in (3.11) and $T_{1}$. The only term needs
modification is
$\|\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega}$. More
precisely, we can apply the ABP estimate as:

$C\|\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h}\|^{2}_{\infty,\Omega}\leq
C\|A:D^{2}(\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h})\|^{2}_{3,\Omega}.$

Notice that if $A$ is piecewise polynomial, the last term in the above
estimate $A:D^{2}(\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h})$ remains to
be piecewise polynomials over ${\mathcal{T}}_{h}$, we can now apply the
inverse inequality (3.1) to obtain:

$\displaystyle\|A:D^{2}(\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h})\|^{2}_{3,\Omega}$
$\displaystyle\leq
Ch^{-1}\|A:D^{2}(\widetilde{\Pi}^{c}_{h}u-\widetilde{u}^{c}_{h})\|^{2}_{\Omega}$
$\displaystyle\leq
Ch^{-1}(\|A:D^{2}(\widetilde{\Pi}^{c}_{h}u-u)\|^{2}_{\Omega}+\|A:D^{2}\widetilde{u}^{c}_{h}-f\|^{2}_{\Omega})$

From here we can use the same argument to control the last term on the right
hand side as in the previous case. Finally, due to the $h^{-1}$ from the
inverse inequality, we have the final estimates as:

$\displaystyle\|u-u_{h}\|^{2}_{{\mathcal{T}}_{h}}$ $\displaystyle\leq
Ch^{2\min\\{k-1,\delta\\}-1}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega}+C\mathcal{C}_{f}h^{-1}\mathcal{J}_{h}(u_{h},{\boldsymbol{q}}_{h})$
$\displaystyle\leq
C\mathcal{C}_{f}h^{2\min\\{k-1,\delta\\}-1}\|u\|^{2}_{\min\\{k+1,2+\delta\\},2,\Omega},$

if the exact solution is smooth enough $u\in H^{2+\delta}(\Omega)$ for some
$\delta>\frac{1}{2}$. This completes the proof.

∎

##  4\. Fixed-Point Proximity Algorithms

In this section, we develop a fast numerical algorithm for optimization
problem (2.3) based on the Fixed-Point Proximity Algorithms (FP2A) [22, 20,
32] and establish a convergence theorem for the algorithm.

A challenge of solving problem (2.3) is that the $L_{1}$ norm involved in its
objective function is non-differentiable, making classical gradient-based
methods invalid. To surmount this difficulty, we first characterize solutions
of the optimization problem (2.3) as fixed-points of a nonlinear map involved
proximity operators of functions appearing in the objective function. Based on
the fixed-point equation, application of a matrix splitting technique to the
fixed-point equation leads to the FP2A with guaranteed convergence.

###  4.1. An Equivalent Formulation of the Basic Optimization Problem.

We first reformulate the optimization problem (2.3) in a form convenient to
apply FP2A. To this end, we choose a basis of the finite dimensional space
$W_{h}\times{\boldsymbol{V}}_{h}$ and denote it by $\mathcal{B}_{W_{h}\times
V_{h}}$. The numerical solution $(u_{h},{\boldsymbol{q}}_{h})$ is represented
by its coefficients $x$ with respect to this basis. Note the choice of the
basis is crucial to obtain a sparse representation of the solution. We will
demonstrate this point with two choices of bases: a standard Lagrange nodal
basis and a multiscale basis.

We associate each element in $(w,{\boldsymbol{v}})\in
W_{h}\times{\boldsymbol{V}}_{h}$ a vector $x\in\mathbb{R}^{N}$, which
represents its coefficients in the basis. This association is one-to-one.
Hence, $W_{h}\times{\boldsymbol{V}}_{h}$ is isomorphic to $\mathbb{R}^{N}$.
Specifically, we use $\\{w_{i}\\}_{i=1}^{N_{W}}$ to denote the basis of
$W_{h}$ and $\\{{\boldsymbol{v}}_{i}\\}_{i=1}^{N_{V}}$ the basis of
${\boldsymbol{V}}_{h}$. Let

$\displaystyle B_{11}$ $\displaystyle:=h^{-2}\begin{bmatrix}{\left(\nabla
w_{i},\nabla w_{j}\right)}_{{\mathcal{T}}_{h}}\end{bmatrix}_{1\leq i,j\leq
N_{W}},$ $\displaystyle B_{12}$
$\displaystyle:=h^{-2}\begin{bmatrix}{\left(\nabla
w_{i},{\boldsymbol{v}}_{j}\right)}_{{\mathcal{T}}_{h}}\end{bmatrix}_{1\leq
i\leq N_{W},1\leq j\leq N_{V}},$ $\displaystyle B_{22}$
$\displaystyle:=\begin{bmatrix}{\left(A:\nabla{\boldsymbol{v}}_{i},A:\nabla{\boldsymbol{v}}_{j}\right)}_{{\mathcal{T}}_{h}}+h^{-2}{\left({\boldsymbol{v}}_{i},{\boldsymbol{v}}_{j}\right)}_{{\mathcal{T}}_{h}}\end{bmatrix}_{1\leq
i,j\leq N_{V}},$

and define

$B:=\begin{bmatrix}B_{11}&B_{12}\\\ B_{12}^{\ast}&B_{22}\end{bmatrix}.$

We also define

$b:=\begin{bmatrix}{\left(-2{\left(A:\nabla{\boldsymbol{v}}_{j},f\right)}_{{\mathcal{T}}_{h}}\right)}_{1\leq
j\leq N_{V}}\\\ 0_{N_{W}}\end{bmatrix},$

where $0_{N_{W}}$ denotes the $N_{W}$-dimensional zero vector. In terms of the
symmetric positive semi-definite matrix $B$ and the vector $b$ defined above,
we write

${\left\|A:\nabla{\boldsymbol{v}}-f\right\|}_{{\mathcal{T}}_{h}}^{2}+h^{-2}{\left\|{\boldsymbol{v}}+\nabla
w\right\|}_{{\mathcal{T}}_{h}}^{2}-{\left\|f\right\|}^{2}_{{\mathcal{T}}_{h}}=x^{\ast}Bx+b^{\ast}x.$

Likewise, we rewrite the boundary terms in (2.3). To this end, we introduce
matrix $L\in\mathbb{R}^{M\times N}$ by

$L:=\begin{bmatrix}\tau h^{-1}L_{1}\\\ \tau h^{-1}L_{2}\\\ \tau
h^{-2}L_{3}\end{bmatrix}$ (4.1)

with $L_{1}\in\mathbb{R}^{2|{\mathcal{E}}_{h}^{I}|\times N}$,
$L_{2}\in\mathbb{R}^{3|{\mathcal{E}}_{h}^{I}|\times N}$ and
$L_{3}\in\mathbb{R}^{3|{\mathcal{E}}_{h}^{B}|\times N}$. Here, $L_{1}x$
contains all endpoint values of
${\left\llbracket{\boldsymbol{v}}\cdot{\boldsymbol{n}}\right\rrbracket}$ on
${\mathcal{E}}_{h}^{I}$, $L_{2}x$ contains all endpoint and midpoint values of
${\left\llbracket w\right\rrbracket}$ on ${\mathcal{E}}_{h}^{I}$ and $L_{3}x$
contains all endpoint and midpoint values of $w$ on ${\mathcal{E}}_{h}^{B}$.
We then introduce $d\in\mathbb{R}^{M}$ by setting

$d:=\begin{bmatrix}0_{2|{\mathcal{E}}_{h}^{I}|}\\\
0_{3|{\mathcal{E}}_{h}^{I}|}\\\ d_{1}\end{bmatrix},$

where $d_{1}\in\mathbb{R}^{{}_{2|{\mathcal{E}}_{h}^{B}|}}$ contains all
endpoint and midpoint values of $g$ on ${\mathcal{E}}_{h}^{B}$, with each
element arranged according to $L_{3}x$. In terms of the matrix $L$ and the
vector $d$, we identify

${\tau
h^{-1}{\left\|{\left\llbracket{\boldsymbol{v}}\cdot{\boldsymbol{n}}\right\rrbracket}\right\|}_{1,{\mathcal{E}}^{I}_{h}}+\tau
h^{-1}{\left\|{\left\llbracket
w\right\rrbracket}\right\|}_{1,{\mathcal{E}}^{I}_{h}}+\tau
h^{-2}{\left\|w-g\right\|}_{1,{\mathcal{E}}^{B}_{h}}}={\left\|Lx-d\right\|}_{1}.$

In this notation, problem (2.3) is reformulated as

$\operatorname*{arg\,min}_{x\in\mathbb{R}^{N}}{\left\\{x^{\ast}Bx+b^{\ast}x+{\left\|Lx-d\right\|}_{1}\right\\}}.$
(4.2)

For convenience of further analysis, we define functions
$g:\mathbb{R}^{N}\to\mathbb{R}$ by

$g(x)\coloneqq x^{\ast}Bx+b^{\ast}x,\ \ \mbox{for}\ \ x\in\mathbb{R}^{N}$

and $h:\mathbb{R}^{M}\to\mathbb{R}$ by

$h(y)\coloneqq\|y-d\|_{1},\ \ \mbox{for}\ \ y\in\mathbb{R}^{M}.$

With these definitions, we identify (4.2) as

$\operatorname*{arg\,min}_{x\in\mathbb{R}^{N}}{\left\\{g{\left(x\right)}+h{\left(Lx\right)}\right\\}}.$
(4.3)

Clearly, (4.3) is a convex optimization problem studied in [20, 22]. Since the
objective function of (4.3) is non-differentiable, the classical gradient-
based method is not applicable.

###  4.2. A Fixed-Point Formulation.

We follow [20, 22] to reformulate (4.3) as a fixed-point problem. To this end,
we recall some terminologies of convex analysis [22, 33]. We denote by
$\Gamma_{0}(\mathbb{R}^{n})$ the class of all lower semi-continuous and proper
convex functions mapping $\mathbb{R}^{N}\to\mathbb{R}\cup\\{+\infty\\}$.
Notice that $g\in\Gamma_{0}(\mathbb{R}^{N})$ and
$h\in\Gamma_{0}(\mathbb{R}^{M})$. By
${\left\langle{\cdot},{\cdot}\right\rangle}$ we denote the Euclidean inner
product on $\mathbb{R}^{N}$. For any function
$f\in\Gamma_{0}(\mathbb{R}^{N})$, we define the convex conjugate of $f$ by

$f^{\ast}{\left(u\right)}\coloneqq\sup_{x\in\mathbb{R}^{N}}{\left\\{{\left\langle{x},{u}\right\rangle}-f{\left(u\right)}\right\\}}.$

Let $\mathbb{S}_{+}^{N}\subset\mathbb{R}^{N\times N}$ denote the collection of
all symmetric positive definite matrices. The proximity operator of
$f\in\Gamma_{0}(\mathbb{R}^{N})$ with respect to $P\in\mathbb{S}_{+}^{N}$ is a
function mapping $\mathbb{R}^{N}\to\mathbb{R}^{N}$, defined by

$\operatorname{prox}_{f,P}{\left(x\right)}\coloneqq\operatorname*{arg\,min}_{u\in\mathbb{R}^{N}}{\left\\{f{\left(u\right)}+\frac{1}{2}{\left\|u-x\right\|}_{P}^{2}\right\\}},\qquad\text{for
all }x\in\mathbb{R}^{N},$ (4.4)

where
$\|{\cdot}\|_{P}^{2}\coloneqq{\left\langle{\cdot},{P\cdot}\right\rangle}$. It
is worth remarking that for $f\in\Gamma_{0}(\mathbb{R}^{N})$, we always have
$f^{\ast}\in\Gamma_{0}(\mathbb{R}^{N})$, and there holds

$\operatorname{prox}_{f^{\ast},P}=\operatorname{{\it
I}}-P^{-1}\circ\operatorname{prox}_{f,P^{-1}}\circ P.$ (4.5)

By a direct application of relation between proximity operator and
subdifferntial of functions in $\Gamma_{0}(\mathbb{R}^{N})$, we can easily
derive the following characterization theorem [20, 22].

######  Theorem 4.1.

If $x_{\ast}\in\mathbb{R}^{N}$ is a solution of (4.3), then there exists
$y_{\ast}\in\mathbb{R}^{M}$ such that for all $\alpha>0$,
$P\in\mathbb{S}^{N}_{+}$ and $Q\in\mathbb{S}^{M}_{+}$, there holds

$\left\\{\begin{aligned} x_{\ast}&=\operatorname{prox}_{\alpha
g,P}{\left(x_{\ast}-P^{-1}L^{\ast}y_{\ast}\right)},\\\
y_{\ast}&=\operatorname{prox}_{(\alpha
h)^{\ast},Q}{\left(Q^{-1}Lx_{\ast}+y_{\ast}\right)}.\end{aligned}\right.$
(4.6)

Conversely, if there exist $P\in\mathbb{S}^{N}_{+}$ and
$Q\in\mathbb{S}^{M}_{+}$ satisfying the fixed-point equations (4.6), then the
corresponding $x_{\ast}$ is a solution of (4.3).

It is worth noticing that due to the smoothness of function $g$, there is an
alternative form of fixed-point characterization of solutions of (4.3)
involving gradient of $g$. This gradient approach simplifies the fixed-point
equation 4.6, but at the same time requires exact estimation of Lipschitz
constant of the gradient for convergence analysis, which is complicated for
problem 2.3. To have a clear convergence condition, we choose 4.6 for further
discuss. For details of the gradient approach, readers are referred to [21].

We can rewrite fixed-point equation (4.6) in a compact form. For $\alpha>0$,
$P\in\mathbb{S}_{+}^{N}$ and $Q\in\mathbb{R}^{M}_{+}$, we define block
matrices

$R\coloneqq\begin{bmatrix}P&0\\\ 0&Q\end{bmatrix},\quad
E\coloneqq\begin{bmatrix}I&-P^{-1}L^{\ast}\\\ Q^{-1}L&I\end{bmatrix},$

and function

$\Phi{\left(z\right)}\coloneqq\alpha g{\left(x\right)}+{\left(\alpha
h{\left(y\right)}\right)}^{\ast}.$

Then, (4.6) is equivalent to

$z_{\ast}=\operatorname{prox}_{\Phi,R}{\left(Ez_{\ast}\right)},$ (4.7)

where $z_{\ast}\coloneqq{\left(x_{\ast},y_{\ast}\right)}\in\mathbb{R}^{N+M}$.
That is, $z_{\ast}$ is a fixed-point of operator
$\operatorname{prox}_{\Phi,R}\circ E$.

The fixed-point characterization (4.7) brings two advantages. First of all,
this formulation provides a guidance for designing convergence-guaranteed
algorithms. Moreover, by this characterization, the non-differentiability of
$\|\cdot\|_{1}$ is no longer an obstacle, as long as the proximity operator
$\operatorname{prox}_{\Phi,R}$ can be easily calculated. Indeed, as we will
show next, operator $\operatorname{prox}_{\Phi,R}$ has a closed-form.

######  Proposition 4.2.

If $\alpha>0$, $P\in\mathbb{S}_{+}^{N}$ and $Q:={\rm
diag}(q_{i})_{i=1}^{M}\in\mathbb{S}_{+}^{M}$, then for any
$x\in\mathbb{R}^{N}$ and $y\in\mathbb{R}^{M}$,

$\operatorname{prox}_{\alpha g,P}{\left(x\right)}={\left(P+2\alpha
B\right)}^{-1}{\left(Px-\alpha b\right)},$

and

$\operatorname{prox}_{(\alpha
h)^{\ast},Q}{\left(y\right)}={\left[\operatorname{Proj}_{{\left[-\alpha,\alpha\right]}}{\left(y_{i}-d_{i}/q_{i}\right)}\right]}_{i=1}^{M},$

where $\operatorname{Proj}_{{\left[-\alpha,\alpha\right]}}$ denotes the
projection operator onto ${\left[-\alpha,\alpha\right]}$.

###### Proof.

Both of these closed-form formulas can be shown by the definition of the
proximity operator. By (4.4),

$\operatorname{prox}_{\alpha
g,P}{\left(x\right)}=\operatorname*{arg\,min}_{u\in\mathbb{R}^{N}}{\left\\{\alpha
u^{\ast}Bu+\alpha
b^{\ast}u+\frac{1}{2}{\left\|u-x\right\|}^{2}_{P}\right\\}},$

whose right-hand side is a minimization of a quadratic function of $u$. Taking
the gradient of the right-hand side and solving for the zero gradient, we
obtain the first formula.

It remains to establish the formula for $\operatorname{prox}_{(\alpha
h)^{\ast},Q}$. By identity (4.5), for any $y\in\mathbb{R}^{M}$ we have that

$\operatorname{prox}_{(\alpha
h)^{\ast},Q}{\left(y\right)}=y-Q^{-1}\operatorname{prox}_{\alpha
h,Q^{-1}}{\left(Qy\right)}.$

By setting $\tilde{v}\coloneqq v-d$, we observe that

$\displaystyle\operatorname{prox}_{\alpha h,Q^{-1}}{\left(Qy\right)}$
$\displaystyle=\operatorname*{arg\,min}_{v\in\mathbb{R}^{M}}{\left\\{\alpha{\left\|v-d\right\|}_{1}+\frac{1}{2}{\left\|v-Qy\right\|}^{2}_{Q^{-1}}\right\\}}$
$\displaystyle=d+\operatorname*{arg\,min}_{\tilde{v}\in\mathbb{R}^{M}}{\left\\{\alpha{\left\|\tilde{v}\right\|}_{1}+\frac{1}{2}{\left\|\tilde{v}-{\left(Qy-d\right)}\right\|}^{2}_{Q^{-1}}\right\\}}$
$\displaystyle=d+\operatorname{prox}_{\alpha{\left\|\cdot\right\|}_{1},Q^{-1}}{\left(Qy-d\right)}.$

Since $Q$ is diagonal, we obtain that

$\operatorname{prox}_{\alpha{\left\|\cdot\right\|}_{1},Q^{-1}}{\left(Qy-d\right)}={\left[{\left({\left(y_{i}-\alpha\right)}q_{i}-d_{i}\right)}_{+}\operatorname{Sgn}{\left({\left(y_{i}-\alpha\right)}q_{i}-d_{i}\right)}\right]}_{i=1}^{M},$

where $(\cdot)_{+}\coloneqq\max\\{\cdot,0\\}$ and $\operatorname{Sgn}$ is the
sign function (for details, see [22]). Finally, combining the equations above
yields the desired formula. ∎

With help of 4.2, we can easily obtain the closed form of
$\operatorname{prox}_{\Phi,R}$ by the following calculation. For any
$z={\left(x,y\right)}\in\mathbb{R}^{N+M}$, by definition of proximity operator
4.4, we have

$\displaystyle\operatorname{prox}_{\Phi,R}{\left(z\right)}$
$\displaystyle=\operatorname*{arg\,min}_{u\in\mathbb{R}^{N},v\in\mathbb{R}^{M}}{\left\\{\alpha
g{\left(u\right)}+{\left(\alpha
h\right)}^{\ast}{\left(v\right)}+\frac{1}{2}{\left\|u-x\right\|}^{2}_{P}+\frac{1}{2}{\left\|v-y\right\|}^{2}_{Q}\right\\}}$
$\displaystyle={\left(\operatorname{prox}_{\alpha
g,P}{\left(x\right)},\operatorname{prox}_{{\left(\alpha
h\right)}^{\ast},Q}{\left(y\right)}\right)}.$

###  4.3. Fixed-Point Proximity Algorithms

In this section we consider the numerical algorithm of solving fixed-point
equation (4.7). The naïve Picard’s iteration might not converge, due to the
expanding nature of the linear mapping $E$ [20, 31, 32]. Fortunately a matrix-
splitting technique will do the trick. Specifically, we split $E$ into

$E={\left(E-M\right)}+M,$

then the Fixed-Point Proximity Algorithm (FP2A) is proposed by

$z_{k+1}=\operatorname{prox}_{\Phi,R}{\left({\left(E-M\right)}z_{k+1}+Mz_{k}\right)}.$
(4.8)

Based on properties of proximity operators, a corollary of [31, 32] provides
the following convergence theorem. Define $\mathbb{S}^{N}$ as the collection
of all positively semi-definite matrices on $\mathbb{R}^{N\times N}$.

######  Theorem 4.3.

Suppose $\alpha>0$. If $RM\in\mathbb{S}^{N+M}$, then for any initial
$z_{0}\in\mathbb{R}^{N+M}$, the FP2A scheme (4.8) converges to a solution of
fixed-point equation (4.7).

Notice that Theorem 4.3 extends results from [20], where
$RM\in\mathbb{S}^{N+M}_{+}$ is required. Although (4.8) is generally an
implicit algorithm, if $E-M$ is chosen to be strictly block-triangular, then
(4.8) is actually explicit. We set

$M\coloneqq\begin{bmatrix}I&P^{-1}L^{\ast}\\\ Q^{-1}L&I\end{bmatrix}$

in (4.8), then the explicit FP2A for (4.2) follows as Algorithm 1.

Algorithm 1 Explicit FP2A for (4.2)

0: $\alpha>0$, $P\in\mathbb{S}_{+}^{N}$,
$\\{q_{i}\\}_{i=1}^{M}\subset(0,+\infty)$, $x_{0}\in\mathbb{R}^{N}$ and
$y_{0}\in\mathbb{R}^{M}$.

1: $k\leftarrow 0$

2: repeat

3:
$y_{k+1}\leftarrow[\operatorname{Proj}_{{\left[-\alpha,\alpha\right]}}{\left((y_{k})_{i}+((Lx_{k})_{i}-d_{i})/q_{i}\right)}]_{i=1}^{M}$

4: $x_{k+1}\leftarrow{\left(P+2\alpha
B\right)}^{-1}{\left(L^{\ast}{\left(y_{k}-2y_{k+1}\right)}+Px_{k}-\alpha
b\right)}$

5: $k\leftarrow k+1$

6: until stop criteria is met

6: $x_{\ast}\leftarrow x_{k}$

###  4.4. Convergence Analysis.

This section provides convergence theorem for FP2A algorithm 1 based on
Theorem 4.3, and discusses one-dimension case for illustration propose. A
direct application of Theorem 4.3 gives the following convergence result of
Algorithm 1.

######  Theorem 4.4.

Suppose $\alpha>0$, $P\in\mathbb{S}_{+}^{N}$ and $Q\in\mathbb{S}_{+}^{M}$ is
diagonal. If ${\left\|Q^{-1/2}LP^{-1/2}\right\|}\leq 1$, then for any
$x_{0}\in\mathbb{R}^{N}$ and $y_{0}\in\mathbb{R}^{M}$, explicit FP2A
Algorithm 1 generates a sequence ${\left\\{x_{k}\right\\}}_{k\in\mathbb{N}}$
converging to a solution of (4.2).

###### Proof.

By setting $R\coloneqq{\rm Diag}(P,Q)$ and

$M\coloneqq\begin{bmatrix}I&P^{-1}L^{\ast}\\\ Q^{-1}L&I\end{bmatrix},$

in (4.8), we have FP2A as an explicit iteration as

$\left\\{\begin{aligned} y_{k+1}&=\operatorname{prox}_{(\alpha
h)^{\ast},Q}{\left(y_{k}+Q^{-1}Lx_{k}\right)},\\\
x_{k+1}&=\operatorname{prox}_{\alpha
g,P}{\left(P^{-1}L^{\ast}{\left(y_{k}-2y_{k+1}\right)}+x_{k}\right)}.\end{aligned}\right.$

The above iteration exactly fits into Algorithm 1 by applying 4.2. It is
direct to check that $RM\in\mathbb{S}$ if and only if
$\|Q^{-1/2}LP^{-1/2}\|\leq 1$. Therefore conditions of Theorem 4.3 obtain, and
the convergence directly follows. ∎

Here we discuss the simple case of Theorem 4.4 with $d=1$ and $g=0$ to
illustrate Theorem 4.4 in greater detail. In this case the second-order
elliptic equation reduces to the following homogeneous boundary value problem

$\displaystyle A:u^{\prime\prime}$
$\displaystyle=f,\quad\text{in}\quad{\left[0,1\right]},$ (4.9) $\displaystyle
u$ $\displaystyle=0,\quad\text{on}\quad{\left\\{0,1\right\\}}.$ (4.10)

The ${\mathcal{T}}_{h}$ consists of intervals
$\\{I_{i}\coloneqq[ih,(i+1)h]\\}_{i=0}^{1/h-1}$, and interior faces
${\mathcal{E}}_{h}^{I}={\left\\{ih\right\\}}_{i=1}^{1/h-1}$ and boundary faces
${\mathcal{E}}_{h}^{B}={\left\\{0,1\right\\}}$. Therefore $W_{h}$ and
${\boldsymbol{V}}_{h}$ are linear spaces of piecewise quadratic and linear
functions on ${\mathcal{T}}_{h}$ respectively.

If we choose Lagrangian bases for $W_{h}$ and ${\boldsymbol{V}}_{h}$, then
each element $(w,{\boldsymbol{v}})\in W_{h}\times{\boldsymbol{V}}_{h}$ is
represented by vector
$(x_{w},x_{v})\in\mathbb{R}^{3/h}\times\mathbb{R}^{2/h}$, where
$((x_{w})_{3i},(x_{w})_{3i+1},(x_{w})_{3i+2})$ are function values of $w$ at
points $(ih,(i+1/2)h,(i+1)h)$ on element $I_{i}$, and
$((x_{v})_{2i},(x_{v})_{2i+1})$ are function values of ${\boldsymbol{v}}$ at
points $(ih,(i+1)h)$ on element $I_{i}$. Then by (4.1), we have
$L_{1}=[0_{(1/h-1)\times(3/h)},L_{12}]$,
$L_{2}=[L_{21},0_{(1/h-1)\times(2/h)}]$ and $L_{3}=[L_{31},0_{2\times(2/h)}]$,
where

$\displaystyle L_{12}$
$\displaystyle=\begin{bmatrix}0&1&-1&0&0&\cdots&0&0&0\\\
0&0&0&1&-1&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&0&\cdots&1&-1&0\end{bmatrix}\in\mathbb{R}^{(1/h-1)\times(2/h)},$
(4.11) $\displaystyle L_{21}$
$\displaystyle=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}0&0&1&-1&0&0&0&\cdots&0&0&0&0\\\
0&0&0&0&0&1&-1&\cdots&0&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&0&0&0&\cdots&1&-1&0&0\end{bmatrix}\in\mathbb{R}^{(1/h-1)\times(3/h)},$
(4.12) $\displaystyle L_{31}$
$\displaystyle=\begin{bmatrix}1&0&0&\cdots&0&0&0\\\
0&0&0&\cdots&0&0&1\end{bmatrix}\in\mathbb{R}^{2\times(3/h)}.$ (4.13)

Therefore we have the following corollary of Theorem 4.4 for problem 4.9 and
4.10.

######  Proposition 4.5.

Suppose $\alpha,\lambda>0$, $\\{q_{i}\\}_{i=1}^{2/h}\subset(0,+\infty)$ and
$P=\lambda I_{5/h}$, $Q={\rm Diag}(q_{i})_{i=1}^{2/h}$, $L$ is defined as in
4.11, 4.12 and 4.13. If

$\lambda\geq\frac{\tau^{2}}{h^{4}}\max{\left\\{2h^{2}\max{\left\\{q_{i}^{-1}\right\\}}_{i=1}^{2/h-2},q_{2/h-1}^{-1},q_{2/h}^{-1}\right\\}},$
(4.14)

then for any $x_{0}\in\mathbb{R}^{5/h}$ and $y_{0}\in\mathbb{R}^{2/h}$,
explicit FP2A Algorithm 1 generates a sequence
${\left\\{x_{k}\right\\}}_{k\in\mathbb{N}}$ converging to a solution of (4.2).

###### Proof.

By Theorem 4.4, it is sufficient to prove $\|Q^{1/2}LP^{-1/2}\|<1$. Observe
that

$LL^{\ast}=\frac{\tau^{2}}{h^{2}}{\rm
Diag}{\left(L_{11}L_{11}^{\ast},L_{22}L_{22}^{\ast},h^{-2}L_{32}L_{32}^{\ast}\right)}=\frac{\tau^{2}}{h^{4}}{\rm
Diag}{\left(2h^{2}I_{2/h-2},I_{2}\right)}.$

Then by condition, we have

$\displaystyle{\left\|Q^{-1/2}LP^{-1/2}\right\|}^{2}$
$\displaystyle={\left\|Q^{-1/2}LP^{-1}L^{\ast}Q^{-1/2}\right\|}$
$\displaystyle=\frac{\tau^{2}}{\lambda
h^{4}}\max{\left\\{2h^{2}\max{\left\\{q_{i}^{-1}\right\\}}_{i=1}^{2/h-2},q_{2/h-1}^{-1},q_{2/h}^{-1}\right\\}}\leq
1.$

Finally, by Theorem 4.4, we finish the proof. ∎

Proposition 4.5 extends results from [20] as previously discussed, where the
strict inequality is required in (4.14). Proposition 4.5 also implies that for
all given $h$, $P=\lambda I$ and diagonal $Q\in\mathbb{S}_{+}$, FP2A
Algorithm 1 always converges with sufficiently large $\lambda>0$. Similar
result can be easily obtained for general $P\in\mathbb{S}_{+}$ with
sufficiently large $\|P\|$. For case $d=2,3$, the matrix $L$ depends on the
adjacency matrix. Nevertheless there are always similar structures as 4.11,
4.12 and 4.13 as long as Lagrangian bases are picked for $W_{n}$ and
${\boldsymbol{V}}_{h}$.

##  5\. Numerical Results

In this section we provide numerical experiments for 1.1 validating the error
estimates obtained in this work via explicit FP2A method. Specially, we test
$k=2$ in numerical experiments, with regular solution on square domain, or
singular solution on L-shaped domain. In each case we test different type of
coefficients $A$, i.e. constant, continuous and discontinuous coefficients. In
Sections 5.1 and 5.2 we choose the Lagrangian basis for $W_{h}$ and
${\boldsymbol{V}}_{h}$ and show the order of various errors as $N=1/h$
doubles. In Section 5.3 we show examples of solving 1.1 with multiscale DG
bases. We demonstrate the numerical solution naturally possess a sparse
representation with such choice of basis.

###  5.1. Numerical Results for Various Coefficients on Square Domain

We first consider problem (2.1) on domain $[0,1]\times[0,1]$ with constant,
continuous and discontinuous coefficient matrices $A$, respectively. Whatever
coefficient matrix $A$ is, the right-hand side functions $f$ and $g$ in (2.1)
are chosen such that the exact solution is always

$u{\left(x,y\right)}=\sin{\left(\pi x\right)}\sin{\left(\pi y\right)}.$

####  5.1.1. Constant Coefficients

Here we choose

$A=\begin{bmatrix}2&1\\\ 1&2\end{bmatrix}$

in (2.1). Table 1 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{0,2}$ | order
---|---|---|---|---|---|---|---|---
4 | 4.92e$-$02 | – | 2.58e$-$01 | – | 3.52e+00 | – | 2.59e$-$01 | –
8 | 1.42e$-$02 | 1.80 | 7.49e$-$02 | 1.78 | 1.76e+00 | 1.00 | 7.59e$-$02 | 1.77
16 | 3.80e$-$03 | 1.90 | 2.01e$-$02 | 1.90 | 8.61e$-$01 | 1.04 | 2.02e$-$02 | 1.91
32 | 9.75e$-$04 | 1.96 | 5.16e$-$03 | 1.96 | 4.23e$-$01 | 1.02 | 5.14e$-$03 | 1.98
64 | 2.46e$-$04 | 1.99 | 1.30e$-$03 | 1.99 | 2.10e$-$01 | 1.01 | 1.29e$-$03 | 2.00
Table 1. Various errors and convergence rates for regular solution with
constant coefficients.

####  5.1.2. Continuous Coefficients

Here we choose

$A=\begin{bmatrix}1+x&\sqrt{xy}\\\ \sqrt{xy}&1+y\end{bmatrix}$

in (2.1). Table 2 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{2}$ | order
---|---|---|---|---|---|---|---|---
4 | 4.48e$-$02 | – | 2.44e$-$01 | – | 3.55e+00 | – | 2.51e$-$01 | –
8 | 1.21e$-$02 | 1.89 | 6.59e$-$02 | 1.89 | 1.75e+00 | 1.02 | 6.68e$-$02 | 1.91
16 | 3.18e$-$03 | 1.93 | 1.72e$-$02 | 1.94 | 8.56e$-$01 | 1.03 | 1.72e$-$02 | 1.96
32 | 8.08e$-$04 | 1.98 | 4.36e$-$03 | 1.98 | 4.23e$-$01 | 1.02 | 4.30e$-$03 | 2.00
64 | 2.02e$-$04 | 2.00 | 1.09e$-$03 | 2.00 | 2.10e$-$01 | 1.01 | 1.07e$-$03 | 2.01
Table 2. Various errors and convergence rates for regular solution with
continuous coefficients.

####  5.1.3. Discontinuous Coefficients

Here we choose

$A=\begin{bmatrix}2&{\rm
sgn}{\left({\left(x-0.5\right)}{\left(y-0.5\right)}\right)}\\\ {\rm
sgn}{\left({\left(x-0.5\right)}{\left(y-0.5\right)}\right)}&2\end{bmatrix}$

in (2.1). Table 3 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{2}$ | order
---|---|---|---|---|---|---|---|---
4 | 4.89e$-$02 | – | 2.56e$-$01 | – | 3.50e+00 | – | 2.52e$-$01 | –
8 | 1.33e$-$02 | 1.88 | 7.04e$-$02 | 1.86 | 1.75e+00 | 1.00 | 7.02e$-$02 | 1.84
16 | 3.58e$-$03 | 1.89 | 1.88e$-$02 | 1.90 | 8.57e$-$01 | 1.03 | 1.88e$-$02 | 1.90
32 | 9.38e$-$04 | 1.93 | 4.88e$-$03 | 1.95 | 4.23e$-$01 | 1.02 | 4.84e$-$03 | 1.96
64 | 2.40e$-$04 | 1.96 | 1.24e$-$03 | 1.97 | 2.10e$-$01 | 1.01 | 1.23e$-$03 | 1.98
Table 3. Various errors and convergence rates for regular solution with
discontinuous coefficients.

###  5.2. Numerical Results for Various Coefficients on L-Shaped Domain

We then consider singular solution of (2.1) on L-shaped domain with constant,
continuous and discontinuous coefficient matrices $A$, respectively. Whatever
coefficient matrix $A$ is, the right-hand side functions $f$ and $g$ in (2.1)
are chosen such that the exact solution is always

$u{\left(r,\theta\right)}=r^{2/3}\sin{\left(2\theta/3\right)}$

over the L-shaped domain
${\left[-1,1\right]}\times{\left[-1,1\right]}\setminus(0,1]\times[-1,0)$.

####  5.2.1. Constant Coefficients

Here we choose

$A=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}$

in (2.1). Table 4 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{0,2}$ | order
---|---|---|---|---|---|---|---|---
4 | 9.32e$-$03 | – | 7.93e$-$02 | – | 1.08e+00 | – | 7.99e$-$02 | –
8 | 6.86e$-$03 | 0.44 | 3.45e$-$02 | 1.20 | 8.67e$-$01 | 0.32 | 3.64e$-$02 | 1.13
16 | 5.29e$-$03 | 0.38 | 1.99e$-$02 | 0.79 | 7.08e$-$01 | 0.29 | 2.08e$-$02 | 0.81
32 | 3.46e$-$03 | 0.61 | 1.21e$-$02 | 0.72 | 5.89e$-$01 | 0.27 | 1.23e$-$02 | 0.75
64 | 2.34e$-$03 | 0.57 | 7.80e$-$03 | 0.63 | 4.85e$-$01 | 0.28 | 7.88e$-$03 | 0.65
Table 4. Various errors and convergence rates for singular solution with
constant coefficients.

####  5.2.2. Continuous Coefficients

Here we choose

$A=\begin{bmatrix}1+{\left|x\right|}&\sqrt{{\left|xy\right|}}\\\
\sqrt{{\left|xy\right|}}&1+{\left|y\right|}\end{bmatrix}$

in (2.1). Table 5 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{2}$ | order
---|---|---|---|---|---|---|---|---
4 | 1.06e$-$01 | – | 8.47e$-$01 | – | 5.02e+00 | – | 8.50e$-$01 | –
8 | 7.78e$-$03 | 3.76 | 3.73e$-$02 | 4.51 | 8.69e$-$01 | 2.53 | 3.89e$-$02 | 4.45
16 | 6.31e$-$03 | 0.30 | 2.28e$-$02 | 0.71 | 7.08e$-$01 | 0.30 | 2.36e$-$02 | 0.72
32 | 4.28e$-$03 | 0.56 | 1.44e$-$02 | 0.66 | 5.90e$-$01 | 0.26 | 1.47e$-$02 | 0.68
64 | 3.04e$-$03 | 0.50 | 9.85e$-$03 | 0.55 | 4.88e$-$01 | 0.27 | 9.93e$-$03 | 0.57
Table 5. Various errors and convergence rates for singular solution with
continuous coefficients.

####  5.2.3. Discontinuous Coefficients

Here we choose

$A=\begin{bmatrix}2&{\rm sgn}{\left(xy\right)}\\\ {\rm
sgn}{\left(xy\right)}&2\end{bmatrix}$

in (2.1). Table 6 shows various error measures and convergence rates.

$N$ | ${\left\|u-u_{h}\right\|}_{2}$ | order | ${\left|u-u_{h}\right|}_{1,2}$ | order | ${\left|u-u_{h}\right|}_{2,2}$ | order | ${\left\|{\boldsymbol{q}}-{\boldsymbol{q}}_{h}\right\|}_{2}$ | order
---|---|---|---|---|---|---|---|---
4 | 9.96e$-$03 | – | 7.94e$-$02 | – | 1.08e+00 | – | 8.01e$-$02 | –
8 | 4.76e$-$03 | 1.06 | 3.16e$-$02 | 1.33 | 8.67e$-$01 | 0.32 | 3.39e$-$02 | 1.24
16 | 4.81e$-$03 | -0.01 | 1.87e$-$02 | 0.76 | 7.11e$-$01 | 0.29 | 1.99e$-$02 | 0.77
32 | 4.17e$-$03 | 0.21 | 1.41e$-$02 | 0.41 | 5.96e$-$01 | 0.26 | 1.44e$-$02 | 0.47
64 | 3.60e$-$03 | 0.21 | 1.16e$-$02 | 0.27 | 4.93e$-$01 | 0.27 | 1.17e$-$02 | 0.30
Table 6. Various errors and convergence rates for singular solution with
discontinuous coefficients.

###  5.3. Numerical Results Using a Multiscale Basis

In this section we perform examples of solving 1.1 using multiscale bases to
illustrate the adaptivity/sparsity feature of $L^{1}$ stabilization. We test
4.9 and 4.10 with following example

$u(x)=\begin{cases}t{\left({\left(\frac{x}{t}\right)}^{n+1}-\frac{x}{t}\right)},&0\leq
x\leq t,\\\
{\left(1-t\right)}{\left(\frac{1-x}{1-t}-{\left(\frac{1-x}{1-t}\right)}^{n+1}\right)},&t<x\leq
1,\\\ \end{cases}$ (5.1)

where $0<t<1$ and $n>0$. The corresponding $q$ has a kink at $x=t$:

$q(x)=u^{\prime}(x)=\begin{cases}{\left(n+1\right)}{\left(\frac{x}{t}\right)}^{n}-1,&0\leq
x\leq t,\\\ {\left(n+1\right)}{\left(\frac{1-x}{1-t}\right)}^{n}-1,&t<x\leq
1.\\\ \end{cases}$

Now we introduce the multiscale piecewise polynomial space and basis on
$[0,1]$ once presented in [7]. Define $\mathbb{X}_{n}\coloneqq W_{2^{-n}}$ for
$n\in\mathbb{N}$ and transformations of $f\in L^{2}{\left[0,1\right]}$ as

${\left(\tau_{0}f\right)}{\left(x\right)}\coloneqq
f{\left(2x\right)}\chi_{{\left[0,1\right]}}{\left(2x\right)},\qquad{\left(\tau_{1}f\right)}{\left(x\right)}\coloneqq
f{\left(2x-1\right)}\chi_{{\left[0,1\right]}}{\left(2x-1\right)},$

where $\chi_{{\left[0,1\right]}}:\mathbb{R}\to\mathbb{R}$ takes $1$ on
${\left[0,1\right]}$ and $0$ otherwise. Then it is proven in [7] that for all
$n\in\mathbb{N}$,

$\mathbb{X}_{n}=\mathbb{X}_{0}\oplus^{\perp}\mathbb{W}_{1}\oplus^{\perp}\cdots\oplus^{\perp}\mathbb{W}_{n},$

where $\mathbb{W}_{1}$ is defined by
$\mathbb{X}_{2}=\mathbb{X}_{1}\oplus^{\perp}\mathbb{W}_{1}$ and for $n\geq 2$,

$\mathbb{W}_{n}\coloneqq\tau_{0}\mathbb{W}_{n-1}\oplus\tau_{1}\mathbb{W}_{n-1}.$

Here ‘$\oplus^{\perp}$’ indicates addition between two orthogonal sets in
$L^{2}$.

In [7] the multiscale piecewise quadratic polynomial basis is proposed as
follows. There are three basis functions for $\mathbb{X}_{0}$, as

$w_{00}{\left(x\right)}\coloneqq 1,\quad
w_{01}{\left(x\right)}\coloneqq\sqrt{3}{\left(2x-1\right)},\quad
w_{02}{\left(x\right)}\coloneqq\sqrt{5}{\left(6x^{2}-6x+1\right)},$

and for $\mathbb{W}_{1}$,

$\displaystyle w_{10}{\left(x\right)}$
$\displaystyle\coloneqq\begin{cases}1-6x,&x\in{\left[0,1/2\right]},\\\
5-6x,&x\in(1/2,1],\end{cases}$ $\displaystyle w_{11}{\left(x\right)}$
$\displaystyle\coloneqq\begin{cases}\frac{\sqrt{91}}{31}{\left(240x^{2}-116x+9\right)},&x\in{\left[0,1/2\right]},\\\
\frac{\sqrt{91}}{31}{\left(3-4x\right)},&x\in(1/2,1],\end{cases}$
$\displaystyle w_{12}{\left(x\right)}$
$\displaystyle\coloneqq\begin{cases}\frac{\sqrt{91}}{31}{\left(4x-1\right)},&x\in{\left[0,1/2\right]},\\\
\frac{\sqrt{91}}{31}{\left(240x^{2}-364x+133\right)},&x\in(1/2,1].\end{cases}$

Define $M_{1}\coloneqq{\left\\{w_{1i}\right\\}}_{i=1}^{3}$ and
$M_{n}\coloneqq\tau_{0}M_{n-1}\cup\tau_{1}M_{n-1}$ for all $n\geq 2$. Then [7]
proves that $M_{n}$ is a basis of $\mathbb{W}_{n}$ for $n\geq 1$. Such
multiscale basis has vanishing moment of $4$ since by construction we have
$M_{n}\perp\mathbb{X}_{0}$, and this gives us a sparse representation of
solutions. Finally the elements of representation are truncated under certain
threshold.

For problem (5.1), we set $t=1/6$ and $h=2^{-7}$, and the reconstructed
function and coefficients of single-/multi-scale basis for reconstructed $u$
are shown in Figure 2. Here we truncate all the multi-scale coefficient with
absolute value under $10^{-5}$. One can see that the coefficients of single-
scale basis are dense, while multiscale basis decade rapidly due to its
vanishing moment, but have strokes around kink $x=1/6$. This shows the
adaptivity/sparsity of the proposed multiscale basis comparing to single-scale
basis. Table 7 quantifies the $L^{2}$ error and the coefficient sparsity of
single-/multi-scale solutions with various truncation thresholds.

Figure 1. Ground truth and reconstructed solutions with single- and multi-scale bases. Single- and multi-scale solutions overlap in plot. Figure 2. Coefficients comparison of multi- and single-scale bases. Basis | Truncation Threshold | ${\left\|u-u_{h}\right\|}_{2}$ | Coefficient Sparsity
---|---|---|---
Single-scale | – | $1.097\times 10^{-2}$ | $0\%$
Multi-scale | $10^{-2}$ | $1.181\times 10^{-2}$ | $90.63\%$
$10^{-3}$ | $1.099\times 10^{-2}$ | $85.42\%$
$10^{-4}$ | $1.097\times 10^{-2}$ | $78.13\%$
Table 7. The $L^{2}$-error and the coefficient sparsity of single- and multi-
scale basis solutions with various truncation thresholds. Coefficient sparsity
is defined as the percentage of zero elements.

##  6\. Concluding Remark

In this paper, we propose a novel $DG$ method for the second order elliptic
equations in the non-divergence form with $L^{1}$ stabilization. The
motivation for such an approach is to develop numerical schemes for PDEs such
that the resulting solutions have sparse representation. This work is the
first attempt along this path. Due to the use of $L^{1}$ stabilization in the
minimization formulation, we are unable to use existing analytical tools for
gradient based methods. We believe the analysis used in the error estimates
and the FP2A algorithm can be applied to more challenging problems with low
regularity settings. It is worth to point out that the original $L^{1}$
regularization technique is different from what we proposed in this paper.
Namely, the purpose of regularization is to relax the functional from ill-
posedness. Meanwhile, in our work the $L^{1}$ term is a stablization term
which ensures the convergence and existence of the minimizer of the discrete
problem without breaking the consistency with the PDE. In turn, the sparsity
promoting feature of the multiscale piecewise polynomial basis is not fully
utilized in the current formulation. Nevertheless, we believe that a possible
way to improve the performance of the scheme is to add a regularization term
such as $\|Bx\|_{1}$ in this functional where $B$ is the multiscale piecewise
polynomial decomposition. This is subject to our ongoing work in this project.

## References

  * [1] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, and L. Donatella Marini, _Unified analysis of discontinuous Galerkin methods for elliptic problems_ , SIAM J. Numer. Anal. 39 (2001/02), no. 5, 1749–1779. MR 1885715
  * [2] Susanne C. Brenner and L. Ridgway Scott, _The mathematical theory of finite element methods_ , third ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. MR 2373954
  * [3] Susanne C. Brenner and Li-Yeng Sung, _Virtual enriching operators_ , Calcolo 56 (2019), no. 4, Paper No. 44, 25. MR 4019733
  * [4] Waixiang Cao, Junping Wang, and Yuesheng Xu, _An $l^{p}$-weak galerkin method for second order elliptic equations in non-divergence form_, 2021.
  * [5] Zhongying Chen, Charles A. Micchelli, and Yuesheng Xu, _A construction of interpolating wavelets on invariant sets_ , Math. Comp. 68 (1999), no. 228, 1569–1587. MR 1651746
  * [6] by same author, _Fast collocation methods for second kind integral equations_ , SIAM J. Numer. Anal. 40 (2002), no. 1, 344–375. MR 1921923
  * [7] by same author, _Multiscale methods for Fredholm integral equations_ , Cambridge Monographs on Applied and Computational Mathematics, vol. 28, Cambridge University Press, Cambridge, 2015. MR 3381986
  * [8] Monique Dauge, _Neumann and mixed problems on curvilinear polyhedra_ , Integral Equations Operator Theory 15 (1992), no. 2, 227–261. MR 1147281
  * [9] Andreas Dedner and Tristan Pryer, _Discontinuous Galerkin methods for a class of nonvariational problems_ , Commun. Appl. Math. Comput. 4 (2022), no. 2, 634–656. MR 4412338
  * [10] Jim Douglas, Jr., Todd Dupont, Peter Percell, and Ridgway Scott, _A family of $C^{1}$ finite elements with optimal approximation properties for various Galerkin methods for 2nd and 4th order problems_, RAIRO Anal. Numér. 13 (1979), no. 3, 227–255. MR 543934
  * [11] Xiaobing Feng, Lauren Hennings, and Michael Neilan, _Finite element methods for second order linear elliptic partial differential equations in non-divergence form_ , Math. Comp. 86 (2017), no. 307, 2025–2051. MR 3647950
  * [12] Xiaobing Feng, Michael Neilan, and Stefan Schnake, _Interior penalty discontinuous Galerkin methods for second order linear non-divergence form elliptic PDEs_ , J. Sci. Comput. 74 (2018), no. 3, 1651–1676. MR 3767822
  * [13] Wendell H. Fleming and H. Mete Soner, _Controlled Markov processes and viscosity solutions_ , second ed., Stochastic Modelling and Applied Probability, vol. 25, Springer, New York, 2006. MR 2179357
  * [14] Dietmar Gallistl, _Variational formulation and numerical analysis of linear elliptic equations in nondivergence form with Cordes coefficients_ , SIAM J. Numer. Anal. 55 (2017), no. 2, 737–757. MR 3628316
  * [15] David Gilbarg and Neil S. Trudinger, _Elliptic partial differential equations of second order_ , Classics in Mathematics, Springer-Verlag, Berlin, 2001, Reprint of the 1998 edition. MR 1814364
  * [16] J. L. Guermond, _A finite element technique for solving first-order PDEs in $L^{P}$_, SIAM J. Numer. Anal. 42 (2004), no. 2, 714–737. MR 2084233
  * [17] Max Jensen and Iain Smears, _On the convergence of finite element methods for Hamilton-Jacobi-Bellman equations_ , SIAM J. Numer. Anal. 51 (2013), no. 1, 137–162. MR 3033005
  * [18] Ohannes A. Karakashian and Frederic Pascal, _Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic problems_ , SIAM J. Numer. Anal. 45 (2007), no. 2, 641–665. MR 2300291
  * [19] Omar Lakkis and Tristan Pryer, _A finite element method for second order nonvariational elliptic problems_ , SIAM J. Sci. Comput. 33 (2011), no. 2, 786–801. MR 2801189
  * [20] Qia Li, Lixin Shen, Yuesheng Xu, and Na Zhang, _Multi-step fixed-point proximity algorithms for solving a class of optimization problems arising from image processing_ , Adv. Comput. Math. 41 (2015), no. 2, 387–422. MR 3337498
  * [21] Qia Li and Na Zhang, _Fast proximity-gradient algorithms for structured convex optimization problems_ , Appl. Comput. Harmon. Anal. 41 (2016), no. 2, 491–517. MR 3534448
  * [22] Charles A. Micchelli, Lixin Shen, and Yuesheng Xu, _Proximity algorithms for image models: denoising_ , Inverse Problems 27 (2011), no. 4, 045009, 30. MR 2781033
  * [23] Charles A. Micchelli, Lixin Shen, Yuesheng Xu, and Xueying Zeng, _Proximity algorithms for the L1/TV image denoising model_ , Adv. Comput. Math. 38 (2013), no. 2, 401–426. MR 3019155
  * [24] Charles A. Micchelli and Yuesheng Xu, _Using the matrix refinement equation for the construction of wavelets on invariant sets_ , Appl. Comput. Harmon. Anal. 1 (1994), no. 4, 391–401. MR 1310661
  * [25] by same author, _Reconstruction and decomposition algorithms for biorthogonal multiwavelets_ , Multidimens. Systems Signal Process. 8 (1997), no. 1-2, 31–69. MR 1433530
  * [26] Michael Neilan, _Convergence analysis of a finite element method for second order non-variational elliptic problems_ , J. Numer. Math. 25 (2017), no. 3, 169–184. MR 3707104
  * [27] Michael Neilan, Abner J. Salgado, and Wujun Zhang, _The Monge-Ampére equation_ , Geometric partial differential equations. Part I, Handb. Numer. Anal., vol. 21, Elsevier/North-Holland, Amsterdam, [2020] ©2020, pp. 105–219. MR 4378427
  * [28] Michael Neilan and Mohan Wu, _Discrete Miranda-Talenti estimates and applications to linear and nonlinear PDEs_ , J. Comput. Appl. Math. 356 (2019), 358–376. MR 3921147
  * [29] Ricardo H. Nochetto and Wujun Zhang, _Discrete ABP estimate and convergence rates for linear elliptic equations in non-divergence form_ , Found. Comput. Math. 18 (2018), no. 3, 537–593. MR 3807356
  * [30] Weifeng Qiu and Shun Zhang, _Adaptive first-order system least-squares finite element methods for second-order elliptic equations in nondivergence form_ , SIAM J. Numer. Anal. 58 (2020), no. 6, 3286–3308. MR 4173220
  * [31] Jin Ren, _Inexact fixed-point proximity algorithms for nonsmooth convex optimization_ , PhD dissertation, Old Dominion University, 2022.
  * [32] Jin Ren and Yuesheng Xu, _Inexact fixed-point iterations with applications to convex optimization_ , Preprint.
  * [33] R. Tyrrell Rockafellar and Roger J.-B. Wets, _Variational analysis_ , Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], vol. 317, Springer-Verlag, Berlin, 1998. MR 1491362
  * [34] Iain Smears and Endre Süli, _Discontinuous Galerkin finite element approximation of nondivergence form elliptic equations with Cordès coefficients_ , SIAM J. Numer. Anal. 51 (2013), no. 4, 2088–2106. MR 3077903
  * [35] Chunmei Wang and Junping Wang, _A primal-dual weak Galerkin finite element method for second order elliptic equations in non-divergence form_ , Math. Comp. 87 (2018), no. 310, 515–545. MR 3739209