# Computation of multipole expansion of potential not converging

According to Beatson and Greengard's short course on FMM: ( Eq. 5.15 & 5.16 setting k=1, q=1 )

We can approximate a potential $\phi = 1/(r-R)$ using:

$${1\over |\vec{r}-\vec{R}|} = \sum_{n=0}^{\infty}{r^n\over R^{n+1}} {4\pi\over(2n+1)} \sum_{m=-n}^{n} Y_n^{-m}(\theta, \phi) Y_l^m(\theta', \phi')$$

I tried this in Python, the error is such that:

$$\left| \phi(P)-\sum_{n=0}^p\sum_{m=-n}^{n} ( \cdots ) \right| \le {1 \over R-r}\left( r\over R \right)^{p+1}$$

def potential_expansion( p,
r1, theta1, phi1,
r2, theta2, phi2 ):
"""
Return inverse r potential expansion upto the pth term.

Input
======
p - terms in the expansion to return
r1, theta1, phi1 - position vector components for r
r2, theta2, phi2 - position vector components for R

"""
coefficients = np.zeros( (p+1, p+1), complex )
for n in xrange(p+1):
for m in xrange(-n, n+1):
coefficient = sph_harm( m, n, theta1, phi1 )*sph_harm( -m, n, theta2, phi2)
coefficients[n][m] = 4*pi/(2*n+1)*(r2/r1)**n/r1*coefficient
return coefficients.sum().real


Using this method I am getting wrong results, taking a simple example

p = 5

r1 = 100.
theta1 = 0.
phi1 = 0.

r2 = 1.
theta2 = pi/2.
phi2 = 0.

cosgamma = cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)*cos(phi1-phi2)
potential = 1/root( r1**2 + r2**2 - 2*r1*r2*cosgamma)
print "Direct calculation: %s" % potential

approx_potential = potential_expansion(p, r1, theta1, phi1, r2, theta2, phi2)
print "Approximation: %s" % approx_potential
print
print "Error: %s" % np.abs((approx_potential - potential))
print "Upper bound on error: %s" % (1/(r1-r2)*(r2/r1)**(p+1))


This outputs completely wrong results:

Approximation:0.010101010101
Direct calculation:0.0099995000375
Error:0.000101510063503
Upper bound on error:1.0101010101e-14


Should I investigate whether or not this is a floating point error? If so, how can I go about this?

• This might be stupid, but you should check that Greengard's notes and the scipy sph_harm function don't swap what they mean by $\theta$ and $\phi$ for the spherical angles. From the scipy documentation: "There are different conventions for the meaning of input arguments theta and phi. We take theta to be the azimuthal angle and phi to be the polar angle. It is common to see the opposite convention - that is theta as the polar angle and phi as the azimuthal angle." Jan 25, 2015 at 20:52
• That seems to be it! I now have error: 3.12250225676e-15 with upper bound on error: 1.0101010101e-14, feel free to put an answer for the question. Jan 25, 2015 at 21:02

In order to inconvenience as many people as possible, long ago, mathematicians and physicists decided to use two different conventions on whether $\theta$ or $\phi$ is the latitude angle. Greengard's notes use the physicists' convention that $\theta$ is latitude and $\phi$ is longitude, whereas Scipy uses $\theta$ for longitude and $\phi$ for latitude, so switching the order of the arguments theta and phi fixed your code.

Just a small comment: The sequence of summation may cause divergence as well. I learned in school always to sum up the coefficients from small to large magnitudes. How is the line

return coefficients.sum().real


sums up, is not specified. You may consider that.

• I believe it is more complicated than just always summing from small to large. See the very readable paper: Higham, N. J. (1993). The accuracy of floating point summation. SIAM Journal on Scientific Computing, 14(4), 783–799. Jan 27, 2015 at 21:45