#### Abstract

We consider a Poisson boundary value problem and its functional a posteriori error estimate derived by S. Repin in 1999. The estimate majorizes the seminorm of the error of the discrete solution computed by FEM method and contains a free ux variable from the (div) space. In order to keep the estimate sharp, a procedure for the minimization of the majorant term with respect to the ux variable is introduced, computing the free ux variable from a global linear system of equations. Since the linear system is symmetric and positive definite, few iterations of a conjugate gradient method with a geometrical multigrid preconditioner are applied. Numerical techniques are demonstated on one benchmark example with a smooth solution on a unit square domain including the computation of the approximate value of the constant in Friedrichs' inequality.

#### 1. Introduction

A priori rate convergence estimates for finite element approximations of elliptic boundary value problems have been investigated in the 70 s–80 s (see, e.g., [1]). However, adaptive multilevel algorithms require a posteriori estimates able to (a) provide a reliable and directly computable estimate of the approximation error, and (b) efficient error indicator that detects the regions with excessively high errors.

In the recent decades, a posteriori error estimates for linear elliptic and parabolic problems were intensively investigated. A reader will find a systematic exposition of the main approaches to a posteriori error estimation of finite element approximations (such as residual or gradient averaging methods) in the monographs [2–5] and papers [6–8] and in literature cited therein.

In this paper, a posteriori estimates that majorate the difference between exact solution of a linear elliptic equation and any function in the admissible (energy) class are studied. For the class of uniformly convex variational problems computable error majorants (for any conforming approximation) were derived by the variational techniques in the mid 90 s using duality theory of the calculus of variations. Key publications related to this subject are in [9–11]. Another “nonvariational" method was introduced in [12]. In this paper it was stated that for linear elliptic problems both methods lead to the same error majorants. Later it was applied to many problems, including parabolic equations and nonliner problems [13–17].

As an example of demonstration, let us consider a scalar boundary value (Poisson's) problem for a searched function from a Sobolev space . The right hand side and an open bounded domain , where denotes a domain dimension, that is, are given. Assume that is an approximation of . Then, a functional error estimate from [9] holds for all functions , denotes the norm. Note that no mesh-dependent constants or any assumptions on regularity of an exact solution are contained in this estimate. The only global constant included represents a constant from the Friedrichs' inequality which holds for all . Thus the constant depends on the domain only and can be precomputed (it is demonstated in Subsection 3.1). For a given , the quality of the estimate (1.2) is measured by a ratio of its right hand and left hand side known as an index of efficiency. It obviously holds and the equality is valid for the choice , that is, if is chosen as the flux of an exact solution . Having this interpretation of in mind, there are known ways [18] how to compute a reasonable flux approximation from the discrete solution .

(1)Averaging on the mesh of the discrete solution . In this case, the flux approximation is computed as , where represents an averaging gradient operator, see, for example, [8] for more details. This is a cheap method providing some preliminary knowledge on the upper bound of the error.(2)Averaging on a refined mesh. It is similar to (1), only with the difference that the averaging is done for the solution calculated on once more (or more times) refined mesh. This method can be regarded as a quantitative form of the Runge's rule. It is more expensive, but provides (generally not always) sharper results.(3)Using partially equilibrated fluxes. By postprocessing of , a function is constructed such that and are sufficiently close to in norm. Then, the substitution of into (1.2) provides an estimate where is arbitrary. (4)Minimization of the right hand side of (1.2) with respect to free variable on the mesh of the discrete solution . This is the most expensive method for a detailed knowledge of the error.A comparison of methods (1), (2), and (4) for a class of problems with nonlinear boundary conditions can be found, for example, in [17]. We should mention that there are many advanced forms of error bounds for the Poisson's equation (1.1). They are based, for example, on decomposition of the domain or on partial equilibrated fluxes which compute as a solution of small local problems. For more information see [18, section 3.5] and papers [19–22].

This paper focuses merely on the method (4), that is, the minimization of the right hand side of (1.2). We are interested in fast computation of the reliable estimate rather than the indication of regions with high error (adaptivity). The main motivation is to provide people working with “functional a posteriori estimates" concepts and software to speed up their computation and go for larger size problems. In particular, the paper

(i)formulates a majorant minimization problem in (div) space and an algorithm for its computation on continuous level;(ii)demonstrates a numerical computation of the Friedrich's constant;(iii)discretizes the minimization algorithm and applies RT0 elements to obtain a linear system of equations. So far, only vector nodal linear elements were applied [23, 24];(iv)introduces a multigrid-preconditioned conjugate graduate method as an iterative solver for resulting linear system and demonstrates its optimality on one benchmark example.We would also like to attract attention of different groups and evoke cooperations. Since the majorant minimization problem (discretized by RT0 elements) is about 3 times larger than the Poisson problem (discretized by linear nodal elements), the other error estimates specially developed for linear problems might perform faster for the benchmark problem with a smooth solution (the exact comparison is not done here). However, functional a posteriori error estimates are the only tool to provide the guaranteed error of approximation of nonlinear problems [17, 25, 26].

#### 2. Majorant Minimization Problem

*Problem 2.1 (Minimization problem). *Given a discrete solution , a right-hand side of the Poisson problem , the Friedrich's constant belonging to the domain . Find a function satisfying the condition

In order to avoid complications with the nondifferentiability of norm terms in (2.1), we apply the inequality valid for all to obtain where the upper bound in (2.2) denotes a functional majorant The majorant arguments are known and and are free parameters. For a fixed choice of parameters , the majorant represents a quadratic functional in . On the other hand, for a fixed , the parameter minimizes amongst all positive . It suggests the following solution algorithm to Problem 2.1.

*Algorithm 2.2 (Majorant minimization algorithm). *Given .(a)Compute from the minimization of the quadratic problem
(b)Update from using (2.4). If the convergence in is not achieved then go to step (a).

For a detailed analysis of step (a) it is convenient to decompose in its -independent and -dependent parts as where Here, denotes the scalar product. Instead of considering (2.5) in Algorithm 2.2, is solved from the minimization problem for given . Since for , the minimal value of must be nonpositive. Besides, it holds and if and only if and .

The finite element method (FEM) is used for the discretization of the minimization problem (2.7). The domain is divided by a regular triangulation in triangles in the sense of Ciarlet [1], that is, is a finite partition of into closed triangles; two distinct elements and are either disjoint, or is a complete edge or a common node of both and . Let us assume a finite element basis in the space and a vector representing in this basis. Norm terms in the definition of are read then after the discretization where and represent the ”mass” and “div div” matrices. Similarly, linear functionals are discretized as with some vectors and . It allows to express a discrete form of in the provided finite element basis, The minimization of with respect to the vector leads for a given value to a linear system of equations for a minimizing vector . The linear system is represented by a matrix which is symmetric and positive definite. Symmetry is the result of symmetries of matrices and . The positive definiteness is proved by a contradiction: if there was a nonzero vector for which the quadratic form in (2.10) was not positive, we could construct a sequence of arguments , , in which the values of would converge (for ) to . Consequently, since is independent of , sequence of values of the majorant would also converge to , which is not possible due to the reliability of the estimate (2.2).

The next benchmark example explains application issues of the majorant estimate (1.2) and of Algorithm 2.2.

#### 3. Benchmark Example

Let us assume the right-hand side in the unit square domain . The square geometry is discretized using a sequence of nested uniform triangular meshes as displayed in Figure 1. A discrete solution is computed from by using FEM with the nodal linear (Courant) ansatz functions on each triangular mesh. The discrete solutions are displayed as the left column pictures of Figure 2. For this particular right-hand side, the exact solution of (1.1) reads for all and its flux is for all . Therefore, the exact solution of the Benchmark example is a polynomial and, thus, integration error can be avoided.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

##### 3.1. Friedrich's Constant and Its Computation

The computation of the functional majorant (2.3) requires the knowledge of the constant from the Friedrich's inequality (1.3). Under the assumption of the finite element basis we can represent the function as the vector , where . In the same basis, the bilinear forms read where and are stiffness and mass matrices of the discretized Poisson problem (1.1). If the same triangular meshes and nodal linear (Courant) ansatz functions are used, the matrix is already at our disposal from the computation of the discrete solution and the mass matrix is generated additively. Then the discretization of the Friedrichs' inequality reads for all vectors respecting zero Dirichlet boundary conditions. The minimal value to satisfy the last inequality represents the maximum Rayleigh quotient, that is, and it is also equal to the maximal eigenvalue of a generalized eigenvalue problem

Matrices and were assembled in MATLAB and the default function “eigs” was applied for the computation of the approximate values of the Friedrichs' constant . The MATLAB code can be downloaded at http://www.mathworks.com/matlabcentral/fileexchange/authors/37756 and it is easily extensible to any 2D geometry. Table 1 reports on the values of computed on the meshes and compares them with a theoretical value known for the unit square domain and zero Dirichlet boundary conditions. Note that the discrete approximations provide a lower bound only, but converge fast to the exact value of the Friedrich's constant. Here, we observe a quadratic convergence with respect to the mesh size . In general, the convergence typically depends on the shape of the domain boundary [27].

##### 3.2. Majorant Computation

One step (a) and one step (b) of Algorithm 2.2 assuming an initial value were applied for the computation of the flux . For the discretization of , Raviart-Thomas spaces of the zero degree (known as elements [28]) defined of the uniform meshes were considered.

*Remark 3.1 (Majorant problem size). *Note that the number of elements and number of edges of the triangulation satisfy recurrences
subject to the conditions , related to the coarse triangulation . These recurrences are solved to provide exact formulas , . By the known Euler formula for planar graphs, the number of vertices reads and it holds
for large . It implies that the matrix in Algorithm 2.2 using RT0 basis is asymptotically 3 times larger than the stiffness matrix from the discretization of Poisson problem by using linear nodal elements.

A MATLAB implementation is based on [29] with some modification with respect to the performance and extension towards a multigrid solver. Quadrature rules exact for polynomials up to the order two were taken for the computation of integrals on triangulations. The right column of Figure 2 displays computed fluxes (only one component due to symmetry reasons) and the exact flux (3.3). It can be observed that both discrete solutions and fluxes (at least visually) converge to exact solution and its flux. By comparing the values of Table 2 or Figure 3, both exact errors and majorant values converge linearly with respect to corresponding degrees of freedom used for their computation. The index of efficiency which remains bounded by the value for all mesh levels.

#### 4. Improving Linear Solver

The highest computation costs of Algorithm 2.2 are caused due to the solution of the linear system of (2.11) in step (a) of Algorithm 2.2. Let us consider an iterative method for its solving. The advantage of iterative over direct methods is apparent in this context, since each iteration vector for as an approximation the solution can be inserted into the majorant (2.10), without the need for solving the linear system very accurately. Since the linear system is represented by a symmetric and positive definite matrix, a preconditioned conjugate gradient method (PCG) is considered.

*Algorithm 4.1 (PCG method for a system of equations including energy computation). *Let an initial iteration vector be given. Compute an initial residual , an initial energy and , where is a given preconditioning matrix. For the iterations do the loop(1)(2)(3)(4)(5)(6)(7)(8)If a given stopping criterion is fulfilled, leave the loop and output the solution and the energy .

This algorithm recalls Algorithm 4.1 from [30] with a modification for the computation of an energy in the step (2). The energy is defined as and the formula in step (2) provides a cheap update of from without an extra matrix-vector multiplication in (4.1). It can be directly derived from the combination of two formulae. The first one is the known relation between the energy error of the PCG-iterations and their energy where denotes the exact solution of the linear system . The second one is a special version of the formula (3.6) for from [30] The knowledge of the energy is required for the computation of the flux-dependent functional majorant part , since it holds (cf. (2.10) where and are the matrix and the right-hand side of the linear system (2.11) and the vector substitutes . For the initial iteration , it holds and we obtain an estimate where is defined in (2.6). Since in (4.3), the PCG method reduces (or at least does not increase) the energy and consequently the value in each iteration to sharpen the estimate (4.5).

##### 4.1. Performance of PCG for Benchmark Example

Let us apply PCG method to one step (a) of Algorithm 2.2 to linear system (2.11), in which is considered. PCG method is terminated in step (8) of Algorithm 4.1 if the stopping criterion is fulfilled for some given tolerance , here we choose .

*Remark 4.2. *If the preconditioner well approximates , that is, , it holds , and (4.6) is equivalent to a stopping criterion based on the relative -norm of the error
In the case of no preconditioning (then we speak of CG method instead of PCG), is an identity matrix and (4.6) is equivalent to a stopping criterion
which is a default stopping criterion implemented in MATLAB.

As an operation of the preconditioner in Algorithm 4.1, we apply a simple V-cycle of a geometrical multigrid method [31] based on provided hierarchy of nested triangulations . As the linear system (2.11) arises from a discretization in space, a special smoother as a part of the multigrid method is required. Our choice is the additive version of the smoother of Arnold, Falk and Winther [32] using one presmoothing and one postsmoothing steps.

Table 3 compares numbers of iterations of non-preconditioned (CG) and multigrid-preconditioned (MPCG) method for various levels of triangulation. Single CG or MPCG iterations and the corresponding majorant values are displayed on Figure 4 for mesh levels 6 and 7. The number of iterations reflects typical properties of conjugate gradients and system matrices arising in elliptic partial boundary value problems. For shape regular triangulations, the condition number of is known to be proportional to , where is the mesh-size parameter, that is, Furthermore, the number of CG-iterations (with respect to the same stopping tolerance ) satisfies . Together, it holds The mesh size is halfened after each uniform refinement and therefore the number of CG-iterations is according to (4.9) expected to be doubled in the non-preconditioned case. For the multigrid-preconditioned CG method (MPCG) we observe that the number of iterations remains bounded, in our case is valid for all mesh levels. This property is so-called mesh independence and it demonstrates the optimality of the chosen multigrid preconditioner. A detailed observation of the right column pictures of Figure 4 indicates that 4 iterations () already provide a very sharp majorant value without the need for additional iterations.

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

The minimization of the majorant term in the functional a posteriori estimate is done by solving a sequence of systems of linear equations for an unknown approximation of the flux of an exact solution. The solution of the first linear system is efficiently obtained by the multigrid-preconditioned conjugate gradient method. The considered Benchmark example shows that few iterations already provide very accurate flux approximation for the computation of the functional majorant. Therefore, an optimal strategy for the termination process of the preconditioned conjugate gradient method in connection to the majorant computation remains an interesting open question.

#### Acknowledgments

The author acknowledges support from the Austrian Science Fund “Fonds zur Förderung der wissenschaftlichen Forschung (FWF)” for support under Grant SFB F013/F1306 in Linz, Austria, as well as the support of YFF Project “Modelling Transport in Porous Media over Multiple Scales” at the University Bergen, Norway. The comments of S. Repin (St. Petersburg) improved an overview on a posteriori estimates, “the Linz multigrid comunity” (U. Langer, R. Simon, and S. Zaglmayr) provided practical hints concerning multigrid application and A. Janka (Fribourg) assisted at the MATLAB vectorization of the smoother of Arnold, Falk, and Winther.