# Boundary element collocation method for solving the exterior Neumann problem for Helmholtz's equation in three dimensions.

AMS subject classifications. 35J05, 45B05, 65N35, 65N38

1. Introduction. Many applications in physics deal with the Helmholtz equation in three dimensions. One specific example is the exterior Neumann problem. There are different approaches to solve this partial differential equation. Two commonly used approaches are finite differences and finite elements. However, the given domain is of infinite extent and the Sommerfeld radiation condition has to be satisfied. One can avoid these problems using a boundary integral equation. In addition, the dimensionality is reduced by one. The integral equation approach is the most widely used method to solve the Helmholtz equation. However, a boundary integral equation based on Green's representation theorem or based on a layer approach will lack uniqueness for certain wave numbers.

Fortunately, there exist different variations and modifications of the boundary integral equation to overcome this problem. The combined Helmholtz integral equation formulation (CHIEF) due to Schenck [36] overdetermines the integral equation with the interior Helmholtz integral formulation by choosing strategically as few interior points as possible. For numerical results and CHIEF point selection refer to Seybert et al. [39], and Seybert and Rengarajan [40], respectively. However, in general the choice of those interior points is not clear.

Another boundary integral equation formulation is due to Burton and Miller [10, 11,12]. They cleverly combine the Helmholtz representation formula with its normal derivative and give an idea for the existence and uniqueness proof. A complete proof with appropriate space settings is given by Lin [24]. However, one of the integral operators is hypersingular and usually Holder spaces have to be considered, which complicates the analysis of the boundary element collocation method. The first attempt to solve the boundary integral equation numerically has been made by Burton [11] in 1976. He used the Maue and Mitzner transformation to deal with the hypersingular operator. It also can be removed by regularization which results in a product of two surface integrals, where the kernels are now weakly singular. A mini and Wilton [1] presented numerical results for a sphere and an ellipsoid in 1984 and Liu and Rizzo [28] illustrated numerical results for a sphere in 1992.

Jones [20] and Ursell [37, 38] introduced the theory of modifying the fundamental solution. They added radiating spherical wave functions to the fundamental solution to ensure the unique solvability of the boundary integral equation. Various articles derived coefficients of these added terms to ensure different criteria for a perturbation of a sphere. Two of them are due to Kleinman and Roach [23] and Kleinman and Kress [22], respectively. However, the choice of the coefficients for general surfaces is still in question. For some numerical results we refer to the article by Lin and Warnapala-Yehiya [27].

Numerical results for the T-matrix method developed by Waterman in 1969 [44] have been given by Tobocman [42]. However, this method has some numerical difficulties; see [42] for a discussion. Numerical results for prolate spheroids are presented in [43]. Related work is given by Martin [29].

The boundary integral equation derived in 1965 by Panich [34] uses a combination of a single-double layer including a regularization technique. His formulation also results in a product of two surface integrals, and Holder spaces have to be considered; see also Silva, Power and Wrobel [41] for the smoothness requirements. Panich did not get the desired attention, since his article is written in Russian.

An extension of Panich's method is known as the "modified acoustic single-double layer approach", which we call the "modified Panich method" (MPM). This method has been stated in [15] and [35]. No analysis and numerical results have yet been reported for this extension. The major advantage of the MPM is that we can use the space of continuous functions; that is, this approach does not require any Holder space settings, which would be more restrictive.

For the MPM we use a boundary element collocation method, since superconvergence is observed at the collocation nodes. Atkinson and Chien [5] prove superconvergence at the collocation points for a boundary integral equation solving the Laplace equation using quadratic interpolation. Chien and Lin [13] extend their idea and prove superconvergence at the collocation points for a boundary integral equation for all wave numbers that solves the exterior Dirichlet problem for the Helmholtz equation using quadratic interpolation. They also prove superconvergence for a boundary integral equation based on Green's formula that solves the exterior Neumann problem for the Helmholtz equation using quadratic interpolation. However, their integral equation will break down for certain wave numbers.

Based on these results, we first prove superconvergence at the collocation points for an integral equation based on a single layer formulation that solves the exterior Neumann problem for the Helmholtz equation, distinguishing the cases of even and odd interpolation. Although, we are able to prove superconvergence, this integral equation breaks down if the wave number is an interior Dirichlet eigenvalue. Unlike Atkinson and Chien [5] and Chien and Lin [13], we use interior collocation nodes as in Atkinson and Chandler [4]. As a byproduct, we also obtain superconvergence of an integral equation based on a single layer formulation that solves the exterior Neumann problem for the Laplace equation by choosing the wave number to be zero. Finally, we conjecture superconvergence at the collocation points for the MPM (note that we prove convergence for the MPM in Corollary 3.2) that solves the exterior Neumann problem for the Helmholtz equation which we can confirm with numerical results distinguishing the cases of even and odd interpolation. Note that we are able to prove superconvergence under a strong condition. This condition is needed in our proof due to the nature of the composition of three integral operators, which is needed to regularize the normal derivative of the double layer--a hypersingular integral operator. We are not able to remove that condition, but we give an observation in the paragraph after Theorem 4.13 as to why the condition might hold. Our proofs are based on ideas of Atkinson and Chandler [4] and Micula [31] who proved superconvergence for the radiosity equation, whose kernel has only a bounded singularity. In addition, Micula also proved superconvergence for the exterior Neumann problem solving the Laplace equation using only constant interpolation; see also [32].

Note that there exist other methods, such as spectral methods, to solve the Helmholtz equation. We mention here the work of Graham and Sloan [19].

Finally note that this overview is by far not complete. Some of the most recent results are given by Antoine and Darbas [2] for example, who use an alternative integral equation for smooth surfaces which can be viewed as a generalization of the usual Burton and Miller approach.

Our new approach is limited to smooth surfaces, but we would like to emphasize that, surprisingly, our approach still works for a cube (a polyhedral domain) using constant interpolation. However, there are some more recent publications dealing with uniquely solvable integral equations for the exterior Dirichlet and Neumann problem for Lipschitz boundaries in appropriate Sobolev spaces given by Buffa, Hiptmair and Sauter in [8, 9] and some theoretical results by Betcke et al. in [7]. Further results are presented by Engleder and Steinbach in [17, 18] and by Meury in [30].

The outline of this article is as follows. Section 2 gives the problem formulation of the exterior Neumann problem for solving Helmholtz's equation. The integral equation based on the MPM is reviewed as well as an existence and uniqueness result. Section 3 explains the boundary element collocation method and we give a convergence and error analysis; that is, we review the consistency, stability and convergence order of the boundary element collocation method. In the next section we first prove superconvergence for an integral equation based on a single layer formulation that solves the exterior Neumann problem for the Helmholtz equation. Although we are able to prove superconvergence, the integral equation breaks down if the wave number is an interior Dirichlet eigenvalue. In addition, we are able to prove superconvergence for the integral equation based on a single layer formulation that solves the exterior Neumann problem for the Laplace equation. Then we prove superconvergence at the collocation points distinguishing the cases of even and odd interpolation under a strong condition. In Section 5 numerical results for several smooth surfaces are presented which are in agreement with the theoretical results. A short summary concludes this article.

2. The exterior Neumann problem for solving Helmholtz's equation. Let D be a bounded open region in [R.sup.3]. The boundary of is denoted by [GAMMA] and is assumed to consist of a finite number of disjoint, closed, bounded surfaces belonging to class [C.sup.2], and we assume that the complement [R.sup.3]\[bar.D] is connected; see [14, p. 32].

The mathematical formulation of the exterior Neumann problem consists of finding a complex-valued solution u [member of] [C.sup.1] ([R.sup.3]\D)[intersection] [C.sup.2] ([R.sup.3]\D) solving the Helmholtz equation

[DELTA]u(A) + [k.sup.2]u(A) = 0, A [member of] [R.sup.3]\[bar.D], Im k [greater than or equal to] 0

with the Neumann boundary condition

[absolute value of u]/[absolute value of v](x) = f(x), x [member of] [GAMMA],

where f is a given continuous function on the surface [GAMMA] and u(x) satisfies the Sommerfeld radiation condition

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where r = [absolute value of x] and the limit holds uniformly in all directions x/[absolute value of x]. First, define the acoustic single layer integral operator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the acoustic double layer integral operator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[PHI].sub.k](x,y) = exp(ikr)/4[pi]r, with r = [absolute value of x- y] for x,y [member of] [R.sup.3] x [not equal to] y is the fundamental solution of the Helmholtz equation, and [sigma] [member of] C([GAMMA]). Next, define the acoustic single and double layer integral operators acting on the boundary

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Both operators are compact from C([GAMMA]) to [C.sup.0,[lambda]]([GAMMA]). In fact, the first operator is compact from [C.sup.0,[lambda]]([GAMMA]) to [C.sup.1,[lambda]] ([GAMMA]). Their normal derivatives are defined, respectively, by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first operator is compact from C([GAMMA]) to [C.sup.0,[lambda]] ([GAMMA]), whereas the second operator is bounded from [C.sup.1,[lambda]]([GAMMA]) to [C.sup.0,[lambda]] ([GAMMA]); see [45] or for more general [C.sup.1,[lambda]] ([GAMMA]) [greater than or equal to] 0, see [26].

The problem at hand can be solved with the aid of integral equations. We use the "modified Panich method" (MPM) to derive an integral equation that solves the problem at hand. This approach has been stated in [15] and is an extension of Panich's method [34].

We can write u(A) as a combination of a single and double layer combination in the form

(2.1) u(A) = ([V.sub.k]+ i[eta][W.sub.k][V.sup.2.sub.0] [[sigma]](A), A [member of] [R.sup.3]\[bar.D]

Take the normal derivative of (2.1); let A [right arrow] x [member of] [GAMMA] and use the jump relations to obtain

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which has to be solved for the unknown density function [sigma](x) on the surface [GAMMA]. The parameter [eta] [member of] R [eta] [not equal to] 0 such that [eta] Re k [greater than or equal to] 0 ensures uniqueness for every wave number satisfying Im k [greater than or equal to] 0. The existence and uniqueness proof is given in [15]. Note that the operator [M.sup.T.sub.k] + i[eta][N.sub.k][L.sup.2.sub.0] is compact from C[GAMMA] to C[GAMMA]. In fact, it is also compact from [L.sup.[infinity]]([GAMMA]) to C([GAMMA]); see [6] for the details.

To remove the hypersingularity of the operator [N.sub.k], we use the identity (see [15, p. 43])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and therefore rewrite the Fredholm integral equation of the second kind (2.2) in the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(2.3)

where the kernels of the operators [L.sub.0], [M.sup.T.sub.k], and [N.sub.k] - [N.sub.0] contain only weak singularities for which numerical approximations can be constructed. Finally, solve

u(A) = ([V.sub.k] + i[eta][W.sub.k][V.sup.2.sub.0])[[sigma]](A),

to obtain u(A) for any point A in the exterior domain.

REMARK 2.1. Panich [34] seeks a solution in the form

u(A) = ([V.sub.k] + i[eta][W.sub.k][V.sub.0])[[sigma]](A), A [member of] [R.sup.3]\D,

where the density is found by solving the Fredholm integral equation of the second kind

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

However, we need [sigma] [member of] [C.sup.1,[lambda]] ([GAMMA]), which would be more restrictive.

3. The boundary element collocation method. The boundary element method is discussed extensively in [4]. We briefly summarize the important parts.

Assume that [GAMMA] is a connected smooth surface of class [C.sup.2]; that is, [GAMMA] can be written as

(3.1) [GAMMA]= [[GAMMA].sub.1] [union] ... [union][[GAMMA].sub.J],

where each [[GAMMA].sub.j] is divided into a triangular mesh and the collection of those is denoted by

(3.2) [T.sub.n] = [[DELTA].sub.k]|

Let the unit simplex in the st-plane be defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For a given constant [alpha], with 0 < [alpha] < 1/3, let

(3.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

be the uniform grid inside [sigma] with [f.sub.r]= (r + )(r+2)/2 nodes. We use interior points to avoid the problem of defining the normal [v.sub.x] at the collocation points which are common to more than one face [[DELTA].sub.k]; see [4, p. 280]. The ordering of this grid is denoted by the nodes [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The interior nodes for constant, linear and quadratic interpolation are illustrated in Figure 3.1 and we explain later why we choose such [alpha]'s.

For each [[DELTA].sub.k], we assume there is a map

(3.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which is used for interpolation and integration on [[DELTA].sub.k]. Define the node points of [[DELTA].sub.k] by

[v.sub.k,j] = [m.sub.k]([q.sub.j]), j = , ... ,[f.sub.r].

To obtain a triangulation (3.2) and the mapping (3.4), we use a parametric representation for each region [[GAMMA].sub.j] of(3.1). Assume that for each [[GAMMA].sub.j], there is a map

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [R.sub.j] is a polygonal region in the plane and [F.sub.j] is sufficiently smooth. That means, a triangulation of [R.sub.j] is mapped onto a triangulation [[GAMMA].sub.j]. Let [[??].sub.k,j] be an element of the triangulation of [R.sub.j] with vertices [[??].sub.1,k], [[??].sub.2,k] and [[??].sub.3,k]. Then the map (3.4) is given by

(3.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with the obvious definition of the map [T.sub.k]. We collect all triangles of [R.sub.j] for all j together and denote the triangulation of the parametrization plane by

(3.7) [[??].sub.n] = {[[??].sub.k]| 1 [less than or equal to] k [less than or equal to] n}

and the mesh size by

(3.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which satisfies [[??].sub.n] [right arrow] 0 as [right arrow][infinity].

[FIGURE 3.1 OMITTED]

Most smooth surfaces can be decomposed as in (3.1). In the sequel we consider conforming triangulations satisfying T3; see [3, p. 188]. That is, if two triangles in [??] have a nonempty intersection, then that intersection consists of either (i) a single vertex, or (ii) all of a common edge. Note that T1 and T2 are automatically satisfied, since our surface is assumed to be smooth. The refinement of [[??].sub.k] [member of] [[??].sub.n] is done by connecting the midpoints of the three sides of [[??].sub.k] yielding four new triangles. Thus, T3 is automatically satisfied and this also leads to symmetry in the triangulation and cancelation of errors occurs; see [3, p. 173].

For interpolation of degree r on [sigma], let u = 1 -s -t and the corresponding Lagrange basis functions of degree r on [sigma] are obtained by the usual condition

(3.9) [l.sub.i]([q.sub.i]) = 1, and [l.sub.i]([q.sub.j])= 0, if i [not equal to] j.

In Table 3.1 we state the nodes and Lagrange basis functions over [sigma] for constant, linear and quadratic interpolation.

The interpolation operator is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where y = [m.sub.k] (s,t). The interpolation polynomial of degree r is denoted by [u.sub.n]. Note that [P.sub.n] defines a family of bounded projections on [L.sup.[infinity]] ([GAMMA]) with the pointwise convergence

(3.10) [P.sub.n]u [right arrow]u as n [right arrow] [infinity], u [member of] [X.sub.n],

since [[??].sub.n] [right arrow] 0. Here [X.sub.n] denotes a finite dimensional subspace of [L.sup.[infinity]].

Recall that we have to solve a Fredholm integral equation of the second kind

(3.11) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using the map (3.6), equation (3.11) is equivalent to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Define the collocation nodes by

(3.12) {[v.sub.k,j] = {[m.sub.k]([q.sub.j]), k = 1, ... , n,j = 1 , ... , [f.sub.r].

Collectively, we refer to the collocation nodes by [v.sub.k,j] by {[v.sub.i], where i = 1 , ... , [nf.sub.r]. Then substitute the approximated solution [u.sub.n] in (3.11) and force the residual

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

to be zero at the collocation nodes. Thus, we have to solve the linear system of size [f.sub.r] x [nf.sub.r] given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where i = 1, ... , [nf.sub.r].

This can be written abstractly in the following form. To solve the Fredholm integral equation of the second kind

(3.13) ([lambda]I - K)u = f,

we approximate it by solving

(3.14) [P.sub.n]([lambda]I - K)[u.sub.n] = [P.sub.n]f, [u.sub.n] [member of] [X.sub.n].

This will lead to equivalent linear systems. We reformulate (3.14) in the equivalent form

(3.15) ([lambda]I - [P.sub.n]K)[u.sub.n] = [P.sub.n]f, [u.sub.n] [member of] [X.sub.n],

where [u.sub.n] is the solution of (3.14). Note that the iterated collocation solution [[??].sub.n] = 1/[lambda] [f + [Ku.sub.n]] satisfies

(3.16) u - [[??].sub.n] = [([lambda]I - [KP.sub.n]).sup.-1] K(I - [P.sub.n])u

and

(3.17) [P.sub.n][[??].sub.n] = [u.sub.n] [right arrow] un([v.sub.i]) = [[??].sub.n]([v.sub.i]),

where the collocation nodes [v.sub.i] are given in (3.12); see [3, Eq. 3.4.101 on p. 78, Eq. 3.4.81 on p. 72 and p. 82] for details.

Note that we have consistency. That is, we have

(3.18) [parallel]K - [P.sub.n]K[parallel] [right arrow] 0 as n [right arrow] [infinity],

since K. : [L.sup.[infinity]]([GAMMA]) C([GAMMA])) is compact and since (3.10) holds. With (3.18) we can prove stability and convergence. That is, for all sufficiently large n [greater than or equal to] N the operator [([lambda] -[P.sub.n]K).sup.-1] exists as a bounded operator from to [L.sup.[infinity]]([GAMMA]) to C([GAMMA]). Moreover, it is uniformly bounded

(3.19) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the solution of (3.15) and (3.13) we have

(3.20)

Now, we can establish the following result which is an easy extension of [3, Theorem 9.2.1].

THEOREM 3.1. Let [GAMMA] be a smooth surface. Further assume that [GAMMA] is parametrized as in (3.1) and (3.4), where each [F.sub.j] [member of] [C.sup.r+2]. Let K be a compact integral operator from [L.sup.[infinity]] ([GAMMA]) to C([GAMMA])) and assume the equation (3.13) is uniquely solvable for all functions f [member of] C([GAMMA])). Let [P.sub.n] be the interpolation operator of degree r and consider the approximate solution of ([lambda]I - K)u = f by means of the collocation approximation (3.15). Then we obtain

Stability: The inverse operators [([lambda]I - [P.sub.n]K.).sup.-1] exist and are uniformly bounded for all sufficiently large n [greater than or equal to] N.

Convergence: The approximation [u.sub.n] has error

(3.21) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and therefore [u.sub.n] [right arrow]u as [right arrow] [infinity].

Convergence order: Assume u [member of] [C.sup.r+1] ([GAMMA]). Then

(3.22) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[??].sub.n] is the mesh size of the parametrization domain given by (3.8).

Proof. Consider [P.sub.n] as a projection operator from [L.sup.[infinity]([GAMMA])] into itself. By (3.10) and the assumption that [[??].sub.n], we have

[P.sub.n]u [right arrow] u as n [right arrow] [infinity]

for all u [member of] C([GAMMA]). Since K is compact from to [L.sup.[infinity]([GAMMA])] to C([GAMMA]), we have

[parallel]K - [P.sub.n]K[parallel] [right arrow] 0 as n [right arrow] [infinity].

The existence and stability of [([lambda]- [P.sub.n]K).sup.-1] is based on (3.19), and the error formula (3.21) is simply (3.20). The formula (3.22) is a consequence of [3, p. 165].

As a consequence we directly have the following corollary.

COROLLARY 3.2. Let the parametrization function [F.sub.j] [member of] [C.sup.r+2] and [sigma] [member of] [C.sup.r+1] ([GAMMA]). Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (2.3) obtained through the MPM.

Next, we will show that we can derive a better result; that is, superconvergence at the collocation nodes.

4. Superconvergence of the boundary element collocation method. Here we show that we can improve the result (3.22), distinguishing the cases of even and odd interpolation.

We confine ourselves to triangles in the parametrization plane and then use the map [F.sub.j]. Therefore, let [tau] [subset] [R.sup.2] an arbitrary triangle with vertices {[[??].sub.1], [[??].sub.2], [[??].sub.3]}. If f [member of] C([tau]), then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is a polynomial of degree r in the parametrization variables s and t that interpolates f at the nodes [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] where [q.sub.i] and [l.sub.i] are given in (3.3) and (3.9) and

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which corresponds to the map [T.sub.k] given in (3.6) by suppressing the index k. We will write explicitly [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] if necessary. The operator norm is given by

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The integration formula over [tau] given by

(4.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

has degree of precision of at least r. If r is even, this implies that whenever [[tau].sub.1] and [[tau].sub.1] are triangles for which [[tau].sub.1] [union] [[tau].sub.2] is a parallelogram, then (4.3) has degree of precision r+1; see [31, pp. 22-24]. If r is odd, then (4.3) has degree of precision r. Suppose we can find [alpha] = [[alpha].sub.0] such that (4.3) has degree of precision r+1, then (4.3) has degree of precision r+2 over a parallelogram. For example using [alpha] + 1/6 for r =1 yields degree of precision two over a triangle and degree of precision three over a parallelogram; see [4, p. 271]. For further discussion regarding this matter refer to [31, pp. 58-67].

For differentiable functions f, we define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In our case the kernel is given by K(P,Q) with points P = ([??],[??]) and Q = (x,y). For simplicity we write instead [K.sub.P](x,y) instead of K(P,Q) = K(P, (x,y)).

The following lemma has been used in [31, Proof of Theorem 3.3.16], although it has not been stated or proved.

LEMMA 4.1. Let [tau] be a planar right triangle. Further, assume that the two sides which form the right angle have length h. Let f [member of] [C.sup.r+1]([tau]) and [K.sub.P] [member of] [L.sup.1]([tau]). Then

(4.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where P [not member of] [tau].

Proof. Let [p.sub.r] (x,y) be a Taylor polynomial of f with degree r over [tau]. We have

(4.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for a suitable constant c. Then, we use (4.5) to get the estimate (4.4); see [21, Lemma 3.4.22] for details.

The result of Lemma 4.1 can be extended to general triangles. However, the derivatives of f and [K.sub.P] will involve the mapping [m.sub.[tau]]. The bound of (4.4) will depend on a term proportional to some power of

Y([tau]) = h([tau])/h*([tau]),

where h([tau]) denotes the diameter of [tau] and h*([tau]) denotes the radius of the circle inscribed in [tau] and tangent to its sides. Our triangulation [[??].sub.n] = {[[??].sub.n,k]}, n [greater than or equal to] 1 satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

that is, it is uniformly bounded in and therefore, it prevents the triangles [[??].sub.n,k] from having angles which approach 0 as n [right arrow ][infinity]; see [4, p. 276]. Hence, we have the following corollary.

COROLLARY 4.2. Let [tau] be a planar triangle of diameter h, f [member of] [C.sup.r+1]([tau]) and [K.sub.P] [member of] [L.sup.1]([tau]). Then

(4.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where is c(Y([tau])) some multiple of a power of Y([tau])) and P [not member of] [tau].

Proof. Let [tau] be a planar triangle of diameter h with vertices [v.sub.1], [v.sub.2] and [v.sub.3] and let [??] be a planar right triangle with vertices [[??].sub.1], [[??].sub.2] and [[??].sub.3]. Further, assume that the two sides which form the right angle have length h. The map [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus, using a change of variables, we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

One can easily check that the Jacobian of this transformation is simply Area([tau])/Area([tau]), since the Jacobian of the maps T and [??] are 2- Area([tau]) and 2- Area([??]), respectively. That is, the Jacobian is a constant. Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and hence this case can be reduced to the right triangle case. However, [K.sub.P] and f, as well as their derivatives, will depend on the mapping [m.sub.[??]] . In addition, the constant c(Y([tau])) is some multiple of a power of Y([tau]).

Before we prove the next lemma we need the following assumption.

ASSUMPTION 4.3. Let i = 0, be an integer and let [GAMMA] be a smooth [C.sup.2] surface. Let K(P,Q) be the kernel of our integral operator given in the form

(4.7) K(P,Q)= [PSI](P,Q)/[absolute value of P-Q],

where [PSI] is smooth and bounded. Assume

(4.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where c denotes a generic constant independent of P and Q.

Note that we obtain an additional order in the next theorem. The proof is based on the Duffy transformation; see [16].

LEMMA 4.4. Let [tau] be a planar triangle of diameter h. Let f [member of] [C.sup.r+1]([tau]) and the kernel K satisfy Assumption 4.3 for i = 0. In addition, assume that there is a singularity inside of [tau] or on the boundary; that is, P = Q inside of [tau] or on the boundary. Then

(4.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where c(Y([tau])) is some multiple of a power of Y([tau]).

Proof. We can assume without loss of generality the right triangular case, otherwise proceed as in Corollary 4.2. Assume that the singularity occurs inside of [tau], say at P. Connect the vertices of [tau] with P. We obtain three triangles [[tau].sub.1], [[tau].sub.2] and [[tau].sub.3]. The singularity occurs at one of the vertices of each triangle. Without loss of generality we can assume that we deal with [tau] where the singularity sits at the origin; that is, P = (0,0) (use a linear transformation for each triangle [[tau].sub.1], [[tau].sub.2] and [[tau].sub.3]). If P happens to be on the boundary, then we can use the same procedure except that we deal only with two triangles. If P is on a vertex, then we deal only with one triangle.

Let [p.sub.r](x,y) be a Taylor polynomial of f with degree r over [tau]. Using (4.5), and the notation [??]= I - [L.sub.[tau]],and K = [K.sub.P], we have

(4.10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Next, define the transformation from [tau] to the unit square [] = [0,1] x [0,1] by

(4.11) x = (1 - v)uh, y= uvh with Jacobian [uh.sup.2].

Note that this is the composition of the map [tau] to [sigma] given by (x,y) = (sh,th) with Jacobian [h.sup.2] and the Duffy transformation mapping [sigma] to [] given by (s,t) = ((1 -v) u, uv) with Jacobian u. Hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

g(u,v) = 1/[square root of [(1 -v).sup.2] +[v.sup.2]]

is nonsingular together with all its derivatives on []. Therefore g(u,v), [member of] C[infinity]([]), and hence it is bounded on the compact set []. Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This can be done for all triangles [[tau].sub.i], and thus we obtain the assertion (4.9).

Remark 4.5. Note that in [4, Theorem 3.7] it is assumed that all integrals over a triangle containing the singularity are evaluated with an error of [h.sup.4]. A similar assumption is stated in [31, Theorem 3.3.5].

4.1. Interpolation of even degree. In this section, we assume that r is even. This implies whenever [[tau].sub.1] and [[tau].sub.2] are triangles for which [[tau].sub.1] [union] [[tau].sub.2] is a parallelogram, then (4.3) has degree of precision r + 1.

The next lemma is a general statement of [31, Lemma 3.3.15].

LEMMA 4.6. Let [[tau].sub.1] and [[tau].sub.2] be two planar triangles with diameter h such that R = [[tau].sub.1][union] [[tau].sub.2] is a parallelogram. Let f [member of] [C.sup.r+2] (R) and let [K.sub.P] [member of] [L.sup.1] (R) be differentiable with first derivatives [D.sub.x][K.sub.P] and belonging to [L.sup.1](R). Then

(4.12) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where Y(R) = [max.sub.i=1,2] Y([[tau].sub.i]), c(Y(R)) is some multiple of a power of Y(R), and P [not member of] R.

We now consider a symmetric triangulation. We can impose symmetry in the triangulation [[??].sub.n]. We refine a triangle in the parametrization plane by dividing it into four new smaller triangles by connecting the midpoints of the three sides by straight lines. Hence, the refinement of [[??].sub.n] increases by a factor of 4. In addition, most of the triangles can be grouped as parallelograms. The number of such triangles is O(n). O([[??].sup.-2.sub.n]) The number of remaining triangles is O([square root of n]) = O([[??].sup.-1.sub.n]). The symmetric pairs can be chosen such that the remaining elements are at a bounded distance from P (independent of n); see [3, 5].

We want to apply these results to the individual subintegrals

(4.13) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where i = 1,..., [f.sub.r]n and which are now defined in the parametrization plane for some j depending on [v.sub.i]. We define

(4.14) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which are defined in the parametrization plane.

In the following, by f [member of] [C.sup.k]([GAMMA]) we mean and f [member of] C[GAMMA]) and f [member of] [C.sup.k]([[GAMMA].sub.j]) (that is, f o [F.sub.j] [member of] [C.sup.k]([R.sub.j])), j = 1, ... , J. Note that in [31, Theorem 3.3.16] the kernel is smoother than ours; that is, the radiosity equation is considered by Micula, whereas we consider single/double layer potentials which do not fit in his framework. He obtained a rate of O([[??].sup.r+2.sub.n]). Furthermore, we have removed an additional assumption.

In the sequel, we assume that the conforming triangulation [[??].sub.n] is symmetric and satisfies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, we are in position to state the following theorem.

THEOREM 4.7. Assume the conditions of Theorem 3.1 with each parametrization function [F.sub.j] [member of] [C.sup.r+2] and u [member of] [C.sup.r+2] ([GAMMA]). Assume that the kernel K satisfies Assumption 4.3. Then

(4.15) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Using (3.16) and (3.17) we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] see also [3, p. 450]. Thus, we will bound

(4.16) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

to prove (4.15), where we used (4.13) and (4.14) in the last step. By assumption, f is [C.sup.r+2] in the parametrization plane.

For a given collocation point [v.sub.i], denote by [DELTA]* the curved triangle containing this point and denote by [??]* the triangle in the parametrization plane containing the point [P.sub.i] [member of] [??]* satisfying [v.sub.i] = [F.sub.j]([P.sub.i]); see (4.14). Define

[[??]*.sub.n]= [[??].sub.n]- [??]*

and subdivide [[??]*.sub.n] into two disjoint classes [[??].sub.n.sup.(1)] and [[??].sub.n.sup.(2)] such that [[??].sub.n.sup.(1)] [union] [[??].sub.n.sup.(2)] = [[??]*.sub.n], where [[??].sub.n.sup.(1)] denotes the set of triangles making up parallelograms to the maximum extent possible (actually parallelograms in the parametrization plane) and [[??].sub.n.sup.(2)] denotes the set of the remaining triangles. Hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Recall that the number of triangles in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Moreover, all but a finite number of the triangles in [[??].sub.n.sup.(2)] bounded independently of n, will be at a minimum distance from d > 0 with [P.sub.i] with d independent of n and i. Hence the function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] uniformly bounded for the point (x,y) being in a triangle in [[??].sub.n.sup.(2)]. Thus, we can split the integral in (4.16) into three parts

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The bound for the three terms on the right-hand side can be done in the following three steps, respectively.

1. By Lemma 4.4, the error in evaluating the integral

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is O([h.sup.r+2]) where h is the diameter of [??]*. Thus, we also have O([[??].sup.r+2.sub.n])by the definition (3.8).

2. Next, consider the error from triangles in [[??].sub.n.sup.(1)]. By Lemma 4.6 we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

According to Assumption 4.3, the last quantity is bounded by

(4.17) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A use of a local representation of the surface and polar coordinates shows that the expression in (4.17) is of order O([[??].sup.r+2.sub.n])ln ([[??].sup.-1.sub.n])); see Appendix. Therefore, the error arising from triangles in [[??].sub.n.sup.(1)] is O([[??].sup.r+2.sub.n])ln ([[??].sup.-1.sub.n])).

3. Lastly, consider the error over each such triangle in [[??].sub.n.sup.(2)]. Applying Corollary 4.2 yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the last step follows, since the area of the k-th triangle is O([h.sup.2.sub.k])and since [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is uniformly bounded as mentioned above. Therefore, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since we have ([[??].sup.-1.sub.n]) such triangles. Thus, the total error coming from triangles in [[??].sub.n.sup.(2)] is O([[??].sup.r+2.sub.n]).

Combining the errors from the integrals over [??]*, [[??].sub.n.sup.(1)], and [[??].sub.n.sup.(2)] gives the result (4.15).

COROLLARY 4.8. Assume the conditions of Theorem 4.7. Then

(4.18) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where P [member of] [GAMMA].

Proof. The proof is almost identical to the proof given in Theorem 4.7. For an arbitrary point with P [member of] [GAMMA], with P = [F.sub.j](x,y), denote by [[??].sub.n.sup.(0)] the set of triangles in the parametrization plane that contain the point (x,y). There is only a finite number of such triangles (see [3, p. 451]), but not more than six triangles. Define

[[??]*.sub.n] = [[??].sub.n]- [[??].sub.n.sup.(0)]

and subdivide [[??].sub.n] as in the previous proof to obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By Lemma 4.4, the error in evaluating the integral over [[??].sub.n.sup.(0)] is O[[??].sup.r+2.sub.n]. The remaining proof is the same as the proof of Theorem 4.7.

First we prove superconvergence at the collocation nodes for the integral equation

(4.19) (- 1/2 I + [K.sub.1])[[sigma]](x)= f(x),

where [K.sub.1] = [M.sup.T.sub.k]. This integral equation solves the exterior Neumann problem and is a special case of (2.3) by choosing [eta] = 0. It is based on a single layer formulation, but it will break down if k is an eigenvalue of the interior Dirichlet problem; see [12, p. 205].

We now show that Assumption 4.3 is satisfied for the kernel of the operator [K.sub.1] = [M.sup.T.sub.k] and [K.sub.2] = [L.sub.0].

LEMMA 4.9. Assumption 4.3 is satisfied for the kernel of the operator [K.sub.1] and [K.sub.2].

Proof. We only prove that Assumption 4.3 is satisfied for the kernel of the operator [K.sub.1]. For [K.sub.2] refer to [21, Lemma 3.4.16]. First rewrite the kernel of the operator [K.sub.1] as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [v.sub.P] is the outer normal at P. Note that [K.sub.1] (P,Q) is a smooth function. First we consider,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Clearly,

(4.20) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since [absolute value of cos [[theta].sub.p]] _ c[absolute value of P - Q] ;see [33]. Note that the kernel [K.sub.2] has been studied before in [31]. Next, we claim that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We mention that we have to show this statement for a given fixed point P only for Q in a neighborhood [GAMMA] [intersection] B(P; R) of P with 0 < [absolute value of P - Q] < R, since for Q [member of] [GAMMA] not belonging to [GAMMA][intersection]B(P; R) we can always use the estimates [absolute value of P - Q] [less than or equal to] d= diam([GAMMA]), [absolute value of [v.sub.P],Q - P)] [less than or equal to] d and 1/[absolute value of P - Q] [less than or equal to] 1/R, because [absolute value of P - Q] [greater than or equal to] R.

According to the assumptions on the surface, there exists an R > 0 such that [GAMMA][intersection]B(P; R) can be projected onto the tangent plane in a one-to-one fashion. Thus, we can assume that the surface [GAMMA] can be represented locally by z = f (x,y) with f [member of] [C.sup.2]. Without loss of generality we let P be the origin of the coordinate system and Q be an arbitrary point in [GAMMA] [intersection] B(P; R). Therefore, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and implicitly we have f(0,0) = [f.sub.x](0,0) = [f.sub.y](0,0)= 0. Hence

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that [v.sub.P] = [N.sup.P] is independent of x and y. Clearly, we have

(4.21) [absolute value of Q/[[absolute value of Q].sup.2] x [N.sup.P]] c, P [not equal to] Q.

Now, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first term is bounded using Taylor expansion about (0,0) to obtain

(4.22) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first factor of the second term is bounded due to (4.21). The second factor is bounded since [absolute value of [f.sub.x]] is bounded and

(4.23) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Hence,

[absolute value of [partial derivative]/[partial derivative]x [K.sub.2]] [less than or equal to] c/[[absolute value of P - Q].sup.2], P [not member of] Q.

Similarly, we can prove

[absolute value of [partial derivative]/[partial derivative]y [K.sub.2]] [less than or equal to] c/[[absolute value of P - Q].sup.2], P [not member of] Q.

and therefore

[absolute value of [D.sub.Q][K.sub.2]] [less than or equal to] c/[[absolute value of P - Q].sup.2], P [not member of] Q.

Lastly, we show that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Clearly, we have

(4.24) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

because the first factor is bounded by (4.21). Using Taylor series for [e.sup.z], one can easily see that the last factor is bounded, too. Next, consider

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first term is bounded, since we can write it as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and then use (4.22). The second term can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and clearly it is bounded, because of(4.23) and (4.24). The last term is bounded due to (4.21). Hence,

[absolute value of [partial derivative]/[partial derivative]x [K.sub.1]]- c/[[absolute value of P - Q], P [not member of] Q.

Similarly, we can prove

[absolute value of [partial derivative]/[partial derivative]y [K.sub.1]] c/[[absolute value of P - Q], P [not member of] Q.

and therefore

[absolute value of [D.sub.Q][K.sub.1]] c/[[absolute value of P - Q], P [not member of] Q.

In summary, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and thus proves Assumption 4.3.

Now, we state the following superconvergence theorem for the integral equation (4.19).

THEOREM 4.10. Let the parametrization function and [F.sub.j] [member of] [C.sup.r+2] and [sigma] [member of] [C.sup.r+2]([GAMMA]). Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.19), where K = [M.sup.T.sub.k].

Proof. Since the kernel of [K.sub.1] satisfies Assumption 4.3 as shown in Lemma 4.9, we have by Theorem 4.7

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and thus, we have proved the theorem.

We note that Chien and Lin [13] prove superconvergence at the collocation nodes for an integral equation based on Green's formula. They use the vertices and midpoints of a triangle as collocation nodes for piecewise quadratic polynomials. However, we use interior nodes and therefore our analysis is different from theirs.

Note also that if we choose k = 0 in (4.19), then the integral equation

(4.25) (- 1/2 I + [M.sup.T.sub.0])[[sigma]](x) = f(x,)

solves the exterior Neumann problem for the Laplace equation. Thus, we have the following superconvergence result for the Laplace equation.

COROLLARY 4.11. Assume the conditions of Theorem 4.10. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.25), where K = [M.sup.T.sub.0].

Note that the case for constant interpolation with [alpha] = 1/3 has already been considered in [31] and agrees with our results.

Finally, we show superconvergence under a strong assumption for the MPM method. First, denote

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then, equation (2.3) can be written as

(- 1/2 I + K)[[sigma]](x) = f(x),

where

(4.26) K = [K.sub.1] + i[eta][K.sub.4] + i[eta][K.sub.3]- i[eta]1/4 [K.sub.2].

For the integral equation obtained from the MPM given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(4.27)

we have the following conjecture.

CONJECTURE 4.12. Let the parametrization function [F.sub.j] [member of] [C.sup.r+2] and [sigma] [member of] [C.sup.r+2]([GAMMA]).

Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.27), where K. is given by (4.26).

Note that the theoretical rate in Conjecture 4.12 can be confirmed with several numerical results. We can prove the conjecture above under the strong assumption that is stated below in the next theorem.

THEOREM 4.13. Let the parametrization function [F.sub.j] [member of] [C.sup.r+2] and [sigma] [member of] [C.sup.r+2] ([GAMMA]). In addition, assume that there is a constant c such that [absolute value of [c.sub.R]][less than or equal to]c for all points R in [GAMMA] in, where [c.sub.R] is defined by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], which is from Corollary 4.8. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.27), where K is given by (4.26).

Proof. Since the kernel of [K.sub.1] and [K.sub.2] satisfy Assumption 4.3 as shown in Lemma 4.9, we have by Theorem 4.7

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We also have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since we can establish

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The integral with respect to T can be written in the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we used Corollary 4.8, i.e.,[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. By assumption, we have [absolute value of [c.sub.R]] [less than or equal to] c. Thus, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A similar argument yields the result for [K.sub.4].

Note that the assumption on the boundedness of [c.sub.R] for all points R in [GAMMA] seems a little strong. However, we know that K maps from [L.sup.[infinity]] ([GAMMA]) into C([GAMMA]). That is, for any fixed [sigma], K(I - [P.sub.n])[sigma] is a continuous function on [GAMMA] and on each triangular element [[DELTA].sub.k]. Therefore, K(I - [P.sub.n])[sigma] is a bounded function on [GAMMA] and on each triangular element [[DELTA].sub.k]. This seems to suggest that [c.sub.R] is bounded on [GAMMA].

4.2. Interpolation of odd degree. We now assume r is odd. Then (4.3) has degree of precision r. Suppose we can find [alpha] = [[alpha].sub.0] such that (4.3) has degree of precision r + 1. Then (4.3) has degree of precision r+2 over a parallelogram. For example using [alpha] = 1/6 for r=1 yields degree of precision two over a triangle and degree of precision three over a parallelogram.

The following lemma is a general statement of [31, Lemma 3.3.12].

LEMMA 4.14. Let [[tau].sub.1] and [[tau].sub.2] be two planar triangles with diameter h such that R = [[tau].sub.1] [union] [[tau].sub.2] is a parallelogram. Let f [member of] [C.sup.r+3](R) and [K.sub.P] [member of] [L.sup.1] (R) be twice differentiable with derivatives of order 1 and 2 belonging to [L.sup.1](R). In addition, assume [alpha] = [[alpha].sub.0]. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(4.28) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where Y(R)= [max.sub.i=1,2] Y ([[tau].sub.i]), c (Y (R)) is some multiple of a power of Y(R), and P [not member of] R.

Before we prove the next theorem we need the following assumption.

ASSUMPTION 4.15. Let [GAMMA] be a smooth [C.sup.3] surface. Let K(P,Q) be the kernel of our integral operator. Assume

[absolute value of [D.sup.2.sub.Q]K(P,Q)] c/[[absolute value of P - Q].sup.3], P [not member of] Q.

where c denotes a generic constant independent of P and Q.

Note that in [31, Theorem 3.3.14] a rate of O([[??].sup.r+3.sub.n])ln O([[??].sup.-1.sub.n]) has been obtained for the radiosity equation.

Now, we are in position to state the following theorem.

THEOREM 4.16. Assume the conditions of Theorem 3.1 with each parametrization function and [F.sub.j] [member of] [C.sup.r+3] and u [member of] [C.sup.r+3]([GAMMA]). Assume [alpha] = [[alpha].sub.0]. Moreover, assume that the kernel K satisfies Assumption 4.3 and 4.15. Then

(4.29) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. We will bound

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

to prove (4.29). Using the same approaches as done in Theorem 4.7, we have

1. By Lemma 4.4, the error in evaluating the integral over [??]* is O([[??].sup.r+2.sub.n]).

2. Next, consider the error from triangles in [[??].sup.(1).sub.n]. By Lemma 4.14 and using the same argument as in item 2 of Theorem 4.7, we have

(4.30) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

According to Assumption 4.3 and Assumption 4.15 the quantity (4.30) is bounded by

(4.31) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A use of a local representation of the surface and polar coordinates shows that the expression in (4.31) is of order [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Therefore, the error arising from triangles in ([[??].sub.n.sup.(1)]) is O([[??].sup.r+2.sub.n]).

3. Lastly, consider the error over each such triangle in ([[??].sub.n.sup.(2)]). As in Theorem 4.7 item 3, we apply Lemma 4.2. Thus, the total error coming from triangles in ([[??].sub.n.sup.(2)]) is O([[??].sup.r+2.sub.n]). Combining the errors from the integrals over [??]*,([[??].sub.n.sup.(2)]), and ([[??].sub.n.sup.(2)]) gives the result (4.29).

COROLLARY 4.17. Assume the conditions of Theorem 4.16. Then

(4.32) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where P [member of] [GAMMA].

Proof. The proof is almost identical to the proof given in Theorem 4.16. The reasoning is similar to the one given in Corollary 4.8.

First, we will prove superconvergence at the collocation nodes for the integral equation

(4.33) (- 1/2 I + [K.sub.1])[[sigma]](x)= f(x),

where [K.sub.1][M.sup.T.sub.k]. This integral equation solves the exterior Neumann problem and is a special case of (2.3) by choosing [eta]= 0. It is based on a single layer formulation, but it will break down if k is an eigenvalue of the interior Dirichlet problem; see [12, p. 205].

We now show that Assumption 4.15 is satisfied for the kernels of the operator [K.sub.1]= [M.sup.T.sub.k] and [K.sub.2] = [L.sub.0].

LEMMA 4.18. Assumption 4.15 is satisfied for the kernel of the operator [K.sub.1] and [K.sub.2].

Proof. We only prove that Assumption 4.15 is satisfied for the kernel of the operator [K.sub.1]. For [K.sub.2] refer to [21, Lemma 3.4.41]. Recall that we have (see Lemma 4.9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. We have to show that for P [not equal to] Q the absolute value of all second order partial derivatives of [K.sub.1] and [K.sub.2] with respect to x and y are bounded by c/[[absolute value of and P-Q].sup.2], and c/[[absolute value of P- Q].sup.3]respectively. Then, we would obtain

[absolute value of [D.sup.2.sub.Q]K] c/[[absolute value of P - Q].sup.3], P [not member of] Q.

Using a local representation of the surface, we obtain (see Lemma 4.9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(4.34)

Thus, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the remainder of the proof refer to [21, Lemma 3.4.12].

Now, we state the following superconvergence theorem for the integral equation (4.33).

THEOREM 4.19. Let the parametrization function and [F.sub.j] [member of] [C.sup.r+3] and [sigma] [member of] [C.sup.r+3]([GAMMA]). Assume [alpha] = [[alpha].sub.0]. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.33), where K = [M.sup.T.sub.k].

Proof. Since the kernel of [K.sub.1] satisfies Assumption 4.3 and Assumption 4.15 as shown in Lemma 4.9 and 4.18, we have by Theorem 4.16

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and thus, we have proved the theorem.

Note also that if we choose k = 0 in (4.19), then the integral equation

(4.35) (- 1/2 I + [M.sup.T.sub.0])[[sigma]](x)= f(x),

solves the exterior Neumann problem for the Laplace equation. Thus, we have the following superconvergence result for the Laplace equation.

COROLLARY 4.20. Assume the conditions of Theorem 4.19. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.35), where K = [M.sup.T.sub.0]

For the integral equation obtained from the MPM given by

(4.36) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we have the following conjecture.

CONJECTURE 4.21. Let the parametrization function [F.sub.j] [member of] [C.sup.r+3] and [sigma] [member of] [C.sup.r+3] ([GAMMA]). Assume [alpha] = [[alpha].sub.0]. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.36).

Note that the theoretical rate in Conjecture 4.21 can be confirmed with several numerical results. We can prove the conjecture above under the strong assumption that is stated in the next theorem.

THEOREM 4.22. Let the parametrization function [F.sub.j] [member of] [C.sup.r+3] and [sigma] [member of] [C.sup.r+3] ([GAMMA]). Assume [alpha] = [[alpha].sub.0]. In addition, assume that there is a constant c such that [absolute value of [c.sub.R]] [less than or equal to] c for all points R in [GAMMA], where [c.sub.R] is defined by [absolute value of [L.sub.0] (I - [P.sub.n])[sigma](R)] [less than or equal to] [c.sub.R] [[??].sup.r+2.sub.n], which is from Corollary 4.17. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the integral equation (4.36).

Proof. The proof is similar to the proof of Theorem 4.13, since Assumption 4.15 is satisfied for the kernels of the operator [K.sub.1] and [K.sub.2] as shown in Lemma 4.18.

5. Numerical result. First, we illustrate the accuracy of the integral equation (2.3) for constant, linear, and quadratic interpolation. The surface under consideration is the unit sphere. The true solution (see [25]) is given by

(5.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where r = [square root of [x.sup.2] + [y.sup.2] + [z.sup.2]]. The wave number k equals 1. We used [N.sub.S] = 128, [N.sub.NS] = 4 and [eta] = k/2. We denote with n the number of faces of the triangulation and with [n.sub.v] the number of node points of the triangulation.

As we can see in Tables 5.1, 5.2, and 5.3 the closer the point to the boundary the worse the error. However, using constant interpolation we obtain two to three digits accuracy. With the linear interpolation we obtain three to four digits accuracy, whereas for the quadratic interpolation the accuracy is four to five digits. For points situated further away we get even five to six digits.

Note that the integral equation (2.3) does not break down for critical wave numbers. The first zero of the spherical Bessel function of order one is 4.4934094579 and is such a critical wave number. A stable solution is obtained with the integral equation (2.3) for [eta] = k/2 and illustrated in Table 5.4. If we choose [eta]= 0, the uniqueness is not guaranteed which can be seen in Table 5.5. No digit is calculated correctly.

Increasing the wave numbers means that the kernel is more oscillatory and the accuracy depends crucially on the integration routines. In Table 5.6 we present the accuracy for an ellipsoid with a = 1.0,b = 1.2 and c = 1.5 for the wave number k = 7.

Next, we present numerical results to illustrate the superconvergence of the collocation method for smooth surfaces without surface approximation (see [21, Appendix A]) for the integral equation (2.3) and compare them with the theoretical results.

Since we do not know the exact density [sigma] on the surface, we define the estimated order of convergence (EOC) at the collocation point [P.sub.i] by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

To compare the estimated order of convergence at all the collocation points for different refinements, we define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where I is the number of collocation points of the initial triangulation. Note that this kind of comparison can only be done for constant interpolation with [alpha] = 1/3, since the set of collocation nodes for the initial triangulation is a subset of any refined triangulation. Table 5.7 shows the agreement with the theoretical results for 4 and 8 collocation nodes belonging to all different refinements.

To verify the superconvergence also for linear and quadratic interpolation, we pick some points in the exterior domain and calculate the rates, since we know the true solution in the exterior domain.

Let [P.sub.1,](3,3,3), [P.sub.2](2,5, 6) and [P.sub.3](10,11,12) be points in the exterior domain. The true solution is given by (5.1). Denote the error between the calculated solution [u.sub.n] and the true solution u at the point [P.sub.i] by [E.sub.n]([P.sub.i]); that is,

[E.sub.n]([P.sub.i]) = [absolute value of u([P.sub.i])- [u.sub.n]([P.sub.i])].

Define the estimated order of convergence (EOC) at the point [P.sub.i] by

EOC ([P.sub.i]) = [log.sub.2] ([E.sub.n]([P.sub.i])/ [E.sub.4n] ([P.sub.i])).

We consider two different smooth surfaces. As a first example we consider an ellipsoid with a = 1.0,b = 1.2 and c = 1.5. The second surface is peanut-shaped. Its surface in spherical coordinates is given through x, and [??]sin([phi]) cos([theta]),y = [??]sin([phi])sin([theta])and z = [??]cos([phi]), where [[??].sup.2]=9 {[cos.sup.2]([phi] + [sin.sup.2]([phi])/4}/4. As a third surface we consider the acorn surface. Its surface is given by x = [??]sin([phi]) cos([theta]),y = [??]sin([phi])sin([theta])and z = [??]cos([phi]), where [[??].sup.2]= 9 {17/4 + 2 cos(3([phi]))/25.

The estimated error of convergence illustrated in Table 5.8 is in agreement with the superconvergence rates predicted by Theorems 4.13 and 4.22; that is we obtain the theoretical order of convergence , O([[??].sup.2.sub.n]ln ([[??].sup.-1.sub.n]),O([[??].sup.3.sub.n])and O([[??].sup.4.sub.n])ln ([[??].sup.-1.sub.n]), and , respectively.

For the linear and quadratic interpolation case we almost obtain an additional order of convergence compared to the theoretical results predicted by Theorems 4.13 and 4.22. This is possibly due to the smoothing effect of the integral of the density function. We get similar result for points closer to the boundary; see [21].

The numerical results shown in Table 5.9 for a peanut are in agreement with the theoretical results predicted by Theorems 4.13 and 4.22. Note that we do not have an additional order of convergence for the linear interpolation as in the ellipsoidal case. The same is true for an acorn as illustrated in Table 5.10. Finally, note that different wave numbers give similar convergence rates.

Lastly, we consider a cube with edge length two centered at the origin, although it contradicts the assumption that the surface has to be of class [C.sup.2]. The points in the exterior are given by [P.sub.1] = (1.1,1.1,1.1,), [P.sub.2]= (1.5,1.5,1.5,) and [P.sub.3] = (3.0,3.0.,3.0).

Surprisingly, we can obtain quite good results for this polyhedral domain provided we use constant interpolation as illustrated in Table 5.11. It seems like we can almost achieve the order of convergence O ([h.sup.2] ln ([h.sup.-1]))and thus, this shows the applicability of our method even for such polyhedral domain, but further investigation is necessary from the theoretical point of view. This will be possible future research which could be based on ideas of[2, 8, 9, 17, 18]. However, note that we do not have the desired superconvergence for the linear and quadratic case.

6. Summary. We use the boundary element collocation method to solve a Fredholm integral equation of the second kind, where we use interpolation at interior points. We prove superconvergence at the collocation nodes, distinguishing the cases of even and odd interpolation, and illustrate it with numerical results for several smooth surfaces. We also demonstrate numerical results for a cube (a polyhedral surface) not satisfying the smoothness assumption which in theory is necessary to obtain superconvergence. But, surprisingly, our approach still works for a cube using constant interpolation and investigation regarding this observation might be possible future research.

Appendix. In this appendix, we prove the following theorem. THEOREM A.1. We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Without loss of generality, assume [v.sub.i] = (0,0,0). Define

[S.sub.0,R] = {y [member of] [T.sub.n] : [absolute value of y] < R},

such that [DELTA]* [subset] [S.sub.0,R], where R > 0 is sufficiently small (since we can assume diam([DELTA]*) [right arrow]0 as [[??].sub.n])[right arrow] 0. For the integral

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Hence,

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now, consider the integral

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with the constant k > 0 sufficiently small such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is contained in [DELTA]*. Note that we can find such k, since the diameter of the curved triangle [DELTA]* is smaller or equal to the maximum mesh size [[delta].sub.n] of the surface. Also note that the projection of [S.sub.0,R] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] into the parametrization plane are circles centered at zero with radius R and k[[??].sub.n], since the maximum diameter of the triangles in the parametrization plane is by definition [[??].sub.n]. Using polar coordinates we get the following (see [14, p. 40])

(6.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Combining (6.1) and (6.2) yields the result.

Acknowledgment. The authors thank Professor K. E. Atkinson for his source code that solves the Laplace equation by a boundary element method.

* Received February 8, 2011. Accepted for publication January 4, 2012. Published online on April 27, 2012. Recommended by O. Widlund.

REFERENCES

[1] S. Amini and D. T. Wilton, An investigation of boundary element methods for the exterior acoustic problem, Comput. Methods Appl. Mech. Engrg., 54 (1986), pp. 49-65.

[2] X. Antoine and M. Darbas, Alternative integral equations for the iterative solution of acoustic scattering problems, Quart. J. Mech. Appl. Math., 58 (2005), pp. 107-128.

[3] K. E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997.

[4] K. E. Atkinson and G. Chandler, The collocation method for solving the radiosity equation for unoc cluded surfaces, J. Integral Equations Appl., 10 (1998), pp. 253-290.

[5] K. E. Atkinson and D. Chien, Piecewise polynomial collocation for boundary integral equations, SIAM J. Sci. Comput., 16 (1995), pp. 651-681.

[6] K. E. Atkinson, I. Graham, and I. Sloan, Piecewise continuous collocation for integral equations, SIAM J. Numer. Anal., 20 (1983), pp. 172-186.

[7] T. Betcke, S. N. Chandler-Wilde, I. G. Graham, S. Langdon, and M. Lindner, Condition number estimates for combined potential integral operators in acoustics and their boundary element discretisation, Numer. Methods Partial Differential Equations, 27 (2011), pp. 31-69.

[8] A. Buffa and R. Hiptmair, Regularized combined field integral equations, Numer. Math., 100 (2005), pp. 1-19.

[9] A. Buffa and S. Sauter, On the acoustic single layer potential: stabilization and Fourier analysis, SIAM J. Sci. Comput., 28 (2006), pp. 1974-1999.

[10] A. J. Burton, The solution of Helmholtz' equation in exterior domains using integral equations, National physical laboratory report NAC 30, Division of numerical analysis and computing, NPL, Teddington, January 1973.

[11] --, Numerical solution of acoustic radiation problems, National physical laboratory contract OC5/535, Division of numerical analysis and computing, NPL, Teddington, February 1976.

[12] A. J. Burton and G. F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proc. R. Soc. Lond. Ser. A, 323 (1971), pp. 201-210.

[13] D. Chien and T.-C. Lin, Piecewise polynomial collocation for boundary integral equations for solving Helmholtz's equation, in preparation.

[14] D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, John Wiley & Sons, New York, 1983.

[15] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, second ed., Springer-Verlag, New York, 1998.

[16] M. G. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAM J. Numer. Anal., 19 (1982), pp. 1260-1262.

[17] S. Englederand O. Steinbach, Modified boundary integral formulations for the Helmholtz equation, J. Math. Anal. Appl., 331 (2007), pp. 396-407.

[18] --, Stabilized boundary element methods for exterior Helmholtz problems, Numer. Math., 110 (2008), pp. 145-160.

[19] I. G. Graham and I. H. Sloan, Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in , Numer. Math., 92 (2002), pp. 289-323.

[20] D. S. Jones, Integral equations for the exterior acoustic problem, Quart. J. Mech. appl. Math., 27 (1974), pp. 129-142.

[21] A. Kleefeld, Direct and Inverse Acoustic Scattering for Three-Dimensional Surfaces, Ph.D. Thesis, University of Wisconsin - Milwaukee, December 2009.

[22] R. E. Kleinman and R. Kress, On the condition number of integral equations in acoustics using modified fundamental solutions, IMA J. Appl. Math., 31 (1983), pp. 79-90.

[23] R. E. Kleinman and G. F. Roach, On modified Green functions in exterior problems for the Helmholtz equation, Proc. R. Soc. Lond. Ser. A, 383 (1982), pp. 313-332.

[24] T.-C. Lin, A proof for the Burton and Miller integral equation approach for the Helmholtz formula, J. Math. Anal. Appl., 103 (1984), pp. 565-574.

[25] --, The numerical solution of Helmholtz's equation for the exterior Dirichlet problem in three dimensions, SIAM J. Numer. Anal., 22 (1985), pp. 670-686.

[26] --, Smoothness results of single and double layer solutions of the Helmholtz equations, J. Integral Equations Appl., 1 (1988), pp. 83-121.

[27] T.-C. Lin and Y. Warnapala-Yehiya, The numerical solution of exterior Neumann problem for Helmholtz's equation via modified Green's function approach, Comput. Math. Appl., 47 (2004), pp. 593-609.

[28] Y. Liu and F. J. Rizzo, A weakly singular form of the hypersingular boundary integral equation applied to 3-D acoustic wave problems, Comput. Methods Appl. Mech. Engrg., 96 (1992), pp. 271-287.

[29] P. Martin, On the null-field equations for the exterior problems of acoustics, Quart. J. Mech. Appl. Math., 33 (1980), pp. 385-396.

[30] P. Meury, Stable Finite Element Boundary Element Galerkin Schemes for Acoustic and Electromagnetic Scattering, Ph.D. Thesis, ETH Zurich, 2007.

[31] S. Micula, Numerical Methods for the Radiosity Equation and Related Problems, Ph.D. Thesis, University of Iowa, December 1997.

[32] --, A collocation method for solving the exterior Neumann problem, Stud. Univ. Babe[section]-Bolyai Math., 48 (2003), pp. 105-113.

[33] S. G. Mikhlin, Mathematical Physics, an Advanced Course, North-Holland, Amsterdam, 1970.

[34] O. I. Panich, On the question of the solvability of the exterior boundary problem for the wave equation and Maxwell's equation, Uspekhi Mat. Nauk, 20 (1965), pp. 221-226.

[35] R. Potthast, Frechet differentiability of boundary integral operators in inverse acoustic scattering, Inverse Problems, 10 (1994), pp. 431-447.

[36] H. A. Schenck, Improved integral formulation for acoustic radiation problems, J. Acoust. Soc. Am., 44 (1968), pp. 41-58.

[37] F. Ursell, On the exterior problems of acoustics, Proc. Camb. Phil. Soc., 74 (1973), pp. 117-125.

[38] --, On the exterior problems of acoustics: II, Proc. Camb. Phil. Soc, 84 (1978), pp. 545-548.

[39] A. F. Seybert and T. K. Rengarajan, The use of CHIEF to obtain unique solutions for acoustic radiation using boundary integral equations, J. Acoust. Soc. Am., 81 (1987), pp. 1299-1306.

[40] A. F. Seybert, B. Soenarko, F. J. Rizzo, and D. J. Shippy, An advanced computational method for radiation and scattering ofacoustic waves in three dimensions, J. Acoust. Soc. Am., 77 (1985), pp. 362-368.

[41] J. J. Silva, H. Power, and L. C. Wrobel, The use ofC[degrees]>a boundary elements in an improved numerical formulation for three-dimensional acoustic radiation problems, J. Acoust. Soc. Am., 95 (1994), pp. 2388-2394.

[42] W. Tobocman, Comparison ofthet-matrix and Helmholtz integral equation methods for wave scattering calculations, J. Acoust. Soc. Am., 77 (1985), pp. 369-374.

[43] V. K. Varadan, V. V. Varadan, L. R. Dragonette, and L. Flax, Computation of rigid body scattering by prolate spheroids using thet-matrix approach , J. Acoust. Soc. Am., 71 (1982), pp. 22-25.

[44] P. C. Waterman, New formulation of acoustic scattering, J. Acoust. Soc. Am., 45 (1969), pp. 1417-1429.

[45] P. Werner, Beugungsprobleme der mathematischen Akustik, Arch. Ration. Mech. Anal., 12 (1963), pp. 155-184.

ANDREAS KLEEFELD ([dagger]) AND TZU-CHU LIN ([double dagger])

([dagger]) Brandenburgische Technische Universitat Cottbus, Konrad-Wachsmann-Allee 1, 03046 Cottbus, Germany (kleefeld@tu-cottbus.de).

([double dagger]) Department of Mathematical Sciences, University of Wisconsin-Milwaukee, PO Box 413, Milwaukee, WI 53201 (lin@uwm.edu).
```Table 3.1

Nodes and Lagrange basis functions over [sigma] for constant,

Constant

i        [q.sub.i]        [l.sub.i](s,t)

1    ([alpha], [alpha])         1

2

3

4

5

6

Linear

i           [q.sub.i]             [l.sub.i](s,t)

1      ([alpha], [alpha])      u-[alpha]/l-3[alpha]

2    ([alpha], 1 - 2[alpha])   t-[alpha]/l-3[alpha]

3    (1 - 2[alpha], [alpha])   s-[alpha]/l-3[alpha]

4

5

6

i       [q.sub.i]               [l.sub.i](s,t)

1       ([alpha],        [MATHEMATICAL EXPRESSION NOT
[alpha])            REPRODUCIBLE IN ASCII]

2     ([alpha], 1 -      [MATHEMATICAL EXPRESSION NOT
2[alpha])           REPRODUCIBLE IN ASCII]

3    (1 - 2[alpha],      [MATHEMATICAL EXPRESSION NOT
[alpha])            REPRODUCIBLE IN ASCII]

4     ([alpha], 1 -      [MATHEMATICAL EXPRESSION NOT
[alpha]/2)           REPRODUCIBLE IN ASCII]

5    (1 - [alpha]/2,     [MATHEMATICAL EXPRESSION NOT
1 - [alpha]/2)         REPRODUCIBLE IN ASCII]

6    (1 - [alpha]/2,     [MATHEMATICAL EXPRESSION NOT
[alpha])            REPRODUCIBLE IN ASCII]

Table 5.1

Accuracy for constant interpolation, n = 1024 and [n.sub.v] = 1024
with wave number K = l for a sphere.

Approximated solution
Absolute
Point                real part     imag. part      error

(10.0,11.0,12.0)    +3.1334D-02    +9.9500D-03   4.6443D-05
(5.0, 6.0, 7.0)     -2.5597D-02    -5.8477D-02   9.1254D-05
(1.0, 2.0, 3.0)     -1.4428D-01    -1.6800D-01   3.5436D-04
(1.0, 1.0, 1.0)     -2.4302D-01    +2.9767D-01   6.3241D-04
(0.0, 0.0, 2.5)     -4.1544D-01    +1.1101D-01   7.9922D-04
(1.0, 0.5, 0.5)     -1.4272D-01    +4.0560D-01   3.8040D-04
(1.0, 0.2, 0.3)     -8.9025D-02    +3.5382D-01   4.1385D-04

Table 5.2

Accuracy for linear interpolation, n = 256 and [n.sub.v] = 768 with
wave number k = 1 for a sphere.

Approximated solution
Absolute
Point                real part    imag. part      error

(10.0,11.0,12.0)    +3.1371D-02   +9.9690D-03   4.8251D-06
(5.0, 6.0, 7.0)     -2.5622D-02   -5.8556D-02   9.2227D-06
(1.0, 2.0, 3.0)     -1.4450D-01   -1.6824D-01   2.8937D-05
(1.0, 1.0, 1.0)     -2.4343D-01   +2.9806D-01   6.2783D-05
(0.0, 0.0, 2.5)     -4.1617D-01   +1.1122D-01   4.6545D-05
(1.0, 0.5, 0.5)     -1.4294D-01   +4.0581D-01   8.7465D-05
(1.0, 0.2, 0.3)     -8.9132D-02   +3.5340D-01   2.9184D-05

Table 5.3

Accuracy for quadratic interpolation, n = 256 and [n.sub.v] = 1536
with wave number k = 1 for a sphere.

Approximated solution
Absolute
Point               real part    imag. part      error

(10.0,11-0,12.0)   +3.1375D-02   +9.9712D-03   9.2773D-08
(5.0, 6.0, 7.0)    -2.5623D-02   -5.8565D-02   1.7742D-07
(1.0, 2.0, 3.0)    -1.4451D-01   -1.6827D-01   5.6770D-07
(1.0, 1.0, 1.0)    -2.4347D-01   +2.9811D-01   1.2875D-06
(0.0, 0.0, 2.5)    -4.1621D-01   +1.1121D-01   9.6095D-07
(1.0, 0.5, 0.5)    -1.4297D-01   +4.0589D-01   4.8034D-07
(1.0, 0.2, 0.3)    -8.9143D-02   +3.5342D-01   6.7632D-06

Table 5.4

Accuracy for quadratic interpolation, n = 256 and [n.sub.v] = 1536
with critical wave number k = 4.4934094579 and [eta] = k/2 for
a sphere.

Approximated solution
Absolute
Point               real part     imag. part      error

(10.0,11.0,12.0)   -1.6783D-02    -2.8273D-02   2.6769D-07
(5.0, 6.0, 7.0)    -6.3632D-02    -1.5657D-03   6.2821D-07
(1.0, 2.0, 3.0)    -8.4887D-02    -1.9717D-01   2.3636D-06
(1.0, 1.0, 1.0)    -1.9022D-02    +3.3554D-01   4.2893D-06
(0.0, 0.0, 2.5)    +1.2889D-01    -3.8034D-01   1.7094D-06
(1.0, 0.5, 0.5)    +2.7959D-01    -1.9134D-01   4.2313D-06
(1.0, 0.2, 0.3)    +7.2496D-02    -2.6136D-01   2.0580D-05

Table 5.5

Accuracy for quadratic interpolation, n = 256 and [n.sub.v] = 1536
with critical wave number k = 4.4934094579 and [eta] = 0 for
a sphere.

Approximated solution
Absolute
Point               real part    imag. part      error

(10.0,11.0,12.0)   -7.5714D-02   -9.0756D-02   8.5889D-02
(5.0, 6.0, 7.0)    -2.6487D-01   +3.7716D-02   2.0504D-01
(1.0, 2.0, 3.0)    -2.4736D-01   -9.5641D-01   7.7643D-01
(1.0, 1.0, 1.0)    -1.5694D+00   +7.4617D-01   1.6039D+00
(0.0, 0.0, 2.5)    -1.5952D+00   +4.7979D-01   1.7270D+00
(1.0, 0.5,0.5)     -7.1122D+01   +4.1024D-02   1.0177D+00
(1.0, 0.2, 0.3)    +2.7680D+00   ^t.5044D-02   2.7042D+00

Table 5.6

Accuracy for quadratic interpolation, n = 256 and [n.sub.v] with
wave number k = 7 for an ellipsoid.

Approximated solution
Absolute
Point               real part     imag. part      error

(10.0,11.0,12.0)   -7.2828D-03    +3.2065D-02   4.1804D-05
(5.0, 6.0, 7.0)    -2.4711D-02    -5.8659D-02   9.7101D-05
(1.0, 2.0, 3.0)    +9.8086D-02    +1.9089D-01   2.5669D-04
(1.0, 1.0, 1.0)    +3.1361D-01    -1.1837D-01   8.3860D-04
(0.0, 0.0, 2.5)    +1.0974D-01    -3.8551D-01   4.3519D-04
(1.0, 0.5, 0.5)    -2.4944D-01    +2.2520D-01   5.9115D-04
(1.0, 0.2, 0.3)    +7.4938D-02    +2.5711D-01   1.1569D-03

Table 5.7

Constant interpolation with [alpha] = 1/3 for a sphere.

n ([n.sub.v])   [[??].sub.4]   n ([n.sub.v])   [[??].sub.8]

4(4)                               8(8)
16 (16)                           32 (32)
64 (64)             1.25         128(128)          1.53
256 (256)           2.10         512(512)          1.73
1024 (1024)         1.51        2048 (2048)        1.68

Table 5.8

Constant, linear and quadratic interpolation for an ellipsoid.

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.1])        EOC

4(4)                7.7155D-02
16 (16)             3.2332D-02       1.25
64 (64)             3.0312D-03       3.42
256 (256)           9.7283D-04       1.64
1024 (1024)         2.2671D-04       2.10

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.1])        EOC

4(12)               8.9214D-02
16 (48)             1.0039D-02       3.15
64 (192)            1.2702D-03       2.98
256(768)           7.4601 D-05       4.09

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.1])        EOC

4(24)               1.8457D-02
16 (96)             5.7908D-03       1.67
64 (384)            1.1266D-04       5.68
256 (1536)          7.7467D-07       7.18

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(4)                5.6560D-02
16 (16)             1.9795D-02       1.51
64 (64)             1.9718D-03       3.33
256 (256)           7.0460D-04       1.48
1024 (1024)         1.6716D-04       2.08

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(12)               5.8612D-02
16 (48)             6.2687D-03       3.22
64 (192)            1.2597D-03       2.32
256(768)            7.6031D-05       4.05

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(24)               8.2448D-03
16 (96)             4.4997D-03       0.87
64 (384)            9.7819D-05       5.52
256 (1536)          1.0026D-06       6.61

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(4)                2.3263D-02
16 (16)             9.3200D-03       1.32
64 (64)             9.5112D-04       3.29
256 (256)           2.9977D-04       1.67
1024 (1024)         6.9379D-05       2.11

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(12)               2.3930D-02
16 (48)             3.0255D-03       2.98
64 (192)            4.0138D-04       2.91
256(768)            2.5589D-05       3.97

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(24)               4.6017D-03
16 (96)             1.6453D-03       1.48
64 (384)            3.5941D-05       5.52
256 (1536)          3.0691D-07       6.87

Table 5.9

Constant, linear and quadratic interpolation for a peanut.

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(4)                1.8968D-01
16 (16)             2.2830D-02       3.05
64 (64)             5.3670D-03       2.09
256 (256)           1.8429D-03       1.54
1024 (1024)         4.0915D-04       2.17

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(12)               1.5041D-01
16 (48)             1.3041D-02       3.53
64 (192)            8.6959D-04       3.91
256(768)            1.2514D-04        2.8

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(24)               2.0368D-02
16 (96)             3.9080D-03       2.38
64 (384)            2.9870D-04       3.71
256 (1536)          7.4099D-06       5.33

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(4)                1.2105D-01
16 (16)             1.4970D-02       3.02
64 (64)             4.1890D-03       1.84
256 (256)           1.2316D-03       1.77
1024 (1024)         2.7237D-04       2.18

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(12)               1.0084D-01
16 (48)             1.1493D-02       3.13
64 (192)            7.0722D-04       4.02
256(768)            8.8471D-05       3.00

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(24)               1.7600D-02
16 (96)             3.2544D-03       2.44
64 (384)            2.3553D-04       3.79
256 (1536)          4.7706D-06       5.63

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(4)                5.0542D-02
16 (16)             6.4053D-03       2.98
64 (64)             1.5020D-03       2.09
256 (256)           4.9754D-04       1.59
1024 (1024)         1.1026D-04       2.17

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(12)               4.1227D-02
16 (48)             4.2706D-03       3.27
64 (192)            2.4639D-04       4.12
256(768)            3.6109D-05       2.77

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(24)               6.1177D-03
16 (96)             1.2040D-03       2.35
64 (384)            8.6645D-05       3.80
256 (1536)          2.0222D-06       5.42

Table 5.10

Constant, linear and quadratic interpolation for an acorn.

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(4)                2.5980D-01
16 (16)             2.2626D-01       0.20
64 (64)             5.4086D-02       2.06
256 (256)           6.5394D-03       3.05
1024 (1024)         1.2956D-03       2.34

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(12)               6.8475D-02
16 (48)             5.7888D-02       0.24
64 (192)            2.9618D-03       4.29
256 (768)           2.3629D-04       3.65

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

4(24)               1.8245D-01
16 (96)             1.1798D-02       3.95
64 (384)            1.5205D-03       2.96
256 (1536)          1.7514D-05       6.44

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(4)                2.8900D-01
16 (16)             1.6782D-01       0.78
64 (64)             4.2960D-02       1.97
256 (256)           5.1552D-03       3.06
1024 (1024)         1.0109D-03       2.35

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(12)               5.2062D-02
16 (48)             3.9641D-02       0.39
64 (192)            2.1820D-03       4.18
256 (768)           1.3344D-04       4.03

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

4(24)               1.2662D-01
16 (96)             9.2453D-03       3.78
64 (384)            1.1285D-03       3.03
256 (1536)          1.0973D-05       6.68

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(4)                8.8502D-02
16 (16)             7.6075D-02       0.22
64 (64)             1.6292D-02       2.22
256 (256)           1.9068D-03       3.09
1024 (1024)         3.7004D-04       2.37

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(12)               1.9049D-02
16 (48)             1.5944D-02       0.26
64 (192)            8.9133D-04       4.16
256 (768)           8.0285D-05       3.47

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

4(24)               4.9425D-02
16 (96)             3.3440D-03       3.89
64 (384)            4.3950D-04       2.93
256 (1536)          5.1072D-06       6.43

Table 5.11

Constant, linear and quadratic interpolation for the unit cube.

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

12 (12)             2.1878D-01
48 (48)             8.0459D-02       1.44
192 (192)           1.8311D-02       2.14
768 (768)           4.4321D-03       2.05
3072 (3072)         1.2767D-03        1.8

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

12 (36)             1.1097D-01
48(144)             4.6919D-03       4.56
192 (576)           7.3286D-04       2.69
768 (2304)          5.0750D-04       0.53

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.l])        EOC

12 (72)             2.1166D-02
48(288)             2.0360D-03       3.38
192(1152)           5.9166D-04       1.78

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

12 (12)             1.5350D-01
48 (48)             3.0520D-02       2.33
192 (192)           8.9033D-03       1.78
768 (768)           2.3155D-03       1.94
3072 (3072)         6.9779D-04       1.73

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

12 (36)             4.4018D-02
48(144)             2.8649D-03       3.94
192 (576)           4.4041D-04       2.70
768 (2304)          2.8441D-04       0.63

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.2])        EOC

12 (72)             1.0229D-02
48(288)             9.2309D-04       3.47
192(1152)           3.2087D-04       1.52

Constant interpolation with [alpha] = 1/3

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

12 (12)             8.2077D-02
48 (48)             9.3621D-03       3.13
192 (192)           3.5481D-03       1.40
768 (768)           9.4360D-04       1.91
3072 (3072)         2.9047D-04       1.70

Linear interpolation with [alpha] = 1/6

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

12 (36)             1.7682D-02       3.32
48(144)             1.7645D-03       3.08
192 (576)           2.0926D-04       0.72
768 (2304)          1.2697D-04

Quadratic interpolation with [alpha] = 1/10

[[epsilon].sub.n]
n ([n.sub.v])      ([P.sub.3])        EOC

12 (72)             5.4511D-03
48(288)             4.0937D-04       3.74
192(1152)           1.4189D-04       1.53
```
COPYRIGHT 2012 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.