See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/2691423

Numerical Integration for Multivariable Functions with Point Singularities

Article · November 1994

Source: CiteSeer

CITATIONS READS
8 59

2 authors, including:

Kendall E. Atkinson
University of Iowa
159 PUBLICATIONS 9,738 CITATIONS

SEE PROFILE

All content following this page was uploaded by Kendall E. Atkinson on 30 October 2013.

The user has requested enhancement of the downloaded file.

 Numerical Integration for Multivariable
Functions with Point Singularities
Yajun Yang  and Kendall E. Atkinson y
October 1992

Abstract
We consider the numerical integration of functions with point singularities over
a planar wedge S using isoparametric piecewise polynomial interpolation of the
function and the wedge. Such integrals often occur in solving boundary integral
equations using the collocation method. In order to obtain the same order of
convergence as is true with uniform meshes for smooth functions, we introduce an
adaptive re nement of the triangulation of S . Error analyses and several examples
are given for a certain type of adaptive re nement.

Key Words: numerical integration, quadratic interpolation, adaptive re nement.

AMS(MOS) Subject Classi cation: 65D30, 65D32

 Department of Mathematics, The University of Iowa, Iowa City, IA 52242

y Department of Mathematics, The University of Iowa, Iowa City, IA 52242

1
1 Introduction
We consider the problem of approximating an integral of the form
Z

f (x; y)dS
S

where S is a closed bounded connected set in R2. When the integrand has singularities
within the integration region, the use of a standard quadrature method may be very inef-
cient. There are several ways to deal with this problem. One approach is to use a change
of variables [10] to transform a singular integrand into a well behaved function over a new
region. A second approach is to use adaptive numerical integration, to place more node
points near where the integrand is badly behaved in order to improve the performance
of a standard quadrature method. A third approach is to use extrapolation methods to
construct new and more accurate integration formulas based on the asymptotic expan-
sion for the quadrature error in the original quadrature rule. The generalization of the
classical Euler-Maclaurin expansion to functions having a particular type of singularity,
as obtained by Lyness [8], [9], provides a basis for extrapolation methods in the same way
as the Euler-Maclaurin expansion is used as a basis for Romberg integration. For the
one dimensional case, a standard discussion of these methods and others can be found in
Atkinson [3].
In this paper we discuss adaptive numerical integration for the two dimensional case.
The process of placing node points with variable spacing so as to better re ect the inte-
grand is called adaptive re nement or grading the mesh. We propose a type of adaptive
re nement for which the order of convergence is the same as for smooth functions, but
with the integrand function having a particular type of singularity, speci ed below.
In the case of smooth integrand functions, Chien [6],[7] obtained that for integrals
over a piecewise smooth surface in R3, the numerical integration using isoparametric
piecewise quadratic interpolation for both the surface and the integrand leads to the
order of convergence O(4) . In this,  = max1KN (K ); K = maxp;q24K j p ? q j

2
and f4K j K = 1;


; N g is a quasi-uniform triangulation of the surface. In terms
of the number N of triangles, the order of convergence is O(1=N 2 ). We extend these
results to singular integrands.
To simplify our discussion, we study only the following class of problems. The inte-
gration region S is a wedge, i.e.
S = f(x; y) 2 R2 j 0  r  1; 0    g ; 0 <  < :
The integration function is singular at only the origin, and it is of the form

f (x; y) = r '()h(r)g(x; y); > ?2 (1:1)

or the form
f (x; y) = r ln r '()h(r)g(x; y); > ?2 (1:2)
Here (x; y) are the Cartesian coordinates with the corresponding polar coordinates
(r; ). The functions '(); h(r); g(x; y) are assumed to be suciently smooth. This is
also the class of functions considered in Lyness [8] .
In section 2, we describe the triangulation of the wedge S and the adaptive re ne-
ment scheme we use. The interpolation-based quadrature formulas are given in section 3.
Section 4 contains the error analyses and a discussion of the results. Numerical examples
are given in section 5.
This paper presents detailed results for only the use of quadratic interpolation. The
method being used generalizes to other degrees of piecewise polynomial interpolation,
and the results are consistent with the kind of results we have obtained for the quadratic
case. This generalization for other degrees of interpolation is given in section 6.
Because integrals of non-smooth functions often occur in solving boundary integral
equations using the numerical method, we expect that the error analysis of the present
numerical integration methods eventually will lead to better numerical methods for the
solution of boundary integral equations.

3
2 The triangulation and adaptive re nement
In this section, we describe the triangulation of the wedge S and discuss its re nement
to a ner mesh. Let
 n = f4K;n j 1  K  Nng (2:1)
be the triangulation mesh for a sequence n = 1; 2; : : :. In referring to the element 4K;n,
the reference to n will be omitted, but understood implicitly.
The initial triangulation  1 of S is obtained by connecting the midpoints of the sides
of S using straight lines. The sequence of triangulations  n of (2.1) will be obtained
by successive adaptive and uniform re nements based on the initial triangulation. We
construct this sequence as follows, and we call it an `La + u' re nement with L a
positive integer. Given f4K;n j 1  K  Nng , we divide the triangle containing the
origin into four new triangular elements. For the resulting triangulation, repeat the
preceding subdivision. After doing this L times, we divide simultaneously every triangle
into four new triangles. The nal triangulation is denoted by f4K;n+1 j 1  K  Nn+1 g.
In other words, at level n, we perform L times an adaptive subdivision on the triangle
containing the origin, and then we do one simultaneous subdivision of all triangles. An
advantage of this form of re nement is that each set of mesh points contains those at the
preceding level.
As an example, we illustrate the `2a + u' re nement for n = 1, 2 in Figures 1, 2.
When n = 2, there are three di erent triangular elements in size. Let  = 1?cos sin  , and

de ne
B0 = f(x; y) 2 S j x + y  =2g;
B1 = f(x; y) 2 S j =4  x + y  =2g;
B2 = f(x; y) 2 S j 0  x + y  =4g:
The set Bl is the union of the triangles of the same size in S , and therefore it is
uniformly divided by the triangulation. The diameter of triangles in Bl is O(2?(2+l)).
4
Figure 1: n = 1

Figure 2: n = 2

5
Moreover, functions in (1.1) and (1.2) are smooth on B0 [ B1.
More generally, by examining the structure of the `La + u' re nement, we can
calculate the total number Nn of triangles at level n: this is (L +1)22n ? 4L = O(22n ).
There are Ln ? (L ? 1) di erent triangular elements in size. The closer the triangle is
to the origin, the smaller it is. As the triangles vary in size from large to small, we name
the region containing triangles of the same size to be B0; B1; : : : ; BLn?L , respectively.
The diameter of triangles in Bl, denoted by l , is O(2?(n+l) ). Let Nl be the number of
triangles in Bl. Then Nl is proportional to 4n?i where l = iL + i1 for 0  i1  L ? 1.
The distance from the origin to Bl, denoted by rl, is O( 2 L1 i ).
( +1)

In each Bl (l = 1; : : : ; Ln ? L), the triangular elements 4K are true triangles and all
are congruent. The triangular elements in B0 are nearly congruent. More importantly,
any symmetric pair of triangles, as shown in Figure 3, have the following property:
v1 ? v2 = ?(v1 ? v4);
v1 ? v3 = ?(v1 ? v5):
The total number of symmetric pairs of triangles in Bl is O(Nl ) and the remaining
p
triangles in Bl is O( Nl ).
If let L = 0, then the re nement is quasi-uniform. The analysis given in [6] indicates
that a quasi-uniform re nement is a better scheme to use with smooth integrands.

3 Interpolation
Let  denote the unit simplex in the st ? plane

 = f(s; t) j 0  s; t; s + t  1g ;
Let 1; : : :; 6 denote the three vertices and three midpoints of the sides of , numbered
according to Figure 4.

6
 v3 2 v
q

 q





  1
v q



 

q q

v4 v5

Figure 3: A symmetric pair of triangles

t
6
2
@@ r

@@
@@
@@
4 @5
@@
r r

@@
@@
@@ -
1
r r
s r

6 3

Figure 4: The unit simplex

7
 To de ne interpolation, introduce the basis functions for quadratic interpolation on
. Letting u = 1 ? (s + t), de ne
l1(s; t) = u(2u ? 1); l2(s; t) = t(2t ? 1); l3(s; t) = s(2s ? 1);
l4(s; t) = 4tu; l5(s; t) = 4st; l6(s; t) = 4su:
We give the corresponding set of basis functions flj;K (q)g on 4K by using its
parametrization over . As a special case of piecewise smooth surfaces in R3, discussed
in [5]-[7], there is a mapping
1{1- 4
mK :  onto (3:1)
K

with mK 2 C 6(). Introduce the node points for 4K by vj;K = mK (j ); j = 1;


; 6.

The rst three are the vertices and the last three are the midpoints of the sides of 4K .
De ne
lj;K (mK (s; t)) = lj (s; t) ; 1  j  6; 1  K  N:
Given a function f , de ne
6
PN f (q) = f (vj;K )lj;K (q) ; q 2 4K ;
X

j =1
for K = 1; : : :; N . This is called the piecewise quadratic isoparametric function inter-
polating f on the nodes of the mesh f4K g for S .

4 Numerical integration and error analyses

1{1- 4 ; we have
With the triangulation f4K g and the mapping mK :  onto K
Z Z

f (q)dS = f (mK (s; t)) j Ds mK (s; t)  Dt mK (s; t) j dsdt (4:1)

4K 

Ds and Dt denote di erentiation with respect to s; t, respectively. The quantity

j Ds mK (s; t)  DtmK (s; t) j is the Jacobian determinant of the mapping mK (s; t) used
in transforming surface integrals over 4K into integrals over . When 4K is a
triangle, the Jacobian is twice the area of 4K .
8
 The numerical integration formula used here is
g(s; t)dsdt  16 [ g(4) + g(5) + g(6) ]
Z

(4:2)

which is based on integrating the quadratic polynomial interpolating g on  at
1; : : :; 6. This integration has degree of precision 2. Applying (4.2) to the right side of
(4.1), we have
6
f (q)dS  61 f (mK (j )) j Ds mK (s; t)  Dt mK (s; t) jj
Z

(4:3)
X

4K j =4
A major problem with (4.3) is that Ds mK and DtmK are inconvenient to compute
for some elements 4K on many surfaces S . Therefore, we approximate mK (s; t) in
terms of only its values at 1; : : : ; 6. De ne
X6 X6
mK (s; t) =
f mK (j )lj (s; t) = vj;K lj (s; t):
j =1 j =1
For the case of the wedge of a circle, in this paper, only the outer triangles will be a ected
by this approximation. Then
6
f (q)dS  16 f (mK (j )) j Ds mK (s; t)  DtmK (s; t) jj
Z
X
f f
4K j =4
6
!j;K f (vj;K )
X
=
j =4
where
!j;K = 61 j Ds mK (s; t)  DtmK (s; t) jj :
f f

LEMMA 1. Let f4K g be the irregular triangulation de ned by the `La + u' re nement,
and let N be its total number of triangles. Assume that f is of form (1.1). Then
6
!j;K f (vj;K ) j  O( N1p )
Z

j B f (q)dS ?
X X
(4:4)
Ln?L 4K BLn?L j =4 1

where p1 = (L+1)(2 +2) :

Proof. Let
1 :
Ln?L = 2n+(Ln?L)

9
The partition of BLn?L is shown in Figure 4. Let f4K j K = 1; : : :; 16g denote the
sixteen triangles in BLn?L , and let 41 contain the origin. Then the area of 4K is
2
Ln?L sin(=2).

y
6
?L
?
? L? L
?
? L ?? LL
L ?
?L ? ?LL? LL
?
?? L ?? LL ?? LL
? ? L?L L?L LL
L ?
?? LL ??? LL ??? LL ??? LL
? L L L L -x
0 Ln?L 2Ln?L 3Ln?L 4Ln?L

Figure 5: The partition of BLn?L

(a) We rst estimate the error of

Z 6
j f (q)dS ? !j;1f (vj;1) j :
X

41 j =4
De ne
g1(x; y) = '()h(r)g(x; y);
and note that g1(x; y) is integrable over 41 . Since r  0, by the integral mean value
theorem we have Z Z

f (q)dS =  r dS
4 1 4 1

where
inf
4 1
g (x; y)    sup g1(x; y):
1 41

10
Also,  is bounded:
j  j  sup j g1(x; y) j
S
 sup j '() j sup j h(r) j sup j g(x; y) j
S S S
= 0max

j '() j 0max
r1
j h(r) j max
S
j g(x; y) j
 M <1
Therefore,
Z Z

j 41
f (q)dS j  M
41
r dS
Z 
Ln?L Z  +1
= M 0 0
r ddr
= M   +2 :
+ 2 Ln?L
On the other hand,
6
j !j;1 f (vj;1) j 16 Ln+2?L sin  (21? + j cos(=2) j ) M
X

j =4

This inequality comes from the following:

!j;1 = 16 j Ds m1(s; t)  Dt m1(s; t) jj :
f f

= 16 2 (Area of 41)
= 61 Ln+2?L sin ; for j = 4; 5; 6
and
j f (vj;1) j  j vj;1 j M 8
>
< ( Ln2?L ) M if j = 4; 6
= >
: Ln?L j cos(=2) j M if j = 5
Hence, Z 6
j f (q)dS ? !j;1f (vj;1) j  O(Ln+2?L ):
X

41 j =4

11
Notice that Ln?L = O( N L1 ). It follows that the error over 41 is O( Nlp ).
+1
2
1

(b) The error over 4K (K = 2; : : : ; 16) can be obtained by using Taylor's error
formula. Since f is smooth on 4K and mK :  onto 1{1- 4 is also smooth, we have
K
6
f (mK (s; t)) ?
X
f (mK (j ))lj (s; t) = Hf;K (s; t; ; ) (4:5)
j =1

where
2 3

4
@ + t @ )3f (m (; )) ? 6 (s @ + t @ )3f (m ( ;  ))l (s; t) :
Hf;K (s; t; ; ) = 3!1 (s @s
X
5
K j j K j j j
@t j =1 @s @t
In this, j = (sj ; tj ), (; ) is on the line segment from (0,0) to (s; t), and (j ; j ) is on
the line segment from (0,0) to (sj ; tj ). Notice that (; ) and (j ; j ) belong to . For
(s; t) 2 ; mK (s; t) 2 4K , and

0 < cos( 2 )Ln?L  r(mK (s; t))  4Ln?L

where r(mK (s; t)) = j mK (s; t) j; the distance from the point mK (s; t) to (0,0).
We would like to examine one term in Hf;K , that associated with s3, and it will
show the general behavior of Hf;K .
@ 3 f (m (; )) = @ 3 (r (m (s; t))
g (m (s; t)))
@s3 K @s3 K 1 K s=
t =
@ 3 (r (m (s; t)))
= g1(mK (; )) @s (4.6)
3 K
s=
t =
@
(
@ 2 )

+ 3 @s g1(mK (s; t)) @s2 (r (mK (s; t))) s =

(4.7)
t =
@ 2 g (m (s; t)) @ (r (m (s; t)))
( )

+ 3 @s 2 1 K @s K
s =
(4.8)
t =
@ 3 g (m (s; t))
+ r (mK (; )) @s (4.9)
3 1 K s=
t =
12
The magnitudes of (4.6), : : : , (4.9) vary from O(Ln?L ) to O(Ln+3?L ), and O(Ln?L )
is the dominant term in (4.6), : : : , (4.9). It follows that the coecient of s3 in Hf;K
is of O(Ln?L ). Consequently,
j Hf;K (s; t; ; ) j  O(Ln?L )
for any (s; t) 2 .
Therefore,
Z 6
j f (q)dS ? !j;K f (vj;K ) j
X

4K j =4
Z 6
 O(Ln
2 ) j f (m (s; t)) ? f (mK (j ))lj (s; t) j dS
X
?L  K
j =1
Z

 Ln?L ) j Hf;K (s; t; ;  ) j dS

O(2

 O(Ln?L ):
+2

Hence, (4.4) holds.

LEMMA 2. Under the assumptions in Lemma 1, let = + 2 + ?2 . Then
L
8
Z 6 >
< O(2?4n?l ) if < 3
j f (q)dS ? !j;K f (vj;K ) j
X X
(4:10)
Bl 4K Bl j =4 >
:
O(2?4n? Ll ?5l) if  3
for l = 1; : : : ; Ln ? L ? 1.
Proof. There are two types of triangles in Bl . Those triangles that are part of symmetric
pairs of triangles (cf. Figure 3) are of the rst type and remaining triangles are of the
second type. By analysis of the `La + u' re nement in section 2, the number of triangles
of the rst type is O(Nl ), where Nl is the number of triangles in Bl. Then Nl is
proportional to 4n?i , where l is decomposed as l = iL + i1 for 0  i1  L ? 1.
Since f is smooth on Bl , then by Taylor's error formula we have
6
f (mK (s; t)) ?
X
f (mK (j ))lj (s; t) = Hf;K (s; t) + Gf;K (s; t; ; ) (4:11)
j =1

13
where
2 3

4
@ + t @ )3f (m (0; 0)) ? 6 (s @ + t @ )3f (m (0; 0))l (s; t) ;
Hf;K (s; t) = 3!1 (s @s
X
5
K j j K j
@t j =1 @s @t
2 3
@ + t @ )4f (m (; )) ? 6 (s @ + t @ )4f (m ( ;  ))l (s; t) :
Gf;K (s; t; ; ) = 4!1 (s @s
X
4
K j j K j j j 5
@t j =1 @s @t
We examine (4.11) and we can nd the following.
First, Hf;K (s; t) is a polynomial of degree three. Second, the coecients of Hf;K (s; t)
are combinations of (v2;K ? v1;K ) and (v3;K ? v1;K ). For instance, the coecient of s3
is
1 @ 3 f (m (0; 0))
3! @s3 K
@ 3 f (m ( ))(v ? v )3 + 3 @ 3 f (m ( ))(v ? v )2(v ? v )
= 3!1 @x
"

K 1 3;x 1;x
3 @x2@y K 1 3;x 1;x 3;y 1;y
@ 3 f (m ( ))(v ? v )(v ? v )2 + @ 3 f (m ( ))(v ? v )3) ; (4.12)
#

+ 3 @x@ 2y K 1 3;x 1;x 3;y 1;y

@y3 K 1 3;y 1;y
where vj;K = (vj;x; vj;y ) for j = 1; 2; 3: For every symmetric pair of triangles, say 41
and 42, in Bl (see gure 3), let
m1(s; t) = (v3 ? v1)s + (v2 ? v1)t + v1
m2(s; t) = (v5 ? v1)s + (v4 ? v1)t + v1
vj = (vj;x; vj;y ):
Then
v1 ? v2 = ?(v1 ? v4)
v1 ? v3 = ?(v1 ? v5) :
We now have for Hf;1 and Hf;2, that the coecient of s3 in Hf;1 is
1 @ 3 f (m ( ))(v ? v )3 + 3 @ 3 f (m ( ))(v ? v )2(v ? v )
"

3! @x3 1 1 3;x 1;x @x2@y 1 1 3;x 1;x 3;y 1;y

14
 @ 3 f (m ( ))(v ? v )(v ? v )2 + @ 3 f (m ( ))(v ? v )3 ; (4.13)
#

+ 3 @x@ 2y 1 1 3;x 1;x 3;y 1;y

@y3 1 1 3;y 1;y
and the coecient of s3 in Hf;2 is
1 @ 3 f (m ( ))(v ? v )3 + 3 @ 3 f (m ( ))(v ? v )2(v ? v )
"

3! @x3 2 1 5;x 1;x @x2@y 2 1 5;x 1;x 5;y 1;y

@ 3 f (m ( ))(v ? v )(v ? v )2 + @ 3 f (m ( ))(v ? v )3 : (4.14)
#

+ 3 @x@ 2y 2 1 5;x 1;x 5;y 1;y

@y3 2 1 5;y 1;y
Adding (4.13) and (4.14) gives us zero, and an analogous argument holds for the
remaining coecients. This means that cancellation happens on any symmetric pair of
triangles. It follows that
Z 2 6 2 Z

j f (q)dS ? !j;K f (vj;K ) j  O(l2) j Gf;K (s; t; ; ) j dsdt

X X X

41 [42 K =1 j =4 K =1 

where l is the diameter of 41 and 42. This argument is based on Chien [6].
We now bound Gf;K (s; t; ; ) on 4K  Bl, for K = 1; 2. Since 1  r(mK (s; t)) 
rl > 0 for (s; t) 2 , and rl = O(2?(L+1)i ), we have for < 4,
6
j Gf;K (s; t; ; ) j  O(4)fO(r(m ?4 ) + X O(r(m ( ;  )) ?4 )g
l K (;  )) K j j
j =1
 O(l4)O(rl ?4)
 O(l4)O(2(4? )(L+1)i):
The rst inequality arises from bounding the individual terms of @ as@ b t ; a + b = 4 , based
@4f

on the kind of expansion done in (4.12). Otherwise, for  4,

j Gf;K (s; t; ; ) j O(l4):

So, 8
Z 6 O(l6)O(2(4?
>
< )(L+1)i) if < 4
j f (q)dS ? !j;K f (vj;K ) j 
X

41 [42 j =4 O(l6) >

: otherwise
15
 The error contributed by triangles of the rst type is as follows. If < 4,
Z 6
j f (q)dS ? !j;K f (vj;K ) ) j  O(Nl)O(l6)O(2(4? )(L+1)i )
X X
(
4K offirst type 4K j =4

 O(22n?2i 26n1+6l )O(2(4? )(L+1)i)

 O(2?4n?l ): (4.15)
For  4,
Z 6
j f (q)dS ? !j;K f (vj;K ) ) j  O(Nl )O(l6)
X X
(
4K offirst type 4K j =4
 O(2?4n? Ll ?6l):
2
(4.16)
For the remaining triangles in Bl, those of the second type, use (4.5), the fact that
p
rl  r(mK (s; t))  1 for 4K  Bl, and that their number is O( Nl). Then
8
>
< O(l3)O(rl ?3 ) if < 3
j Hf;K (s; t; ; ) j >
: O(l3) otherwise
and 8
Z 6 O(2?4n?l ) if < 3 >
<
j ( f (q)dS ? !j;K f (vj;K ) ) j (4:17)
X X

4K ofsecond type 4K j =4 O(2?4n? Ll ?5l) otherwise

>
:

The total error over Bl is given by (4.10).

LEMMA3. Under the assumptions in Lemma 1,
6
!j;K f (vj;K ) j O( N12 ):
Z

j B f (q)dS ? (4:18)
X X

0 4K B j =4 0

Proof. The function f is smooth on B0, and B0 is uniformly divided by triangular

elements. By the results in [6], (4.18) follows.
Combining the above lemmas, we get the following result, which gives the total error
of integrating over S .
THEOREM 1. Let f be of the form (1.1). Then
8
Z N 6 O( lnN
N ) if p1 = 2
>
<
j S f (q)dS ? !j;K f (vj;K ) j 
X X 2

K =1 j =4 O( N1p ) otherwise >

:

16
where N is the total number of triangles in the triangulation, p1 = ( +2)(2 L+1) , and
p = minfp1; 2g.
Proof. We rst add all errors contributed by each 4K  B1 [


[ BLn?L?1 . For < 3,
if p1 6= 2 (i.e. 6= 0),
?L?1
LnX Z 6
j f (q)dS ? !j;K f (vj;K ) j
X X

l=1 Bl 4K Bl j =4
?L?1
LnX
 O(2?4n?l )
l=1

= O 2?4n 2
? ? 2? (Ln?L) !

1 ? 2?
 O(2?2pn )
= O( N1p )
while
?L?1 6
!j;K f (vj;K ) j O( 2n2n )  O( lnN
LnX Z

j f (q)dS ?
X X

l=1 Bl 4K Bl j =4 N 2 ) for = 0.

When  3, is nonzero. By a similar argument,
?L?1 6
!j;K f (vj;K ) j O( N1p ) for  3.
LnX Z

j f (q)dS ?
X X

l=1 Bl 4K Bl j =4
Combining with Lemma 1 and Lemma 3, we have
Z N 6
j S f (q)dS ? !j;K f (vj;K ) j
X X

K =1 j =4
Z 6
 j f (q)dS ? !j;K f (vj;K ) j
X X

BLn?L 4K BLn?L j =4
?L?1
LnX Z 6
j f (q)dS ? !j;K f (vj;K ) j
X X
+
l=1 Bl 4K Bl j =4
Z 6
+j f (q)dS ? !j;K f (vj;K ) j
X X

B0 4K B0 j =4
 O( N1p ) + O( N1p ) + O( N12 )
1

 O( N1p ) for 6= 0:
17
And it is O( lnN
N ) for
2 = 0. The quantity p1 is de ned in Lemma 1, and p =
minfp1; 2g.
COROLLARY Let L be any positive integer greater than +2 4 ? 1. Then the error for

evaluating the integral over S is O(1=N 2 ). In particular, a `1a + u' re nement gives an
error of O(1=N 2 ) , for any > 0.
THEOREM 2. Let f be of the form (1.2). Then
6
!j;K f (vj;K ) j  O( lnN 1)
Z N
j S f (q)dS ?
X X

N p ) + O ( N2
K =1 j =4
1

where N is the total number of triangles in the triangulation, p1 = ( +2)(2 L+1) .

Proof. The proof is analogous to that of Theorem 1.
COROLLARY Let L be any positive integer greater than +2 4 ? 1. Then the error for

evaluating the integral over S is O( N1 ). For > 0, a `1a + u' re nement still gives an
2

error of O(1=N 2 ).

5 Numerical examples
We give numerical examples using the method analyzed in section 4. The method was
implemented with a package of programs written by K. Atkinson, which is described
in [1], [4]. All examples were computed on a Hewlett-Packard workstation in double
precision arithmetic.
EXAMPLE 1. Let S = f(x; y) 2 R2 j 0  r  1; 0    =2g,

f (x; y) = r ; > ?2: (5:1)

The results are given in Table 1 for = 0:1; 0:5 with L = 1. The column labeled
Order gives the value
ln jEn=En+1 j
pn = ln( N =N ) n+1 n

18
where En in the error at level n. Since the theoretical result shows that the error is
O(1=N 2 ), we expect that pn will converge to 2.
Table 2 gives the results for = ?1 with L = 1 and L = 3. The empirical orders
of convergence pn approach 1 and 2, respectively, as expected.

Table 1: Evaluation with L = 1

n N = :1 = :5
Error Order Error Order
1 4 -2.36E-4 -9.47E-4
2 28 -9.00E-6 1.68 2.90E-6 3.10
3 124 3.97E-7 2.10 3.22E-6 -.23
4 508 7.13E-8 1.22 2.91E-7 1.71
5 2044 6.95E-9 1.67 2.10E-8 1.89
6 8188 5.70E-10 1.81 1.40E-9 1.95

Table 2: Evaluation = ?1
n N L=1 N L=3
Error Order Error Order
1 4 1.49E-1 4 1.49E-1
2 28 4.01E-2 0.68 52 9.51E-3 1.08
3 124 1.07E-2 0.93 244 5.61E-4 1.83
4 508 2.53E-3 0.98 1012 3.35E-5 1.98
5 2044 6.32E-4 1.00 4084 2.00E-6 2.02
6 8188 1.58E-4 1.00

19
6 Generalization
We have presented results for only the planar wedge, while using polynomial interpolation
of degree two to approximate the integrand and the integration region. Any other degree
of interpolation could also have been used. In such cases, the de nition of the nodes will
change appropriately. But the de nition of the triangulation will remain the same, and
we will use the `La + u' re nement.
In addition to using quadratic interpolation, we have also examined the use of linear,
cubic and quartic interpolation. For linear and cubic interpolation, the function value at
the origin is needed in the numerical integration. We simply let f (0; 0) = 0. The results
are consistent with the kind of results we have obtained for the quadratic case, and they
are as follows.
Suppose that we use interpolation of degree d to approximate both the integrand
and the wedge S . Let fq1; : : : ; qv g be the node points in the unit simplex and let
fl1; : : : ; lvg be the basis functions in the Lagrange form, where v = (d+1)(2 d+2) . The
points fq1; : : :; qv g will be equally spaced over  and of the form ( di ; dj ); 0  i + j  d:
De ne the interpolating operation
v
PN h(s; t) = h(qj )lj (s; t)
X

j =1
The numerical integration we use is based on integrating this interpolation polynomial,
i.e.
Z Z

h(s; t)dsdt  PN h(s; t)dsdt

 
v
!j(d)h(qj )
X
=
j =1

Assume that the integrand is of the form (1.1). Then with the `La + u' re nement,
we have 8

O( lnN 
Nn v
dn ) if p1 = d
Z >

j S f (q)dS ?
X X
(d)
!j;K f (vj;K ) j  N n
<

K =1 j =1 O( N1np ) otherwise
>
:

20
where Nn the number of triangular elements in the triangulation. The order p =
minfp1; dg, p1 = ( +2)(2 L+1) , d = d+2  d+1
2 when d is an even number, and d = 2 when
d is an odd number. The proof is completely analogous to that given earlier for the
quadratic case. In addition, we need the results for smooth integrands as stated in [6].
The results can be generalized to other integration regions S , for example, a wedge
with the central angle larger than  , triangles, squares, regions containing the origin
as an interior point, etc. We also can generalize this to curved surfaces in R3, with the
integrand having a point singularity of the type given in (1.1). We omit statements of
these results, as they are straightforward.

21
References
[1] K. Atkinson(1985) Piecewise polynomial collocation for boundary integral equations
on surfaces in three dimensions, Journal of Integral Equations. 9(suppl.), pp. 25-48.
[2] K. Atkinson(1985) Solving integral equations on surfaces in space, in Constructive
Methods for the Practical Treatment of Integral Equations, ed. by G. Hammerlin and
K. Ho man, Birkhauser, Basel, pp. 20-43.
[3] K. Atkinson(1989) An Introduction to Numerical Analysis, 2nd ed., John Wiley and
Sons.
[4] K. Atkinson(1992) User's guide to a boundary element package for solving integral
equations on piecewise smooth surfaces, Dept of Mathematics, Univ. of Iowa.
[5] K. Atkinson and D. Chien(1992) Piecewise polynomial collocation for boundary inte-
gral equations, Reports on Computational Mathematics, #29, Dept of Mathematics,
Univ. of Iowa, Iowa City, Iowa.
[6] D. Chien(1991) Piecewise Polynomial Collocation for Integral Equations on Surfaces
in Three Dimensions, PhD thesis, Univ. of Iowa, Iowa City, Iowa.
[7] D. Chien(1992) Piecewise polynomial collocation for integral equations with a
smooth kernel on surfaces in three dimensions, Journal of Integral Equations and
Applications, to appear.
[8] J. Lyness(1976) An error functional expansion for n-dimensional quadrature with
an integrand function singular at a point, Math. Comp. 30, pp. 1-23.
[9] J. Lyness(1990) Extrapolation-based boundary element quadrature, Rend. Semin.
Mat.,Torino., to appear.

22
 [10] C. Schwab and W. Wendland(1992) On numerical cubatures of singular surface
integrals in boundary element methods, Numerische Math. 62, pp. 343-369.

23

View publication stats