计算科学 - 不收敛势的多极展开计算 - 吾爱随笔录

根据Beatson 和 Greengard 关于 FMM 的短期课程：（等式 5.15 和 5.16 设置 k=1，q=1）

我们可以近似一个潜在的 $\phi = 1/(r-R)$ 使用：

\frac{1}{| \vec{r} - \vec{R} |} = \sum_{n = 0}^{\infty} \frac{r^{n}}{R^{n + 1}} \frac{4 π}{(2 n + 1)} \sum_{m = - n}^{n} Y_{n}^{- m} (θ, ϕ) Y_{l}^{m} (θ^{'}, ϕ^{'})

${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')$

我在Python中试过这个，错误是这样的：

| ϕ (P) - \sum_{n = 0}^{p} \sum_{m = - n}^{n} (\dots) | \leq \frac{1}{R - r} {(\frac{r}{R})}^{p + 1}

$\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

使用这种方法我得到错误的结果，举一个简单的例子

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))

这会输出完全错误的结果：

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

我是否应该调查这是否是浮点错误？如果是这样，我该怎么做？