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 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
The 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?