Boundary Integral Equations and A Posteriori Error Estimates

TSINGHUA SCIENCE AND TECHNOLOGY ISSN 1007-0214 05/18 pp35-42 Volume 10, Number 1, February 2005

Boundary Integral Equations and A Posteriori Error Estimates * YU Dehao (余德浩)**,

ZHAO Longhua (赵龙花)

Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and System Science, Chinese Academy of Sciences, P. O. Box 2719, Beijing 100080, China Abstract:

Adaptive methods have been rapidly developed and applied in many fields of scientific and engi-

neering computing. Reliable and efficient a posteriori error estimates play key roles for both adaptive finite element and boundary element methods. The aim of this paper is to develop a posteriori error estimates for boundary element methods. The standard a posteriori error estimates for boundary element methods are obtained from the classical boundary integral equations. This paper presents hyper-singular a posteriori error estimates based on the hyper-singular integral equations. Three kinds of residuals are used as the estimates for boundary element errors. The theoretical analysis and numerical examples show that the hypersingular residuals are good a posteriori error indicators in many adaptive boundary element computations. Key words:

boundary integral equation; natural boundary reduction; a posteriori error estimate; hypersingular residual; pseudo-differential operator

Introduction In recent years, adaptive techniques with a posteriori error estimates have been widely used in many fields of scientific and engineering computing. These techniques have been applied to the finite element methods first by Babuska and his coworkers[1-4]. Then a posteriori error estimates were also applied to the adaptive boundary element method[5-7]. Various reliable and efficient a posteriori error estimates play key roles in both adaptive finite element and boundary element methods. There are many ways to obtain a posteriori error estimates with residuals usually used as a posteriori error estimates. The residual is defined as the difference between two sides of the boundary integral Received: 2004-08-02

﹡ Supported by the National Key Basic Research and Development (973) Program of China (No. G19990328) and the Knowledge Innovation Program of the Chinese Academy of Sciences

﹡﹡ To whom correspondence should be addressed. E-mail: [email protected]; Tel: 86-10-62553888

equations when the exact solution is replaced by the approximate solution. Since boundary value problems have many kinds of boundary integral equations, different residuals can be defined for each problem. This paper introduces boundary reduction for the potential problem with the comparison of the residuals for Dirichlet, Neumann, and mixed problems analyzed using boundary element methods.

1

Boundary Integral Equations

Scientific and engineering problems are described by mathematical formulations such as partial differential equations inside a domain or integral equations on a boundary. These formulations are equivalent in theory, but usually have different computational effects. A typical problem is the potential problem described by the Poisson equation in some domain Ω which is a boundary value problem of an elliptic partial differential equation. Let E(P, Q) be the fundamental solution

36

Tsinghua Science and Technology, February 2005, 10(1): 35–42

of the Laplace equation and Γ be the boundary of Ω. The classical boundary integral equation for the potential problem is φ 0 = Kφ 0 + V φ n (1) where the linear operators are defined as: Vu( P) = 2 Ku( P) = −2

∫Γ E ( P, Q)u(Q)ds

Q,

∂ E ( P, Q)u(Q)dsQ . ∂nQ

∫Γ

A boundary integral equation with a hyper-singular operator is φ n = Lφ 0 + M φ n (2) where Mu( P) = 2

∂

∫Γ ∂n

E ( P, Q)u(Q)dsQ ,

P

Lu( P) = −2

∂

∫Γ ∂n

P

∂ E ( P, Q)u(Q )dsQ . ∂nQ

Here V, K, M, and L are all pseudo-differential operators with integer orders. The weakly-singular operator V of order −1 is a smoothing operator. For smooth domains, K and M are compact. Therefore, I−K and I−M are operators of order 0, where I is the identity operator. The hyper-singular operator L is of order +1. Let G(P,Q) be the Green function for the problem. Then by the theory of natural boundary reduction given by Yu[8,9], the natural boundary integral equation is φ n = Sφ 0 (3)

this case, Eq. (1) is a Fredholm integral equation of the first kind, Eq. (2) is a singular integral equation of the Cauchy type which is a Fredholm integral equation of the second kind, and Eq. (3) is a hyper-singular integral which eliminates the need to solve any equations. In Neumann boundary value problems, φ n is given, and φ 0 should be found. Then Eq. (1) is a Fredholm integral equation of the second kind and Eqs. (2) and (3) are hyper-singular integral equations. Even though Eqs. (2) and (3) are both hyper-singular integral equations, Eq. (3) is much simpler than Eq. (2) which simplifies theoretical analyses and numerical computations. Furthermore, Eq. (3) can be coupled with the finite element method directly, so it plays an important role in the domain decomposition methods[8,9]. In mixed boundary value problems, φ n1 is to be found on Γ 1 where φ 01 is given and φ 02 is to be found on Γ 2 where φ n 2 is given, and Γ= Γ1+ Γ2. Define the operators: Vij u( P) = 2 K ij u( P) = −2 M ij u( P) = 2 Lij u( P) = −2

∫Γ

∫Γ

∫Γ

∫Γ j

j

j

E ( P, Q)u(Q)dsQ , j

∂ E ( P, Q)u(Q)dsQ , ∂nQ ∂ E ( P, Q)u(Q )dsQ , ∂nP

∂ ∂ E ( P, Q)u(Q )dsQ , ∂nP ∂nQ

for P ∈ Γ i . Similarly, we can also define Sij , i, j=1,2:

where Su( P) = −

∂

∫Γ ∂n

P

∂ G ( P, Q )u(Q)dsQ， ∂nQ

Sij u( P) = −

∫Γ

∂ ∂ G ( P, Q)u(Q)dsQ j ∂nP ∂nQ

which relates the Dirichlet and Neumann conditions. The hyper-singular integral operator S is the Dirichlet to Neumann (DtN) operator, which is a pseudodifferential operator of order +1. Particularly, on a circular boundary with radius R, the natural boundary integral operator for harmonic problems is 1 *: H α (Γ ) → H α −1 (Γ ) S =− (4) θ 2 4πR sin 2 α where H (Γ) is a Sobolev space on the boundary Γ. For the Dirichlet boundary value problems, φ 0 is

for P ∈ Γ i . By using Eq. (1) ,

given while φ n is unknown and should be found. In

This hyper-singular integral equation has the same defect as Eq. (5).

⎡ V11 ⎢ −V ⎣ 21

K12 ⎤ ⎡φ n1 ⎤ ⎡ I − K11 ⎢ ⎥= I − K 22 ⎥⎦ ⎣⎢φ 02 ⎦⎥ ⎢⎣ K 21

−V12 ⎤ ⎡φ 01 ⎤ ⎢ ⎥ V22 ⎥⎦ ⎣⎢φ n 2 ⎦⎥

(5)

which is a Cauchy type integral equation. The coefficient matrix of the linear system is not symmetric positive definite. By using Eq. (2) , M12 ⎤ ⎡φ 01 ⎤ ⎡ I − M11 − L12 ⎤ ⎡φ n1 ⎤ ⎡ L11 = ⎢ ⎥ ⎢ ⎥ (6) ⎢ −M − L22 ⎥⎦ ⎢⎣φ 02 ⎥⎦ ⎢⎣ − L21 M 22 − I ⎥⎦ ⎢⎣φ n 2 ⎥⎦ 21 ⎣

YU Dehao (余德浩) et al：Boundary Integral Equations and A Posteriori Error Estimates

Using Eq. (3),

37

Using Eq. (12), the residual is ⎡φ n1 ⎤ ⎡ S11 ⎢ ⎥=⎢ ⎣⎢φ n 2 ⎦⎥ ⎣ S21

S12 ⎤ ⎡φ 01 ⎤ ⎢ ⎥ S22 ⎥⎦ ⎣⎢φ 02 ⎦⎥

h h r (2) h = φ n − ( Lφ 0 + M φ n ) =

(7)

(φ n − Lφ 0 − M φ n ) + eh − Me h =

eh − Meh

φ 02 is obtained by solving S22φ 02 = φ n 2 − S21φ 01

(8)

i.e., rh(2) = ( I − M )eh

and then using

φ n1 = S11φ 01+ S12φ 02

(13)

(9)

to get the unknown φ n1 . This method has several advantages. First, the scale of the equation has been reduced since Eq. (8) can be solved alone and the operator S22 is symmetric posi-

(14)

This hyper-singular residual is computable and leads to the following result. Theorem 1 There are two positive constants C2 ≥ C1 > 0 such that C1 rh(2)

0

≤ eh

0

≤ C2 rh(2)

0

(15)

tive definite, so it is easy to be solved. Secondly, by using the solution φ 02 of Eq. (8), only the integral

where the error, eh, is defined by Eq. (10) and residual, rh(2) , is defined by Eq.(13).

must be evaluated to obtain φ n1 , which is much easier

2.3

than solving the equation. The calculation can be easily done in parallel.

Using Eq. (3), the residual is defined by

Residual using Eq. (3)

rh(3) = φ hn − Sφ 0 = (φ n − Sφ 0 ) + eh = eh

2 Residuals for the Dirichlet Problem Any potential problem has several equivalent integral equations[10,11] which can be used to define different residuals from the computed boundary element solution[12]. If the boundary integral in Eq. (1) is used for the initial solution, the error in the approximate solution φ nh is

(16)

This hyper-singular residual is just the error in the primary approximation. Equation (16) shows that eh = rh(3) for both global and local norms. Therefore, very accurate integrals of Sφ 0 give the best a posteriori error estimates as the residual r (3) h .

The residuals can then be calculated using Eqs. (1), (2), or (3).

The hyper-singular integrals in Eqs. (13) and (16) are both difficult to calculate. When the expression of S is known, Eq. (16) is much easier and more effective to apply than Eq. (13). In fact, the natural boundary integral method can be easily used to obtain the approximate solution of Eq. (3).

2.1

3

eh = φ hn − φ n

(10)

Residual using Eq. (1)

Using Eq. (1), the classical residual is defined as h r (1) h = φ 0 − ( K φ 0 + V φ n ) = (φ 0 − K φ 0 − V φ n ) − Veh = −Veh

(11)

Residuals for the Neumann Problem

Since V is a weakly singular integral operator, it smoothes the error eh , the properties of the residual

For the Neumann problems, φ n = g is known and φ 0 is unknown. If the solution is smooth enough and the approximate solution is given by φ h0 , the error is de-

rh(1) cannot accurately reflect the error properties, so it

noted by eh = φ 0h − φ 0

is not a proper measure of the error, eh . 2.2

which also has three kinds of residuals.

Residual using Eq. (2)

3.1

With Eq. (2), the iterative sequence is defined by

φ hn ' = Lφ 0 + M φ hn

(17)

(12)

Residual using Eq. (1)

Using the boundary integral Eq. (1), the approximate

38

Tsinghua Science and Technology, February 2005, 10(1): 35–42

C1 rh(3)

solution is

φ 0h

'=

K φ 0h

+ Vφn

(18)

≤ C2 rh(3)

−1/ 2

(25)

φ 0 + eh − K (φ 0 + eh ) − V φ n =

required precision, rh(3) is a good a posteriori error (19)

The operator K is compact (for C1 boundary). The error estimate is then: Theorem 2 There are two constants C1, C2, C2 ≥ C1>0, such that

C1 rh(1)

0

≤ eh

0

≤

C2 rh(1)

0

(20)

where the error is defined by Eqs. (17) and (19). Using Eq. (18), the relation between the errors is

K (φ 0 + eh ) − V φ h − φ 0 = Keh

(21)

The compact operator K filters out all the oscillatory components; therefore, the iteration result is superconvergent. Residual using Eq. (2)

Using Eq. (2), the hyper-singular residual is defined as

rh(2) = φn − ( Lφ 0h + M φ n ) =

φ n − L(φ 0 + eh ) − M φ n = − Leh

estimate.

4

(22)

The pseudo-differential operator L is of order +1. L tends to amplify the high frequencies of eh . Thus,

enough and that the error is given by h eh = [e1h , e2h ]T = Wh − W = [φ nh1 − φ n1 , φ 02 − φ 02 ]T

C1 rh(2)

−1/ 2

≤ eh

1/ 2

≤ C2 rh(2)

−1/ 2

(23)

where the error is defined by Eq. (17) and the residual is defined by Eq. (22).

Residual using Eq. (1)

Using Eq. (1), the residual is given by ⎡ r (1) ⎤ ⎡ V11 K12 ⎤ ⎡e1h ⎤ h rh(1) = ⎢ 1(1) ⎥=⎢ ⎥ ⎢ h ⎥ ≡ Ae h ⎣⎢ r2 h ⎦⎥ ⎣ −V21 I − K 22 ⎦ ⎣⎢e2 ⎦⎥

a posteriori error estimate. 4.2

Residual using Eq. (2)

Using Eq. (2), the residual is ⎡ r (2) ⎤ ⎡ I − M11 h rh(2) = ⎢ 1(2) ⎥=⎢ ⎣ r2 h ⎦ ⎣ − M 21

− L12 ⎤ ⎡ e1h ⎤ ⎢ ⎥ ≡ Be h − L22 ⎦⎥ ⎣e2h ⎦

(28)

which requires the computation of hyper-singular integrals. ⎡0 1⎤ The operator B is of order ⎢ ⎥ . If we assume ⎣0 1⎦ that the inverse of B exists, then C1 rh(2)

H −1/ 2 × H −1/ 2

≤ eh

C2 rh(2)

into Eq. (3), the hyper-singular re-

(27)

⎡ −1 0 ⎤ where the operator A is of order ⎢ ⎥ . Since this ⎣ −1 0 ⎦ operator smoothes the error eh , Eq. (27) is not a good

Residual using Eq. (3)

φ nh

(26)

Then there are three kinds of residuals as follows.

rh(2) should be a good a posteriori estimate of the error. Theorem 3 There exist two constants C1, C2, C2 ≥ C1>0, such that

Residuals for Mixed Problems

Assume that the solution W = [φ n1 , φ 02 ]T is smooth

4.1

eh ' = φ 0h '− φ 0 = K φ 0h + V φ n − φ 0 =

Substituting

1/ 2

rh(1) = φ 0h − ( K φ 0h + V φ n ) =

φ 0 − ( K φ 0 + V φ n ) + eh − Keh = ( I − K )eh

3.3

≤ eh

where the error is defined by Eqs. (17) and the residual is defined by Eq. (24). When the computed integral Sφ 0 reaches the

The residual is again defined as in Eq. (11) as

3.2

−1/ 2

H −1/ 2 × H 1/ 2 / P0

≤

H −1/ 2 × H −1/ 2

(29)

where C1, C2 are constants, C2 ≥ C1>0.

sidual is defined as rh(3) = φ n − Sφ 0h = (φ n − Sφ 0 ) − Seh = − Seh

(24)

The hyper-singular operator S is of order +1. Theorem 4 There exist two constants C1, C2, C2 ≥ C1>0, such that

4.3

Residual using Eq. (3)

Substituting Wh into Eq. (3), the residual is defined as

⎡ r (3) ⎤ ⎡φ nh1 ⎤ ⎡ S11 h rh(3) = ⎢ 1(3) ⎥=⎢ ⎥−⎢ ⎣ r2 h ⎦ ⎣⎢φ n 2 ⎦⎥ ⎣ S21

S12 ⎤ ⎡ φ1 ⎤ ⎢ h⎥ = S22 ⎥⎦ ⎣φ 2 ⎦

YU Dehao (余德浩) et al：Boundary Integral Equations and A Posteriori Error Estimates

⎡ I − S12 ⎤ ⎢ 0 − S ⎥ ≡ Deh 22 ⎦ ⎣

e

(30)

elements be defined on a family of K-meshes. Let u ∈ H k ( I 2 ) ∩ H l (Γ ), max(α , 0) ≤ 1 ≤ k , I1 ⊂ I 20 ⊂ Γ, and then for 0 < h ≤ h0 , S h ⊂ H k (Γ ), the local error of the Galerkin boundary element method satisfies the following estimates: e

≤ C{ r

0, I1

+ h k +1+ l u l + h k r

−α , I 2

−α

}

(31)

for −2 ≤ α ≤ 0, and

e

α

0, I1

≤ C{h r

+h

0, I 2

min( k +α ,t )

e e

e

0, I1 0, I1

≤ C{ r

1, I 2

≤ C{ r e

0, I1

0, I 2

+ h k +1 u

(32)

+ h k +1 u

≤ C{h r

h k +1 r

0, I 2

0, Γ

+ hk r

0, Γ 0, Γ

+h

+ hk r

k +3 / 2

1, Γ

}, α = 0 ;

u 1/ 2, Γ +

α =1

},

}, α = −1 ;

0, Γ

(33)

Note that V is a pseudo-differential operator of order −1; I −M, I−K, and I are all pseudo-differential operators of order 0, and L and S are pseudodifferential operators of order +1. The following results can be obtained using Theorem 5: For the Dirichlet problem, eh = φ nh − φ n , e

0, I1

≤ C r (1)

1, I 2

+ O(h k ),

e

0, I1

0, I1

0, I1

≤ C r (1)

(34)

0, I 2

≤ Ch r (2) ≤ Ch r (3)

+ O (h k ),

0, I 2

+ O(h k +1 ),

0, I 2

0, I 2

+ O (h k +1 )

(35)

The main terms on the right-hand sides of Eqs. (34) and (35) can be used as local error estimators.

6

Numerical Examples

Numerical examples will illustrate that the hypersingular residuals, calculated from the natural boundary integral equation, track the local error more effectively than the other error estimates. Therefore, the hyper-singular residual can be easily used as the basis for adaptive mesh refinement. The hyper-singular integrals can be computed using the method due to Yu[8,13]. The circular boundary [0,2π] is partitioned uniformly into N elements. The basis function is the “hat function” Li defined in Ref.[8]. The approximations for N

∑φ L and i i

i =1

for 0 ≤ α ≤ 2 . For the special cases,

e

= r (3)

the potential and flux are φ =

ul+

h min( k +α ,t ) r 0 }

0, I1

+ O(h k ),

0, I 2

For the Neumann problem, eh = φ 0h − φ 0 ,

A Posteriori Local Error Estimates

Suppose that the relationship between the residual, r, and the error, e, can be formulated as r = Ae. The following theorem for a posteriori local error estimates was given by Wendland and Yu[7]. Theorem 5 Assume A is a strongly elliptic pseudodifferential operator of order α with | α |≤ 2. Let | α | + min(0,α / 2) ≤ k ≤ t − 1 and the boundary

≤ C r (2) e

⎡0 1⎤ The operator D is of order ⎢ ⎥ which is much ⎣0 1⎦ simpler than B. Therefore, the residual in Eq. (30) is valuable as an a posteriori estimate for the mixed boundary value problem.

5

0, I1

39

φn =

N

∑φ

ni Li

.

i =1

Example 1 Consider a circular domain and the har1 monic function φ = r 2 sin 2θ , where (r ,θ ) denote 5 ∂φ 2 = r sin 2θ . Asthe polar coordinates. The flux is ∂n 5 suming r =5, then φ0 = 5sin 2θ and φn = 2sin 2θ . Con-

forming piecewise linear elements were used. All potentials, fluxes, and errors were normalized with respect to their maximums. For the Dirichlet problem described in Example 1, the integrals Lφ0 and Sφ0 were computed using the method due to Yu[8,9,13]. Figure 1 shows the exact and computed solutions. Figures 2 and 3 compare the error and residuals for 8 nodes and 16 nodes. Figure 4 compares the errors and residuals on a log-scale for 32 nodes. The hypersingular residuals accurately track the primary errors as shown in Figs.1-3.

40

Tsinghua Science and Technology, February 2005, 10(1): 35–42

For the Neumann problem,

φ n is given and satis-

fies the consistency condition: ∫ φn ds = 0 . The inteΓ

grals Lφ and Sφ h 0

h 0

were computed using the method

[8,9,13]

. due to Yu Figure 5 compares the computed and exact solutions. Figures 6-8 show how the errors were tracked by the different residuals. Fig. 1 Exact and computed solutions for the Dirichlet problem

Fig. 5 Computed and exact solutions for the Neumann problem Fig. 2 Errors and residuals for the Dirichlet problem using 8 nodes

Fig. 3 Errors and residuals for the Dirichlet problem with 16 nodes

Fig. 4 Errors and residuals for the Dirichlet problem with 32 nodes

Fig. 6 Errors and residuals for the Neumann problem using 8 nodes

Fig. 7 Errors and residuals for the Neumann problem using 16 nodes

YU Dehao (余德浩) et al：Boundary Integral Equations and A Posteriori Error Estimates

41

Eq. (1) with 32 nodes. Figure 10 compares the exact flux and the computed flux from the natural boundary integral iteration. Figure 11 compares the residuals from the error and the three different schemes.

Fig. 8 Errors and residuals for the Neumann problem using 32 nodes

Example 2 Consider a potential problem with a physical singularity on the boundary. Specifically, consider the Neumann problem whose solution with the singularity is

φn =

⎛θ π ⎞ sin ⎜ + ⎟ − cosθ ⎝ 2 4⎠

Fig. 10 Exact and computed flux for the Neumann problem with a singularity

， 0 ≤ θ < 2π .

1

⎡ ⎤2 ⎛θ π ⎞ 2 ⎢ 4 sin ⎜ + ⎟ − 2cosθ ⎥ ⎝ 2 4⎠ ⎣ ⎦

The exact solution is 1

⎧ ⎛θ π ⎞ 1 ⎫2 φ0 = ⎨ sin ⎜ + ⎟ − cosθ ⎬ , on Γ . ⎩ ⎝ 2 4⎠ 2 ⎭ Because e ∉ L2 (Γ ) , the norm e

L2 ( Γ j 0 )

cannot be

calculated, where Γ j 0 is an interval containing the singularity, (e, N j ) , j=1, …, N, is calculated instead of e

Fig. 11 Errors and residuals for the Neumann problem with a singularity

(2N )

1 − ,Γ j 2

, where N j = L2 j −1 , j=1, …, N, are the

piecewise linear functions for the mesh with 2N nodes. Figure 9 compares the exact solution and the computed solution using the classical boundary integral

Fig. 9 Exact and computed solutions for the Neumann problem with a singularity

The computational results for adaptive meshes, starting with 8 nodes or 16 nodes, and using |(r, Nj)| as error indicators, are compared with uniform meshes in Fig. 12. In the adaptive mesh, elements with error

Fig. 12 Comparison of results using adaptive and uniform meshes

42

Tsinghua Science and Technology, February 2005, 10(1): 35–42

indicators higher than some threshold value were refined to reduce the errors. The results in Fig. 12 show that the adaptive refinement leads to higher convergence rate than the uniform grids. Many other numerical examples have shown that the hyper-singular residuals are good a posteriori error estimates with many potential applications in adaptive boundary element methods.

[5]

Yu D H. A posteriori error estimates and adaptive approaches for some boundary element methods. In: Brebbia C A, Wendland W L, Kuhn G, eds. Boundary Elements. Springer-Verlag, 1987: 241-256.

[6]

Wendland W L, Yu D H. Adaptive boundary element methods for strongly elliptic integral equations. Numer. Math., 1988, 53: 539-558.

[7]

Wendland W L, Yu D H. A posteriori local error estimates of boundary element methods with some pseudo-

7

Conclusions

Residuals based on different boundary integral equations can be used as a posteriori error estimates in the boundary element methods. This paper describes estimates for Dirichlet, Neumann, and mixed boundary value problems. Theoretical analyses and numerical examples illustrate the important features of these error estimates. The hyper-singular residuals obtained from the natural boundary integral equation provide the best a posteriori error estimates, which can be used efficiently to improve adaptive boundary element calculations. References

differential equations on closed curves. J. Comput. Math., 1992, 10(3): 273-289. [8]

Yu Dehao. The Mathematical Theory of the Natural Boundary

Method.

Beijing:

Science

Press,

1993.

(in Chinese) [9]

Yu D H. Natural Boundary Integral Method and Its Applications. Dordrecht: Kluwer Academic Publishers, 2002.

[10] Wendland W L. Strongly elliptic boundary integral equations. In: Iserles A, Powell M J D, eds. The State of the Art in Numerical Analysis. London: Oxford University Press, 1987, 9: 511-652. [11] Zhu Jialin. Analysis of the Boundary Element Method on Elliptical Boundary Problems. Beijing: Science Press, 1991. (in Chinese) [12] Menon G, Paulino G H, Mrkherjee S. Analysis of

[1] [2]

Babuska I, Strouboulis T. The Finite Element Method and

hypersingular residual error estimates in boundary element

Its Reliability. Oxford: Clarendon Press, 2001.

methods for potential problems. Comput. Methods Appl.

Babuska I, Yu D H. Asymptotically exact a posteriori error estimator for bi-quadratic elements. Finite Elements

[3]

in Analysis and Design, 1987, 3: 341-354.

canonical integral equations in interior or exterior circular

Yu D H. Asymptotically exact a-posteriori error estimator

domains. J. Comput. Math., 1983, 1(1): 52-62.

for elements of bi-even degree. Math. Numer. Sinica, 1991, 13(1): 89-101. [4]

Mech. Engrg., 1999, 173: 449-473. [13] Yu D H. Numerical solutions of harmonic and biharmonic

Yu D H. Asymptotically exact a-posteriori error estimator for elements of bi-odd degree. Math. Numer. Sinica, 1991, 13(3): 307-314.