4
$\begingroup$

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?

||cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Sumedh Joshi Jun 16 '14 at 18:51
  • $\begingroup$ 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. $\endgroup$ – Grieverheart Jun 16 '14 at 19:17
  • 2
    $\begingroup$ 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. $\endgroup$ – Sumedh Joshi Jun 16 '14 at 19:31
  • $\begingroup$ 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. $\endgroup$ – Grieverheart Jun 16 '14 at 20:21
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter Brune Jun 16 '14 at 22:11
4
$\begingroup$

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.

||cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.