PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Commun Comput Phys. Author manuscript; available in PMC 2010 April 21.
Published in final edited form as:
Commun Comput Phys. 2009; 6: 625–638.
PMCID: PMC2857787
NIHMSID: NIHMS121781

Fast Spectral Collocation Method for Surface Integral Equations of Potential Problems in a Spheroid

Abstract

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 O(M2N4) to O(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.

Keywords: Boundary integral equations, three-dimensional potential problems, collocation spectral methods, spherical harmonics, Fourier series, hypergeometric functions

1 Introduction

In this paper we consider a spectral collocation method for the surface integral equations arising from three-dimensional potential problems

Sμ(r)1|rr|dr=f(r),
(1.1)

Sρ(r)nr(1|rr|)dr=g(r),
(1.2)

where f and g are given, r,r′ [set membership] S and S is the surface of a three-dimensional spheroid defined by

S={r=(x,y,z)|x2+y2R2+z2L2=1}.
(1.3)

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 [16] 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 [7] and the free energy [2].

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 [6]. For the boundary integral equations considered in this paper, the spherical harmonics have been used by Atkinson [3,4], Graham and Sloan [12], Chen [8], and Ganesh [10] to construct spectral methods for three-dimensional potential problems. Also, double Fourier series have found applications in solving potential problems [11].

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 O(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 O(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.

2 Spectral approximations in spherical coordinates

Let us consider Eq. (1.1) with a spherical coordinate representation for the spheroidal surface,

r=(x,y,z)=(Rcosϕsinθ,Rsinϕsinθ,Lcosθ),
(2.1)

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 θ,

J(θ)=|rϕ×rθ|=(LR)2sin4θ+R4sin2θcos2θ.
(2.2)

Eq. (1.1) can then be re-written as a two-dimensional integral as follows:

02π0πμ(r)|rr|J(θ)dθdϕ=f(θ,ϕ).
(2.3)

The grid collocation points along the longitudinal and latitudinal directions are selected as

ϕl=2πlN+1,andθj=π(2j1)2(N+1),  for  l,j=0,,N,
(2.4)

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,

Ynm(θ,ϕ)=2n+14π(n|m|)!(n+|m|)!Pn|m|(cosθ)eimϕ,
(2.5)

for −nmn and 0≤nN, where Pnm 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

Pnm(x)=(1)m(1x2)m/2dmdxmPn(x),
(2.6)

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

μ(θ,ϕ)n=0mN=nnμnmYnm(θ,ϕ),
(2.7)

in which μnm 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 [17].

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 μnm,

n=0mN=nnAnmμnm=f(θj,ϕl),
(2.8)

with j,l=0,(...),N, where the collocation coefficients are

Anm(θj,ϕl)=02π0πYnm(θ,ϕ)J(θ)|rjlr|dθdϕ.
(2.9)

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.

3 Fast calculations of the collocation coefficients Anm

An important issue for the spectral collocation is how to evaluate the collocation matrix entry Anm accurately and efficiently. Naively, a 2M-trapezoidal rule in ϕ and an M-Gaussian quadrature in θ will require O(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 O(M). For this purpose, let us write the distance between the field point and the collocation point as

|rrjl|=R2(cosϕsinθcosϕlsinθj)2+R2(sinϕsinθsinϕlsinθj)2+L2(cosθcosθj)2=[R2(sin2θ+sin2θj)+L2(cosθcosθj)2]2R2sinθsinθjcos(ϕϕl)=a(θ)1+b(θ)cos(ϕϕl),
(3.1)

where a(θ) and b(θ) are defined by

a(θ)=R2(sin2θ+sin2θj)+L2(cosθcosθj)2,
(3.2a)

b(θ)=2R2sinθsinθj/a(θ).
(3.2b)

Re-arranging integrations in (2.9), we reach

Anm(θj,ϕl)=0πΘnm(θ)J(θ)a(θ)02πeimϕ1+b(θ)cos(ϕϕl)dϕdθ,
(3.3)

with

Θnm(θ)=2n+14π(n|m|)!(n+|m|)!Pn|m|(cosθ).
(3.4)

Noticing that the functions of ϕ are 2π-periodic and using the oddness of function sinmϕ, we have

02πeimϕ1+b(θ)cos(ϕϕl)dϕ=eimϕl02πcosmϕ1+b(θ)cosϕdϕ.
(3.5)

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

F12(α,β;γ;x)=k=0(α)k(β)k(γ)kxkk!,

with (α)k=α(α+1)(...)(α+k−1) being the rising factorial, and we have

Figure 1
Illustration of formula (3.6) for m = 0,(...),4. The integral is singular at b = ±1.

Theorem 3.1

For b satisfying |b| < 1, we then have

02πcosmϕ1+bcosϕdϕ=212m(b)mπ(2m1)!!m!F12(2m+14,2m+34;m+1;b2).
(3.6)

Proof

Due to |b| < 1, we have the Taylor series expansion

11+bcosϕ=k=0(14)k(2k)!(k!)2(bcosϕ)k.
(3.7)

We can write the integral as

I=k=0(b4)k(2k)!(k!)202πcosmϕcoskϕdϕ.
(3.8)

Since

02πcosmϕcoskϕdϕ={21k(kkm2)π,for0mk,and(km  mod2)=0,0otherwise.
(3.9)

The integral becomes

I=i=0(b4)m+2i(m+2ii)21m2iπ(2m+4i)![(m+2i)!]2=212m(b)mπ(2m1)!!m!F12(2m+14,2m+34;m+1;b2).
(3.10)

This completes the proof.

Therefore, we obtain from Eqs. (3.3)(3.6) the following method to compute the collocation coefficients:

Anm(θj,ϕl)=212mπ(2m1)!!eimϕlm!0πΘnm(θ)J(θ)(b)maF12(2m+14,2m+34;m+1;b2)dθ,
(3.11)

where a,b and J are functions of θ defined in Eq. (3.2) and Eq. (2.2), respectively.

Remark 3.1

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 [1] (pp. 559):

F12(α,β;α+β;x)=Γ(α+β)Γ(α)Γ(β)n=0(α)n(β)n(n!)2[2ψ(n+1)ψ(α+n)ψ(β+n)log(1x)](1x)n,

where ψ(z) is the logarithmic derivative of the Gamma function Γ,

ψ(z)=d[logΓ(z)]dz=Γ(z)Γ(z).

Table 1
Truncation terms of the hypergeometric series F12(2m+14,2m+34;m+1;b2) needed for a 1×10−6 error tolerance for varying b and m = 0,(...),7. For large b, the series can be approximated by the linear transformation formulas.

Furthermore, for m = 1,2,(...), the following linear transformation formulas will be used in the double layer potential in the next section:

F12(α,β;α+βm;x)=Γ(m)Γ(α+βm)Γ(α)Γ(β)(x)mn=0m1(αm)n(βm)nn!(1m)n(1x)n+(1)mΓ(α+βm)Γ(αm)Γ(βm)n=0(α)n(β)nn!(n+m)![ψ(n+1)+ψ(n+m+1)ψ(α+n)ψ(β+n)log(1x)](1x)n.

Remark 3.2

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).

4 Double layer potential

Extension to the double layer potential problem is straightforward with some more involved manipulations. First note that

nr(1|rjlr|)=nr(1|rjlr|)=nr(rrjl|rjlr|3),
(4.1)

and, for the surface of the spheroid defined by (1.3), the unit normal direction inward S is

nr=rϕ×rθ|rϕ×rθ|=(Lcosϕsinθ,Lsinϕsinθ,Rcosθ)L2sin2θ+R2cos2θ.
(4.2)

Using the spherical coordinate representations of r and rjl yields

nr(rrjl)=LR[1cosθcosθjsinθsinθjcos(ϕϕl)]L2sin2θ+R2cos2θ;
(4.3)

thus,

nr(1|rjlr|)=LR(1cosθcosθj)LRsinθsinθjcos(ϕϕl)a(θ)3/2[1+b(θ)cos(ϕϕl)]3/2L2sin2θ+R2cos2θ=1a(θ)3/2L2sin2θ+R2cos2θA(θ)+B(θ)cos(ϕϕl)[1+b(θ)cos(ϕϕl)]3/2,
(4.4)

where a(θ) and b(θ) is the same as in the definition for the single layer potential, and A(θ) and B(θ) are defined by

A(θ)=LR(1cosθcosθj),
(4.5a)

B(θ)=LRsinθsinθj.
(4.5b)

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

II=02πcosmϕ[A(θ)+B(θ)cosϕ][1+b(θ)cosϕ]3/2dϕ
(4.6)

for |b| < 1, we have the following identities.

Theorem 4.1

If |b| < 1, then

02πcosmϕ(1+bcosϕ)3/2dϕ=22m(b)mπΓ(32+m)Γ(m+1)F12(2m+34,2m+55;m+1;b2),
(4.7)

and, furthermore, for m ≥ 1,

02πcosmϕcosϕ(1+bcosϕ)3/2dϕ=21m(b)m1π{2Γ(m+12)Γ(m)F12(2m+14,2m+34;m+1;b2)+b2Γ(m+52)Γ(m+2)F12(2m+54,2m+74;m+2;b2)}.
(4.8)

The proof is similar to that of Theorem 3.1 by using the Taylor series expansion, thus omitted.

5 Quadrature for singular functions

Finally, we are ready to evaluate the two dimensional integral in the definition for the collocation coefficient Anm. Note that the integral of ϕ was done by the analytical formula above, which leaves us just one-dimensional singular integrals to handle. We rewrite Anm in the form of

I(G)=0πG(θ)dθ.
(5.1)

Note that b(θ) → −1 as θ → θj. For the single layer potential problem we have a logarithmic singularity due to the following fact that

limb11log(1+b)02πcosmϕ1+bcosϕdϕ=212mπ(2m1)!!Γ(m+1)m!Γ(2m+14)Γ(2m+34).
(5.2)

Similarly, for the double layer potential, we have

limθθj1log(1+b(θ))02πcosmϕ[A(θ)+B(θ)cosϕ][1+b(θ)cosϕ]3/2dϕ=LRsin2θjlimb11log(1+b)02πcosmϕ1+bcosϕdϕ.
(5.3)

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

limθθjG(θ)log(|θθj|)=𝒢(θj),
(5.4)

then we can write I(G) into the following form

I(G)=0π[G(θ)𝒢(θj)log(|θθj|)]dθ+𝒢(θj)0πlog(|θθj|)dθ.
(5.5)

Here the second term can be integrated analytically as

𝒢(θj)0πlog(|θθj|)dθ=𝒢(θj)[π+θjlogθj+(πθj)log(πθj)].
(5.6)

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.

6 Numerical examples

6.1 Single and double layer potentials

To test the method, we first consider the Dirichlet problem of the Laplace equation inside a spheroid,

Δu(r)=0,rΩ,
(6.1)

u(r)=f(r),whenrS=Ω.
(6.2)

The solution of this problem can be expressed either as a single layer potential

u(r)=14πSμ(r)|rr|dr,rΩ,
(6.3)

or as a double layer potential

u(r)=Sρ(r)nr(1|rr|)dr,rΩ.
(6.4)

The density functions μ(r) and ρ(r) are then determined by letting r′ tend to the surface to yield the following surface integral equations

14πSμ(r)|rr|dr=f(r),rS,
(6.5)

and

2πρ(r)+Sρ(r)nr(1|rr|)dr=f(r),rS,
(6.6)

for the single and double layer potentials, respectively.

To test the performance, we use two analytical solutions of the Laplace equation

u(r)=x,r=(x,y,z),
(6.7)

and

u(r)=1ep(epxcospy+epzsinpz),forp=1,2,
(6.8)

and the surface S defined by (1.3) with L = 1.5 and R = 1. In Table 2Table 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.

Table 2
Errors |uuN| at eight points inside the spheroid for the analytical solution u(r)=x.
Table 4
Errors |uuN| at eight points inside the spheroid for the analytical solution u(r)=1e2(e2xcos2y+e2zsin2z).

6.2 Solving boundary integral formulations for the Poisson equation

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:

ϵiΔΦ=δ(rrs),rΩ.
(6.9)

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:

Φi(r)=Φe(r+),ϵiΦir(r)=ϵeΦer(r+),
(6.10)

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 [15], can be transformed into boundary integral equations by using the fundamental solution of the Poisson equation,

G(r,r)=14π|rr|.
(6.11)

Using the Green’s second theorem, we can obtain

ϕi(r)=Ω[G(r,r)Φi(r)nG(r,r)nΦi(r)]dr+1ϵiG(r,rs),
(6.12)

Φe(r)=Ω[G(r,r)Φe(r)n+G(r,r)nΦe(r)]dr.
(6.13)

When r approaches the surface S=[partial differential]Ω, by using limiting properties of double-layer potentials at boundary, the above two equations reduce to

12Φi(r)=Ω[G(r,r)Φi(r)nG(r,r)nΦi(r)]dr+1ϵiG(r,rs),
(6.14)

12Φi(r)=Ω[G(r,r)ϵiϵeΦi(r)n+G(r,r)nΦi(r)]dr,
(6.15)

for r [set membership] 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 [partial differential]Φi(r)/[partial differential]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.

Figure 2
Relative L2 errors for the potential problem of the Poisson equation for charge locations rs=0.4, 0.6 and 0.8.

7 Conclusions

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.

Table 3
Errors |uuN| at eight points inside the spheroid for the analytical solution u(r)=1e(excosy+ezsinz).

Acknowledgments

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.

Footnotes

AMS subject classifications: 31B10, 33C90, 65R20, 76M22

References

1. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover; 1964.
2. Ambia-Garrido J, Pettitt BM. Free energy calculations for DNA near surfaces using an ellipsoidal geometry. Commun. Comput. Phys. 2008;3:1117–1131. [PMC free article] [PubMed]
3. Atkinson KE. The numerical solution of Laplace’s equation in three dimensions. SIAM J. Numer. Anal. 1982;19:263–274.
4. Atkinson KE. Algorithm 629: An integral equation program for Laplace’s equation in three dimensions. ACM Trans. Math. Soft. 1985;11:85–96.
5. Atkinson KE. The numerical solution of boundary integral equations. In: Duff I, Watson G, editors. The State of the Art in Numerical Analysis. Clarendon Press; 1997. pp. 223–259.
6. Boyd JP. Chebyshev and Fourier Spectral Methods. 2nd Edition. New York: Dover Publications; 2001.
7. Cai W, Deng S, Jacobs D. Extending the fast multipole method to charges inside or outside a dielectric sphere. J. Comput. Phys. 2007;223:846–864.
8. Chen Y. Ph.D. thesis. University of Iowa; 1994. Galerkin methods for solving single layer integral equations in three dimensions.
9. Cheong HB. Double Fourier series on a sphere: Applications to elliptic and vorticity equations. J. Comput. Phys. 2000;157:327–349.
10. Ganesh M. Fully discrete spectral methods for boundary integral equations on slender spheroids. J. Comput. Appl. Math. 2004;164–165:307–322.
11. Ganesh M, Graham IG, Sivaloganathan J. A new spectral boundary integral collocation method for three-dimensional potential problems. SIAM J. Numer. Anal. 1998;35:778–805.
12. Graham IG, Sloan IH. Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in R3. Numer. Math. 2002;92:289–323.
13. Huang HY, Yu DH. Natural boundary elementmethod for three dimensional exterior harmonic problem with an inner prolate spheroid boundary. J. Comput. Math. 2006;24:193–208.
14. Kirkwood JG. Theory of solutions of molecules containing widely separated charges with special applications to awitterions. J. Chem. Phys. 1934;2:351–361.
15. Lu BZ, Zhou YC, Holst MJ, McCammon JA. Recent progress in numerical methods for the Poisson-Boltzmann equation in biophysical applications. Commun. Comput. Phys. 2008;3:973–1009.
16. Nedelec JC. New York: Springer-Verlag; 2001. Acoustic and Electromagnetic Equations – Integral Representations for Harmonic problems.
17. Orszag SA. Fourier series on spheres. Mon. Weather Rev. 1974;102:56–75.
18. Shen J, Tang T. Spectral and High-Order Methods with Applications. Beijing: Science Press; 2006.
19. Ying L, Biros G, Zorin D. A high-order 3D boundary integral equation solver for elliptic PDEs in smooth domains. J. Comput. Phys. 2006;219:247–275.
20. Yu D. Mathematical Theory of Natural Boundary Element Methods. Beijing: Science Press; 1993.