I'm trying to solve the Poisson equation with pure Neumann boundary conditions,$$\nabla^2\phi = \rho \quad in \quad \Omega\\\mathbf{\nabla}\phi \cdot \mathbf{n} = 0 \quad on \quad \partial \Omega$$using a Fourier transform method I found in Numerical Recipes. The method uses a discrete cosine transform, if you don't have access to the book, you can find a derivation here. I tried implementing the algorithm in Python; my code is listed below:import numpy as npimport scipy.sparse as sparseimport scipy.fftpack as fftif __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 potentialThe charge density is the following,[[ 0.11111111  0.11111111  0.11111111] [ 0.11111111 -0.88888889  0.11111111] [ 0.11111111  0.11111111  0.11111111]]while, multiplying the solution with the Laplacian $L$ gives,[[ 0.25  0.25  0.25] [ 0.25 -0.75  0.25] [ 0.25  0.25  0.25]]instead of the same results as above.I've been staring to the code for some time now and am unable to understand what I'm doing wrong. I've even checked to see if scipy's discrete cosine transform gives the correct results and that seems fine too.If anyone could point out my mistake, I would be really grateful!EDIT: I found out that if I multiply the solution by $L - J$, where $J$ is a matrix filled with ones, instead of just $L$, I do get the expected charge density. Why is that?
Poisson equation with Neumann boundary conditions
finite difference;python;poisson