\begin{document}

\title{\sc Functional A Posteriori Error Equalities\\ for Conforming Mixed Approximations\\ of Elliptic Problems}

\begin{dedi}{0}
Dedicated to Sergey Igorevich Repin on the occasion of his $60^{\textrm{th}}$ birthday
\end{dedi}

\begin{abstr}{12}
In this paper we show how to find the {\it exact error}
(not just an estimate of the error)
of a conforming mixed approximation by using the functional type
a posteriori error estimates in the spirit of Repin \cite{repinbookone}. 
The error is measured in a mixed norm which takes into account both
the primal and dual variables. We derive this result 
for all elliptic partial differential equations of the class
$$\A^{*}\A\,x+x=f,$$
where $\A$ is a linear, densely defined and closed (usually a differential) operator 
and $\A^{*}$ its adjoint. 
We first derive a special version of our main result by using a simplified reaction-diffusion problem 
to demonstrate the strong connection to the classical functional a posteriori error estimates 
of Repin \cite{repinbookone}. After this we derive the main result in an abstract setting. 
Our main result states that in order to obtain the {\it exact global error} value 
of a conforming mixed approximation with primal variable $x$ 
and dual variable $y$, i.e.,
$$\A^{*} y+x=f,\quad\A x=y,$$ 
one only needs the problem data and the mixed approximation $(\tilde{x},\tilde{y})\in D(\A)\times D(\A^{*})$ 
of the exact solution $(x,y)\in D(\A)\times\big(D(\A^{*})c_{\A}p R(\A)\big)$, i.e.,
the {\it equality}
$$\norm{x-\tilde{x}}^2
+\norm{\A(x-\tilde{x})}^2
+\norm{y-\tilde{y}}^2
+\norm{\A^{*}(y-\tilde{y})}^2
=\norm{f-\tilde{x}-\A^{*}\tilde{y}}^2
+\norm{\tilde{y}-\A\tilde{x}}^2$$
holds. There is no need for calculating any auxiliary data. 
The calculation of the exact error consists of simply calculating two (usually integral) quantities 
where all the quantities are known after the approximate solution has been obtained
by any conforming method guaranteeing $(\tilde{x},\tilde{y})\in D(\A)\times D(\A^{*})$. 
We also show some numerical computations to confirm the results.
\end{abstr}

\begin{keyw}{0}
functional a posteriori error estimate, error equality, 
elliptic boundary value problem, mixed formulation, combined norm
\end{keyw}

\begin{amscl}{0}
65N15
\end{amscl}

\tableofcontents

\section{Introduction}

The results presented in this paper are based on the conception of functional type a posteriori error estimates. 
These type estimates are valid for any conforming approximation and contain only global constants. 
We note that estimates for nonconforming approximations are known as well but will not be discussed in this paper.
In the case of the class of PDEs studied in this paper, the estimates do not contain even global constants. 
For a detailed exposition of the theory see the books \cite{repinbookone} 
by Repin and \cite{NeittaanmakiRepin2004} by Repin and Neittaanm\"aki
or for a more computational point of view \cite{MaliRepinNeittaanmaki2014}
by Mali, Repin, and Neittaanm\"aki.

We will measure the error of our approximations in a combined norm, which includes the error of both, 
the primal and the dual variable. This is especially useful for mixed methods 
where one calculates an approximation for both the primal and dual variables, see e.g.
the book of Brezzi and Fortin \cite{brezzifortinbookone}.

In this paper, we study the linear equation 
$$(\A^{*}\alpha_{2}\A+\alpha_{1})x=f$$
presented in the mixed form
$$\A^{*} y+\alpha_{1} x=f,\quad\alpha_{2}\A=y,$$
where $\alpha_{1},\alpha_{2}$ are linear and self adjoint topological isomorphisms
on two Hilbert spaces $\hilbert_{1}$ and $\hilbert_{2}$
and $\A:D(\A)\subset\hilbert_{1}\to\hilbert_{2}$ is a linear, densely defined and closed operator
with adjoint operator $\A^{*}:D(\A^{*})\subset\hilbert_{2}\to\hilbert_{1}$.
Our main result is Theorem \rgen{}{}ef{thm:Gmain} and it shortly reads 
as the {\it functional a posterior error equality}
\begin{align}
\label{erroreq}
\begin{split}
&\qquad\normhioao{x-\tilde{x}}^2
+\normhitat{\A(x-\tilde{x})}^2
+\normhitatmo{y-\tilde{y}}^2
+\normhioaomo{\A^{*}(y-\tilde{y})}^2\\
&=\normhioaomo{f-\alpha_{1}\tilde{x}-\A^{*}\tilde{y}}^2
+\normhitatmo{\tilde{y}-\alpha_{2}\A\tilde{x}}^2
\end{split}
\end{align}
being valid for any conforming mixed approximation $(\tilde{x},\tilde{y})\in D(\A)\times D(\A^{*})$ 
of the exact solution $(x,y)\in D(\A)\times D(\A^{*})$.

Functional a posteriori error estimates for combined norms 
were first exposed in the paper \cite{repinsautersmolianskiaposttwosideell},
where the authors present two-sided estimates bounding the error
by the same quantity from below and from above aside from multiplicative constants. 
Unlike in other estimates, these constants are $1$ and $\sqrt{3}$.
In \cite{repinsautersmolianskiaposttwosideell} the authors studied problems of the type 
\begin{align}
\label{probone}
\A^{*}\alpha\A x=f ,
\end{align}
i.e., the case $\alpha=\alpha_{2}$, $\alpha_{1}=0$.

The paper is organized as follows. 
In Section \rgen{}{}ef{sec:M} we prove our main results for a simple model problem 
and show the strong connection to the classical functional a posteriori error estimates. 
In Section \rgen{}{}ef{sec:G} we derive our main results in an abstract Hilbert space setting and
in Section \rgen{}{}ef{sec:A} we show applications of the general results to several classical problems. 
Section \rgen{}{}ef{sec:inhomobc} is devoted to inhomogeneous boundary conditions and
finally in Section \rgen{}{}ef{sec:N} we present some numerical experiments to confirm our theoretical results.

\section{Results for a Model Problem} \label{sec:M}

Let $\Omega\subset\mathbb{R}^d$, $d\geq1$, be open and without loss of generality connected,
so let $\Omega$ be a domain with boundary $\Gamma:=\partial\Omega$. 
We emphasize that $\Omega$ may be bounded or unbounded, like an exterior domain, or non of both.
Moreover, $\Gamma$ does not need to have any smoothness.
We denote by $\scplt{\,\cdot\,}{\,\cdot\,}$ and $\normlt{\,\cdot\,}$ the inner product and the norm 
in $\lgen{2}{}$ for scalar-, vector- and matrix-valued functions. 
Throughout the paper we will not indicate the dependence on $\Omega$
in our notations of the functional spaces.
Moreover, we define the usual Sobolev spaces
$$\hgen{1}{}:=\set{\varphi\in\lgen{2}{}}{\nabla\varphi\in\lgen{2}{}},\quad\dgen{}{}:=\set{\partialsi\in\lgen{2}{}}{\opdiv\partialsi\in\lgen{2}{}}$$
and as the closure of smooth and compactly supported test functions\footnote{The 
spaces $\cgen{\infty}{\ga}$ and $\hgen{1}{}ga$ are often denoted by $\cgen{\infty}{\circ}$ and $\hgen{1}{}c$.}
$$\hgen{1}{}ga:=\overline{\cgen{\infty}{\ga}}^{\hgen{1}{}}.$$
These are Hilbert spaces equipped with the respective graph norms denoted by
$\normho{\,\cdot\,}$, $\normd{\,\cdot\,}$.

Our simple model reaction-diffusion problem reads as follows: 
Find the potential $u\in\hgen{1}{}ga$, i.e., the primal variable, such that
\begin{align}
-\Delta u+u=-\opdiv\nabla u+u=f,\label{eq:Ipde}
\end{align}
where $f \in \lgen{2}{}$ is the source term. The variational formulation
of this problem consists of finding $u\in\hgen{1}{}ga$ such that
\begin{align} 
\label{eq:Igen}
\forall\,\varphi\in\hgen{1}{}ga\qquad
\scplt{\nabla u}{\nabla\varphi}+\scplt{u}{\varphi}=\scplt{f}{\varphi}.
\end{align}
The natural energy norm for this problem is $\normho{\,\cdot\,}$.
Of course, by the Lax-Milgram lemma or Riesz' representation theorem 
\eqref{eq:Igen} has a unique solution $u\in\hgen{1}{}ga$ satisfying 
$$\normho{u}\leq\normlt{f}.$$
Often, a variable of interest is also the flux, i.e., the dual variable,
$$p:=\nabla u\in\dgen{}{},$$
leading to the mixed formulation
$$-\opdiv p+u=f,\quad\nabla u=p.$$ 
We note that indeed by \eqref{eq:Igen} the flux $p$ belongs to $\dgen{}{}$ and $\opdiv p=u-f$ holds.
Let us further emphasize that even
$$p\in\dgen{}{}c_{\A}p\,\nabla\hgen{1}{}ga$$
holds, this is, $p$ is also irrotational, has got vanishing tangential trace 
and is $\lgen{2}{}$-perpendicular to the so-called Dirichlet fields.

We will understand a pair $(\tilde{u},\partialt)\in\hgen{1}{}ga\times\dgen{}{}$ 
without further requirements as an approximation 
of the exact solution pair $(u,p)\in\hgen{1}{}ga\times\dgen{}{}$. 
For the convenience of the reader, 
we first present the classical functional error upper bounds,
frequently called \emph{error majorants}, for the approximations of $u$ and $p$. 

\begin{theo} 
\label{thm:Iest}
For any approximation $\tilde{u}\in\hgen{1}{}ga$ of the exact potential $u$ 
\begin{align}
\normho{u-\tilde{u}}^2=\min_{\partialsi\in\dgen{}{}}\mathcal{M}na(\tilde{u},\partialsi)=\mathcal{M}na(\tilde{u},p),\label{eq:Iest}
\end{align}
holds, where
\begin{align} 
\label{eq:Imaj}
\mathcal{M}na(\tilde{u},\partialsi):=\normlt{f-\tilde{u}+\opdiv\partialsi}^2+\normlt{\partialsi-\nabla\tilde{u}}^2.
\end{align}
\end{theo}

\begin{proof}
To derive the upper bound, we subtract $\scplt{\nabla\tilde{u}}{\nabla\varphi}+\scplt{\tilde{u}}{\varphi}$ 
from both sides of the generalized form \eqref{eq:Igen}, and obtain for all $\varphi\in\hgen{1}{}ga$
\begin{align} 
\label{eq:Imaj_1}
\scplt{\nabla(u-\tilde{u})}{\nabla \varphi}+\scplt{u-\tilde{u}}{\varphi}
=\scplt{f-\tilde{u}}{\varphi}-\scplt{\nabla\tilde{u}}{\nabla\varphi}.
\end{align}
For an arbitrary function $\partialsi\in\dgen{}{}$ and any $\varphi\in\hgen{1}{}ga$  
we have $\scplt{\opdiv\partialsi}{\varphi}+\scplt{\partialsi}{\nabla\varphi}=0$.
By adding this to the right hand side of \eqref{eq:Imaj_1} it becomes
\begin{align} 
\label{eq:Imaj_2}
\scplt{\nabla(u-\tilde{u})}{\nabla \varphi}+\scplt{u-\tilde{u}}{\varphi}
&=\scplt{f-\tilde{u}+\opdiv\partialsi}{\varphi}+\scplt{\partialsi-\nabla\tilde{u}}{\nabla\varphi}\nonumber\\
& \leq\normlt{f-\tilde{u}+\opdiv\partialsi}\normlt{\varphi}
+\normlt{\partialsi-\nabla\tilde{u}}\normlt{\nabla\varphi}\\
&	\leq\mathcal{M}na(\tilde{u},\partialsi)^{\nicefrac{1}{2}}\normho{\varphi}.\nonumber
\end{align}
By choosing $\varphi:=u-\tilde{u}\in\hgen{1}{}ga$  we obtain 
$\normho{u-\tilde{u}}^2\leq\mathcal{M}na(\tilde{u},\partialsi)$.
Since $p\in\dgen{}{}$, we see that $\mathcal{M}na(\tilde{u},p)=\normho{u-\tilde{u}}^2$.
\end{proof}

As the majorant $\mathcal{M}na$ is sharp,
it immediately provides a technique to obtain approximations for the exact flux $p$. 
Minimizing $M_{\nabla}(\partialsi):=\mathcal{M}na(\tilde{u},\partialsi)$ with respect to $\partialsi$
yields by differentiation for all $\partialsi\in\dgen{}{}$
\begin{align*}
0\overset{!}{=}M_{\nabla}'(p)\partialsi&=2\scplt{f-\tilde{u}+\opdiv p}{\opdiv\partialsi}+2\scplt{p-\nabla\tilde{u}}{\partialsi}\\
&=2\scplt{f+\opdiv p}{\opdiv\partialsi}+2\scplt{p}{\partialsi}
\end{align*}
since $\scplt{\tilde{u}}{\opdiv\partialsi}=-\scplt{\nabla\tilde{u}}{\partialsi}$ because $\tilde{u}\in\hgen{1}{}ga$.
Hence the following problem occurs: Find $p\in\dgen{}{}$ such that
\begin{align} 
\label{eq:Igenp}
\forall\partialsi\in\dgen{}{}\qquad
\scplt{\opdiv p}{\opdiv\partialsi}+\scplt{p}{\partialsi}=-\scplt{f}{\opdiv\partialsi}.
\end{align}
Note that $\tilde{u}$ is not present here and the natural energy norm for this problem is $\normd{\,\cdot\,}$.
Once again, by the Lax-Milgram lemma \eqref{eq:Igenp} has a unique solution $p\in\dgen{}{}$ satisfying 
$$\normd{p}\leq\normlt{f}.$$
Since $\nabla u\in\dgen{}{}$ solves \eqref{eq:Igenp}, i.e., with \eqref{eq:Ipde}
$$\scplt{\opdiv\nabla u}{\opdiv\partialsi}
=\scplt{u}{\opdiv\partialsi}
-\scplt{f}{\opdiv\partialsi}
=-\scplt{\nabla u}{\partialsi}
-\scplt{f}{\opdiv\partialsi},$$
we get indeed $p=\nabla u$.

\begin{rem}
\mbox{}
\begin{itemize}
\item[\bf(i)]
The variational formulation \eqref{eq:Igenp} for $p$ 
can also be achieved by testing \eqref{eq:Ipde} with $\opdiv\partialsi$ for all $\partialsi\in\dgen{}{}$ since
\begin{align*}
-\scplt{f}{\opdiv\partialsi}
&=\scplt{\opdiv\nabla u}{\opdiv\partialsi}-\scplt{u}{\opdiv\partialsi}\\
&=\scplt{\opdiv\nabla u}{\opdiv\partialsi}+\scplt{\nabla u}{\partialsi}
=\scplt{\opdiv p}{\opdiv\partialsi}+\scplt{p}{\partialsi}.
\end{align*}
\item[\bf(ii)]
By \eqref{eq:Igenp}
$$p\,\bot\,\dgen{}{}z:=\set{v\in\dgen{}{}}{\opdiv v=0}$$
holds. Thus, by the Helmholtz decomposition, i.e., $\lgen{2}{}=\overline{\nabla\hgen{1}{}ga}\oplus\dgen{}{}z$,
we get $p\in\overline{\nabla\hgen{1}{}ga}$. 
Here, $\bot$ and $\oplus$ denote orthogonality and the orthogonal sum in $\lgen{2}{}$.
\item[\bf(iii)]
\eqref{eq:Igenp} is the dual problem to \eqref{eq:Igen} and its strong formulation 
in duality to \eqref{eq:Ipde} is
\begin{align}
\label{nadivpeq}
-\nabla\opdiv p+p=\nabla f
\end{align}
with mixed formulation
$$\nabla v+p=\nabla f,\quad-\opdiv p=v.$$
We note that in general $\opdiv p$ does not belong to $\hgen{1}{}ga$, not even to $\hgen{1}{}$.
On the other hand, by \eqref{eq:Igenp} we see $\opdiv p+f\in\hgen{1}{}ga$ with $\nabla(\opdiv p+f)=p$
and the natural Neumann boundary condition $\opdiv p+f=0$ at $\Gamma$ appears.
Hence $f\in\hgen{1}{}ga$, if and only if $v:=-\opdiv p\in\hgen{1}{}ga$,
and $f\in\hgen{1}{}$, if and only if $v\in\hgen{1}{}$. In both cases \eqref{nadivpeq} holds
and moreover for all $\varphi\in\hgen{1}{}ga$
$$\scplt{\nabla v}{\nabla\varphi}+\scplt{v}{\varphi}
=-\scplt{p}{\nabla\varphi}+\scplt{v}{\varphi}+\scplt{\nabla f}{\nabla\varphi}
=\scplt{\nabla f}{\nabla\varphi},$$
thus $v\in\hgen{1}{}$ solves in the strong sense $-\Delta v+v=-\Delta f$ and $v=f$ at $\Gamma$
if $\Delta f\in\lgen{2}{}$.
\end{itemize}
\end{rem}

\begin{theo} 
\label{thm:Iestp}
For any approximation $\partialt\in\dgen{}{}$ of the exact flux $p$
\begin{align}
\label{eq:Iestp}
\normd{p-\partialt}^2=\min_{\varphi\in\hgen{1}{}ga}\mathcal{M}div(\partialt,\varphi)=\mathcal{M}div(\partialt,u),
\end{align}
holds, where
\begin{align} 
\label{eq:Imajp}
\mathcal{M}div(\partialt,\varphi):=\normlt{f-\varphi+\opdiv\partialt}^2+\normlt{\partialt-\nabla\varphi}^2.
\end{align}
\end{theo}

\begin{proof}
We add $-\scplt{\opdiv\partialt}{\opdiv\partialsi}-\scplt{\partialt}{\partialsi}$ to the both sides 
of the variational formulation \eqref{eq:Igenp} and obtain for all $\partialsi\in\dgen{}{}$
\begin{align}\label{eq:Imajp_1}
\scplt{\opdiv(p-\partialt)}{\opdiv\partialsi}+\scplt{p-\partialt}{\partialsi}=
-\scplt{f+\opdiv\partialt}{\opdiv\partialsi}-\scplt{\partialt}{\partialsi}.
\end{align}
For any $\varphi\in\hgen{1}{}ga$ we have again $\scplt{\nabla\varphi}{\partialsi}+\scplt{\varphi}{\opdiv\partialsi}=0$. 
By adding this to the right hand side of \eqref{eq:Imajp_1} it becomes
\begin{align} 
\label{eq:Imajp_2}
\scplt{\opdiv(p-\partialt)}{\opdiv \partialsi}+\scplt{p-\partialt}{\partialsi}
&=-\scplt{f-\varphi+\opdiv\partialt}{\opdiv\partialsi}-\scplt{\partialt-\nabla\varphi}{\partialsi}\nonumber\\
&\leq\normlt{f-\varphi+\opdiv\partialt}\normlt{\opdiv\partialsi}
+\normlt{\partialt-\nabla\varphi}\normlt{\partialsi}\\
&\leq\mathcal{M}div(\partialt,\varphi)^{\nicefrac{1}{2}}\normd{\partialsi}.\nonumber
\end{align}
Choosing $\partialsi=p-\partialt\in\dgen{}{}$ yields $\normd{p-\partialt}^2\leq\mathcal{M}div(\partialt,\varphi)$.
Finally $\mathcal{M}div(\partialt,u)=\normd{p-\partialt}^2$ follows by $u\in\hgen{1}{}ga$.
\end{proof}

As before, the sharpness of the majorant $\mathcal{M}div$ gives us a technique
to obtain approximations of the potential $u$. 
In fact, now global minimization of $M_{\opdiv}(\varphi):=\mathcal{M}div(\partialt,\varphi)$ 
with respect to $\varphi$ would lead to the variational formulation \eqref{eq:Igen}
for finding $u$, since for all $\varphi\in\hgen{1}{}ga$
\begin{align*}
0\overset{!}{=}M_{\opdiv}'(u)\varphi&=-2\scplt{f-u+\opdiv\partialt}{\varphi}-2\scplt{\partialt-\nabla u}{\nabla\varphi}\\
&=2\scplt{u-f}{\varphi}+2\scplt{\nabla u}{\nabla\varphi}
\end{align*}
since $\scplt{\opdiv\partialt}{\varphi}=-\scplt{\partialt}{\nabla\varphi}$ by $\partialt\in\dgen{}{}$.
 
Finally, we note that the functional a posteriori error majorants $\mathcal{M}na$ and $\mathcal{M}div$ 
contain only the problem data, conforming numerical approximations
and the free functions $\partialsi$ and $\varphi$.

Now, we define the combined norm for the reaction-diffusion problem 
in a canonical way as the sum of the energy norms for the potential and the flux:
\begin{equation*}
\tnorm{(\varphi,\partialsi)}^2 :=
\normho{\varphi}^2+\normd{\partialsi}^2
=\normlt{\varphi}^2+\normlt{\nabla\varphi}^2 
+\normlt{\partialsi}^2+\normlt{\opdiv\partialsi}^2
\end{equation*}

\begin{rem}
We know $\normho{u}\leq\normlt{f}$ and $\normd{p}\leq\normlt{f}$.
It is indeed notable that 
$$\tnorm{(u,p)}=\normlt{f}$$ 
holds, which follows immediately by $f=-\opdiv p+u$ and $p=\nabla u$ since
$$\normlt{f}^2
=\normlt{\opdiv p}^2
+\normlt{u}^2
-2\scplt{\opdiv p}{u}\\
=\normlt{\opdiv p}^2
+\normlt{u}^2
+2\normlt{p}^2
=\tnorm{(u,p)}^2.$$
Hence the solution operator
$$L:\lgen{2}{}\to\hgen{1}{}ga\times\dgen{}{};f\mapsto(u,p)$$
has norm $\norm{L}=1$, i.e., $L$ is an isometry.
\end{rem}

Our main result for this simple reaction-diffusion problem 
basically combines Theorems \rgen{}{}ef{thm:Iest} and \rgen{}{}ef{thm:Iestp}.
However, we outline that the resulting right hand side does not contain $u$ or $p$ anymore
and is even an equality.

\begin{theo} 
\label{thm:Imain}
For any approximation $(\tilde{u},\partialt)\in\hgen{1}{}ga\times\dgen{}{}$ 
of the exact solution $(u,p)$
\begin{align}
\tnorm{(u,p)-(\tilde{u},\partialt)}^2
&=\mathcal{M}mix(\tilde{u},\partialt)
\label{eq:Imain1}
\intertext{and the normalized counterpart}
\frac{\tnorm{(u,p)-(\tilde{u},\partialt)}^2}{\tnorm{(u,p)}^2}
&=\frac{\mathcal{M}mix(\tilde{u},\partialt)}{\normlt{f}^2}
\label{eq:Imain2}
\end{align}
hold, where
\begin{align} 
\label{eq:Imain}
\begin{split}
\mathcal{M}mix(\tilde{u},\partialt)
&:=\mathcal{M}na(\tilde{u},\partialt)=\mathcal{M}div(\partialt,\tilde{u})
=\normlt{f-\tilde{u}+\opdiv\partialt}^2+\normlt{\partialt-\nablabla \tilde{u}}^2.
\end{split}
\end{align}
\end{theo}

The error in the combined norm can thus be exactly computed by quantities we already know: 
the given problem data $f$ and the conforming approximation $(\tilde{u},\partialt)$.

\begin{proof}
Set $\partialsi=\partialt$ in \eqref{eq:Imaj_2} and $\varphi=\tilde{u}$ in \eqref{eq:Imajp_2}.
Then, for any $\varphi\in\hgen{1}{}ga$ and any $\partialsi\in\dgen{}{}$ we have
\begin{align} 
\label{eq:Imain_1}
\scplt{\nabla(u-\tilde{u})}{\nabla \varphi}+\scplt{u-\tilde{u}}{\varphi}
&=\scplt{f-\tilde{u}+\opdiv\partialt}{\varphi}+\scplt{\partialt-\nabla\tilde{u}}{\nabla\varphi},\\
\label{eq:Imain_2}
\scplt{\opdiv (p-\partialt)}{\opdiv\partialsi}+\scplt{p-\partialt}{\partialsi}
&=-\scplt{f-\tilde{u}+\opdiv\partialt}{\opdiv\partialsi}-\scplt{\partialt-\nabla\tilde{u}}{\partialsi}.
\end{align}
Adding \eqref{eq:Imain_1} and \eqref{eq:Imain_2} we obtain
\begin{align}
\label{eq:main_3}
\begin{split}
&\qquad\scplt{\nabla(u-\tilde{u})}{\nabla\varphi}+\scplt{u-\tilde{u}}{\varphi}
+\scplt{\opdiv(p-\partialt)}{\opdiv\partialsi}+\scplt{p-\partialt}{\partialsi}\\
&=\scplt{f-\tilde{u}+\opdiv\partialt}{\varphi-\opdiv\partialsi}+\scplt{\partialt-\nabla\tilde{u}}{\nabla\varphi-\partialsi}.
\end{split}
\end{align}
By choosing $\varphi:=u-\tilde{u}\in\hgen{1}{}ga$ and $\partialsi:=p-\partialt\in\dgen{}{}$, 
the left hand side of \eqref{eq:main_3} turns to the combined norm 
of the error of the approximation. Since we have
\begin{align*}
\varphi-\opdiv\partialsi
&=u-\tilde{u}-\opdiv p+\opdiv\partialt
=f-\tilde{u}+\opdiv\partialt,\\
\nabla\varphi-\partialsi
&=\nabla u-\nabla\tilde{u}-p+\partialt
=\partialt-\nabla\tilde{u},
\end{align*}
\eqref{eq:main_3} becomes \eqref{eq:Imain1}.
Putting $\tilde{u}=0$, $\partialt=0$ in \eqref{eq:Imain1} shows $\tnorm{(u,p)}=\normlt{f}$
and thus \eqref{eq:Imain2}.
\end{proof}

\begin{rem}
\label{thm:Imainrem}
\mbox{}
\begin{itemize}
\item[\bf(i)]
We note the similarity of the error majorants in 
Theorems \rgen{}{}ef{thm:Iest}, \rgen{}{}ef{thm:Iestp} and \rgen{}{}ef{thm:Imain}.
\item[\bf(ii)] 
It is clear that Theorem \rgen{}{}ef{thm:Imain} generalizes Theorems \rgen{}{}ef{thm:Iest} and \rgen{}{}ef{thm:Iestp}
since these two can be recovered from Theorem \rgen{}{}ef{thm:Imain}.
We just estimate
\begin{align*}
\mathcal{M}na(\tilde{u},p)
=\normho{u-\tilde{u}}^2
&\leq\tnorm{(u,p)-(\tilde{u},\partialt)}^2
=\mathcal{M}mix(\tilde{u},\partialt)=\mathcal{M}_\nabla(\tilde{u},\partialt)
\intertext{and note that the left hand side does not depend on $\partialsi:=\partialt\in\dgen{}{}$.
Analogously we estimate}
\mathcal{M}div(\partialt,u)
=\normd{p-\partialt}^2
&\leq\tnorm{(u,p)-(\tilde{u},\partialt)}^2
=\mathcal{M}mix(\tilde{u},\partialt)=\mathcal{M}div(\partialt,\tilde{u})
\end{align*}
and note that the left hand side does not depend on $\varphi:=\tilde{u}\in\hgen{1}{}ga$.
\end{itemize}
\end{rem}

\begin{rem} 
\label{rem:Iproof}
There is a simple proof of Theorem \rgen{}{}ef{thm:Imain}
using just \eqref{eq:Ipde} and $p=\nabla u$:
\begin{align*}
\mathcal{M}mix(\tilde{u},\partialt)
&=\normlt{f-\tilde{u}+\opdiv\partialt}^2+\normlt{\partialt-\nabla\tilde{u}}^2\\
&=\normlt{u-\tilde{u}+\opdiv\partialt-\opdiv p}^2+\normlt{\partialt-p+\nabla u-\nabla\tilde{u}}^2\\
&=\normlt{u-\tilde{u}}^2+\normlt{\opdiv(\partialt-p)}^2+2\scplt{u-\tilde{u}}{\opdiv(\partialt-p)}\\
&\qquad+\normlt{\partialt-p}^2+\normlt{\nabla(u-\tilde{u})}^2+2\scplt{\partialt-p}{\nabla(u-\tilde{u})}\\
&=\tnorm{(u,p)-(\tilde{u},\partialt)}^2
\end{align*}
In the last line we have used as before 
$\scplt{u - \tilde{u}}{\opdiv(\partialt - p)}=-\scplt{\nabla(u-\tilde{u})}{\partialt - p}$
since $u-\tilde{u}\in\hgen{1}{}ga$.
This shows immediately, that Theorem \rgen{}{}ef{thm:Imain} 
extends to more general situations as well.
E.g. inhomogeneous boundary conditions can be treated since
only $u-\tilde{u}\in\hgen{1}{}ga$ is needed. 
\end{rem}

\section{Results for the General Case} \label{sec:G}

In this section we derive our main result in an abstract setting which
allows for mixed boundary conditions as well as coefficients for the PDEs. 
We will prove the main result by using the simple approach presented in Remark \rgen{}{}ef{rem:Iproof}.

Let ${\hilbert_{1}}$ and ${\hilbert_{2}}$ be two Hilbert spaces with inner products 
$\scphio{\,\cdot\,}{\,\cdot\,}$ and $\scphit{\,\cdot\,}{\,\cdot\,}$, respectively. 
Moreover, let $\A:D(\A)\subset\hilbert_{1}\to\hilbert_{2}$ be a densely defined and closed linear operator
and $\A^{*}:D(\A^{*})\subset\hilbert_{2}\to\hilbert_{1}$ its adjoint. We note $\A^{**}=\bar{\A}=\A$ and 
\begin{align}
\label{partint}
\forall\,\varphi\in D(\A)\quad\forall\,\partialsi\in D(\A^{*})\qquad
\scphit{\A\varphi}{\partialsi}=\scphio{\varphi}{\A^{*}\partialsi}.
\end{align}
Equipped with the natural graph norms $D(\A)$ and $D(\A^{*})$ are Hilbert spaces. 
Furthermore, we introduce two linear, self adjoint and positive topological isomorphisms 
$\alpha_{1}:{\hilbert_{1}}\to{\hilbert_{1}}$ and $\alpha_{2}:{\hilbert_{2}}\to{\hilbert_{2}}$. Especially we have
$$\exists\,c>0\quad\forall\,\varphi\in\hilbert_{1}\qquad 
c^{-1}\normhio{\varphi}^2\leq\scphio{\alpha_{1}\varphi}{\varphi}\leq c\normhio{\varphi}^2$$
and the corresponding holds for $\alpha_{2}$.
For any inner product and corresponding norm we introduce weighted counterparts with sub-index notation. 
For example, for elements from ${\hilbert_{1}}$ we define a new inner product 
$\scphioao{\,\cdot\,}{\,\cdot\,}:=\scphio{\alpha_{1}\,\cdot\,}{\,\cdot\,}$ 
and a new induced norm $\normhioao{\,\cdot\,}$.
Using this notation we can define for $\varphi\in D(\A)$ and $\partialsi\in D(\A^{*})$ new weighted norms
on $D(\A)$, $D(\A^{*})$ as well as on the product space $D(\A)\times D(\A^{*})$ by
\begin{align*}
\normdaidaoat{\varphi}^2
&:=\normhioao{\varphi}^2
+\normhitat{\A\varphi}^2,\\
\normdasaomoatmo{\partialsi}^2
&:=\normhitatmo{\partialsi}^2
+\normhioaomo{\A^{*}\partialsi}^2,\\
\tnorm{(\varphi,\partialsi)}^2
&:=\normdaidaoat{\varphi}^2
+\normdasaomoatmo{\partialsi}^2.
\end{align*}

Let $f \in {\hilbert_{1}}$. By the Lax-Milgram lemma
(or by Riesz' representation theorem) we get immediately:

\begin{lem}
\label{laxmilgramA}
The (primal) variational problem
\begin{align}
\label{varA}
\forall\varphi\in D(\A)\qquad
\scphitat{\A x}{\A\varphi}+\scphioao{x}{\varphi}=\scphio{f}{\varphi}
\end{align}
admits a unique solution $x\in D(\A)$ satisfying $\normdaidaoat{x}\leq\normhioaomo{f}$. 
Moreover, $y_{x}:=\alpha_{2}\A x$ belongs to $D(\A^{*})$ and $\A^{*} y_{x}=f-\alpha_{1} x$.
Hence, the strong and mixed formulations
\begin{align} 
\label{eq:Gstrong}
\A^{*}\alpha_{2}\A x+\alpha_{1} x&=f,\\
\label{mixedformulation}
\A^{*} y_{x}+\alpha_{1} x&=f,\quad\alpha_{2}\A x=y_{x}
\end{align}
hold with $(x,y_{x})\in D(\A)\times \big(D(\A^{*})\times\alpha_{2} R(\A)\big)$.
\end{lem}

To get the dual problem, we multiply the first equation of \eqref{mixedformulation}
by $\A^{*}\partialsi$ with $\partialsi\in D(\A^{*})$ 
taking the right weighted scalar product and use $y_{x}=\alpha_{2}\A x\in D(\A^{*})$.
We obtain
$$\scphioaomo{\A^{*} y_{x}}{\A^{*}\partialsi}
+\scphioaomo{\alpha_{1} x}{\A^{*}\partialsi}
=\scphioaomo{f}{\A^{*}\partialsi}.$$
Since $x\in D(\A)$
$$\scphioaomo{\alpha_{1} x}{\A^{*}\partialsi}
=\scphio{x}{\A^{*}\partialsi}
=\scphit{\A x}{\partialsi}
=\scphitatmo{y_{x}}{\partialsi}$$
holds, we get again by the Lax-Milgram's lemma

\begin{lem}
\label{laxmilgramAs}
The (dual) variational problem
\begin{align}
\label{varAs}
\forall\partialsi\in D(\A^{*})\qquad
\scphioaomo{\A^{*} y}{\A^{*}\partialsi}
+\scphitatmo{y}{\partialsi}
=\scphioaomo{f}{\A^{*}\partialsi}
\end{align}
admits a unique solution $y\in D(\A^{*})$ satisfying $\normdasaomoatmo{y}\leq\normhioaomo{f}$.
Moreover, $y=y_{x}$ holds and thus $y$ even belongs to 
$D(\A^{*})c_{\A}p\alpha_{2} R(\A)$ with $x$ and $y_{x}$ from Lemma \rgen{}{}ef{laxmilgramA}.
Furthermore, $\alpha_{1}^{-1}(\A^{*} y-f)\in D(\A)$ with
$A\alpha_{1}^{-1}(\A^{*} y-f)=-\alpha_{2}^{-1}y$.
\end{lem}

\begin{proof}
We just have to show that $y_{x}\in D(\A^{*})$ solves \eqref{varAs}.
But this follows directly since for all $\partialsi\in D(\A^{*})$
\begin{align*}
\scphioaomo{\A^{*} y_{x}}{\A^{*}\partialsi}
&=-\scphio{x}{\A^{*}\partialsi}
+\scphioaomo{f}{\A^{*}\partialsi}\\
&=-\scphit{\A x}{\partialsi}
+\scphioaomo{f}{\A^{*}\partialsi}
=-\scphitatmo{y_{x}}{\partialsi}
+\scphioaomo{f}{\A^{*}\partialsi}.
\end{align*}
Hence $y_{x}=y$ and $\A^{**}=\A$ completes the proof.
\end{proof}

\begin{rem}
\label{rem:isometrygen}
We know $\normdaidaoat{x}\leq\normhioaomo{f}$ and $\normdasaomoatmo{y}\leq\normhioaomo{f}$.
It is indeed notable that 
$$\tnorm{(x,y)}=\normhioaomo{f}$$ 
holds, which follows immediately by $y=\alpha_{2}\A x$ and
\begin{align*}
\normhioaomo{f}^2
=\normhioaomo{\A^{*}\alpha_{2}\A x+\alpha_{1} x}^2
&=\normhioaomo{\A^{*} y}^2
+\normhioaomo{\alpha_{1} x}^2
+2\underbrace{\scphioaomo{\A^{*}\alpha_{2}\A x}{\alpha_{1} x}}_{\dgen{}{}isplaystyle=\scphio{\A^{*}\alpha_{2}\A x}{x}}\\
&=\normhioaomo{\A^{*} y}^2
+\normhioao{x}^2
+2\underbrace{\scphit{\alpha_{2}\A x}{\A x}}_{\dgen{}{}isplaystyle=\normhitat{\A x}^2}
=\tnorm{(x,y)}^2.
\end{align*}
Thus the solution operator
$$L:\hilbert_{1}\to D(\A)\times D(\A^{*});f\mapsto(x,y)$$
(equipped with the proper weighted norms)
has norm $\norm{L}=1$, i.e., $L$ is an isometry.
\end{rem}

By the latter remark the mixed norm on $D(\A)\times D(\A^{*})$ yields an isomtery.
This motivates to use the mixed norm also for error estimates.
As it turns out, we even obtain an error equality.
We present our main result of the paper.

\begin{theo} 
\label{thm:Gmain}
Let $(x,y),(\tilde{x},\tilde{y})\in D(\A)\times D(\A^{*})$ be the exact solution of \eqref{mixedformulation}
and any conforming approximation, respectively. Then
\begin{align} 
\label{eq:Gmain1}
\tnorm{(x,y)-(\tilde{x},\tilde{y})}^2=\mathcal{M}(\tilde{x},\tilde{y})
\end{align}
and the normalized counterpart
\begin{align} 
\label{eq:Gmain2}
\frac{\tnorm{(x,y)-(\tilde{x},\tilde{y})}^2}{\tnorm{(x,y)}^2}
=\frac{\mathcal{M}(\tilde{x},\tilde{y})}{\normhioaomo{f}^2} 
\end{align}
hold, where
\begin{align} 
\label{eq:Gmain}
\mathcal{M}(\tilde{x},\tilde{y}):=\normhioaomo{f-\alpha_{1}\tilde{x}-\A^{*}\tilde{y}}^2
+\normhitatmo{\tilde{y}-\alpha_{2}\A\tilde{x}}^2.
\end{align}
\end{theo}

\begin{proof}
Using \eqref{eq:Gstrong} and inserting $0=\alpha_{2}\A x-y$ we get by \eqref{partint}
\begin{align*}
\mathcal{M}(\tilde{x},\tilde{y})
&=\normhioaomo{\alpha_{1} x-\alpha_{1}\tilde{x}+\A^{*} y-\A^{*}\tilde{y}}^2
+\normhitatmo{\tilde{y}-y+\alpha_{2}\A x-\alpha_{2}\A\tilde{x}}^2\\
&=\normhioao{x-\tilde{x}}^2
+\normhioaomo{\A^{*}(y-\tilde{y})}^2
+2\scphioaomo{\alpha_{1}(x-\tilde{x})}{\A^{*}(y-\tilde{y})}\\
&\qquad+\normhitatmo{\tilde{y}-y}^2
+\normhitat{\A(x-\tilde{x})}^2
+2\scphitatmo{\tilde{y}-y}{\alpha_{2}\A(x-\tilde{x})}\\
&=\normdaidaoat{x-\tilde{x}}^2+\normdasaomoatmo{y-\tilde{y}}^2\\
&\qquad+2\scphio{x-\tilde{x}}{\A^{*}(y-\tilde{y})}
-2\scphit{\A(x-\tilde{x})}{y-\tilde{y}}\\
&=\tnorm{(x,y)-(\tilde{x},\tilde{y})}^2.
\end{align*}
\eqref{eq:Gmain2} follows by the isometry property in Remark \rgen{}{}ef{rem:isometrygen}, completing the proof.
\end{proof}

We note that the isometry property, i.e., $\tnorm{(x,y)}=\normhioaomo{f}$,
can be seen by inserting $(\tilde{x},\tilde{y})=(0,0)$ into \eqref{eq:Gmain1} as well.

\begin{rem}
Theorem \rgen{}{}ef{thm:Gmain} can also be deduced as a special case of the equation \cite[(7.2.14)]{NeittaanmakiRepin2004} in the book of Neittaam\"aki and Repin.
\end{rem}

\begin{rem} 
\label{rem:Gbehavior}
Of course, the majorant $\mathcal{M}$ is continuous. Especially we have
\begin{align*}
\mathcal{M}(\tilde{x},\tilde{y})\xrightarrow{\tilde{x}\rgen{}{}ightarrow x\text{ in }D(\A)}
&\normdasaomoatmo{y-\tilde{y}}^2=\mathcal{M}(x,\tilde{y}),\\
\mathcal{M}(\tilde{x},\tilde{y})\xrightarrow{\tilde{y}\rgen{}{}ightarrow y\text{ in }D(\A^{*})}
&\normdaidaoat{x-\tilde{x}}^2=\mathcal{M}(\tilde{x},y)
\end{align*}
and $\mathcal{M}(\tilde{x},\tilde{y})\to\mathcal{M}(x,y)=0$ if $(\tilde{x},\tilde{y})\to(x,y)$ in $D(\A)\times D(\A^{*})$.
This suggests that the majorant $\mathcal{M}$ can also be used as an error indicator for adaptive computations, 
even though the equality \eqref{eq:Gmain1} is global.
\end{rem}

\begin{cor}
\label{cor:maintheocor}
Theorem \rgen{}{}ef{thm:Gmain} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{x}\in D(\A)$ it holds
$\dgen{}{}isplaystyle\normdaidaoat{x-\tilde{x}}^2
=\min_{\partialsi\in D(\A^{*})}\mathcal{M}(\tilde{x},\partialsi)
=\mathcal{M}(\tilde{x},y)$.
\item[\bf(ii)]
For any $\tilde{y}\in D(\A^{*})$ it holds
$\dgen{}{}isplaystyle\normdasaomoatmo{y-\tilde{y}}^2
=\min_{\varphi\in D(\A)}\mathcal{M}(\varphi,\tilde{y})
=\mathcal{M}(x,\tilde{y})$.
\end{itemize}
\end{cor}

\begin{proof}
We just have to estimate
$$\normdaidaoat{x-\tilde{x}}^2
\leq\tnorm{(x,y)-(\tilde{x},\tilde{y})}^2
=\mathcal{M}(\tilde{x},\tilde{y})$$
and note that the left hand side does not depend on $\tilde{y}\in D(\A^{*})$.
Setting $\partialsi:=\tilde{y}\in D(\A^{*})$ we get
$$\normdaidaoat{x-\tilde{x}}^2
\leq\inf_{\partialsi\in D(\A^{*})}\mathcal{M}(\tilde{x},\partialsi).$$
But for $\partialsi=y\in D(\A^{*})$ we see $\mathcal{M}(\tilde{x},y)=\normdaidaoat{x-\tilde{x}}^2$,
which proves (i). Analogously, we estimate
$$\normdasaomoatmo{y-\tilde{y}}^2
\leq\tnorm{(x,y)-(\tilde{x},\tilde{y})}^2
=\mathcal{M}(\tilde{x},\tilde{y})$$
and note that the left hand side does not depend on $\tilde{x}\in D(\A)$.
Setting $\varphi:=\tilde{x}\in D(\A)$ we get
$$\normdasaomoatmo{y-\tilde{y}}^2
\leq\inf_{\varphi\in D(\A)}\mathcal{M}(\varphi,\tilde{y}).$$
But for $\varphi=x\in D(\A)$ we see $\mathcal{M}(x,\tilde{y})=\normdasaomoatmo{y-\tilde{y}}^2$,
which shows (ii).
\end{proof}

\begin{rem}
\label{laxmilgramAsrem}
\mbox{}
\begin{itemize}
\item[\bf(i)]
Since $y\,\bot_{\alpha_{2}^{-1}}\,N(\A^{*})$ by \eqref{varAs} we get immediately $y\in\alpha_{2}\overline{R(\A)}$
by the Helmholtz decomposition $\hilbert_{2}=N(\A^{*})\oplus_{\alpha_{2}^{-1}}\alpha_{2}\overline{R(\A)}$.
\item[\bf(ii)] 
If $\alpha_{1}^{-1}f\in D(\A)$ we have $z:=\alpha_{1}^{-1}\A^{*} y\in D(\A)$ 
and the strong and mixed formulations of \eqref{varAs} read
\begin{align*}
\A\alpha_{1}^{-1}\A^{*} y+\alpha_{2}^{-1}y&=\A\alpha_{1}^{-1}f,\\
\A z+\alpha_{2}^{-1}y&=\A\alpha_{1}^{-1}f,\quad \alpha_{1}^{-1}\A^{*} y=z.
\end{align*}
Then for all $\varphi\in D(\A)$ we have
\begin{align*}
\scphitat{\A z}{\A\varphi}+\scphioao{z}{\varphi}
&=-\scphit{y}{\A\varphi}+\scphioao{z}{\varphi}+\scphitat{\A\alpha_{1}^{-1}f}{A\varphi}\\
&=\scphitat{\A\alpha_{1}^{-1}f}{A\varphi}
\end{align*}
and hence $z\in\big(D(\A)c_{\A}p\alpha_{1}^{-1} R(\A^{*})\big)\subset D(\A)$ 
is the unique solution of this variational problem.
Furthermore, $\alpha_{2}(\A z-\A\alpha_{1}^{-1}f)\in D(\A^{*})$ 
and $\A^{*}\alpha_{2}(\A z-\A\alpha_{1}^{-1}f)=-\alpha_{1} z$.
If $\alpha_{2}\A\alpha_{1}^{-1}f$ belongs to $D(\A^{*})$ this yields 
$\alpha_{2}\A z\in D(\A^{*})$ and the strong equation
$$\A^{*}\alpha_{2}\A z+\alpha_{1} z=\A^{*}\alpha_{2}\A\alpha_{1}^{-1}f.$$
\end{itemize}
\end{rem}

Our error equalities may also be used to compute 
the radius of the indeterminacy set of solutions 
in terms of the radius of the indeterminacy set of right hand sides.
Often the right hand $f$ of a problem is not known exactly  
but known to belong to an indeterminacy ball around some known mean data $\hat{f}$.
Let us write $f=\hat{f}+f_{\mathtt{osc}}$.
Since the solution operator $L$ from Remark \rgen{}{}ef{rem:isometrygen} is an isometry, 
we have for the solutions $(x,y)=(\hat{x},\hat{y})+(x_{\mathtt{osc}},y_{\mathtt{osc}})$
$$\tnorm{(x_{\mathtt{osc}},y_{\mathtt{osc}})}
=\tnorm{Lf_{\mathtt{osc}}}
=\normhioaomo{f_{\mathtt{osc}}}.$$
Hence, the solutions belong to a ball of the same radius as the data.
In other words, any modeling error is mapped to an error of same size.
If the magnitude of the oscillating part $f_{\mathtt{osc}}$ is known, we also
know the magnitude of variations of the solution set.

\subsection{Application to Time Discretization}

One main application of our error equalities might be that equations of the type
\begin{align}
\label{AsAeqapp}
\A^{*}\alpha_{2}\A x+\alpha_{1} x=f
\end{align}
naturally occur in many types of time discretizations for plenty of linear wave propagation models.
A large class of wave propagation models, like electro-magnetics, acoustics or elasticity, have the structure 
$$(\partial_{t}\Lambdaambda^{-1}+\mathcal{M}ax)\begin{bmatrix}x\\y\end{bmatrix}
=\begin{bmatrix}g\\h\end{bmatrix},\quad
\mathcal{M}ax=\begin{bmatrix}0&-\A^{*}\\\A&0\end{bmatrix},\quad
\Lambdaambda=\begin{bmatrix}\lambda_{1}&0\\0&\lambda_{2}\end{bmatrix}$$
or
\begin{align}
\label{dteq}
\partial_{t}\lambda_{1}^{-1}x-\A^{*} y=g,\quad
\partial_{t}\lambda_{2}^{-1}y+\A x=h
\end{align}
with initial condition $(x,y)(0)=(x_{0},y_{0})$.
Often the material is assumed to be time-independent, i.e., $\Lambdaambda$ does not depend on time.
In this case $i\Lambdaambda\mathcal{M}ax$ is selfadjoint in the proper Hilbert spaces 
and the solution theory follows immediately by the spectral theorem.
We note that formally the second order wave equation 
$$\big(\partial_{t}^2-(\Lambdaambda\mathcal{M}ax)^2\big)\begin{bmatrix}x\\y\end{bmatrix}
=(\partial_{t}-\Lambdaambda\mathcal{M}ax)\Lambdaambda\begin{bmatrix}g\\h\end{bmatrix},\quad
(\Lambdaambda\mathcal{M}ax)^2
=\begin{bmatrix}-\lambda_{1}\A^{*}\lambda_{2}\A&0\\0&-\lambda_{2}\A\lambda_{1}\A^{*}\end{bmatrix}$$
holds. A standard implizit time discretization for \eqref{dteq} is e.g. the backward Euler scheme, i.e.,
$$\dgen{}{}elta_{n}^{-1}\lambda_{1}^{-1}(x_{n}-x_{n-1})-\A^{*} y_{n}=g_{n},\quad
\dgen{}{}elta_{n}^{-1}(y_{n}-y_{n-1})+\lambda_{2}\A x_{n}=\lambda_{2}h_{n},\quad
\dgen{}{}elta_{n}:=t_{n}-t_{n-1}.$$
Hence, we obtain e.g. for $x_{n}$
$$\A^{*}\lambda_{2}\A x_{n}+\dgen{}{}elta_{n}^{-2}\lambda_{1}^{-1}x_{n}
=f_{n}:=\A^{*}(\lambda_{2}h_{n}+\dgen{}{}elta_{n}^{-1}y_{n-1})+\dgen{}{}elta_{n}^{-2}\lambda_{1}^{-1}x_{n-1} + \dgen{}{}elta_n^{-1} g_n$$
provided that $\lambda_{2}h_{n}\in D(\A^{*})$.
Therefore \eqref{AsAeqapp} holds for $x_{n}$ with e.g. 
$\alpha_{1}=\dgen{}{}elta_{n}^{-2}\lambda_{1}^{-1}$ and $\alpha_{2}=\lambda_{2}$.
Of course, a similar equation holds for $y_{n}$ as well.
We note that our arguments extend to `all' practically used time discretizations.

Functional a posteriori error estimates for wave equations can be found 
in \cite{repinapostwave,paulyrepinrossihypmax}.

\section{Applications} \label{sec:A}

We will discuss some standard applications.
Let $\Omega\subset\mathbb{R}^d$, $d\geq1$.
Since we want to handle mixed boundary conditions, 
let us assume for simplicity, that $\Omega$ is a bounded or an exterior domain 
with (compact) Lipschitz continuous boundary $\Gamma$.
Moreover, let $\Gammad$ be an open subset of $\Gamma$ and $\Gamman:=\Gamma\setminus\overline{\Gammad}$
its complement. We will denote by $n$ the outward unit normal of the boundary. 
The results presented in this section are direct consequences of Theorem \rgen{}{}ef{thm:Gmain}
and, of course, Lemmas \rgen{}{}ef{laxmilgramA}, \rgen{}{}ef{laxmilgramAs}
and Remarks \rgen{}{}ef{rem:Gbehavior}, \rgen{}{}ef{rem:isometrygen}
as well as Corollary \rgen{}{}ef{cor:maintheocor} hold for all special applications.

\subsection{Reaction-Diffusion} \label{subsec:RD}

Find the scalar potential $u\in\hgen{1}{}$, such that
\begin{align}
-\opdiv\alpha\nabla u+\rgen{}{}ho\,u&=f&
\textrm{in }&\Omega,\nonumber\\
u&=0&
\textrm{on }&\Gammad,\label{eq:RD}\\
n\cdot\alpha\nabla u&=0&
\textrm{on }&\Gamman\nonumber.
\end{align}
The quadratic diffusion matrix $\alpha\in\lgen{\infty}{}$ is symmetric, real valued and uniformly positive definite. 
The reaction coefficient $\rgen{}{}ho\ge\rgen{}{}ho_0>0$ belongs to $\lgen{\infty}{}$
and the source $f$ to $\lgen{2}{}$. The dual variable for this problem is the flux $p=\alpha\nabla u\in\dgen{}{}$.
We need more Sobolev spaces
$$\hgen{1}{}gad:=\overline{\cgen{\infty}{\ga}d}^{\hgen{1}{}},\quad
\dgen{}{}gan:=\overline{\cgen{\infty}{\ga}n}^{\dgen{}{}},\quad
\dgen{}{}zgan:=\set{\partialsi\in\dgen{}{}gan}{\opdiv\partialsi=0},$$
where $\cgen{\infty}{\ga}d$ resp. $\cgen{\infty}{\ga}n$ are smooth test functions resp. vector fields 
having supports bounded away from $\Gammad$ resp. $\Gamman$.
The following table shows the relation to the notation of Section \rgen{}{}ef{sec:G}.
\begin{center}\begin{tabular}{c|c||c|c||c|c||c|c}
$\alpha_{1}$ & $\alpha_{2}$ & $\A$ & $\A^{*}$ & $\hilbert_{1}$ & $\hilbert_{2}$ & $D(\A)$ & $D(\A^{*})$ \\
\hline
$\rgen{}{}ho$ & $\alpha$ & $\nabla$ & $-\opdiv$ & $\lgen{2}{}$ & $\lgen{2}{}$ & $\hgen{1}{}gad$ & $\dgen{}{}gan$
\end{tabular}\end{center}
We note that indeed $D(\A^{*})=\dgen{}{}gan$ holds for Lipschitz domains, see e.g. \cite{jochmanncompembmaxmixbc},
which is not trivial at all. The relation \eqref{partint} reads now
$$\forall\,\varphi\in\hgen{1}{}gad\quad\forall\,\partialsi\in\dgen{}{}gan\qquad
\scplt{\nabla\varphi}{\partialsi}=-\scplt{\varphi}{\opdiv\partialsi}.$$
Considering the norms we have
\begin{align*}
\norm{u}_{\hgen{1}{},\rgen{}{}ho,\alpha}^2
&=\norm{u}_{\lgen{2}{},\rgen{}{}ho}^2
+\norm{\nabla u}_{\lgen{2}{},\alpha}^2,\\
\norm{p}_{\dgen{}{},\rgen{}{}ho^{-1},\alpha^{-1}}^2
&=\norm{p}_{\lgen{2}{},\alpha^{-1}}^2
+\norm{\opdiv p}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2,\\
\tnorm{(u,p)}^2
&=\norm{u}_{\hgen{1}{},\rgen{}{}ho,\alpha}^2
+\norm{p}_{\dgen{}{},\rgen{}{}ho^{-1},\alpha^{-1}}^2.
\end{align*}
Now \eqref{eq:RD} reads: Find $u\in\hgen{1}{}gad$ with $\alpha\nabla u\in\dgen{}{}gan$ such that
\begin{align}
\label{formapprd}
-\opdiv\alpha\nabla u+\rgen{}{}ho\,u=f.
\end{align}
Equivalently, in mixed formulation we have: Find $(u,p)\in\hgen{1}{}gad\times\dgen{}{}gan$ such that
\begin{align}
\label{mixedformapprd}
-\opdiv p+\rgen{}{}ho\,u=f,\quad\alpha\nabla u=p.
\end{align}
The primal and dual variational problems are: Find $(u,p)\in\hgen{1}{}gad\times\dgen{}{}gan$ such that
\begin{align*}
\forall\,\varphi&\in\hgen{1}{}gad&
\scp{\nabla u}{\nabla\varphi}_{\lgen{2}{},\alpha}+\scp{u}{\varphi}_{\lgen{2}{},\rgen{}{}ho}
&=\scplt{f}{\varphi},\\
\forall\,\partialsi&\in\dgen{}{}gan&
\scp{\opdiv p}{\opdiv\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}+\scp{p}{\partialsi}_{\lgen{2}{},\alpha^{-1}}
&=-\scp{f}{\opdiv\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}.
\end{align*}

\begin{theo} 
\label{thm:RD}
Let $(u,p),(\tilde{u},\partialt)\in\hgen{1}{}gad\times\dgen{}{}gan$ 
be the exact solution of \eqref{mixedformapprd} 
and any approximation, respectively. Then
$$\tnorm{(u,p)-(\tilde{u},\partialt)}^2
=\mathcal{M}rd(\tilde{u},\partialt),\quad
\frac{\tnorm{(u,p)-(\tilde{u},\partialt)}^2}{\tnorm{(u,p)}^2}
=\frac{\mathcal{M}rd(\tilde{u},\partialt)}{\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2}$$
hold, where $\mathcal{M}rd(\tilde{u},\partialt)=\norm{f-\rgen{}{}ho\tilde{u}+\opdiv\partialt}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2
+\norm{\partialt-\alpha\nabla\tilde{u}}_{\lgen{2}{},\alpha^{-1}}^2$.
\end{theo}

\begin{rem}
We note $\norm{u}_{\hgen{1}{},\rgen{}{}ho,\alpha}\leq\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}$ 
and $\norm{p}_{\dgen{}{},\rgen{}{}ho^{-1},\alpha^{-1}}\leq\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}$ and indeed
$$\tnorm{(u,p)}=\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}.$$
The solution operator $L:\lgen{2}{}\to\hgen{1}{}gad\times\dgen{}{}gan;f\mapsto(u,p)$ is an isometry, i.e. $\norm{L}=1$.
\end{rem}

\begin{cor}
Theorem \rgen{}{}ef{thm:RD} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{u}\in\hgen{1}{}gad$ it holds
$\dgen{}{}isplaystyle\norm{u-\tilde{u}}_{\hgen{1}{},\rgen{}{}ho,\alpha}^2
=\min_{\partialsi\in\dgen{}{}gan}\mathcal{M}rd(\tilde{u},\partialsi)
=\mathcal{M}rd(\tilde{u},p)$.
\item[\bf(ii)]
For any $\partialt\in\dgen{}{}gan$ it holds
$\dgen{}{}isplaystyle\norm{p-\partialt}_{\dgen{}{},\rgen{}{}ho^{-1},\alpha^{-1}}^2
=\min_{\varphi\in\hgen{1}{}gad}\mathcal{M}rd(\varphi,\partialt)
=\mathcal{M}rd(u,\partialt)$.
\end{itemize}
\end{cor}

\begin{rem}
We have $p=\alpha\nabla u\in\dgen{}{}ganc_{\A}p\,\alpha\nabla\hgen{1}{}gad$ and $u$ and $(u,p)$ 
solve \eqref{formapprd} and \eqref{mixedformapprd}, respectively.
Moreover, $\opdiv p+f\in\rgen{}{}ho\hgen{1}{}gad$
with $\nabla\rgen{}{}ho^{-1}(\opdiv p+f)=\alpha^{-1}p\in\nabla\hgen{1}{}gad=\rgen{}{}zgadc_{\A}p\,\mathcal{H}_{\gad,\gan}^{\bot}$.
Hence, for $f\in\rgen{}{}ho\hgen{1}{}$ we have $\opdiv p\in\rgen{}{}ho\hgen{1}{}$ and therefore
the strong and mixed formulations of the dual problem
\begin{align*}
-\nabla\rgen{}{}ho^{-1}\opdiv p+\alpha^{-1}p&=\nabla\rgen{}{}ho^{-1}f&
&&
\textrm{in }&\Omega,\\
\nabla v+\alpha^{-1}p&=\nabla\rgen{}{}ho^{-1}f,&
-\rgen{}{}ho^{-1}\opdiv p&=v&
\textrm{in }&\Omega
\intertext{hold, which are completed by the equations}
\opdiv p+f&=0&
&&
\textrm{on }&\Gammad,\\
n\cdot p&=0&
&&
\textrm{on }&\Gamman,\\
\rgen{}{}ot\alpha^{-1}p&=0&
&&
\textrm{in }&\Omega,\\
n\times\alpha^{-1}p&=0&
&&
\textrm{on }&\Gammad,\\
\alpha^{-1}p&\,\,\bot\,\,\mathcal{H}_{\gad,\gan}.
\end{align*}
Here the Dirichlet-Neumann fields $\mathcal{H}_{\gad,\gan}$ and
the space $\rgen{}{}zgad$ will be defined in Section \rgen{}{}ef{subsec:EC}.
Of course, $\rgen{}{}ho v=f$ on $\Gammad$
and by $\rgen{}{}ho v\in\opdiv\dgen{}{}gan$ we also have $\rgen{}{}ho v\bot\,\mathbb{R}$ if $\Gamma=\Gamman$.
\end{rem}

Related results and numerical tests for exterior domains can be found in e.g.
\cite{paulyrepinell,malimuzalevskiypaulyextell}.

\subsection{Eddy-Current (3D)} \label{subsec:EC}

Let $d=3$. The problem reads: Find the electric field $E\in\rgen{}{}$ such that
\begin{align}
\rgen{}{}ot\mu^{-1}\rgen{}{}ot E+\epsilon E&=J&
\textrm{in }&\Omega,\nonumber\\
n\times E&=0&
\textrm{on }&\Gammad,\label{eq:EC}\\
n\times\mu^{-1}\rgen{}{}ot E&=0&
\textrm{on }&\Gamman,\nonumber
\end{align}
where
$$\rgen{}{}:=\set{\Phi\in\lgen{2}{}}{\rgen{}{}ot\Phi\in\lgen{2}{}},\quad
\rgen{}{}z:=\set{\Phi\in\rgen{}{}}{\rgen{}{}ot\Phi=0}.$$
We assume that the magnetic permeability $\mu$ and the electric permittivity $\epsilon$ 
are symmetric, real valued and uniformly positive definite matrices 
from $\lgen{\infty}{}$. Of course, the extension to complex valued matrices is straight forward.
The electric current $J$ belongs to $\lgen{2}{}$. 
The dual variable for this problem is the magnetic field $H=\mu^{-1}\rgen{}{}ot E\in\rgen{}{}$.
We define the Sobolev spaces
$$\rgen{}{}gad:=\overline{\cgen{\infty}{\ga}d}^{\rgen{}{}},\quad
\rgen{}{}zgad:=\set{\Phi\in\rgen{}{}gad}{\rgen{}{}ot\Phi=0}$$
and analogously $\rgen{}{}gan$ and $\rgen{}{}zgan$.
Moreover, we introduce the co-called Dirichlet-Neumann and Neumann-Dirichlet fields by
\begin{align*}
\mathcal{H}_{\gad,\gan}&:=\rgen{}{}zgadc_{\A}p\dgen{}{}zgan=\set{\Psi\in\rgen{}{}gadc_{\A}p\dgen{}{}gan}{\rgen{}{}ot\Psi=0\,\wedge\,\opdiv\Psi=0},\\
\mathcal{H}_{\gan,\gad}&:=\rgen{}{}zganc_{\A}p\dgen{}{}zgad=\set{\Psi\in\rgen{}{}ganc_{\A}p\dgen{}{}gad}{\rgen{}{}ot\Psi=0\,\wedge\,\opdiv\Psi=0},
\end{align*}
respectively. The following table shows the relation to the notation of Section \rgen{}{}ef{sec:G}.
\begin{center}\begin{tabular}{c|c||c|c||c|c||c|c}
$\alpha_{1}$ & $\alpha_{2}$ & $\A$ & $\A^{*}$ & $\hilbert_{1}$ & $\hilbert_{2}$ & $D(\A)$ & $D(\A^{*})$\\
\hline
$\epsilon$ & $\mu^{-1}$ & $\rgen{}{}ot$ & $\rgen{}{}ot$ & $\lgen{2}{}$ & $\lgen{2}{}$ & $\rgen{}{}gad$ & $\rgen{}{}gan$
\end{tabular}\end{center}
We note that indeed $D(\A^{*})=\rgen{}{}gan$ holds for Lipschitz domains, see e.g. \cite{jochmanncompembmaxmixbc},
which is not trivial at all. The relation \eqref{partint} reads now
$$\forall\,\Phi\in\rgen{}{}gad\quad\forall\,\Psi\in\rgen{}{}gan\qquad
\scplt{\rgen{}{}ot\Phi}{\Psi}=\scplt{\Phi}{\rgen{}{}ot\Psi}.$$
Considering the norms we have
\begin{align*}
\norm{E}_{\rgen{}{},\epsilon,\mu^{-1}}^2
&=\norm{E}_{\lgen{2}{},\epsilon}^2
+\norm{\rgen{}{}ot E}_{\lgen{2}{},\mu^{-1}}^2,\\
\norm{H}_{\rgen{}{},\epsilon^{-1},\mu}^2
&=\norm{H}_{\lgen{2}{},\mu}^2
+\norm{\rgen{}{}ot H}_{\lgen{2}{},\epsilon^{-1}}^2,\\
\tnorm{(E,H)}^2
&=\norm{E}_{\rgen{}{},\epsilon,\mu^{-1}}^2
+\norm{H}_{\rgen{}{},\epsilon^{-1},\mu}^2.
\end{align*}
Now \eqref{eq:EC} reads: Find $E\in\rgen{}{}gad$ with $\mu^{-1}\rgen{}{}ot E\in\rgen{}{}gan$ such that
$$\rgen{}{}ot\mu^{-1}\rgen{}{}ot E+\epsilon E=J.$$
In mixed formulation we have: Find $(E,H)\in\rgen{}{}gad\times\rgen{}{}gan$ such that
$$\rgen{}{}ot H+\epsilon E=J,\quad\mu^{-1}\rgen{}{}ot E=H.$$
The primal and dual variational problems are: Find $(E,H)\in\rgen{}{}gad\times\rgen{}{}gan$ such that
\begin{align*}
\forall\,\Phi&\in\rgen{}{}gad&
\scp{\rgen{}{}ot E}{\rgen{}{}ot\Phi}_{\lgen{2}{},\mu^{-1}}+\scp{E}{\Phi}_{\lgen{2}{},\epsilon}
&=\scplt{J}{\Phi},\\
\forall\,\Psi&\in\rgen{}{}gan&
\scp{\rgen{}{}ot H}{\rgen{}{}ot\Psi}_{\lgen{2}{},\epsilon^{-1}}+\scp{H}{\Psi}_{\lgen{2}{},\mu}
&=\scp{J}{\rgen{}{}ot\Psi}_{\lgen{2}{},\epsilon^{-1}}.
\end{align*}

\begin{theo} 
\label{thm:EC}
For any approximation $(\tilde{E},\tilde{H})\in\rgen{}{}gad\times\rgen{}{}gan$ 
$$\tnorm{(E,H)-(\tilde{E},\tilde{H})}^2
=\mathcal{M}ec(\tilde{E},\tilde{H}),\quad
\frac{\tnorm{(E,H)-(\tilde{E},\tilde{H})}^2}{\tnorm{(E,H)}^2}
=\frac{\mathcal{M}ec(\tilde{E},\tilde{H})}{\norm{J}_{\lgen{2}{},\epsilon^{-1}}^2}$$
hold, where $\mathcal{M}ec(\tilde{E},\tilde{H})=\norm{J-\epsilon\tilde{E}-\rgen{}{}ot\tilde{H}}_{\lgen{2}{},\epsilon^{-1}}^2
+\norm{\tilde{H}-\mu^{-1}\rgen{}{}ot\tilde{E}}_{\lgen{2}{},\mu}^2$.
\end{theo}

\begin{rem}
\label{rem:EC}
We note $\norm{E}_{\rgen{}{},\epsilon,\mu^{-1}}\leq\norm{J}_{\lgen{2}{},\epsilon^{-1}}$ 
and $\norm{H}_{\rgen{}{},\epsilon^{-1},\mu}\leq\norm{J}_{\lgen{2}{},\epsilon^{-1}}$ and indeed
$$\tnorm{(E,H)}=\norm{J}_{\lgen{2}{},\epsilon^{-1}}.$$
The solution operator $L:\lgen{2}{}\to\rgen{}{}gad\times\rgen{}{}gan;f\mapsto(E,H)$ is an isometry, i.e. $\norm{L}=1$.
\end{rem}

\begin{cor}
\label{cor:EC}
Theorem \rgen{}{}ef{thm:EC} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{E}\in\rgen{}{}gad$ it holds
$\dgen{}{}isplaystyle\norm{E-\tilde{E}}_{\rgen{}{},\epsilon,\mu^{-1}}^2
=\min_{\Psi\in\rgen{}{}gan}\mathcal{M}ec(\tilde{E},\Psi)
=\mathcal{M}ec(\tilde{E},H)$.
\item[\bf(ii)]
For any $\tilde{H}\in\rgen{}{}gan$ it holds
$\dgen{}{}isplaystyle\norm{H-\tilde{H}}_{\rgen{}{},\epsilon^{-1},\mu}^2
=\min_{\Phi\in\rgen{}{}gad}\mathcal{M}ec(\Phi,\tilde{H})
=\mathcal{M}ec(E,\tilde{H})$.
\end{itemize}
\end{cor}

\begin{rem}
We have $H=\mu^{-1}\rgen{}{}ot E\in\rgen{}{}ganc_{\A}p\,\mu^{-1}\rgen{}{}ot\rgen{}{}gad$ and $E$ and $(E,H)$ 
solve the strong and mixed formulation, respectively.
Moreover, $\rgen{}{}ot H-J\in\epsilon\rgen{}{}gad$
with $\rgen{}{}ot\epsilon^{-1}(\rgen{}{}ot H-J)=-\mu H$ 
belonging to $\rgen{}{}ot\rgen{}{}gad=\dgen{}{}zgadc_{\A}p\,\mathcal{H}_{\gan,\gad}^{\bot}$.
Hence, for $J\in\epsilon\rgen{}{}$ we have $\rgen{}{}ot H\in\epsilon\rgen{}{}$ and therefore
the strong and mixed formulations of the dual problem
\begin{align*}
\rgen{}{}ot\epsilon^{-1}\rgen{}{}ot H+\mu H&=\rgen{}{}ot\epsilon^{-1}J&
&&
\textrm{in }&\Omega,\\
\rgen{}{}ot D+\mu H&=\rgen{}{}ot\epsilon^{-1}J,&
\epsilon^{-1}\rgen{}{}ot H&=D&
\textrm{in }&\Omega
\intertext{hold, which are completed by the equations}
n\times\epsilon^{-1}(\rgen{}{}ot H-J)&=0&
&&
\textrm{on }&\Gammad,\\
n\times H&=0&
&&
\textrm{on }&\Gamman,\\
\opdiv\mu H&=0&
&&
\textrm{in }&\Omega,\\
n\cdot\mu H&=0&
&&
\textrm{on }&\Gammad,\\
\mu H&\,\,\bot\,\,\mathcal{H}_{\gan,\gad}.
\end{align*}
Of course, $n\times D=n\times\epsilon^{-1}J$ on $\Gammad$
and by $\epsilon D\in\rgen{}{}ot\rgen{}{}gan$ we also have
$\opdiv\epsilon D=0$ in $\Omega$ and $n\cdot\epsilon D=0$ on $\Gamman$
as well as $\epsilon D\bot\mathcal{H}_{\gad,\gan}$.
\end{rem}

Earlier results for eddy current and static Maxwell problems can be found in 
\cite{anjammalimuzalevskiyneittaanmakirepinmaxtype,paulyrepinmaxst}.

\subsection{Eddy-Current (2D)} \label{subsec:EC2D}

Let $d=2$. We just indicate the changes compared to the latter section.
First, we have to understand the double $\rgen{}{}ot$ as $\nabla^\bot\rgen{}{}ot$, where
$$\rgen{}{}ot E:=\opdiv\R E=\partial_1E_2-\partial_2E_1,\quad
\nabla^\bot H:=\R\nabla H=\begin{bmatrix}\partial_2H\\-\partial_1H\end{bmatrix},\quad
\R:=\begin{bmatrix}0&1\\-1&0\end{bmatrix}$$
and $E\in\rgen{}{}$ is a vector field and $H\in\hgen{1}{}$ a scalar function. 
In the literature, the operator $\nabla^\bot$ is often called 
co-gradient or vector rotation $\vec{\rgen{}{}ot}$ as well. 
Also $\mu$ is scalar. \eqref{eq:EC} reads: Find the electric field $E\in\rgen{}{}$ such that
\begin{align*}
\nabla^{\bot}\mu^{-1}\rgen{}{}ot E+\epsilon E&=J&
\textrm{in }&\Omega,\\
n\times E&=0&
\textrm{on }&\Gammad,\\
\mu^{-1}\rgen{}{}ot E&=0&
\textrm{on }&\Gamman.
\end{align*}
We have
\begin{center}\begin{tabular}{c|c||c|c||c|c||c|c}
$\alpha_{1}$ & $\alpha_{2}$ & $\A$ & $\A^{*}$ & $\hilbert_{1}$ & $\hilbert_{2}$ & $D(\A)$ & $D(\A^{*})$\\
\hline
$\epsilon$ & $\mu^{-1}$ & $\rgen{}{}ot$ & $\nabla^{\bot}$ & $\lgen{2}{}$ & $\lgen{2}{}$ & $\rgen{}{}gad$ & $\hgen{1}{}gan$
\end{tabular}\end{center}
and \eqref{partint} turns to
$$\forall\,\Phi\in\rgen{}{}gad\quad\forall\,\partialsi\in\hgen{1}{}gan\qquad
\scplt{\rgen{}{}ot\Phi}{\partialsi}=\scplt{\Phi}{\nabla^{\bot}\partialsi}.$$
The norm for $H$ is
$$\norm{H}_{\hgen{1}{},\epsilon^{-1},\mu}^2
=\norm{H}_{\lgen{2}{},\mu}^2
+\norm{\nabla^{\bot}H}_{\lgen{2}{},\epsilon^{-1}}^2.$$
The strong formulation of the problem is: 
Find $E\in\rgen{}{}gad$ with $\mu^{-1}\rgen{}{}ot E\in\hgen{1}{}gan$ such that
$$\nabla^{\bot}\mu^{-1}\rgen{}{}ot E+\epsilon E=J.$$
The mixed formulation is: Find $(E,H)\in\rgen{}{}gad\times\hgen{1}{}gan$ such that
$$\nabla^{\bot}H+\epsilon E=J,\quad\mu^{-1}\rgen{}{}ot E=H.$$
The primal and dual variational problems are: Find $(E,H)\in\rgen{}{}gad\times\hgen{1}{}gan$ such that
\begin{align*}
\forall\,\Phi&\in\rgen{}{}gad&
\scp{\rgen{}{}ot E}{\rgen{}{}ot\Phi}_{\lgen{2}{},\mu^{-1}}+\scp{E}{\Phi}_{\lgen{2}{},\epsilon}
&=\scplt{J}{\Phi},\\
\forall\,\partialsi&\in\hgen{1}{}gan&
\scp{\nabla^{\bot}H}{\nabla^{\bot}\partialsi}_{\lgen{2}{},\epsilon^{-1}}+\scp{H}{\partialsi}_{\lgen{2}{},\mu}
&=\scp{J}{\nabla^{\bot}\partialsi}_{\lgen{2}{},\epsilon^{-1}}.
\end{align*}

Theorem \rgen{}{}ef{thm:EC} reads:

\begin{theo} 
\label{thm:EC2D}
For any approximation $(\tilde{E},\tilde{H})\in\rgen{}{}gad\times\hgen{1}{}gan$ 
$$\tnorm{(E,H)-(\tilde{E},\tilde{H})}^2
=\mathcal{M}ec(\tilde{E},\tilde{H}),\quad
\frac{\tnorm{(E,H)-(\tilde{E},\tilde{H})}^2}{\tnorm{(E,H)}^2}
=\frac{\mathcal{M}ec(\tilde{E},\tilde{H})}{\norm{J}_{\lgen{2}{},\epsilon^{-1}}^2}$$
hold, where $\mathcal{M}ec(\tilde{E},\tilde{H})=\norm{J-\epsilon\tilde{E}-\nabla^{\bot}\tilde{H}}_{\lgen{2}{},\epsilon^{-1}}^2
+\norm{\tilde{H}-\mu^{-1}\rgen{}{}ot\tilde{E}}_{\lgen{2}{},\mu}^2$.
\end{theo} 

\begin{rem}
\label{rem:EC2D}
We note $\norm{E}_{\rgen{}{},\epsilon,\mu^{-1}}\leq\norm{J}_{\lgen{2}{},\epsilon^{-1}}$ 
and $\norm{H}_{\hgen{1}{},\epsilon^{-1},\mu}\leq\norm{J}_{\lgen{2}{},\epsilon^{-1}}$ and indeed
$$\tnorm{(E,H)}=\norm{J}_{\lgen{2}{},\epsilon^{-1}}.$$
The solution operator $L:\lgen{2}{}\to\rgen{}{}gad\times\hgen{1}{}gan;f\mapsto(E,H)$ is an isometry, i.e. $\norm{L}=1$.
\end{rem}

\begin{cor}
\label{cor:EC2D}
Theorem \rgen{}{}ef{thm:EC} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{E}\in\rgen{}{}gad$ it holds
$\dgen{}{}isplaystyle\norm{E-\tilde{E}}_{\rgen{}{},\epsilon,\mu^{-1}}^2
=\min_{\partialsi\in\hgen{1}{}gan}\mathcal{M}ec(\tilde{E},\partialsi)
=\mathcal{M}ec(\tilde{E},H)$.
\item[\bf(ii)]
For any $\tilde{H}\in\hgen{1}{}gan$ it holds
$\dgen{}{}isplaystyle\norm{H-\tilde{H}}_{\hgen{1}{},\epsilon^{-1},\mu}^2
=\min_{\Phi\in\rgen{}{}gad}\mathcal{M}ec(\Phi,\tilde{H})
=\mathcal{M}ec(E,\tilde{H})$.
\end{itemize}
\end{cor}

\begin{rem}
We have again $H=\mu^{-1}\rgen{}{}ot E\in\hgen{1}{}ganc_{\A}p\,\mu^{-1}\rgen{}{}ot\rgen{}{}gad$ and
as in the 3D case $E$ and $(E,H)$ 
solve the strong and mixed formulation, respectively.
Moreover, $\nabla^{\bot}H-J\in\epsilon\rgen{}{}gad$
with $\rgen{}{}ot\epsilon^{-1}(\nabla^{\bot}H-J)=-\mu H$.
Hence, for $J\in\epsilon\rgen{}{}$ we have $\nabla^{\bot}H\in\epsilon\rgen{}{}$ and therefore
the strong and mixed formulations of the dual problem
\begin{align*}
\rgen{}{}ot\epsilon^{-1}\nabla^{\bot}H+\mu H&=\rgen{}{}ot\epsilon^{-1}J&
&&
\textrm{in }&\Omega,\\
\rgen{}{}ot D+\mu H&=\rgen{}{}ot\epsilon^{-1}J,&
\epsilon^{-1}\nabla^{\bot}H&=D&
\textrm{in }&\Omega
\intertext{hold, which are completed by the equations}
n\times\epsilon^{-1}(\nabla^{\bot}H-J)&=0&
&&
\textrm{on }&\Gammad,\\
H&=0&
&&
\textrm{on }&\Gamman,\\
\mu H&\,\,\bot\,\,\mathbb{R}\quad(\text{if }\Gammad=\Gamma).
\end{align*}
Of course, $n\times D=n\times\epsilon^{-1}J$ on $\Gammad$
and by $\epsilon D\in\nabla^{\bot}\hgen{1}{}gan$ we also have
$\opdiv\epsilon D=0$ in $\Omega$ and $n\cdot\epsilon D=0$ on $\Gamman$
as well as $\epsilon D\bot\mathcal{H}_{\gad,\gan}$.
\end{rem}

\subsection{Linear Elasticity} \label{subsec:LE}

Find the displacement vector field $u\in\hgen{1}{}$ such that
\begin{align}
-\Div\Lambda\nablas u+\rgen{}{}ho\,u&=f&
\textrm{in }&\Omega,\nonumber\\
u&=0&
\textrm{on }&\Gammad,\label{eq:LE}\\
n\cdot\Lambda\nablas u&=0&
\textrm{on }&\Gamman.\nonumber
\end{align}
Here $\nablas$ is the symmetric part of the gradient\footnote{Here, 
as usual in elasticity the gradient $\nabla u$
is to be understood as the Jacobian of the vector field $u$.}
$$\nablas u:=\sym\nabla u=\frac{1}{2}\big(\nabla u+(\nabla u)^\top\big),$$
where ${}^\top$ denotes the transpose. 
$\nablas u$, often denoted by $\epsilon(u)$, is also called the infinitesimal strain tensor.
The fourth order stiffness tensor of elastic moduli $\Lambda\in\lgen{\infty}{}$,
mapping symmetric matrices to symmetric matrices point-wise,
and the second order tensor (quadratic matrix) of reaction $\rgen{}{}ho$ 
are assumed to be symmetric, real valued and uniformly positive definite.
The vector field $f$ (body force) belongs to $\lgen{2}{}$ and
the dual variable for this problem is the Cauchy stress tensor $\sigma=\Lambda\nablas u\in\dgen{}{}$,
where the application of $\Div$ to $\sigma$ and the notation $\sigma\in\dgen{}{}$ 
is to be understood row-wise as the usual divergence $\opdiv$.
We note that the first equation can also be written as
$$-\Div_{\mathrm{s}}\,\Lambda\nablas u+\rgen{}{}ho\,u=f,\quad\Div_{\mathrm{s}}\,:=\Div\sym.$$
We have:
\begin{center}\begin{tabular}{c|c||c|c||c|c||c|c}
$\alpha_{1}$ & $\alpha_{2}$ & $\A$ & $\A^{*}$ & $\hilbert_{1}$ & $\hilbert_{2}$ & $D(\A)$ & $D(\A^{*})$\\
\hline
$\rgen{}{}ho$ & $\Lambda$ & $\nablas$ & $-\Div_{\mathrm{s}}\,$ & $\lgen{2}{}$ & $\lgen{2}{}$ & $\hgen{1}{}gad$ & $\sym^{-1}\dgen{}{}gan$
\end{tabular}\end{center}
The notation $\sigma\in\sym^{-1}\dgen{}{}gan$ means $\sym\sigma\in\dgen{}{}gan$.
More precisely, $\partialsi\in D(\A^{*})$ if and only if
$$\forall\,\varphi\in D(\A)=\hgen{1}{}gad\qquad
\scplt{\nablas\varphi}{\partialsi}=\scplt{\varphi}{\A^{*}\partialsi}.$$
Since $\scplt{\nablas\varphi}{\partialsi}=\scplt{\nabla\varphi}{\sym\partialsi}$
we see that this holds if and only if $\sym\partialsi\in\dgen{}{}gan$ and $\A^{*}\partialsi=-\Div\sym\partialsi$.
Equation \eqref{partint} turns into
$$\forall\,\varphi\in\hgen{1}{}gad\quad\forall\,\partialsi\in\sym^{-1}\dgen{}{}gan
\qquad\scplt{\nablas\varphi}{\partialsi}=-\scplt{\varphi}{\Div_{\mathrm{s}}\,\partialsi}.$$
For the norms we have
\begin{align*}
\norm{u}_{\hgen{1}{},\rgen{}{}ho,\Lambda}^2
&=\norm{u}_{\lgen{2}{},\rgen{}{}ho}^2
+\norm{\nablas u}_{\lgen{2}{},\Lambda}^2,\\
\norm{\sigma}_{\sym^{-1}\dgen{}{},\rgen{}{}ho^{-1},\Lambda^{-1}}^2
&=\norm{\sigma}_{\lgen{2}{},\Lambda^{-1}}^2
+\norm{\Div_{\mathrm{s}}\,\sigma}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2,\\
\tnorm{(u,\sigma)}^2
&=\norm{u}_{\hgen{1}{},\rgen{}{}ho,\Lambda}^2
+\norm{\sigma}_{\sym^{-1}\dgen{}{},\rgen{}{}ho^{-1},\Lambda^{-1}}^2.
\end{align*}
Now \eqref{eq:LE} reads: Find $u\in\hgen{1}{}gad$ with $\sym\Lambda\nablas u=\Lambda\nablas u\in\dgen{}{}gan$ such that
$$-\Div\Lambda\nablas u+\rgen{}{}ho\,u=f.$$
In mixed formulation we have: Find $(u,\sigma)\in\hgen{1}{}gad\times\dgen{}{}gan$ such that
$$-\Div\sigma+\rgen{}{}ho\,u=f,\quad\Lambda\nablas u=\sigma.$$
Note that then $\sigma$ is automatically symmetric.
The primal and dual variational problems are: Find $(u,\sigma)\in\hgen{1}{}gad\times\sym^{-1}\dgen{}{}gan$ such that
\begin{align*}
\forall\,\varphi&\in\hgen{1}{}gad&
\scp{\nablas u}{\nablas\varphi}_{\lgen{2}{},\Lambda}+\scp{u}{\varphi}_{\lgen{2}{},\rgen{}{}ho}
&=\scplt{f}{\varphi},\\
\forall\,\partialsi&\in\sym^{-1}\dgen{}{}gan&
\scp{\Div_{\mathrm{s}}\,\sigma}{\Div_{\mathrm{s}}\,\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}+\scp{\sigma}{\partialsi}_{\lgen{2}{},\Lambda^{-1}}
&=-\scp{f}{\Div_{\mathrm{s}}\,\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}.
\intertext{Since $\sigma\in\dgen{}{}gan$ must be symmetric, we can formulate the dual problem also as}
\forall\,\partialsi&\in\dgen{}{}gan,\;\partialsi\text{ symmetric}&
\scp{\Div\sigma}{\Div\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}+\scp{\sigma}{\partialsi}_{\lgen{2}{},\Lambda^{-1}}
&=-\scp{f}{\Div\partialsi}_{\lgen{2}{},\rgen{}{}ho^{-1}}.
\end{align*}
Then, the norms reduce to
$$\tnorm{(u,\sigma)}^2
=\norm{u}_{\hgen{1}{},\rgen{}{}ho,\Lambda}^2
+\norm{\sigma}_{\dgen{}{},\rgen{}{}ho^{-1},\Lambda^{-1}}^2,\quad
\norm{\sigma}_{\dgen{}{},\rgen{}{}ho^{-1},\Lambda^{-1}}^2
=\norm{\sigma}_{\lgen{2}{},\Lambda^{-1}}^2
+\norm{\Div\sigma}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2.$$

\begin{theo} 
\label{thm:LE}
For any approximation $(\tilde{u},\tilde{\sigma})\in\hgen{1}{}gad\times\sym^{-1}\dgen{}{}gan$ 
\begin{align}
\label{linelatheoeq}
\tnorm{(u,\sigma)-(\tilde{u},\tilde{\sigma})}^2
=\mathcal{M}le(\tilde{u},\tilde{\sigma}),\quad
\frac{\tnorm{(u,\sigma)-(\tilde{u},\tilde{\sigma})}^2}{\tnorm{(u,\sigma)}^2}
=\frac{\mathcal{M}le(\tilde{u},\tilde{\sigma})}{\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2}
\end{align}
hold, where $\mathcal{M}le(\tilde{u},\tilde{\sigma})=\norm{f-\rgen{}{}ho\tilde{u}+\Div_{\mathrm{s}}\,\tilde{\sigma}}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2
+\norm{\tilde{\sigma}-\Lambda\nablas\tilde{u}}_{\lgen{2}{},\Lambda^{-1}}^2$.
Moreover, since $\sigma$ is automatically symmetric we have \eqref{linelatheoeq} for all
$(\tilde{u},\tilde{\sigma})\in\hgen{1}{}gad\times\dgen{}{}gan$ with $\tilde{\sigma}$ symmetric
and the right hand side simplifies to
$\mathcal{M}le(\tilde{u},\tilde{\sigma})=\norm{f-\rgen{}{}ho\tilde{u}+\Div\tilde{\sigma}}_{\lgen{2}{},\rgen{}{}ho^{-1}}^2
+\norm{\tilde{\sigma}-\Lambda\nablas\tilde{u}}_{\lgen{2}{},\Lambda^{-1}}^2$.
\end{theo}

\begin{rem}
We note $\norm{u}_{\hgen{1}{},\rgen{}{}ho,\Lambda}\leq\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}$ 
and $\norm{\sigma}_{\dgen{}{},\rgen{}{}ho^{-1},\Lambda^{-1}}\leq\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}$ and indeed
$$\tnorm{(u,\sigma)}=\norm{f}_{\lgen{2}{},\rgen{}{}ho^{-1}}.$$
The solution operator $L:\lgen{2}{}\to\hgen{1}{}gad\times\dgen{}{}gan;f\mapsto(u,\sigma)$ is an isometry, i.e. $\norm{L}=1$.
\end{rem}

\begin{cor}
Theorem \rgen{}{}ef{thm:LE} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{u}\in\hgen{1}{}gad$ it holds
$\dgen{}{}isplaystyle\norm{u-\tilde{u}}_{\hgen{1}{},\rgen{}{}ho,\alpha}^2
=\min_{\partialsi\in\sym^{-1}\dgen{}{}gan}\mathcal{M}le(\tilde{u},\partialsi)
=\mathcal{M}le(\tilde{u},\sigma)$.
\item[\bf(ii)]
For any $\tilde{\sigma}\in\sym^{-1}\dgen{}{}gan$ it holds
$\dgen{}{}isplaystyle\norm{\sigma-\tilde{\sigma}}_{\sym^{-1}\dgen{}{},\rgen{}{}ho^{-1},\alpha^{-1}}^2
=\min_{\varphi\in\hgen{1}{}gad}\mathcal{M}le(\varphi,\tilde{\sigma})
=\mathcal{M}le(u,\tilde{\sigma})$.
\end{itemize}
If $\tilde{\sigma}$ and $\partialsi$ are already symmetric we can skip the $\sym^{-1}$
and replace $\Div_{\mathrm{s}}\,$ by $\Div$.
\end{cor}

\begin{rem}
We have $\sigma=\Lambda\nablas u\in\dgen{}{}ganc_{\A}p\,\Lambda\nablas\hgen{1}{}gad$ is symmetric with $\Div_{\mathrm{s}}\,\sigma=\Div\sigma$
and $u$ and $(u,\sigma)$ solve the strong and mixed formulation, respectively.
Moreover, $\Div\sigma+f\in\rgen{}{}ho\hgen{1}{}gad$
with $\nablas\rgen{}{}ho^{-1}(\Div\sigma+f)=\Lambda^{-1}\sigma\in\nablas\hgen{1}{}gad$.
Hence, for $f\in\rgen{}{}ho\hgen{1}{}$ we have $\Div\sigma\in\rgen{}{}ho\hgen{1}{}$ and therefore
strong and mixed formulations of the dual problem hold, i.e.,
\begin{align*}
-\nablas\rgen{}{}ho^{-1}\Div\sigma+\Lambda^{-1}\sigma&=\nablas\rgen{}{}ho^{-1}f&
&&
\textrm{in }&\Omega,\\
\nablas v+\Lambda^{-1}\sigma&=\nablas\rgen{}{}ho^{-1}f,&
-\rgen{}{}ho^{-1}\Div\sigma&=v&
\textrm{in }&\Omega.
\end{align*}
\end{rem}

\subsection{Generalized Reaction-Diffusion, Linear Accoustics and Eddy-Current} \label{subsec:diff}

Let $\Omega$ be a $d$-dimensional smooth Riemannian manifold 
with compact Lipschitz boundary $\Gamma$. If $\Omega$ is unbounded, 
we assume that outside of some compact set, $\Omega$ is isomorphic 
to the exterior unit domain $\set{x\in\mathbb{R}^{d}}{|x|>1}$.
Moreover, let $\Gammad$ be an open subset of $\Gamma$ and $\Gamman:=\Gamma\setminus\overline{\Gammad}$ its complement.
The problem reads:
For $f\in\lgen{2}{}q$ find the differential form potential ($q$-form) $u\in\dgen{}{}q$, such that
\begin{align}
-\cd\alpha\ed u+\rgen{}{}ho\,u&=f&\textrm{in }&\Omega,\nonumber\\
\ttr{\ga}d u&=0&\textrm{on }&\Gammad,\label{eq:diff}\\
\ntr{\ga}n \alpha\ed u&=0&\textrm{on }&\Gamman\nonumber.
\end{align}
Here, $\ed$ denotes exterior derivative, $\cd=\partialm*\ed*$ the co-derivative
and $\ttr{\ga}d$ resp. $\ntr{\ga}n$ the restrictions of the tangential resp. normal
traces $\ttr{\ga}$ resp. $\ntr{\ga}$ to the proper subspaces.
We also introduce the Sobolev spaces
$$\dgen{}{}q:=\set{\varphi\in\lgen{2}{}q}{\ed\varphi\in\lgen{2}{}qpo},\quad
\dgen{}{}eq:=\set{\partialsi\in\lgen{2}{}q}{\cd\partialsi\in\lgen{2}{}qmo}$$
and $\dgen{}{}qgad:=\overline{\cgen{\infty,q}{\gad}}^{\dgen{}{}q}$, $\dgen{}{}eqgan:=\overline{\cgen{\infty,q}{\gan}}^{\dgen{}{}eq}$,
where $\cgen{\infty,q}{\gad}$ resp. $\cgen{\infty,q}{\gan}$ are smooth test $q$-forms
having supports bounded away from $\Gammad$ resp. $\Gamman$.
Moreover, $\lgen{2}{}q$ denotes the Lebesgue space of all square integrable $q$-forms on $\Omega$
equipped with the inner or scalar product 
$$\scpltq{u}{\varphi}:=\int_{\Omega}u\wedge*\varphi$$
and corresponding norm $\normltq{\,\cdot\,}$. Of course, $\dgen{}{}q$ and $\dgen{}{}eq$
are equipped with the respective graph norms, making them Hilbert spaces.
Finally, $\rgen{}{}ho$ and $\alpha$ denote linear, symmetric, real valued, bounded 
and uniformly positive definite transformations on $q$- resp. $(q+1)$-forms. 
It is again straight forward to discuss complex valued transformations.
We also need the spaces
$$\dgen{}{}qz:=\set{\varphi\in\dgen{}{}q}{\ed\varphi=0},\quad
\dgen{}{}qzgad:=\set{\varphi\in\dgen{}{}qgad}{\ed\varphi=0}$$
and the corresponding spaces for the co-derivative as well as 
the space of harmonic Dirichlet-Neumann forms
$$\mathcal{H}^q_{\gad,\gan}:=\dgen{}{}qzgadc_{\A}p\dgen{}{}eqzgan.$$
The dual variable for this problem is the `flux' $p=\alpha\ed u\in\dgen{}{}eqpo$.
The next table shows the relation to the notations of Section \rgen{}{}ef{sec:G}.
\begin{center}\begin{tabular}{c|c||c|c||c|c||c|c}
$\alpha_{1}$ & $\alpha_{2}$ & $\A$ & $\A^{*}$ & $\hilbert_{1}$ & $\hilbert_{2}$ & $D(\A)$ & $D(\A^{*})$\\
\hline
$\rgen{}{}ho$ & $\alpha$ & $\ed$ & $-\cd$ & $\lgen{2}{}q$ & $\lgen{2}{}qpo$ & $\dgen{}{}qgad$ & $\dgen{}{}eqpogan$
\end{tabular}\end{center}
Also here indeed $D(\A^{*})=\dgen{}{}eqpogan$ holds, see 
e.g. \cite{goldshteinmitreairinamariushodgedecomixedbc,jakabmitreairinamariusfinensolhodgedeco,kuhndiss}. 
The relation \eqref{partint} turns into
$$\forall\,\varphi\in\dgen{}{}qgad\quad\forall\,\partialsi\in\dgen{}{}eqpogan\qquad
\scpltqpo{\ed\varphi}{\partialsi}=-\scpltq{\varphi}{\cd\partialsi}.$$
Considering the norms we have
\begin{align*}
\norm{u}_{\dgen{}{}q,\rgen{}{}ho,\alpha}^2
&=\norm{u}_{\lgen{2}{}q,\rgen{}{}ho}^2
+\norm{\ed u}_{\lgen{2}{}qpo,\alpha}^2,\\
\norm{p}_{\dgen{}{}eqpo,\rgen{}{}ho^{-1},\alpha^{-1}}^2
&=\norm{p}_{\lgen{2}{}qpo,\alpha^{-1}}^2
+\norm{\cd p}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}^2,\\
\tnorm{(u,p)}^2
&=\norm{u}_{\dgen{}{}q,\rgen{}{}ho,\alpha}^2
+\norm{p}_{\dgen{}{}eqpo,\rgen{}{}ho^{-1},\alpha^{-1}}^2.
\end{align*}
Now \eqref{eq:diff} reads: Find $u\in\dgen{}{}qgad$ with $\alpha\ed u\in\dgen{}{}eqpogan$ such that
$$-\cd\alpha\ed u+\rgen{}{}ho\,u=f.$$
In mixed formulation we have: Find $(u,p)\in\dgen{}{}qgad\times\dgen{}{}eqpogan$ such that
$$-\cd p+\rgen{}{}ho\,u=f,\quad\alpha\ed u=p.$$
The primal and dual variational problems are: Find $(u,p)\in\dgen{}{}qgad\times\dgen{}{}eqpogan$ such that
\begin{align*}
\forall\,\varphi&\in\dgen{}{}qgad&
\scp{\ed u}{\ed\varphi}_{\lgen{2}{}qpo,\alpha}+\scp{u}{\varphi}_{\lgen{2}{}q,\rgen{}{}ho}
&=\scpltq{f}{\varphi},\\
\forall\,\partialsi&\in\dgen{}{}eqpogan&
\scp{\cd p}{\cd\partialsi}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}+\scp{p}{\partialsi}_{\lgen{2}{}qpo,\alpha^{-1}}
&=-\scp{f}{\cd\partialsi}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}.
\end{align*}

\begin{theo} 
\label{thm:diff}
For any approximation $(\tilde{u},\partialt)\in\dgen{}{}qgad\times\dgen{}{}eqpogan$ 
$$\tnorm{(u,p)-(\tilde{u},\partialt)}^2
=\mathcal{M}diff(\tilde{u},\partialt),\quad
\frac{\tnorm{(u,p)-(\tilde{u},\partialt)}^2}{\tnorm{(u,p)}^2}
=\frac{\mathcal{M}diff(\tilde{u},\partialt)}{\norm{f}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}^2}$$
hold, where $\mathcal{M}diff(\tilde{u},\partialt)=\norm{f-\rgen{}{}ho\tilde{u}+\cd\partialt}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}^2
+\norm{\partialt-\alpha\ed\tilde{u}}_{\lgen{2}{}qpo,\alpha^{-1}}^2$.
\end{theo}

\begin{rem}
We note $\norm{u}_{\dgen{}{}q,\rgen{}{}ho,\alpha}\leq\norm{f}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}$ 
and $\norm{p}_{\dgen{}{}eqpo,\rgen{}{}ho^{-1},\alpha^{-1}}\leq\norm{f}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}$ and indeed
$$\tnorm{(u,p)}=\norm{f}_{\lgen{2}{}q,\rgen{}{}ho^{-1}}.$$
The solution operator $L:\lgen{2}{}q\to\dgen{}{}qgad\times\dgen{}{}eqpogan;f\mapsto(u,p)$ is an isometry, i.e. $\norm{L}=1$.
\end{rem}

\begin{cor}
Theorem \rgen{}{}ef{thm:diff} provides the a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{u}\in\dgen{}{}qgad$ it holds
$\dgen{}{}isplaystyle\norm{u-\tilde{u}}_{\dgen{}{}q,\rgen{}{}ho,\alpha}^2
=\min_{\partialsi\in\dgen{}{}eqpogan}\mathcal{M}diff(\tilde{u},\partialsi)
=\mathcal{M}diff(\tilde{u},p)$.
\item[\bf(ii)]
For any $\partialt\in\dgen{}{}eqpogan$ it holds
$\dgen{}{}isplaystyle\norm{p-\partialt}_{\dgen{}{}eqpo,\rgen{}{}ho^{-1},\alpha^{-1}}^2
=\min_{\varphi\in\dgen{}{}qgad}\mathcal{M}diff(\varphi,\partialt)
=\mathcal{M}diff(u,\partialt)$.
\end{itemize}
\end{cor}

We note that for $q=0$ we get back the reaction-diffusion problem from Section \rgen{}{}ef{subsec:RD}
and for $d=3$ or $d=2$ and $q=1$ we obtain the eddy-current problems 
from Sections \rgen{}{}ef{subsec:EC} and \rgen{}{}ef{subsec:EC2D},
identifying $\Omega\subset\mathbb{R}^{d}$ with a proper domain and
$0$-forms with functions and $1$- and $2$-forms with vector fields
by Riesz' representation theorem and Hodge's star operator.

\begin{rem}
It holds $p=\alpha\ed u\in\dgen{}{}eqpoganc_{\A}p\,\alpha\ed\dgen{}{}qgad$ and $u$ and $(u,p)$ 
solve the strong and mixed formulations, respectively.
Moreover, $\cd p+f$ belongs to $\rgen{}{}ho\dgen{}{}qgad$ and we see immediately 
$\ed\rgen{}{}ho^{-1}(\cd p+f)=\alpha^{-1}p\in\ed\dgen{}{}qgad=\dgen{}{}qpozgadc_{\A}p\,(\mathcal{H}^{q+1}_{\gad,\gan})^{\bot}$.
Hence, for $f\in\rgen{}{}ho\dgen{}{}q$ we have $\cd p\in\rgen{}{}ho\dgen{}{}q$ and therefore
the strong and mixed formulations of the dual problem
\begin{align*}
-\ed\rgen{}{}ho^{-1}\cd p+\alpha^{-1}p&=\ed\rgen{}{}ho^{-1}f&
&&
\textrm{in }&\Omega,\\
\ed v+\alpha^{-1}p&=\ed\rgen{}{}ho^{-1}f,&
-\rgen{}{}ho^{-1}\cd p&=v&
\textrm{in }&\Omega
\intertext{hold, which are completed by the equations}
\ttr{\ga}d\rgen{}{}ho^{-1}(\cd p+f)&=0&
&&
\textrm{on }&\Gammad,\\
\ttr{\ga}n p&=0&
&&
\textrm{on }&\Gamman,\\
\ed\alpha^{-1}p&=0&
&&
\textrm{in }&\Omega,\\
\ttr{\ga}d\alpha^{-1}p&=0&
&&
\textrm{on }&\Gammad,\\
\alpha^{-1}p&\,\,\bot\,\,\mathcal{H}^{q+1}_{\gad,\gan}.
\end{align*}
Of course, there are also more equations for $v$
following from $\rgen{}{}ho v\in\rgen{}{}ho\dgen{}{}qc_{\A}p\cd\dgen{}{}eqpogan$, e.g. $\cd\rgen{}{}ho v=0$,
which we will not list here explicitly.
\end{rem}

\section{Inhomogeneous and More Boundary Conditions}
\label{sec:inhomobc}

In this section we will demonstrate that our error equalities also hold 
for Robin type boundary conditions, which means that our error equalities are true
for many commonly used boundary conditions. Moreover, we emphasize that we can also handle
inhomogeneous boundary conditions.
Since it is clear that this method works in the general setting as well
we will discuss it here just for the simple reaction-diffusion model problem from the introduction.

Let $\Omega$ be as in the latter section and now the boundary $\Gamma$ 
be decomposed into three disjoint parts $\Gammad$, $\Gamman$ and $\Gammar$.
The model problem is:
Find the scalar potential $u\in\hgen{1}{}$ such that
\begin{align*}
-\opdiv\nabla u+u&=f&\textrm{in }&\Omega,\\
u&=g_1&\textrm{on }&\Gammad,\\
n\cdot\nabla u&=g_2&\textrm{on }&\Gamman,\\
n\cdot\nabla u+\Gammamma u&=g_3&\textrm{on }&\Gammar
\end{align*}
hold. Hence, on $\Gammad, \Gamman$ and $\Gammar$ we impose Dirichlet, Neumann 
and Robin type boundary conditions, respectively. 
In the Robin boundary condition, we assume that the coefficient $\Gammamma\ge\Gammamma_0>0$ 
belongs to $\lgen{\infty}{}$. The dual variable for this problem 
is the flux $p:=\nabla u\in\dgen{}{}$.
Furthermore, as long as $\Gammar\neq\emptyset$ and
to avoid tricky discussions about traces 
and the corresponding $\hgen{-\nicefrac{1}{2}}{}$-spaces of $\Gamma$, $\Gammad, \Gamman$ and $\Gammar$,
which can be quite complicated, we assume for simplicity that
$u\in\hgen{2}{}$. Then, $p\in\hgen{1}{}$ and all $g_{i}$ belong to $\lgen{2}{}$ even to $\hgen{1}{}h$ of $\Gamma$.
For the norms we simply have
$$\tnorm{(u,p)}^2
=\normho{u}^2
+\normd{p}^2.$$

\begin{theo} 
\label{thm:dnr}
For any approximation $(\tilde{u},\partialt)\in\hgen{2}{}\times\hgen{1}{}$ with
$u-\tilde{u}\in\hgen{1}{}gad$ and $p-\partialt\in\dgen{}{}gan$ as well as
$n\cdot(p-\partialt)+\Gammamma(u-\tilde{u})=0$ on $\Gammar$
$$\tnorm{(u,p)-(\tilde{u},\partialt)}^2
+\norm{u-\tilde{u}}_{\lgen{2}{}(\Gammar),\Gammamma}^2
+\norm{n\cdot(p-\partialt)}_{\lgen{2}{}(\Gammar),\Gammamma^{-1}}^2
=\mathcal{M}mix(\tilde{u},\partialt)$$
holds with $\mathcal{M}mix$ from Theorem \rgen{}{}ef{thm:Imain}.
Moreover, $\norm{u-\tilde{u}}_{\lgen{2}{}(\Gammar),\Gammamma}=\norm{n\cdot(p-\partialt)}_{\lgen{2}{}(\Gammar),\Gammamma^{-1}}$.
\end{theo}

\begin{proof}
Following Remark \rgen{}{}ef{rem:Iproof} we have
\begin{align*}
\mathcal{M}mix(\tilde{u},\partialt)
&=\underbrace{\normho{u-\tilde{u}}^2
+\normd{p-\partialt}^2}_{\dgen{}{}isplaystyle=\tnorm{(u,p)-(\tilde{u},\partialt)}^2}
+2\scplt{\nabla(u-\tilde{u})}{\partialt-p}
+2\scplt{u-\tilde{u}}{\opdiv(\partialt-p)}.
\end{align*}
Moreover, since $n\cdot(\partialt-p)$ and $u-\tilde{u}$ belong to $\lgen{2}{}(\Gamma)$ we have 
\begin{align*}
&\qquad\scplt{\nabla(u-\tilde{u})}{\partialt-p}
+\scplt{u-\tilde{u}}{\opdiv(\partialt-p)}\\
&=\scp{n\cdot(\partialt-p)}{u-\tilde{u}}_{\lgen{2}{}(\Gamma)}
=\scp{n\cdot(\partialt-p)}{u-\tilde{u}}_{\lgen{2}{}(\Gammar)}
=\scp{\Gammamma(u-\tilde{u})}{u-\tilde{u}}_{\lgen{2}{}(\Gammar)}.
\end{align*}
As $\scp{\Gammamma(u-\tilde{u})}{u-\tilde{u}}_{\lgen{2}{}(\Gammar)}=\scp{\Gammamma^{-1}n\cdot(p-\partialt)}{n\cdot(p-\partialt)}_{\lgen{2}{}(\Gammar)}$
we get the assertion.
\end{proof}

\begin{rem}
If all $g_{i}=0$, we can set $(\tilde{u},\partialt)=(0,0)$ and get
$$\tnorm{(u,p)}^2
+\norm{u}_{\lgen{2}{}(\Gammar),\Gammamma}^2
+\norm{n\cdot p}_{\lgen{2}{}(\Gammar),\Gammamma^{-1}}^2
=\normlt{f}^2,$$
which follows also directly from Remark \rgen{}{}ef{thm:Imainrem} (ii'),
$p=\nabla u$ and $n\cdot p=-\Gammamma u$ on $\Gammar$ as well as
\begin{align*}
\normlt{f}^2
&=\normlt{\opdiv p}^2
+\normlt{u}^2
-2\scplt{\opdiv\nabla u}{u}\\
&=\normlt{\opdiv p}^2
+\normlt{u}^2
+2\normlt{\nabla u}
-2\scp{n\cdot\nabla u}{u}_{\lgen{2}{}(\Gamma)}\\
&=\normlt{\opdiv p}^2
+\normlt{u}^2
+2\normlt{\nabla u}
-2\underbrace{\scp{n\cdot\nabla u}{u}_{\lgen{2}{}(\Gammar)}}_{\dgen{}{}isplaystyle=-\norm{u}_{\lgen{2}{}(\Gammar),\Gammamma}^2}.
\end{align*}
Thus, in this case the assertion of Theorem \rgen{}{}ef{thm:dnr} has a normalized counterpart as well.
\end{rem}

If $\Gammar=\emptyset$ we have a pure mixed Dirichlet and Neumann boundary. 

\begin{theo} 
\label{thm:dnrnorobin}
Let $\Gammar=\emptyset$.
For any approximation $(\tilde{u},\partialt)\in\hgen{1}{}\times\dgen{}{}$ with
$u-\tilde{u}\in\hgen{1}{}gad$ and $p-\partialt\in\dgen{}{}gan$
$$\tnorm{(u,p)-(\tilde{u},\partialt)}^2
=\mathcal{M}mix(\tilde{u},\partialt)$$
holds with $\mathcal{M}mix$ from Theorem \rgen{}{}ef{thm:Imain}.
\end{theo}

\begin{cor}
Let $\Gammar=\emptyset$.
Theorem \rgen{}{}ef{thm:dnrnorobin} provides the well known a posteriori error estimates 
for the primal and dual problems. 
\begin{itemize}
\item[\bf(i)]
For any $\tilde{u}\in\hgen{1}{}$ with $u-\tilde{u}\in\hgen{1}{}gad$ it holds
$\dgen{}{}isplaystyle\normho{u-\tilde{u}}^2
=\min_{\substack{\partialsi\in\dgen{}{}\\partial-\partialsi\in\dgen{}{}gan}}\mathcal{M}mix(\tilde{u},\partialsi)
=\mathcal{M}mix(\tilde{u},p)$.
\item[\bf(ii)]
For any $\partialt\in\dgen{}{}$ with $p-\partialt\in\dgen{}{}gan$ it holds
$\dgen{}{}isplaystyle\normd{p-\partialt}^2
=\min_{\substack{\varphi\in\hgen{1}{}\\u-\varphi\in\hgen{1}{}gad}}\mathcal{M}mix(\varphi,\partialt)
=\mathcal{M}mix(u,\partialt)$.
\end{itemize}
\end{cor}

\section{Numerical Examples} \label{sec:N}

In this section we show by some academic test cases 
the numerical performance of our error equalities.
All the calculations have been done using MATLAB, 
and the reported values in the tables have not been rounded, 
but are simply cut-offs of values reported by MATLAB. 
The main quantity of interest is the difference between the exact error 
and the value given by the majorant for a certain approximation $(\tilde{u},\partialt)$, i.e.,
$$\dgen{}{}elta:=\big|\tnorm{(u,p)-(\tilde{u},\partialt)}-\sqrts{\mathcal{M}_{\cdots}(\tilde{u},\partialt)}\big|,$$
where the test problems are either from the reaction-diffusion problems from subsection \rgen{}{}ef{subsec:RD} 
or from the eddy-current problems from subsections \rgen{}{}ef{subsec:EC} and \rgen{}{}ef{subsec:EC2D}. 
Where the finite element method (FEM) has been used, 
we have employed only linear triangular elements in 2D and linear tetrahedral elements in 3D. 
In all the examples below we calculated the approximations $\tilde{u}$ and $\partialt$ (or $\tilde{E}$ and $\tilde{H}$)
in the same mesh only for the sake of convenience. 
Using different meshes for the primal and dual approximations is allowed. 
We also used only regular meshes, but irregular meshes can be used as well. 
The only requirement is that the approximations must be conforming, 
meaning that they belong to the appropriate Sobolev spaces 
and fulfill the boundary conditions exactly.
All finite element solvers were implemented in the vectorized manner explained in \cite{rahmanvaldmanfastmatlab}.

\begin{ex}
\label{ex1}
We take the 3D-reaction-diffusion problem from Section \rgen{}{}ef{subsec:RD}
and choose the unit cube $\Omega:=(0,1)^3$ with exact solution 
$$u(x):=\partialrod_{i=1}^{3}x_{i}(1-x_{i}),$$
where $u$ satisfies the zero Dirichlet boundary conditions on the whole boundary, i.e.,
$\Gammad=\Gamma$ and $\Gamman=\emptyset$,
and the following data

$$\alpha(x):=\alpha:=\begin{bmatrix}1&0&0\\0&5&0\\0&0&10\end{bmatrix},\quad
\rgen{}{}ho(x):=\begin{cases}
1 & \textrm{if} \quad 0 < x_1 < \nicefrac{1}{4} \\
10 & \textrm{if} \quad \nicefrac{1}{4} < x_1 < \nicefrac{3}{4} \\
25 & \textrm{if} \quad \nicefrac{3}{4} < x_1 < 1
\end{cases}.$$
This means that the approximation 
of the dual variable does not have any boundary condition. 
We calculated the approximation globally by solving 
the primal and dual problem
with standard linear Courant elements and linear Raviart-Thomas elements, respectively. 
We will denote this finite element approximation pair by $(u_{\mathtt{h}},\partialh)$. 
The resulting linear systems were solved directly in MATLAB. 
The approximations were calculated in uniformly refined regular meshes, 
where the jumps in the reaction coefficient $\rgen{}{}ho$ coincide with element boundaries. 
For each mesh we computed the exact combined error and the majorant $\mathcal{M}rd(u_{\mathtt{h}},\partialh)$. 
The results are displayed in Table \rgen{}{}ef{tbl:Ex1}. 
The first column shows the number of elements $N_{\mathtt{elem}}$ of the mesh. 
The second and third column show the exact error and the value given by the majorant. 
The fourth column shows the difference $\dgen{}{}elta$ between the exact error and the value given by the majorant.
\end{ex}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex1} (3D-reaction-diffusion)}
\begin{center}
\begin{tabular}{r|c|c|c}
$N_{\mathtt{elem}}$ & $\tnorm{(u,p)-(u_{\mathtt{h}},\partialh)}$ & $\sqrts{\mathcal{M}rd(u_{\mathtt{h}},\partialh)}$ & difference $\dgen{}{}elta$ \\
\hline\hline
384    & 0.12803218100 & 0.12803218100 & 5.551115123e-17 \\
3072   & 0.06736516349 & 0.06736516349 & 4.163336342e-17 \\
24576  & 0.03433600867 & 0.03433600867 & 9.714451465e-17  \\
196608 & 0.01728806289 & 0.01728806289 & 3.469446952e-18
\end{tabular}
\end{center}
\label{tbl:Ex1}
\end{table}

\begin{ex}
\label{ex2}
This test is similar to the Example 1 except that the linear systems resulting 
from the finite element computations were not solved directly, but with an iterative method,
where the stopping tolerance was set to the crude value of $10^{-4}$.
The approximation pair obtained by this method is denoted by $(u_{\mathtt{iter}},\partialiter)$. 
No preconditioning was done. The iterative solver of the linear system of the dual problem converged 
only for the smallest mesh, and the error actually grows between the two last meshes. 
With this stopping tolerance this is expected
and was purposefully done so in order to obtain approximations 
which are relatively far from having the Galerkin orthogonality property. 
We did this test simply to demonstrate that Galerkin orthogonality 
is not a requirement for the equality to hold. 
The results are displayed in Table \rgen{}{}ef{tbl:Ex2}.
\end{ex}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex2} (3D-reaction-diffusion)}
\begin{center}
\begin{tabular}{r|c|c|c}
$N_{\mathtt{elem}}$ & $\tnorm{(u,p)-(u_{\mathtt{iter}},\partialiter)}$ & $\sqrts{\mathcal{M}rd(u_{\mathtt{iter}},\partialiter)}$ & difference $\dgen{}{}elta$ \\
\hline\hline
384    & 0.12803483290 & 0.12803483290 & 2.775557562e-17 \\
3072   & 0.06868358511 & 0.06868358511 & 6.938893904e-17 \\
24576  & 0.05294561599 & 0.05294561599 & 6.245004514e-17 \\
196608 & 0.09166231565 & 0.09166231565 & 9.714451465e-17
\end{tabular}
\end{center}
\label{tbl:Ex2}
\end{table}

\begin{ex}
\label{ex3}
We ran the problem data of Example \rgen{}{}ef{ex1} with subsequently refined regular meshes, 
where the approximation of the primal variable $u_{\mathtt{h}}$ was again obtained by the linear Courant finite elements. 
The resulting linear system was solved directly. 
The approximation of the dual variable was calculated by averaging the values 
$\alpha\nablau_{\mathtt{h}}$ to the nodes of the mesh. This procedure is often called the gradient averaging method
and we will denote the resulting function by $\partialav$. 
The results can be seen in Table \rgen{}{}ef{tbl:Ex3}.
\end{ex}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex3} (3D-reaction-diffusion)}
\begin{center}
\begin{tabular}{r|c|c|c}
$N_{\mathtt{elem}}$ & $\tnorm{(u,p)-(u_{\mathtt{h}},\partialav)}$ & $\sqrts{\mathcal{M}rd(u_{\mathtt{h}},\partialav)}$ & difference $\dgen{}{}elta$ \\
\hline\hline
384    & 0.2698605861 & 0.2698605861 & 0 \\
3072   & 0.2285323585 & 0.2285323585 & 0 \\
24576  & 0.1831121412 & 0.1831121412 & 6.106226635e-16 \\
196608 & 0.1333268308 & 0.1333268308 & 1.693090113e-15
\end{tabular}
\end{center}
\label{tbl:Ex3}
\end{table}

\begin{ex}
\label{ex4}
We take the 2D-eddy-current problem from Section \rgen{}{}ef{subsec:EC2D} 
and choose the unit square $\Omega:=(0,1)^2$ 
with $\epsilon=\id$ and $\mu=1$. We split the domain in two parts
$\Omega_1:=\set{x\in\Omega}{x_1>x_2}$ and $\Omega_2=\Omega\setminus\overline{\Omega_1}$ 
in order to define the following discontinuous solution
$$E|_{\Omega_1}(x):=
\begin{bmatrix}
\sin(2\partiali x_1)+2\partiali\cos(2\partiali x_1)(x_1-x_2)\\
\sin\big((x_1-x_2)^2(x_1-1)^2x_2\big)-\sin(2\partiali x_1)
\end{bmatrix},\quad
E|_{\Omega_2}(x):=0.$$
Note that indeed $E\in\rgen{}{}\setminus\hgen{1}{}$ and $\rgen{}{}ot E\in\hgen{1}{}$ with
$$\rgen{}{}ot E|_{\Omega_1}(x)=2x_{2}(x_1-x_2)(x_1-1)(2x_{1}-x_{2}-1)\cos(2\partiali x_1).$$
We set zero Neumann boundary conditions on the whole boundary, i.e.,
$\Gammad=\emptyset$ and $\Gamman=\Gamma$. 
The exact solution and its rotation is visualized in Figure \rgen{}{}ef{fig:ex4}.
We calculated the approximation globally by solving 
the primal and dual problem with linear N\'ed\'elec elements and linear Courant elements, respectively. 
This finite element approximation pair will be denoted by $(E_{\mathtt{h}},H_{\mathtt{h}})$. 
The resulting linear systems were solved directly. 
The approximations were calculated in uniformly refined regular meshes, 
where the jumps in the exact solution and in the right hand side $J$ coincide with element boundaries. 
For each mesh we calculated the exact combined error and the majorant $\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})$. 
The results are displayed in Table \rgen{}{}ef{tbl:Ex4}.
\end{ex}

\begin{figure}\label{fig:ex4}
\end{figure}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex4} (2D-eddy-current)}
\begin{center}
\begin{tabular}{r|c|c|c}
$N_{\mathtt{elem}}$ & $\tnorm{(E,H)-(E_{\mathtt{h}},H_{\mathtt{h}})}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}$ & difference $\dgen{}{}elta$ \\
\hline\hline
800    & 0.151485078300 & 0.151485078300 & 2.220446049e-16 \\
3200   & 0.075877018950 & 0.075877018950 & 0 \\
12800  & 0.037956449900 & 0.037956449900 & 7.632783294e-17 \\
51200  & 0.018980590110 & 0.018980590110 & 6.938893904e-17 \\
204800 & 0.009490605462 & 0.009490605462 & 2.602085214e-17
\end{tabular}
\end{center}
\label{tbl:Ex4}
\end{table}

\begin{ex}
\label{ex5}
We take the 3D-eddy-current problem from Section \rgen{}{}ef{subsec:EC} 
and choose the unit cube $\Omega:=(0,1)^3$ 
with $\epsilon=\mu=\id$. Again we split the domain in two parts
$\Omega_1:=\set{x\in\Omega}{x_1>x_2}$ and $\Omega_2=\Omega\setminus\overline{\Omega_1}$ 
in order to define the following discontinuous solution
$$E(x):=
\chi_{\Omega_{1}}(x)
\begin{bmatrix}
\sin(2\partiali x_1)+2\partiali\cos(2\partiali x_1)(x_1-x_2)\\
\sin\big((x_1-x_2)^2(x_1-1)^2x_2\big)-\sin(2\partiali x_1)\\
0
\end{bmatrix}
+\xi(x)\begin{bmatrix}0\\0\\1\end{bmatrix},\quad
\xi(x):=\partialrod_{i=1}^{3}x_{i}^2(1-x_{i})^2.$$
Thus, we extended the discontinuous vector field of Example \rgen{}{}ef{ex4} 
by zero in the third component and added a smooth bubble in the third component.
Hence, $E\in\rgen{}{}\setminus\hgen{1}{}$ and $\rgen{}{}ot E\in\rgen{}{}$ with
$$\rgen{}{}ot E(x)=
\chi_{\Omega_{1}}(x)\big(2x_{2}(x_1-x_2)(x_1-1)(2x_{1}-x_{2}-1)\cos(2\partiali x_1)\big)
\begin{bmatrix}0\\0\\1\end{bmatrix}
+\begin{bmatrix}\partial_{2}\xi\\-\partial_{1}\xi\\0\end{bmatrix}(x).$$
Note that even $\rgen{}{}ot E\in\hgen{1}{}$ holds.
We set zero Neumann boundary conditions on the whole boundary, i.e.,
$\Gammad=\emptyset$ and $\Gamman=\Gamma$. 
We calculated the approximation globally by solving 
the primal and dual problem with linear N\'ed\'elec elements.
This finite element approximation pair will be denoted by $(E_{\mathtt{h}},H_{\mathtt{h}})$.
The resulting linear systems were solved directly.
The approximations were calculated in uniformly refined regular meshes, 
where the jumps in the exact solution and in the right hand side $J$ coincide with element boundaries.
For each mesh we calculated the exact combined error and the majorant $\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})$. 
The results are displayed in Table \rgen{}{}ef{tbl:Ex5}.
\end{ex}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex5} (3D-eddy-current)}
\begin{center}
\begin{tabular}{r|c|c|c}
$N_{\mathtt{elem}}$ & $\tnorm{(E,H)-(E_{\mathtt{h}},H_{\mathtt{h}})}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}$ & difference $\dgen{}{}elta$ \\
\hline\hline
384    & 0.7228185218 & 0.7228185218 & 3.330669074e-16 \\
3072   & 0.3717887807 & 0.3717887807 & 6.106226635e-16 \\
24576  & 0.1883612515 & 0.1883612515 & 2.775557562e-16 \\
196608 & 0.0945757836 & 0.0945757836 & 8.604228441e-16
\end{tabular}
\end{center}
\label{tbl:Ex5}
\end{table}

\begin{ex}
\label{ex6}
We take the problem data of Example \rgen{}{}ef{ex4} and solve the primal and dual problems 
in adaptively refined meshes with linear N\'ed\'elec elements and linear Courant elements, respectively. 
This finite element approximation pair will be denoted by $(E_{\mathtt{h}},H_{\mathtt{h}})$ and the linear systems are solved directly. 
We compare optimal refinement achieved by using the exact error distribution $e_T$ 
to the refinement provided by the distribution of the majorant $\eta_T$, where 
\begin{align*}
e_T^2
&:=\tnorm{(E,H)-(E_{\mathtt{h}},H_{\mathtt{h}})}^2_T
:=\norm{E-E_{\mathtt{h}}}_{\rgen{}{}(T)}^2
+\norm{H-H_{\mathtt{h}}}_{\hgen{1}{}(T)}^2,\\
\eta_T^2
&:=\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})_{T}
:=\norm{J-E_{\mathtt{h}}-\nabla^{\bot}H_{\mathtt{h}}}_{\lgen{2}{}(T)}^2
+\norm{H_{\mathtt{h}}-\rgen{}{}otE_{\mathtt{h}}}_{\lgen{2}{}(T)}^2
\end{align*}
and $T$ denotes an element (triangle) of the mesh discretization. 
We start from a regular mesh with $200$ elements, and perform nine refinement iterations, 
where on each iteration $30\%$ of elements with the highest amount of error are refined. 
The refinement of element meshes is done by regular refinement such that
the resulting mesh does not contain hanging nodes. 
The results of Figure \rgen{}{}ef{fig:ex6} show that even though the equality is \emph{global}, 
the majorant can still be used to perform reliable adaptive computations. 
We see from Table \rgen{}{}ef{tbl:Ex6} that the number of elements in the optimal meshes and
the meshes produced using $\eta_T$ are very close to each other.
In Figure \rgen{}{}ef{fig:ex6mesh} we have depicted the meshes after the fourth refinement.
Figure \rgen{}{}ef{fig:ex6mesh_final} depicts one of the finest parts of the final meshes.
In fact, the adaptive refinement using $\eta_T$ is very close to optimal in each step, 
and the resulting approximation after the last refinement is practically the same.
\end{ex}

\begin{figure}\label{fig:ex6}
\end{figure}

\begin{figure}\label{fig:ex6mesh}
\end{figure}

\begin{figure}\label{fig:ex6mesh_final}
\end{figure}

\begin{table} [htdp]
c_{\A}ption{Adaptive computation of Example \rgen{}{}ef{ex6}. 
The number of elements in the optimal meshes and the meshes generated by the help of $\eta_T$.}
\begin{center}
\begin{tabular}{r|r|r|r|l}
Ref. & optimal & with $\eta_T$ & difference & difference $\%$ \\
\hline\hline
- &    200 &    200 &    0 & 0 \\
1 &    434 &    434 &    0 & 0 \\
2 &    998 &   1002 &    4 & 0.40 \\
3 &   2240 &   2252 &   12 & 0.53 \\
4 &   4823 &   4878 &   55 & 1.14 \\
5 &  10378 &  10446 &  68  & 0.65 \\
6 &  22116 &  22337 &  221 & 0.99 \\
7 &  46388 &  46768 &  380 & 0.81 \\
8 &  96859 &  97832 &  973 & 1.00 \\
9 & 198704 & 200970 & 2266 & 1.14
\end{tabular}
\end{center}
\label{tbl:Ex6}
\end{table}

\begin{ex}
\label{ex7}
We take the 2D-eddy-current problem of Section \rgen{}{}ef{subsec:EC2D} in the $L$-shaped domain 
$\Omega:=(0,1)^2\setminus\big([\nicefrac{1}{2},1]\times[0,\nicefrac{1}{2}]\big)$ 
with $\epsilon=\id$, $\mu=1000$ and $J=[1,0]^\top$. 
We set zero Dirichlet boundary conditions on the whole boundary, i.e., $\Gammad=\Gamma$ and $\Gamman=\emptyset$. 
The exact solution of this problem is unknown. 
However, since the majorant gives indeed the exact error in the combined norm, 
we will use this information in this example. Therefore, all the error values 
in Figure \rgen{}{}ef{fig:ex7} and Table \rgen{}{}ef{tbl:Ex7} are actually the values of the majorant. 
We compare uniform refinement and adaptive refinement using $\eta_T$ with
$$\eta_T^2
=\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})_{T}
=\norm{J-E_{\mathtt{h}}-\nabla^{\bot}H_{\mathtt{h}}}_{\lgen{2}{}(T)}^2
+\norm{H_{\mathtt{h}}-\mu^{-1}\rgen{}{}otE_{\mathtt{h}}}_{\lgen{2}{}(T),\mu}^2,$$
refining $30\%$ of elements on each refinement iteration as before. 
We solve the primal and dual problems with linear N\'ed\'elec elements and linear Courant elements, respectively.
The resulting linear systems are solved directly.
We see from Figure \rgen{}{}ef{fig:ex7} that the adaptive procedure is beneficial in this example.
We have also depicted the approximation in Figure \rgen{}{}ef{fig:ex7approx} 
and the mesh in Figure \rgen{}{}ef{fig:ex7mesh} after the fifth refinement.
It can be concluded that in addition to providing the exact error, 
the majorant also provides a good error indicator without any additional computational expenditures.
\end{ex}

\begin{figure}\label{fig:ex7}
\end{figure}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex7} (2D-eddy-current) 
Adaptively refined meshes.}
\begin{center}
\begin{tabular}{r|c|c}
$N_{\mathtt{elem}}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}/\norm{J}_{\lgen{2}{}}$ \\
\hline\hline
96     & 0.2534 & 0.2926 \\
230    & 0.1534 & 0.1771 \\
541    & 0.0842 & 0.0973 \\
1204   & 0.0467 & 0.0539 \\
2623   & 0.0309 & 0.0357 \\
6082   & 0.0203 & 0.0234 \\
13514  & 0.0135 & 0.0155 \\
29530  & 0.0093 & 0.0107 \\
63363  & 0.0062 & 0.0072 \\
134205 & 0.0043 & 0.0050
\end{tabular}
\end{center}
\label{tbl:Ex7}
\end{table}

\begin{figure}\label{fig:ex7approx}
\end{figure}

\begin{figure}\label{fig:ex7mesh}
\end{figure}

\begin{ex}
\label{ex8}
We take the 2D-eddy-current problem of Section \rgen{}{}ef{subsec:EC2D} in $\Omega:=(0,1)^2$.
In order to define discontinuous data, we define with $\xi(x):=\ln(2+x_2)$ and
\begin{align*}
\Omega_1
&:=\big((0,1)\times(0.4,0.6)\big)\cup\big((0.3,0.5)\times(0,1)\big),&
\epsilon|_{\Omega_1}
&:=\id,&
\epsilon|_{\Omega\setminus\overline{\Omega}_1}
&:=100\cdot\id,\\
& &
\mu|_{\Omega_1}
&:=1000,&
\mu|_{\Omega\setminus\overline{\Omega}_1}
&:=1,\\
\Omega_2
&:=(0,1)\times(0.35,0.65),&
J|_{\Omega_2}
&:=\xi\begin{bmatrix}1\\0\end{bmatrix},&
J|_{\Omega\setminus\overline{\Omega}_2}
&:=-\xi\begin{bmatrix}0\\1\end{bmatrix}.
\end{align*}
We set zero Dirichlet boundary conditions on the right side of the rectanglular boundary
and zero Neumann boundary condition on the remaining part, i.e., 
$\Gammad=\set{x\in\Omega}{x_1=1}$. 
As in Example \rgen{}{}ef{ex7}, the exact solution of this problem is unknown, so the error values 
in Figure \rgen{}{}ef{fig:ex8} and Table \rgen{}{}ef{tbl:Ex8} are the values of the majorant. 
We compare uniform refinement and adaptive refinement using $\eta_T$ with
$$\eta_T^2
=\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})_{T}
=\norm{J-\epsilonE_{\mathtt{h}}-\nabla^{\bot}H_{\mathtt{h}}}_{\lgen{2}{}(T),\epsilon^{-1}}^2
+\norm{H_{\mathtt{h}}-\mu^{-1}\rgen{}{}otE_{\mathtt{h}}}_{\lgen{2}{}(T),\mu}^2,$$
refining $30\%$ of elements on each refinement iteration as before. 
We solve the primal and dual problems with linear N\'ed\'elec elements and linear Courant elements, respectively.
The resulting linear systems are solved directly.
Again, we see from Figure \rgen{}{}ef{fig:ex8} that the adaptive procedure is beneficial in this example.
We have also depicted the approximation in Figure \rgen{}{}ef{fig:ex8approx} 
and the mesh in Figure \rgen{}{}ef{fig:ex8mesh} after the third refinement.
\end{ex}

\begin{figure}\label{fig:ex8}
\end{figure}

\begin{table} [htdp]
c_{\A}ption{Example \rgen{}{}ef{ex8} (2D-eddy-current) 
Adaptively refined meshes.}
\begin{center}
\begin{tabular}{r|c|c}
$N_{\mathtt{elem}}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}$ & $\sqrts{\mathcal{M}ec(E_{\mathtt{h}},H_{\mathtt{h}})}/\norm{J}_{\lgen{2}{},\epsilon^{-1}}$ \\
\hline\hline
800    & 0.1632 & 0.2941 \\
1827   & 0.0921 & 0.1659 \\
4367   & 0.0513 & 0.0924 \\
10214  & 0.0307 & 0.0554 \\
23657  & 0.0199 & 0.0359 \\
51429  & 0.0128 & 0.0231 \\
113073 & 0.0085 & 0.0153
\end{tabular}
\end{center}
\label{tbl:Ex8}
\end{table}

\begin{figure}\label{fig:ex8approx}
\end{figure}

\begin{figure}\label{fig:ex8mesh}
\end{figure}

To conclude, in all the tests performed, nonzero values of $\dgen{}{}elta$ were of magnitude $10^{-18}$-$10^{-15}$. 
This is within the limit of machine precision, 
so numerically these numbers are considered zero. 
In addition to verifying the equality, 
we also performed three simple examples to show that the majorant 
can be used to perform refinement of element meshes
without any additional computational expenditures.

\begin{acknow}{10}
We are deeply indebted to Sergey Repin for so many interesting
and encouraging discussions. We thank him gratefully not only 
for being an academic colleague and teacher but also a good friend.

The first author also thanks Jan Valdman for showing how to program vectorized FEM solvers.
The first author is funded by the Finnish foundations KAUTE Foundation and
V\"ais\"al\"a Foundation of the Finnish Academy of Science and Letters.

This contribution has been worked out mainly while the first author 
was visiting the Fakult\"at f\"ur Mathematik of the Universit\"at Duisburg-Essen during 2013.
\end{acknow}

 

\end{document}