我正在尝试用纯诺依曼边界条件求解泊松方程，

$$\nabla^2\phi = \rho \quad in \quad \Omega\\ \mathbf{\nabla}\phi \cdot \mathbf{n} = 0 \quad on \quad \partial \Omega$$ 使用我在Numerical Recipes中找到的傅里叶变换方法。该方法使用了离散余弦变换，如果你没有接触到这本书，你可以在这里找到一个推导。我尝试在 Python 中实现该算法；我的代码如下：

import numpy as np
import scipy.sparse as sparse
import scipy.fftpack as fft

if __name__ == '__main__':
    shape = (3, 3)
    nx, ny = shape

    charges = np.zeros(shape)
    charges[:] = 1.0 / (nx * ny)
    charges[nx / 2, ny / 2] = 1.0 / (nx * ny) - 1.0
    print charges

    charges = charges.flatten()

#Build Laplacian
    ex  = np.append(np.ones(nx - 2), [2, 2])
    ey  = np.append(np.ones(ny - 2), [2, 2])
    Dxx = sparse.spdiags([ex, -2 * np.ones(nx), ex[::-1]], [-1, 0, 1], nx, nx)
    Dyy = sparse.spdiags([ey, -2 * np.ones(ny), ey[::-1]], [-1, 0, 1], ny, ny)

    L = sparse.kronsum(Dxx, Dyy).todense()
###############

#Fourier method
    rhofft = np.zeros(shape, dtype = float)
    for i in range(shape[0]):
        rhofft[i,:] = fft.dct(charges.reshape(shape)[i,:], type = 1) / (shape[1] - 1.0)
    for j in range(shape[1]):
        rhofft[:,j] = fft.dct(rhofft[:,j], type = 1) / (shape[0] - 1.0)

    for i in range(shape[0]):
        for j in range(shape[1]):
            factor = 2.0 * (np.cos((np.pi * i) / (shape[0] - 1)) + np.cos((np.pi * j) / (shape[1] - 1)) - 2.0)
            if factor != 0.0:
                rhofft[i, j] /= factor
            else:
                rhofft[i, j] = 0.0

    potential = np.zeros(shape, dtype = float)
    for i in range(shape[0]):
        potential[i,:] = 0.5 * fft.dct(rhofft[i,:], type = 1)
    for j in range(shape[1]):
        potential[:,j] = 0.5 * fft.dct(potential[:,j], type = 1)

################
    print np.dot(L, potential.flatten()).reshape(shape)
    print potential

电荷密度如下，

[[ 0.11111111 0.11111111 0.11111111]
 [ 0.11111111 -0.88888889 0.11111111]
 [0.11111111 0.11111111 0.11111111]]

同时，将解与拉普拉斯算子相乘$L$给，

[[ 0.25 0.25 0.25]
 [ 0.25 -0.75 0.25]
 [ 0.25 0.25 0.25]]

而不是与上面相同的结果。

我已经盯着代码一段时间了，无法理解我做错了什么。我什至检查了 scipy 的离散余弦变换是否给出了正确的结果，这似乎也很好。

如果有人能指出我的错误，我将非常感激！

编辑：我发现如果我将解决方案乘以$L - J$， 在哪里$J$是一个用 1 填充的矩阵，而不仅仅是$L$，我确实得到了预期的电荷密度。这是为什么？