Exact a posteriori error analysis of the least squares finite element method

Exact a posteriori error analysis of the least

squares ®nite element method

1

Jinn-Liang Liu

Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan

Abstract

A residual type a posteriori error estimator is presented for the least squares ®nite element method. The estimator is proved to equal the exact error in a norm induced by some least squares functional. The error indicator of each element is equal to the exact error norm restricted to the element as well. In other words, the estimator is perfectly e€ective and reliable for error control and for adaptive mesh re®nement. The exactness property requires virtually no assumptions on the regularity of the solution and on the ®nite element order in the approximation or in the estimation. The least squares method is in a very general setting that applies to various linear boundary-value problems such as the elliptic systems of ®rst-order and of even-order and the mixed type partial dif-ferential equations. Numerical results are given to demonstrate the exactness. Ó 2000 Elsevier Science Inc. All rights reserved.

Keywords: Exact error estimator; Least squares ®nite elements

1. Introduction

A posteriori error estimation is one of the most important components in adaptive methodology [19,22,28]. It seems that there are no exact a posteriori error analyses available for all the ®nite element, ®nite volume, and boundary element methods in the literature. An estimator for the least squares ®nite

www.elsevier.com/locate/amc

*_{E-mail address: jinnliu@math.nctu.edu.tw; http://www.math.nctu.edu.tw/~jinliu (J.-L. Liu).} 1_{This work was supported by NSC under grant NSC87-2115-M-009-005, Taiwan.}

0096-3003/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S0096-3003(99)00153-8

element method (LSFEM) is presented here and is proved to be not only globally but also locally exact.

Since the pioneering work of Babuska and Rheinboldt [6], it has been an important subject to study the reliability of error estimators for adaptive nu-merical methods in various applications. The reliability is usually quanti®ed by the so-called e€ectivity index h de®ned as

h ˆ_keke ;

where e denotes an error estimator and the denominator is the exact error e measured in a suitable norm k
k. Theoretically as well as practically, the quantity is often desired to be either bounded below and above by some positive constants, say C1 and C2, which are independent of the mesh

param-eter h of the approximation, i.e., C16 h 6 C2

or asymptotically exact, i.e., lim

h!0h ˆ 1

or both, see e.g., [1,2,6,7,10,15,16,21,24,26,27].

The error estimator described here is for LSFEM in a very general setting. It is obtained by calculating the residual in L2 _{norm and is proved not only}

globally exact, i.e., h ˆ 1;

but also locally exact, i.e., hi:ˆ_kekei

iˆ 1;

where ei denotes the residual norm (an error indicator) of an element i in a

particular triangulation of the domain associated with the parameter h and kek_i is the exact error norm restricted to that element. Here, the norm k
k is in-duced by the bilinear form derived from a least squares functional which in turn is associated with a given boundary value problem. In other words, this is a perfect error estimator if the error of numerical integration is not taken into account. Even more surprisingly, the proof of the exactness in the context of least squares formulation is almost trivial when compared with the previous a posteriori error analysis. Furthermore, no extra assumptions on the regularity of the solution or on the ®nite element spaces used in the approximation are required for this exactness.

The residual type error estimation is, of course, not new [6,15,18,19,22,24] and is presented in various formulations for various applications. The present analysis shows that the exactness of the estimator may be an additional out-standing feature for LSFEM which has been recognized as an attractive method in many applications in recent years, see e.g., [3,5,8,9,11±14,18,20].

2. An exact estimator

The boundary value problems considered herein are in a very general setting and are expressed in the form

Lju ˆ fj; j ˆ 1; 2; . . . ; m; in X; …1†

Bku ˆ gk; k ˆ 1; 2; . . . ; n; on oX; …2†

where X  Rd_{, d ˆ 1; 2; 3, is a bounded domain with the boundary oX,}

u ˆ …u1; . . . ; um†t, Lj are linear di€erential operators and Bj are linear

boundary operators. We always assume that the problem (1) and (2) has a unique solution u in some (m-tuple) Cartesian product of function spaces de-noted by H…X† with the given functions fj2 L2…X†, gk 2 L2…oX†. The setting is

so general that it applies to a wide class of linear boundary value problems such as the elliptic systems of Refs. [3,25], the ®rst-order systems of [8,9,11±14] and the mixed type PDEs of [4,5,17,23] in the ®rst-order formulation, etc.

The problem (1) and (2) will be approximated by the LSFEM. The method itself will also be presented in a very general setting by this we mean that the associated least squares functional contains equation residuals and boundary residuals all weighted by positive parameters. It is convenient to weight (1) and (2) before de®ning the functional. Let (1) and (2) be weighted as follows:

Lwu ˆ fw; in iX; …3† Bwu ˆ gw; on oX; …4† where Lw:ˆ …l1L1; . . . ; lmLm†t; Bw:ˆ …b1B1; . . . ; bnBn†t; fw:ˆ …l1f1; . . . ; lmfm†t; gw:ˆ …b1g1; . . . ; bngn†t;

and lj and bk are the weighting parameters.

The functional, denoted by J: H…X† ! R, is de®ned as J…v† ˆXm jˆ1 Z iX l2 j…Ljv ÿ fj†2dX ‡ Xn kˆ1 Z oX b2 k…Bkv ÿ gk†2ds ˆ: kLwv ÿ fwk20;X‡ kBwv ÿ gwk20;oX: …5†

Here, the norms k
k_0;X and k
k_0;oX are associated with the inner products …
;
†_0;Xand …
;
†_0;oX, respectively. Evidently, the solution u 2 H…X† of (1) and (2) minimizes the functional and vice versa, i.e.,

J…u† ˆ min

Taking the ®rst variation, the solution equivalently satis®es the equation

B…u; v† ˆ F …v† 8v 2 H…X†; …6†

where

B…u; v† :ˆ …Lwu; Lwv†0;X‡ …Bwu; Bwv†0;oX;

F …v† :ˆ …fw; Lwv†0;X‡ …gw; Bwv†0;oX:

The LSFEM for (1) and (2) is to solve (6) in a ®nite element subspace Shof

H…X†, that is, to ®nd uh2 Shsuch that

B…uh; vh† ˆ F …vh† 8vh2 Sh; …7†

where h is the mesh size of a regular triangulation Th of X. Note that if the

boundary condition (2) is computed exactly, i.e., if the ®nite element space Sh

consists of functions that satisfy (2), the weighting parameters bk are

indi€er-ent. Hence, the least squares functional (5) leads to a very general setting of LSFEMs with weighted or non-weighted di€erential and/or boundary residu-als. It thus applies to all the above mentioned least squares problems.

Evidently, the uniqueness assumption of the solution of (1) and (2) implies that the bilinear form B…
;
† induces a norm by which we denote k
k_B. For most ®rst-order systems, the norm is equivalent to the H1 _{norm [8,11±14]. It}

also implies the following result for which a proof can be found in [5]. Theorem 1. If the problem (1) and (2) has a unique solution u 2 H…X†, then there exists a unique function uh2 Sh H…X† satisfying Eq. (7).

Once the approximate solution uhis available, one of the major concerns in

practice is to assess the reliability of this approximation, i.e., to estimate the exact error u ÿ uh in some suitable norm. A posteriori error estimators

repre-sent an important means towards the assessment. For LSFE solution, we in-troduce a residual error estimator e de®ned as

e2_{ˆ kL}

wuhÿ fwk20;X‡ kBwuhÿ gwk20;oX: …8†

The estimator is readily computable. In fact, if we compute the residual in each element si of the current triangulation Th, we obtain an error indicator eifor

that element, i.e., e2

i ˆ kLwuhÿ fwk0;s2 i‡ kBwuhÿ gwk

2

0;oX\si: …9†

The square of the estimator is thus the sum of all squares of the indicators. De®ne the e€ectivity index

h ˆ e

We also de®ne local e€ectivity indices by hiˆ_kekei

B;si

; …11†

where kek_B;si is the restriction of kek_B on si.

Theorem 2. Let u 2 H…X† and uh 2 Sh H…X† be the solutions of (6) and (7),

respectively. Then h ˆ hiˆ 1 …12† for all si2 Th. Proof. e2 i ˆ kLwuhÿ fwk0;s2 i‡ kBwuhÿ gwk 2 0;oX\si ˆ kLwuhÿ Lwuk0;s2 i‡ kBwuhÿ Bwuk 2 0;oX\si

ˆ …Lwe; Lwe†0;si‡ …Bwe; Bwe†0;oX\si

ˆ B…e; e†j_si:

Therefore, we have the exactness of the local e€ectivity indices. The global exactness is again immediate. 

Several remarks are in order.

Remark 1. The global as well as local exactness of (12) appears to be ®rst presented in the literature. The proof of this statement is almost trivial when compared with the previous a posteriori error analysis, see e.g. [1,2,6,7,10, 16,21,24,26,27]. Furthermore, there are virtually no assumptions for this exactness, i.e., no saturation assumptions like those of [7,21], no extra regularity assumptions on the exact solution like those of [1,2,6,10,16,24,26] and no restrictions on the ®nite element orders used in the approximation. In fact, we even do not require the approximate solutions uh of (7) to be

convergent at all (cf. Theorem 1).

Remark 2. It is well-known that the LSFEM inherently provides very attractive properties in applications. For example, the trial and test functions are not required to satisfy the boundary conditions, its discretization results in symmetric and positive de®nite algebraic system, a single piecewise polynomial ®nite element space may be used for all test and trial functions and it does not require the inf±sup condition to be satis®ed when compared with the mixed ®nite element method etc., see loc. cit. The exactness of the error estimator may yet provide an additional outstanding feature of the LSFEM versus other methods, since it is perfectly reliable and e€ective.

Remark 3. Obviously, the implementation of the residual estimator is simple. The exactness is also preserved in terms of the numerical integration. The error indicators can be computed on parallel processors without any communication cost since no jump terms across element boundaries and no local boundary conditions are involved. Together with the symmetric property of the algebraic system, the entire adaptive procedure of least squares computations can be parallel and distributed if a conjugate gradient solver is used since, in this case, there is no need for a global assembly and the iterative process can be done locally [18].

3. Numerical example

By introducing the vorticity x ˆ ov=ox ÿ ou=oy, the 2-D dimensionless Stokes equations can be written as [13]

Lu ˆ 0 0 mo=oy o=ox 0 0 ÿ mo=ox o=oy o=ox o=oy 0 0 o=oy ÿ o=ox 1 0 2 6 6 6 4 3 7 7 7 5 u v x p 2 6 6 6 4 3 7 7 7 5ˆ f1 f2 0 0 2 6 6 6 4 3 7 7 7 5ˆ f in X; …13† where u and v are the x and y components of velocity, p the total head of pressure, m the inverse of Reynolds number and f1 and f2 the given body

forces. We consider a model problem [8] of (13) to which the exact solution is given by

u1ˆ …x h ÿ a†2_{‡ …y ÿ b†}2ic=2_; u2ˆ …x h ÿ a†2_{‡ …y ÿ b†}2ic=2_; x ˆ …xh ÿ a†2‡ …y ÿ b†2ic=2; p ˆ …xh ÿ a†2_{‡ …y ÿ b†}2ic=2_;

in a unit square f…x; y† j 0 6 x 6 1; 0 6 y 6 1g, where a ˆ b ˆ 0:1234 and c ˆ 0:9. Singularity appears at the point …a; b†. The appropriate boundary condition (2) can be constructed via the exact solution. For the least squares

Table 1

NN kekB hmin hmax h

21 0.4416 1 1 1 30 0.3549 1 1 1 39 0.0962 1 1 1 137 0.0456 1 1 1 165 0.0313 1 1 1 243 0.0221 1 1 1 289 0.0180 1 1 1 294 0.0174 1 1 1 353 0.0131 1 1 1 823 0.0054 1 1 1 1071 0.0038 1 1 1 1522 0.0026 1 1 1 2382 0.0016 1 1 1

approximation (7), we impose the boundary condition in biquadratic ®nite element spaces Shand set all weighting parameters to one. An adaptive process

using the residual error estimator (8) begins with the initial mesh Fig. 1 and ends with the ®nal mesh Fig. 2. All of the (global and local) e€ectivity indices of the adaptive computation are equal to one and are shown in Table 1 where NN denotes the number of nodes, hminˆ minihiand hmaxˆ maxihi.

References

[1] S. Adjerid, J.E. Flaherty, Y.J. Wang, A posteriori error estimation with ®nite element methods of lines for 1-D parabolic systems, Numer. Math. 65 (1993) 1±21.

[2] M. Ainsworth, J.T. Oden, A uni®ed approach to a posteriori error estimation using element residual methods, Numer. Math. 65 (1993) 23±50.

[3] A.K. Aziz, R.B. Kellogg, A.B. Stephens, Least squares methods for elliptic systems, Math. Comput. 44 (1985) 53±70.

[4] A.K. Aziz, S. Leventhal, Finite lement approximation for ®rst order systems, SIAM J. Numer. Anal. 6 (1978) 1103±1111.

[5] A.K. Aziz, J.-L. Liu, A weighted least squares method for the backward-forward heat equation, SIAM J. Numer. Anal. 28 (1991) 156±167.

[6] I. Babuska, W.C. Rheinboldt, Error estimates for adaptive ®nite element computations, SIAM J. Numer. Anal. 15 (1978) 736±754.

[7] R.E. Bank, A. Weiser, Some a posteriori error estimators for elliptic partial di€erential equations, Math. Comput. 44 (1985) 283±301.

[8] P.B. Bochev, M.D. Gunzburger, Analysis of least-squares ®nite element methods for the Stokes equations, Math. Comput. 63 (1994) 479±506.

[9] Z. Cai, R. Lazarov, T.A. Manteu€el, S.F. McCormick, First-order system least squares for second-order partial di€erential equations: Part I, SIAM J. Numer. Anal. 31 (1994) 1785±1799. [10] C. Carstensen, E.P. Stephan, A posteriori error estimates for boundary element method, Math.

Comput. 64 (1995) 483±500.

[11] C.L. Chang, Finite element approximation for grad-div type systems in the plane, SIAM J. Numer. Anal. 29 (1992) 452±461.

[12] C.L. Chang, An error estimate of the least squares ®nite element method for the Stokes problem in 3-D, Math. Comput. 63 (1994) 41±50.

[13] C.L. Chang, B.-N. Jiang, An error analysis of least-squares ®nite element method of velocity± pressure±vorticity formulation for Stokes problem, Comput. Meth. Appl. Mech. Engrg. 84 (1990) 247±255.

[14] C.L. Chang, S.-Y. Yang, C.-H. Hsu, A least-squares ®nite element method for incompressible ¯ow in stress±velocity±pressure version, Comput. Meth. Appl. Mech. Engrg. 128 (1995) 1±9. [15] L. Demkowicz, J.T. Oden, M. Ainsworth, P. Geng, Solution of elastic scattering problems in linear acoustics using h±p boundary element method, J. Comput. Appl. Math. 36 (1991) 29±63. [16] R. Duran, R. Rodriguez, On the asymptotic exactness of Bank±Weiser's estimator, Numer.

Math. 62 (1992) 297±303.

[17] K.O. Friedrichs, Symmetric positive di€erential equations, Comm. Pure Appl. Math. 11 (1958) 333±418.

[18] B.-N. Jiang, G.F. Carey, Adaptive re®nement for least-squares ®nite elements with element-by-element conjugate gradient solution, Int. J. Numer. Meth. Engrg. 24 (1987) 569±580. [19] J.-L. Liu, On weak residual error estimation, SIAM J. Sci. Comput. 14 (1996) 1249±1268. [20] J.-L. Liu, I.-J. Lin, M.-Z. Shih, R.-C. Chen, M.-C. Hsieh, Object oriented programming of

[21] J.-L. Liu, W.C. Rheinboldt, A posteriori ®nite element error estimators for inde®nite elliptic boundary value problems, Numer. Func. Anal. Optimiz. 15 (1994) 335±356.

[22] J.T. Oden, L. Demkowicz, W. Rachowicz, T.A. Westermann, Toward a universal h±p adaptive ®nite element strategy, Part 2. A posteriori error estimation, Comput. Meth. Appl. Mech. Engrg. 77 (1989) 113±180.

[23] P. Sermer, R. Mathon, Least-squares methods for mixed-type equations, SIAM J. Numer. Anal. 18 (1981) 705±723.

[24] R. Verfurth, A posteriori error estimates for the Stokes equations, Numer. Math. 55 (1989) 309±325.

[25] W.L. Wendland, Elliptic Systems in the Plane, Pitman, London, 1979.

[26] W.L. Wendland, D. Yu, Adaptive boundary element methods for strongly elliptic integral equations, Numer. Math. 53 (1988) 539±558.

[27] J.Z. Zhu, O.C. Zienkiewicz, Superconvergence recovery technique and a posteriori error estimators, Int. J. Numer. Meth. Engrg. 30 (1990) 1321±1339.

[28] O.C. Zienkiewicz, Computational mechanics today, Int. J. Numer. Meth. Engrg. 34 (1992) 9±33.