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.

    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 "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:

Direct calculation:0.0099995000375
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?

  • 1
    $\begingroup$ 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." $\endgroup$ – Daniel Shapero Jan 25 '15 at 20:52
  • $\begingroup$ 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. $\endgroup$ – 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.

| cite | improve this answer | |

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.

| cite | improve this answer | |
  • $\begingroup$ 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. $\endgroup$ – OscarB Jan 27 '15 at 21:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.