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().realUsing this method I am getting wrong results, taking a simple example p = 5r1 = 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 % potentialapprox_potential = potential_expansion(p, r1, theta1, phi1, r2, theta2, phi2)print Approximation: %s % approx_potentialprintprint 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.010101010101Direct calculation:0.0099995000375Error:0.000101510063503Upper bound on error:1.0101010101e-14Should I investigate whether or not this is a floating point error? If so, how can I go about this?
Computation of multipole expansion of potential not converging
python;floating point;approximation algorithms
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.