A New Spectral Boundary Integral Collocation Method for Three-Dimensional Potential Problems.
In this work we propose and analyze a new global collocation method for classical second-kind boundary integral equations of potential theory on smooth simple closed surfaces $\Gamma \subset {\Bbb R}^3$.  Under the assumption that $\Gamma$ is diffeomorphic to the unit sphere $\partial B$, the original equation is transferred to an equivalent one on $\partial B$ which is solved using collocation onto a nonstandard set of basis functions. The collocation points are situated on lines of constant latitude and longitude.  The interpolation operator used in the collocation method is equivalent to a certain  discrete orthogonal (pseudospectral) projection, and this equivalence allows us to establish the fundamental properties of the interpolation process and subsequently to prove that our collocation method is stable and super-algebraically convergent.  In addition, we describe a fast method for computing the weakly singular collocation integrals and present some numerical experiments illustrating the use of the method. These show that at least for model problems the method attains an exponential rate of convergence and exhibits a good accuracy for very small numbers of degrees of freedom.
Introduction
In this work we present a new spectrally convergent collocation method for second-kind
boundary integral equations of potential theory. Our method is the result of a search
for a three-dimensional analogue of well-known discrete global Galerkin or "pseudospec-
tral" schemes for integral equations on planar contours (see, e.g., [15], [19], [4] or [14]).
In these two-dimensional methods the contour is parametrized by a periodic univariate
function and the basis functions used to approximate the solution are constructed from
trigonometric polynomials with respect to this parameter. With an appropriate choice of
quadrature rule the discrete orthogonal projection is then also an interpolation operator,
the discrete Galerkin method is equivalent to a collocation scheme and hence is relatively
cheap to implement. Stability and spectral orders of convergence are then proved using
reasonably standard arguments.
It is interesting to consider the extension of such ideas to three dimensions. In that
case a spectral method may be quite attractive since it has the potential to achieve an
acceptable accuracy with a significantly smaller linear algebra overhead than that required
by conventional (piecewise) methods. In line with spectral methods in general, we consider
here only problems with smooth solutions and hence (in common with the two-dimensional
case) we assume that our boundary integral equation is posed on a smooth surface.
Although most boundaries of interest in structural or mechanical engineering are not
smooth, there are many applications of smooth boundaries in other branches of science
- such as biology and medicine, where problems with corners and edges are rare.
We make the further assumption that our boundary is diffeomorphic to the unit sphere,
which means that the sphere may be used as a parametrization of the boundary. Although
this does restrict the range of problems which can be treated, it is in some sense a
natural generalisation of the fact that in two dimensions all simple closed contours may
be parametrized by the unit circle. Moreover there exists quite a lot of interest already
in the literature on problems of this form.
Applications of boundary value problems on spherical or nearly spherical geometries
are found for example in global weather prediction and in geodesy see, e.g. [21], [11].
Numerical analyses of various methods for such problems and further applications can be
found in [2], [1], [5], [7, x3.6], [8], [9], [12], [16], [23], [24], [25] and the associated references.
For a recent review of spectral methods in general see [10].
For boundary integral equations of the type considered in this paper, nondiscrete (and
therefore not fully practical) Galerkin methods using spherical harmonic basis functions
are the only global approximation schemes whose convergence has been proved ([2], [1],
[16]).
Two questions then naturally arise:
(a) whether these (or related) Galerkin methods still converge when quadrature rules
are used to compute the Galerkin integrals, and
(b) whether any such discrete Galerkin method is equivalent to a collocation scheme.
A discrete method using spherical harmonics and closely related to those used for implementation
(but not covered by the proofs) in [2], [3], [1] and [16] is analysed completely
in [12] for boundaries which are exactly spherical.
However the convergence proof used there fails when the boundary is merely diffeomorphic
to the sphere. This is essentially because the corresponding discrete Galerkin
projection has a bound which grows so fast that the usual stability proof fails. Thus the
answer to (a) is still unknown in this case. With regard to (b), it is known from [22] that
the answer is "no" when the underlying basis functions are spherical harmonics.
The very interesting unpublished Ph.D. thesis of Wienert [25] (see also [7]) suggested
some fully discrete variants of Galerkin's method with both spherical harmonics and other
basis sets and proved some new approximation theory (yielding exponential convergence
rates) for the discrete integral operators. However Wienert did not give a convergence
analysis for the corresponding approximate integral equations and thus does not provide
an answer to (a). Nevertheless one of his methods suggests the use of a discrete orthogonal
projection, with basis functions chosen as Fourier modes on R 2 (mapped to the unit sphere
with polar coordinates). This happens also to be an interpolatory projection and hence
satisfies (b).
Our method takes Wienert's idea as starting point. We construct a basis set by
choosing a distinguished set of Fourier modes on the plane and mapping them to the unit
sphere, again with spherical polar coordinates. However unlike Wienert, we only choose
those Fourier modes which map to continuous functions on the sphere and which have
a certain reflectional symmetry. (Similar ideas are in [21].) Our choice ensures that the
spherical polar coordinate transformation forms a bijection between our basis set and the
distinguished set of Fourier modes and that the usual discrete orthogonal projection onto
these Fourier modes corresponds to a discrete orthogonal projection on the sphere. Both
are interpolatory and have norm growing with O(log 2 (N)), when the number of degrees
of freedom is O(N 2 ). This growth is slow enough to allow stability of the corresponding
discrete Galerkin method to be proved. So the answer to both (a) and (b) is "yes" for
this new method. These results and the implementation of the method are the principal
aim of this paper.
The layout of the paper is as follows. In Section 2 we describe our basis set on
the sphere and prove the required properties of the corresponding discrete orthogonal
projection. The space spanned by our basis functions has dimension 2.
It contains, as a proper subset, spherical polynomials (which are smooth on the sphere)
as well as extra functions which are continuous on the sphere but lack smoothness at
the poles. Once Section 2 is done the convergence of the collocation method is obtained
very easily in Section 3. The method is proved to be super-algebraically convergent. It
should be possible also to show exponential convergence in an appropriate analytic setting
(similar to that in [25]), but we do not attempt that here.
Computation of the collocation integrals is considered in Section 4. Here were are
again influenced by [25]. If the basis functions were smooth on the sphere there is a very
nice way of evaluating the required weakly singular integrals using a rotational change
of variables and then polar coordinates, which eliminates the weak singularity in the
integrand (see also [7, x3.6]). In Section 4 we consider this transformation applied to
our collocation integrals and we show that it yields integrands on polar coordinate space
('; OE) 2 [0; -] \Theta [0; 2-] with a mild singularity in the second derivative at certain points.
Fortunately this singularity seems not to affect the empirical convergence of standard
quadrature rules for the collocation integrals. These perform as if the kernel were
smooth, as is seen in Section 5 where we use the tensor product Gauss-Legendre rule for
the integrals (after having transformed the singularities to the edges of the domain of in-
tegration) and observe exponential convergence. This would be difficult to prove because
of the mild singularity in the integrands. We also observe exponential convergence of the
solutions to the integral equation and to the corresponding solution to the underlying
harmonic boundary value problem, provided the number of quadrature points increases
appropriately with the number of degrees of freedom. For regular problems the method is
extremely accurate: relative errors of about 10 \Gamma3 are obtained by inverting linear systems
of dimension about 40. Because of the relatively small linear algebra requirements of this
method we are able to code it using MATLAB ([18]). On modern machines this has been
found to yield accurate solutions of potential problems quickly enough to be usable as an
interactive tool.
We finish this section by describing in more detail the problem to be solved. Let
\Omega ae R 3 be a bounded domain with simple smooth closed boundary \Gamma. For any j 2 \Gamma let
denote the outward unit normal to \Gamma and let d\Sigma(j) denote surface measure on \Gamma.
denote the Euclidean norm on R 3 . We are primarily concerned in this paper with
fast (spectrally accurate) solutions of the boundary integral equation
Z
\Gamma2-
This is the classical second kind equation which arises when the indirect boundary integral
method with a double layer potential ansatz is used to solve the Dirichlet problem for
Laplace's equation on \Omega\Gamma The method of this paper can also be used to solve the Neumann
problem
on\Omega or the corresponding (exterior) problems
Our method for solving (1.1) is based upon expanding the solution u in terms of a
certain set of global basis functions, with the coefficients in the expansion determined by
collocation. We make the following assumptions which are fundamental to the construction
of this collocation method.
There exists a smooth map q : D ae R 3 (where D is a domain containing the
unit sphere
bijective and its Jacobian dq is assumed to satisfy
det dq(x) 6= 0; x
Since @B ae D is compact it follows from (A.1) that:
All derivatives of q are bounded on @B : (1.2)
With the change of variable the directed measures on \Gamma and @B are
related by Nanson's formula ([20, p.88]):
where doe(y) denotes surface measure on @B at y 2 @B, y is the unit normal to
@B at y and Adj(M) denotes the adjugate of the 3 \Theta 3 matrix M (i.e. the transposed
matrix of its cofactors).
Now, with the substitution using (1.3), we can transplant
onto @B to obtain
Z
where
We will be concerned with the solution of (1.4), which is equivalent to (1.1). Thus we
rewrite u ffi q as u and f ffi q as f and we abbreviate (1.4) as
where K is the integral operator on @B induced by the kernel k. It is well known that
(1.6) has a unique solution u which is smooth when f is smooth.
series and interpolation on the
unit sphere
2.1 Functions on the unit sphere
Let C r (@B); r - 0, denote the usual space of r-times continuously differentiable functions
on the unit sphere @B with norm jj:jj 1;r;@B . Let
denote the space of r-times continuously differentiable functions on R 2 which are 2-
periodic in each argument with norm jj:jj 1;r;D . In addition, for r
we let C r;ff (@B) denote the functions in C r (@B) whose rth partial derivatives (defined with
respect to a suitable chart) satisfy a H-older condition of order ff. (C r;0 (@B) j C r (@B).)
Also, let C r;ff (D) denote the analogous space on the Euclidean domain D. Let k:k 1;r;ff;@B
and k:k 1;r;ff;D denote the corresponding norms. When we drop the index r in the
above notations.
The standard mapping from R 2 to @B is the spherical polar coordinate transformation
This is invertible when considered as a map from (0; -) \Theta [0; 2-) onto the punctured sphere
@Bnfn; sg, where
(for any OE 2 R) are the north and south poles, respectively. On @Bnfn; sg, p has a
continuous inverse given by p is the solution
of the equations
The function p is smooth and 2-periodic in ('; OE) and any v 2 C(@B) thus induces a
function J v 2 C(D) given by
Analogously, for all r
C r;ff (D) and
It is clear that J is an isometry from C(@B) into C(D), i.e.
Since p has reflectional symmetry:
and is independent of OE 2 R at it is natural to introduce the
proper subspaces of C(D):
are independent of OE 2 Rg:
Then it is easy to show that J : C(@B) ! S(D) is an isometric isomorphism with inverse
given by J
w(0;
(w
w(-; s:
We shall use (2.6) to define our approximation procedure in C(@B). We do this by
first recalling the discrete Fourier projection on C(D). Provided this is carefully defined
it turns out also to be a projection operator in S(D). Using the isometry (2.6) this gives
us an analogue of discrete Fourier approximation in C(@B). The details are explored in
the following subsections.
2.2 Discrete Fourier series in C(D)
The Fourier modes
are orthonormal with respect to the usual inner product
We shall use a corresponding discrete inner product defined using the quadrature points
Throughout the paper we assume that N is odd, so that -=2 is never a quadrature point.
The commonly used discrete analogue of (2.8) is then obtained by applying the tensor
product trapezoidal rule based at the nodes (2.9) to obtain
for Throughout the paper
means the usual summation with only half
the first and last terms included and C denotes a generic constant which is independent
of N .
The truncated discrete Fourier series of a function w 2 C(D) is then defined by
m=\GammaN
The following lemma discusses various properties of EN . Before we prove it we recall the
following useful formula2N
Lemma 2.1
(iii)
log 2 N
Proof Let EN denote the set defined on the right hand side of (2.12) and let w 2 C(D).
Note first that
Hence, for m;
\GammaN thus it follows that hw; e m
\GammaN
Hence from the
definition of EN , we get
Moreover, using the identity
he m
we obtain, for
together with
EN (e N
Hence a short calculation shows that
Thus (i) and (ii) follow from (2.15) and (2.16). Also observe that for
m=\GammaN
Then (iii) follows on using (2.12). For (iv), let w 2 C(D) and ('; OE) 2 D. Then
where D
N is the modified Dirichlet kernel [26, Vol I, p.50] given by
Hence
sup
We show in the Appendix that sup
proof is obtained by modifying the argument in [26, Vol.II, p.19]). Hence (iv) follows from
(2.19).
Finally, to obtain (v), let w 2 C r;ff (D); r - 0; ff 2 (0; 1]: Then recall that for each
there exists a trigonometric polynomial T n ('; OE) of degree n in each variable such
that
[17, p.89, Theorem 7]. Since EN TN
Using (iv) and (2.20) in (2.21), we get (v).
2.3 Discrete Fourier series in R(D) and C(D)
The next result shows that the discrete orthogonal projection EN maps R(D) into itself
with range a subspace of EN (C(D)). To show this, it is convenient to introduce functions
\Gamman
Observe that r m
Lemma 2.2
m=\GammaN
and hence EN
Proof Let RN denote the set defined by the right hand side of (2.23). We prove first
that EN (R(D)) ' RN , using essentially the method of [25, p.26] but in a slightly different
setting.
Let w 2 R(D). We write
where
m=\GammaN
and
w (2)
m=\GammaN
m=\GammaN
Using the fact that w 2 R(D), and the 2- periodicity of w and e m
n we have
w (2)
m=\GammaN
m=\GammaN
Using this together with (2.24) and (2.25), we have
m=\GammaN
n , this yields
m=\GammaNN
where
By easily that r N
Z. Hence we have
. Also for 0 - jmj - N with m odd we have
which implies ff m
Thus we have shown
The fact that RN ' EN (R(D)) follows from the following four identities, which can
be proved using standard arguments.
EN r q
EN (r \GammaN
EN r q
r q
EN (r \GammaN
r \GammaN

The next result shows that in fact EN leaves S(D) invariant.
Lemma 2.3 For all w 2 S(D), and all OE 2 R, we have
Proof Let w 2 S(D) and set w(0;
(2.12) we have, for all OE 2 R,
m=\GammaN
m=\GammaN
m=\GammaN
Similarly Lemma 2.2, we have
The previous lemmas show that EN is an interpolatory projection operator on S(D). Its
image EN (S(D)) constitutes a distinguished subspace of the space spanned by the Fourier
modes. It is from this subspace that a basis of interpolating trigonometric polynomials in
S(D) can be constructed. Using J \Gamma1 given by (2.6) this leads directly to an interpolatory
basis in C(@B).
To define the basis in S(D), define a function
\Theta R! R by
f0g. Moreover the sets fs m
are orthogonal with respect to the inner product
(2.8). The spaces
are mutually orthogonal and have dimension
Hence the space
is a subspace of S(D) with dimension 2N 2. We shall show in Theorem 2.6
below that for each w 2 S(D), EN w is the unique element of   N which interpolates w at
a certain set of 2N points. Thus EN is an interpolation operator. The proof is
obtained with the help of the following two lemmas.
Lemma 2.4
where
Proof Observe first that, since the basis functions for   N are independent of those for
f0g. Now, by examining each of the basis functions for   N \Phi \Delta N
in turn and using the characterisation (2.13) and the definition of R(D) we can show that
Hence if wN 2   N \Phi \Delta N then by Lemma 2.1 (ii),
Conversely, suppose wN 2 EN (R(D)). Then from Lemma 2.2 there exist constants ff m
with ff \GammaN
n and ff m
when jmj is odd such that
m=\GammaN
Thus, recalling (2.22),
wN ('; OE) =-
m=\GammaN
even
m=\GammaN
sin n' exp imOE:
where
and
w (2)
Since w (1)
it remains to show
w (2)
To obtain (2.35) let m 2 f0; :::; Ng be even. Then observe that for n 2 f1; :::; Ng with n
odd,
cos n' cos
sin j' sin ' cos mOE
In addition if n 2 f0; :::; Ng with n even then
cos n' cos
sin j' sin ' cos mOE
even and n 2 f0; :::; Ng
cos n' sin mOE 2   N \Phi
Now (2.36)-(2.38) imply (2.35) and wN 2   N \Phi \Delta N . Hence
and, together with (2.34) this proves the result.
Lemma 2.5 EN
Proof Let wN 2   N . It is easily checked that wN 2 S(D). Since by (2.33), wN 2
EN (C(D)), Lemma 2.1 (ii) yields wN 2 EN (S(D)) and so
On the other hand suppose wN 2 EN (S(D)). Then also wN 2 EN (R(D)) and by
Lemma 2.4 there exist unique w (1)
N . Also,
by Lemma 2.3, wN 2 S(D) and, since   N ' S(D), we have w (2)
Hence there exist C
w (2)
for all OE 2 R. From the definition of \Delta N it follows easily that w (2)
Hence EN (S(D)) '   N and the result follows from this and (2.39).
These lemmas then lead to the following theorem describing the interpolatory properties
of EN .
Theorem 2.6 For all w 2 S(D), EN w is the unique element of   N which satisfies
Proof Let w 2 S(D). From Lemmas 2.1 and 2.5, EN w 2   N has the interpolation
property (2.40) -(2.42). To establish uniqueness, suppose wN 2   N satisfies wN (0;
Then since   N ae S(D) we have in fact wN thus
we have (using Lemma 2.1 (ii)),
as required.
2.4 A discrete orthogonal projection on C(@B)
Finally in this section we study the interpolatory operator on C(@B) induced by the
operator EN on S(D). That is we set
is the isometric isomorphism introduced in (2.3). We introduce
the interpolation points on the sphere
Together with the north pole n and the south pole s these form a set of 2N
points on @B. Similarly we introduce the interpolation space
Theorem 2.7
(ii) For v 2 C(@B); EN v is the unique element of SN with the property
(iv) For
log 2 N
Proof This is immediate from Lemma 2.1, Theorem 2.6 and the properties of J .
3 The Collocation (Pseudospectral) Method
In this section we introduce our pseudospectral method for (1.6) and prove its convergence.
First we give an algorithmic description of the method. Let
the dimension of the space introduced in (2.30) and (2.43). A basis for this
space is
For convenience we denote this basis f/ dg. It is natural to choose the
ordering implicit in (3.1), i.e. to choose the first 2N basis functions 2N to be
are
and so on with / being the last two basis functions J
N+1 .
The collocation points are the set
These have a corresponding natural ordering x 1k
and so on, with n and s being the last two collocation points.
Our method for (1.6) then consists of seeking a numerical solution
d
a
where the fa p g are scalars chosen so that
or equivalently
d
a
We discuss the calculation of the matrix entries
in the next section.
The proof of convergence of (3.3) is straightforward from the results of x2.
Theorem 3.1 Let f 2 C(@B). For N sufficiently large (3.3) has a unique solution
. If, in addition, f 2 C r (@B) with r - 1 then, for N sufficiently large,
log 2 N
Proof Observe that by Theorem 2.7, uN satisfies (3.3) if and only if
It is well known e.g. [13, x6.4] that exists and is bounded on C(@B). Moreover
[6, is bounded for ff 2 (0; 1). Hence, for all u 2 C(@B),
Theorem 2.7 (iv) yields
log 2 N
log 2 N
for any ff 2 (0; 1). Hence perturbation theory shows
that exists for N sufficiently large and is uniformly bounded. Then, since
then since I \Gamma K is also invertible on C r (@B); u 2 C r (@B) and the result
follows from Theorem 2.7 (iv).
4 Collocation Integrals
To implement the collocation method we must calculate integrals of the form
Z
where v is one of the basis functions for SN introduced in Section 2 and k is given by
(1.5).
To compute the weakly singular integrals (4.1) we follow the approach of Wienert
([25]) and introduce the Householder matrix defined for x 2 @B by
where
and
and
Moreover, since H(x) has eigenvalues \Gamma1 and 1 in directions u(x) and fu(x)g ? , it follows
that
geometrically, multiplication by H(x) represents reflection
in the plane perpendicular to the line joining x to (0; 0; - (x)) T , passing through the mid-point
of that line. Since surface measure on @B is unaffected by an orthogonal change of
variable, we may rewrite (4.1) as
Z
By virtue of (4.5) the integrand in (4.7) now has a singularity at We
handle this by taking polar coordinates to obtain
Z -Z 2--(x; '; OE)v(H(x)p('; OE))dOEd' (4.8)
with
We shall approximate the integrals (4.8) by quadrature. As we shall see in Lemmas 4.1-
4.3 below, -(x; '; OE) is smooth on ('; OE) 2 D with all its derivatives uniformly bounded
in x 2 @B. This is obtained using the arguments of Wienert [25] (see also [7]) but in the
context of C 1 rather than analytic function spaces.
We shall also examine below the behaviour of v(H(x)p('; OE)). As a function of ('; OE) 2
D this may have a very mild singularity when p('; We characterise this behaviour
in results from Lemma 4.4 below onwards. Together these results show that the integrand
in (4.8) is smooth except for a mild singularity at a single point. We devise a suitable
method for approximating these integrals in the context of a model example in Section 5.
When examining the properties of - it is useful to define the function
We can then write
with (for x 3 -
sin '
and
then we replace 0 by - in (4.11) and (4.12).)
Before we analyse - 1 and - 2 we point out a useful simpler expression for the last factor
on the right-hand side of (4.12) This is obtained by recalling the formula (valid for any
n \Theta n matrix A and any two vectors v; w 2 R n
Then (as is easily verified),
and using (4.4), (4.6) we have
Hence (4.4) and (4.13)-(4.15) yield
sin '
('; OE)
Using (4.13) again together with the chain rule, we have
sin '
('; OE)
To examine the smoothness of the first factor in the integrand of (4.8) it is convenient
to make the following definition.
Definition Given v : @B \Theta D ! R, we say v(x; '; OE) is C 1 (D) uniformly on @B if the
function ('; OE) 7! v(x; '; OE) is in C 1 (D) and all its derivatives are uniformly bounded in
Our first lemma is the following.
Lemma 4.1 Q(x; '; OE) is C 1 (D) uniformly on @B.
Proof From (1.2) it is clear that the function ('; OE) 7! Q(x; '; OE) is in C 1 (D). By the
chain rule, its ('; OE) derivatives are products of the Euclidean derivatives of q evaluated
at H(x)p('; OE) 2 @B and ('; OE) derivatives of H(x)p('; OE). The former are bounded
independently of x (by (1.2)). For the latter, observe that for all x 2 @B, we have, using
(4.4),
d k+l
d k+l p
('; OE)
d k+l p
('; OE)
and the results follow.
The next two lemmas treat the two factors on the right-hand side of (4.10).
Lemma 4.2 - 1 (x; '; OE) is C 1 (D) uniformly on @B.
Proof Let x 2 @B and without loss of generality assume x 3 - 0. Then using (4.11),
we have
The result then follows from the following properties of -
Q which we shall establish
below
uniformly in x ; (4.18)
and proceed analogously.)
To obtain (4.18), observe that
Z 1d
dt
Z 1@Q
Hence
Z 1@Q
Z 1@Q
and (4.18) follows from Lemma 4.1.
To obtain (4.19) observe that -
Since D is compact we simply have to show that
lim
However by (4.21), -
Q(x; 0; OE) is well-defined and non-negative. Thus, if (4.22) is not true,
then (4.21) implies
for some OE 2 [0; 2-]. By (4.9) and (4.5) this implies
cos OE
sin OE3
cos OE
sin OE3
Hence dq(x) has a zero eigenvalue. However this is ruled out by assumption (A.2) of the
Introduction and so (4.19) (and hence the result) follows.
Lemma 4.3 - 2 (x; '; OE) is C 1 (D) uniformly on @B.
Proof As in Lemma 4.2, without loss of generality, let x 2 @B with x 3 - 0. By (4.20),
we have
Z 1(
Z 1Z 1d
ds
ds
ds
Thus, in view of (4.16),
ds
\Theta
Hence dividing (4.12) top and bottom by ' 2 we have a numerator which is C 1 (D) uniformly
on @B and a denominator which is precisely -
complete the proof.
Our next set of results concerns the second factor in the integrand of (4.8), namely
the function
where v is one of the basis functions for the space SN defined in (2.43), i.e.
where w is one of the basis functions for   N listed in (2.29). It turns out to be convenient
to study the function v(x; '; OE) as ('; OE) ranges over all of R 2 . Then, recalling (2.6)
motivates us to introduce the set
('; OE) 2 P(x) if and only if p('; OE) 2 f\Gammax; xg. If, in particular,
then it is easy to see that
Then by (2.43) and (2.6), we have, for ('; OE) 2 R 2 nP(x),
Moreover if ('; OE) 2 P(x), then either H(x)p('; or s in which case
It follows easily that ('; OE) 7! v(x; '; OE) defines a function in S(D) for each fixed x 2 @B.
It is convenient to write (4.27) as
where
The first step in characterising the smoothness of v is to examine the properties of
G. Here some care must be taken, since (as well as being undefined at n and s), the
second component of p \Gamma1 suffers a jump discontinuity of amplitude 2- along the line
(0; -)g on @B, This motivates the introduction of the set
Then G(x; '; OE) is well-defined for ('; OE) 2 R 2 nZ(x), where
and in fact is smooth there, as is shown in the following lemma.
Lemma 4.4 Let x 2 @B. Then the function ('; OE) 7! G(x; '; OE) is in (C 1 (R 2 nZ(x))) 2 .
Proof If ('; OE) a solution of the parameter-dependent
problem:
F is infinitely continuously differentiable in y and ('; OE). The Jacobian of F with respect
to y 2 R 2 is the 3 \Theta 2 matrix
which may easily be shown to have rank 2 when
evaluated at G(x; '; OE) for all ('; OE) 2 R 2 nZ(x). Hence the result follows from the implicit
function theorem.
Our next result obtains formulae for the first partial derivatives of G.
Notation If a; b are vectors in R n then [a; b] denotes the n \Theta 2 matrix which has a as
its first column and b as its second column.
Lemma 4.5 Let x 2 @B. Then, for ('; OE) 2 R 2 nZ(x), we have
sin
@' sin g 1
cos
sin
('; OE)
where the components of G and their derivatives are evaluated at (x; '; OE). Moreover
extend as continuous functions to ('; OE) 2 R 2 nP(x).
Proof Since p(G(x; '; we can differentiate this with respect to '
first and then OE to obtainB @
cos
cos
sin
('; OE)
Then, multiplying both sides of (4.32) by the matrix
cos
\Gammasin
yields (4.31). Since the right-hand side of (4.31) is continuous across L(x) the result
follows.
Our next result concerns the smoothness of the function v(x; '; OE) given by (4.29).
We use the following standard multi-index notation : Given
denotes the partial derivative @ jffj /=@' ff 1 @OE ff 2 where
Theorem 4.6 (i) For fixed x 2 @B the function ('; OE) 7! v(x; '; OE) is infinitely continuously
differentiable on R 2 nZ(x).
(ii) If x 2 fn; sg all the partial derivatives of v(x; \Delta) on R 2 nZ(x) have a continuous
extension to R 2 nP(x) and are bounded on R 2 nP(x).
(iii) If x 2 @Bnfn; sg then for each ff with jffj - 1, D ff v(x; \Delta) has a continuous
extension to R 2 nP(x) and there exists a function Rx;ff 2 C 1 (R 4 ) which is 2-periodic in
each argument such that
where v; g 1 and g 2 are evaluated at (x; '; OE).
Proof By (4.29), Lemma 4.4 and the fact that w 2 C 1 (D), part (i) follows directly.
To obtain (ii), observe that if x 2 fn; sg then, for ('; OE) 2 (0; -) \Theta [0; 2-), H(x)p(';
Hence, we have v(x; '; which is
smooth and 2-periodic on (0; -) \Theta R. Also, since v(x; \Delta) 2 S(D), if ('; OE) 2 (\Gamma-; 0) \Theta R,
then
Then v(x; '; OE) is defined by extending these formulae 2-periodically and part (ii) follows.
Part (iii) is proved by induction on first examine the case
simple analysis of (2.29) shows that any basis function w 2   N satisfies either
or
where functions which are 2-periodic in each argument. We shall consider
only the case (4.34). The case (4.35) is similar but slightly simpler. Substituting (4.34)
into (4.29) and differentiating with respect to ' gives
cos
where v; g 1 and g 2 are evaluated at (x; '; OE). Now substituting the expressions for @g 1
and
sin
found in (4.31) yields an expression of the form (4.33) in the case
identical argument can be used to obtain an analogous expression for @v
. This completes
the proof of (4.33) when
Now suppose (4.33) holds for all
By the inductive hypothesis we have
Differentiating with respect to ' yields
sin
sin
sin
with Rx;~ ff evaluated at (g 1 substituting the expressions for @g 1 =@' and
sin found in (4.31) yields an expression for D -
ff v of the form (4.33).
If -
then we set ~
proceed analogously with ' replaced by OE.
Hence (4.33) holds when and the assertion follows.
Theorem 4.6 (iii) allows the higher derivatives of v(x; \Delta) to blow up as ('; OE) approaches
any point in P(x). In the following corollary we give more detail of the blow-up behaviour
which may arise. Recall that P(x) is characterised by (4.26) when x satisfies (4.25).
Corollary 4.7 Let x 2 @Bnfn; sg. For any multi-index ff, D ff v(x; '; OE) has a continuous
extension to ('; OE) 2 R 2 nP(x). Moreover there exists a constant M ff such that, if
for all ('; OE) 2 R 2 sufficiently close (but not equal) to ( -
OE). (Here k \Delta k denotes the
Euclidean norm in R 2 .)
Proof Since Rx;ff defined in Theorem 4.6(iii) is 2- periodic in its second argument and
2 has a jump discontinuity of amplitude 2- along the curves in L(x) the first assertion
follows directly from (4.33). This yields the estimate
for ('; OE) sufficiently close to ( -
OE), where
To obtain the second assertion from the bound (4.37), we observe that when ('; OE) is
near to but not equal to ( -
OE) then (';
2 P(x) and g 1 2 (0; -) is uniquely defined by
cos
Suppose x is given by (4.25) and recall (4.26). It is then a simple calculation to show
that
A Taylor expansion about ( -
for some constant C and all ('; OE) sufficiently close to ( -
OE). Hence from (4.38)
sin
which, together with (4.37), yields (4.36).
Remark 4.8 It is not difficult to find an illustration of the sharpness of the estimates
in Corollary 4.7. For example
OE) with ( ~
sin ' sin OE
sin ' cos
When examining the smoothness of G and v we can consider the limit ('; OE) ! ( ~ '; ~
along infinitely many paths. Two such paths are
Then a simple calculation shows that G(x; ';
(ffi; -=2) on Path 2. Making use of this and (4.31) shows that as ('; OE) ! ( ~
OE) on Path 1,
OE) on Path 2 then
Now let us consider, for example, the basis function
Using (4.39), (4.40), it is then straightforward to check that v(x; ';
OE) on Path 1
whereas
OE) on Path 2:
So @v=@' is discontinuous at (';
OE). In this case the limiting values of @ 2 v=@' 2 are
actually equal along each of Paths 1 and 2. However if we consider instead v(x; ';
then @v
has the same limiting value along Paths 1 and 2 but an elementary calculation
shows
OE) along Path 1
and
OE) along Path 2:
So the estimate O((sin g 1 ) 1\Gammajffj ) of Corollary 4.7 cannot, in general, be strengthened.

5 Numerical Experiments
As a model case we solve (1.1) where \Gamma is the ellipsoid
and\Omega is the interior of \Gamma. This was converted to (1.6) using the mapping
We used the method in x3 to solve (1.6) yielding the solution
uN on @B which is identified with a function on \Gamma (also called uN ), defined using the
inverse of q. Having found this we can compute the potential in \Omega\Gamma
Z
which approximates the solution U of the boundary value problem
on \Gamma. We consider two cases
2.
In Case A, (1.1) has the exact solution u and the solution to the boundary value
problem is U j 1. Because of the quasioptimal estimate (3.5) and since the constant
functions lie in the basis space for the collocation method we have uN j 1 for all N
provided the collocation integrals are computed exactly. Thus Case A can be used as
a test of the accuracy of the quadrature method which we shall use for the collocation
integrals.
In Case B the solution of (1.1) is not known analytically but
which can be used for comparison with the computed value UN (-). In all the
tables of results given below kuN \Gamma 1k1 is computed by finding the maximum of
at the collocation points on @B, and -  is the (randomly chosen)
point x  on the unit sphere. For any table of results containing
a sequence (a N ) which tends to 0 as N ! 1, exponential convergence is tested by
conjecturing that
and then computing ff from results for two different values of N by
oe
log(a
It remains to explain how we use the information derived in x4 to devise an accurate
method of computing the collocation integrals. Since the collocation method has a spectral
convergence rate we do not want to use any quadrature rule which would destroy that
fast rate and so this question is vitally important.
Recall that the collocation integrals are given by (4.8). The integrands have singularities
(described by Corollary 4.7) if x 6= n or s. So suppose
OE) with ~
and ~
Then we are concerned with singularities in the integrand of (4.8) at any
point in the region [0; -] \Theta [0; 2-]. By Corollary 4.7 these can only happen at points
and at their translates through 2- in the OE direction.
Since the integrand in (4.8) is 2-periodic in OE we can write
where the integrand now has singularities at three points in the domain of integration:
We split (5.2) into the sum of integrals over [0; -] \Theta [0; -] and [0; -] \Theta [-; 2-] and apply
the tensor product Gauss-Legendre rule with M points in each coordinate direction to
each integral. Since the singularities lie on the boundaries of the integration domains we
expect this to converge reasonably well.
We solved (1.1) with above. The results are given in

Table

1 and they indicate that in fact the quadrature rule yields solutions to the integral
equation and approximate potentials which converge with order roughly O(exp(\Gamma1:2M)).
The convergence is less regular for the potential at the arbitrarily chosen point -   . These
results are surprisingly good (we have no proof of the exponential convergence), and they
indicate that we have no problems from the weak singularities in (4.8).

Table

A.
Next, to check the performance of the collocation method in the absence of quadrature
errors we solved (1.1) with From Table 1, we
expect quadrature error to pollute at worst (about) the 11th decimal place of our solution.
The results are in Table 2 and show a convergence rate for approximate potentials which
is roughly of order O(exp(\Gamma2:5N)).
9

Table


Table

2 shows that good solutions to the boundary value problem are obtained by
solving a system of very small dimension (e.g. when are 86 degrees of
freedom). A similar accuracy is reported in [2]. However since the method used there is
a Galerkin method more quadrature is needed to assemble the stiffness matrix (for the
same number of degrees of freedom). The program used to compute these solutions is
very simple and compact, was written in MATLAB ([18]) and runs interactively on a Sun
4 workstation with many concurrent users and also on stand alone PC-486.
Finally it is interesting to try to reduce the overall cost of the method by balancing the
quadrature error with the underlying error in the collocation method. To do this let u e
denote the (theoretical) collocation solution (if the collocation integrals are done exactly)
and let uN denote the corresponding solution in the presence of quadrature. Then
We know that the first term on the right-hand side of (5.3) converges exponentially. We
want to ensure the second does also. Recall that from x3
suppose that denotes the quadrature approximation
to ENK (using the M point Gauss rule described above). Then a simple manipulation
shows
Thus, assuming stability for using the fact that u e
u, we expect that
where SN is as in x3. The right-hand side of (5.4) is the uniform norm of a matrix with
dimension O(N 2 ), each entry of which is converging exponentially. Thus it is reasonable
to conjecture
for some fi. We checked the sharpness of the N-dependence in this estimate by solving
Case A. The results are in Table 3 and confirm the
growth conjectured in (5.5). (Note u e
u in this example.)
9

Table

A.
To ensure exponential convergence for u e
for some ff 0 ? 0) it is then clear from (5.5) that we must choose
In principle we would choose ff 0 equal to ff where
and then M by (5.6), but in practice ff, and fi are unknown in general, and ff 0 , and fi in
would have to be chosen by experiment. Here we can verify this analysis since Table
In (5.6) we put 2. The results
are in

Table

4 and exhibit a very similar convergence and accuracy to Table 2, but much
less quadrature is needed.
9 23 7.917(-10) 2.7

Table

Finally, we tried the method on the long slender ellipsoid 10. The results
are in

Table

5 and exhibit similar convergence to Table 4 but with a somewhat larger
asymptotic constant. The accuracy of 3 \Theta 10 \Gamma6 with 86 degrees of freedom is still very
good.
9 23 1.521(-7) 1.6

Table




--R

The numerical solution of Laplace's equation in three dimensions II.
The numerical solution of Laplace's equation in three dimensions.
Algorithm 629: An integral equation program for Laplace's equation in Three dimensions.
A discrete Galerkin method for first kind integral equations with a logarithmic kernel.
Galerkin methods for solving single layer integral equations in three dimen- sions
Integral Equation Methods in Scattering Theory.
Inverse acoustic and electromagnetic scattering theory.
Spectral algorithms for vector elliptic equations in a spherical gap.
A pseudospectral approach for polar and spherical geometries.
A review of pseudospectral methods for solving partial differential equations.
Geophysical data inversion methods and applications.
A pseudospectral 3D boundary integral method applied to a nonlinear model problem from finite elasticity.
Linear Integral Equations.
On the numerical solution of a logarithmic integral equation of the first kind for the Helmholtz equation.
The fast Fourier transform and the numerical solution of one-dimensional boundary integral equations
The numerical solution of Helmholtz's equation for the exterior Dirichlet problem in three dimensions.
Approximation of Functions.
Mathworks Inc.
A spectral Galerkin method for a boundary integral equation.
Nonlinear Elastic Deformations.
Fourier series on spheres.
Polynomial interpolation and hyperinterpolation over general regions.
On the spectral approximation of discrete scalar and vector functions on the sphere.
Cubature for the sphere and the discrete spherical harmonic transform.
Die numerische approximation von randintegraloperatoren f?
Trigonometric Series.
--TR

--CTR
M. Ganesh , I. G. Graham, A high-order algorithm for obstacle scattering in three dimensions, Journal of Computational Physics, v.198 n.1, p.211-242, 20 July 2004
Lexing Ying , George Biros , Denis Zorin, A high-order 3D boundary integral equation solver for elliptic PDEs in smooth domains, Journal of Computational Physics, v.219 n.1, p.247-275, 20 November 2006
Piotr Boronski, Spectral method for matching exterior and interior elliptic problems, Journal of Computational Physics, v.225 n.1, p.449-463, July, 2007
fourier approximation;pseudospectral method;boundary integral equation;three-dimensional potential problems;collocation