On a residual-based a posteriori error estimator for the total error

On a residual-based a posteriori error estimator for the total error Abstract A posteriori error analysis in numerical partial differential equations aims at providing sufficiently accurate information about the distance of the numerically computed approximation to the true solution. Besides estimating the total error, a posteriori analysis should also provide information about its discretization and (inexact) algebraic computation parts. This issue has been addressed by many authors using different approaches. Historically, probably the first and practically very important approach is based on combination of the classical residual-based bound on the discretization error with the adaptive hierarchy of discretizations and computations that allow to incorporate, using various heuristic arguments, the algebraic error. Motivated by some recent publications, this text uses a complementary approach and examines subtleties of the (generalized) residual-based a posteriori error estimator for the total error that rigorously accounts for the algebraic part of the error. The aim is to show on the standard Poisson model problem example, which is used here as a case study, that a rigorous incorporation of the algebraic error represents an intriguing problem that is not yet completely resolved. That should be of concern in $$h$$-adaptivity approaches, where the refinement of the mesh is determined using the residual-based a posteriori error estimator assuming Galerkin orthogonality. The commonly used terminology such as ‘guaranteed computable upper bounds’ should be in the presence of algebraic error cautiously examined. 1. Introduction Historically, most a posteriori analysis in numerical partial differential equations (PDEs) focuses on estimating the discretization error, i.e., on the discrepancy between the solution of the original infinite-dimensional formulation of the problem and the exact solution of its discretized counterpart. This information is crucial for $$h$$-adaptivity, which refines discretization in the parts of the domain where the estimator indicates a large discretization error with the goal of achieving its close-to-uniform spatial distribution over the domain. Estimation of the discretization error, however, has to deal with the principal difficulty: the exact solution of the original problem is unknown, and unless the algebraic computations providing the coordinates of the discrete solution in the discretization basis are performed exactly or with a negligible algebraic error, the exact solution of the discretized problem is also unknown. Near-to-exact algebraic computations can be prohibitive due to extensive computational cost. When solving practical problems, one typically needs to estimate the error even for the computed approximation far from the exact solution of the discretized problem (see, e.g., Nordbotten & Bjørstad, 2008, Conclusions; Keilegavlen & Nordbotten, 2015; Nissen et al., 2015). Reaching exact algebraic results can even be theoretically prohibitive. The matrix eigenvalues, e.g., are (in general) in principle uncomputable by any finite formula as proved by the Abel–Galois theorem, and they can only be approximated iteratively. Moreover, in case of highly nonnormal matrices, there is no guaranteed forward estimate of the accuracy of the computed eigenvalue approximations and we can guarantee the backward error only. As a consequence, because of the inexactness of algebraic computations, the a posteriori error estimates should be based, from their derivation to their application, on the available computed approximations to the solution of the discrete problem. When considering simple model problems, the previous point is seemingly unimportant. The need for solving adaptively large-scale problems, however, requires abandoning the concept of highly accurate algebraic solutions. Incorporating multilevel discretization structure and the associated preconditioned iterative algebraic solvers with reliable stopping criteria is becoming a prerequisite for efficient very large-scale numerical PDE solvers. This has been clearly formulated as a program by Becker et al. (1995); see also other relevant references below. Since then, there is a growing number of work in this direction. The recent survey with many references can be found, e.g., in Arioli et al. (2013a, Section 4); see also the introductory parts of Jiránek et al. (2010), Arioli et al. (2013b), Papež et al. (2016) and (Málek & Strakoš, 2015, Chapter 12). Multilevel discretization structure is often obtained using discretization mesh adaptivity that requires, as mentioned above, local estimation of the discretization error in order to identify mesh elements that need to be refined. For that purpose a standard residual-based a posteriori bound on the discretization error is used (see, e.g., Babuška & Rheinboldt, 1978; Verfürth, 1996, Section 1.2; Ainsworth & Oden, 2000, Section 2.2, Brenner & Scott, 2008, Section 9.2). Since the exact solution of the algebraic problem is not available, the computed approximation that does not satisfy Galerkin orthogonality is often used in the bound, which violates the assumptions under which the bound has been derived. This raises a question to which extent the standard a posteriori error bounds, assuming Galerkin orthogonality can be modified to estimate the total error and, at the same time, to allow comparison of the discretization and algebraic part of the error and, consequently, construction of reliable stopping criteria for algebraic iterative solvers. The standard residual-based a posteriori bound is also a key ingredient of the works that use hierarchy of approximation spaces for estimating the total error and its algebraic part; for the early examples we refer, in particular, to Bramble et al. (1990), Xu (1992), Oswald (1993), Rüde (1993a,b), Becker et al. (1995) and Griebel & Oswald (1995). Therefore, we focus in this article on the residual-based a posteriori bound as a basic building block used elsewhere and investigate its extension for estimating the total error. We do not share the belief that the matter has been fully resolved, apart from simple technicalities. We describe difficulties that still remain open and have to be taken into account in various circumstances. Using, in particular, the articles by Becker & Mao (2009), Carstensen (1999) and Arioli et al. (2013b), we discuss the subtleties one has to deal with while estimating the total and discretization error using a residual-based a posteriori error estimator. We show that removing the standard Galerkin orthogonality assumption, which cannot be used in a large-scale practical application of the bound, requires a nontrivial revision of the known estimator. Even for the simple model problem, the derived extension of the bound contains multiplicative factors that are potentially very large, as shown also in Arioli et al. (2013b), and that cannot be, in general, easily and accurately determined. Moreover, despite the claims published in literature, there exists no a posteriori estimator for the algebraic part of the error that is cheap, easily computable and that gives in practice a tight guaranteed upper bound. As pointed out below, a rigorous a posteriori analysis that incorporates algebraic errors is for realistic problems substantially more difficult than the analysis that assumes exact algebraic computations. In Section 2 we set the notation, discuss some methodological questions and recall several approaches that use a residual-based a posteriori error estimator as a building block. Section 3 recalls the results from Carstensen (1999) on quasi-interpolation. Section 4 presents the revision of the upper bound on the total error from Becker & Mao (2009) and gives its detailed proof that abandons the Galerkin orthogonality assumption. Section 5 comments on the upper bound on the total error presented in Arioli et al. (2013b). Numerical illustrations of the difficulties associated with the multiplicative factors are present in Section 6. Section 7 addresses estimates of the algebraic error and explains misunderstandings present in literature. Finally, the conclusions are drawn. 2. Model problem and comments on methodology We will use the following standard model problem. Let $${\it {\Omega}} \subset \mathbb{R}^2$$ be a polygonal domain (open, bounded and connected set with a polygonal boundary). We consider the Poisson problem with the homogeneous Dirichlet boundary condition   \begin{equation} \text{find }u:{\it {\Omega}} \rightarrow \mathbb{R}\,{:}\qquad- {\it {\Delta}} u = f \quad \textrm{in} \quad {\it {\Omega}}, \qquad u = 0 \quad \textrm{on} \quad \partial {\it {\Omega}}, \label{eq:modelproblem} \end{equation} (2.1) where $$f:{\it {\Omega}} \rightarrow \mathbb{R}$$ is the source term. Hereafter, we use the standard notation for the Sobolev spaces. For $$D \subset {\it {\Omega}}$$, $$L^1(D)$$ denotes the space of the (Lebesgue) integrable functions in $$D$$, $$L^2(D)$$ denotes the space of the square integrable functions in $$D$$, $$(w,v)_D = \int_{D}{v\, w}$$ denotes the $$L^2$$-inner product on $$L^2(D)$$ and $$\left\| {{w}} \right\|_D = (w,w)^{1/2}_D$$ denotes the associated $$L^2$$-norm. We omit the subscripts for $$D = {\it {\Omega}}$$. $$H^k({\it {\Omega}})$$ denotes the Hilbert space of functions in $$L^2({\it {\Omega}})$$ whose weak derivatives up to the order $$k$$ belong to $$L^2({\it {\Omega}})$$. $$H^1_0({\it {\Omega}})$$ denotes the space of functions in $$H^1({\it {\Omega}})$$ with vanishing trace on the boundary $$\partial {\it {\Omega}}$$. Assuming $$f \in L^2({\it {\Omega}})$$, the problem (2.1) can be written in the following weak form   \begin{equation} \text{find } u\in V \equiv H^1_0({\it {\Omega}})\,{:}\qquad(\nabla u, \nabla v) = (\,f,v)\qquad \text{for}\,\text{all}~ v \in V. \label{eq:weakform} \end{equation} (2.2) Let $$\mathcal{T}$$ be a conforming triangulation of the domain $${\it {\Omega}}$$, i.e., two distinct and intersecting elements $$T_1, T_2 \in \mathcal{T}$$ share a common face, edge or vertex. Let $${\mathcal{N}}$$ denote the set of all nodes (i.e., the vertices of the elements of $$\mathcal{T}$$), whereas $${\mathcal{N}_\mathrm{int}} \equiv {\mathcal{N}} \backslash \partial {\it {\Omega}}$$ denotes the set of the free nodes. By $$\mathcal{E}$$ we denote the set of all edges of the elements of $$\mathcal{T}$$ and, similarly, $${\mathcal{E}_\mathrm{int}} \equiv \mathcal{E} \backslash \partial {\it {\Omega}}$$. For any node $$z \in {\mathcal{N}}$$, let $$\varphi_z$$ be the corresponding hat-function, i.e., the piecewise linear function that takes value 1 at the node $$z$$ and vanishes at all other nodes. By $$\omega_z$$ we denote the support of $$\varphi_z$$, which is equal to the patch $$\omega_z = \cup \left\{ T \in \mathcal{T} | z \in T \right\}$$. For an element $$T \in \mathcal{T}$$ we denote $$h_T \equiv {\text{diam}}(T)$$, similarly $$h_z \equiv {\text{diam}}(\omega_z)$$ denotes the diameter of $$\omega_z$$, $$z \in {\mathcal{N}}$$. By $$V_h \subset V$$ we denote the space of the continuous piecewise linear functions on the triangulation $$\mathcal{T}$$ vanishing on the boundary $$\partial {\it {\Omega}}$$, i.e., $$V_h \equiv \textrm{span}\left\{\varphi_z|z \in {\mathcal{N}_\mathrm{int}}\right\}$$. The discrete formulation of (2.2) then reads   \begin{equation} \text{find }u_h \in V_h\,{:}\qquad (\nabla u_h, \nabla v_h) = (\,f,v_h)\qquad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:discreteweakform} \end{equation} (2.3) The solution $$u_h$$ of (2.3) is called the Galerkin solution. Subtracting (2.3) from (2.2) and using $$V_h \subset V$$, we get the Galerkin orthogonality  \begin{equation} (\nabla (u - u_h), \nabla v_h) = 0 \qquad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:galerkinorthogonality} \end{equation} (2.4) The difficulty in estimating the discretization error $$u-u_h$$ mentioned above can be formulated as follows. Consider any estimator $$\text{EST}(\cdot)$$ that provides an upper bound   \begin{equation} \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(u_h), \label{eq:theorbound} \end{equation} (2.5) where $$\left\| {\left| {{\cdot}} \right|} \right\|$$ denotes an appropriate norm (for the model problem (2.1) typically the energy norm $$\left\| {\left| {{w}} \right|} \right\| = \left\| {{\nabla w}} \right\| = (\nabla w, \nabla w)^{1/2}$$). To evaluate the right-hand side of (2.5), we need $$u_h$$ that is not available. The common practice is then replacing $$u_h$$ by the computed approximation $$u_h^C$$, giving the seemingly easy solution   \[ \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(u_h^C). \] This inequality is, however, not guaranteed to hold without further justification that can be highly nontrivial or even impossible to achieve. Provided that   \begin{equation} \text{EST}(u_h) = \inf_{v_h \in V_h} \text{EST}(v_h), \label{eq:boundinfimum} \end{equation} (2.6) the bound (2.5) does indeed lead to a guaranteed upper bound   \begin{equation} \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(v_h)\quad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:theorbound2} \end{equation} (2.7) Proving (2.6) can, however, represent a challenge, and the authors are unaware of any applicable results of this kind published in the literature. Another option is simply writing   \begin{equation} \text{EST}(u_h) = \text{EST}(u_h^C) + \big(\text{EST}(u_h)-\text{EST}(u_h^C)\big)\label{eq:estdiff} \end{equation} (2.8) or using a variant of this based on a specific form of the estimator. Then the evaluation of the first term $$\text{EST}(u_h^C)$$ does not require any assumptions (it should not be confused with estimating the total error). Provided that the second term $$\text{EST}(u_h)-\text{EST}(u_h^C)$$ could be bounded using the computed quantity $$u_h^C$$, the relation (2.8) would give a rigorous bound on the discretization error. This consideration has been used in combination with heuristics in various approaches. Before focusing on the residual-based error estimator itself, we recall several ideas from the approaches combining this estimator with the hierarchy of discretizations. In the prototype paper (Becker et al., 1995), the error $$u - v_h$$ of any approximation $$v_h \in H^1_0({\it {\Omega}})$$ is expressed1 as   \begin{align*} \left\| {(u-v_h)} \right\|^2 &= (\,f,u-v_h) - (\nabla v_h, \nabla (u-v_h))\\ &=: \langle r(v_h), u-v_h \rangle, \end{align*} where $$r(v_h) \in H^{-1}({\it {\Omega}})$$ is the associated residual and $$\langle \cdot, \cdot \rangle : H^{-1}({\it {\Omega}}) \times H^1_0({\it {\Omega}}) \to \mathbb{R}$$ represents the duality pairing. Using the hierarchy of meshes and the associated discrete approximation subspaces $$V_0 \subset V_1 \subset \cdots \subset V_L \subset H^1_0({\it {\Omega}})$$, the article considers a multigrid algorithm with the Galerkin projection property and the operators $$I_j$$, $$j = 0,1, \ldots, L$$, where   \[ I_j: H^1_0({\it {\Omega}}) \to V_j \quad j = 0,1,\ldots,L. \] A straightforward substitution then gives   \begin{align} \langle r(v_h), u-v_h \rangle &= \langle r(v_h), (u-v_h) - I_L (u-v_h) \rangle \label{eq:BJR_bound1}\\ \end{align} (2.9)  \begin{align} &\quad{} + \sum^L_{j=1} \langle r(v_h), (I_j - I_{j-1})(u-v_h) \rangle \label{eq:BJR_bound2}\\ \end{align} (2.10)  \begin{align} &\quad{} + \langle r(v_h), I_0(u-v_h) \rangle. \label{eq:BJR_bound3} \end{align} (2.11) The first term (2.9) is then bounded using the standard residual-based a posteriori error estimator with the (non-Galerkin) input quantity $$v_h$$. The second term (2.10) is bounded using the algebraic residuals on the individual levels $$j = L, L-1, \ldots, 1$$. This will bring in nontrivial multiplicative factors analogous to those presented later in our article. Finally, the last term (2.11) is assumed to vanish because of the exact solution of the problem on the coarsest mesh. This assumption is substantial, as demonstrated also by the numerical experiment in Section 7 of the quoted article. The authors also suggest heuristics for stopping criteria in an adaptive algorithm and for approximation of the unknown multiplicative factors. The derivation does not consider roundoff error. There is a large amount of work that in principle can be put into the framework of (2.8), where $$\text{EST}(\cdot)$$ is again the residual-based error estimator and the difference $$(\text{EST}(u_h)-\text{EST}(u_h^C))$$ that reflects the algebraic error is estimated using the hierarchy of splittings of the approximation space $$H^1_0({\it {\Omega}})$$ or of its appropriate discretization; see, e.g., Bramble et al. (1990), Xu (1992), Oswald (1993), Rüde (1993a,b), Griebel & Oswald (1995), Becker et al. (1995) and Harbrecht & Schneider (2016), as well as further references in Arioli et al. (2013a, Section 4) (Málek & Strakoš, 2015, Chapter 12). The instructive article by Stevenson (2007) extends the approach of another remarkable article by Morin et al. (2002) on convergence of adaptive FEM, where the adaptivity is based on the residual-based a posteriori error bound on the discretization error. Stevenson’s rigorously presented results account for inexact algebraic computations, but they show that such extension is indeed highly intriguing. The main result in Stevenson (2007) relies on the continuity argument, i.e., it assumes that the algebraic solver deviates from the exact result in a sufficiently small way. For the algebraic solver, this article refers to the work of Wu & Chen (2006) on uniform convergence of multigrid V-cycle algorithm. The article by Wu & Chen (2006) assumes exact arithmetic and, in particular, the exact solution of the problem on the coarsest mesh. The recalled results underline the importance of understanding the extension of the residual-based a posteriori error bounds to inexact algebraic computations (non-Galerkin solutions). In In Becker & Ma (2009, Lemma 3.1) the bound on the total error is given in the form   \begin{equation} \left\| {{\nabla(u - v_h)}} \right\|^2 \leq \widetilde{C} \cdot\text{EST}(v_h) + 2 \,\left\| {{\nabla(u_h - v_h)}} \right\|^2, \label{eq:BecMaobound} \end{equation} (2.12) where $$\widetilde{C}$$ is stated to depend only on the minimal angle of the triangulation $$\mathcal{T}$$, and $$v_h\in V_h$$ is arbitrary, i.e., it can account for inexact algebraic computations, where $$v_h = u_h^C$$. Here, the first term $$\widetilde{C} \cdot\text{EST}(v_h)$$ represents the first term in (2.8) given by the standard residual-based a posteriori error estimator for the discretization error (with replacing the Galerkin solution $$u_h$$ by $$v_h$$). The second term, $$2 \left\| {{\nabla(u_h - v_h)}} \right\|^2$$ does not have the meaning of the second term in (2.8). The proof refers for the case $$v_h=u_h$$, i.e., for estimating the discretization error, to the article by Carstensen (1999). The proof is completed by arguing, without detailed explanation, that the general case (2.12) follows via the triangle inequality. To be valid and applicable in practical computations as an upper bound, the total error estimator of the form (2.12) must resolve two challenges. (1) It must rigorously justify using the arbitrary (non-Galerkin) $$v_h \in V_h$$ in the first term on the right-hand side of (2.12) and give the necessary information on the value of the factor $$\widetilde{C}$$. (2) It must be proved that in practical computations, we can indeed provide a meaningful (i.e., inexpensive and tight) upper bound on the norm $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ of the algebraic error. In this article, we examine the derivation of (2.12) and prove that an estimator of the analogous form can indeed be used also for approximations $$u_h^C$$ that do not solve the discretized problem exactly. Although the derivation is simple, and in comparison with the standard residual-based a posteriori error estimator assuming Galerkin orthogonality, the resulting bound seemingly only adds the term bounding the algebraic error, the whole matter is, in our opinion, not simple. First, as explained below, there is no proof that the given estimator provides in practice a guaranteed computable upper bound due to subtleties determining the multiplicative factors and due to effects of roundoff on estimating the algebraic error. Second, here we consider a simple model problem and point out difficulties that are not technical but methodological. It is not clear whether for more complicated problems an extension of the estimator preserves the same form or whether an estimator based on the related methodology can meaningfully incorporate algebraic errors. 3. Quasi-interpolation results This section presents results from Carstensen (1999) used further in the text. We include them here for completeness and self-consistency of the text. Denote by $$\psi$$ a piecewise linear function taking value 1 at the inner nodes $$z \in {\mathcal{N}_\mathrm{int}}$$ and vanishing on the boundary $$\partial {\it {\Omega}}$$, $$\psi \equiv \sum_{z\in{\mathcal{N}_\mathrm{int}}} {\varphi_z}$$. Then $$\varphi_z/\psi$$, $$z\in{\mathcal{N}_\mathrm{int}},$$ represents in $${\it {\Omega}}$$ a partition of unity. Indeed, since $$\varphi_z$$, $$z\in{\mathcal{N}_\mathrm{int}}$$, sum up to $$\psi$$, we have $$\sum_{z\in{\mathcal{N}_\mathrm{int}}}{\varphi_z/\psi} = 1$$ in $${\it {\Omega}}$$ (see Carstensen, 1999, Proposition 2.1). The quasi-interpolation operator $${\mathcal{I}}: L^1({\it {\Omega}}) \rightarrow V_h$$ is then defined in the following way. For a given $$w \in L^1({\it {\Omega}})$$,   \[ {\mathcal{I}} w \equiv \sum_{z \in {\mathcal{N}_\mathrm{int}}} {w_z \varphi_z}, \quad \textrm{where} \quad w_z \equiv \frac{(w, \varphi_z/\psi)}{(1, \varphi_z)}. \] The error $$w - {\mathcal{I}} w$$ has a vanishing weighted average. Namely, for $$w$$, $$R \in L^2({\it {\Omega}})$$ and arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{equation} \int_{{\it {\Omega}}} {R \, (w - {\mathcal{I}} w)} = \sum_{z\in {\mathcal{N}_\mathrm{int}}}{\int_{{\it {\Omega}}}{(R - R_z) (w - w_z \psi)(\varphi_z/\psi)}}, \label{eq:weightedaverage} \end{equation} (3.1) see Carstensen (1999, Remark 2.4). Since   \[ {\int_{{\it {\Omega}}}{(w - w_z \psi)(\varphi_z/\psi)}} = 0 \qquad \text{for}\,\text{all}~ z \in {\mathcal{N}_\mathrm{int}}, \] the numbers $$R_z \in \mathbb{R}$$ can be chosen arbitrarily. In particular, $$R_z$$ can be chosen as the mean value of $$R$$ on $$\omega_z$$. Then $$\int_{\omega_z}{|R-R_z|^2}$$ is minimal among all $$R_z \in \mathbb{R}$$. The following lemmas are stated and proved in Carstensen (1999) for a more general case. Considering the model problem (2.1), we restrict ourselves to the case $$w \in H^1_0({\it {\Omega}})$$. The multiplicative factors in the lemmas then depend on (I) The shape of $$\omega_z$$. (II) The shapes of $$\omega_{z\partial{\it {\Omega}}} \equiv (\omega_z \cup \omega_{\xi}\ |\ z \in {\mathcal{N}_\mathrm{int}}, \xi \in {\mathcal{N}} \backslash {\mathcal{N}_\mathrm{int}}, \omega_z \cap \omega_{\xi} \neq 0)$$. (III) The shape coefficients $$(\int_{\omega_z}{\varphi_z}/|\omega_z| \ |\ z\in {\mathcal{N}_\mathrm{int}})$$, where $$|\omega_z|$$ stands for the Lebesgue measure of $$\omega_z$$. (IV) The overlap   \[ M_1 \equiv \max_{z \in {\mathcal{N}_\mathrm{int}}} \mathrm{card} \{ \xi \in {\mathcal{N}} \backslash {\mathcal{N}_\mathrm{int}} \ |\ \omega_z \cap \omega_{\xi} \neq 0\}. \] (V) The shape of the elements $$T \in \mathcal{T}$$. (VI) The value $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\varphi_z}} \right\|_{\infty}}$$, where $$\left\| {{\cdot}_{\infty}} \right\|$$ denotes the $$L^{\infty}({\it {\Omega}})$$-norm and $$h_z = \text{diam}(\omega_z)$$. (VII) The value   \[ M_2 \equiv \underset {x \in {\it {\Omega}}}{\text{ess sup}}{ \{ h(x)/h_T\ |\ x \in T \in \mathcal{T} \} }, \] where $$h(x) \equiv \max\{h_z\ |\ \varphi_z(x) > 0, z \in {\mathcal{N}_\mathrm{int}}\}$$, $$h_T = \text{diam}(T)$$. The proofs of the lemmas use the Poincaré inequality on $$\omega_z$$ and the Friedrichs inequality on $$\omega_{z\partial{\it {\Omega}}}$$ defined in II. To prove Lemma 3.2, the so-called trace theorem (see, e.g., Carstensen, 1999, Proposition 4.1) is used on each element of the triangulation $$T \in \mathcal{T}$$; the multiplicative factor then depends on the shape of the elements; see V. The quantities $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\nabla \varphi_z}} \right\|_{\infty}}$$ and $$M_2$$ (see VI and VII) are of the order 1 on a shape-regular mesh, where $$\left\| {{\nabla \varphi_z|_T}} \right\|_{\infty} \approx h^{-1}_T$$ and $$h_z \approx h_T, T\in \omega_z$$. They deteriorate on a mesh consisting of triangles with small inner angles, where small and large elements (in the sense of their diameter) adjoin. To see the development of the argument, for clarity we present the following two lemmas. Lemma 3.1 (see Carstensen, 1999, Theorem 3.1, statement 1) There exists a multiplicative factor $$C>0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VI), but not on the size of the elements $$h_T$$, such that, $$\text{for}\,\text{all}~ R \in L^2({\it {\Omega}})$$, $$\text{for}\,\text{all}~ w \in H^{1}_0 ({\it {\Omega}})$$ and for arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{eqnarray} \int_{{\it {\Omega}}} {R\, (w - {\mathcal{I}} w)} \leq C \left\| {{\nabla w}} \right\| \left\{ \sum_{z \in {\mathcal{N}_\mathrm{int}}}{h^2_z \int_{\omega_z}{\varphi_z/\psi \ |R-R_z|^2}} \right\}^{1/2}. \nonumber \end{eqnarray} Lemma 3.1 is a consequence of the definition of the quasi-interpolation operator $${\mathcal{I}}$$; see (3.1). Lemma 3.2 (see Carstensen, 1999, Theorem 3.2) Let $$S \subset \mathcal{E}$$. There exists a multiplicative factor $$C > 0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VII), but not on the size of the elements $$h_T$$, such that $$\text{for}\,\text{all}~ J \in L^2(S)$$ and $$\text{for}\,\text{all}~ w \in H^1_0({\it {\Omega}})$$,   \[ \int_S{J \, (w - {\mathcal{I}} w)} \leq C \left\| {{\nabla w}} \right\| \left( \sum_{T\in\mathcal{T}}{h_T \left\| {{J}} \right\|^2_{S \cap \partial T}} \right)^{1/2}. \] Combining Lemmas 3.1 and 3.2 we get the final inequality. Lemma 3.3 (see Carstensen, 1999, Corollary 3.1) Let $$S \subset \mathcal{E}$$. There exists a multiplicative factor $$C_1>0$$ depending on I–VII such that, for all $$J \in L^2(S)$$, $$\text{for}\,\text{all}~ R \in L^2({\it {\Omega}})$$, for all $$w \in H^{1}_0 ({\it {\Omega}})$$ and for arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{eqnarray*} &&\int_{{\it {\Omega}}} {R\, (w - {\mathcal{I}} w)} + \int_{S} {J\, (w - {\mathcal{I}} w)} \\ &&\qquad \leq C_1 \left\| {{\nabla w}} \right\| \left\{ \sum_{z \in {\mathcal{N}_\mathrm{int}}}{h^2_z \left\| {{R - R_z}} \right\|^2_{\omega_z}} + \sum_{T \in \mathcal{T}} {h_T \left\| {{J}} \right\|^2_{S \cap \partial T}} \right\}^{1/2}. \end{eqnarray*} The following lemma introduces a positive multiplicative factor $$C_{\mathrm{intp}}$$ that plays a key role in our discussion on incorporating the algebraic error into the a posteriori bound on the total error; see Section 4 and the numerical experiments in Section 6. Lemma 3.4 (see Carstensen, 1999, Theorem 3.1, statement 3) There exists a multiplicative factor $$C_{\mathrm{intp}} > 0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VI and VI) such that, for all $$\text{for}\,\text{all}~ w \in H^{1}_0({\it {\Omega}})$$,   \begin{equation} \left\| {{{\mathcal{I}} w}} \right\| \leq C_{\mathrm{intp}} \left\| {{\nabla w}} \right\| . \label{eq:Cintp} \end{equation} (3.2) Remark 3.5 The factor $$C_{\mathrm{intp}}$$ satisfies   \[ \sup_{w \in H^{1}_0({\it {\Omega}})} \frac{\left\| {{{\mathcal{I}} w}} \right\|}{\left\| {{\nabla w}} \right\|} \leq C_{\mathrm{intp}}. \] Obtaining a value of $$C_{\mathrm{intp}}$$ such that the above inequality is tight represents a nontrivial issue. Using the proof of Carstensen (2006, Theorem 2.4) and the discussion in Carstensen (2006, Example 2.3),, we can get a better idea about the size of $$C_{\mathrm{intp}}$$. For a shape-regular mesh with $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\nabla \varphi_z}} \right\|_{\infty}} \approx 2$$ (see VI), there holds $$C_{\mathrm{intp}} \approx 6$$. In general, as stated in Carstensen (2006), it may be very large for small angles in the triangulation. For $$f \in L^1({\it {\Omega}})$$, define the mean-value operator $$\pi_{\omega_z}(\,f) \equiv \int_{\omega_z}{f}/|\omega_z|$$. We denote, for $$z \in {\mathcal{N}}$$ and for any subset $$Z \subset {\mathcal{N}}$$,   \[ {\text{osc}}_z \equiv |\omega_z|^{1/2} \left\| {{ f - \pi_{\omega_z} f}} \right\|_{\omega_z}, \qquad {\text{osc}}(Z) \equiv \left( \sum_{z \in Z}{{\text{osc}}^2_z} \right)^{1/2}, \] measuring the oscillations of $$f$$, i.e., the variations of the function $$f$$ from the mean value $$\pi_{\omega_z} f$$ on the subdomains $$\omega_z$$. Given $$v_h \in V_h$$, define for $$E \in {\mathcal{E}_\mathrm{int}}$$ and any subset $$F \subset {\mathcal{E}_\mathrm{int}}$$ the edge residuals   \[ J_E(v_h) \equiv |E|^{1/2} \left\| \left[ \frac{\partial v_h}{\partial n_E} \right] \right\|_{E}, \qquad J(v_h, F) \equiv \left( \sum_{E \in F}{J^2_E(v_h)} \right)^{1/2}, \] where $$[\cdot]$$ denotes the jump of a piecewise continuous function and $$n_E$$ denotes the unit normal to $$E$$ (for each $$E \in {\mathcal{E}_\mathrm{int}}$$ the orientation of the unit normal is set arbitrarily but fixed). The edge residual $$J_E(v_h)$$, $$v_h \in V_h,$$ measures the jump of the piecewise constant gradient $$\nabla v_h$$ over the inner edge $$E$$. We set for brevity $${\text{osc}} \equiv {\text{osc}}({\mathcal{N}})$$ and $$J(v_h) \equiv J(v_h, {\mathcal{E}_\mathrm{int}})$$. For a given $$v_h \in V_h$$, we define the jump function $$\mathcal{J}(v_h) \in L^2({\mathcal{E}_\mathrm{int}})$$ on the inner edges   \begin{equation} \mathcal{J}(v_h)|_E \equiv \left[ \frac{\partial v_h}{\partial n_E} \right]\quad E \in {\mathcal{E}_\mathrm{int}}. \end{equation} (3.3) The Green’s formula (see, e.g., Ciarlet, 2002, p. 14) gives for a domain $$D$$ with a Lipschitz continuous boundary $$\partial D$$ and for $$v \in H^2(D)$$, $$w \in H^1(D)$$  \begin{equation} \int_{D} \nabla v \cdot \nabla w = - \int_{D} {\it {\Delta}} v \, w + \int_{\partial D} \left(\frac{\partial v}{\partial n_{\partial D}} \right)\, w, \label{eq:Green} \end{equation} (3.4) where $$n_{\partial D}$$ denotes the unit normal to $$\partial D$$ pointing outward. Let $$v_h \in V_h$$ and $$T \in \mathcal{T}$$. Then $$v_h |_T$$ is a linear function, $$v_h |_T \in H^2(T)$$ and $${\it {\Delta}} v_h|_T = 0$$. Then, applying the Green’s formula (3.4) elementwise yields, for any $$v_h \in V_h$$, $$w \in H^1_0({\it {\Omega}})$$,   \begin{eqnarray} \int_{{\it {\Omega}}}{\nabla v_h \cdot \nabla w} &=& \sum_{T \in \mathcal{T}}{\int_T{\nabla v_h \cdot \nabla w}} = \sum_{T \in \mathcal{T}} \left( - {\int_T{{\it {\Delta}} v_h \, w}} + {\int_{\partial T}{\left(\frac{\partial v_h}{\partial n_{\partial T}} \right)\, w} } \right) \nonumber \\ &=& \sum_{E \in {\mathcal{E}_\mathrm{int}}}{\int_E {\left[ \frac{\partial v_h}{\partial n_E} \right] w}} = \int_{{\mathcal{E}_\mathrm{int}}}{\mathcal{J}(v_h)\, w} . \label{eq:perpartes} \end{eqnarray} (3.5) The results recalled in this section are used to prove the upper bound on the total error in the next section. 4. Upper bound on the total error Now we state the upper bound on the energy norm of the total error using the residual-based a posteriori error estimator. Theorem 4.1 There exist triangulation-dependent positive multiplicative factors $$C_1$$, $$C_{\mathrm{intp}}$$ and $$C_2$$ such that for the solution $$u$$ of (2.2), the Galerkin solution $$u_h$$ of (2.3) and an arbitrary $$v_h \in V_h$$,   \begin{equation} \label{eq:lemma31} \left\| {{\nabla (u - v_h)}} \right\|^2 \leq 2\, C^2_1\, C^2_2 \left( J^2(v_h) + {\text{osc}}^2 \right) + 2\, C^2_{\mathrm{intp}} \left\| {{ \nabla (u_h - v_h) }} \right\|^2. \end{equation} (4.1) In particular, $$C_1$$ depends on I–VII (see Lemma 3.3), $$C_{\mathrm{intp}}$$ depends on I–VI and VI (see Lemma 3.4) and the factor $$C_2$$ depends on the ratios $$h_z^2/|\omega_z|$$, $$z \in {\mathcal{N}_\mathrm{int}},$$ and $$h_T/|E|$$, $$T\in\mathcal{T}$$, $$E\in\partial T\cap {\mathcal{E}_\mathrm{int}}$$. Proof. We will use the standard expression for the norm   \begin{equation} \left\| {{\nabla (u - v_h) }} \right\| = \sup_{0 \neq w \in H^1_0({\it {\Omega}})} {\frac{1}{\left\| {{\nabla w}} \right\|} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w }}. \label{eq:dualnorm} \end{equation} (4.2) Let $$v_h \in V_h$$ and $$w \in H^1_0({\it {\Omega}}), w \neq 0,$$ be arbitrary,   \begin{eqnarray*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla {\mathcal{I}} w}\\ &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u - u_h) \cdot \nabla {\mathcal{I}} w} \\ && + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{eqnarray*} It follows from the definition of the interpolation operator that $${\mathcal{I}} w \in V_h$$. The Galerkin orthogonality (2.4) gives   \[ \int_{{\it {\Omega}}}{\nabla(u - u_h) \cdot \nabla {\mathcal{I}} w} = 0. \] Then   \begin{eqnarray*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w} \\ &=& \int_{{\it {\Omega}}}{\nabla u \cdot \nabla(w - {\mathcal{I}} w)} - \int_{{\it {\Omega}}}{\nabla v_h \cdot \nabla(w - {\mathcal{I}} w)}\\ && + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{eqnarray*} Using the weak formulation (2.2) and the equality (3.5),   \[ \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} = \int_{{\it {\Omega}}}{f ( w - {\mathcal{I}} w)} - \int_{{\mathcal{E}_\mathrm{int}}}{\mathcal{J}(v_h) ( w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \] Then Lemma 3.3 with $$S = {\mathcal{E}_\mathrm{int}}$$, $$R = f$$, $$R_z = \pi_{\omega_z}f,\ z \in {\mathcal{N}_\mathrm{int}}$$, $$J = -\mathcal{J}(v_h)$$ gives   \begin{align*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &\leq C_1 \left\| {{\nabla w}} \right\| \left\{ \sum_{T \in \mathcal{T}} {h_T \left\| {{\mathcal{J}(v_h)}} \right\|^2_{{\mathcal{E}_\mathrm{int}} \cap \partial T}} + \sum_{z \in {\mathcal{N}_\mathrm{int}}}{\! h^2_z \left\| {{f - \pi_{\omega_z}f}} \right\|^2_{\omega_z}} \right\}^{\frac12} \\ & \quad + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}\\ &\leq C_1 C_2 \left\| {{\nabla w}} \right\| \left\{ \sum_{E \in {\mathcal{E}_\mathrm{int}}} {|E| \, \left\| {{ \mathcal{J}(v_h) }} \right\|^2_{E}} + \sum_{z \in {\mathcal{N}_\mathrm{int}}}{|\omega_z|\, \left\| {{f - \pi_{\omega_z}f}} \right\|^2_{\omega_z}} \right\}^{\frac12} \\ & \quad + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}\\ &= C_1 C_2 \left\| {{\nabla w}} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{align*} Using the Cauchy–Schwarz inequality,   \begin{equation} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} \leq C_1 C_2 \left\| {{\nabla w}} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \left\| {{\nabla {\mathcal{I}} w}} \right\| \, \left\| {{ \nabla (u_h - v_h)}} \right\|. \label{eq:proofI} \end{equation} (4.3) Dividing (4.3) by $$\left\| {{\nabla w}} \right\|$$ and using Lemma 3.4,   \begin{eqnarray*} \frac{1}{\left\| {{\nabla w}} \right\|}\int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w } &\leq& C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \frac{\left\| {{\nabla {\mathcal{I}} w}} \right\|}{\left\| {{\nabla w}} \right\|} \left\| {{\nabla (u_h - v_h)}} \right\| \\ &\leq& C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + C_{\mathrm{intp}} \left\| {{\nabla(u_h - v_h)}} \right\|. \end{eqnarray*} Using the representation (4.2) of the energy norm, recall that $$w \in H^1_0({\it {\Omega}})$$ was chosen arbitrarily,   \begin{equation} \label{eq:totalboundCintp} \left\| {{\nabla (u - v_h)}} \right\| \leq C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + C_{\mathrm{intp}} \left\| {{ \nabla (u_h - v_h) }} \right\|. \end{equation} (4.4) Finally, we deduce (4.1) using the inequality $$(a + b)^2 \leq 2a^2 + 2b^2$$. The multiplicative factor $$C_{\mathrm{intp}}$$ is independent of the solution $$u \in H^1_0$$ of (2.2)2 and of its approximation $$v_h \in V_h$$. We will now give another bound on the total error that will contain information about the exact solution $$u$$ of (2.2). It is useful for illustration of the role of the factor $$C_{\mathrm{intp}}$$ in (4.1). Setting $$w \equiv u - v_h \in H^1_0({\it {\Omega}})$$ in (4.3), we get   \[ \left\| {{\nabla (u - v_h)}} \right\|^2 \leq C_1 C_2 \left\| {{ \nabla(u - v_h) }} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + {\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|} \, \left\| {{\nabla (u_h - v_h) }} \right\| . \] Dividing both sides by $$\left\| {{ \nabla(u - v_h) }} \right\|$$ gives   \begin{equation} \left\| {{ \nabla(u - v_h) }} \right\| \leq C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \frac {\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla (u - v_h) }} \right\|} \ \left\| {{ \nabla(u_h - v_h) }} \right\|, \label{eq:lemma31Cinttilde} \end{equation} (4.5) which is formulated as the following corollary. Corollary 4.2 Using the notation of Theorem 4.1,   \begin{equation} \label{eq:lemma31Ctilde} \left\| {{\nabla(u - v_h)}} \right\|^2 \leq 2\, C^2_1\, C^2_2 \left( J^2(v_h) + {\text{osc}}^2 \right) + 2\, \widetilde{C}^{\,2}_{\mathrm{intp}}(u,v_h) \,\left\| {{ \nabla(u_h - v_h) }} \right\|^2, \end{equation} (4.6) where   \begin{equation} \widetilde{C}_{\mathrm{intp}}(u, v_h) \equiv \frac{\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla (u - v_h) }} \right\|}\leq C_{\mathrm{intp}}. \label{eq:Cint} \end{equation} (4.7) Since $$C_{\mathrm{intp}}$$ is independent of $$u \in H^1_0$$ (independent of the source term $$f$$) and of $$v_h \in V_h$$, we must have   \begin{equation} \sup_{u\in H^1_0,\ u \text{ solves (2.2)},\ f \in L^2({\it {\Omega}})}\hspace{.5em} \sup_{v_h \in V_h} \frac{\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla(u - v_h) }} \right\|}\leq C_{\mathrm{intp}}. \label{eq:worstcase} \end{equation} (4.8) Therefore, $$C_{\mathrm{intp}}$$ represents a worst-case scenario factor and one may expect that most likely $$\widetilde{C}_{\mathrm{intp}}(u,v_h) \ll C_{\mathrm{intp}}$$. 5. A related result In Arioli et al. (2013b), the authors consider elliptic self-adjoint problems, and they use a residual-based error estimator for setting the stopping criterion for the conjugate gradient method (CG). Following their approach, one can easily get a theoretical upper bound on the total error. Although the bound is not stated explicitly in Arioli et al. (2013b), it appears in the proof of Theorem 3.3; see the inequality (Arioli et al., 2013b, (3.22)). The derivation proceeds differently from the proof of the bound (4.1); it again demonstrates that the price to be paid for removing the assumption on exact algebraic computations in terms of including unknown and possibly large multiplicative factors can be high. Moreover, in contrast to the statements in Arioli et al. (2013b), the resulting estimator cannot be considered a guaranteed practical upper bound due to the difficulties in estimating the algebraic part of the error; see Section 7 below. First, Arioli et al. (2013b, Theorem 2.2) recalls the bound on the discretization error: there exists a multiplicative factor $$C_{2.2}>0$$ that is dependent on the minimal angle of the triangulation $$\mathcal{T}$$, but independent of $$h$$, $$u$$ and $$u_h$$, such that   \begin{equation} \|{\nabla (u - u_h) }\|^2 \leq C_{2.2}\, \eta^2(u_h), \qquad \eta(u_h) \equiv \left( \sum_{T\in{\mathcal{T}}} |T| \, \|{ f + {\it{\Delta}} u_h}\|^2_T + (J(u_h))^2 \right)^{1/2},\label{eq:resbounddiscr} \end{equation} (5.1) see, e.g., Verfürth (1996, Section 1.2) and Ainsworth & Oden (2000, Section 2.2). Using inverse estimates for piecewise polynomial functions and Young’s inequality, Arioli et al. (2013b, Lemma 3.1) yields the inequality   \[ \eta^2(w_h) \leq (1+\gamma)\, \eta^2(v_h) + C_{3.1} (1+\gamma^{-1}) \left\| {{\nabla (v_h - w_h) }} \right\|^2 \quad \text{for}~\text{all}~ w_h, v_h \in V_h,\ \gamma > 0, \] where the positive factor $$C_{3.1}$$ depends on the minimal angle of the triangulation $$\mathcal{T}$$. Combining these bounds and the equality   \[ \left\| {{\nabla (u - v_h) }} \right\|^2 = \left\| {{\nabla (u - u_h) }} \right\|^2 + \left\| {{ \nabla(u_h - v_h) }} \right\|^2 \] that follows from the Galerkin orthogonality (2.4), we get the upper bound on the total error   \begin{eqnarray*} \left\| {{\nabla (u - v_h) }} \right\|^2 &\leq& C_{2.2}\, \eta^2(u_h) + \left\| {{\nabla (u_h - v_h) }} \right\|^2 \\ &\leq& C_{2.2} \, (1+\gamma)\, \eta^2(v_h) + \left(1 + C_{2.2}\, C_{3.1} \, (1+\gamma^{-1}) \right) \left\| {{\nabla (u_h - v_h) }} \right\|^2. \end{eqnarray*} Finally, by setting $$\gamma \equiv 1$$,   \begin{equation} \left\| {{ \nabla(u - v_h) }} \right\|^2 \leq 2\, C_{2.2} \, \eta^2(v_h) + \left(1 + 2\, C_{2.2}\, C_{3.1}\right) \left\| {{\nabla (u_h - v_h) }} \right\|^2. \label{eq:Ariolietalbound} \end{equation} (5.2) Comparing (5.1) with (5.2), we see that replacing the Galerkin solution $$u_h$$ by an arbitrary $$v_h \in V_h$$ results in the additional term $$\left\| {{ \nabla(u_h - v_h) }} \right\|^2$$, which is, however, multiplied by an unknown and potentially very large multiplicative factor $$(1 + 2\, C_{2.2}\, C_{3.1})$$. In the criteria proposed in Arioli et al. (2013b, Section 5) for numerical experiments, the factors are empirically set to $$C_{2.2} \equiv 40$$, $$C_{3.1} \equiv 10$$, giving $$(1 + 2\, C_{2.2}\, C_{3.1}) = 801$$; cf. (2.12). This underlines the subtleties of the residual-based bounds discussed above. In addition, as mentioned above and explained in detail in Section 7 below, getting a tight practical upper bound on $$\left\| {{\nabla (u_h - v_h) }} \right\|^2$$ represents an unresolved challenge. 6. Numerical illustrations We use, on purpose, very simple problems to illustrate the possible difference in the values of $$\widetilde{C}_{\mathrm{intp}}(u,v_h)$$ and $$C_{\mathrm{intp}}$$. Although $$\widetilde{C}_{\mathrm{intp}}(u,v_h)$$ can be, assuming the knowledge of the exact solution $$u$$, evaluated up to a negligible quadrature error, for the factor $$C_{\mathrm{intp}}$$, we present a lower bound given by plugging a chosen function into (3.2). The derivation of a more accurate estimate for $$C_{\mathrm{intp}}$$ (see also the discussion in Remark 3.5) is beyond the scope of this article. 6.1 Numerical illustration in one dimension We first consider a one-dimensional analog of the Clément-type quasi-interpolation operator $${\mathcal{I}}$$ to illustrate that $$C_{\mathrm{intp}}$$ can be significantly larger than one. Consider the domain $${\it {\Omega}} = (0,1)$$ with the (nonuniform) partition   \[ [0, \beta, {1}/{3} \pm \beta, {2}/{3} \pm \beta, 1-\beta, 1]\quad \beta = 0.01. \] This partition is adapted to the one-dimensional Laplace problem with the solution   \begin{eqnarray} \label{eq:atan-solution} u(x) &=& {\tan^{-1}}(cx) - {\tan^{-1}}(c(x-{1}/{3})) + {\tan^{-1}}(c(x-{2}/{3})) \\ && \qquad - {\tan^{-1}}(c(x-1)) - {\tan^{-1}}({c}/{3}) + {\tan^{-1}}({2c}/{3}) - {\tan^{-1}}(c), \nonumber \end{eqnarray} (6.1) with $$c = 1000$$. The left part of Fig. 1 depicts the solution $$u$$ and the Clément-type quasi-interpolant $${\mathcal{I}} u$$. For a zero approximate vector, we have Fig. 1. View largeDownload slide Left: the solution $$u$$ (6.1) (solid line) and the interpolant $${\mathcal{I}} u$$ (dotted line). Right: the function $$w(x) = x(1-x)$$ (solid line) and $${\mathcal{I}} w$$ (dotted line). Fig. 1. View largeDownload slide Left: the solution $$u$$ (6.1) (solid line) and the interpolant $${\mathcal{I}} u$$ (dotted line). Right: the function $$w(x) = x(1-x)$$ (solid line) and $${\mathcal{I}} w$$ (dotted line).   \begin{equation} \widetilde{C}_{\mathrm{intp}}(u,0) = \frac{\left\| {{({\mathcal{I}} u)'}} \right\|}{\left\| {{u'}} \right\|} = 0.77. \label{eq:Ctilde} \end{equation} (6.2) For a quadratic function $$w(x) = x(1-x)$$, $$w\in H^1_0({\it {\Omega}}),$$ we have   \[ \frac{\left\| {{({\mathcal{I}} w)'}} \right\|}{\left\| {{w'}} \right\|} = 3.70; \] see the right part of Fig. 1 for the plot of the function $$w$$ and the interpolant $${\mathcal{I}} w$$. Consequently, $$C_{\mathrm{intp}} \geq 3.70$$. 6.2 Two-dimensional numerical illustration For two-dimensional numerical illustration, we consider the square domain $${\it {\Omega}} \equiv (-1,1) \times (-1,1)$$ and the triangulation $$\mathcal{T}$$ generated by MATLAB3 command initmesh (‘squareg’, ‘Hmax’, 0.1)+ that provides a Delaunay triangulation consisting of 1 368 elements with the maximal diameter less than or equal to $$0.1$$. The minimal angle of the mesh is 35.9° and the average of the minimal angles of the elements is 50.3°. Consider the solution of problem (2.1):   \begin{equation} u^{(1)}(x,y) = (x-1)(x+1)(y-1)(y+1). \label{eq:u1} \end{equation} (6.3) For the zero approximate solution and the Galerkin solution $$u^{(1)}_h$$ corresponding to $$u^{(1)}$$, we have   \[ \widetilde{C}_{\mathrm{intp}}(u^{(1)},0) = 1.02, \qquad \widetilde{C}_{\mathrm{intp}}(u^{(1)},u^{(1)}_h) = 0.16. \] Similarly, for the exact solution   \begin{equation} u^{(2)}(x,y) = (x-1)(x+1)(y-1)(y+1)\exp(-100(x^2 + y^2)), \label{eq:u2} \end{equation} (6.4) we have   \[ \widetilde{C}_{\mathrm{intp}}(u^{(2)},0) = 0.76, \qquad \widetilde{C}_{\mathrm{intp}}(u^{(2)}, u^{(2)}_h) = 0.28. \] In Fig. 2 we show the values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)},v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$ defined above. Fig. 2. View largeDownload slide The values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)}, v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$; see (6.3) and (6.4), respectively. Fig. 2. View largeDownload slide The values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)}, v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$; see (6.3) and (6.4), respectively. To bound the constant $$C_{\mathrm{intp}}$$ from below, we consider $$w_h \in V_h$$ such that   \begin{equation} w_h(z) = 1 \quad z \in {\mathcal{N}_\mathrm{int}}, \qquad w_h = 0 \quad \text{on } \partial {\it {\Omega}}. \label{eq:wh} \end{equation} (6.5) For this function   \[ 1.10 = \frac{\left\| {{\nabla {\mathcal{I}} w_h }} \right\|}{\left\| {{\nabla w_h }} \right\|} \leq C_{\mathrm{intp}}. \] Figure 3 gives the difference $$w_h-{\mathcal{I}} w_h$$. This is a piecewise linear function that is on the machine precision level in most of the domain except patches around the inner nodes adjacent to the boundary $$\partial {\it {\Omega}}$$. We recall that for this simple problem and a shape-regular mesh $$C_{\mathrm{intp}} \approx 6$$, see Remark 3.5 and the original article by Carstensen (2006). It can therefore indeed hold $$C_{\mathrm{intp}} \gg \widetilde{C}_{\mathrm{intp}}(u, v_h)$$. Fig. 3. View largeDownload slide The difference $$w_h-{\mathcal{I}} w_h$$ for $$w_h$$ given by (6.5). Fig. 3. View largeDownload slide The difference $$w_h-{\mathcal{I}} w_h$$ for $$w_h$$ given by (6.5). 7. Algebraic error estimates Various published articles that include algebraic errors in a posteriori error analysis consider preconditioned iterative methods (such as CG) for solving the associate symmetric positive definite algebraic problem. Estimating the algebraic error $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ is then often considered to be resolved, apart from some seemingly technical issues such as the choice of a number of extra CG iterations, or the approximation of the smallest eigenvalue of the system matrix (see, e.g., Arioli et al., 2013b, Section 4), where the presented Gauss–Radau quadrature-based bound is considered “the only guaranteed upper bound for the $$A$$-norm of the CG error.” Here, we present mathematical arguments that challenge this opinion. In short, as above, theoretical results with assumptions that do not hold in practical computations cannot be applied without an appropriate theoretical justification, here the numerical stability analysis. In practice, the Gauss–Radau quadrature-based estimate does not give, in general, a guaranteed and tight upper bound. The Gauss–Radau quadrature requires an approximation of the smallest eigenvalue of the system matrix from below. A close approximation to the smallest eigenvalue from above can be in theory obtained in the process of (exact arithmetic) CG computations. A closer investigation reveals, however, a serious difficulty. If the approximation to the smallest eigenvalue is not accurate enough, the estimate for the norm of the algebraic error can be inaccurate. If, on the other hand, the eigenvalue approximation approaches the smallest eigenvalue, then the computation of the Gauss–Radau estimator must be cautiously done in a numerically stable way, since it (implicitly) inverts a matrix that becomes in such case close to numerically singular. Moreover, if a second approximation of the smallest eigenvalue is formed due to loss of linear independence of the computed vectors generating the associated Krylov subspaces, then the Gauss–Radau quadrature estimates exhibit instabilities that are not yet understood. We can see again the generally valid principle: In evaluation of a posteriori error estimates, all sources of errors, including roundoff, must be taken into consideration. This is also illustrated by the following point, which is valid even if the difficulty with a tight approximation of the smallest eigenvalue from below is resolved. In numerical computations, we cannot guarantee that Gauss–Radau quadrature estimates give an upper bound due to effects of roundoff errors to the underlying short recurrences in CG computations. This is a principal issue as explained next. Because of the loss of orthogonality caused by roundoff errors, the formulas that express the error norms in CG iterations under the assumption of exact arithmetic and, consequently, orthogonality among the computed basis vectors are no longer valid. In practical computations, we neither have orthogonal bases nor we have Krylov subspaces in the strict mathematical meaning. Unless the problem is really computationally simple (as in the Poisson model problem),4 orthogonality and also linear independence among the vectors computed using short recurrences is, in general, ratherquickly lost and convergence of the computed iterations is substantially delayed. In terms of the computed Jacobi matrices that are substantial in the derivation of the Gauss–Radau quadrature bounds, their entries will quickly become far away (orders of magnitude) from their theoretical counterparts. Recent description of the related issues can be found in the article by Gergelits & Strakoš(2014) with the references to many earlier articles (see, e.g., Paige, 1980; Greenbaum, 1989; Strakoš, 1991; Greenbaum & Strakoš, 1992; Notay, 1993; Strakoš & Tichý, 2002; Strakoš & Tichý, 2005; Meurant & Strakoš, 2006; O’Leary et al., 2007). The matter has been comprehensively discussed already in the survey article by Strakoš & Liesen (2005) and it is also covered in the recent monograph (Liesen & Strakoš, 2013, Section 5.9). In summary, the following principal question is omitted in works assuming exact arithmetic: Considering the effects of roundoff recalled above, how do we know that the formulas derived under the assumptions that are so drastically violated give anything meaningful in practical computations? This fundamental issue cannot be resolved by heuristic arguments; it requires thorough analysis. The question on the effects of roundoff errors to error estimation in iterative methods such as CG has already been raised in the seminal article by Dahlquist et al. (1979). It has been investigated in Section 5 (called ‘Rounding error analysis’) of the article by Golub & Strakoš (1994), and again very thoroughly (in the context of different estimates) in the article by Strakoš & Tichý (2002) with the title ‘On error estimation in the conjugate gradient method and why it works in finite precision computations’; see also the survey article by Meurant & Strakoš (2006). This underlines the point. Without thorough rounding error analysis, we may say that we have observed some behaviour on some examples and nothing more. With thorough rounding error analysis, we can explain why and under which conditions some estimates work and prove them numerically safe. Perhaps even more important, we can prove that other (mathematically equivalent) estimates can behave in a numerically unstable way and they should not be used. Estimates that have been proved numerically unstable (see, e.g., Strakoš & Tichž, 2002) are unfortunately indeed frequently used in practice. The facts presented, e.g., in Golub & Strakoš (1994, Section 5), Strakoš & Tichý(2002, Section 7), Meurant (2006, Chapter 7), Meurant & Strakoš (2006, Section 5), O’Leary et al. (2007) and Liesen & Strakoš (2013, Section 5.9), show that without a thorough and rigorous analysis, application of the error bounds derived under the assumption of exact computation, to the results of finite precision computations, is not only methodologically wrong, but it can also indeed lead to false conclusions. A partial theoretical justification of the Gauss and Gauss–Radau quadrature estimates is provided (based on the earlier works of Paige, Greenbaum and others) in Section 5 of the article by Golub & Strakoš (1994). In short, justification of the quadrature-based bounds must be based on the Riemann–Stieltjes integration using the distribution function, with possibly many more points of increase than the original distribution function determined by the data of the problem; see also the associated arguments in the article by O’Leary et al. (2007) on sensitivity of the Gauss quadrature. The article by Golub & Strakoš (1994) mentioned earlier justifies using Gauss and Gauss–Radau quadrature estimates in finite precision computations (with limitations specified in the article). The estimates based on the Gauss–Radau quadrature can be useful, but they cannot be proved to give a tight guaranteed upper bound on the error norms. Summarizing, accurate and computationally efficient estimation of the algebraic error $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ still represents a challenge. This challenge is not resolved by deriving estimators assuming exact computations and then plugging in the computed quantities. The fact that some estimators can be used in the same form for input entries computed using finite precision arithmetic is not at all obvious and it cannot be guessed a priori by any heuristics. This fact can only result from a rigorous analysis that is (in these cases) highly nontrivial. 8. Conclusion The presented article investigates changes in derivation and application of the standard residual-based a posteriori error bound on the discretization error needed to get an estimator (not necessarily a practically applicable guaranteed upper bound) for the total error. Technically, this article provides a detailed proof of the residual-based upper bound on the total approximation error that requires knowledge of the associated multiplicative factors. As published previously in Becker & Mao (2009) and Arioli et al. (2013b), abandoning the Galerkin orthogonality assumption in the derivation leads naturally to an additional term accounting for the algebraic part of the error. We show that the contribution of the algebraic error is scaled by a nontrivial multiplicative factor; see (4.1) and (4.5)–(4.7). This multiplicative factor $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$ depends, besides the computed approximation $$v_h$$, also on the unknown infinite-dimensional solution $$u$$ of (2.2). Therefore, it is generally uncomputable. It can be bounded, using the a priori information, by the factor $$C_{\mathrm{intp}}$$ independent of $$u$$ and $$v_h$$ given by (3.2). The value of $$C_{\mathrm{intp}}$$ can therefore overestimate the value of $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$. Moreover, as another substantial point, a guaranteed computationally efficient upper bound on the algebraic part of the error that provides, in general, sufficiently tight results is not available. The main message of this article does not concern technical details about multiplicative factors in the residual-type a posteriori error estimators for the total error developed for the given model problem. It shows that handling inexact algebraic computations still needs, despite many remarkable results, further work. In practical computations, we cannot avoid using various heuristics. This article supports using heuristics (they are used in many articles coauthored by the authors of this article). It points out, however, a need for supporting analysis. It shows that even for the simple model problem and the standard residual-based a posteriori error estimator, the matter is not easy and the unresolved questions can be practically important. As, e.g., numerically illustrated in in Papež (2016, Section 4.2) and as mentioned in the Introduction, application of the residual-based error estimator for the mesh refinement adaptivity remains in the presence of algebraic errors an open problem. When the standard estimator for the discretization error is evaluated using computed quantities, there is no guarantee that its local contributions provide a meaningful indication of the spatial distribution of the discretization error over the domain. The matter is practically very important because it can affect efficiency of $$h$$-adaptive computations. Extrapolation of the observations obtained for simple model problems cannot be considered the final and generally valid justification. The fact that in (4.1) and (4.6) the algebraic error is estimated globally, and the multiplicative factors $$C_{\mathrm{intp}}$$ and $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$ are not easy to estimate, suggests that the matter is intriguing and it requires further work. This article cannot survey the previously published and recently developed results towards robust stopping criteria that balance the algebraic and discretization error. This is addressed, e.g., in the works using the hierarchy of subspace splittings recalled in Section 2 (see also Huber et al., 2017) or in Arioli et al. (2013a, Section 4), Jirànek et al. (2010), Papež et al. (2014), Carstensen et al. (2014, Section 7.1), and Papež et al. (2016). In Papež et al. (2016), it is shown, using the methodology based on flux reconstruction, that such stopping criteria could indeed be constructed and rigorously supported by analysis. It also shows, however, that in practical applications, there is a substantial computational cost to be paid, and this cost can even become in some cases excessive. There is still a work to be done. As shown at the example of the residual-based a posteriori error estimator in this article, such work should go hand in hand with analysis. Practical importance of Theorem 4.1 is not in application of the bound (4.1). It is in the warning that such application can be difficult and unjustified heuristics can be misleading. In particular, any use of the adjective ‘guaranteed’ should carefully examine the assumptions under which this adjective holds in practice. Acknowledgements The authors are grateful to David Silvester, Endre Süli and an anonymous referee for valuable comments, which improved the text. Funding ERC-CZ project LL1202 financed by the Ministry of Education, Youth and Sports of the Czech Republic. References Ainsworth M. & Oden J. T. ( 2000) A Posteriori Error Estimation in Finite Element Analysis . Pure and Applied Mathematics (New York). New York: Wiley-Interscience [John Wiley & Sons] pp. xx+ 240. Google Scholar CrossRef Search ADS   Arioli M., Liesen J., Miȩdlar A. & Strakoš Z. ( 2013a) Interplay between discretization and algebraic computation in adaptive numerical solution of elliptic PDE problems. GAMM-Mitt. , 36, 102– 129. Google Scholar CrossRef Search ADS   Arioli M., Georgoulis E. H. & Loghin D. ( 2013b) Stopping criteria for adaptive finite element solvers. SIAM J. Sci. Comput. , 35, A1537– A1559. Google Scholar CrossRef Search ADS   Babuška I. & Rheinboldt W. C. ( 1978) Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. , 15, 736– 754. Google Scholar CrossRef Search ADS   Becker R. & Mao S. ( 2009) Convergence and quasi-optimal complexity of a simple adaptive finite element method. M2AN Math. Model. Numer. Anal. , 43, 1203– 1219. Google Scholar CrossRef Search ADS   Becker R., Johnson C. & Rannacher R. ( 1995) Adaptive error control for multigrid finite element methods. Computing , 55, 271– 288. Google Scholar CrossRef Search ADS   Bramble J. H., Pasciak J. E. & Xu J. ( 1990) Parallel multilevel preconditioners. Numerical Analysis 1989 (Dundee, 1989) . Pitman Res. Notes Math. Ser., vol. 228. Harlow: Longman Sci. Tech., pp. 23– 39. Google Scholar CrossRef Search ADS   Brenner S. C. & Scott L. R. ( 2008) The Mathematical Theory of Finite Element Methods . Texts in Applied Mathematics, vol. 15, 3rd edn. New York: Springer, pp. xviii+ 397. Google Scholar CrossRef Search ADS   Carstensen C. ( 1999) Quasi-interpolation and a posteriori error analysis in finite element methods. M2AN Math. Model. Numer. Anal. , 33, 1187– 1202. Google Scholar CrossRef Search ADS   Carstensen C. ( 2006) Clément interpolation and its role in adaptive finite element error control. Partial Differential Equations and Functional Analysis . Oper. Theory Adv. Appl., vol. 168. Basel: Birkhäuser, pp. 27– 43. Google Scholar CrossRef Search ADS   Carstensen C., Feischl M., Page M. & Praetorius D. ( 2014) Axioms of adaptivity. Comput. Math. Appl. , 67, 1195– 1253. Google Scholar CrossRef Search ADS PubMed  Ciarlet P. G. ( 2002) The Finite Element Method for Elliptic Problems . Classics in Applied Mathematics, vol. 40. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xxviii+ 530 (Reprint of the 1978 original [North-Holland, Amsterdam]). Google Scholar CrossRef Search ADS   Dahlquist G., Golub G. H. & Nash S. G. ( 1979) Bounds for the error in linear systems. Semi-infinite programming (Proc. Workshop, Bad Honnef, 1978) . Lecture Notes in Control and Information Science vol. 15. Berlin: Springer, pp. 154– 172. Google Scholar CrossRef Search ADS   Gergelits T. & Strakoš Z. ( 2014) Composite convergence bounds based on Chebyshev polynomials and finite precision conjugate gradient computations. Numer. Algorithms , 65, 759– 782. Google Scholar CrossRef Search ADS   Golub G. H. & Strakoš Z. ( 1994) Estimates in quadratic formulas. Numer. Algorithms , 8, 241– 268. Google Scholar CrossRef Search ADS   Greenbaum A. ( 1989) Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences. Linear Algebra Appl. , 113, 7– 63. Google Scholar CrossRef Search ADS   Greenbaum A. & Strakoš Z. ( 1992) Predicting the behavior of finite precision Lanczos and conjugate gradient computations. SIAM J. Matrix Anal. Appl. , 13, 121– 137. Google Scholar CrossRef Search ADS   Griebel M. & Oswald P. ( 1995) On the abstract theory of additive and multiplicative Schwarz algorithms. Numer. Math. , 70, 163– 180. Google Scholar CrossRef Search ADS   Harbrecht H. & Schneider R. ( 2016) A note on multilevel based error estimation. Comput. Methods Appl. Math. , 16, 447– 458. Google Scholar CrossRef Search ADS   Huber M., Rüde U., Strakoš Z. & Wohlmuth B. ( 2017) Error estimators for highly parallel multigrid solver (in preparation). Jiránek P., Strakoš Z. & Vohralík M. ( 2010) A posteriori error estimates including algebraic error and stopping criteria for iterative solvers. SIAM J. Sci. Comput. , 32, 1567– 1590. Google Scholar CrossRef Search ADS   Keilegavlen E. & Nordbotten J. M. ( 2015) Inexact linear solvers for control volume discretizations in porous media. Comput. Geosci. , 19, 159– 176. Google Scholar CrossRef Search ADS   Liesen J. & Strakoš Z. ( 2013) Krylov Subspace Methods: Principles and Analysis . Numerical Mathematics and Scientific Computation. Oxford: Oxford University Press. Málek J. & Strakoš Z. ( 2015) Preconditioning and the Conjugate Gradient Method in the Context of Solving PDEs.  SIAM Spotlights, vol. 1. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. x+ 104. Meurant G. ( 2006) The Lanczos and Conjugate Gradient Algorithms. From Theory to Finite Precision Computations.  Software, Environments, and Tools, vol. 19. Philadelphia, PA: SIAM, pp. xvi+ 365. Google Scholar CrossRef Search ADS   Meurant G. & Strakoš Z. ( 2006) The Lanczos and conjugate gradient algorithms in finite precision arithmetic. Acta Numer. , 15, 471– 542. Google Scholar CrossRef Search ADS   Morin P., Nochetto R. H. & Siebert K. G. ( 2002) Convergence of adaptive finite element methods. SIAM Rev. , 44, 631– 658 (2003) (Revised reprint of ‘Data oscillation and convergence of adaptive FEM’ [SIAM J. Numer. Anal., 38 (2000), no. 2, 466–488]). Google Scholar CrossRef Search ADS   Nissen A., Pettersson P., Keilegavlen E. & Nordbotten J. M. ( 2015) Incorporating geological uncertainty in error control for linear solvers. SPE Reservoir Simulation Symposium . Society of Petroleum Engineers. DOI: http://doi.org/10.2118/173272-MS Nordbotten J. M. & Bjørstad P. E. ( 2008) On the relationship between the multiscale finite-volume method and domain decomposition preconditioners. Comput. Geosci. , 12, 367– 376. Google Scholar CrossRef Search ADS   Notay Y. ( 1993) On the convergence rate of the conjugate gradients in presence of rounding errors. Numer. Math. , 65, 301– 317. Google Scholar CrossRef Search ADS   O’Leary D. P., Strakoš Z. & Tichý P. ( 2007) On sensitivity of Gauss-Christoffel quadrature. Numer. Math. , 107, 147– 174. Google Scholar CrossRef Search ADS   Oswald P. ( 1993) On a BPX-preconditioner for $$\rm P1$$ elements. Computing , 51, 125– 133. Google Scholar CrossRef Search ADS   Paige C. C. ( 1980) Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem. Linear Algebra Appl. , 34, 235– 258. Google Scholar CrossRef Search ADS   Papež J. ( 2016) Algebraic error in matrix computations in the context of numerical solution of partial differential equations. Ph.D. thesis . Prague: Charles University. Papež J., Liesen J. & Strakoš Z. ( 2014) Distribution of the discretization and algebraic error in numerical solution of partial differential equations. Linear Algebra Appl. , 449, 89– 114. Google Scholar CrossRef Search ADS   Papež J., Strakoí Z. & Vohralík M. ( 2016) Estimating and localizing the algebraic and total numerical errors using flux reconstructions. Preprint MORE/2016/12 (accepted for publication in Numer. Math.). Rüde U. ( 1993a) Fully adaptive multigrid methods. SIAM J. Numer. Anal. , 30, 230– 248. Google Scholar CrossRef Search ADS   Rüde U. ( 1993b) Mathematical and Computational Techniques for Multilevel Adaptive Methods . Frontiers in Applied Mathematics, Vol. 13. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xii+ 140. Google Scholar CrossRef Search ADS   Stevenson R. ( 2007) Optimality of a standard adaptive finite element method. Found. Comput. Math. , 7, 245– 269. Google Scholar CrossRef Search ADS   Strakoš Z. ( 1991) On the real convergence rate of the conjugate gradient method. Linear Algebra Appl. , 154–156, 535– 549. Google Scholar CrossRef Search ADS   Strakoš Z. & Liesen J. ( 2005) On numerical stability in large scale linear algebraic computations. ZAMM Z. Angew. Math. Mech. , 85, 307– 325. Google Scholar CrossRef Search ADS   Strakoš Z. & Tichý P. ( 2002) On error estimation in the conjugate gradient method and why it works in finite precision computations. Electron. Trans. Numer. Anal. , 13, 56– 80. Strakoš Z. & Tichý P. ( 2005) Error estimation in preconditioned conjugate gradients. BIT , 45, 789– 817. Google Scholar CrossRef Search ADS   Verfürth R. ( 1996) A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques . Advances in Numerical Mathematics Series. Wiley-Teubner. Wu H. & Chen Z. ( 2006) Uniform convergence of multigrid V-cycle on adaptively refined finite element meshes for second order elliptic problems. Sci. China Ser. A , 49, 1405– 1429. Google Scholar CrossRef Search ADS   Xu J. ( 1992) Iterative methods by space decomposition and subspace correction. SIAM Rev. , 34, 581– 613. Google Scholar CrossRef Search ADS   Footnotes 1 Here, we use the notation of our article and consider the estimate for the energy norm of the error. 2 In the setting of this article it means independent of the source term $$f \in L^2({\it {\Omega}})$$. 3 Using the Partial Differential Equation Toolbox. 4 Simple model problems cannot be used for verification of computational efficiency and numerical behaviour of methods and algorithms. Extrapolation of the results observed on simple model problems to difficult practical problems can lead to false conclusions, because some phenomena (such as significant loss of orthogonality and delay of convergence in CG) are on model problems (such as the Poisson boundary value problem) simply not observable. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

On a residual-based a posteriori error estimator for the total error

Loading next page...
 
/lp/ou_press/on-a-residual-based-a-posteriori-error-estimator-for-the-total-error-0rENBYLt0s
Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx037
Publisher site
See Article on Publisher Site

Abstract

Abstract A posteriori error analysis in numerical partial differential equations aims at providing sufficiently accurate information about the distance of the numerically computed approximation to the true solution. Besides estimating the total error, a posteriori analysis should also provide information about its discretization and (inexact) algebraic computation parts. This issue has been addressed by many authors using different approaches. Historically, probably the first and practically very important approach is based on combination of the classical residual-based bound on the discretization error with the adaptive hierarchy of discretizations and computations that allow to incorporate, using various heuristic arguments, the algebraic error. Motivated by some recent publications, this text uses a complementary approach and examines subtleties of the (generalized) residual-based a posteriori error estimator for the total error that rigorously accounts for the algebraic part of the error. The aim is to show on the standard Poisson model problem example, which is used here as a case study, that a rigorous incorporation of the algebraic error represents an intriguing problem that is not yet completely resolved. That should be of concern in $$h$$-adaptivity approaches, where the refinement of the mesh is determined using the residual-based a posteriori error estimator assuming Galerkin orthogonality. The commonly used terminology such as ‘guaranteed computable upper bounds’ should be in the presence of algebraic error cautiously examined. 1. Introduction Historically, most a posteriori analysis in numerical partial differential equations (PDEs) focuses on estimating the discretization error, i.e., on the discrepancy between the solution of the original infinite-dimensional formulation of the problem and the exact solution of its discretized counterpart. This information is crucial for $$h$$-adaptivity, which refines discretization in the parts of the domain where the estimator indicates a large discretization error with the goal of achieving its close-to-uniform spatial distribution over the domain. Estimation of the discretization error, however, has to deal with the principal difficulty: the exact solution of the original problem is unknown, and unless the algebraic computations providing the coordinates of the discrete solution in the discretization basis are performed exactly or with a negligible algebraic error, the exact solution of the discretized problem is also unknown. Near-to-exact algebraic computations can be prohibitive due to extensive computational cost. When solving practical problems, one typically needs to estimate the error even for the computed approximation far from the exact solution of the discretized problem (see, e.g., Nordbotten & Bjørstad, 2008, Conclusions; Keilegavlen & Nordbotten, 2015; Nissen et al., 2015). Reaching exact algebraic results can even be theoretically prohibitive. The matrix eigenvalues, e.g., are (in general) in principle uncomputable by any finite formula as proved by the Abel–Galois theorem, and they can only be approximated iteratively. Moreover, in case of highly nonnormal matrices, there is no guaranteed forward estimate of the accuracy of the computed eigenvalue approximations and we can guarantee the backward error only. As a consequence, because of the inexactness of algebraic computations, the a posteriori error estimates should be based, from their derivation to their application, on the available computed approximations to the solution of the discrete problem. When considering simple model problems, the previous point is seemingly unimportant. The need for solving adaptively large-scale problems, however, requires abandoning the concept of highly accurate algebraic solutions. Incorporating multilevel discretization structure and the associated preconditioned iterative algebraic solvers with reliable stopping criteria is becoming a prerequisite for efficient very large-scale numerical PDE solvers. This has been clearly formulated as a program by Becker et al. (1995); see also other relevant references below. Since then, there is a growing number of work in this direction. The recent survey with many references can be found, e.g., in Arioli et al. (2013a, Section 4); see also the introductory parts of Jiránek et al. (2010), Arioli et al. (2013b), Papež et al. (2016) and (Málek & Strakoš, 2015, Chapter 12). Multilevel discretization structure is often obtained using discretization mesh adaptivity that requires, as mentioned above, local estimation of the discretization error in order to identify mesh elements that need to be refined. For that purpose a standard residual-based a posteriori bound on the discretization error is used (see, e.g., Babuška & Rheinboldt, 1978; Verfürth, 1996, Section 1.2; Ainsworth & Oden, 2000, Section 2.2, Brenner & Scott, 2008, Section 9.2). Since the exact solution of the algebraic problem is not available, the computed approximation that does not satisfy Galerkin orthogonality is often used in the bound, which violates the assumptions under which the bound has been derived. This raises a question to which extent the standard a posteriori error bounds, assuming Galerkin orthogonality can be modified to estimate the total error and, at the same time, to allow comparison of the discretization and algebraic part of the error and, consequently, construction of reliable stopping criteria for algebraic iterative solvers. The standard residual-based a posteriori bound is also a key ingredient of the works that use hierarchy of approximation spaces for estimating the total error and its algebraic part; for the early examples we refer, in particular, to Bramble et al. (1990), Xu (1992), Oswald (1993), Rüde (1993a,b), Becker et al. (1995) and Griebel & Oswald (1995). Therefore, we focus in this article on the residual-based a posteriori bound as a basic building block used elsewhere and investigate its extension for estimating the total error. We do not share the belief that the matter has been fully resolved, apart from simple technicalities. We describe difficulties that still remain open and have to be taken into account in various circumstances. Using, in particular, the articles by Becker & Mao (2009), Carstensen (1999) and Arioli et al. (2013b), we discuss the subtleties one has to deal with while estimating the total and discretization error using a residual-based a posteriori error estimator. We show that removing the standard Galerkin orthogonality assumption, which cannot be used in a large-scale practical application of the bound, requires a nontrivial revision of the known estimator. Even for the simple model problem, the derived extension of the bound contains multiplicative factors that are potentially very large, as shown also in Arioli et al. (2013b), and that cannot be, in general, easily and accurately determined. Moreover, despite the claims published in literature, there exists no a posteriori estimator for the algebraic part of the error that is cheap, easily computable and that gives in practice a tight guaranteed upper bound. As pointed out below, a rigorous a posteriori analysis that incorporates algebraic errors is for realistic problems substantially more difficult than the analysis that assumes exact algebraic computations. In Section 2 we set the notation, discuss some methodological questions and recall several approaches that use a residual-based a posteriori error estimator as a building block. Section 3 recalls the results from Carstensen (1999) on quasi-interpolation. Section 4 presents the revision of the upper bound on the total error from Becker & Mao (2009) and gives its detailed proof that abandons the Galerkin orthogonality assumption. Section 5 comments on the upper bound on the total error presented in Arioli et al. (2013b). Numerical illustrations of the difficulties associated with the multiplicative factors are present in Section 6. Section 7 addresses estimates of the algebraic error and explains misunderstandings present in literature. Finally, the conclusions are drawn. 2. Model problem and comments on methodology We will use the following standard model problem. Let $${\it {\Omega}} \subset \mathbb{R}^2$$ be a polygonal domain (open, bounded and connected set with a polygonal boundary). We consider the Poisson problem with the homogeneous Dirichlet boundary condition   \begin{equation} \text{find }u:{\it {\Omega}} \rightarrow \mathbb{R}\,{:}\qquad- {\it {\Delta}} u = f \quad \textrm{in} \quad {\it {\Omega}}, \qquad u = 0 \quad \textrm{on} \quad \partial {\it {\Omega}}, \label{eq:modelproblem} \end{equation} (2.1) where $$f:{\it {\Omega}} \rightarrow \mathbb{R}$$ is the source term. Hereafter, we use the standard notation for the Sobolev spaces. For $$D \subset {\it {\Omega}}$$, $$L^1(D)$$ denotes the space of the (Lebesgue) integrable functions in $$D$$, $$L^2(D)$$ denotes the space of the square integrable functions in $$D$$, $$(w,v)_D = \int_{D}{v\, w}$$ denotes the $$L^2$$-inner product on $$L^2(D)$$ and $$\left\| {{w}} \right\|_D = (w,w)^{1/2}_D$$ denotes the associated $$L^2$$-norm. We omit the subscripts for $$D = {\it {\Omega}}$$. $$H^k({\it {\Omega}})$$ denotes the Hilbert space of functions in $$L^2({\it {\Omega}})$$ whose weak derivatives up to the order $$k$$ belong to $$L^2({\it {\Omega}})$$. $$H^1_0({\it {\Omega}})$$ denotes the space of functions in $$H^1({\it {\Omega}})$$ with vanishing trace on the boundary $$\partial {\it {\Omega}}$$. Assuming $$f \in L^2({\it {\Omega}})$$, the problem (2.1) can be written in the following weak form   \begin{equation} \text{find } u\in V \equiv H^1_0({\it {\Omega}})\,{:}\qquad(\nabla u, \nabla v) = (\,f,v)\qquad \text{for}\,\text{all}~ v \in V. \label{eq:weakform} \end{equation} (2.2) Let $$\mathcal{T}$$ be a conforming triangulation of the domain $${\it {\Omega}}$$, i.e., two distinct and intersecting elements $$T_1, T_2 \in \mathcal{T}$$ share a common face, edge or vertex. Let $${\mathcal{N}}$$ denote the set of all nodes (i.e., the vertices of the elements of $$\mathcal{T}$$), whereas $${\mathcal{N}_\mathrm{int}} \equiv {\mathcal{N}} \backslash \partial {\it {\Omega}}$$ denotes the set of the free nodes. By $$\mathcal{E}$$ we denote the set of all edges of the elements of $$\mathcal{T}$$ and, similarly, $${\mathcal{E}_\mathrm{int}} \equiv \mathcal{E} \backslash \partial {\it {\Omega}}$$. For any node $$z \in {\mathcal{N}}$$, let $$\varphi_z$$ be the corresponding hat-function, i.e., the piecewise linear function that takes value 1 at the node $$z$$ and vanishes at all other nodes. By $$\omega_z$$ we denote the support of $$\varphi_z$$, which is equal to the patch $$\omega_z = \cup \left\{ T \in \mathcal{T} | z \in T \right\}$$. For an element $$T \in \mathcal{T}$$ we denote $$h_T \equiv {\text{diam}}(T)$$, similarly $$h_z \equiv {\text{diam}}(\omega_z)$$ denotes the diameter of $$\omega_z$$, $$z \in {\mathcal{N}}$$. By $$V_h \subset V$$ we denote the space of the continuous piecewise linear functions on the triangulation $$\mathcal{T}$$ vanishing on the boundary $$\partial {\it {\Omega}}$$, i.e., $$V_h \equiv \textrm{span}\left\{\varphi_z|z \in {\mathcal{N}_\mathrm{int}}\right\}$$. The discrete formulation of (2.2) then reads   \begin{equation} \text{find }u_h \in V_h\,{:}\qquad (\nabla u_h, \nabla v_h) = (\,f,v_h)\qquad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:discreteweakform} \end{equation} (2.3) The solution $$u_h$$ of (2.3) is called the Galerkin solution. Subtracting (2.3) from (2.2) and using $$V_h \subset V$$, we get the Galerkin orthogonality  \begin{equation} (\nabla (u - u_h), \nabla v_h) = 0 \qquad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:galerkinorthogonality} \end{equation} (2.4) The difficulty in estimating the discretization error $$u-u_h$$ mentioned above can be formulated as follows. Consider any estimator $$\text{EST}(\cdot)$$ that provides an upper bound   \begin{equation} \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(u_h), \label{eq:theorbound} \end{equation} (2.5) where $$\left\| {\left| {{\cdot}} \right|} \right\|$$ denotes an appropriate norm (for the model problem (2.1) typically the energy norm $$\left\| {\left| {{w}} \right|} \right\| = \left\| {{\nabla w}} \right\| = (\nabla w, \nabla w)^{1/2}$$). To evaluate the right-hand side of (2.5), we need $$u_h$$ that is not available. The common practice is then replacing $$u_h$$ by the computed approximation $$u_h^C$$, giving the seemingly easy solution   \[ \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(u_h^C). \] This inequality is, however, not guaranteed to hold without further justification that can be highly nontrivial or even impossible to achieve. Provided that   \begin{equation} \text{EST}(u_h) = \inf_{v_h \in V_h} \text{EST}(v_h), \label{eq:boundinfimum} \end{equation} (2.6) the bound (2.5) does indeed lead to a guaranteed upper bound   \begin{equation} \left\| {\left| {{u-u_h}} \right|} \right\| \leq \text{EST}(v_h)\quad \text{for}\,\text{all}~ v_h \in V_h. \label{eq:theorbound2} \end{equation} (2.7) Proving (2.6) can, however, represent a challenge, and the authors are unaware of any applicable results of this kind published in the literature. Another option is simply writing   \begin{equation} \text{EST}(u_h) = \text{EST}(u_h^C) + \big(\text{EST}(u_h)-\text{EST}(u_h^C)\big)\label{eq:estdiff} \end{equation} (2.8) or using a variant of this based on a specific form of the estimator. Then the evaluation of the first term $$\text{EST}(u_h^C)$$ does not require any assumptions (it should not be confused with estimating the total error). Provided that the second term $$\text{EST}(u_h)-\text{EST}(u_h^C)$$ could be bounded using the computed quantity $$u_h^C$$, the relation (2.8) would give a rigorous bound on the discretization error. This consideration has been used in combination with heuristics in various approaches. Before focusing on the residual-based error estimator itself, we recall several ideas from the approaches combining this estimator with the hierarchy of discretizations. In the prototype paper (Becker et al., 1995), the error $$u - v_h$$ of any approximation $$v_h \in H^1_0({\it {\Omega}})$$ is expressed1 as   \begin{align*} \left\| {(u-v_h)} \right\|^2 &= (\,f,u-v_h) - (\nabla v_h, \nabla (u-v_h))\\ &=: \langle r(v_h), u-v_h \rangle, \end{align*} where $$r(v_h) \in H^{-1}({\it {\Omega}})$$ is the associated residual and $$\langle \cdot, \cdot \rangle : H^{-1}({\it {\Omega}}) \times H^1_0({\it {\Omega}}) \to \mathbb{R}$$ represents the duality pairing. Using the hierarchy of meshes and the associated discrete approximation subspaces $$V_0 \subset V_1 \subset \cdots \subset V_L \subset H^1_0({\it {\Omega}})$$, the article considers a multigrid algorithm with the Galerkin projection property and the operators $$I_j$$, $$j = 0,1, \ldots, L$$, where   \[ I_j: H^1_0({\it {\Omega}}) \to V_j \quad j = 0,1,\ldots,L. \] A straightforward substitution then gives   \begin{align} \langle r(v_h), u-v_h \rangle &= \langle r(v_h), (u-v_h) - I_L (u-v_h) \rangle \label{eq:BJR_bound1}\\ \end{align} (2.9)  \begin{align} &\quad{} + \sum^L_{j=1} \langle r(v_h), (I_j - I_{j-1})(u-v_h) \rangle \label{eq:BJR_bound2}\\ \end{align} (2.10)  \begin{align} &\quad{} + \langle r(v_h), I_0(u-v_h) \rangle. \label{eq:BJR_bound3} \end{align} (2.11) The first term (2.9) is then bounded using the standard residual-based a posteriori error estimator with the (non-Galerkin) input quantity $$v_h$$. The second term (2.10) is bounded using the algebraic residuals on the individual levels $$j = L, L-1, \ldots, 1$$. This will bring in nontrivial multiplicative factors analogous to those presented later in our article. Finally, the last term (2.11) is assumed to vanish because of the exact solution of the problem on the coarsest mesh. This assumption is substantial, as demonstrated also by the numerical experiment in Section 7 of the quoted article. The authors also suggest heuristics for stopping criteria in an adaptive algorithm and for approximation of the unknown multiplicative factors. The derivation does not consider roundoff error. There is a large amount of work that in principle can be put into the framework of (2.8), where $$\text{EST}(\cdot)$$ is again the residual-based error estimator and the difference $$(\text{EST}(u_h)-\text{EST}(u_h^C))$$ that reflects the algebraic error is estimated using the hierarchy of splittings of the approximation space $$H^1_0({\it {\Omega}})$$ or of its appropriate discretization; see, e.g., Bramble et al. (1990), Xu (1992), Oswald (1993), Rüde (1993a,b), Griebel & Oswald (1995), Becker et al. (1995) and Harbrecht & Schneider (2016), as well as further references in Arioli et al. (2013a, Section 4) (Málek & Strakoš, 2015, Chapter 12). The instructive article by Stevenson (2007) extends the approach of another remarkable article by Morin et al. (2002) on convergence of adaptive FEM, where the adaptivity is based on the residual-based a posteriori error bound on the discretization error. Stevenson’s rigorously presented results account for inexact algebraic computations, but they show that such extension is indeed highly intriguing. The main result in Stevenson (2007) relies on the continuity argument, i.e., it assumes that the algebraic solver deviates from the exact result in a sufficiently small way. For the algebraic solver, this article refers to the work of Wu & Chen (2006) on uniform convergence of multigrid V-cycle algorithm. The article by Wu & Chen (2006) assumes exact arithmetic and, in particular, the exact solution of the problem on the coarsest mesh. The recalled results underline the importance of understanding the extension of the residual-based a posteriori error bounds to inexact algebraic computations (non-Galerkin solutions). In In Becker & Ma (2009, Lemma 3.1) the bound on the total error is given in the form   \begin{equation} \left\| {{\nabla(u - v_h)}} \right\|^2 \leq \widetilde{C} \cdot\text{EST}(v_h) + 2 \,\left\| {{\nabla(u_h - v_h)}} \right\|^2, \label{eq:BecMaobound} \end{equation} (2.12) where $$\widetilde{C}$$ is stated to depend only on the minimal angle of the triangulation $$\mathcal{T}$$, and $$v_h\in V_h$$ is arbitrary, i.e., it can account for inexact algebraic computations, where $$v_h = u_h^C$$. Here, the first term $$\widetilde{C} \cdot\text{EST}(v_h)$$ represents the first term in (2.8) given by the standard residual-based a posteriori error estimator for the discretization error (with replacing the Galerkin solution $$u_h$$ by $$v_h$$). The second term, $$2 \left\| {{\nabla(u_h - v_h)}} \right\|^2$$ does not have the meaning of the second term in (2.8). The proof refers for the case $$v_h=u_h$$, i.e., for estimating the discretization error, to the article by Carstensen (1999). The proof is completed by arguing, without detailed explanation, that the general case (2.12) follows via the triangle inequality. To be valid and applicable in practical computations as an upper bound, the total error estimator of the form (2.12) must resolve two challenges. (1) It must rigorously justify using the arbitrary (non-Galerkin) $$v_h \in V_h$$ in the first term on the right-hand side of (2.12) and give the necessary information on the value of the factor $$\widetilde{C}$$. (2) It must be proved that in practical computations, we can indeed provide a meaningful (i.e., inexpensive and tight) upper bound on the norm $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ of the algebraic error. In this article, we examine the derivation of (2.12) and prove that an estimator of the analogous form can indeed be used also for approximations $$u_h^C$$ that do not solve the discretized problem exactly. Although the derivation is simple, and in comparison with the standard residual-based a posteriori error estimator assuming Galerkin orthogonality, the resulting bound seemingly only adds the term bounding the algebraic error, the whole matter is, in our opinion, not simple. First, as explained below, there is no proof that the given estimator provides in practice a guaranteed computable upper bound due to subtleties determining the multiplicative factors and due to effects of roundoff on estimating the algebraic error. Second, here we consider a simple model problem and point out difficulties that are not technical but methodological. It is not clear whether for more complicated problems an extension of the estimator preserves the same form or whether an estimator based on the related methodology can meaningfully incorporate algebraic errors. 3. Quasi-interpolation results This section presents results from Carstensen (1999) used further in the text. We include them here for completeness and self-consistency of the text. Denote by $$\psi$$ a piecewise linear function taking value 1 at the inner nodes $$z \in {\mathcal{N}_\mathrm{int}}$$ and vanishing on the boundary $$\partial {\it {\Omega}}$$, $$\psi \equiv \sum_{z\in{\mathcal{N}_\mathrm{int}}} {\varphi_z}$$. Then $$\varphi_z/\psi$$, $$z\in{\mathcal{N}_\mathrm{int}},$$ represents in $${\it {\Omega}}$$ a partition of unity. Indeed, since $$\varphi_z$$, $$z\in{\mathcal{N}_\mathrm{int}}$$, sum up to $$\psi$$, we have $$\sum_{z\in{\mathcal{N}_\mathrm{int}}}{\varphi_z/\psi} = 1$$ in $${\it {\Omega}}$$ (see Carstensen, 1999, Proposition 2.1). The quasi-interpolation operator $${\mathcal{I}}: L^1({\it {\Omega}}) \rightarrow V_h$$ is then defined in the following way. For a given $$w \in L^1({\it {\Omega}})$$,   \[ {\mathcal{I}} w \equiv \sum_{z \in {\mathcal{N}_\mathrm{int}}} {w_z \varphi_z}, \quad \textrm{where} \quad w_z \equiv \frac{(w, \varphi_z/\psi)}{(1, \varphi_z)}. \] The error $$w - {\mathcal{I}} w$$ has a vanishing weighted average. Namely, for $$w$$, $$R \in L^2({\it {\Omega}})$$ and arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{equation} \int_{{\it {\Omega}}} {R \, (w - {\mathcal{I}} w)} = \sum_{z\in {\mathcal{N}_\mathrm{int}}}{\int_{{\it {\Omega}}}{(R - R_z) (w - w_z \psi)(\varphi_z/\psi)}}, \label{eq:weightedaverage} \end{equation} (3.1) see Carstensen (1999, Remark 2.4). Since   \[ {\int_{{\it {\Omega}}}{(w - w_z \psi)(\varphi_z/\psi)}} = 0 \qquad \text{for}\,\text{all}~ z \in {\mathcal{N}_\mathrm{int}}, \] the numbers $$R_z \in \mathbb{R}$$ can be chosen arbitrarily. In particular, $$R_z$$ can be chosen as the mean value of $$R$$ on $$\omega_z$$. Then $$\int_{\omega_z}{|R-R_z|^2}$$ is minimal among all $$R_z \in \mathbb{R}$$. The following lemmas are stated and proved in Carstensen (1999) for a more general case. Considering the model problem (2.1), we restrict ourselves to the case $$w \in H^1_0({\it {\Omega}})$$. The multiplicative factors in the lemmas then depend on (I) The shape of $$\omega_z$$. (II) The shapes of $$\omega_{z\partial{\it {\Omega}}} \equiv (\omega_z \cup \omega_{\xi}\ |\ z \in {\mathcal{N}_\mathrm{int}}, \xi \in {\mathcal{N}} \backslash {\mathcal{N}_\mathrm{int}}, \omega_z \cap \omega_{\xi} \neq 0)$$. (III) The shape coefficients $$(\int_{\omega_z}{\varphi_z}/|\omega_z| \ |\ z\in {\mathcal{N}_\mathrm{int}})$$, where $$|\omega_z|$$ stands for the Lebesgue measure of $$\omega_z$$. (IV) The overlap   \[ M_1 \equiv \max_{z \in {\mathcal{N}_\mathrm{int}}} \mathrm{card} \{ \xi \in {\mathcal{N}} \backslash {\mathcal{N}_\mathrm{int}} \ |\ \omega_z \cap \omega_{\xi} \neq 0\}. \] (V) The shape of the elements $$T \in \mathcal{T}$$. (VI) The value $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\varphi_z}} \right\|_{\infty}}$$, where $$\left\| {{\cdot}_{\infty}} \right\|$$ denotes the $$L^{\infty}({\it {\Omega}})$$-norm and $$h_z = \text{diam}(\omega_z)$$. (VII) The value   \[ M_2 \equiv \underset {x \in {\it {\Omega}}}{\text{ess sup}}{ \{ h(x)/h_T\ |\ x \in T \in \mathcal{T} \} }, \] where $$h(x) \equiv \max\{h_z\ |\ \varphi_z(x) > 0, z \in {\mathcal{N}_\mathrm{int}}\}$$, $$h_T = \text{diam}(T)$$. The proofs of the lemmas use the Poincaré inequality on $$\omega_z$$ and the Friedrichs inequality on $$\omega_{z\partial{\it {\Omega}}}$$ defined in II. To prove Lemma 3.2, the so-called trace theorem (see, e.g., Carstensen, 1999, Proposition 4.1) is used on each element of the triangulation $$T \in \mathcal{T}$$; the multiplicative factor then depends on the shape of the elements; see V. The quantities $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\nabla \varphi_z}} \right\|_{\infty}}$$ and $$M_2$$ (see VI and VII) are of the order 1 on a shape-regular mesh, where $$\left\| {{\nabla \varphi_z|_T}} \right\|_{\infty} \approx h^{-1}_T$$ and $$h_z \approx h_T, T\in \omega_z$$. They deteriorate on a mesh consisting of triangles with small inner angles, where small and large elements (in the sense of their diameter) adjoin. To see the development of the argument, for clarity we present the following two lemmas. Lemma 3.1 (see Carstensen, 1999, Theorem 3.1, statement 1) There exists a multiplicative factor $$C>0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VI), but not on the size of the elements $$h_T$$, such that, $$\text{for}\,\text{all}~ R \in L^2({\it {\Omega}})$$, $$\text{for}\,\text{all}~ w \in H^{1}_0 ({\it {\Omega}})$$ and for arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{eqnarray} \int_{{\it {\Omega}}} {R\, (w - {\mathcal{I}} w)} \leq C \left\| {{\nabla w}} \right\| \left\{ \sum_{z \in {\mathcal{N}_\mathrm{int}}}{h^2_z \int_{\omega_z}{\varphi_z/\psi \ |R-R_z|^2}} \right\}^{1/2}. \nonumber \end{eqnarray} Lemma 3.1 is a consequence of the definition of the quasi-interpolation operator $${\mathcal{I}}$$; see (3.1). Lemma 3.2 (see Carstensen, 1999, Theorem 3.2) Let $$S \subset \mathcal{E}$$. There exists a multiplicative factor $$C > 0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VII), but not on the size of the elements $$h_T$$, such that $$\text{for}\,\text{all}~ J \in L^2(S)$$ and $$\text{for}\,\text{all}~ w \in H^1_0({\it {\Omega}})$$,   \[ \int_S{J \, (w - {\mathcal{I}} w)} \leq C \left\| {{\nabla w}} \right\| \left( \sum_{T\in\mathcal{T}}{h_T \left\| {{J}} \right\|^2_{S \cap \partial T}} \right)^{1/2}. \] Combining Lemmas 3.1 and 3.2 we get the final inequality. Lemma 3.3 (see Carstensen, 1999, Corollary 3.1) Let $$S \subset \mathcal{E}$$. There exists a multiplicative factor $$C_1>0$$ depending on I–VII such that, for all $$J \in L^2(S)$$, $$\text{for}\,\text{all}~ R \in L^2({\it {\Omega}})$$, for all $$w \in H^{1}_0 ({\it {\Omega}})$$ and for arbitrary numbers $$R_z \in \mathbb{R}$$, $$z \in {\mathcal{N}_\mathrm{int}}$$,   \begin{eqnarray*} &&\int_{{\it {\Omega}}} {R\, (w - {\mathcal{I}} w)} + \int_{S} {J\, (w - {\mathcal{I}} w)} \\ &&\qquad \leq C_1 \left\| {{\nabla w}} \right\| \left\{ \sum_{z \in {\mathcal{N}_\mathrm{int}}}{h^2_z \left\| {{R - R_z}} \right\|^2_{\omega_z}} + \sum_{T \in \mathcal{T}} {h_T \left\| {{J}} \right\|^2_{S \cap \partial T}} \right\}^{1/2}. \end{eqnarray*} The following lemma introduces a positive multiplicative factor $$C_{\mathrm{intp}}$$ that plays a key role in our discussion on incorporating the algebraic error into the a posteriori bound on the total error; see Section 4 and the numerical experiments in Section 6. Lemma 3.4 (see Carstensen, 1999, Theorem 3.1, statement 3) There exists a multiplicative factor $$C_{\mathrm{intp}} > 0$$ depending on the triangulation $$\mathcal{T}$$ (more precisely on I–VI and VI) such that, for all $$\text{for}\,\text{all}~ w \in H^{1}_0({\it {\Omega}})$$,   \begin{equation} \left\| {{{\mathcal{I}} w}} \right\| \leq C_{\mathrm{intp}} \left\| {{\nabla w}} \right\| . \label{eq:Cintp} \end{equation} (3.2) Remark 3.5 The factor $$C_{\mathrm{intp}}$$ satisfies   \[ \sup_{w \in H^{1}_0({\it {\Omega}})} \frac{\left\| {{{\mathcal{I}} w}} \right\|}{\left\| {{\nabla w}} \right\|} \leq C_{\mathrm{intp}}. \] Obtaining a value of $$C_{\mathrm{intp}}$$ such that the above inequality is tight represents a nontrivial issue. Using the proof of Carstensen (2006, Theorem 2.4) and the discussion in Carstensen (2006, Example 2.3),, we can get a better idea about the size of $$C_{\mathrm{intp}}$$. For a shape-regular mesh with $$\max_{z \in {\mathcal{N}}}{h_z \left\| {{\nabla \varphi_z}} \right\|_{\infty}} \approx 2$$ (see VI), there holds $$C_{\mathrm{intp}} \approx 6$$. In general, as stated in Carstensen (2006), it may be very large for small angles in the triangulation. For $$f \in L^1({\it {\Omega}})$$, define the mean-value operator $$\pi_{\omega_z}(\,f) \equiv \int_{\omega_z}{f}/|\omega_z|$$. We denote, for $$z \in {\mathcal{N}}$$ and for any subset $$Z \subset {\mathcal{N}}$$,   \[ {\text{osc}}_z \equiv |\omega_z|^{1/2} \left\| {{ f - \pi_{\omega_z} f}} \right\|_{\omega_z}, \qquad {\text{osc}}(Z) \equiv \left( \sum_{z \in Z}{{\text{osc}}^2_z} \right)^{1/2}, \] measuring the oscillations of $$f$$, i.e., the variations of the function $$f$$ from the mean value $$\pi_{\omega_z} f$$ on the subdomains $$\omega_z$$. Given $$v_h \in V_h$$, define for $$E \in {\mathcal{E}_\mathrm{int}}$$ and any subset $$F \subset {\mathcal{E}_\mathrm{int}}$$ the edge residuals   \[ J_E(v_h) \equiv |E|^{1/2} \left\| \left[ \frac{\partial v_h}{\partial n_E} \right] \right\|_{E}, \qquad J(v_h, F) \equiv \left( \sum_{E \in F}{J^2_E(v_h)} \right)^{1/2}, \] where $$[\cdot]$$ denotes the jump of a piecewise continuous function and $$n_E$$ denotes the unit normal to $$E$$ (for each $$E \in {\mathcal{E}_\mathrm{int}}$$ the orientation of the unit normal is set arbitrarily but fixed). The edge residual $$J_E(v_h)$$, $$v_h \in V_h,$$ measures the jump of the piecewise constant gradient $$\nabla v_h$$ over the inner edge $$E$$. We set for brevity $${\text{osc}} \equiv {\text{osc}}({\mathcal{N}})$$ and $$J(v_h) \equiv J(v_h, {\mathcal{E}_\mathrm{int}})$$. For a given $$v_h \in V_h$$, we define the jump function $$\mathcal{J}(v_h) \in L^2({\mathcal{E}_\mathrm{int}})$$ on the inner edges   \begin{equation} \mathcal{J}(v_h)|_E \equiv \left[ \frac{\partial v_h}{\partial n_E} \right]\quad E \in {\mathcal{E}_\mathrm{int}}. \end{equation} (3.3) The Green’s formula (see, e.g., Ciarlet, 2002, p. 14) gives for a domain $$D$$ with a Lipschitz continuous boundary $$\partial D$$ and for $$v \in H^2(D)$$, $$w \in H^1(D)$$  \begin{equation} \int_{D} \nabla v \cdot \nabla w = - \int_{D} {\it {\Delta}} v \, w + \int_{\partial D} \left(\frac{\partial v}{\partial n_{\partial D}} \right)\, w, \label{eq:Green} \end{equation} (3.4) where $$n_{\partial D}$$ denotes the unit normal to $$\partial D$$ pointing outward. Let $$v_h \in V_h$$ and $$T \in \mathcal{T}$$. Then $$v_h |_T$$ is a linear function, $$v_h |_T \in H^2(T)$$ and $${\it {\Delta}} v_h|_T = 0$$. Then, applying the Green’s formula (3.4) elementwise yields, for any $$v_h \in V_h$$, $$w \in H^1_0({\it {\Omega}})$$,   \begin{eqnarray} \int_{{\it {\Omega}}}{\nabla v_h \cdot \nabla w} &=& \sum_{T \in \mathcal{T}}{\int_T{\nabla v_h \cdot \nabla w}} = \sum_{T \in \mathcal{T}} \left( - {\int_T{{\it {\Delta}} v_h \, w}} + {\int_{\partial T}{\left(\frac{\partial v_h}{\partial n_{\partial T}} \right)\, w} } \right) \nonumber \\ &=& \sum_{E \in {\mathcal{E}_\mathrm{int}}}{\int_E {\left[ \frac{\partial v_h}{\partial n_E} \right] w}} = \int_{{\mathcal{E}_\mathrm{int}}}{\mathcal{J}(v_h)\, w} . \label{eq:perpartes} \end{eqnarray} (3.5) The results recalled in this section are used to prove the upper bound on the total error in the next section. 4. Upper bound on the total error Now we state the upper bound on the energy norm of the total error using the residual-based a posteriori error estimator. Theorem 4.1 There exist triangulation-dependent positive multiplicative factors $$C_1$$, $$C_{\mathrm{intp}}$$ and $$C_2$$ such that for the solution $$u$$ of (2.2), the Galerkin solution $$u_h$$ of (2.3) and an arbitrary $$v_h \in V_h$$,   \begin{equation} \label{eq:lemma31} \left\| {{\nabla (u - v_h)}} \right\|^2 \leq 2\, C^2_1\, C^2_2 \left( J^2(v_h) + {\text{osc}}^2 \right) + 2\, C^2_{\mathrm{intp}} \left\| {{ \nabla (u_h - v_h) }} \right\|^2. \end{equation} (4.1) In particular, $$C_1$$ depends on I–VII (see Lemma 3.3), $$C_{\mathrm{intp}}$$ depends on I–VI and VI (see Lemma 3.4) and the factor $$C_2$$ depends on the ratios $$h_z^2/|\omega_z|$$, $$z \in {\mathcal{N}_\mathrm{int}},$$ and $$h_T/|E|$$, $$T\in\mathcal{T}$$, $$E\in\partial T\cap {\mathcal{E}_\mathrm{int}}$$. Proof. We will use the standard expression for the norm   \begin{equation} \left\| {{\nabla (u - v_h) }} \right\| = \sup_{0 \neq w \in H^1_0({\it {\Omega}})} {\frac{1}{\left\| {{\nabla w}} \right\|} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w }}. \label{eq:dualnorm} \end{equation} (4.2) Let $$v_h \in V_h$$ and $$w \in H^1_0({\it {\Omega}}), w \neq 0,$$ be arbitrary,   \begin{eqnarray*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla {\mathcal{I}} w}\\ &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u - u_h) \cdot \nabla {\mathcal{I}} w} \\ && + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{eqnarray*} It follows from the definition of the interpolation operator that $${\mathcal{I}} w \in V_h$$. The Galerkin orthogonality (2.4) gives   \[ \int_{{\it {\Omega}}}{\nabla(u - u_h) \cdot \nabla {\mathcal{I}} w} = 0. \] Then   \begin{eqnarray*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &=& \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla(w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w} \\ &=& \int_{{\it {\Omega}}}{\nabla u \cdot \nabla(w - {\mathcal{I}} w)} - \int_{{\it {\Omega}}}{\nabla v_h \cdot \nabla(w - {\mathcal{I}} w)}\\ && + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{eqnarray*} Using the weak formulation (2.2) and the equality (3.5),   \[ \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} = \int_{{\it {\Omega}}}{f ( w - {\mathcal{I}} w)} - \int_{{\mathcal{E}_\mathrm{int}}}{\mathcal{J}(v_h) ( w - {\mathcal{I}} w)} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \] Then Lemma 3.3 with $$S = {\mathcal{E}_\mathrm{int}}$$, $$R = f$$, $$R_z = \pi_{\omega_z}f,\ z \in {\mathcal{N}_\mathrm{int}}$$, $$J = -\mathcal{J}(v_h)$$ gives   \begin{align*} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} &\leq C_1 \left\| {{\nabla w}} \right\| \left\{ \sum_{T \in \mathcal{T}} {h_T \left\| {{\mathcal{J}(v_h)}} \right\|^2_{{\mathcal{E}_\mathrm{int}} \cap \partial T}} + \sum_{z \in {\mathcal{N}_\mathrm{int}}}{\! h^2_z \left\| {{f - \pi_{\omega_z}f}} \right\|^2_{\omega_z}} \right\}^{\frac12} \\ & \quad + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}\\ &\leq C_1 C_2 \left\| {{\nabla w}} \right\| \left\{ \sum_{E \in {\mathcal{E}_\mathrm{int}}} {|E| \, \left\| {{ \mathcal{J}(v_h) }} \right\|^2_{E}} + \sum_{z \in {\mathcal{N}_\mathrm{int}}}{|\omega_z|\, \left\| {{f - \pi_{\omega_z}f}} \right\|^2_{\omega_z}} \right\}^{\frac12} \\ & \quad + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}\\ &= C_1 C_2 \left\| {{\nabla w}} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \int_{{\it {\Omega}}}{\nabla(u_h - v_h) \cdot \nabla {\mathcal{I}} w}. \end{align*} Using the Cauchy–Schwarz inequality,   \begin{equation} \int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w} \leq C_1 C_2 \left\| {{\nabla w}} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \left\| {{\nabla {\mathcal{I}} w}} \right\| \, \left\| {{ \nabla (u_h - v_h)}} \right\|. \label{eq:proofI} \end{equation} (4.3) Dividing (4.3) by $$\left\| {{\nabla w}} \right\|$$ and using Lemma 3.4,   \begin{eqnarray*} \frac{1}{\left\| {{\nabla w}} \right\|}\int_{{\it {\Omega}}}{\nabla(u - v_h) \cdot \nabla w } &\leq& C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \frac{\left\| {{\nabla {\mathcal{I}} w}} \right\|}{\left\| {{\nabla w}} \right\|} \left\| {{\nabla (u_h - v_h)}} \right\| \\ &\leq& C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + C_{\mathrm{intp}} \left\| {{\nabla(u_h - v_h)}} \right\|. \end{eqnarray*} Using the representation (4.2) of the energy norm, recall that $$w \in H^1_0({\it {\Omega}})$$ was chosen arbitrarily,   \begin{equation} \label{eq:totalboundCintp} \left\| {{\nabla (u - v_h)}} \right\| \leq C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + C_{\mathrm{intp}} \left\| {{ \nabla (u_h - v_h) }} \right\|. \end{equation} (4.4) Finally, we deduce (4.1) using the inequality $$(a + b)^2 \leq 2a^2 + 2b^2$$. The multiplicative factor $$C_{\mathrm{intp}}$$ is independent of the solution $$u \in H^1_0$$ of (2.2)2 and of its approximation $$v_h \in V_h$$. We will now give another bound on the total error that will contain information about the exact solution $$u$$ of (2.2). It is useful for illustration of the role of the factor $$C_{\mathrm{intp}}$$ in (4.1). Setting $$w \equiv u - v_h \in H^1_0({\it {\Omega}})$$ in (4.3), we get   \[ \left\| {{\nabla (u - v_h)}} \right\|^2 \leq C_1 C_2 \left\| {{ \nabla(u - v_h) }} \right\| \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + {\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|} \, \left\| {{\nabla (u_h - v_h) }} \right\| . \] Dividing both sides by $$\left\| {{ \nabla(u - v_h) }} \right\|$$ gives   \begin{equation} \left\| {{ \nabla(u - v_h) }} \right\| \leq C_1 C_2 \left( J^2(v_h) + {\text{osc}}^2 \right)^{1/2} + \frac {\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla (u - v_h) }} \right\|} \ \left\| {{ \nabla(u_h - v_h) }} \right\|, \label{eq:lemma31Cinttilde} \end{equation} (4.5) which is formulated as the following corollary. Corollary 4.2 Using the notation of Theorem 4.1,   \begin{equation} \label{eq:lemma31Ctilde} \left\| {{\nabla(u - v_h)}} \right\|^2 \leq 2\, C^2_1\, C^2_2 \left( J^2(v_h) + {\text{osc}}^2 \right) + 2\, \widetilde{C}^{\,2}_{\mathrm{intp}}(u,v_h) \,\left\| {{ \nabla(u_h - v_h) }} \right\|^2, \end{equation} (4.6) where   \begin{equation} \widetilde{C}_{\mathrm{intp}}(u, v_h) \equiv \frac{\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla (u - v_h) }} \right\|}\leq C_{\mathrm{intp}}. \label{eq:Cint} \end{equation} (4.7) Since $$C_{\mathrm{intp}}$$ is independent of $$u \in H^1_0$$ (independent of the source term $$f$$) and of $$v_h \in V_h$$, we must have   \begin{equation} \sup_{u\in H^1_0,\ u \text{ solves (2.2)},\ f \in L^2({\it {\Omega}})}\hspace{.5em} \sup_{v_h \in V_h} \frac{\left\| {{ \nabla({\mathcal{I}} u - {\mathcal{I}} v_h) }} \right\|}{\left\| {{\nabla(u - v_h) }} \right\|}\leq C_{\mathrm{intp}}. \label{eq:worstcase} \end{equation} (4.8) Therefore, $$C_{\mathrm{intp}}$$ represents a worst-case scenario factor and one may expect that most likely $$\widetilde{C}_{\mathrm{intp}}(u,v_h) \ll C_{\mathrm{intp}}$$. 5. A related result In Arioli et al. (2013b), the authors consider elliptic self-adjoint problems, and they use a residual-based error estimator for setting the stopping criterion for the conjugate gradient method (CG). Following their approach, one can easily get a theoretical upper bound on the total error. Although the bound is not stated explicitly in Arioli et al. (2013b), it appears in the proof of Theorem 3.3; see the inequality (Arioli et al., 2013b, (3.22)). The derivation proceeds differently from the proof of the bound (4.1); it again demonstrates that the price to be paid for removing the assumption on exact algebraic computations in terms of including unknown and possibly large multiplicative factors can be high. Moreover, in contrast to the statements in Arioli et al. (2013b), the resulting estimator cannot be considered a guaranteed practical upper bound due to the difficulties in estimating the algebraic part of the error; see Section 7 below. First, Arioli et al. (2013b, Theorem 2.2) recalls the bound on the discretization error: there exists a multiplicative factor $$C_{2.2}>0$$ that is dependent on the minimal angle of the triangulation $$\mathcal{T}$$, but independent of $$h$$, $$u$$ and $$u_h$$, such that   \begin{equation} \|{\nabla (u - u_h) }\|^2 \leq C_{2.2}\, \eta^2(u_h), \qquad \eta(u_h) \equiv \left( \sum_{T\in{\mathcal{T}}} |T| \, \|{ f + {\it{\Delta}} u_h}\|^2_T + (J(u_h))^2 \right)^{1/2},\label{eq:resbounddiscr} \end{equation} (5.1) see, e.g., Verfürth (1996, Section 1.2) and Ainsworth & Oden (2000, Section 2.2). Using inverse estimates for piecewise polynomial functions and Young’s inequality, Arioli et al. (2013b, Lemma 3.1) yields the inequality   \[ \eta^2(w_h) \leq (1+\gamma)\, \eta^2(v_h) + C_{3.1} (1+\gamma^{-1}) \left\| {{\nabla (v_h - w_h) }} \right\|^2 \quad \text{for}~\text{all}~ w_h, v_h \in V_h,\ \gamma > 0, \] where the positive factor $$C_{3.1}$$ depends on the minimal angle of the triangulation $$\mathcal{T}$$. Combining these bounds and the equality   \[ \left\| {{\nabla (u - v_h) }} \right\|^2 = \left\| {{\nabla (u - u_h) }} \right\|^2 + \left\| {{ \nabla(u_h - v_h) }} \right\|^2 \] that follows from the Galerkin orthogonality (2.4), we get the upper bound on the total error   \begin{eqnarray*} \left\| {{\nabla (u - v_h) }} \right\|^2 &\leq& C_{2.2}\, \eta^2(u_h) + \left\| {{\nabla (u_h - v_h) }} \right\|^2 \\ &\leq& C_{2.2} \, (1+\gamma)\, \eta^2(v_h) + \left(1 + C_{2.2}\, C_{3.1} \, (1+\gamma^{-1}) \right) \left\| {{\nabla (u_h - v_h) }} \right\|^2. \end{eqnarray*} Finally, by setting $$\gamma \equiv 1$$,   \begin{equation} \left\| {{ \nabla(u - v_h) }} \right\|^2 \leq 2\, C_{2.2} \, \eta^2(v_h) + \left(1 + 2\, C_{2.2}\, C_{3.1}\right) \left\| {{\nabla (u_h - v_h) }} \right\|^2. \label{eq:Ariolietalbound} \end{equation} (5.2) Comparing (5.1) with (5.2), we see that replacing the Galerkin solution $$u_h$$ by an arbitrary $$v_h \in V_h$$ results in the additional term $$\left\| {{ \nabla(u_h - v_h) }} \right\|^2$$, which is, however, multiplied by an unknown and potentially very large multiplicative factor $$(1 + 2\, C_{2.2}\, C_{3.1})$$. In the criteria proposed in Arioli et al. (2013b, Section 5) for numerical experiments, the factors are empirically set to $$C_{2.2} \equiv 40$$, $$C_{3.1} \equiv 10$$, giving $$(1 + 2\, C_{2.2}\, C_{3.1}) = 801$$; cf. (2.12). This underlines the subtleties of the residual-based bounds discussed above. In addition, as mentioned above and explained in detail in Section 7 below, getting a tight practical upper bound on $$\left\| {{\nabla (u_h - v_h) }} \right\|^2$$ represents an unresolved challenge. 6. Numerical illustrations We use, on purpose, very simple problems to illustrate the possible difference in the values of $$\widetilde{C}_{\mathrm{intp}}(u,v_h)$$ and $$C_{\mathrm{intp}}$$. Although $$\widetilde{C}_{\mathrm{intp}}(u,v_h)$$ can be, assuming the knowledge of the exact solution $$u$$, evaluated up to a negligible quadrature error, for the factor $$C_{\mathrm{intp}}$$, we present a lower bound given by plugging a chosen function into (3.2). The derivation of a more accurate estimate for $$C_{\mathrm{intp}}$$ (see also the discussion in Remark 3.5) is beyond the scope of this article. 6.1 Numerical illustration in one dimension We first consider a one-dimensional analog of the Clément-type quasi-interpolation operator $${\mathcal{I}}$$ to illustrate that $$C_{\mathrm{intp}}$$ can be significantly larger than one. Consider the domain $${\it {\Omega}} = (0,1)$$ with the (nonuniform) partition   \[ [0, \beta, {1}/{3} \pm \beta, {2}/{3} \pm \beta, 1-\beta, 1]\quad \beta = 0.01. \] This partition is adapted to the one-dimensional Laplace problem with the solution   \begin{eqnarray} \label{eq:atan-solution} u(x) &=& {\tan^{-1}}(cx) - {\tan^{-1}}(c(x-{1}/{3})) + {\tan^{-1}}(c(x-{2}/{3})) \\ && \qquad - {\tan^{-1}}(c(x-1)) - {\tan^{-1}}({c}/{3}) + {\tan^{-1}}({2c}/{3}) - {\tan^{-1}}(c), \nonumber \end{eqnarray} (6.1) with $$c = 1000$$. The left part of Fig. 1 depicts the solution $$u$$ and the Clément-type quasi-interpolant $${\mathcal{I}} u$$. For a zero approximate vector, we have Fig. 1. View largeDownload slide Left: the solution $$u$$ (6.1) (solid line) and the interpolant $${\mathcal{I}} u$$ (dotted line). Right: the function $$w(x) = x(1-x)$$ (solid line) and $${\mathcal{I}} w$$ (dotted line). Fig. 1. View largeDownload slide Left: the solution $$u$$ (6.1) (solid line) and the interpolant $${\mathcal{I}} u$$ (dotted line). Right: the function $$w(x) = x(1-x)$$ (solid line) and $${\mathcal{I}} w$$ (dotted line).   \begin{equation} \widetilde{C}_{\mathrm{intp}}(u,0) = \frac{\left\| {{({\mathcal{I}} u)'}} \right\|}{\left\| {{u'}} \right\|} = 0.77. \label{eq:Ctilde} \end{equation} (6.2) For a quadratic function $$w(x) = x(1-x)$$, $$w\in H^1_0({\it {\Omega}}),$$ we have   \[ \frac{\left\| {{({\mathcal{I}} w)'}} \right\|}{\left\| {{w'}} \right\|} = 3.70; \] see the right part of Fig. 1 for the plot of the function $$w$$ and the interpolant $${\mathcal{I}} w$$. Consequently, $$C_{\mathrm{intp}} \geq 3.70$$. 6.2 Two-dimensional numerical illustration For two-dimensional numerical illustration, we consider the square domain $${\it {\Omega}} \equiv (-1,1) \times (-1,1)$$ and the triangulation $$\mathcal{T}$$ generated by MATLAB3 command initmesh (‘squareg’, ‘Hmax’, 0.1)+ that provides a Delaunay triangulation consisting of 1 368 elements with the maximal diameter less than or equal to $$0.1$$. The minimal angle of the mesh is 35.9° and the average of the minimal angles of the elements is 50.3°. Consider the solution of problem (2.1):   \begin{equation} u^{(1)}(x,y) = (x-1)(x+1)(y-1)(y+1). \label{eq:u1} \end{equation} (6.3) For the zero approximate solution and the Galerkin solution $$u^{(1)}_h$$ corresponding to $$u^{(1)}$$, we have   \[ \widetilde{C}_{\mathrm{intp}}(u^{(1)},0) = 1.02, \qquad \widetilde{C}_{\mathrm{intp}}(u^{(1)},u^{(1)}_h) = 0.16. \] Similarly, for the exact solution   \begin{equation} u^{(2)}(x,y) = (x-1)(x+1)(y-1)(y+1)\exp(-100(x^2 + y^2)), \label{eq:u2} \end{equation} (6.4) we have   \[ \widetilde{C}_{\mathrm{intp}}(u^{(2)},0) = 0.76, \qquad \widetilde{C}_{\mathrm{intp}}(u^{(2)}, u^{(2)}_h) = 0.28. \] In Fig. 2 we show the values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)},v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$ defined above. Fig. 2. View largeDownload slide The values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)}, v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$; see (6.3) and (6.4), respectively. Fig. 2. View largeDownload slide The values of $$\widetilde{C}_{\mathrm{intp}}(u^{(j)}, v_h)$$ for $$v_h$$ generated in CG iterations with zero initial vector for solving the linear algebraic systems corresponding to the discretization of (2.2) with the solutions $$u^{(1)}$$, $$u^{(2)}$$; see (6.3) and (6.4), respectively. To bound the constant $$C_{\mathrm{intp}}$$ from below, we consider $$w_h \in V_h$$ such that   \begin{equation} w_h(z) = 1 \quad z \in {\mathcal{N}_\mathrm{int}}, \qquad w_h = 0 \quad \text{on } \partial {\it {\Omega}}. \label{eq:wh} \end{equation} (6.5) For this function   \[ 1.10 = \frac{\left\| {{\nabla {\mathcal{I}} w_h }} \right\|}{\left\| {{\nabla w_h }} \right\|} \leq C_{\mathrm{intp}}. \] Figure 3 gives the difference $$w_h-{\mathcal{I}} w_h$$. This is a piecewise linear function that is on the machine precision level in most of the domain except patches around the inner nodes adjacent to the boundary $$\partial {\it {\Omega}}$$. We recall that for this simple problem and a shape-regular mesh $$C_{\mathrm{intp}} \approx 6$$, see Remark 3.5 and the original article by Carstensen (2006). It can therefore indeed hold $$C_{\mathrm{intp}} \gg \widetilde{C}_{\mathrm{intp}}(u, v_h)$$. Fig. 3. View largeDownload slide The difference $$w_h-{\mathcal{I}} w_h$$ for $$w_h$$ given by (6.5). Fig. 3. View largeDownload slide The difference $$w_h-{\mathcal{I}} w_h$$ for $$w_h$$ given by (6.5). 7. Algebraic error estimates Various published articles that include algebraic errors in a posteriori error analysis consider preconditioned iterative methods (such as CG) for solving the associate symmetric positive definite algebraic problem. Estimating the algebraic error $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ is then often considered to be resolved, apart from some seemingly technical issues such as the choice of a number of extra CG iterations, or the approximation of the smallest eigenvalue of the system matrix (see, e.g., Arioli et al., 2013b, Section 4), where the presented Gauss–Radau quadrature-based bound is considered “the only guaranteed upper bound for the $$A$$-norm of the CG error.” Here, we present mathematical arguments that challenge this opinion. In short, as above, theoretical results with assumptions that do not hold in practical computations cannot be applied without an appropriate theoretical justification, here the numerical stability analysis. In practice, the Gauss–Radau quadrature-based estimate does not give, in general, a guaranteed and tight upper bound. The Gauss–Radau quadrature requires an approximation of the smallest eigenvalue of the system matrix from below. A close approximation to the smallest eigenvalue from above can be in theory obtained in the process of (exact arithmetic) CG computations. A closer investigation reveals, however, a serious difficulty. If the approximation to the smallest eigenvalue is not accurate enough, the estimate for the norm of the algebraic error can be inaccurate. If, on the other hand, the eigenvalue approximation approaches the smallest eigenvalue, then the computation of the Gauss–Radau estimator must be cautiously done in a numerically stable way, since it (implicitly) inverts a matrix that becomes in such case close to numerically singular. Moreover, if a second approximation of the smallest eigenvalue is formed due to loss of linear independence of the computed vectors generating the associated Krylov subspaces, then the Gauss–Radau quadrature estimates exhibit instabilities that are not yet understood. We can see again the generally valid principle: In evaluation of a posteriori error estimates, all sources of errors, including roundoff, must be taken into consideration. This is also illustrated by the following point, which is valid even if the difficulty with a tight approximation of the smallest eigenvalue from below is resolved. In numerical computations, we cannot guarantee that Gauss–Radau quadrature estimates give an upper bound due to effects of roundoff errors to the underlying short recurrences in CG computations. This is a principal issue as explained next. Because of the loss of orthogonality caused by roundoff errors, the formulas that express the error norms in CG iterations under the assumption of exact arithmetic and, consequently, orthogonality among the computed basis vectors are no longer valid. In practical computations, we neither have orthogonal bases nor we have Krylov subspaces in the strict mathematical meaning. Unless the problem is really computationally simple (as in the Poisson model problem),4 orthogonality and also linear independence among the vectors computed using short recurrences is, in general, ratherquickly lost and convergence of the computed iterations is substantially delayed. In terms of the computed Jacobi matrices that are substantial in the derivation of the Gauss–Radau quadrature bounds, their entries will quickly become far away (orders of magnitude) from their theoretical counterparts. Recent description of the related issues can be found in the article by Gergelits & Strakoš(2014) with the references to many earlier articles (see, e.g., Paige, 1980; Greenbaum, 1989; Strakoš, 1991; Greenbaum & Strakoš, 1992; Notay, 1993; Strakoš & Tichý, 2002; Strakoš & Tichý, 2005; Meurant & Strakoš, 2006; O’Leary et al., 2007). The matter has been comprehensively discussed already in the survey article by Strakoš & Liesen (2005) and it is also covered in the recent monograph (Liesen & Strakoš, 2013, Section 5.9). In summary, the following principal question is omitted in works assuming exact arithmetic: Considering the effects of roundoff recalled above, how do we know that the formulas derived under the assumptions that are so drastically violated give anything meaningful in practical computations? This fundamental issue cannot be resolved by heuristic arguments; it requires thorough analysis. The question on the effects of roundoff errors to error estimation in iterative methods such as CG has already been raised in the seminal article by Dahlquist et al. (1979). It has been investigated in Section 5 (called ‘Rounding error analysis’) of the article by Golub & Strakoš (1994), and again very thoroughly (in the context of different estimates) in the article by Strakoš & Tichý (2002) with the title ‘On error estimation in the conjugate gradient method and why it works in finite precision computations’; see also the survey article by Meurant & Strakoš (2006). This underlines the point. Without thorough rounding error analysis, we may say that we have observed some behaviour on some examples and nothing more. With thorough rounding error analysis, we can explain why and under which conditions some estimates work and prove them numerically safe. Perhaps even more important, we can prove that other (mathematically equivalent) estimates can behave in a numerically unstable way and they should not be used. Estimates that have been proved numerically unstable (see, e.g., Strakoš & Tichž, 2002) are unfortunately indeed frequently used in practice. The facts presented, e.g., in Golub & Strakoš (1994, Section 5), Strakoš & Tichý(2002, Section 7), Meurant (2006, Chapter 7), Meurant & Strakoš (2006, Section 5), O’Leary et al. (2007) and Liesen & Strakoš (2013, Section 5.9), show that without a thorough and rigorous analysis, application of the error bounds derived under the assumption of exact computation, to the results of finite precision computations, is not only methodologically wrong, but it can also indeed lead to false conclusions. A partial theoretical justification of the Gauss and Gauss–Radau quadrature estimates is provided (based on the earlier works of Paige, Greenbaum and others) in Section 5 of the article by Golub & Strakoš (1994). In short, justification of the quadrature-based bounds must be based on the Riemann–Stieltjes integration using the distribution function, with possibly many more points of increase than the original distribution function determined by the data of the problem; see also the associated arguments in the article by O’Leary et al. (2007) on sensitivity of the Gauss quadrature. The article by Golub & Strakoš (1994) mentioned earlier justifies using Gauss and Gauss–Radau quadrature estimates in finite precision computations (with limitations specified in the article). The estimates based on the Gauss–Radau quadrature can be useful, but they cannot be proved to give a tight guaranteed upper bound on the error norms. Summarizing, accurate and computationally efficient estimation of the algebraic error $$\left\| {{\nabla(u_h - u_h^C)}} \right\|$$ still represents a challenge. This challenge is not resolved by deriving estimators assuming exact computations and then plugging in the computed quantities. The fact that some estimators can be used in the same form for input entries computed using finite precision arithmetic is not at all obvious and it cannot be guessed a priori by any heuristics. This fact can only result from a rigorous analysis that is (in these cases) highly nontrivial. 8. Conclusion The presented article investigates changes in derivation and application of the standard residual-based a posteriori error bound on the discretization error needed to get an estimator (not necessarily a practically applicable guaranteed upper bound) for the total error. Technically, this article provides a detailed proof of the residual-based upper bound on the total approximation error that requires knowledge of the associated multiplicative factors. As published previously in Becker & Mao (2009) and Arioli et al. (2013b), abandoning the Galerkin orthogonality assumption in the derivation leads naturally to an additional term accounting for the algebraic part of the error. We show that the contribution of the algebraic error is scaled by a nontrivial multiplicative factor; see (4.1) and (4.5)–(4.7). This multiplicative factor $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$ depends, besides the computed approximation $$v_h$$, also on the unknown infinite-dimensional solution $$u$$ of (2.2). Therefore, it is generally uncomputable. It can be bounded, using the a priori information, by the factor $$C_{\mathrm{intp}}$$ independent of $$u$$ and $$v_h$$ given by (3.2). The value of $$C_{\mathrm{intp}}$$ can therefore overestimate the value of $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$. Moreover, as another substantial point, a guaranteed computationally efficient upper bound on the algebraic part of the error that provides, in general, sufficiently tight results is not available. The main message of this article does not concern technical details about multiplicative factors in the residual-type a posteriori error estimators for the total error developed for the given model problem. It shows that handling inexact algebraic computations still needs, despite many remarkable results, further work. In practical computations, we cannot avoid using various heuristics. This article supports using heuristics (they are used in many articles coauthored by the authors of this article). It points out, however, a need for supporting analysis. It shows that even for the simple model problem and the standard residual-based a posteriori error estimator, the matter is not easy and the unresolved questions can be practically important. As, e.g., numerically illustrated in in Papež (2016, Section 4.2) and as mentioned in the Introduction, application of the residual-based error estimator for the mesh refinement adaptivity remains in the presence of algebraic errors an open problem. When the standard estimator for the discretization error is evaluated using computed quantities, there is no guarantee that its local contributions provide a meaningful indication of the spatial distribution of the discretization error over the domain. The matter is practically very important because it can affect efficiency of $$h$$-adaptive computations. Extrapolation of the observations obtained for simple model problems cannot be considered the final and generally valid justification. The fact that in (4.1) and (4.6) the algebraic error is estimated globally, and the multiplicative factors $$C_{\mathrm{intp}}$$ and $$\widetilde{C}_{\mathrm{intp}}(u, v_h)$$ are not easy to estimate, suggests that the matter is intriguing and it requires further work. This article cannot survey the previously published and recently developed results towards robust stopping criteria that balance the algebraic and discretization error. This is addressed, e.g., in the works using the hierarchy of subspace splittings recalled in Section 2 (see also Huber et al., 2017) or in Arioli et al. (2013a, Section 4), Jirànek et al. (2010), Papež et al. (2014), Carstensen et al. (2014, Section 7.1), and Papež et al. (2016). In Papež et al. (2016), it is shown, using the methodology based on flux reconstruction, that such stopping criteria could indeed be constructed and rigorously supported by analysis. It also shows, however, that in practical applications, there is a substantial computational cost to be paid, and this cost can even become in some cases excessive. There is still a work to be done. As shown at the example of the residual-based a posteriori error estimator in this article, such work should go hand in hand with analysis. Practical importance of Theorem 4.1 is not in application of the bound (4.1). It is in the warning that such application can be difficult and unjustified heuristics can be misleading. In particular, any use of the adjective ‘guaranteed’ should carefully examine the assumptions under which this adjective holds in practice. Acknowledgements The authors are grateful to David Silvester, Endre Süli and an anonymous referee for valuable comments, which improved the text. Funding ERC-CZ project LL1202 financed by the Ministry of Education, Youth and Sports of the Czech Republic. References Ainsworth M. & Oden J. T. ( 2000) A Posteriori Error Estimation in Finite Element Analysis . Pure and Applied Mathematics (New York). New York: Wiley-Interscience [John Wiley & Sons] pp. xx+ 240. Google Scholar CrossRef Search ADS   Arioli M., Liesen J., Miȩdlar A. & Strakoš Z. ( 2013a) Interplay between discretization and algebraic computation in adaptive numerical solution of elliptic PDE problems. GAMM-Mitt. , 36, 102– 129. Google Scholar CrossRef Search ADS   Arioli M., Georgoulis E. H. & Loghin D. ( 2013b) Stopping criteria for adaptive finite element solvers. SIAM J. Sci. Comput. , 35, A1537– A1559. Google Scholar CrossRef Search ADS   Babuška I. & Rheinboldt W. C. ( 1978) Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. , 15, 736– 754. Google Scholar CrossRef Search ADS   Becker R. & Mao S. ( 2009) Convergence and quasi-optimal complexity of a simple adaptive finite element method. M2AN Math. Model. Numer. Anal. , 43, 1203– 1219. Google Scholar CrossRef Search ADS   Becker R., Johnson C. & Rannacher R. ( 1995) Adaptive error control for multigrid finite element methods. Computing , 55, 271– 288. Google Scholar CrossRef Search ADS   Bramble J. H., Pasciak J. E. & Xu J. ( 1990) Parallel multilevel preconditioners. Numerical Analysis 1989 (Dundee, 1989) . Pitman Res. Notes Math. Ser., vol. 228. Harlow: Longman Sci. Tech., pp. 23– 39. Google Scholar CrossRef Search ADS   Brenner S. C. & Scott L. R. ( 2008) The Mathematical Theory of Finite Element Methods . Texts in Applied Mathematics, vol. 15, 3rd edn. New York: Springer, pp. xviii+ 397. Google Scholar CrossRef Search ADS   Carstensen C. ( 1999) Quasi-interpolation and a posteriori error analysis in finite element methods. M2AN Math. Model. Numer. Anal. , 33, 1187– 1202. Google Scholar CrossRef Search ADS   Carstensen C. ( 2006) Clément interpolation and its role in adaptive finite element error control. Partial Differential Equations and Functional Analysis . Oper. Theory Adv. Appl., vol. 168. Basel: Birkhäuser, pp. 27– 43. Google Scholar CrossRef Search ADS   Carstensen C., Feischl M., Page M. & Praetorius D. ( 2014) Axioms of adaptivity. Comput. Math. Appl. , 67, 1195– 1253. Google Scholar CrossRef Search ADS PubMed  Ciarlet P. G. ( 2002) The Finite Element Method for Elliptic Problems . Classics in Applied Mathematics, vol. 40. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xxviii+ 530 (Reprint of the 1978 original [North-Holland, Amsterdam]). Google Scholar CrossRef Search ADS   Dahlquist G., Golub G. H. & Nash S. G. ( 1979) Bounds for the error in linear systems. Semi-infinite programming (Proc. Workshop, Bad Honnef, 1978) . Lecture Notes in Control and Information Science vol. 15. Berlin: Springer, pp. 154– 172. Google Scholar CrossRef Search ADS   Gergelits T. & Strakoš Z. ( 2014) Composite convergence bounds based on Chebyshev polynomials and finite precision conjugate gradient computations. Numer. Algorithms , 65, 759– 782. Google Scholar CrossRef Search ADS   Golub G. H. & Strakoš Z. ( 1994) Estimates in quadratic formulas. Numer. Algorithms , 8, 241– 268. Google Scholar CrossRef Search ADS   Greenbaum A. ( 1989) Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences. Linear Algebra Appl. , 113, 7– 63. Google Scholar CrossRef Search ADS   Greenbaum A. & Strakoš Z. ( 1992) Predicting the behavior of finite precision Lanczos and conjugate gradient computations. SIAM J. Matrix Anal. Appl. , 13, 121– 137. Google Scholar CrossRef Search ADS   Griebel M. & Oswald P. ( 1995) On the abstract theory of additive and multiplicative Schwarz algorithms. Numer. Math. , 70, 163– 180. Google Scholar CrossRef Search ADS   Harbrecht H. & Schneider R. ( 2016) A note on multilevel based error estimation. Comput. Methods Appl. Math. , 16, 447– 458. Google Scholar CrossRef Search ADS   Huber M., Rüde U., Strakoš Z. & Wohlmuth B. ( 2017) Error estimators for highly parallel multigrid solver (in preparation). Jiránek P., Strakoš Z. & Vohralík M. ( 2010) A posteriori error estimates including algebraic error and stopping criteria for iterative solvers. SIAM J. Sci. Comput. , 32, 1567– 1590. Google Scholar CrossRef Search ADS   Keilegavlen E. & Nordbotten J. M. ( 2015) Inexact linear solvers for control volume discretizations in porous media. Comput. Geosci. , 19, 159– 176. Google Scholar CrossRef Search ADS   Liesen J. & Strakoš Z. ( 2013) Krylov Subspace Methods: Principles and Analysis . Numerical Mathematics and Scientific Computation. Oxford: Oxford University Press. Málek J. & Strakoš Z. ( 2015) Preconditioning and the Conjugate Gradient Method in the Context of Solving PDEs.  SIAM Spotlights, vol. 1. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. x+ 104. Meurant G. ( 2006) The Lanczos and Conjugate Gradient Algorithms. From Theory to Finite Precision Computations.  Software, Environments, and Tools, vol. 19. Philadelphia, PA: SIAM, pp. xvi+ 365. Google Scholar CrossRef Search ADS   Meurant G. & Strakoš Z. ( 2006) The Lanczos and conjugate gradient algorithms in finite precision arithmetic. Acta Numer. , 15, 471– 542. Google Scholar CrossRef Search ADS   Morin P., Nochetto R. H. & Siebert K. G. ( 2002) Convergence of adaptive finite element methods. SIAM Rev. , 44, 631– 658 (2003) (Revised reprint of ‘Data oscillation and convergence of adaptive FEM’ [SIAM J. Numer. Anal., 38 (2000), no. 2, 466–488]). Google Scholar CrossRef Search ADS   Nissen A., Pettersson P., Keilegavlen E. & Nordbotten J. M. ( 2015) Incorporating geological uncertainty in error control for linear solvers. SPE Reservoir Simulation Symposium . Society of Petroleum Engineers. DOI: http://doi.org/10.2118/173272-MS Nordbotten J. M. & Bjørstad P. E. ( 2008) On the relationship between the multiscale finite-volume method and domain decomposition preconditioners. Comput. Geosci. , 12, 367– 376. Google Scholar CrossRef Search ADS   Notay Y. ( 1993) On the convergence rate of the conjugate gradients in presence of rounding errors. Numer. Math. , 65, 301– 317. Google Scholar CrossRef Search ADS   O’Leary D. P., Strakoš Z. & Tichý P. ( 2007) On sensitivity of Gauss-Christoffel quadrature. Numer. Math. , 107, 147– 174. Google Scholar CrossRef Search ADS   Oswald P. ( 1993) On a BPX-preconditioner for $$\rm P1$$ elements. Computing , 51, 125– 133. Google Scholar CrossRef Search ADS   Paige C. C. ( 1980) Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem. Linear Algebra Appl. , 34, 235– 258. Google Scholar CrossRef Search ADS   Papež J. ( 2016) Algebraic error in matrix computations in the context of numerical solution of partial differential equations. Ph.D. thesis . Prague: Charles University. Papež J., Liesen J. & Strakoš Z. ( 2014) Distribution of the discretization and algebraic error in numerical solution of partial differential equations. Linear Algebra Appl. , 449, 89– 114. Google Scholar CrossRef Search ADS   Papež J., Strakoí Z. & Vohralík M. ( 2016) Estimating and localizing the algebraic and total numerical errors using flux reconstructions. Preprint MORE/2016/12 (accepted for publication in Numer. Math.). Rüde U. ( 1993a) Fully adaptive multigrid methods. SIAM J. Numer. Anal. , 30, 230– 248. Google Scholar CrossRef Search ADS   Rüde U. ( 1993b) Mathematical and Computational Techniques for Multilevel Adaptive Methods . Frontiers in Applied Mathematics, Vol. 13. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), pp. xii+ 140. Google Scholar CrossRef Search ADS   Stevenson R. ( 2007) Optimality of a standard adaptive finite element method. Found. Comput. Math. , 7, 245– 269. Google Scholar CrossRef Search ADS   Strakoš Z. ( 1991) On the real convergence rate of the conjugate gradient method. Linear Algebra Appl. , 154–156, 535– 549. Google Scholar CrossRef Search ADS   Strakoš Z. & Liesen J. ( 2005) On numerical stability in large scale linear algebraic computations. ZAMM Z. Angew. Math. Mech. , 85, 307– 325. Google Scholar CrossRef Search ADS   Strakoš Z. & Tichý P. ( 2002) On error estimation in the conjugate gradient method and why it works in finite precision computations. Electron. Trans. Numer. Anal. , 13, 56– 80. Strakoš Z. & Tichý P. ( 2005) Error estimation in preconditioned conjugate gradients. BIT , 45, 789– 817. Google Scholar CrossRef Search ADS   Verfürth R. ( 1996) A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques . Advances in Numerical Mathematics Series. Wiley-Teubner. Wu H. & Chen Z. ( 2006) Uniform convergence of multigrid V-cycle on adaptively refined finite element meshes for second order elliptic problems. Sci. China Ser. A , 49, 1405– 1429. Google Scholar CrossRef Search ADS   Xu J. ( 1992) Iterative methods by space decomposition and subspace correction. SIAM Rev. , 34, 581– 613. Google Scholar CrossRef Search ADS   Footnotes 1 Here, we use the notation of our article and consider the estimate for the energy norm of the error. 2 In the setting of this article it means independent of the source term $$f \in L^2({\it {\Omega}})$$. 3 Using the Partial Differential Equation Toolbox. 4 Simple model problems cannot be used for verification of computational efficiency and numerical behaviour of methods and algorithms. Extrapolation of the results observed on simple model problems to difficult practical problems can lead to false conclusions, because some phenomena (such as significant loss of orthogonality and delay of convergence in CG) are on model problems (such as the Poisson boundary value problem) simply not observable. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Sep 11, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off