# Minimization of Functional Majorant in a Posteriori Error Analysis Based onH(div) Multigrid-Preconditioned CG Method

Authors: Jan Valdman
Journal: Advances in Numerical Analysis
 (2009)

---

## Abstract

We consider a Poisson boundary value problem and its functional a posteriori error estimate derived by S. Repin in 1999. The estimate majorizes theH1 seminorm of the error of the discrete solution computed by FEM method and contains a free ux variable from the H(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.


---

## Body

## 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(1.1)Δu+f=0inΩ,u=0on∂Ω
for a searched function u from a Sobolev space H∘1(Ω). The right hand side f∈L2(Ω) and an open bounded domain Ω∈ℝd, where d denotes a domain dimension, that is, d=1,2,3 are given. Assume that v∈H∘1(Ω) is an approximation of u. Then, a functional error estimate from [9]
(1.2)∥∇(u-v)∥Ω≤∥∇v-y∥Ω+CΩ∥divy+f∥Ω
holds for all functions y∈H(Ω,div), ∥·∥Ω denotes the L2(Ω) 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 CΩ included represents a constant from the Friedrichs' inequality
(1.3)∥w∥Ω≤CΩ∥∇w∥Ω
which holds for all w∈H∘1(Ω). Thus the constant CΩ depends on the domain Ω only and can be precomputed (it is demonstated in Subsection 3.1). For a given y∈H(Ω,div), the quality of the estimate (1.2) is measured by a ratio of its right hand and left hand side
(1.4)Ieff:=∥∇v-y∥Ω+CΩ∥divy+f∥Ω∥∇(u-v)∥Ω
known as an index of efficiency. It obviously holds Ieff≥1 and the equality Ieff=1 is valid for the choice y=∇u, that is, if y is chosen as the flux of an exact solution u. Having this interpretation of y in mind, there are known ways [18] how to compute a reasonable flux approximation y from the discrete solution v.(1)
Averaging on the mesh of the discrete solutionv. In this case, the flux approximation is computed as y=𝒢v, 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 solutionv 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 ofv, a function yf̅ is constructed such that divyf̅+f̅=0 and f̅ are sufficiently close to f in L2(Ω) norm. Then, the substitution of y=yf̅+τ into (1.2) provides an estimate
(1.5)∥∇(u-v)∥Ω≤∥∇v+curlτ-yf̅∥Ω+CΩ∥f-f̅∥Ω,
where τ∈H(Ω,curl) is arbitrary.(4)
Minimization of the right hand side of (1.2) with respect to free variable y on the mesh of the discrete solution v. 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 y 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 inH(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 solutionv∈H∘1(Ω), a right-hand side of the Poisson problem f∈L2(Ω), the Friedrich's constant belonging to the domain Ω. Find a function y∈H(Ω,div) satisfying the condition
(2.1)∥∇v-y∥Ω+CΩ∥divy+f∥Ω→min.In order to avoid complications with the nondifferentiability of norm terms in (2.1), we apply the inequality (a+b)2≤(1+β)a2+(1+(1/β))b2 valid for all β>0 to obtain
(2.2)(∥∇v-y∥Ω+CΩ∥divy+f∥Ω)2≤ℳ(∇v,f,β,CΩ,y),
where the upper bound in (2.2) denotes a functional majorant
(2.3)ℳ(∇v,f,β,CΩ,y):=(1+β)∥∇v-y∥Ω2+(1+1β)CΩ2∥divy+f∥Ω2.
The majorant arguments f,v,CΩ are known and β>0 and y∈H(Ω,div) are free parameters. For a fixed choice of parameters β, the majorant represents a quadratic functional in y. On the other hand, for a fixed y, the parameter
(2.4)β=CΩ∥divy+f∥Ω∥∇v-y∥Ω
minimizes ℳ amongst all positive β. It suggests the following solution algorithm to Problem 2.1.Algorithm 2.2 (Majorant minimization algorithm).
Givenβ>0.(a)
Computey from the minimization of the quadratic problem
(2.5)ℳ(∇v,f,β,CΩ,y)→min.(b)
Updateβ from y using (2.4). If the convergence in y is not achieved then go to step (a).For a detailed analysis of step (a) it is convenient to decomposeℳ in its y-independent and y-dependent parts as ℳ=ℳ1+ℳ2, where
(2.6)ℳ1:=(1+β)∥∇v∥Ω2+(1+1β)CΩ2∥f∥Ω2,ℳ2:=(1+β)(∥y∥Ω2-2(∇v,y))+(1+1β)CΩ2(∥divy∥Ω2+2(f,divy)).
Here, (·,·) denotes the L2(Ω) scalar product. Instead of considering (2.5) in Algorithm 2.2, y is solved from the minimization problem
(2.7)ℳ2(∇v,f,β,CΩ,y)→min
for given ∇v,f,β,CΩ. Since ℳ2=0 for y=0, the minimal value of ℳ2 must be nonpositive. Besides, it holds ℳ1≥0 and ℳ1=0 if and only if f=0 and v=0.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 T1 and T2 are either disjoint, or T1∩T2 is a complete edge or a common node of both T1 and T2. Let us assume a finite element basis in the space H(Ω,div) and a vector ŷ representing y∈H(Ω,div) in this basis. Norm terms in the definition of ℳ2 are read then after the discretization
(2.8)∥y∥Ω2=ŷTMŷ,∥divy∥Ω2=ŷTDIVDIVŷ,
where M and DIVDIV represent the ”mass” and “div div” matrices. Similarly, linear functionals are discretized as
(2.9)(∇v,y)=l1Tŷ,(f,divy)=l2Tŷ,
with some vectors l1 and l2. It allows to express a discrete form of ℳ2 in the provided finite element basis,
(2.10)ℳ2=ŷT[(1+β)M+(1+1β)CΩ2DIVDIV]ŷ-2ŷT[(1+β)l1-(1+1β)CΩ2l2].
The minimization of ℳ2 with respect to the vector ŷ leads for a given value β to a linear system of equations
(2.11)[(1+β)M+(1+1β)CΩ2DIVDIV]ŷ=(1+β)l1-(1+1β)CΩ2l2
for a minimizing vector ŷ. The linear system is represented by a matrix which is symmetric and positive definite. Symmetry is the result of symmetries of matrices M and DIVDIV. The positive definiteness is proved by a contradiction: if there was a nonzero vector ŷ for which the quadratic form in (2.10) was not positive, we could construct a sequence of arguments αiŷ, i=1,2,…, in which the values of ℳ2 would converge (for i→∞) to -∞. Consequently, since ℳ1 is independent of ŷ, 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(3.1)f(x,y)=2x(1-x)+2y(1-y)
in the unit square domain Ω=(0,1)×(0,1). The square geometry is discretized using a sequence of nested uniform triangular meshes 𝒯0,𝒯1,… as displayed in Figure 1. A discrete solution v is computed from by using FEM with the nodal linear (Courant) ansatz functions on each triangular mesh. The discrete solutions v are displayed as the left column pictures of Figure 2. For this particular right-hand side, the exact solution of (1.1) reads
(3.2)u=x(1-x)y(1-y)
for all (x,y)∈Ω and its flux is
(3.3)∂u∂x=(1-2x)y(1-y),∂u∂y=x(1-x)(1-2y)
for all (x,y)∈Ω. Therefore, the exact solution of the Benchmark example is a polynomial and, thus, integration error can be avoided.Uniform triangulations (a)𝒯0, (b) 𝒯1, and (c) 𝒯2.
(a)(b)(c)Discrete solutionsv and y-component of the flux y computed from Algorithm 2.2 on mesh triangulations 𝒯2,𝒯3,𝒯4, and the exact solution u and its flux y-component ∂u/∂y (in the bottom row denoted as 𝒯∞ triangulation).
(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 CΩ2 from the Friedrich's inequality (1.3). Under the assumption of the finite element basis ϕ1,…,ϕN we can represent the function w∈H∘1(Ω) as the vector ŵ∈ℝN, where w=∑i=1Nŵiϕi. In the same basis, the bilinear forms read
(3.4)(∇w,∇w)=ŵTKΔŵ,(w,w)=ŵTMΔŵ,
where KΔ and MΔ 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 KΔ is already at our disposal from the computation of the discrete solution v and the mass matrix is generated additively. Then the discretization of the Friedrichs' inequality reads
(3.5)ŵTMΔŵ≤CΩ2ŵTKΔŵ
for all vectors ŵ∈ℝN respecting zero Dirichlet boundary conditions. The minimal value CΩ2 to satisfy the last inequality represents the maximum Rayleigh quotient, that is,
(3.6)CΩ2=maxŵŵTMΔŵŵTKΔŵ,
and it is also equal to the maximal eigenvalue λ of a generalized eigenvalue problem
(3.7)MΔŵ=λKΔŵ.MatricesMΔ and KΔ were assembled in MATLAB and the default function “eigs” was applied for the computation of the approximate values of the Friedrichs' constant CΩ. 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 CΩ computed on the meshes 𝒯1,…,𝒯9 and compares them with a theoretical value
(3.8)CΩ=12π≈0.22507907903928
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 h. In general, the convergence typically depends on the shape of the domain boundary [27].Table 1
Friedrich's constant for the unit square domain computed from the generalized eigenvalue problem (3.7) on nested uniform meshes.
LevelMesh sizehErrorCΩ12-10.048302383742640.1767766952966322-20.015953561742050.2091255172972232-30.004245890661530.2208331883777442-40.001078758172360.2240003208669052-50.000270799181720.2248082798575562-60.000067769672600.2250113093666772-70.000016946794000.2250621322452782-80.000004236972100.2250748420671792-90.000001059260030.22507801977924Theory000.22507907903928


### 3.2. Majorant Computation

One step (a) and one step (b) of Algorithm2.2 assuming an initial value β=1 were applied for the computation of the flux y∈H(Ω,div). For the discretization of H(Ω,div), Raviart-Thomas spaces of the zero degree (known as RT0 elements [28]) defined of the uniform meshes 𝒯0,…,𝒯8 were considered.Remark 3.1 (Majorant problem size).
Note that the number of elements|𝒯n| and number of edges |ℰn| of the triangulation 𝒯n satisfy recurrences
(3.9)|𝒯n+1|=4|𝒯n|,|ℰn+1|=2|ℰn|+3|𝒯n|
subject to the conditions |𝒯0|=2, |ℰ0|=5 related to the coarse triangulation T0. These recurrences are solved to provide exact formulas |𝒯n|=2·4n, |ℰn|=2·2n+3·4n. By the known Euler formula for planar graphs, the number of vertices reads |𝒩n|=1+2·2n+4n and it holds
(3.10)|ℰn|≈3|𝒩n|≈32|𝒯n|
for large n. 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 Ieff≈1.73 for all mesh levels.Table 2
Comparison of the square of exact error∥∇(u−v)∥Ω2 and majorant values ℳ and on various uniform meshes levels.
LevelDOFv∥∇(u-v)∥Ω2DOFyMajorantℳIeff198.873457e-03162.135435e-021.552253.240138e-03569.050160e-031.673818.950033e-042082.635564e-031.7242892.295153e-048006.858397e-041.73510895.774695e-0531361.732119e-041.73642251.445987e-05124164.341359e-051.737166413.616415e-06494081.086033e-051.738660499.041944e-071971202.715516e-061.73Figure 3
Values of exact error∥∇(u-v)∥Ω  and the square root of the majorant ℳ displayed versus degrees of freedom (displayed as DOF v and DOF y in Table 2) of their calculation on various uniform mesh levels. Both lines indicate linear convergence with respect to the mesh size.


## 3.1. Friedrich's Constant and Its Computation

The computation of the functional majorant (2.3) requires the knowledge of the constant CΩ2 from the Friedrich's inequality (1.3). Under the assumption of the finite element basis ϕ1,…,ϕN we can represent the function w∈H∘1(Ω) as the vector ŵ∈ℝN, where w=∑i=1Nŵiϕi. In the same basis, the bilinear forms read
(3.4)(∇w,∇w)=ŵTKΔŵ,(w,w)=ŵTMΔŵ,
where KΔ and MΔ 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 KΔ is already at our disposal from the computation of the discrete solution v and the mass matrix is generated additively. Then the discretization of the Friedrichs' inequality reads
(3.5)ŵTMΔŵ≤CΩ2ŵTKΔŵ
for all vectors ŵ∈ℝN respecting zero Dirichlet boundary conditions. The minimal value CΩ2 to satisfy the last inequality represents the maximum Rayleigh quotient, that is,
(3.6)CΩ2=maxŵŵTMΔŵŵTKΔŵ,
and it is also equal to the maximal eigenvalue λ of a generalized eigenvalue problem
(3.7)MΔŵ=λKΔŵ.MatricesMΔ and KΔ were assembled in MATLAB and the default function “eigs” was applied for the computation of the approximate values of the Friedrichs' constant CΩ. 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 CΩ computed on the meshes 𝒯1,…,𝒯9 and compares them with a theoretical value
(3.8)CΩ=12π≈0.22507907903928
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 h. In general, the convergence typically depends on the shape of the domain boundary [27].Table 1
Friedrich's constant for the unit square domain computed from the generalized eigenvalue problem (3.7) on nested uniform meshes.
LevelMesh sizehErrorCΩ12-10.048302383742640.1767766952966322-20.015953561742050.2091255172972232-30.004245890661530.2208331883777442-40.001078758172360.2240003208669052-50.000270799181720.2248082798575562-60.000067769672600.2250113093666772-70.000016946794000.2250621322452782-80.000004236972100.2250748420671792-90.000001059260030.22507801977924Theory000.22507907903928


## 3.2. Majorant Computation

One step (a) and one step (b) of Algorithm2.2 assuming an initial value β=1 were applied for the computation of the flux y∈H(Ω,div). For the discretization of H(Ω,div), Raviart-Thomas spaces of the zero degree (known as RT0 elements [28]) defined of the uniform meshes 𝒯0,…,𝒯8 were considered.Remark 3.1 (Majorant problem size).
Note that the number of elements|𝒯n| and number of edges |ℰn| of the triangulation 𝒯n satisfy recurrences
(3.9)|𝒯n+1|=4|𝒯n|,|ℰn+1|=2|ℰn|+3|𝒯n|
subject to the conditions |𝒯0|=2, |ℰ0|=5 related to the coarse triangulation T0. These recurrences are solved to provide exact formulas |𝒯n|=2·4n, |ℰn|=2·2n+3·4n. By the known Euler formula for planar graphs, the number of vertices reads |𝒩n|=1+2·2n+4n and it holds
(3.10)|ℰn|≈3|𝒩n|≈32|𝒯n|
for large n. 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 Ieff≈1.73 for all mesh levels.Table 2
Comparison of the square of exact error∥∇(u−v)∥Ω2 and majorant values ℳ and on various uniform meshes levels.
LevelDOFv∥∇(u-v)∥Ω2DOFyMajorantℳIeff198.873457e-03162.135435e-021.552253.240138e-03569.050160e-031.673818.950033e-042082.635564e-031.7242892.295153e-048006.858397e-041.73510895.774695e-0531361.732119e-041.73642251.445987e-05124164.341359e-051.737166413.616415e-06494081.086033e-051.738660499.041944e-071971202.715516e-061.73Figure 3
Values of exact error∥∇(u-v)∥Ω  and the square root of the majorant ℳ displayed versus degrees of freedom (displayed as DOF v and DOF y in Table 2) of their calculation on various uniform mesh levels. Both lines indicate linear convergence with respect to the mesh size.


## 4. Improving Linear Solver

The highest computation costs of Algorithm2.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 ŷk for k=0,1,… as an approximation the solution ŷ 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 equationsAx=b including energy computation).
Let an initial iteration vectorx0 be given. Compute an initial residual r0=b-Ax0, an initial energy E0=(1/2)(Ax0,x0)-(b,x0) and s0=C-1r0,p0=s0,   where C is a given preconditioning matrix. For the iterations j=0,1,… do the loop(1)
γj=(rj,sj)/(pj,Apj)(2)
Ej+1=Ej-(γj(rj,sj)/2)(3)
xj+1=xj+γjpj(4)
rj+1=rj-γjApj(5)
sj+1=C-1rj+1(6)
δj+1=(rj+1,sj+1)/(rj,sj)(7)
pj+1=sj+1+δj+1pj(8)
If a given stopping criterion is fulfilled, leave the loop and output the solutionxj+1 and the energy Ej+1.This algorithm recalls Algorithm4.1 from [30] with a modification for the computation of an energy in the step (2). The energy is defined as
(4.1)Ej:=12(Axj,xj)-(b,xj),
and the formula in step (2) provides a cheap update of Ej+1 from Ej without an extra matrix-vector multiplication Axj+1 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
(4.2)Ej-E=12(Axj,xj)-12(Ax,x)-(b,xj)+(b,x)=12((Axj,xj)+(Ax,x)-2(Ax,xj))=12(A(x-xj),x-xj):=12∥xj-x∥A2,
where x denotes the exact solution of the linear system Ax=b. The second one is a special version of the formula (3.6) for d=1 from [30]
(4.3)∥x-xj+1∥A2=∥x-xj∥A2-γj(rj,sj).
The knowledge of the energy Ej is required for the computation of the flux-dependent functional majorant part ℳ2, since it holds (cf. (2.10)
(4.4)ℳ2(∇v,f,β,CΩ,ŷj)=2Ej,
where A and b are the matrix and the right-hand side of the linear system (2.11) and the vector x substitutes ŷ. For the initial iteration ŷ0=0, it holds ℳ2(∇v,f,β,CΩ,ŷ0)=0 and we obtain an estimate
(4.5)∥∇v-y∥Ω2≤ℳ1(∇v,f,β,CΩ),
where ℳ1 is defined in (2.6). Since γj(rj,sj)≥0 in (4.3), the PCG method reduces (or at least does not increase) the energy Ej and consequently the value ℳ(∇v,f,β,CΩ,ŷj) in each iteration j=1,2,… to sharpen the estimate (4.5).


### 4.1. Performance of PCG for Benchmark Example

Let us apply PCG method to one step (a) of Algorithm2.2 to linear system (2.11), in which β=1 is considered. PCG method is terminated in step (8) of Algorithm 4.1 if the stopping criterion
(4.6)(rj,sj)(r0,s0)≤ε
is fulfilled for some given tolerance ε, here we choose ε=10-6.Remark 4.2.
If the preconditionerC well approximates A, that is, C≈A, it holds (rj,sj)=(A(x-xj), C-1A(x-xj))≈∥x-xj∥A2 and (4.6) is equivalent to a stopping criterion based on the relative A-norm of the error
(4.7)∥x-xj∥A∥x-x0∥A≤ε.
In the case of no preconditioning (then we speak of CG method instead of PCG), C is an identity matrix and (4.6) is equivalent to a stopping criterion
(4.8)∥b-Axj∥2∥b∥2≤ε,
which is a default stopping criterion implemented in MATLAB.As an operation of the preconditionerC in Algorithm 4.1, we apply a simple V-cycle of a geometrical multigrid method [31] based on provided hierarchy of nested triangulations 𝒯0,𝒯1,𝒯2,…. As the linear system (2.11) arises from a discretization in H(Ω,div) 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.Table3 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 A is known to be proportional to h-2, where h is the mesh-size parameter, that is, cond(A)≈h-2.  Furthermore, the number of CG-iterations (with respect to the same stopping tolerance ε) satisfies ITER≈cond(A). Together, it holds
(4.9)ITER≈h-1.
The mesh size h 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
(4.10)ITER≤11
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 (ITER=4) already provide a very sharp majorant value without the need for additional iterations.Table 3
Number of CG/MPCG iterations (using stopping criterion (4.6) for ε=10−6) for various nested uniform mesh levels.
LevelMesh sizehProblem sizeCG iterationsMPCG iterations0152212-1165522-25615832-3208451042-48001081152-531362191162-6124163491172-7494088771182-819712017731192-9787456350611Majorant values displayed versus iterations of CG and MPCG method for𝒯6 and 𝒯7 meshes. For better illustration, the majorant part ℳ1 (the upper bound) from (2.6) and the square of the exact error, that is, ∥∇(u-v)∥Ω2 (the lower bound) are displayed.
(a)(b)(c)(d)


## 4.1. Performance of PCG for Benchmark Example

Let us apply PCG method to one step (a) of Algorithm2.2 to linear system (2.11), in which β=1 is considered. PCG method is terminated in step (8) of Algorithm 4.1 if the stopping criterion
(4.6)(rj,sj)(r0,s0)≤ε
is fulfilled for some given tolerance ε, here we choose ε=10-6.Remark 4.2.
If the preconditionerC well approximates A, that is, C≈A, it holds (rj,sj)=(A(x-xj), C-1A(x-xj))≈∥x-xj∥A2 and (4.6) is equivalent to a stopping criterion based on the relative A-norm of the error
(4.7)∥x-xj∥A∥x-x0∥A≤ε.
In the case of no preconditioning (then we speak of CG method instead of PCG), C is an identity matrix and (4.6) is equivalent to a stopping criterion
(4.8)∥b-Axj∥2∥b∥2≤ε,
which is a default stopping criterion implemented in MATLAB.As an operation of the preconditionerC in Algorithm 4.1, we apply a simple V-cycle of a geometrical multigrid method [31] based on provided hierarchy of nested triangulations 𝒯0,𝒯1,𝒯2,…. As the linear system (2.11) arises from a discretization in H(Ω,div) 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.Table3 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 A is known to be proportional to h-2, where h is the mesh-size parameter, that is, cond(A)≈h-2.  Furthermore, the number of CG-iterations (with respect to the same stopping tolerance ε) satisfies ITER≈cond(A). Together, it holds
(4.9)ITER≈h-1.
The mesh size h 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
(4.10)ITER≤11
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 (ITER=4) already provide a very sharp majorant value without the need for additional iterations.Table 3
Number of CG/MPCG iterations (using stopping criterion (4.6) for ε=10−6) for various nested uniform mesh levels.
LevelMesh sizehProblem sizeCG iterationsMPCG iterations0152212-1165522-25615832-3208451042-48001081152-531362191162-6124163491172-7494088771182-819712017731192-9787456350611Majorant values displayed versus iterations of CG and MPCG method for𝒯6 and 𝒯7 meshes. For better illustration, the majorant part ℳ1 (the upper bound) from (2.6) and the square of the exact error, that is, ∥∇(u-v)∥Ω2 (the lower bound) are displayed.
(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.



---

*Source: 164519-2009-07-12.xml*