|Home | About | Journals | Submit | Contact Us | Français|
This paper proposes a new technique to speed up the computation of the matrix of spectral collocation discretizations of surface single and double layer operators over a spheroid. The layer densities are approximated by a spectral expansion of spherical harmonics and the spectral collocation method is then used to solve surface integral equations of potential problems in a spheroid. With the proposed technique, the computation cost of collocation matrix entries is reduced from (M2N4) to (MN4), where N2 is the number of spherical harmonics (i.e., size of the matrix) and M is the number of one-dimensional integration quadrature points. Numerical results demonstrate the spectral accuracy of the method.
In this paper we consider a spectral collocation method for the surface integral equations arising from three-dimensional potential problems
where f and g are given, r,r′ S and S is the surface of a three-dimensional spheroid defined by
The surface is called a prolate spheroid if L > R, and an oblate spheroid if L < R, and if L = R it is a sphere. In Eqs. (1.1) and (1.2), the unknown scalar functions μ(r) and ρ(r) are called single and double layer density functions, respectively, and n is the unit inner normal direction at a field point r on S. The potential problems arise from many fields of physics and engineering  through boundary integral formulations of differential equations of elliptic-type. In particular, in biological applications, spheroidal geometries are often used to model some types of molecule in the study of electrostatics  and the free energy .
There are many issues [5, 19] in accurate and efficient numerical approximations of potential problems, including the choices of basis functions, the treatment of singularities and fast solvers for the resulting discrete algebraic systems. For problems defined on a smooth boundary, it is natural to consider spectral methods [6, 18] due to its infinite-order convergence, so-called the spectral convergence. The spectral methods have been applied for numerical solutions of problems in spheroidal geometries, for example, in geodesy and in global atmospheric models [9, 17], by using the spherical harmonics as the basis functions. This is due to the fact that the spherical harmonics offer nice properties in approximating spherical functions with exponential convergence and the ease in handling pole singularities . For the boundary integral equations considered in this paper, the spherical harmonics have been used by Atkinson [3,4], Graham and Sloan , Chen , and Ganesh  to construct spectral methods for three-dimensional potential problems. Also, double Fourier series have found applications in solving potential problems .
For spectral methods of surface integral equations, the calculation of collocation coefficients (or Galerkin coefficients in Galerkin methods) remains one of the computational challenges as numerical quadratures over the full boundary surface is required in calculating each entry of the coefficient matrix. For example, if N2 basis functions and M2-point two-dimensional quadratures are used, the total operations of order (M2N4) are required to produce the coefficient matrix of a collocation method. Here M can become large with the increase of basis function order N in order to maintain the high-frequency information from the higher order basis functions. Even with relatively small N, the cost can be still high. In this paper, we will propose a new technique which will reduce the two-dimensional integral to one dimensional one. This is made possible by an analytical formula along the longitudinal direction for the spheroidal surface with the help of a hypergeometric function. The hypergeometric function can be calculated fast even at singular points thanks to linear transformation formulas. As a result, the complexity of filling up the collocation matrix is then reduced to the order of (MN4). With some generalization, this technique can be extended to boundary element methods and artificial boundary methods [13, 20] where the spheroidal boundary is usually adopted.
The rest of this paper is organized as follows. In Section 2, basis functions of spherical harmonics are introduced. In Sections 3 and 4, we discuss the technique for fast calculations of the collocation coefficients for the single and double layer potentials, respectively. In Section 5, the treatment of singularities is discussed. Numerical test examples are then performed in Section 6 with conclusions given in Section 7.
Let us consider Eq. (1.1) with a spherical coordinate representation for the spheroidal surface,
where θ and ϕ represent latitude and longitude of the surface with 0≤θ≤π and 0≤ϕ<2π. It is noted that the Jacobian of the coordinate transformation between the Cartesian and spherical coordinates is a function only of variable θ,
Eq. (1.1) can then be re-written as a two-dimensional integral as follows:
The grid collocation points along the longitudinal and latitudinal directions are selected as
and (N+1)2 is the total number of collocation points, and N is chosen to be an even integer.
A standard finite dimensional orthogonal basis over a spherical surface is the (N+1)2 spherical harmonics,
for −n≤m≤n and 0≤n≤N, where are the associate Legendre functions of degree n and order m, and i the imaginary unit. The associate Legendre polynomials can be defined by the Rodrigue’s formula
where Pn(x) is the ordinary Legendre polynomial of degree n.
For a given basis, a function defined on the surface of the sphere can be approximated by a linear combination of the basis functions as
in which are constant coefficients, and the total number of the basis functions is (N+1)2. If the function μ(θ,ϕ) is infinitely differentiable, then it can be proved that the approximation is spectrally convergent .
Although Galerkin-type methods have been studied more thoroughly, we will adopt a collocation method for its ease of implementation and simpler filling of the discretization matrix. With the (N+1)2 collocation points in (2.4), the spectral collocation method for the single layer surface integral equation yields the following linear system for the unknowns ,
with j,l=0,,N, where the collocation coefficients are
After numerically approximating the integrals of coefficients with quadratures over the spheroid surface, the linear algebraic systemcan be solved through the Gauss elimination direct solver or other iterative solvers.
An important issue for the spectral collocation is how to evaluate the collocation matrix entry accurately and efficiently. Naively, a 2M-trapezoidal rule in ϕ and an M-Gaussian quadrature in θ will require (M2) operations for the calculation of each entry. In this section we will first cast the integral (2.9) analytically with a special function along the latitudinal direction such that the complexity is reduced to (M). For this purpose, let us write the distance between the field point and the collocation point as
where a(θ) and b(θ) are defined by
Re-arranging integrations in (2.9), we reach
Noticing that the functions of ϕ are 2π-periodic and using the oddness of function sinmϕ, we have
The remaining step is to compute the integration on the right hand side (see Fig. 1), which can be done with the following Theorem 3.1. Let us first recall the generalized hypergeometric function 2F1(α,β;γ;x) defined by
with (α)k=α(α+1)(α+k−1) being the rising factorial, and we have
For b satisfying |b| < 1, we then have
Due to |b| < 1, we have the Taylor series expansion
We can write the integral as
The integral becomes
This completes the proof.
The objective of Theorem 3.1 is to express the integral in the form of a series which can be identified with the hypergeometric special function. The calculation of the series is convenient due to its geometric convergence and can be done by an efficient algorithm. Table 1 contains the required truncations of the hypergeometric series for a 10−6 error tolerance. Even if the parameter b is close to the singularity, for example b ≥ 0.9, the hypergeometric function can be evaluated by the following linear transformation formulas  (pp. 559):
where ψ(z) is the logarithmic derivative of the Gamma function Γ,
Furthermore, for m = 1,2,, the following linear transformation formulas will be used in the double layer potential in the next section:
Once the two-dimensional integral is reduced to a one-dimensional one, the singularity of the potential function is easier to treat. Actually, the induced hypergeometric function is a logarithmic singularity, which can be separated from the integrand (to be discussed in Section 5).
Extension to the double layer potential problem is straightforward with some more involved manipulations. First note that
and, for the surface of the spheroid defined by (1.3), the unit normal direction inward S is
Using the spherical coordinate representations of r and rjl yields
where a(θ) and b(θ) is the same as in the definition for the single layer potential, and A(θ) and B(θ) are defined by
Particularly, if the spheroid is a sphere L = R, then a=2A and b=B/A, Eq. (4.4) reduces to a single layer kernel.
As in the case of the single layer potential, we need to do an analytical integration to reduce the complexity of the two dimensional integrals.
for |b| < 1, we have the following identities.
If |b| < 1, then
and, furthermore, for m ≥ 1,
The proof is similar to that of Theorem 3.1 by using the Taylor series expansion, thus omitted.
Finally, we are ready to evaluate the two dimensional integral in the definition for the collocation coefficient . Note that the integral of ϕ was done by the analytical formula above, which leaves us just one-dimensional singular integrals to handle. We rewrite in the form of
Note that b(θ) → −1 as θ → θj. For the single layer potential problem we have a logarithmic singularity due to the following fact that
Similarly, for the double layer potential, we have
Therefore, the integrals in both the single and double layer potential problems have logarithmic singularities when b tends to 1 (i.e., when θ approaches to θj), which can then be treated by a simple singularity substraction technique. Assuming
then we can write I(G) into the following form
Here the second term can be integrated analytically as
As for the first term, because the integrand becomes zero at the point θj, it is an integral without singularity and can be then approximated accurately by Gaussian quadratures.
To test the method, we first consider the Dirichlet problem of the Laplace equation inside a spheroid,
The solution of this problem can be expressed either as a single layer potential
or as a double layer potential
The density functions μ(r) and ρ(r) are then determined by letting r′ tend to the surface to yield the following surface integral equations
for the single and double layer potentials, respectively.
To test the performance, we use two analytical solutions of the Laplace equation
and the surface S defined by (1.3) with L = 1.5 and R = 1. In Table 2–Table 4, The errors for different number of basis functions are listed. The (N+1)2 one-dimensional collocation coefficients are calculated using a 64-th order Gaussian quadrature. Exponential convergence can be observed from the results, and the single layer problem does have a smaller error than the double layer one.
In the second example, we consider the following electrostatic problem: a unit point charge is located at rs inside a spherical cavity Ω modeled as a uniform dielectric medium with dielectric constant ϵi; the cavity is embedded into a homogeneous medium with dielectric constant ϵe. The potential satisfies the Poisson equation:
The potential is continuous across spherical boundary and the normal derivative has a jump related to the ratio of the dielectric constant. Thus, we have boundary conditions:
where Φi and Φe are the potential functions inside and outside the cavity, and r− and r+ are inner and outer limits, respectively.
This model, being widely used in computing electrostatic contributions of salvation free energy of biomolecules , can be transformed into boundary integral equations by using the fundamental solution of the Poisson equation,
Using the Green’s second theorem, we can obtain
When r approaches the surface S=Ω, by using limiting properties of double-layer potentials at boundary, the above two equations reduce to
for r S, where we have used the two boundary conditions.
We then solve the integral equations by the proposed method for the two unknowns Φi(r) and Φi(r)/n. A unit spherical cavity Ω is used, for which the analytical solution is known by the Kirkwood series expansion [7,14] in terms of the Legendre polynomials. In the calculations, we set ϵi = 2, ϵe = 80, and the charge position rs = (0,0,rs) for rs=0.4, 0.6, and 0.8. The one-dimensional collocation coefficients are again calculated using a 64-th order Gaussian quadrature. We compute the L2 relative errors of the solutions of Φi(r), which are shown in Fig. 2. It can seen the results rapidly converge to the exact solutions with the magnitude of the relative error around 10−5, illustrating the spectral convergence of the proposed algorithm.
A new technique was proposed to reduce the complexity of two dimensional surface single and double layer potential integrals over a spheroid and as a result, a fast spectral collocation method is obtained for the surface integral equations of potential problems in a spheroid. Numerical results shows the spectral accuracy and efficiency of the collocation method.
We thanks some mathematicians on Google math group for helpful discussions on the hypergeometric function. Financial support for this work was provided by the National Institutes of Health (grant number: 1R01GM083600-01). Z. Xu is also partially supported by the Charlotte Research Institute through a Duke Postdoctoral Fellowship.
AMS subject classifications: 31B10, 33C90, 65R20, 76M22