Mixed formulation of Laplacian

The following deals with mixed problem of Dirichelet kind for Poisson equation. To get started, we cast Poisson equation into first order equation.

Recall, original Poisson equation is
\begin{equation}
  -u^{''}(x) = f(x)
\end{equation}
Let's set $\sigma = -\text{grad}\; u$. $grad$ - aka Gradient is a vector field. In its full glory written as
\begin{equation}
  \nabla u = (\frac{\partial u}{\partial x},\frac{\partial u}{\partial y}, \frac{\partial u}{\partial z} )
\end{equation}
clearly $\text{grad}\; u$ is a vector field. For our case, we simply set this to $\sigma= -\text{grad}\; u$.
Our original equation is $-u^{''}(x) = f(x)$. To express this in terms of $\sigma$, we use another operator - div.
Let's write down what div is.
The divergence of a vector field, say $F=(F_x,F_y,F_z)$, is a scalar function that represents the net rate of outward flux per unit volume at each point in the field. It gives a measure of how much the vector field is spreading out or compressing at a given point. The divergence is calculated as follows:
\begin{equation}
  \text{div} F = \frac{\partial F_x}{\partial x}+ \frac{\partial F_x}{\partial y} + \frac{\partial F_x}{\partial z}
\end{equation}
Clear from above definition if we $\text{div}\; \text{grad}\; u$, we get
\begin{equation}
  \text{ div } \text{grad}\; u = \frac{\partial^2 u}{\partial x^2}+ \frac{\partial^2 u}{\partial y^2}+ \frac{\partial^2 u}{\partial z^2}
\end{equation}
which in our original equation is same as $u^{''}$.
Hence, the original equation becomes,
\begin{equation}
  \text{div}\; \text{grad}\; u(x) = f(x)
\end{equation}
This in turn is shortened as,

\begin{equation}
  \sigma = - \text{grad}\; u, \text{div}\; \sigma = f
\end{equation}
The pair $(\sigma,u)$ can be characterized as critical point (unique) of the functional
\begin{equation}
  I(\sigma,u) = \int_\Omega (\frac{1}{2} \sigma.\sigma - u \text{div}  \sigma)dx + \int_\Omega fu dx
\end{equation}
over $H(div:\Omega) \times L^2(\Omega)$ where $H(div:\Omega) = {\sigma \in L^2; div \sigma \in L^2}$.
Equivalently one can solve weak problem
\begin{equation}
  \int_{\Omega}\tau.\sigma dx - \int_{\Omega}u \text{div} \tau dx = 0, \tau \in H(div:\Omega)
\end{equation}
\begin{equation}
  \int_{\Omega} \text{div} \tau v dx = \int_{\Omega} f v dx
\end{equation}
This fits into abstract framework if we define, $V=H(div:\Omega) \times L^(\Omega)$.
\begin{equation}
  B(\sigma,u;\tau,v) = \int_{\Omega}\sigma.\tau dx - \int_{\Omega}u \text{div} \tau dx + \int_{\Omega}\text{div} \sigma v dx, F(\tau,v) = \int_{\Omega}fvdx
\end{equation}


In this case the bilinear form $B$ is not coercive, and so the choice of subspaces an the analysis is not so simple
as for the standard finite element method for Poisson’s equation. Finite element discretizations based on such saddle point variational principles are called mixed finite element methods. Thus a mixed finite element for Poisson’s equation is obtained by choosing subspaces $\Sigma_h \subset H(div;\Omega)$ and $V_h \subset L_2(\Omega)$ and seeking a critical point of $I$ over $\Sigma_h \times V_h$. The resulting Galerkin method has the form: Find $\sigma_h \in \Sigma_h,u_h \in V_h$ satisfying
\begin{equation}
      \int_{\Omega} \sigma_h \cdot \tau \, dx - \int_{\Omega} u_h \, \text{div} \tau \, dx = 0, \quad \forall \tau \in \Sigma_h,
    \int_{\Omega} \text{div} \sigma_h v \, dx = \int_{\Omega} fv \, dx, \quad \forall v \in V_h.
  \end{equation}
Since the bilinear form is not coercive, it is not automatic that the linear system is nonsingular.
If $f=0$, then $\int_{\Omega} fv \, dx=0$. This means $\int_{\Omega} \text{div} \sigma_h v \, dx = 0$. This foces $\sigma_h=0$ as $v$ is trail function. This, in turn, via $\int_{\Omega} \sigma_h \cdot \tau \, dx - \int_{\Omega} u_h \, \text{div} \tau \, dx = 0$ forces $\int_{\Omega} u_h \, \text{div} \tau \, dx = 0$. Then $u_h=0$.
Choosing $\tau = \sigma h$ and $v = uh$ and adding the discretized variational equations, it follows
immediately that when $f = 0, \sigma_h = 0$. However, $u_h$ need not vanish.