Mixed Finite Element Methods for the Poisson Equation Using Biorthogonal and Quasi-Biorthogonal Systems

We introduce two three-field mixed formulations for the Poisson equation and propose finite element methods for their approximation. Both mixed formulations are obtained by introducing a weak equation for the gradient of the solution by means of a Lagrange multiplier space. Two efficient numerical schemes are proposed based on using a pair of bases for the gradient of the solution and the Lagrange multiplier space forming biorthogonal and quasi-biorthogonal systems, respectively. We also establish an optimal a priori error estimate for both finite element approximations.


Introduction
In many practical situations, it is important to compute dual variables of partial differential equations more accurately.For example, the gradient of the solution is the dual variable in case of the Poisson equation, whereas the stress or pressure variable is the dual variable in case of elasticity equation.Working with the standard finite element approach these variables should be obtained a posteriori by differentiation, which will result in a loss of accuracy.In these situations, a mixed method is often preferred as these variables can be directly computed using a mixed method.
In this paper, we introduce two mixed finite element methods for the Poisson equation using biorthogonal or quasi-biorthogonal systems.Both formulations are obtained by introducing the gradient of the solution of Poisson equation as a new unknown and writing an additional variational equation in terms of a Lagrange multiplier.This gives rise to two additional vector unknowns: the gradient of the solution and the Lagrange multiplier.In order to obtain an efficient numerical scheme, we carefully choose a pair of bases for the space of the gradient of the solution and the Lagrange multiplier space in the discrete setting.Choosing the pair of bases forming a biorthogonal or quasi-biorthogonal system for these two spaces, we can eliminate the degrees of freedom associated with the gradient of the solution and the Lagrange multiplier and arrive at a positive definite formulation.The positive definite formulation involves only the degrees of freedom associated with the solution of the Poisson equation.Hence a reduced system is obtained, which is easy to solve.The first formulation is discretized by using a quasi-biorthogonal system, whereas the second one, which is a stabilized version of the first one, is discretized using a biorthogonal system.
There are many mixed finite element methods for the Poisson equation [1][2][3][4][5][6][7][8].However, all of them are based on the two-field formulation of the Poisson equation and hence are not amenable to the application of the biorthogonal and quasi-biorthogonal systems.We need a three-field formulation to apply the biorthogonal and quasi-biorthogonal systems which leads to a symmetric formulation (see [9] for a three-field formulation in linear elasticity).The use of biorthogonal and quasi-biorthogonal systems allows an easy static condensation of the auxiliary variables (gradient of the solution and Lagrange multiplier) leading to a reduced linear system.These variables can be recovered just by inverting a diagonal matrix.Therefore, in this paper we present three-field formulations of the Poisson equation to apply the biorthogonal and quasi-biorthogonal systems.
The structure of the rest of the paper is as follows.In Section 2, we introduce our three-field formulations and

Two Three-Field Formulations of the Poisson Equation
In this section we introduce two three-field mixed formulations of the Poisson problem.Let Ω ⊂ R  ,  ∈ {2, 3}, be a bounded convex domain with polygonal or polyhedral boundary Ω with the outward pointing normal n on Ω.
Our new formulation is obtained by introducing an auxiliary variable  = ∇ such that the minimization problem (4) is rewritten as the following constrained minimization problem: arg min where the natural norm for the product space  ×  is √‖‖ 2 1,Ω + ‖‖ 2 0,Ω for (, ) ∈  × .We write a weak variational equation for equation  = ∇ in terms of the Lagrange multiplier space  to obtain the saddle-point problem of the minimization problem (8).The saddle-point formulation is to find (, , ) ∈  ×  ×  such that where A similar three-field formulation-called the Hu-Washizu formulation-is very popular in designing finite element methods to alleviate locking effect in linear elasticity [11][12][13].
In order to show that the saddle-point problem (9) has a unique solution, we want to apply a standard saddlepoint theory [3,4,10].To this end, we need to show the following three conditions of well-posedness.
The linear form ℓ(⋅) and the two bilinear forms (⋅, ⋅) and (⋅, ⋅) are continuous by the Cauchy-Schwarz inequality.The coercivity of the bilinear form (⋅, ⋅) on the kernel space  follows from Poincaré inequality: It remains to show that the bilinear form (⋅, ⋅) satisfies the inf-sup condition.
Lemma 1.The bilinear form (⋅, ⋅) satisfies the inf-sup condition; that is, Proof.Using V = 0 in the expression on the left above, we have sup To summarize we have proved the following theorem.
Theorem 2. The saddle-point problem (9) has a unique solution (, , ) ∈  ×  ×  and The main difficulty in this formulation is that the bilinear form (⋅, ⋅) is not coercive on the whole space  × .It is only coercive on the subspace  of  × .One has to consider this difficulty in developing a finite element method for this problem.One approach to make the bilinear form (⋅, ⋅) coercive on the whole space  ×  is to stabilize it.There are many approaches to stabilize this formulation (see [2,5,6]).Here we modify the bilinear form (⋅, ⋅) to get a consistent stabilization as proposed in [14] for the Mindlin-Reissner plate so that it is now coercive on the whole space  × .This is obtained by adding the term to the bilinear form (⋅, ⋅).Thus our problem is to find where Theorem 3. The saddle-point problems (9) and (18) yield the same solution.
Proof.By construction the solution to ( 18) is also the solution to (9).Moreover, as the second equation of ( 18) yields  = ∇, we can substitute this into the first equation and set V = 0 to get  = −∇.This proves that the solution to ( 18) is also the solution to (9).
Here the bilinear form (⋅, ⋅) is coercive on the whole space  ×  from the triangle and Poincaré inequalities.

Lemma 4. The bilinear form
Proof.We start with the triangle inequality From Poincaré inequality, there exists a constant  > 0 such that and hence we have Moreover, the other conditions of well-posedness for this new saddle-point formulation (18) can also be shown exactly as for (9).
In the following, we present finite element approximations for both proposed three-field formulations.In the first step, we propose a finite element method for the nonstabilized formulation (9), which is based on quasi-biorthogonal systems.In the second step, we present a finite element method for the stabilized formulation (18), which is based on biorthogonal systems.The first method works only for simplicial elements, whereas the second method works for meshes of parallelotopes and simplices.

Finite Element Approximation and a Priori Error Estimate
Let T ℎ be a quasi-uniform partition of the domain Ω in simplices or parallelotopes having the mesh size ℎ.Let T be a reference simplex or -cube, where the reference simplex is defined as and a -cube T := (−1, 1)  .Note that for  = 2 a simplex is a triangle, and for  = 2 it is a tetrahedron.Similarly, a 2-cube is a square, whereas, 3-cube is a cube.The finite element space is defined by affine maps   from a reference element T to a physical element  ∈ T ℎ .Let Q( T) be the space of bilinear/trilinear polynomials in T if T is a reference square/cube or the space of linear polynomials in T if T is a reference simplex.Then the finite element space based on the mesh T ℎ is defined as the space of continuous functions whose restrictions to an element  are obtained by maps of bilinear, trilinear, or linear functions from the reference element; that is, Advances in Numerical Analysis (See [3,10,15]).We start with some abstract assumptions for two discrete spaces [ ℎ ]  ⊂  and [ ℎ ]  ⊂  to discretize the gradient of the solution and the Lagrange multiplier spaces, where  ℎ and  ℎ both are piecewise polynomial spaces with respect to the mesh T ℎ and are both subsets of  2 (Ω).Explicit construction of local basis functions for these spaces for the mixed finite element method ( 9) is given in Section 3.1 and for the mixed finite element method ( 18) is given in Section 3.2.We assume that  ℎ and  ℎ satisfy the following assumptions (see also [16,17]).
(ii) There is a constant  > 0 independent of the mesh T ℎ such that (iii) The space  ℎ has the approximation property: inf (iv) The space  ℎ has the approximation property: inf Now we define biorthogonality and quasi-biorthogonality (see [16,18,19]).

Definition 5. The pair of bases {𝜇
where D 1 and D 2 are diagonal matrices and R is a sparse rectangular matrix.
We recall that a Gram matrix G of two sets of basis functions {  } 1≤≤ and {  } 1≤≤ is the matrix G whose th entry is It is worth noting that the quasi-diagonal matrix is inverted by inverting two diagonal matrices and scaling the rectangular matrix.

Finite Element
Approximation for the Nonstabilized Formulation (9).We now turn our attention to a finite element approximation for the nonstabilized formulation (9).This approximation works only for simplicial meshes and is based on a quasi-biorthogonal system as in [19].Let be the space of bubble functions, where  is a constant and {   } +1 =1 is the set of barycentric coordinates on .Let  ℎ =  ℎ ⊕  ℎ .Explicit construction of basis functions of  ℎ and  ℎ on a reference element is provided below.
The finite element problem is to find We note that the standard finite element space  ℎ is enriched with elementwise defined bubble functions to obtain the space  ℎ , which is done to obtain the coercivity of the bilinear form (⋅, ⋅) on the kernel space  ℎ defined as (see Lemma 7) To simplify the numerical analysis of the finite element problem, we introduce a quasiprojection operator: This type of operator was introduced in [9] to obtain the finite element interpolation of nonsmooth functions satisfying boundary conditions and is used in [20] in the context of mortar finite elements.The definition of  ℎ allows us to write the weak gradient as where operator  ℎ is applied to the vector ∇ ℎ componentwise.We see that  ℎ is well defined due to Assumption 1(i) and (ii).Furthermore, the restriction of  ℎ to  ℎ is the identity.Hence  ℎ is a projection onto the space  ℎ .We note that  ℎ is not the orthogonal projection onto  ℎ but an oblique projection onto  ℎ see [21,22].The stability and approximation properties of  ℎ in  2 and  1 -norm can be shown as in [16,23].In the following, we will use a generic constant , which will take different values at different places but will be always independent of the mesh size ℎ.

Lemma 6. Under Assumption 1(i)-(iii)
and for 0 <  ≤ 1 and We need to construct another finite element space  ℎ ⊂  2 (Ω) that satisfies Assumption 1(i)-(iii) and also leads to an efficient numerical scheme.The main goal is to choose the basis functions for  ℎ and  ℎ so that the matrix D associated with the bilinear form ∫ Ω  ℎ :  ℎ  is easy to invert.To this end, we consider the algebraic form of the weak equation for the gradient of the solution, which is given by where ⃗  and ⃗  are the solution vectors for  ℎ and  ℎ and D is the Gram matrix between the bases of [ ℎ ]  and [ ℎ ]  .Hence if the matrix D is diagonal or quasi-diagonal, we can easily eliminate the degrees of freedom corresponding to  ℎ and  ℎ .This then leads to a formulation involving only one unknown  ℎ .Statically condensing out variables  ℎ and  ℎ we arrive at the variational formulation of the reduced system given by (68) (see Section 3.3).Note that the matrix D is the Gram matrix between the bases of [ ℎ ]  and [ ℎ ]  .
We cannot use a biorthogonal system in this case as the stability will be lost (see the Appendix).This motivates our interest in a quasi-biorthogonal system.We now show the construction of the local basis functions for  ℎ and  ℎ so that they form a quasi-biorthogonal system.The construction is shown on the reference triangle T = {(, ) : 0 ≤ , 0 ≤ , + ≤ 1} in the two-dimensional case and on the reference tetrahedron T = {(, , ) : 0 ≤ , 0 ≤ , 0 ≤ ,  +  +  ≤ 1} in the three-dimensional case.
Similarly, on the reference tetrahedron, the local basis functions { φ1 , . . ., φ5 } of  ℎ defined as are quasi-biorthogonal.These five basis functions { φ1 , . . ., φ5 } of  ℎ or { μ1 , . . ., μ5 } of  ℎ are associated with four vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1) and the barycenter (1/4, 1/4, and 1/4) of the reference tetrahedron.Using the finite element basis functions for  ℎ , the finite element basis functions for  ℎ are constructed in such a way that the Gram matrix on the reference element T is given by for the two-dimensional case and for the three-dimensional case.After ordering the degrees of freedom in  ℎ and  ℎ so that the degrees of freedom associated with the barycenters of elements come after the degrees of freedom associated with the vertices of elements, the global Gram matrix D is quasi-diagonal.
Thus in the two-dimensional case, the local basis functions μ1 , μ2 , and μ3 are associated with the vertices of the reference triangle, and the function μ4 is associated with the barycenter of the triangle and defined elementwise.That means μ4 is supported only on the reference element.The situation is similar in the three-dimensional case.Hence Advances in Numerical Analysis using a proper ordering, the support of a global basis function   of  ℎ is the same as that of the global basis function   of  ℎ for  = 1, . . ., .However, the global basis functions  1 , . . .,   of  ℎ are not continuous on interelement boundaries.Now we show that  ℎ satisfies Assumption 1(i)-(iii).As the first assumption is satisfied by construction, we consider the second assumption.Let  ℎ = ∑  =1     ∈  ℎ and set  ℎ = ∑  =1     ∈  ℎ .By using the quasibiorthogonality relation (1) and the quasiuniformity assumption, we get where ℎ  denotes the meshsize at th vertex.Taking into account the fact that ‖ ℎ ‖ 2 0,Ω ≡ ‖ ℎ ‖ 2 0,Ω ≡ ∑  =1  2  ℎ   , we find that Assumption 1(ii) is satisfied.Since the sum of the local basis functions of  ℎ is one, Assumption 1(iii) can be proved as in [23,25].Assumption 1(iv) follows by standard arguments.
There are ( + 2) local basis functions of  ℎ , where ( + 1) of them are associated with vertices and one is associated with the element barycenter.Then the local basis functions of  ℎ are also divided into two parts: ( + 2) basis functions associated with the vertices and 1 basis function associated with the element.Let   be the local basis function associated with the element barycenter of  ∈ T ℎ for  ℎ .Let  ℎ be the space of piecewise constant functions with respect to T ℎ , and  ℎ :  2 (Ω) →  ℎ is defined as Then  ℎ is a projection operator onto  ℎ .Moreover, and this implies that This yields the following coercivity result.
Lemma 7. The bilinear form (⋅, ⋅) is coercive on the kernel space  ℎ .In fact, One needs to modify the finite element method for these meshes.
The inf-sup condition follows exactly as in the continuous setting using Assumption 1(ii).Lemma 9.There exists a constant β > 0 such that sup Proof.The proof follows exactly as in the continuous setting.We use V ℎ = 0 in the expression sup to get sup Assumption 1(ii) then yields the final result.
Hence we have the following stability and approximation result.
Owing to Assumption 1(iii)-(iv) the discrete solution converges optimally to the continuous solution.(18).After stabilizing formulation (9), the bilinear form (⋅, ⋅) (18) is coercive on the whole product space  × .Therefore, we can use the standard finite element space  ℎ to discretize the gradient of the solution for the saddle-point problem (18) so that  ℎ =  ℎ .Thus our discrete saddle-point formulation of ( 18) is to find

Finite Element Approximation for the Stabilized Formulation
We can now verify the three conditions of well-posedness for the discrete saddle point problem.The continuity of (⋅, ⋅), (⋅, ⋅), and ℓ(⋅) follows as in the continuous setting, and coercivity is immediate from Lemma 4 as  ℎ × [ ℎ ]  ⊂  × .
Explicit construction of basis functions of  ℎ on a reference element will be given below.Assumption 1(ii) guarantees that the bilinear form (⋅, ⋅) satisfies the inf-sup condition as in the case of the previous finite element method.
We have thus proved the following theorem from the standard theory of saddle-point problem (see [4]).
The optimal convergence of the discrete solution is then guaranteed by the standard approximation property of  ℎ and  ℎ and Assumption 1(iii).
In the following, we give these basis functions for linear simplicial finite elements in two and three dimensions [18].The parallelotope case is handled by using the tensor product construction of the one-dimensional basis functions presented in [20,26].These basis functions are very useful in the mortar finite element method and extensively studied in one-and two-dimensional cases [16,23,26].
Note that we have imposed the condition dim  ℎ = dim  ℎ to get that our Gram matrix is a square matrix.Therefore, the local basis functions for linear/bilinear/trilinear finite spaces are associated with the vertices of triangles, tetrahedra, quadrilaterals, and hexahedra.For example, for the reference triangle T := {(, ) : 0 < , 0 < ,  +  < 1}, we have where the basis functions μ1 , μ2 and μ3 are associated with three vertices (0, 0), (1, 0), and (0, 1) of the reference triangle.
The global basis functions for the test space are constructed by gluing the local basis functions together following exactly the same procedure of constructing global finite element basis functions from the local ones.These global basis functions then satisfy the condition of biorthogonality with global finite element basis functions of  ℎ .These basis functions also satisfy Assumption 1(i)-(iv) (see [17,18]).

Static Condensation.
Here we want to eliminate  ℎ and Lagrange multiplier  ℎ from the saddle-point problem (57).We make use of the operator  ℎ defined previously.Using the biorthogonality relation between the basis functions of  ℎ and  ℎ , the action of operator  ℎ on a function V ∈  2 (Ω) can be written as which tells that the operator  ℎ is local in the sense to be given below (see also [27]).
Using the property of operator  ℎ , we can eliminate the degrees of freedom corresponding to  ℎ and  ℎ so that our problem is to find  ℎ ∈  ℎ such that  ( ℎ ) = min where  (V ℎ ) =     ∇ ( ℎ (∇V ℎ )) for the stabilized formulation (18).Let the bilinear form Ã(⋅, ⋅) be defined as for the non-stabilized formulation (9) and