A POSTERIORI ERROR ESTIMATES FOR hp {BOUNDARY ELEMENT ...

A POSTERIORI ERROR ESTIMATES FOR hp { BOUNDARY ELEMENT METHODS Carsten Carstensen Mathematisches Seminar II, Christian-Albrechts-Universitat Kiel Ludewig-Meyn-Str. 4, 24098 Kiel, Germany Stefan A. Funken, Ernst P. Stephan Institut fur Angewandte Mathematik, Universitat Hannover Welfengarten 1, 30167 Hannover, Germany Abstract: This paper presents a posteriori error estimates for the hp{version of the boundary element method. We discuss two rst kind integral operator equations, namely Symm's integral equation and the integral equation with a hypersingular operator. The computable upper error bounds indicate an algorithm for the automatic hp{adaptive mesh{re nement. The eciency of this method is shown by numerical experiments yielding almost optimal convergence even in the presence of corner singularities. AMS(MOS): 65N35, 65R20, 65D07, 45L10. KEY WORDS: a posteriori error estimates, boundary element methods, adaptive boundary element methods, hp{version

1 Introduction In the numerical treatment of partial dierential equations hp{methods are established for the ecient solution of problems even if they are involving corner singularities. The exponential convergence is proved for the boundary element method (BEM) in [2] where we know a priori how to re ne the mesh: geometric grading with linear (decreasing) distribution of polynomial degrees towards corner singularities. In general, for example treating singular or nearly singular data or in case of a smooth boundary with high curvature, the precise information is missing or is, e.g., in case of elasticity operators with nonhomogeneous coecients, only obtained very dicultly and costly. In these situations, a self{adaptive mesh{re nement is an important tool for the ecient numerical solution. Although hp{methods and adaptive methods are well established for the nite element method, comparably little is known for the BEM: cf. [9, 12, 15, 18] for hp{methods and [7, 8, 10, 14, 15, 16, 20, 21] for h{adaptive methods. Clearly, to gain exponential

2

C. CARSTENSEN; S. A. FUNKEN; E. P. STEPHAN

convergence we need to use the hp{method adaptively; numerical experiments can be found in [14, 15, 16] lacking a rigorous justi cation. In this paper we present a{posteriori estimates for the hp{version of the BEM. The computable error bound is important in itself for the reliability of a numerical computation. Here we use the upper bound to steer the mesh{re nement even if the number of degrees of freedom is small. For coarse meshes, a priori information on the asymptotic behavior are not applicable to balance a global re nement for part of the domains with smooth solutions and a local re nement towards singularities. But for an ecient treatment you have to do this with only a few degrees of freedom. We address this questions considering two prototypes of rst kind integral equations: Symm's integral equation (with a pseudo{dierential operator of order ?1) related to the Dirichlet problem for the Laplacian, and the hypersingular equation (with the normal derivative of the double layer potential; a pseudo{dierential operator of order +1), equivalently related to the Neumann problem for the Laplacian, in x3 and x4. Observation shows that the abstract framework for the derivation of a posteriori estimates, given in [7] and recalled in x2 for convenient reading, is feasible for hp{methods as well. The resulting a posteriori estimates depend on the sum of the local residuals where the latter are weighted with the factor hpjj ; hj denotes the length of the element ?j and pj the polynomial degree of the BEM Galerkin solution on ?j . Based on these estimates we formulate an hp{adaptive algorithm in x5. The implementation of the hp{adaptive BEM is described in x6, numerical experiments are reported in x7. The results clearly indicate the hp{adaptive BEM converges exponentially in a model case where this is already proved for the BEM on geometrically graded meshes [2, 12].

2 Abstract framework of a posteriori error bounds for boundary integral equations In this section, we brie y describe the abstract frame for for BEM established in [7]. The setting of the notation becomes clearer in the applications in x3 and x5. Let X and Y1 Y0 be real Banach spaces and let Y := [Y0 ; Y1] be de ned by interpolation, 0 1, cf. [1]. Let A : X ! Y be linear, bounded, surjective and injective. Then, x a right hand side f 2 Y1 and the solution u 2 X of Au = f: (1) We apply the Galerkin method to approximate u. Let Sh X and Th Y0 be nite dimensional subspaces such that there exists some uh 2 S h satisfying

th (Auh ) = th (f ) for all th 2 Th:

(2)

The residual Rh := f ? Auh 2 Y1 and the error eh := u ? uh are related by the following estimate.

A POSTERIORI ERROR ESTIMATES

Theorem 1 [7] and there holds

3

There exists 2 Y0 satisfying

kk2Y0 = kRhk2Y0 = (Rh)

(3)

kehkX C kRhkY1 : t inf k ? th k1Y?0 : 2T

(4)

h

h

The constant C := c kA?1 kL(Y ;X ) depends on A and the interpolation constant c . Proof: To be self contained we give a brief outline of the proof. The existence of 2 Y0 satisfying (3) follows from a well{known corollary of Hahn{Banach's theorem. Since A : X ! Y is bijective, eh = A?1 Rh and

keh kX kA?1kL(Y;X ) kRhkY : By interpolation, kRhkY c k Rh k1Y?0 k Rh kY1 [1]. The proof is concluded by standard duality arguments using (2) and (3)

k Rh k2Y0 = (Rh) = ( ? th)(Rh) k Rh kY0 k ? th kY0 (th 2 Th): 2 The crucial point is that can explicitly be determined in applications.

3 A posteriori error bounds for hp{BEM for a weakly singular equation Given a bounded Lipschitz domain IRn with boundary ? the Dirichlet problem for the Laplacian is equivalently related with the so{called Symm's integral equation for the unknown density ( is the normal derivative of the displacement) on ?

V (x) = g(x) (x 2 ?)

(5)

with the weakly singular operator V of the single layer potential,

V (x) :=

Z

?

(y)G(x; y)dsy:

Here G(x; y ) := ? 1 log jx ? y j for n = 2, G(x; y ) := cn jx ? y j?n for n 3 with some constant cn > 0. In (5) g is some right hand side determined by the given Dirichlet data uj? . Let H s (?) := fuj? : u 2 H s+1=2(IRn )g, 0 s 1, be the usual Sobolev space on ?, H 0(?) = L2 (?) and H ?s (?) is the dual space of H s (?), with respect to the duality which is de ned for smooth functions u; v by

< u; v >=

Z

?

u v ds:

4


V is a Fredholm operator of index zero, it is injective in L1 (?) for n 3 or for n = 2 if cap(?) 6= 1 [4, 17]. Thus V is bijective from H s?1 (?) onto H s(?) provided that n 3 or n = 2 and cap(?) = 6 1 which is assumed throughout this paper (cap(?) denotes the capacity or conformal radius or trans nite diameter of ?). Let Thp := Shp (?) L2 (?), (p 0) denote the vector space of piecewise polynomials of degree p with respect to a "triangulation" of the boundary ? = [Nj=1 ?j where the "elements" ?1 ; : : :; ?N are either equal or have at most one common point or side, respectively. We assume ?j is an interval for n=2 or a triangle or parallelogram for n=3. If n = 3, we assume the interior angles being uniformly bounded below, i.e., there exists a global mesh{independent constant c > 0 such that meas(?j ) ch2j . We de ne h 2 Sh0(?) as the local mesh size, h(x) := hj := diam(?j ) if x 2 ?j , diam(?j ) > 0 is the element size. The polynomial degrees are described by p 2 Sh0 (?) the (integer valued) degree function de ned such that each hp 2 Shp (?) is a polynomial on ?j of degree p(x) = pj , x 2 ?j . Let ph 2 Thp denote the Galerkin solution of (5), i.e. 0 =< tph ; R >

for all tph 2 Shp(?)

(6)

where R := g ? V ph is the residual. Assume g 2 L2 (?) so that, according to Thp L2 (?), R = V ( ? ph ) 2 H 1(?) due to the mapping properties of the single layer potential [4]. Let r? denote the gradient with respect to ?. Theorem 2 Under the above assumptions, we have for 0 s 1

k ? ph kH ?s(?) C kRk1H?1s(?) k p +h 1 r?RksL2(?): The constant C > 0 depends not on h or p. The proof relies on Theorem 2 and the following approximation estimate. Lemma 1 [3, 18] On each element ?j and for each f 2 H k (?j ) there is a polynomial fj of degree pj such that

k f ? fj kH q(?j ) Ck

hj ?q

pjk?q + 1

k f kH k(?j )

for all q with 0 q k and := minfpj + 1; kg. The constant Ck > 0 depends only on k but not on h or p. Proof: The assertion is proved in [3, Lemma 4.5] for n = 3 and in [18] for n = 2; we remark that we enlarge the nominator by 1 changing the constant Ck to allow pj = 0 (in which case the assertion is well{known). 2


5

Proof of Theorem 3. Apply Theorem 2 with A = V , X = H ?s (?), Y0 = L2(?), Y1 = H 1(?), = 1 ? s, Y = H 1?s (?) so that (4) holds (with c = 1). Note that = R

satis es (3) and it remains to verify

infp kR ? tph kL2 (?) C k p +h 1 r? RkL2(?) :

t 2Sh (?) p h

(7)

Since there is no continuity required across inter element boundaries it remains to prove for each element ?j that inf kR ? tph kL2 (?j ) C k p +h 1 r? RkL2(?j ) :

tph 2Shp (?)

(8)

From Lemma 3 we obtain with, k = 1 and q = 0,

infp kR ? tph kL2 (?j ) C p +h 1 k R kH 1 (?j ) :

t 2Sh (?) p h

(9)

It is known that we can replace the H 1{norm in the right{hand side of (9) by its seminorm which then proves (8). To see this, we may use Lemma 3 and obtain for some constant polynomial cj (pj = 0) (de ned as the integral mean of R on ?j )

k R ? cj k2L2(?j ) (C hj )2 k R ? cj k2H 1(?j ) and so, for h small enough,

k R ? cj kL2(?j ) 2C hj k r?R kL2(?j ): Finally, we use Lemma 3 again for R ? cj and obtain from this (owing to pj 0) inf kR ? tph kL2 (?j ) p infp kR ? cj ? tph kL2 (?j ) th 2Sh (?) tph 2Shp (?) C p h+j 1 k R ? cj kH 1(?j ) j 2C k p +h 1 r?R kL2(?j ): This indicates (9) and veri es (8). Therefore, Theorem 2 concludes the proof.

2

4 A posteriori error bounds for hp{BEM for a hypersingular equation Given a two dimensional bounded domain with polygonal boundary ? the Neumann problem for the Laplacian is equivalently related to the hypersingular integral equation,

Wv(x) = f (x) (x 2 ?)

(10)

6


for the unknown displacement v (= uj? ) on ? with the hypersingular operator Z Wv(x) := ? 1 @ v(y) @ log jx ? yjds :

@nx ? @ny Let H0s (?) := f 2 H s (?) :< 1; >= 0g. It is known that W : H0s(?) ! H0s?1 (?)

y

is linear, bounded and bijective [7]. Let Thp := S ph (?) H01(?) be the vector space of continuous and piecewise polynomials of degree p, p 1 (with respect to the mesh as in the previous section) with integral mean zero. Given f 2 H00(?) let vhp 2 Thp denote the Galerkin solution, i.e. for all tph 2 Thp ; (11) 0 =< tph ; R > where R := f ? Wvhp 2 H00 (?). Theorem 3 Under the above assumptions, we have for 0 s 1 ()

s : kv ? vhpkH s(?) C kRksL2(?) k hp Rk1L?2(?)

Proof: Apply Theorem 2 with A = W , X = H0s (?), Y0 = H0?1 (?), Y1 = H00(?), = s, Y = H0s?1 (?). Note that, by construction of the Sobolev spaces H s(?), the duality in (11) is the duality between L2(?) and L2 (?) as well as between H01(?) and H0?1 (?). Hence (2) is satis ed. Note that integrating ?R twice with respect to the arc length gives , = ?I I (R) satis es (3) where I denotes the operator of integration. Therefore inf p k ? tph kH01 (?) = pinf p kI (R) ? tph 0 kL2 (?) k hp RkL2 (?)

t 2Th th 2Th tph 0 as the L2{projection of I (R) and p h

by choosing using standard interpolation/approximation arguments as in Lemma 3. Therefore, Theorem 2 yields the assertion. 2

5 Adaptive feedback algorithm For a given "triangulation" of the boundary ? = [Nj=1 ?j and a corresponding p{ distribution we can compute an approximation of the contribution aj of one element ?j by numerical integration. For Symm's integral equation aj := kR0kL2 (?j ) and for the integral equation with the hypersingular operator aj := kRkL2(?j ) . Then, with some constant c > 0, we have the error in the energy norm (s = 1=2) is bounded above by

c

N

X

j =1

a2j

=4

1

N

h2j 21=4 2 aj j =1 pj

X

(12)


7

where hj := j?j j is the length of the element ?j and pj is the polynomial degree on ?j . This a posteriori error estimate is not very useful for absolute error control unless the constant c > 0 (or at least an upper bound) is known. But it can be used to compare the contributions to the local error and to control the relative error. Hence, we may steer the process of mesh re ning such q that (12) becomes small and the PN 2 terms in the sums are uniformly distributed. Let a := j =1 aj , then the contribution 2 of any element ?j to the error is c4 hp2jj a2j a2 . Thus we are lead to consider hpjj aj and to organize the mesh re nement so that all terms hp11 a1 , . . . , hpNN aN will be of the same magnitude. The mesh degree distribution in our numerical examples are steered by the following algorithm. Algorithm (A) Given some coarse e.g. uniform mesh, re ne it successively by halving some of the elements or increasing the polynomial degree due to the following rule. For any triangulation de ne a1 ; : : :; aN as above and "re ne" the element (?j ; pj ) if and only if P k a2 hj a2 Nk=1 pkh+1 k : (13) j p +1 N j

After the decision to "re ne element ?j " according to (13) we have to decide whether to perform an h{ or a p{re nement, i.e. whether to halve the element length j?j j or to increase the polynomial degree pj by one. This is decided as follows. First, one increases pj by one for each element ?j which might be re ned and computes a new Galerkin solution ph . Now one compares the old residual R(?j ; pj ) with the new residual R(?j ; pj + 1). Thus, if R(?j ; pj + 1) R(?j ; pj ) one performs a p{re nement whereas, if R(?j ; pj + 1) > R(?j ; pj ) one performs an h{re nement, where is given.

6 Implementation In this section we explain the implementation of the method used in the numerical experiments. On a standard element (?1; 1) we de ne as trial functions antiderivatives of Legendre q R 1 ? 1+ 2 j ? 3 polynomials Pi , namely N1( ) := 2 , N2( ) := 2 and Nj ( ) := 2 ?1 Pj ?2 (t)dt with j = 3; : : :; p + 1. We consider the single layer potential V kj on an element ?j with end points ; 2 IR2. With the linear mapping ( [; ] IR2 Q : [?1t; 1] ! 7! 12 [( ? )t + ( + )]

8


we obtain the integral

lj () log j ? jds = j?2j j

Z

?j

1

Z

?1

p

xl log ax2 + bx + cdx;

where a; b; c depend on and ?j . With the two conjugate complex roots z0; z0 of x2 + b c a x + a we get j?j j Z 1 xl log pax2 + bx + cdx 2 ?1 Z 1 Z 1 j ? = 4j j xl log jajdx + j?2j j < xl log(x ? z0 )dx : ?1

?1

For the last integral we have the relation Z

xl log(x ? z)dx = l +1 1 xl+1 ? z l+1 log(x ? z) ? l +1 1 h

i

l

xl?k+2 z k?1 k=1 l ? k + 2 +1 X

(14)

[13, formula 2.729.1]. Hence, Z 1 xldx V lj (x) = ? j?4j j log 41 [( 1 ? 1 )2 + ( 2 ? 2)2 ] ?1 Z 1 ? j?2j j < xl log(x ? z0)dx

?1

with

z0 = c3 + ic4 c1 = 12 x1 ? ( 1 + 1); c2 = 12 x2 ? ( 2 + 2); ( 1 ? 1 )c2 + ( 2 ? 2 )c1 ( 1 ? 1 )c1 + ( 2 ? 2 )c2 c3 = ( ? ) + ( ? ) c4 = ( ? ) + ( ? ) : 1 1 2 2 1 1 2 2 If we neglect the appropriate mappings Q and Q , it is obvious that the main problem is to solve

Z

1

Z

1

xl log(x ? z0 )dxdy (15) where we have to pay attention to the fact, that z0 depends on y . Next we deal with the

where ki , lj are monomials which occur in the Galerkin matrix. The term < W ; > satis es the relation < Wu; v >=< V dsd u; dsd v > and since the single layer potential is linear, we can restrict ourselves to describe the calculation of < V ; > for monomials kj (degree k, supp(kj ) = ?j ). In the case of the hypersingular equation p p we use instead of vh 2 Sh the additional restriction < uph ; 1 >= 0. The Galerkin schemes for (10) is now: Find (uph ; a) 2 Shp (?) IR with < Wuph ; vhp > + < a; vhp > = < f; vhp > < uph ; 1 > = 0 for all vhp 2 Shp(?) [5]. Here the second term of the rst equation and the second equation can be calculated exactly. Given a harmonic function u, then functions f and g on the right hand side of the integral equations (5) and (10) are calculated by the relations (I ? K 0) @u @n = f and (I + K )u = g , @u j p p ? respectively. The expressions < f; vh >=< @n ; (I ? K )vh > and < g; ph >=< uj? ; (I + K 0)ph > are calculated adaptively in the way, that the inner expressions (I ? K )vhp and (I + K 0)ph are computed exactly and the outer integrals are calculated adaptively: We subdivide the interval ?j until the accuracy (we used a parameter of 10?10 ) for each subinterval is reached by the used quadrature rule (4-point Gauss-Legendre quadrature rule). The decision whether or not to subdivide the interval ?j is made by adopting an idea of Gander [11]. The residual corresponding to Symm's equation := V ph ? g is approximated by a 16 point Gaussian quadrature formula. ph and g are known exactly and since ph is a piecewise polynomial function V ph can be calculated exactly. We compute k0 kL2 (?k ) with the same Gaussian quadrature as above. The L2(?)-norm of ()

Wuph ? f

is approximated by a 16 point Gaussian quadrature formula.p Here, uph and f are known exactly and since uph is a piecewise ppolynomial function V dudsh can be calculated exactly. p du du p Using the relation Wuh = ? dsd V dsh , the remaining dierentiation of V dsh is computed with replacing dsd by a symmetric dierence quotient with step size 10?5. The H ?1=2(?)-norm ( H01=2(?)-norm) is equivalent to the \energy norm" p p kkV = < V ; > (kukW = < Wu; u >) which is used in the sequel. For x 2 ?j and = ? hp we compute M Z X 1 V (x) = ? (y) log jx ? yjds

k=1

?k

y


11

by numerical quadrature rules. For j 6= k we apply a 16 point Gaussian quadrature formula. For j = k we divide and transform the integral such that the \singular point" x lies at the end of the unit interval. Then we apply an 8 point Gaussian quadrature rule with logarithmic weights [19]. This explains the approximation of the \energy norm" we use.

7 Numerical experiments We consider two numerical examples, namely the hp{version of the Galerkin method for the weakly singular integral equation and the hypersingular integral equation. In both examples we choose the parameter = 0:5 in the adaptive algorithm to steer the p{re nement.

7.1 Example 1: Symm's integral equation

For ? we take the L{shaped polygon with vertices (0; 0), (1; 0), (1; 1), (?1; 1), (?1; ?1), (0; ?1). We consider (5) with g = (I + K )f , where f = r2=3 sin( 32 )j? with polar coordinates (r; ) concentrated at the origin. In Fig. 1 we present the numerical results for the above described adaptive hp{version. For comparison we also show the error curve for the Galerkin solution obtained via both hp{version and various geometric meshes (with mesh parameters) and the h{ and p{version.

Figure 1, 2: The error of the Galerkin solutions in Example 1. Figure 3{5: Adapted meshes for Symm's equation

7.2 Example 2: A hypersingular integral equation

On the polygon ? with vertices (0; 0), (1; 1), (?1; 1), (?1; ?1), (1; ?1) we consider the hypersingular integral equation (10) with the solution v = r2=3 sin( 32 ) in polar coordinates. Again, the convergence rates for uniform meshes are non{optimal due to the singular behavior of the solution.

Figure 6, 7: The error of the Galerkin solutions in Example 2. Figure 8{10: Adapted meshes for hypersingular integral equation

7.3 Conclusions

The conclusions in the two examples are analogous and can be described as follows. The hp{adaptive BEM re nes the mesh towards the singularity and increases the polynomial degree on elements not close to the singularity. Note carefully, that only the origin is a singular point but the other corner points (potentially yielding singularities as well) are treated correctly, i.e., no h{re nement but p{increase. So far the meshes indicate that

12


the hp{self adaptive scheme mimics a geometric mesh grading which is regarded as an optimal strategy. However, the moderate use of p{increase is perhaps surprising but is caused by the fact that the approximation error is very high near the singularity and comparably low in the remaining part of the boundary. So the main increase of degrees of freedom is caused by a local h{re nement to approximate the singularity. The plots in Figure 1 and 6 represent the relative errors in the energy norm as a function of degrees of freedom for various experiments connected by straight lines. Note the logarithmic scaling on both axis so that an algebraic convergence results in a straight line with a slope which is the experimental convergence rate. One observes convergence rates for the uniform h{ and p{method which are poor according to the singularity (see [18] for a proof of that). The hp{method used with various gradings gives a better convergence, the curve is concave not a straight line. To check if we have exponential convergence or not we consider Figure 2 and 7 where the error is plotted in a logarithmic scale while thep number n of degrees of freedom is increased with a square root, i.e., the x axis shows n. A straight line in this picture indicates exponentially fast convergence and this is observed for the hp{version with a geometrically graded mesh as well as for the adaptive hp{method introduced in x5. This proves numerically, that our adaptive hp{ BEM is of exponential convergence and hence a powerfull tool in the ecient treatment of integral equations of the rst kind.

Acknowledgement The work is partly supported by DFG research group at the University of Hannover.

References

[1] J. Bergh, J. Lofstrom: Interpolation spaces. Springer Berlin 1976. [2] I. Babuska, B.Q. Guo and E.P. Stephan, On the exponential convergence of the h{p version for boundary element Galerkin methods on polygons, Math. Meth. Appl. Sci. 12 (1990) 413{427. [3] I. Babuska, M. Suri: The h{p version of the nite element method with quasiuniform meshes. Mathematical Modelling and Numerical Analysis 21 (1987) 199|238. [4] M. Costabel: Boundary integral operators on Lipschitz domains: Elementary results. SIAM J. Math. Anal. 19 (1988) 613{626. [5] M. Costabel, E.P. Stephan: The normal derivative of the double layer potential on polygons and Galerkin approximation. Appl. Anal. 16 (1983) 205{228. [6] M. Costabel, E.P. Stephan: Boundary integral equations for mixed boundary value problems in polygonal domains and Galerkin approximation. Banach Center Publ. 15 (1985) 175{251. [7] C. Carstensen, E.P. Stephan: A posteriori error estimates for boundary element methods. Math. Comp. (1994) in press.


13

[8] C. Carstensen, E.P. Stephan: Adaptive boundary element methods for some rst kind integral equations. SIAM J. Numer. Anal. (1994) accepted for publication. [9] V. Ervin, N. Heuer, E.P. Stephan: On the h ? p version of the boundary element method for Symm's integral equation on polygons. Comput. Meth. Appl. Mech. Engin. 110 (1993) 25{38. [10] B. Faehrmann: Lokale a{posteriori{Fehlerschatzer bei der Diskretisierung von Randintegralgleichungen. PhD{thesis, University of Kiel, FRG (1993). [11] W. Gander: Computermathematik. Birkhauser, Basel (1985) [12] B. Guo, N. Heuer and E.P. Stephan, The h{p version of the boundary element method for transmission problems with piecewise analytic data, to appear. [13] I.S. Gradstein, I.M. Ryshik: Summen-, Produkt und Integral-Tafeln. Band 1, Verlag Harri Deutsch, Frankfurt/Main (1981) [14] N. Heuer: hp{Versionen der Randelementemethode. PhD{thesis, University of Hannover, FRG (1992). [15] F.V. Postell, E.P. Stephan: On the h?, p? and h ? p versions of the boundary element method { numerical results. Computer Meth. in Appl. Mechanics and Egin. 83 (1990) 69| 89. [16] E. Rank: Adaptive boundary element methods. in: C.A. Brebbia, W.L. Wendland and G. Kuhn, eds., Boundary Elements 9, Vol. 1, 259|273. Springer Verlag Heidelberg 1987. [17] I.H. Sloan, A. Spence: The Galerkin Method for Integral Equations of the rst kind with Logarithmic Kernel: Theory. IMA J. Numer. Anal. 8 (1988) 105|122. [18] E.P. Stephan, M. Suri: The hp{version of the boundary element method on polygonal domains with quasiuniform meshes. Mathematical Modelling and Numerical Analysis 25 (1991) 783{807 [19] H.H. Stroud, D. Secrest: Gaussian quadrature formulas. Prentice Hall, Englewood Cli 1966. [20] W.L. Wendland, De-hao Yu: Adaptive boundary element methods for strongly elliptic integral equations. Numer. Math. 53 (1988) 539|558. [21] W.L. Wendland and De-hao Yu: A posteriori local error estimates of boundary element methods with some pseudo{dierential equations on closed curves. Journal for Computational Mathematics 10 (1992) 273|289.

14


h-uniform p-unifom h-adaptiv hp-adaptiv hp-geometric, s=0.15 hp-geometric, s=0.20

Relative Error in Energy Norm

0.1

0.01

0.001

10

100 Number of Unknowns

1000

Figure 2: The error of the Galerkin solutions in Example 1. hp-adaptiv hp-geometric, s=0.15 hp-geometric, s=0.20


0.1

0.01

0.001

2

3

4

5 6 7 ( Number of Unknowns )^(1/2)

8

9

Figure 3: The error of the Galerkin solutions in Example 1.

10

A POSTERIORI ERROR ESTIMATES 0

0

15

0

0

0

0

0

0

0

0

1

0

0

1

0

0

0

0

0

0

0

1

0

0

0

1

1

1

1

1 1 0

0 1

1

0

0

1

0

1

0

1 11 1 1

1

1

1

1

1

1

1

1

1

0

0

Figure 4: Adapted meshes for Symm's equation

1

16

C. CARSTENSEN; S. A. FUNKEN; E. P. STEPHAN 1

0

1

1

1

1

1

1

1

1

1 1

1

1

1 1

1

1

1

1 1

1

1

1

0

1

1

0

1

1

1

1

1

1

1 1

1

1

1

1

0

1

1

0

1

1

1

1

1

1

1 1

1

1

1

1

0

1



17

0

1

1

1

1

1

1

1 1

1

1

1

1

1 1

1

1

1 1

1

1

1 1

1

1

0

1

1

0

0

1

1

1 0

0

1

1

1

1 1

1

1

1

1

1

1

0

1

1

2

1 0

1

1

1

1

1 1

1

1

1

1

2


18


h-uniform p-unifom h-adaptiv hp-adaptiv hp-geometric 0.2


0.1

0.01

0.001

0.0001 10

100

1000

10000

Number of Unknowns

Figure 7: The error of the Galerkin solutions in Example 2. hp-adaptiv hp-geometric 0.2


0.1

0.01

0.001

0.0001 2

3

4

5

6 7 8 ( Number of Unknowns )^(1/2)

9

10

Figure 8: The error of the Galerkin solutions in Example 2.

11


1

19

1

1

1 1

1

1

1

1 1

1

1 1

1

1

1

1

1

1

1

1

1 1

1 1

1

1 1

1

1 1

1

1 1

1

1

1 1

1

1

1

2

1

1

1

1 1

1 1

1

1

1

1

1

1

1 1

1

1 2

1 1

1

1

Figure 9: Adapted meshes for the hypersingular integral equation

20

C. CARSTENSEN; S. A. FUNKEN; E. P. STEPHAN 1

2

1 1

1 1

1 1

1

1

1 1

1

1

1 1

1 2

1 1

1

1

2

2

1 1

2 1

1 1

1

1

1 1

1

1

1 1

1 2

1 1

2

2

2

2 1 1

1

2 1

1 1

1 1 1 1

1 1

1 2

1 2

2

2



21

2 1 1

1

2 2

1 1

1

1

1 1

1

1

1 1

2 2

1 2

2

2

3

2 1 1

1

3 2

1 1

2

1

1 1

1

2

1 1

2 2

1 2

2

2

3

2 2 2

1

3 2

1 1

2

1

2 2

1

2

1 1

2 2

1 2

3

3
