Poisson equation with Neumann boundary conditions

Question

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?

It looks like your two solutions may differ by a constant. If that's true, then the discrepancy makes sense, because the Poisson equation + Neumann boundary conditions has a non-trivial kernel. If you have a solution $u(x)$ to your problem, $u(x) + c$ is also a solution, so you can only solve your problem up to an undetermined constant. — Sumedh Joshi, Jun 16 '14 at 18:51
That is true, but you still expect $\nabla^2 u(x) = \rho(x)$. The charge density is the right-hand side, which is known, i.e. the solution I get doesn't satisfy the equation I want to solve. — Grieverheart, Jun 16 '14 at 19:17
I'm not terribly familiar with Fourier methods, but the problem from a numerical standpoint is that the discrete Laplacian $L$ you'll construct will be rank-deficient, which means that some component of $\rho(x)$ might lie outside of the range space of $L$. So what you'll end up solving is a linear system $Lx = \tilde\rho$, where $\tilde\rho$ is the projection of $\rho$ onto the range space of $L$. In your case, it looks like $\int \rho(x)dx \neq 0$; try this same problem with a RHS vector $\rho$ that has zero mean, since all vectors in the range $y = Lx$ will have zero mean themselves. — Sumedh Joshi, Jun 16 '14 at 19:31
Conversely, I'm not very familiar with advanced linear-algebra concepts. From my understanding, I can somehow solve the equation using 'null space projection', but all these terms are beyond my knowledge. Btw, it happens to be the case here that $\int dx \rho(x) = 0$, if you sum all elements in the first matrix. — Grieverheart, Jun 16 '14 at 20:21
Nullspace projection in a Fourier solution method for Poisson is merely neglecting the singular constant-mode term. You do that with your "if factor != 0" statement in the above algorithm. — Peter Brune, Jun 16 '14 at 22:11
@Grieverheart Sorry I confused the RHS with the Laplacian times the solution in your original post. This is strange and I'm not sure I know why; does your Poisson matrix applied to the vector of ones give you 0? — Sumedh Joshi, Jun 17 '14 at 19:49
Do you properly norm the FFT? Try if just forward and backward gives you the original. — Vladimir F Героям слава, Jun 18 '14 at 18:09
I know this is quite old now but did you ever find a solution to this? i am also currently experiencing the same issue thanks — weddle_32, Oct 07 '14 at 18:31

score 5 · Answer 1 · edited Mar 24 '18 at 16:39

5

I think this happened because $L$ is a finite-difference operator that approximates your equation with some error (it is $\mathcal O(dx^2)$). To get small error your should increase nx and ny. If you set in your example nx = ny = 101 ($dx = 0.01$) you will get error about 1e-4.

edited Mar 24 '18 at 16:39

Anton Menshov

8,672
7
38
94

answered Oct 08 '14 at 13:06

Kirill Belyaev

91
1
4

Poisson equation with Neumann boundary conditions

1 Answers1

Linked