# 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." – Daniel Shapero Jan 25 '15 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. – shilov Jan 25 '15 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.
return coefficients.sum().real