Computation of multipole expansion of potential not converging

Question

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?

Answer 1

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.

Answer 2

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.