# On the necessity of the inf-sup condition for a mixed finite element
formulation

Fleurianne Bertrand  University of Twente, The Netherlands  and  Daniele Boffi
King Abdullah University of Science and Technology (KAUST), Saudi Arabia and
University of Pavia, Italy

###### Abstract.

We study a non standard mixed formulation of the Poisson problem, sometimes
known as dual mixed formulation. For reasons related to the equilibration of
the flux, we use finite elements that are conforming in
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ for the approximation of the
gradients, even if the formulation would allow for discontinuous finite
elements. The scheme is not uniformly inf-sup stable, but we can show
existence and uniqueness of the solution, as well as optimal error estimates
for the gradient variable when suitable regularity assumptions are made.
Several additional remarks complete the paper, shedding some light on the
sources of instability for mixed formulations.

###### 2000 Mathematics Subject Classification:

Primary 65N30. Secondary 65N12

##  1\. Introduction

In this paper we discuss the numerical approximation of saddle point problems
of the following form: given two Hilbert spaces $V$ and $Q$, two continuous
bilinear forms $a(\cdot,\cdot):V\times V\to\mathbb{R}$ and
$b(\cdot,\cdot):V\times Q\to\mathbb{R}$, and two functionals $f\in V^{\prime}$
and $g\in Q^{\prime}$, find $u\in V$ and $p\in Q$ such that

(1) $\left\\{\begin{aligned} &a(u,v)+b(v,p)={}_{V^{\prime}}\langle
f,v\rangle_{V}&&\forall v\in V\\\ &b(u,q)={}_{Q^{\prime}}\langle
g,q\rangle_{Q}&&\forall q\in Q.\end{aligned}\right.$

It is well-known that a conforming approximation of the problem relies on
suitable stability conditions, which are usually referred to as _inf-sup
conditions_ (see Section 2). While it is universally understood that the inf-
sup conditions are sufficient for the quasi-optimal convergence of any
Galerkin discretization, the question whether such conditions are also
necessary is less studied; nevertheless it is a common belief that in general
the inf-sup conditions are essentially sufficient and necessary for the
optimal behavior of a numerical scheme and everybody agrees that inf-sup
unstable formulations should be avoided unless special tricks are adopted
(stabilizations, filtering of spurious modes, special meshes, etc.).

In this paper we discuss the necessity of the inf-sup conditions by studying
the approximation of a non standard mixed formulation for the Poisson
equation. In particular, we present a scheme which, under suitable conditions,
is optimally convergent even if the inf-sup constant goes to zero as the mesh
is refined. This (counter-)example can be consider as an extension of a one
dimensional toy problem that has been studied in [1]. The interested reader is
also referred to the abstract setting sketched in [6, Section 5.6.2] where it
is shown what can happen when the inf-sup condition goes wrong.

After an introductory section about the inf-sup conditions, in Section 3 we
present the dual mixed formulation of the Poisson equation. Section 4 deals
with discrete inf-sup constant, including some numerical tests showing the
mesh dependent behavior of the stability condition. In Section 5 we show how
the finite element spaces can be spit into a stable part and an unstable one.
Sections 6 and 7 present the convergence theoretical results and some
numerical tests confirming the theory. Finally, two appendices conclude the
paper with some links between the considered problem and a flux equilibration
strategy.

##  2\. Generalities about the inf-sup conditions

In this section and in the sequel of this paper we follow the framework of
[6]; we recall some relevant and well known results for completeness and for
setting our notation.

The conforming Galerkin approximation of the mixed formulation presented in
(1) consists in choosing appropriate finite element subspaces $V_{h}\subset V$
and $Q_{h}\subset Q$, and in finding $u_{h}\in V_{h}$ and $p_{h}\in Q_{h}$
such that

(2) $\left\\{\begin{aligned} &a(u_{h},v)+b(v,p_{h})={}_{V^{\prime}}\langle
f,v\rangle_{V}&&\forall v\in V_{h}\\\ &b(u_{h},q)={}_{Q^{\prime}}\langle
g,q\rangle_{Q}&&\forall q\in Q_{h}.\end{aligned}\right.$

Let $N_{V}$ be the dimension of $V_{h}$ and $N_{Q}$ the one of $Q_{h}$, the
matrix form of the discrete problem expressed in (2) is given by

(3) $\left(\begin{matrix}\mathsf{A}&\mathsf{B}^{\top}\\\
\mathsf{B}&\mathsf{0}\end{matrix}\right)\left(\begin{matrix}\mathsf{u}\\\
\mathsf{p}\end{matrix}\right)=\left(\begin{matrix}\mathsf{f}\\\
\mathsf{g}\end{matrix}\right),$

where $\mathsf{A}$ is a square matrix of size $N_{V}\times N_{V}$,
$\mathsf{B}$ is a rectangular matrix of size $N_{Q}\times N_{V}$,
$\mathsf{u}\in\mathbb{R}^{N_{V}}$ and $\mathsf{p}\in\mathbb{R}^{N_{Q}}$ are
column vector representations of $u\in V_{h}$ and $p\in Q_{h}$, respectively,
and $\mathsf{f}\in\mathbb{R}^{N_{V}}$ and $\mathsf{g}\in\mathbb{R}^{N_{Q}}$
are column vector realizations of the right hand sides $f$ and $g$,
respectively.

The necessary and sufficient conditions for the solvability of (3) are
summarized in [6, Theorem 3.2.1], which we now recall for the reader’s
convenience. We denote by $\mathsf{K}$ the kernel of $\mathsf{B}$

$\mathsf{K}=\ker\mathsf{B}.$

The restriction of $\mathsf{A}$ to $\mathsf{K}$ is denoted by
$\mathsf{A}_{\mathsf{K}\mathsf{K}}$

$\mathsf{A}_{\mathsf{K}\mathsf{K}}=\Pi_{\mathsf{K}}\mathsf{A}E_{\mathsf{K}},$

where $\Pi_{\mathsf{K}}$ and $E_{\mathsf{K}}$ are the projection from
$\mathbb{R}^{N_{V}}$ onto $\mathsf{K}$ and the embedding from $\mathsf{K}$
into $\mathbb{R}^{N_{V}}$, respectively. The two conditions equivalent to the
solvability of (3) (for all possible right hand sides) are the following ones.

M1:

The matrix $\mathsf{A}_{\mathsf{K}\mathsf{K}}$ is invertible.

M2:

$N_{V}\geq N_{Q}$ and the matrix $\mathsf{B}$ is full rank.

M1 can be expressed by saying that the operator associated with
$\mathsf{A}_{\mathsf{K}\mathsf{K}}$ is surjective or, equivalently, injective;
analogously, M2 states that the operator associated with $\mathsf{B}$ is
surjective or, equivalently, that the one associated with $\mathsf{B}^{\top}$
is injective.

The essential ideas behind the inf-sup theory is that a uniform stability of
problem (2) with respect to the parameter $h$ requires that conditions M1 and
M2 are made explicit and uniform with respect to $h$. This is done by
introducing suitable inf-sup conditions.

Let $\mathsf{K}_{\mathsf{B}_{h}}$ be the subspace of $V_{h}$ associated with
the kernel of the matrix $\mathsf{B}$

$\mathsf{K}_{\mathsf{B}_{h}}=\\{v_{h}\in V_{h}:b(v_{h},q)=0\ \forall q\in
Q_{h}\\}.$

Then hypotheses M1 and M2 correspond to the following two inf-sup conditions,
respectively (where, as usual, when there is an inf-sup involving fractions,
we understand that the infimum and the supremum are taken over non vanishing
functions).

IS1:

There exists a constant $\alpha_{h}>0$ such that

(4)
$\inf_{v_{h}\in\mathsf{K}_{\mathsf{B}_{h}}}\sup_{w_{h}\in\mathsf{K}_{\mathsf{B}_{h}}}\frac{a(v_{h},w_{h})}{\|v_{h}\|_{V}\|w_{h}\|_{V}}\geq\alpha_{h}.$

IS2:

There exists a constant $\beta_{h}>0$ such that

(5) $\inf_{q_{h}\in Q_{h}}\sup_{v_{h}\in
V_{h}}\frac{b(v_{h},q_{h})}{\|v_{h}\|_{V}\|q_{h}\|_{Q}}\geq\beta_{h}.$

We assume that the continuous problem (1) is stable; in particular, the
following conditions analogue to IS1 and IS2 guarantee this property, where
the continuous kernel is defined as

$\mathsf{K}_{\mathsf{B}}=\\{v\in V:b(v,q)=0\ \forall q\in Q\\}.$

  * •

There exists a constant $\alpha>0$ such that

(6)
$\inf_{v\in\mathsf{K}_{\mathsf{B}}}\sup_{w\in\mathsf{K}_{\mathsf{B}}}\frac{a(v,w)}{\|v\|_{V}\|w\|_{V}}\geq\alpha$

and

(7)
$\inf_{w\in\mathsf{K}_{\mathsf{B}}}\sup_{v\in\mathsf{K}_{\mathsf{B}}}\frac{a(v,w)}{\|v\|_{V}\|w\|_{V}}\geq\alpha.$

  * •

There exists a constant $\beta>0$ such that

(8) $\inf_{q\in Q}\sup_{v\in V}\frac{b(v,q)}{\|v\|_{V}\|q\|_{Q}}\geq\beta.$

The following theorem summarizes the stability and convergence result that is
obtained when the inf-sup conditions are uniform with respect to the parameter
$h$.

######  Theorem 1.

If there exist $\alpha_{0}>0$ and $\beta_{0}>0$ such that the constants in IS1
and IS2 are uniformly bounded below, that is, $\alpha_{h}\geq\alpha_{0}$ and
$\beta_{h}\geq\beta_{0}$ for all $h$, then the following quasi-optimal error
estimate holds true

(9) $\|u-u_{h}\|_{V}+\|p-p_{h}\|_{Q}\leq C\inf_{\begin{subarray}{c}v_{h}\in
V_{h}\\\ q_{h}\in
Q_{h}\end{subarray}}\left(\|u-v_{h}\|_{V}+\|p-q_{h}\|_{Q}\right),$

where $(u,p)$ and $(u_{h},p_{h})$ are the solutions of (1) and (2),
respectively.

######  Remark 2.

The constant $C$ in (9) could be made explicit in terms of $\alpha_{0}$ and
$\beta_{0}$. The interested reader is referred to [6, Theorem 5.2.1].

##  3\. A non-standard mixed formulation for the Laplace equation

In this paper we consider the homogeneous Dirichlet problem for the Laplace
equation on a polygonal domain in $\mathbb{R}^{2}$: given $f$ find $u$ such
that

(10) $\left\\{\begin{aligned} &-\Delta u=f&&\text{in }\Omega\\\ &u=0&&\text{on
}\partial\Omega.\end{aligned}\right.$

More general elliptic equations might be considered, but we believe that this
is the simplest and most effective setting for the presentation of our
results.

We split the second order equation as a system of two first order equations by
introducing the variable $\boldsymbol{\sigma}=-\operatorname{\nabla}u$. As
opposed to the standard mixed formulation, we integrate by parts the
equilibrium equation and not the equation defining $\boldsymbol{\sigma}$, and
obtain the following problem: given $f$ in $H^{-1}(\Omega)$, find
$\boldsymbol{\sigma}\in L^{2}(\Omega)^{2}$ and $u\in H^{1}_{0}(\Omega)$ such
that

(11) $\left\\{\begin{aligned}
&(\boldsymbol{\sigma},\boldsymbol{\tau})+(\boldsymbol{\tau},\operatorname{\nabla}u)=0&&\forall\boldsymbol{\tau}\in
L^{2}(\Omega)^{2}\\\ &(\boldsymbol{\sigma},\operatorname{\nabla}v)=-\langle
f,v\rangle&&\forall v\in H^{1}_{0}(\Omega).\end{aligned}\right.$

Sometimes this formulation is called dual mixed formulation to differentiate
it from the standard primal mixed formulation.

######  Remark 3.

Our interest in problem (11) is related to an analogue formulation used in
elasticity (see [12]). A numerical study of the inf-sup condition in the case
of the Laplace equation was presented in [4].

The mixed formulation is well posed in the chosen functional spaces. For
completeness, this is proved in the next theorem.

######  Theorem 4.

The mixed formulation presented in (11) is well posed in the sense that the
inf-sup conditions (6), (7), and (8) are satisfied.

###### Proof.

The inf-sup conditions (6) and (7) follow from the fact that the bilinear form
$a(\cdot,\cdot)$ corresponds to the identity operator in $L^{2}(\Omega)^{2}$,
which is clearly invertible on the whole space. The inf-sup condition (8) for
the bilinear form $b(\cdot,\cdot)$ follows from the inclusion
$\operatorname{\nabla}H^{1}_{0}(\Omega)\subset L^{2}(\Omega)^{2}$: given $v\in
H^{1}_{0}(\Omega)$, the vectorfield $\boldsymbol{\tau}=\operatorname{\nabla}v$
satisfies
$b(\boldsymbol{\tau},v)=\|\operatorname{\nabla}v\|^{2}_{L^{2}(\Omega)}\geq
C_{1}\|v\|^{2}_{H^{1}_{0}(\Omega)}$ and
$\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}=\|\operatorname{\nabla}v\|_{L^{2}(\Omega)}\leq
C_{2}\|v\|_{H^{1}_{0}(\Omega)}$, that is (8) with $\beta=C_{1}/C_{2}$.

∎

Given two discrete subspaces $\Sigma_{h}\subset L^{2}(\Omega)^{2}$ and
$U_{h}\subset H^{1}_{0}(\Omega)$, the discretization or problem (11) reads:
find $\boldsymbol{\sigma}_{h}\in\Sigma_{h}$ and $u_{h}\in U_{h}$ such that

(12) $\left\\{\begin{aligned}
&(\boldsymbol{\sigma}_{h},\boldsymbol{\tau})+(\boldsymbol{\tau},\operatorname{\nabla}u_{h})=0&&\forall\boldsymbol{\tau}\in\Sigma_{h}\\\
&(\boldsymbol{\sigma}_{h},\operatorname{\nabla}v)=-\langle f,v\rangle&&\forall
v\in U_{h}.\end{aligned}\right.$

######  Remark 5.

Since the inf-sup conditions (6) and (7) are satisfied on the entire space
$L^{2}(\Omega)^{2}$, any choice of discrete spaces $\Sigma_{h}$ and $U_{h}$
will satisfy uniformly the discrete inf-sup condition (4). It follows that the
only condition to be shown for the stability of the discretization, is the
uniform bound of the discrete inf-sup constant in (5) associated with the
bilinear form $b(\cdot,\cdot)$.

A natural choice for the discrete spaces is given by discontinuous piecewise
polynomials of degree $k$ (in each component) for $\Sigma_{h}$ and continuous
piecewise polynomials of degree $k+1$ for $U_{h}$. In the next proposition we
state the stability and the quasi-optimal convergence of the resulting scheme.

######  Proposition 6.

For $k\geq 0$ let $\Sigma_{h}$ be the space of discontinuous piecewise
polynomials of degree $k$ in each component and $U_{h}$ be the space of
continuous piecewise polynomials of degree $k+1$ with zero boundary
conditions. Then the approximation (12) of the mixed formulation (11) is
uniformly stable in the sense of Theorem 1 and the following quasi-optimal
error estimate holds true

$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}+\|u-u_{h}\|_{H^{1}_{0}(\Omega)}\leq
C\inf_{\begin{subarray}{c}\boldsymbol{\tau}\in\Sigma_{h}\\\ v_{h}\in
U_{h}\end{subarray}}\left(\|\boldsymbol{\sigma}-\boldsymbol{\tau}\|_{L^{2}(\Omega)}+\|u-v\|_{H^{1}_{0}(\Omega)}\right).$

In particular, if the solution $(\boldsymbol{\sigma},u)$ is smooth enough, we
have

$\displaystyle\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}+\|u-u_{h}\|_{H^{1}_{0}(\Omega)}$
$\displaystyle\leq
Ch^{k+1}\left(\|\boldsymbol{\sigma}\|_{H^{k+1}(\Omega)}+\|u\|_{H^{k+2}(\Omega)}\right)$
$\displaystyle\leq Ch^{k+1}\|u\|_{H^{k+2}(\Omega)}.$

In general, if $u\in H^{1+s}(\Omega)$ with $s\leq k+1$, we obtain

$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}+\|u-u_{h}\|_{H^{1}_{0}(\Omega)}\leq
Ch^{s}\|u\|_{H^{1+s}(\Omega)}.$

###### Proof.

The stability proof follows the same lines as in Theorem 4. In particular, the
uniform inf-sup condition (4) is guaranteed by the global invertibility of the
bilinear form $a(\cdot,\cdot)$ which corresponds to the identity in
$L^{2}(\Omega)$. The uniform inf-sup condition (5) for the bilinear form
$b(\cdot,\cdot)$ follows from the inclusion
$\operatorname{\nabla}U_{h}\subset\Sigma_{h}$.

∎

It can actually be easily observed that if
$\operatorname{\nabla}U_{h}\subset\Sigma_{h}$ then the mixed formulation (12)
is equivalent to the standard Galerkin formulation where the space $U_{h}$ is
used. We state this result in the following proposition.

######  Proposition 7.

If the inclusion $\operatorname{\nabla}U_{h}\subset\Sigma_{h}$ is satisfied,
then the mixed formulation (12) is well posed in the sense of Theorem 1 and
the component $u_{h}$ of its solution solves the standard Galerkin formulation

(13) $(\operatorname{\nabla}u_{h},\operatorname{\nabla}v)=\langle
f,v\rangle\qquad\forall v\in U_{h}.$

The other component of the solution is given by
$\boldsymbol{\sigma}_{h}=\operatorname{\nabla}u_{h}$.

###### Proof.

It is easily seen that inserting
$\boldsymbol{\sigma}_{h}=\operatorname{\nabla}u_{h}$ into (12), the first
equation is an identity and the second equation corresponds precisely to (13).
∎

The lowest-order element presented in Proposition 6 will be referred to as the
$P_{0}-P_{1}$ scheme.

In a more general context, [12] discusses the approximation of a linear
elasticity problem with a mixed scheme which has some analogies with our
formulation (11). For particular reasons related to some equilibration
properties that will be made more precise later on, it is proposed the use of
Raviart–Thomas elements $RT_{0}$ for the definition of $\Sigma_{h}$. If
_discontinuous_ $RT_{0}$ elements are used, then the same stability proof as
for the $P_{0}-P_{1}$ scheme applies. This follows from the fact that
$\Sigma_{h}$ contains the space of piecewise constants. We state this result
in the following corollary for any degree $k$.

######  Corollary 8.

Let $\Sigma_{h}$ be the space of discontinuous Raviart–Thomas finite elements
of degree $k\geq 0$ and $U_{h}$ be the space of continuous piecewise
polynomials of degree $k+1$ with homogeneous boundary conditions. Then the
formulation (12) is uniformly stable and the following error estimate holds
true if the solution $u$ is smooth enough

$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}+\|u-u_{h}\|_{H^{1}_{0}(\Omega)}\leq
Ch^{s}\|u\|_{H^{1+s}(\Omega)}$

with $s\leq k+1$. Moreover, the mixed problem is equivalent to the standard
Galerkin approximation (13) with continuous polynomials of degree $k+1$ and
$\boldsymbol{\sigma}_{h}=\operatorname{\nabla}u_{h}$.

We now consider the lowest order case, that is $k=0$, so that the
approximation of $u$ is obtained by standard piecewise linear elements.

For reasons related to the equilibration property
$\operatorname{\mathrm{div}}\boldsymbol{\sigma}=f$, in [12] it is proposed to
use $\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$-conforming $RT_{0}$
elements for the definition of $\Sigma_{h}$. We are going to denote this
element by $RT_{0}-P_{1}$. Numerical evidence seems to indicate that for the
elasticity problem this choice provides a uniformly stable scheme [12], while
this is not the case for the Poisson problem [4]. We are going to analyze in
more detail this element in the next sections.

##  4\. The inf-sup condition for the $RT_{0}-P_{1}$ scheme

The first inf-sup condition IS1 (4) is automatically satisfied for our mixed
formulation (see Remark 2), so that we are only discussing the second inf-sup
IS2 (5) which reads:

(14) $\inf_{v\in
U_{h}}\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}v)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}\|v\|_{H^{1}(\Omega)}}\geq\beta_{h}.$

We start with a positive result, showing that for all $h$ the constant
$\beta_{h}$ is strictly greater than zero.

######  Theorem 9.

For all $v\in U_{h}$ with $v\neq 0$ there exists
$\boldsymbol{\tau}\in\Sigma_{h}$ such that

$(\boldsymbol{\tau},\operatorname{\nabla}v)>0,$

that is the inf-sup constant in (14) satisfies $\beta_{h}>0$ for all $h$.

###### Proof.

After integration by parts and taking into account the boundary conditions, we
have

$(\boldsymbol{\tau},\operatorname{\nabla}v)=-(\operatorname{\mathrm{div}}\boldsymbol{\tau},v).$

We denote by $\Pi_{0}$ the $L^{2}(\Omega)$ projection onto the space of
piecewise constant functions; the term
$(\boldsymbol{\tau},\operatorname{\nabla}v)$ is maximized by taking
$\boldsymbol{\tau}\in\Sigma_{h}$ with
$\operatorname{\mathrm{div}}\boldsymbol{\tau}=\Pi_{0}v$, so that we have

$(\boldsymbol{\tau},\operatorname{\nabla}v)=-\|\Pi_{0}v\|^{2}_{L^{2}(\Omega)}.$

Hence the result follows by observing that $\Pi_{0}v=0$ implies $v=0$ if $v$
is vanishing on $\partial\Omega$. Indeed, it is not possible to construct a
function $v\in U_{h}$ that is zero mean valued in each element. This is easily
seen by starting from a boundary element $T$ (with two vertices on
$\partial\Omega$): if $v$ is zero mean valued on $T$, then necessarily it
vanishes on $T$; the same argument can then be applied to the neighboring
elements sharing an edge with $T$ and so on until it is seen that $v$ must
vanish on all elements of the triangulation. ∎

The immediate consequence of the previous theorem is that problem (12) is
solvable.

######  Corollary 10.

For all $f\in H^{-1}(\Omega)$ and all $h$ there exists a unique solution to
problem (12).

We postpone to Appendix B further theoretical investigations about the
behavior of the inf-sup constant. Here we continue this study numerically.

It is well known that an estimate of the inf-sup constant $\beta_{h}$
appearing in (14) can be obtained by solving an algebraic problem. In [6] a
singular value decomposition is used. An essentially equivalent approach was
described in [9] (see Remark 12), based on an idea from [10].

######  Proposition 11.

Let $\Sigma_{h}$ and $U_{h}$ be finite element spaces and consider the
discrete saddle point problem (12). Let $\mathsf{A}$ and $\mathsf{B}$ be the
corresponding matrices appearing in (3) and introduce the matrix $\mathsf{M}$
corresponding to the $H^{1}_{0}(\Omega)$ inner product
$(\operatorname{\nabla}\cdot,\operatorname{\nabla}\cdot)$ in $U_{h}$. Then the
inf-sup constant $\beta_{h}$ in (14) is equal to $\sqrt{\mu_{min}}$, where
$\mu_{min}$ is the smallest eigenvalue $\mu$ of the generalized eigenvalue
problem

(15)
$\mathsf{B}\mathsf{A}^{-1}\mathsf{B}^{\top}\mathsf{x}=\mu\mathsf{M}\mathsf{x}.$

######  Remark 12.

In [6, Section 3.4.3] it is shown that the constant $\beta_{h}$ is equal to
the smallest singular value of the matrix
$\mathsf{M}^{-1/2}\mathsf{B}\mathsf{A}^{-1/2}$. This statement is equivalent
to Proposition 15 since the singular values of
$\mathcal{A}=\mathsf{M}^{-1/2}\mathsf{B}\mathsf{A}^{-1/2}$ are the square
roots of the eigenvalues of
$\mathcal{A}\mathcal{A}^{\top}=\mathsf{M}^{-1/2}\mathsf{B}\mathsf{A}^{-1}\mathsf{B}^{\top}\mathsf{M}^{-1/2}$.
Please note the typo in [6, Proposition 3.4.5], where the involved matrix
should read $S_{Y}^{-1}BS_{X}^{-1}$ instead of $S_{Y}BS_{X}$.

The eigenvalue problem (15) has $N(h)=\dim(U_{h})$ eigensolutions (taking into
account possibly repeated eigenvalues); we number the eigenvalues starting
from the largest one, so that $\mu_{min}=\mu_{N(h)}$,

$\mu_{1}\geq\mu_{2}\geq\dots\geq\mu_{N(h)}\geq 0$

and we denote the corresponding eigenvectors by
$\\{\mathsf{x}_{i}\\}\subset\mathbb{R}^{N(h)}$. Each eigenvector
$\mathsf{x}_{i}$ represents an element $g_{i}$ of $U_{h}$ and we have

$U_{h}=\text{span}\\{g_{1},\dots,g_{N(h)}\\}.$

The eigenvalue problem (15) associated with the inf-sup condition has also a
mixed equivalent formulation:

$\left(\begin{matrix}\mathsf{A}&\mathsf{B}^{\top}\\\
\mathsf{B}&\mathsf{0}\end{matrix}\right)\left(\begin{matrix}\mathsf{y}\\\
\mathsf{x}\end{matrix}\right)=\mu\left(\begin{matrix}\mathsf{0}&\mathsf{0}\\\
\mathsf{0}&-M\end{matrix}\right)\left(\begin{matrix}\mathsf{y}\\\
\mathsf{x}\end{matrix}\right).$

The two formulations are easily shown equivalent to each other by solving for
$\mathsf{y}$ the first equation
$\mathsf{A}\mathsf{y}+\mathsf{B}^{\top}\mathsf{x}=0$ and substituting into the
second equation.

This allows the definition of associated eigenvectors
$\\{\mathsf{y}_{i}\\}\subset\mathbb{R}^{N(h)}$ in addition to the
$\\{\mathsf{x}_{i}\\}$’s satisfying the relation

$\mathsf{A}\mathsf{y}_{i}+\mathsf{B}^{\top}\mathsf{x}_{i}=0\qquad(i=1,\dots,N(h)).$

Translating into the finite element notation, we have constructed two sets of
finite element functions $\\{\boldsymbol{\sigma}_{i,h}\\}\subset\Sigma_{h}$
and $\\{u_{i,h}\\}\subset U_{h}$ that satisfy the following variational
problem

(16) $\left\\{\begin{aligned}
&(\boldsymbol{\sigma}_{i,h},\boldsymbol{\tau})+(\boldsymbol{\tau},\operatorname{\nabla}u_{i,h})=0&&\forall\boldsymbol{\tau}\in\Sigma_{h}\\\
&(\boldsymbol{\sigma}_{i,h},\operatorname{\nabla}v)=-\mu_{i}(\operatorname{\nabla}u_{i,h},\operatorname{\nabla}v)&&\forall
v\in U_{h}\end{aligned}\right.$

and such that

$\displaystyle\text{span}\\{\boldsymbol{\sigma}_{1,h},\dots,\boldsymbol{\sigma}_{N(h),h}\\}\subset\Sigma_{h}$
$\displaystyle\text{span}\\{u_{1,h},\dots,u_{N(h),h}\\}=U_{h}.$

The eigenvectors $\\{u_{i,h}\\}$ can be chosen so that

$\displaystyle(\operatorname{\nabla}u_{i,h},\operatorname{\nabla}u_{j,h})_{L^{2}(\Omega)}=0$
$\displaystyle\text{if }i\neq j$
$\displaystyle\|\operatorname{\nabla}u_{i,h}\|_{L^{2}(\Omega)}=1.$

It follows that also the $\\{\boldsymbol{\sigma}_{i,h}\\}$ are orthogonal,
since

$(\boldsymbol{\sigma}_{i,h},\boldsymbol{\sigma}_{j,h})_{L^{2}(\Omega)}=-(\boldsymbol{\sigma}_{i,h},\operatorname{\nabla}u_{j,h})_{L^{2}(\Omega)}=\mu_{i}(\operatorname{\nabla}u_{i,h},\operatorname{\nabla}u_{j,h})_{L^{2}(\Omega)}=0\qquad\text{if
}i\neq j.$

In [4] it was shown numerically that the inf-sup constant $\beta_{h}$ is not
uniformly bounded from below. More precisely, $\mu_{N(h)}$ (together with
other eigenvalues of (15)) tends to zero when $h$ goes to zero. We report here
again this behavior and show how it may be different depending on the mesh
sequence that we are considering.

We use three sequences of mesh on the unit square: two structured and one
unstructured. With obvious meaning, we call them Crossed, Right, and Non
Structured mesh, respectively. An example of such meshes is shown in Figure 1.

Figure 1. Crossed, Right, and Non Structured meshes

Table 1 shows the last four computed eigenvalues
$\\{\mu_{N(h)-3},\dots,\mu_{N(h)}\\}$.

Crossed mesh

$N(h)$ | $\mu_{N(h)-3}$ | $\mu_{N(h)-2}$ | $\mu_{N(h)-1}$ | $\mu_{N(h)}$
---|---|---|---|---
13 | 0.66666667 | 0.5 | 0.5 | 0.22222222
41 | 0.16521696 | 0.15643855 | 0.15643855 | 0.06604647
145 | 0.04880971 | 0.04191655 | 0.04191655 | 0.01698587
545 | 0.01268672 | 0.01065182 | 0.01065182 | 0.00427448
2113 | 0.00320245 | 0.00267372 | 0.00267372 | 0.00107035

Right mesh

$N(h)$ | $\mu_{N(h)-3}$ | $\mu_{N(h)-2}$ | $\mu_{N(h)-1}$ | $\mu_{N(h)}$
---|---|---|---|---
25 | 0.44698968 | 0.41649077 | 0.23888594 | 0.23720409
81 | 0.14099494 | 0.14089618 | 0.06715927 | 0.06707865
289 | 0.03714468 | 0.03714446 | 0.01720941 | 0.01720741
1089 | 0.00938762 | 0.00938762 | 0.00432346 | 0.00432341
4225 | 0.00235165 | 0.00235165 | 0.00108154 | 0.00108154

Non Structured mesh

$N(h)$ | $\mu_{N(h)-3}$ | $\mu_{N(h)-2}$ | $\mu_{N(h)-1}$ | $\mu_{N(h)}$
---|---|---|---|---
38 | 0.25105003 | 0.22389003 | 0.16903879 | 0.11779969
140 | 0.10876822 | 0.0991606 | 0.08762351 | 0.05880215
531 | 0.0736162 | 0.06851846 | 0.06604323 | 0.05821584
2066 | 0.05617133 | 0.05523537 | 0.05366541 | 0.04984366
8128 | 0.04867193 | 0.04780397 | 0.0451059 | 0.04200516

Table 1. Numerical estimate of the inf-sup constant for $RT_{0}-P_{1}$ scheme

It turns out that in the case of the structured meshes the smallest eigenvalue
$\mu_{N(h)}$ goes to zero quadratically in $h$, thus giving the estimate
$\beta_{h}=O(h)$. The behavior of the inf-sup constant on the unstructured
mesh is less critical, even if not optimal, showing a slower decay of
$\mu_{N(h)}$ as $h$ goes to zero.

It interesting to look at the eigenfunctions corresponding to (16). Figure 2
shows the eigenfunctions $u_{1,h}$ and $u_{N(h),h}$ corresponding to the
maximum and minimum eigenvalue for the three meshes.

Crossed mesh

Right mesh

Non Structured mesh

Figure 2. Eigenfunctions associated with first (left) and last (right)
eigenvalues $\mu_{1}$ and $\mu_{N(h)}$ of (16)

It is apparent that the eigenfunctions corresponding to the smallest
eigenvalue are highly oscillatory, while the ones associated to the largest
one don’t change their sign. This fact can be made more precise by identifying
an appropriate subspace of $U_{h}$ for which a uniform inf-sup condition holds
true. This will be done in the next section where we discuss possible
splittings of the spaces in order to separate the _stable_ part of the
solution from the _unstable_ one.

##  5\. Stable and unstable subspaces

The following discussion identifies special subspaces of $U_{h}$ for which a
uniform inf-sup condition holds true. According to what we have seen in the
previous section, we are expecting that the degeneracy of the inf-sup constant
is associated with highly oscillatory eigenfunctions of problem (16).

Let $w_{k}\in U_{h}$, $k=1,\dots,N(h)$, be the $k$-th discrete eigenfunction
of the Laplace operator approximated by the standard Galerkin method, that is

$(\operatorname{\nabla}w_{k},\operatorname{\nabla}v)=\lambda_{k}(w_{k},v)\quad\forall
v\in U_{h}$

with $0<\lambda_{1}<\lambda_{2}\leq\cdots$.

Our first result shows that if we fix $\bar{N}$ and restrict $U_{h}$ to the
space spanned by $\\{w_{1},w_{2},\dots,w_{\bar{N}}\\}$, then the inf-sup
condition holds with a constant uniform in $h$.

######  Theorem 13.

Let $\bar{W}_{h}$ be the subspace of $U_{h}$ spanned by
$\\{w_{1},w_{2},\dots,w_{\bar{N}}\\}$. Then there exists $\bar{\beta}>0$
independent of $h$ such that

(17)
$\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}w)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}\geq\bar{\beta}\qquad\forall
w\in\bar{W}_{h}.$

###### Proof.

First we show that (17) is satisfied when $w=w_{k}$ for a fixed $k$. We define
$\boldsymbol{\varsigma}_{k}$ as the solution of the following mixed problem:
find $\boldsymbol{\sigma}_{k}\in\Sigma_{h}$ and $p$ in the space of piecewise
constant functions $\mathcal{P}_{0}$ such that

$\left\\{\begin{aligned}
&(\boldsymbol{\varsigma}_{k},\boldsymbol{\tau})+(\operatorname{\mathrm{div}}\boldsymbol{\tau},p)=0&&\forall\boldsymbol{\tau}\in\Sigma_{h}\\\
&(\operatorname{\mathrm{div}}\boldsymbol{\varsigma}_{k},q)=-\lambda_{k}(w_{k},q)&&\forall
q\in\mathcal{P}_{0}.\end{aligned}\right.$

We then have

$\displaystyle\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}\leq
C_{\sigma}\lambda_{k}\|w_{k}\|_{L^{2}(\Omega)}=C_{\sigma}\lambda_{k}^{1/2}\|\operatorname{\nabla}w_{k}\|_{L^{2}(\Omega)}$
$\displaystyle\operatorname{\mathrm{div}}\boldsymbol{\varsigma}_{k}=-\lambda_{k}\Pi_{0}w_{k},$

where $\Pi_{0}$ is the $L^{2}$ projection onto $\mathcal{P}_{0}$.

By choosing $\boldsymbol{\tau}=\boldsymbol{\varsigma}_{k}$ in (17) we can
conclude

(18)
$\displaystyle\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}w_{k})}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}$
$\displaystyle\geq\frac{(\boldsymbol{\varsigma}_{k},\operatorname{\nabla}w_{k})}{\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}}=-\frac{(\operatorname{\mathrm{div}}\boldsymbol{\varsigma}_{k},w_{k})}{\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}}=\frac{\lambda_{k}(\Pi_{0}w_{k},w_{k})}{\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}}$
$\displaystyle=\frac{\lambda_{k}(w_{k},w_{k})-\lambda_{k}(w_{k}-\Pi_{0}w_{k},w_{k})}{\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}}$
$\displaystyle=\frac{\|\operatorname{\nabla}w_{k}\|^{2}_{L^{2}(\Omega)}-\lambda_{k}(w_{k}-\Pi_{0}w_{k},w_{k}-\Pi_{0}w_{k})}{\|\boldsymbol{\varsigma}_{k}\|_{L^{2}(\Omega)}}$
$\displaystyle\geq\frac{\|\operatorname{\nabla}w_{k}\|^{2}_{L^{2}(\Omega)}-\lambda_{k}\|w_{k}-\Pi_{0}w_{k}\|^{2}_{L^{2}(\Omega)}}{C_{\sigma}\lambda_{k}^{1/2}\|\operatorname{\nabla}w_{k}\|_{L^{2}(\Omega)}}$
$\displaystyle\geq\frac{\|\operatorname{\nabla}w_{k}\|^{2}_{L^{2}(\Omega)}-\lambda_{k}C^{2}_{\Pi}h^{2}\|\operatorname{\nabla}w_{k}\|^{2}_{L^{2}(\Omega)}}{C_{\sigma}\lambda_{k}^{1/2}\|\operatorname{\nabla}w_{k}\|_{L^{2}(\Omega)}}$
$\displaystyle=\frac{1-C^{2}_{\Pi}h^{2}\lambda_{k}}{C_{\sigma}\lambda_{k}^{1/2}}\|\operatorname{\nabla}w_{k}\|_{L^{2}(\Omega)}=\beta_{k}\|\operatorname{\nabla}w_{k}\|_{L^{2}(\Omega)},$

where $C_{\Pi}$ denotes the constant appearing in the approximation property
of $\Pi_{0}$

$\|f-\Pi_{0}f\|_{L^{2}(\Omega)}\leq
C_{\Pi}h\|\operatorname{\nabla}f\|_{L^{2}(\Omega)}.$

It follows that $\beta_{k}$ is bounded below uniformly in $h$ for $k$ fixed
and $h$ small enough.

In order to complete the proof it remains to extend the result to a generic
$w\in\bar{W}_{h}$; the restriction of $h$ small enough is removed by comparing
with Theorem 9.

We detail how to deal with $w=w_{i}+w_{j}$; the generic result follows with
similar arguments by considering a finite linear combination of discrete
eigenfunctions. We define
$\boldsymbol{\tau}=\boldsymbol{\varsigma}=\boldsymbol{\varsigma}_{i}+\boldsymbol{\varsigma}_{j}$
in (17), where $\boldsymbol{\varsigma}_{k}$ ($k=i,j$) is defined in the
previous step. We have

(19)
$\displaystyle\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}w)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}$
$\displaystyle\geq\frac{(\boldsymbol{\varsigma},\operatorname{\nabla}w)}{\|\boldsymbol{\varsigma}\|_{L^{2}(\Omega)}}=\frac{(\lambda_{i}\Pi_{0}w_{i}+\lambda_{j}\Pi_{0}w_{j},w)}{\|\boldsymbol{\varsigma}\|_{L^{2}(\Omega)}}$
$\displaystyle\geq\frac{\|\operatorname{\nabla}w\|^{2}_{L^{2}(\Omega)}-(\lambda_{i}(w_{i}-\Pi_{0}w_{i})+\lambda_{j}(w_{j}-\Pi_{0}w_{j}),w)}{C_{\sigma}\|\lambda_{i}w_{i}+\lambda_{j}w_{j}\|_{L^{2}(\Omega)}}.$

Let us study separately the numerator and the denominator of the last
expression by starting with the scalar product appearing in the numerator. In
order to make the notation shorted, we denote by $f_{k}$ ($k=i,j$) the term
$\lambda_{k}w_{k}$.

$\displaystyle(f_{i}-\Pi_{0}f_{i}$
$\displaystyle+f_{j}-\Pi_{0}f_{j},w_{i}+w_{j})$
$\displaystyle=(f_{i}-\Pi_{0}f_{i},w_{i})+(f_{j}-\Pi_{0}f_{j},w_{j})$
$\displaystyle\quad+(f_{i}-\Pi_{0}f_{i},w_{j})+(f_{j}-\Pi_{0}f_{j},w_{i})$
$\displaystyle\leq
C_{\Pi}h\left(\|\operatorname{\nabla}f_{i}\|_{L^{2}(\Omega)}\|w_{i}\|_{L^{2}(\Omega)}+\|\operatorname{\nabla}f_{j}\|_{L^{2}(\Omega)}\|w_{j}\|_{L^{2}(\Omega)}\right)$
$\displaystyle\quad+(f_{i}-\Pi_{0}f_{i},w_{j}-\Pi_{0}w_{j})+(f_{j}-\Pi_{0}f_{j},w_{i}-\Pi_{0}w_{i})$
$\displaystyle\leq
C_{\Pi}h\left(\max(\lambda_{i},\lambda_{j})^{1/2}\|\operatorname{\nabla}w\|_{L^{2}(\Omega)}\right)$
$\displaystyle\quad+C^{2}_{\Pi}h^{2}\left(\lambda_{i}\|\operatorname{\nabla}w_{i}\|_{L^{2}(\Omega)}\|\operatorname{\nabla}w_{j}\|_{L^{2}(\Omega)}+\lambda_{j}\|\operatorname{\nabla}w_{j}\|_{L^{2}(\Omega)}\|\operatorname{\nabla}w_{i}\|_{L^{2}(\Omega)}\right)$
$\displaystyle\leq
C_{\Pi}h\left(\max(\lambda_{i},\lambda_{j})^{1/2}\|\operatorname{\nabla}w\|_{L^{2}(\Omega)}\right)$
$\displaystyle\quad+(C^{2}_{\Pi}h^{2}/2)\max(\lambda_{i},\lambda_{j})\|\operatorname{\nabla}w\|^{2}_{L^{2}(\Omega)},$

where we used twice the orthogonality of $\operatorname{\nabla}w_{i}$ and
$\operatorname{\nabla}w_{j}$. It follows that the numerator in (19) can be
bounded below by a positive constant times
$\|\operatorname{\nabla}w\|^{2}_{L^{2}(\Omega)}$ for $h$ small enough. For the
denominator, we have

$\displaystyle\|\lambda_{i}w_{i}+\lambda_{j}w_{j}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle=\lambda_{i}^{2}\|w_{i}\|^{2}_{L^{2}(\Omega)}+\lambda_{j}^{2}\|w_{j}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle=\lambda_{i}\|\operatorname{\nabla}w_{i}\|^{2}_{L^{2}(\Omega)}+\lambda_{j}\|\operatorname{\nabla}w_{j}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle\leq\max(\lambda_{i},\lambda_{j})\|\operatorname{\nabla}w\|^{2}_{L^{2}(\Omega)},$

where we used again the orthogonality of $\operatorname{\nabla}w_{i}$ and
$\operatorname{\nabla}w_{j}$ together with the orthogonality of $w_{i}$ and
$w_{j}$. It follows that for $h$ small enough there exists a constant $C$
independent of $h$, but dependent on $i$ and $j$, such that

$\frac{(\boldsymbol{\varsigma},\operatorname{\nabla}w)}{\|\boldsymbol{\varsigma}\|_{L^{2}(\Omega)}}\geq
C\|\operatorname{\nabla}w\|_{L^{2}(\Omega)}.$

∎

######  Remark 14.

The inf-sup constant of the previous theorem depends on the dimension of
$\bar{W}_{h}$. In particular, it gives a confirmation that the components of
$U_{h}$ that are source of instability are associated with highly oscillating
functions. We now discuss how it is possible to split the solution of problem
(12) into a stable part and an unstable one in a more abstract way.

It can be easily seen that all eigenvalues of (15) are not larger than $1$, so
that we fix a threshold value $\underline{\mu}$ between $0$ and $1$ and define
an index $N_{\underline{\mu}}$ (depending on $h$) so that

$\displaystyle\mu_{i}\geq\underline{\mu}$ $\displaystyle\forall i\leq
N_{\underline{\mu}}$ $\displaystyle\mu_{i}<\underline{\mu}$
$\displaystyle\forall i>N_{\underline{\mu}}.$

We can then introduce the following spaces

(20)
$\displaystyle\Sigma_{1}=\text{span}\\{\boldsymbol{\sigma}_{1,h},\dots,\boldsymbol{\sigma}_{N_{\underline{\mu}},h}\\}$
$\displaystyle\Sigma_{2}=\text{span}\\{\boldsymbol{\sigma}_{N_{\underline{\mu}}+1,h},\dots,\boldsymbol{\sigma}_{N(h),h}\\}$
$\displaystyle U_{1}=\text{span}\\{u_{1,h},\dots,u_{N_{\underline{\mu}},h}\\}$
$\displaystyle
U_{2}=\text{span}\\{u_{N_{\underline{\mu}}+1,h},\dots,u_{N(h),h}\\}.$

It follows that $\Sigma_{h}=\Sigma_{1}\oplus\Sigma_{2}\oplus\tilde{\Sigma}$
and that $U_{h}=U_{1}\oplus U_{2}$, where $\tilde{\Sigma}$ is a remainder
space whose dimension is equal to $\dim\Sigma_{h}-N(h)\geq 0$. This space
satisfies the orthogonality

$(\boldsymbol{\sigma}_{i,h},\boldsymbol{\tau})_{L^{2}(\Omega)}=0\quad\forall
i\ \forall\boldsymbol{\tau}\in\tilde{\Sigma}.$

######  Remark 15.

It turns out that the subspaces $\Sigma_{1}-U_{1}$ provide a stable
discretization of (12) since the corresponding inf-sup constant is associated
with $(\underline{\mu})^{1/2}>0$ which is independent of $h$.

We can now look at the matrix form of our discrete problem (12)

$\left(\begin{matrix}\mathsf{A}&\mathsf{B}^{\top}\\\
\mathsf{B}&\mathsf{0}\end{matrix}\right)\left(\begin{matrix}\mathsf{y}\\\
\mathsf{x}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
-\mathsf{f}\end{matrix}\right).$

The above problem can be written in the following form by using the splitting
of the spaces $\Sigma_{h}$ and $U_{h}$:

$\left(\begin{matrix}\mathsf{A}_{11}&\mathsf{A}_{21}^{\top}&\mathsf{A}_{31}^{\top}&\mathsf{B}_{11}^{\top}&\mathsf{B}_{21}^{\top}\\\
\mathsf{A}_{21}&\mathsf{A}_{22}&\mathsf{A}_{32}^{\top}&\mathsf{B}_{12}^{\top}&\mathsf{B}_{22}^{\top}\\\
\mathsf{A}_{31}&\mathsf{A}_{32}&\mathsf{A}_{33}&\mathsf{B}_{13}^{\top}&\mathsf{B}_{23}^{\top}\\\
\mathsf{B}_{11}&\mathsf{B}_{12}&\mathsf{B}_{13}&\mathsf{0}&\mathsf{0}\\\
\mathsf{B}_{21}&\mathsf{B}_{22}&\mathsf{B}_{23}&\mathsf{0}&\mathsf{0}\\\
\end{matrix}\right)\left(\begin{matrix}\mathsf{y}_{1}\\\ \mathsf{y}_{2}\\\
\mathsf{y}_{3}\\\ \mathsf{x}_{1}\\\
\mathsf{x}_{2}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
\mathsf{0}\\\ \mathsf{0}\\\ -\mathsf{f}_{1}\\\
-\mathsf{f}_{2}\end{matrix}\right).$

Some of the matrices involved with this formulation are vanishing. In
particular, it is easy to see that $\mathsf{A}_{21}=\mathsf{0}$ due to the
orthogonalities of the $\operatorname{\nabla}u_{i}$’s; analogously,
$\mathsf{B}_{12}=\mathsf{0}$ and $\mathsf{B}_{21}=\mathsf{0}$. Moreover, due
to the orthogonality of $\tilde{\Sigma}$ with the rest of $\Sigma_{h}$, it
follows that $\mathsf{A}_{31}=\mathsf{0}$ and $\mathsf{A}_{32}=\mathsf{0}$,
and also that $\mathsf{B}_{13}=\mathsf{0}$ and $\mathsf{B}_{23}=\mathsf{0}$.

Hence the matrix problem is reduced to

$\left(\begin{matrix}\mathsf{A}_{11}&\mathsf{0}&\mathsf{0}&\mathsf{B}_{11}^{\top}&\mathsf{0}\\\
\mathsf{0}&\mathsf{A}_{22}&\mathsf{0}&\mathsf{0}&\mathsf{B}_{22}^{\top}\\\
\mathsf{0}&\mathsf{0}&\mathsf{A}_{33}&\mathsf{0}&\mathsf{0}\\\
\mathsf{B}_{11}&\mathsf{0}&\mathsf{0}&\mathsf{0}&\mathsf{0}\\\
\mathsf{0}&\mathsf{B}_{22}&\mathsf{0}&\mathsf{0}&\mathsf{0}\\\
\end{matrix}\right)\left(\begin{matrix}\mathsf{y}_{1}\\\ \mathsf{y}_{2}\\\
\mathsf{y}_{3}\\\ \mathsf{x}_{1}\\\
\mathsf{x}_{2}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
\mathsf{0}\\\ \mathsf{0}\\\ -\mathsf{f}_{1}\\\
-\mathsf{f}_{2}\end{matrix}\right).$

Since $\mathsf{A}_{33}$ is invertible, it follows that $\mathsf{y}_{3}=0$ so
that the system can be reduced to the following compact matrix form

$\left(\begin{matrix}\mathsf{A}_{11}&\mathsf{0}&\mathsf{B}_{11}^{\top}&\mathsf{0}\\\
\mathsf{0}&\mathsf{A}_{22}&\mathsf{0}&\mathsf{B}_{22}^{\top}\\\
\mathsf{B}_{11}&\mathsf{0}&\mathsf{0}&\mathsf{0}\\\
\mathsf{0}&\mathsf{B}_{22}&\mathsf{0}&\mathsf{0}\\\
\end{matrix}\right)\left(\begin{matrix}\mathsf{y}_{1}\\\ \mathsf{y}_{2}\\\
\mathsf{x}_{1}\\\
\mathsf{x}_{2}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
\mathsf{0}\\\ -\mathsf{f}_{1}\\\ -\mathsf{f}_{2}\end{matrix}\right)$

where all the non-vanishing blocks are square and diagonal.

It follows that the system decouples into the following two systems:

$\left(\begin{matrix}\mathsf{A}_{11}&\mathsf{B}_{11}^{\top}\\\
\mathsf{B}_{11}&\mathsf{0}\end{matrix}\right)\left(\begin{matrix}\mathsf{y}_{1}\\\
\mathsf{x}_{1}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
-\mathsf{f}_{1}\end{matrix}\right)$

and

$\left(\begin{matrix}\mathsf{A}_{22}&\mathsf{B}_{22}^{\top}\\\
\mathsf{B}_{22}&\mathsf{0}\end{matrix}\right)\left(\begin{matrix}\mathsf{y}_{2}\\\
\mathsf{x}_{2}\end{matrix}\right)=\left(\begin{matrix}\mathsf{0}\\\
-\mathsf{f}_{2}\end{matrix}\right)$

This corresponds to the following two variational problems. Find
$\boldsymbol{\varphi}_{1,h}\in\Sigma_{1}$ and $w_{1,h}\in U_{1}$ such that

(21) $\left\\{\begin{aligned}
&(\boldsymbol{\varphi}_{1,h},\boldsymbol{\tau})+(\boldsymbol{\tau},\operatorname{\nabla}w_{1,h})=0&&\forall\boldsymbol{\tau}\in\Sigma_{1}\\\
&(\boldsymbol{\varphi}_{1,h},\operatorname{\nabla}v)=-\langle
f,v\rangle&&\forall v\in U_{1}\end{aligned}\right.$

and find $\boldsymbol{\varphi}_{2,h}\in\Sigma_{2}$ and $w_{2,h}\in U_{2}$ such
that

(22) $\left\\{\begin{aligned}
&(\boldsymbol{\varphi}_{2,h},\boldsymbol{\tau})+(\boldsymbol{\tau},\operatorname{\nabla}w_{2,h})=0&&\forall\boldsymbol{\tau}\in\Sigma_{2}\\\
&(\boldsymbol{\varphi}_{2,h},\operatorname{\nabla}v)=-\langle
f,v\rangle&&\forall v\in U_{2}\end{aligned}\right.$

The solution of (12) is then given simply by

$\boldsymbol{\sigma}_{h}=\boldsymbol{\varphi}_{1,h}+\boldsymbol{\varphi}_{2,h}\qquad
u_{h}=w_{1,h}+w_{2,h}.$

The heuristic idea of the splitting is that (21) corresponds to a _stable_
problem, where the inf-sup constant is bounded below by
$(\underline{\mu})^{1/2}$, while (22) has an inf-sup constant that goes to
zero as $h$ goes to zero.

We state in the following proposition the results obtained so far.

######  Proposition 16.

Let us consider the solution
$(\boldsymbol{\sigma}_{h},u_{h})\in\Sigma_{h}\times U_{h}$ of (12). Let
$\underline{\mu}$ be a constant between $0$ and $1$ and consider the subspaces
$\Sigma_{1}$, $\Sigma_{2}$, $U_{1}$, and $U_{2}$ introduced in (20). Then it
holds

$\displaystyle\boldsymbol{\sigma}_{h}=\boldsymbol{\varphi}_{1,h}+\boldsymbol{\varphi}_{2,h}$
$\displaystyle u_{h}=w_{1,h}+w_{2,h},$

where $(\boldsymbol{\varphi}_{1,h},w_{1,h})\in\Sigma_{1}\times U_{1}$ and
$(\boldsymbol{\varphi}_{2,h},w_{2,h})\in\Sigma_{2}\times U_{2}$ solve (21) and
(22), respectively.

The characterization of the solution $u_{h}$ can be pushed further by looking
at the matrices involved with (21) and (22) and by testing the systems with
$v=u_{i,h}$. Indeed, it turns out that $\mathsf{A}_{11}=-\mathsf{B}_{11}$ and
$\mathsf{A}_{22}=-\mathsf{B}_{22}$ have diagonal entries equal to the
eigenvalues $\mu_{i}$ if the bases of the spaces are chosen as in (20). The
following theorem, which is an immediate consequence of the previous
observations, shows how $u_{h}$ can be represented in terms of $f$ and the
basis of $U_{h}$.

######  Theorem 17.

The second component $u_{h}\in U_{h}$ of the solution to (12) can be
represented as

(23) $u_{h}=\sum_{i=1}^{N(h)}\frac{\alpha_{i}}{\mu_{i}}u_{i,h}$

where the coefficients $\alpha_{i}$ are defined as

(24) $\alpha_{i}=\langle f,u_{i,h}\rangle.$

######  Remark 18.

It is interesting to observe that, thanks to the orthogonalities of the
$\\{u_{i,h}\\}$, the coefficients $\\{\alpha_{i}\\}$ in Theorem 24 can be used
also to define the solution $u_{G}\in U_{h}$ of the standard Galerkin method

$(\operatorname{\nabla}u_{G},\operatorname{\nabla}v)=\langle
f,v\rangle\quad\forall v\in U_{h}$

as

(25) $u_{G}=\sum_{i=1}^{N(h)}\alpha_{i}u_{i,h}.$

It follows that $u_{G}$ and $u_{h}$ have similar representations and that the
$i$-th coefficient of their representations differs by a factor equal to
$1/\mu_{i}$.

In light of the previous results, it is interesting to see how the
coefficients $\\{\alpha_{i}\\}$ behave. To this aim we consider different
right hand sides $f$ in the unit square and compare their behavior.

We start with a smooth $f$ equal to the first eigenfunction of the Poisson
problem, then we take a constant right and side and we conclude with two
approximations of the Dirac delta function: one centered at $(1/3,1/5)$ and
the other centered at the center of the domain. The corresponding results are
reported in Figure 3.

$f(x,y)=2\pi^{2}\sin(\pi x)\sin(\pi y)$: Crossed mesh (left), Right mesh
(middle), Non Structured mesh (right)

(a) *

$f(x,y)=1$: Crossed mesh (left), Right mesh (middle), Non Structured mesh
(right)

(b) *

$f(x,y)$ equal to an approximation to the Dirac delta function centered at
$(1/3,1/5)$: Crossed mesh (left), Right mesh (middle), Non Structured mesh
(right)

(c) *

$f(x,y)$ equal to an approximation to the Dirac delta function centered at
$(1/2,1/2)$: Crossed mesh (left), Right mesh (middle), Non Structured mesh
(right)

(d) *

Figure 3. Coefficients $\\{\alpha_{i}\\}$ from (24) for various choices of $f$

As it is natural from the oscillatory behavior of the eigenfunctions, the
coefficients $\alpha_{i}$ corresponding to smooth functions are different from
zero only for few low values of $i$ and vanish for larger $i$’s, while non
smooth functions contain nonzero coefficients in the right most part of the
spectrum. Moreover, the eigenvalues $\mu_{i}$ are close to one for small
values of $i$ while become smaller as we move to the right in the spectrum.

Hence, by comparing the two representation of $u_{h}$ and $U_{G}$ in (23) and
(25), it should be expected that their difference is small when $f$ is smooth
and possibly large when $f$ contains large components along the last part of
the spectrum. Indeed in the next section we are going to analyze the
convergence of the scheme when $f$ is smooth enough.

##  6\. Convergence of $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|$ for
smooth data

Let $\mathcal{P}_{0}$ be the finite element space of piecewise constant
functions. We introduce the following two subspaces of $\Sigma_{h}$:

$\displaystyle
Z_{h}^{m}(f)=\\{\boldsymbol{\tau}_{h}\in\Sigma_{h}:(\operatorname{\mathrm{div}}\boldsymbol{\tau}_{h},q_{h})=(f,q_{h})\
\forall q_{h}\in\mathcal{P}_{0}\\}$ $\displaystyle
Z_{h}(f)=\\{\boldsymbol{\tau}_{h}\in\Sigma_{h}:(\boldsymbol{\tau}_{h},\operatorname{\nabla}v_{h})=-(f,v_{h})\
\forall v_{h}\in U_{h}\\}.$

We have the following crucial result.

######  Lemma 19.

If $f$ belongs to $\mathcal{P}_{0}$, then $Z_{h}^{m}(f)\subset Z_{h}(f)$.
Moreover, if
$\boldsymbol{\sigma}\in\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$
satisfies $\operatorname{\mathrm{div}}\boldsymbol{\sigma}=f$, then

$\inf_{\boldsymbol{\tau}_{h}\in
Z_{h}(f)}\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}\leq
C\|\boldsymbol{\sigma}-\boldsymbol{\sigma}^{I}\|_{L^{2}(\Omega)},$

where $\boldsymbol{\sigma}^{I}\in\Sigma_{h}$ is the interpolant of
$\boldsymbol{\sigma}$.

###### Proof.

If $\operatorname{\mathrm{div}}\boldsymbol{\sigma}$ is in $\mathcal{P}_{0}$,
then $\boldsymbol{\sigma}$ is smooth enough so that its interpolant
$\boldsymbol{\sigma}^{I}\in\Sigma_{h}$ is well defined. In this setting the
definition of $Z_{h}(f)$ reads

$Z_{h}(f)=\\{\boldsymbol{\tau}_{h}\in\Sigma_{h}:(\operatorname{\mathrm{div}}\boldsymbol{\tau}_{h},\Pi_{0}v_{h})=(f,\Pi_{0}v_{h})\
\forall v_{h}\in U_{h}\\},$

where $\Pi_{0}$ denotes the $L^{2}(\Omega)$ projection onto $\mathcal{P}_{0}$.
It follows that if $\boldsymbol{\tau}_{h}$ belongs to $Z_{h}^{m}(f)$ then it
is also in $Z_{h}(f)$.

From the standard inf-sup condition that is valid for the spaces $\Sigma_{h}$
and $\mathcal{P}_{0}$, Proposition 5.1.3 of [6] gives

$\inf_{\boldsymbol{\tau}_{h}\in
Z_{h}^{m}(f)}\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{\mathbf{H}(\operatorname{\mathrm{div}};\Omega)}\leq
C\inf_{\boldsymbol{\tau}_{h}\in\Sigma_{h}}\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{\mathbf{H}(\operatorname{\mathrm{div}};\Omega)},$

so that the inclusion $Z_{h}^{m}(f)\subset Z_{h}(f)$ gives, in particular,

$\inf_{\boldsymbol{\tau}_{h}\in
Z_{h}(f)}\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}\leq
C\inf_{\boldsymbol{\tau}_{h}\in\Sigma_{h}}\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{\mathbf{H}(\operatorname{\mathrm{div}};\Omega)}.$

From the commuting diagram property
$\operatorname{\mathrm{div}}\boldsymbol{\sigma}^{I}=\Pi_{0}\operatorname{\mathrm{div}}\boldsymbol{\sigma}$,
observing that $\operatorname{\mathrm{div}}\boldsymbol{\sigma}$ belongs to
$\mathcal{P}_{0}$, we have that
$\operatorname{\mathrm{div}}(\boldsymbol{\sigma}-\boldsymbol{\sigma}^{I})=0$,
so that we can conclude that the right hand side can be estimated by
$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}^{I}\|_{L^{2}(\Omega)}$. ∎

The following theorem shows the convergence of the approximation of
$\boldsymbol{\sigma}$ given by (12) in the case when $f$ belongs to
$\mathcal{P}_{0}$.

######  Theorem 20.

Let $\boldsymbol{\sigma}\in\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ be
the first component of the solution of (11) for $f\in\mathcal{P}_{0}$ and
$\boldsymbol{\sigma}_{h}$ the corresponding approximation given by (12). Then
we have the estimate

$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C\left(\|\boldsymbol{\sigma}-\boldsymbol{\sigma}^{I}\|_{L^{2}(\Omega)}+\inf_{v_{h}\in
U_{h}}\|\operatorname{\nabla}(u-v_{h})\|_{L^{2}(\Omega)}\right),$

where $\boldsymbol{\sigma}^{I}\in\Sigma_{h}$ is the interpolant of
$\boldsymbol{\sigma}$.

###### Proof.

We estimate
$\|\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}$ for
$\boldsymbol{\tau}_{h}\in Z_{h}(f)$. The result will then follow from the
triangular inequality and Lemma 19.

From the error equation and the properties of $Z_{h}(f)$ we have for all
$v_{h}\in U_{h}$

$\displaystyle\|\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=(\boldsymbol{\sigma}_{h}-\boldsymbol{\sigma},\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h})+(\boldsymbol{\sigma}-\boldsymbol{\tau}_{h},\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h})$
$\displaystyle=(\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h},\operatorname{\nabla}(u-u_{h}))+(\boldsymbol{\sigma}-\boldsymbol{\tau}_{h},\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h})$
$\displaystyle=(\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h},\operatorname{\nabla}(u-v_{h}))+(\boldsymbol{\sigma}-\boldsymbol{\tau}_{h},\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h})$
$\displaystyle\leq\|\boldsymbol{\sigma}_{h}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}\left(\|\operatorname{\nabla}(u-v_{h})\|_{L^{2}(\Omega)}+\|\boldsymbol{\sigma}-\boldsymbol{\tau}_{h}\|_{L^{2}(\Omega)}\right),$

from which we obtain the required estimate. ∎

######  Corollary 21.

If $f\in H^{1}(\Omega)$ and if the Poisson problem has $H^{1+s}(\Omega)$
regularity for some $s\in(0,1]$, then the following optimal error estimate
holds true

$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
Ch^{s}\|f\|_{H^{1}(\Omega)}.$

###### Proof.

The result follows by approximating $f$ with $\Pi_{0}f$ and using the previous
theorem together with the standard estimate

$\|f-\Pi_{0}f\|_{L^{2}(\Omega)}\leq
Ch\|\operatorname{\nabla}f\|_{L^{2}(\Omega)}.$

∎

##  7\. Numerical experiments

In this section we report some numerical results related to the solution of
problem (12). As investigated in Section 5 for the variable $u$ and proved in
Section 6 for the variable $\boldsymbol{\sigma}$, we are expecting the
solution to be convergent when the right hand side $f$ is smooth enough.

$N=20$

(a) *

$N=40$

(b) *

Figure 4. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/3,1/5)$ on a coarser and a finer Crossed
mesh

$N=20$

(a) *

$N=40$

(b) *

Figure 5. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/3,1/5)$ on a coarser and a finer Right
mesh

$N=20$

(a) *

$N=40$

(b) *

Figure 6. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/3,1/5)$ on a coarser and a finer Non
Structured mesh

We first consider the case of a singular $f$, equal to an approximation of the
Dirac delta function centered at $(1/3,1/5)$. The results for the Crossed,
Right, and Non Structured mesh sequences are reported in Figures 4, 5, and 6,
respectively. It is clear that the solution is affected by spurious components
that are pretty much related to oscillating eigenfunctions like the one
reported in Figure 2.

$N=20$

(a) *

$N=40$

(b) *

Figure 7. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/2,1/2)$ on a coarser and a finer Crossed
mesh

$N=20$

(a) *

$N=40$

(b) *

Figure 8. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/2,1/2)$ on a coarser and a finer Right
mesh

$N=20$

(a) *

$N=40$

(b) *

Figure 9. Solution corresponding to $f$ equal to the approximation of the
Dirac delta function centered at $(1/2,1/2)$ on a coarser and a finer Non
Structured mesh

The situation is even more apparent if the Dirac delta function is centered at
the point $(1/2,1/2)$, that is at the center of our domain $\Omega$. The
corresponding solutions are reported in Figures 7, 8, and 9.

Crossed mesh

(a) *

Right mesh

(b) *

Non Structured mesh

(c) *

Figure 10. Solution corresponding to the smooth right hand side
$f(x,y)=x-3y+\sin(x)$

Figure 10 shows the solution in a case when $f$ is smooth. We take
$f(x,y)=x-3y+\sin(x)$ and we can observe from the figures that the discrete
solutions do not present any spurious oscillation.

We conclude this section by reporting the rates of convergence in some cases
when the exact solutions are known.

Table 2 shows the rates of convergence with respect of the meshsize in the
case when the solution is one of the Laplace eigenfunctions, namely
$u(x,y)=\sin(\pi x)\sin(2\pi y)$

dof’s | $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{0}$ | $\|u-u_{h}\|_{0}$ | $\|\operatorname{\nabla}(u-u_{h})\|_{0}$
---|---|---|---
3281 | 7.178e-02 | 4.415e-03 | 7.324e-02
12961 | 3.389e-02 (1.09) | 1.215e-03 (1.88) | 4.081e-02 (0.85)
51521 | 1.695e-02 (1.00) | 3.029e-04 (2.01) | 2.036e-02 (1.01)
205441 | 8.341e-03 (1.02) | 7.697e-05 (1.98) | 1.023e-02 (1.00)
820481 | 4.173e-03 (1.00) | 1.920e-05 (2.01) | 5.094e-03 (1.01)
3279361 | 2.083e-03 (1.00) | 4.751e-06 (2.02) | 2.513e-03 (1.02)
13112321 | 1.051e-03 (0.99) | 1.182e-06 (2.01) | 1.263e-03 (0.99)
Table 2. Rate of convergence with respect to $h$: smooth solution on the
Crossed mesh in a square

The last example is related to a case where the solution is singular due to a
reentrant corner on the L-shaped domain
$\Omega=(-1,1)^{2}\setminus(0,1)\times(-1,0)$. We take as exact solution the
harmonic function $u(\rho,\theta)=\rho^{2/3}\sin((2/3)\theta)$, where
$(\rho,\theta)$ are the polar coordinates centered at the origin. We consider
$f=0$ and Dirichlet boundary conditions given by the exact solution. In
particular, the solution is vanishing along the two sides of the reentrant
corner meeting at the origin. It is well known that $u\in
H^{5/3-\epsilon}(\Omega)$ for $\epsilon>0$ but $u\not\in H^{5/3}(\Omega)$. The
rates of convergence are shown in Table 3 and are the expected ones:
approximately order $2/3$ for the energy norm and the suboptimal order $4/3$
for the error $\|u-u_{h}\|$ in $L^{2}(\Omega)$.

dof’s | $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{0}$ | $\|u-u_{h}\|_{0}$ | $\|\operatorname{\nabla}(u-u_{h})\|_{0}$
---|---|---|---
1155 | 7.181e-02 | 3.970e-03 | 1.136e-01
4485 | 4.864e-02 (0.57) | 1.302e-03 (1.64) | 7.063e-02 (0.70)
17673 | 3.096e-02 (0.66) | 4.392e-04 (1.59) | 4.457e-02 (0.67)
70161 | 1.967e-02 (0.66) | 1.518e-04 (1.54) | 2.810e-02 (0.67)
279585 | 1.246e-02 (0.66) | 5.395e-05 (1.50) | 1.779e-02 (0.66)
1116225 | 7.875e-03 (0.66) | 1.963e-05 (1.46) | 1.125e-02 (0.66)
4460673 | 4.972e-03 (0.66) | 7.304e-06 (1.43) | 7.102e-03 (0.66)
Table 3. Rate of convergence with respect to $h$: singular solution on a Non
Structured mesh in the L-shaped domain

Figure 11 shows a typical mesh for this computation and Figure 12 reports the
solution $(\boldsymbol{\sigma}_{h},u_{h})$ computed on the same mesh.

Figure 11. Non Structured mesh of the L-shaped domain

Figure 12. Singular solution approximated on the Non Structured mesh in the
L-shaped domain

##  Appendix A The asymptotic behavior of the inf-sup constant

In this section we discuss theoretically the asymptotic behavior of the inf-
sup constant that has been studied numerically in Section 4.

We start by some remarks related to the continuous inf-sup condition: given
$u\in H^{1}_{0}(\Omega)$ there exists $\boldsymbol{\sigma}\in
L^{2}(\Omega)^{2}$ such that

(26)
$\displaystyle(\boldsymbol{\sigma},\operatorname{\nabla}u)=\|\operatorname{\nabla}u\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\|\boldsymbol{\sigma}\|_{L^{2}(\Omega)}\leq
C\|\operatorname{\nabla}u\|_{L^{2}(\Omega)}.$

This can be easily achieved by defining
$\boldsymbol{\sigma}=\operatorname{\nabla}u$ and the constant $C$ is equal to
one in this case. The same consideration could lead to a uniform discrete inf-
sup condition if we had the inclusion
$\operatorname{\nabla}U_{h}\subset\Sigma_{h}$. In our case, however, the
inclusion is not satisfied because the space $\Sigma_{h}$ is in
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$. The natural question is then
if it is possible, at the continuous level, to find
$\boldsymbol{\sigma}\in\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$
satisfying (26). This is certainly true if $\operatorname{\nabla}u$ belongs to
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ which is the case, for
instance, when $u$ is the solution of a Poisson problem $-\Delta u=g$ for some
$g$ in $L^{2}(\Omega)$. We state this result in the following proposition.

######  Proposition 22.

Let $u\in H^{1}_{0}(\Omega)$ be the solution of the problem

$(\operatorname{\nabla}u,\operatorname{\nabla}v)=(g,v)\quad\forall v\in
H^{1}_{0}(\Omega)$

for some $g\in L^{2}(\Omega)$. Then
$\boldsymbol{\sigma}=\operatorname{\nabla}u$ belongs to
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ and satisfies (26) with
$C=1$.

We now try to mimic this proposition at the discrete level. Since the
divergence of $\Sigma_{h}$ is piecewise constant, it is natural to start by
considering $u_{h}\in U_{h}$ that solves

(27)
$(\operatorname{\nabla}u_{h},\operatorname{\nabla}v)=(g_{h},v)\quad\forall
v\in U_{h}$

for some piecewise constant right hand side $g_{h}$. If we define
$\boldsymbol{\sigma}_{h}=\operatorname{\nabla}u_{h}$ we have that
$\boldsymbol{\sigma}_{h}$ is not in
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ since
$\operatorname{\nabla}u_{h}$ has only the tangential component continuous
across the elements while we would need the normal one. We are then led to
some sort of equilibration strategy, in the spirit of [11, 8, 7, 13, 5]. The
optimal way to achieve the equilibration is to solve a global problem and to
define $\boldsymbol{\sigma}_{h}\in\Sigma_{h}$ as the solution of the mixed
problem: find $\boldsymbol{\sigma}_{h}\in\Sigma_{h}$ and
$p_{h}\in\mathcal{P}_{0}$ such that

(28) $\left\\{\begin{aligned}
&(\boldsymbol{\sigma}_{h},\boldsymbol{\tau})+(\operatorname{\mathrm{div}}\boldsymbol{\tau},p_{h})=0&&\forall\boldsymbol{\tau}\in\Sigma_{h}\\\
&(\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h},q)=-(g_{h},q)&&\forall
q\in\mathcal{P}_{0}\end{aligned}\right.$

Then we have $\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h}=-g_{h}$ which
implies

$(\boldsymbol{\sigma},\operatorname{\nabla}u_{h})=(g_{h},u_{h})=\|\operatorname{\nabla}u_{h}\|_{L^{2}(\Omega)}^{2}$

Hence, we get a bound for the discrete inf-sup condition if we can estimate
the constant $C(h)$ for which it holds

$\|\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C(h)\|\operatorname{\nabla}u_{h}\|_{L^{2}(\Omega)}.$

Clearly, from (28) we have $\|\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C\|g_{h}\|_{L^{2}(\Omega)}$ so that we need a bound of
$\|g_{h}\|_{L^{2}(\Omega)}$ in terms of
$\|\operatorname{\nabla}u_{h}\|_{L^{2}(\Omega)}$. Since $u_{h}$ solves (27),
this bound is related to an inf-sup condition between the spaces $U_{h}$ and
$\mathcal{P}_{0}$. More precisely, the following lemma holds true.

######  Lemma 23.

Let $U_{0}$ be the subspace of $U_{h}$ defined as the solutions of (27) for
some $g_{h}\in\mathcal{P}_{0}$. Let $\zeta(h)>0$ be such that the following
inf-sup condition holds true

(29) $\inf_{w_{h}\in
U_{h}}\sup_{g_{h}\in\mathcal{P}_{0}}\frac{(g_{h},w_{h})}{\|g_{h}\|_{L^{2}(\Omega)}\|w_{h}\|_{L^{2}(\Omega)}}\geq\zeta(h).$

Then the inf-sup condition holds true for the space $U_{0}$

$\inf_{v\in
U_{0}}\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}v)}{\|v\|_{H^{1}(\Omega)}\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}\geq
C\zeta(h).$

The inf-sup constant $\zeta(h)$ in (29) can be investigated theoretically on
special meshes or in the case of a one dimensional domain. It turns out that
its behavior is similar to the one of $\beta_{h}$ in (14). For instance, Table
4 shows the value of the lowest eigenvalue associated with the inf-sup
constant (computed as in Section 4). It turns out that with a uniform mesh
$\zeta(h)$ decays as O(h) while with a non structured mesh the decay is less
pronounced.

h |  Non Structured mesh |  Right mesh
---|---|---
$1/2$ | 0.66666667 | 0.66666667
$1/2^{2}$ | 0.20980372 | 0.33333333
$1/2^{3}$ | 0.09283759 | 0.11409783
$1/2^{4}$ | 0.08061165 | 0.03137791
$1/2^{5}$ | 0.06760737 | 0.00803861
$1/2^{6}$ | 0.0626653 | 0.00202209
Table 4. Lowest eigenvalue associated with the inf-sup constant in (29),
corresponding to the square of $\zeta(h)$

We conclude this appendix by showing the behavior of the _worse_ $w_{h}$ in
(29) as it comes out from the numerical experiments. In one dimension, where
the behavior of $\zeta(h)$ is $O(h^{2})$ in agreement with [1], we get the
highly oscillating function plotted in Figure 13.

Figure 13. The function $w_{h}$ corresponding to the first singular value of
(29) in one dimension

The analogous functions in the two dimensional Right and Non Structured meshes
are plotted in Figures 14 and 15, respectively.

Figure 14. The function $w_{h}$ corresponding to the first singular value of
(29) on the Right mesh Figure 15. The function $w_{h}$ corresponding to the
first singular value of (29) on the Non Structured mesh

##  Appendix B Some connections with flux equilibration

We have seen that the inf-sup constant $\beta_{h}$ in general is not bounded
below independently of $h$. We have also seen in several parts of this paper
some analogies between the inf-sup condition and the well known flux
equilibration strategy used, for instance, in the a posteriori analysis of
standard finite elements (see, for instance, [11, 8, 7, 13, 5]). In this
appendix we explore this connections in more detail.

To this aim, we make use of the following natural modification of the usual
Fortin trick [6].

######  Proposition 24.

Let $U_{0}$ be a subspace of $U_{h}$ and assume that there exists a linear
projection $\Pi:\operatorname{\nabla}(U_{0})\to\Sigma_{h}$ such that for any
$\boldsymbol{\tau}\in\operatorname{\nabla}(U_{0})$

(30)
$\displaystyle(\boldsymbol{\tau}-\Pi\boldsymbol{\tau},\operatorname{\nabla}v)=0$
$\displaystyle\forall v\in U_{0}$
$\displaystyle\|\Pi\boldsymbol{\tau}\|_{L^{2}(\Omega)}\leq
C_{\Pi}\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}.$

Then the following inf-sup condition holds true

(31) $\inf_{v\in
U_{0}}\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}v)}{\||\boldsymbol{\tau}\|_{L^{2}(\Omega)}\|v\|_{H^{1}(\Omega)}}\geq\beta$

with $\beta=1/(C_{P}C_{\Pi})$, being $C_{P}$ the Poincaré constant.

###### Proof.

Take $v\in U_{0}$, then we have

$\displaystyle\|v\|_{H^{1}(\Omega)}$ $\displaystyle\leq
C_{P}\|\operatorname{\nabla}v\|_{L^{2}(\Omega)}=C_{P}\frac{(\operatorname{\nabla}v,\operatorname{\nabla}v)}{\|\operatorname{\nabla}v\|_{L^{2}(\Omega)}}$
$\displaystyle\leq
C_{P}\sup_{\boldsymbol{\tau}\in\operatorname{\nabla}(U_{0})}\frac{(\boldsymbol{\tau},\operatorname{\nabla}v)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}$
$\displaystyle=C_{P}\sup_{\boldsymbol{\tau}\in\operatorname{\nabla}(U_{0})}\frac{(\Pi\boldsymbol{\tau},\operatorname{\nabla}v)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}}\leq{C_{P}}{C_{\Pi}}\sup_{\boldsymbol{\tau}\in\operatorname{\nabla}(U_{0})}\frac{(\Pi\boldsymbol{\tau},\operatorname{\nabla}v)}{\|\Pi\boldsymbol{\tau}\|_{L^{2}(\Omega)}}$
$\displaystyle\leq{C_{P}}{C_{\Pi}}\sup_{\boldsymbol{\tau}\in\Sigma_{h}}\frac{(\boldsymbol{\tau},\operatorname{\nabla}v)}{\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}},$

which implies the inf-sup condition (31) with $\beta=1/(C_{P}C_{\Pi})$.

∎

Looking carefully at the properties of the Fortin projector in (30), we can
see that the construction of $\Pi$ consists in finding a mapping from
gradients of continuous piecewise polynomials into the lowest order
Raviart–Thomas space in $\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$. This
problem has been widely studied, for instance, in the framework of flux
equilibration of finite element spaces.

In this appendix we try to identify a suitable subspace $W_{1}$ of $U_{h}$
that could be used as $U_{0}$ in Proposition 24. We start with the following
heuristic reasoning. Let us take a function $f\in L^{2}(\Omega)$ and for each
$h$ consider the solution $u_{h}\in U_{h}$ of the following standard Galerkin
problem

(32)
$(\operatorname{\nabla}u_{h},\operatorname{\nabla}v_{h})=(\Pi_{0}f,v_{h})\qquad\forall
v_{h}\in U_{h}.$

Then the most natural way to find an equilibrated
$\boldsymbol{\sigma}_{h}\in\Sigma_{h}$ is to solve the following global mixed
problem: find
$(\boldsymbol{\sigma}_{h},p_{h})\in\Sigma_{h}\times\mathcal{P}_{0}$ such that

$\left\\{\begin{aligned}
&(\boldsymbol{\sigma}_{h},\boldsymbol{\tau}_{h})+(\operatorname{\mathrm{div}}\boldsymbol{\tau}_{h},p_{h})=0&&\forall\boldsymbol{\tau}_{h}\in\Sigma_{h}\\\
&(\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h},q_{h})=-(f,q_{h})&&\forall
q_{h}\in\mathcal{P}_{0}.\end{aligned}\right.$

By considering $\boldsymbol{\tau}=\operatorname{\nabla}u_{h}$ and
$\Pi\boldsymbol{\tau}=\boldsymbol{\sigma}_{h}$ and observing that
$\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h}=-\Pi_{0}f$, we then have

$(\boldsymbol{\tau}-\Pi\boldsymbol{\tau},\operatorname{\nabla}v_{h})=(\Pi_{0}f,v)+(\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h},v)=0\qquad\forall
v\in U_{h},$

that is, the first Fortin condition in (30) is satisfied. Moreover, we can
bound $\boldsymbol{\sigma}_{h}$ by triangular inequality as follows

$\|\boldsymbol{\sigma}_{h}\|_{0}\leq\|\boldsymbol{\sigma}_{h}-\operatorname{\nabla}u_{h}\|_{0}+\|\operatorname{\nabla}u_{h}\|_{0}.$

The first term on the right hand side is going to zero and can be bounded by

$\|\boldsymbol{\sigma}_{h}-\operatorname{\nabla}u_{h}\|_{0}\leq
h^{s}\|\Pi_{0}f\|_{0}$

if the solution of the continuous problem corresponding to (32) has regularity
$H^{1+s}(\Omega)$. In order to estimate $\|\boldsymbol{\sigma}_{h}\|_{0}$
uniformly in terms of $\|\operatorname{\nabla}u_{h}\|_{0}$ we then need to
bound $\|\Pi_{0}f\|_{0}$ uniformly by
$h^{-s}\|\operatorname{\nabla}u_{h}\|_{0}$. Unfortunately this cannot be done
in general.

The equilibration technique is used in the next theorem to show how the Fortin
operator restricted to a suitable subspace of $U_{h}$ behaves asymptotically
in $h$.

######  Theorem 25.

Let $W_{1}$ be the following subspace of $U_{h}$:

$W_{1}=\\{w\in U_{h}\,|\,\exists
g_{0}\in\mathcal{P}_{0}:(\operatorname{\nabla}w,\operatorname{\nabla}v)=(g_{0},v)\
\forall v\in U_{h}\\},$

where $\mathcal{P}_{0}$ is the space of piecewise constant functions. Let
$\zeta(h)$ be the inf-sup constant introduced in (29) and $\rho(h)$ a function
of $h$ so that the inverse inequality
$\|\operatorname{\nabla}v\|_{0}\leq\rho(h)\|v\|_{0}$ holds true. Then there
exists a Fortin operator $\Pi$ as in Proposition 24 with $U_{0}=W_{1}$ and
$C_{\Pi}\leq Ch\rho(h)/\zeta(h)$.

###### Proof.

We denote by $\mathcal{T}_{h}$ our triangulation of $\Omega$ and by
$\mathcal{E}_{h}$ the skeleton of the edges. Given $u_{h}\in W_{1}$, we
consider $\boldsymbol{\tau}=\operatorname{\nabla}u_{h}$ so that the flux
reconstruction procedure provided in [7] and described with explicit formulas
in [3] will give a function
$\Pi\boldsymbol{\tau}=\boldsymbol{\sigma}_{h}\in\Sigma_{h}$. The construction
is local and is performed as

$\boldsymbol{\sigma}_{h}=\nabla
u_{h}+\sum\limits_{z\in\mathcal{V}}\boldsymbol{\sigma}_{z}$

where $\boldsymbol{\sigma}_{z}$ is a _discontinuous_ $RT_{0}$ function
($DRT_{0}$) supported on the patch

$\omega_{z}=\\{T\in\mathcal{T}:z\text{ is a vertex of }T\\}.$

As the space $DRT_{0}$ consists of edge basis functions, we will also consider
the set

$\mathcal{E}_{z}=\\{E\in\mathcal{E}:z\text{ is a vertex of }E\\}.$

Let $g_{0}\in\mathcal{P}_{0}$ be such that
$(\operatorname{\nabla}u_{h},\operatorname{\nabla}v)=(g_{0},v)$ $\forall v\in
U_{h}$ as in the definition of $W_{1}$. Following the standard procedure of
flux reconstruction we have that $\boldsymbol{\sigma}_{h}$ belongs to
$\Sigma_{h}$ and that the Fortin property
$(\operatorname{\nabla}u_{h}-\boldsymbol{\sigma}_{h},\operatorname{\nabla}v)=0$
($\forall v\in W_{1}$) is satisfied if the $\boldsymbol{\sigma}_{z}$ are
chosen such that

(33)
$\displaystyle\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{z}|_{T}=-\frac{1}{|T|}(g_{0},\phi_{z})_{T}$
$\displaystyle\forall T\in\omega_{z}$
$\displaystyle\llbracket\boldsymbol{\sigma}_{z}\cdot\mathbf{n}\rrbracket_{E}=-\frac{1}{2}\llbracket\operatorname{\nabla}u_{h}\cdot\mathbf{n}\rrbracket_{E}$
$\displaystyle\forall E\in\mathcal{E}_{z}$
$\displaystyle\boldsymbol{\sigma}_{z}\cdot\mathbf{n}=0$ $\displaystyle\text{on
}\partial\omega_{z}$

where $\phi_{z}$ is the linear nodal basis function corresponding to the node
$z$ and vanishing on the boundary of the patch $\omega_{z}$.

Indeed, the second equation in (33) guarantees that
$\llbracket\boldsymbol{\sigma}_{h}\rrbracket_{E}=0$ for all internal edges
$E$, since each edge belongs to exactly two patches. Hence
$\boldsymbol{\sigma}_{h}$ belongs to
$\mathbf{H}(\operatorname{\mathrm{div}};\Omega)$ and we can evaluate its
divergence element by element. From the first equation in (33) we have

$\left(\operatorname{\mathrm{div}}\sum\limits_{z\in\mathcal{V}}\boldsymbol{\sigma}_{z}\right)\Bigg{|}_{T}=-\frac{1}{|T|}\sum_{i=1}^{3}(g_{0},\phi_{z_{i}})_{T}=-\frac{1}{|T|}(g_{0},1)_{T}=g_{0}|_{T},$

where $z_{i}$ ($i=1,2,3$) are the three vertices of $T$. From
$(\operatorname{\mathrm{div}}\operatorname{\nabla}u_{h})|_{T}=0$ it follows
$\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{h}=-g_{0}$, so that we get
the Fortin property

$(\boldsymbol{\sigma}_{h},\operatorname{\nabla}v)=(g_{0},v)=(\operatorname{\nabla}u_{h},\operatorname{\nabla}v)\qquad\forall
v\in U_{h}.$

The edge basis functions considered for the space $\Sigma_{h}$ are supported
in the two adjacent triangles of the edge $E$. We denote by $T^{-}_{E}$ and
$T^{+}_{E}$ the two triangles adjacent to $E$ and define an edge oriented
basis
$\\{\boldsymbol{\psi}_{E}^{-}\\}_{E\in\mathcal{E}_{h}}\cup\\{\boldsymbol{\psi}_{E}^{+}\\}_{E\in\mathcal{E}_{h},E\not\subset\partial\Omega}$
for $DRT_{0}$, using the basis functions

$\boldsymbol{\psi}_{E}^{-}(\mathbf{x})=\begin{cases}-\frac{1}{2|T^{-}_{E}|}(\mathbf{x}-P^{-}_{E})&\text{on
}T^{-}_{E}\\\ 0&\text{elsewhere}\end{cases}\text{ and
}\boldsymbol{\psi}_{E}^{+}(\mathbf{x})=\begin{cases}\frac{1}{2|T^{+}_{E}|}(\mathbf{x}-P^{+}_{E})&\text{on
}T^{+}_{E}\\\ 0&\text{elsewhere}\end{cases}$

where $P^{-}_{E}$ and $P^{+}_{E}$ are the vertices of $T^{-}_{E}$ and
$T^{+}_{E}$, respectively, not shared by the two triangles. This basis uses a
similar representation of the one presented in [2] although here

$\operatorname{\mathrm{div}}\boldsymbol{\psi}_{E}^{\pm}=\pm\frac{1}{|T^{\pm}_{E}|}$

and thus

$\int_{E}\boldsymbol{\psi}^{\pm}_{E}\cdot\mathbf{n}_{E}=\pm\left(\operatorname{\mathrm{div}}\boldsymbol{\psi}_{E}^{\pm},1\right)_{T^{\pm}_{E}}=1,$

where the normal $\mathbf{n}_{E}$ is pointing from $T^{+}_{E}$ to $T^{-}_{E}$.
It follows

$\boldsymbol{\psi}^{\pm}_{E}\cdot\mathbf{n}_{E}=\frac{1}{|E|}.$

These basis functions allow for the following construction:

$\boldsymbol{\sigma}_{z}=\sum\limits_{E\in\mathcal{E}_{z}}(\tau_{E,z}^{-}\boldsymbol{\psi}_{E}^{-}+\tau_{E,z}^{+}\boldsymbol{\psi}_{E}^{+})$

where the coefficients $\tau_{E,z}^{+}$ and $\tau_{E,z}^{-}$ will be chosen so
that the equilibration conditions (33) hold.

For the computation of the coefficients, we index the triangles in the patch
from $1$ to $n_{z}$ and define $T_{0}:=T_{n_{z}}$ and $T_{n_{z}+1}:=T_{1}$.
Further, we define $E_{i}=T_{i}\cap T_{i+1}$ and $T^{+}_{E_{i}}=T_{i+1}$.
Then,

$\operatorname{\mathrm{div}}\boldsymbol{\sigma}_{z}|_{T_{i}}=\tau_{E_{i-1},z}^{+}\frac{1}{|T_{i}|}-\tau_{E_{i},z}^{-}\frac{1}{|T_{i}|}$

and the first condition in (33) reads

$\tau_{E_{i},z}^{-}=(g_{0},\phi_{z})_{T_{i}}+\tau^{+}_{E_{i-1},z}\quad\forall
i=1,\dots,n_{z}.$

The second condition in (33) implies

$\tau_{E_{i},z}^{+}-\tau_{E_{i},z}^{-}=-\frac{|E_{i}|}{2}\llbracket\operatorname{\nabla}u_{h}\cdot\mathbf{n}\rrbracket_{E_{i}}\quad\forall
i=1,\dots,n_{z}.$

This leads to

$\tau_{E_{i},z}^{+}-\tau^{+}_{E_{i-1},z}=-\frac{|E_{i}|}{2}\llbracket\operatorname{\nabla}u_{h}\cdot\mathbf{n}\rrbracket_{E_{i}}+(g_{0},\phi_{z})_{T_{i}}\quad\forall
i=1,\dots,n_{z}$

and thus to

$\tau_{E_{i},z}^{+}=\tau^{+}_{E_{0},z}+\sum\limits_{j=1}^{i}\left(-\frac{|E_{j}|}{2}\llbracket\operatorname{\nabla}u_{h}\cdot\mathbf{n}\rrbracket_{E_{j}}+(g_{0},\phi_{z})_{T_{j}}\right)\quad\forall
i=1,\dots,n_{z}.$

We can bound the two terms on the right hand side that are involved in the
summation as follows:

$\left|\frac{|E_{j}|}{2}\llbracket\nabla
u_{h}\cdot\mathbf{n}\rrbracket_{E_{j}}\right|\leq
C\|\operatorname{\nabla}u_{h}\|_{0,\omega_{z}}$

and

$|(g_{0},\phi_{z})_{T_{j}}|\leq\|g_{0}\|_{0,T_{j}}\|\phi_{z}\|_{0,T_{j}}\leq
Ch\|g_{0}\|_{0,T_{j}}.$

Choosing $\tau^{+}_{E_{0},z}=0$ we then have

$\displaystyle|\tau_{E_{i},z}^{+}|\leq
C(\|\operatorname{\nabla}u_{h}\|_{0,\omega_{z}}+h\|g_{0}\|_{0,\omega_{z}})$
$\displaystyle|\tau_{E_{i},z}^{-}|\leq
C(\|\operatorname{\nabla}u_{h}\|_{0,\omega_{z}}+h\|g_{0}\|_{0,\omega_{z}})$
$\displaystyle\forall i=1,\dots,n_{z}.$

By a scaling argument or by using a suitable quadrature rule we see that

$\|\boldsymbol{\psi}_{E}^{\pm}\|_{0}^{2}=\frac{1}{4|T_{E}^{\pm}|^{2}}\|\mathbf{x}-P_{E}^{\pm}\|_{0}^{2}\leq
Ch^{2}/|T_{E}^{\pm}|\leq C$

so that it holds $\|\boldsymbol{\sigma}_{z}\|_{0,\omega_{z}}\leq
C(\|\operatorname{\nabla}u_{h}\|_{0,\omega_{z}}+h\|g_{0}\|_{0,\omega_{z}})$.
By putting all the patches together and considering that the intersections
between patches contain a bounded number of elements and that each element
belongs to a bounded number of patches we get

$\|\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C(\|\operatorname{\nabla}u_{h}\|_{L^{2}(\Omega)}+h\|g_{0}\|_{L^{2}(\Omega)}).$

It remains to estimate $g_{h}$, which can be done by considering the inf-sup
constant discussed in Appendix A concerning the $P_{0}-P_{1}$ element (29). By
the definition of $g_{0}$ we have

$(g_{0},v_{h})=(\operatorname{\nabla}u_{h},\operatorname{\nabla}v_{h})\quad\forall
v_{h}\in U_{h}$

and hence

$\displaystyle\zeta(h)\|g_{0}\|_{0}$ $\displaystyle\leq\sup_{v_{h}\in
U_{h}}\frac{(g_{0},v_{h})}{\|v\|_{0}}=\sup_{v_{h}\in
U_{h}}\frac{(\operatorname{\nabla}u_{h},\operatorname{\nabla}v_{h})}{\|v\|_{0}}\leq\sup_{v_{h}\in
U_{h}}\frac{\|\operatorname{\nabla}u_{h}\|_{0}\|\operatorname{\nabla}v_{h}\|_{0}}{\|v_{h}\|_{0}}$
$\displaystyle\leq\rho(h)\|\operatorname{\nabla}u_{h}\|_{0}.$

Finally, we arrive at the final estimate

$\|\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}\leq
C\|\operatorname{\nabla}u_{h}\|_{L^{2}(\Omega)}(1+h\rho(h)/\zeta(h)),$

which leads to the bound $C_{\Pi}\leq Ch\rho(h)/\zeta(h)$. ∎

Although the above considerations do not provide a rigorous proof that the
inf-sup constant is vanishing with $h$, they give a clear indication that we
should not expect the constant $\beta_{h}$ to be uniformly bounded away from
zero. Indeed, the behavior of the inf-sup constant depends on the chosen mesh
as shown in Section 4: if the mesh is quasiuniform then $h\rho(h)$ is bounded
from above and below so that we have a confirmation that the inf-sup constant
cannot be better than $\zeta(h)$.

## Acknowledgments

The first author gratefully acknowledges support by the Deutsche
Forschungsgemeinschaft in the Priority Program SPP 1748 Reliable simulation
techniques in solid mechanics, Development of non standard discretization
methods, mechanical and mathematical analysis under the project number BE
6511/1-1. The second author is member of the INdAM Research group GNCS and his
research is partially supported by IMATI/CNR and by PRIN/MIUR.

## References

  * [1] I. Babuška and R. Narasimhan, _The Babuška-Brezzi condition and the patch test: an example_ , Comput. Methods Appl. Mech. Engrg. 140 (1997), no. 1-2, 183–199. MR 1423460
  * [2] C. Bahriawati and C. Carstensen, _Three matlab implementations of the lowest-order Raviart–Thomas MFEM with a posteriori error control_ , Computational Methods in Applied Mathematics 5.4 (2005), 333–361.
  * [3] F. Bertrand, B. Kober, M. Moldenhauer, and G. Starke, _Equilibrated stress reconstruction and a posteriori error estimation for linear elasticity_ , Novel Finite Element Technologies for Solids and Structures. CISM International Centre for Mechanical Sciences (Courses and Lectures), vol 597. Springer, Cham (de Mattos Pimenta P. Schröder J., ed.), 2020.
  * [4] Daniele Boffi and Fleurianne Bertrand, _A counterexample for the inf-sup stability of the $RT^{0}-P^{1}\subset L^{2}\times H^{1}_{0}$ finite element combination for the mixed Poisson equation_, Proceedings in Applied Mathematics and Mechanics (2019), to appear.
  * [5] by same author, _The Prager–Synge theorem in reconstruction based a posteriori error estimation_ , 75 Years of Mathematics of Computation (Susanne C. Brenner, Igor Shparlinski, Chi-Wang Shu, and Daniel B. Szyld, eds.), Contemporary Mathematics, vol. 754, American Mathematical Society, 2020, to appear.
  * [6] Daniele Boffi, Franco Brezzi, and Michel Fortin, _Mixed finite element methods and applications_ , Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013. MR 3097958
  * [7] D. Braess and J. Schöberl, _Equilibrated residual error estimator for edge elements_ , Math. Comp. 77 (2008), no. 262, 651–672. MR 2373174
  * [8] Dietrich Braess, _Finite elements_ , third ed., Cambridge University Press, Cambridge, 2007, Theory, fast solvers, and applications in elasticity theory, Translated from the German by Larry L. Schumaker. MR 2322235
  * [9] D. Chapelle and K.-J. Bathe, _The inf-sup test_ , Comput. & Structures 47 (1993), no. 4-5, 537–545. MR 1224095
  * [10] D. S. Malkus, _Eigenproblems associated with the discrete LBB condition for incompressible finite elements_ , Internat. J. Engrg. Sci. 19 (1981), no. 10, 1299–1310. MR 660563
  * [11] W. Prager and J. L. Synge, _Approximations in elasticity based on the concept of function space_ , Quart. Appl. Math. 5 (1947), 241–269. MR 25902
  * [12] N. Viebahn, K Steeger, and J. Schröder, _A simple and efficient Hellinger–Reissner type mixed finite element for nearly incompressible elasticity_ , Computer Methods in Applied Mechanics and Engineering 340 (2018), 278–295.
  * [13] Martin Vohralík, _Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods_ , Math. Comp. 79 (2010), no. 272, 2001–2032. MR 2684353