Solve Poisson Equation Using FFT I am trying to solve Poisson equation using FFT. The issue appears at wavenumber $k = 0$ when I want to get inverse Laplacian  which means division by zero.
We have 
${\nabla ^2}\phi  = f$ 
Taking FFT from both side we get:
$-k^2\hat\phi = \hat f $
or 
$\hat\phi = \frac{\hat f}{-k^2} $
Assuming that we want to solve this equation in periodic domain and using DFT using FFTW package, at $k=0$ we have a division by zero. Anybody knows how to deal with this singularity?
 A: THE PROBLEM
Let me first restate the problem in 2 dimensions:
Solve the Poisson equation
$\nabla^2\phi(x,y)=f(x,y)\;\;\;(x,y)\in[0, L_x]\times[0, L_y]$
subject to periodic boundary conditions, namely
$\phi(x,0)=\phi(x,L_y) \;\;\; x\in[0, L_x]$
$\phi(0,y)=\phi(L_x,y) \;\;\; y\in[0, L_y]$
$\frac{\partial \phi}{\partial x}(0, y) = \frac{\partial \phi}{\partial x}(L_x, y)\;\;\; y\in[0, L_y]$
$\frac{\partial \phi}{\partial y}(x, 0) = \frac{\partial \phi}{\partial y}(x, L_y)\;\;\; x\in[0, L_x]$
THE FFT-BASED SOLUTION APPROACH
Let us expand $\phi(x,y)$ in the Fourier series
$$\phi(x,y)=\sum_{p,q=-\infty}^{\infty}\hat{\phi}_{pq}e^{-j2\pi p\frac{x}{L_x}}e^{-j2\pi q\frac{y}{L_y}}$$
where the $\hat{\phi}_{pq}$'s are the expansion coefficients, and let us do the same for $f(x,y)$ with expansion coefficients $\hat{f}_{pq}$. On substituting the Fourier series expansion into the Poisson equation we have
$$\left[\left(-j2\pi \frac{p}{L_x}\right)^2+\left(-j2\pi \frac{q}{L_y}\right)^2\right]\hat{\phi}_{pq}=\hat{f}_{pq}$$
Solving the Poisson equation amounts at solving such an equation for the $\hat{\phi}_{pq}$'s.
On specializing the Fourier series of $\phi(x,y)$ at the points $(x_m, y_n) = (m \Delta x, n \Delta y)$, with $m = 0,1,\ldots,M$ and $n=0,1,\ldots,N$, we have
$$\phi_{mn}=\phi(m\Delta x,n\Delta y)=\sum_{p,q=-\infty}^{\infty}\hat{\phi}_{pq}e^{-j2\pi p\frac{m\Delta x}{L_x}}e^{-j2\pi q\frac{n\Delta y}{L_y}}$$
Finally, if we choose
$\Delta x = \frac{L_x}{M}$ and $\Delta y = \frac{L_y}{N}$ and truncate the Fourier series to $(M+1)\times(N+1)$ terms, then
$$\phi_{mn}=\sum_{p=0}^{M}\sum_{q=0}^{N}\hat{\phi}_{pq}e^{-j2\pi p\frac{m\Delta x}{L_x}}e^{-j2\pi q\frac{n\Delta y}{L_y}}$$
which is the expression of a Discrete Fourier Transform (DFT). The same holds true for $f(x,y)$, i.e.
$$f_{mn}=\sum_{p=0}^{M}\sum_{q=0}^{N}\hat{f}_{pq}e^{-j2\pi p\frac{m\Delta x}{L_x}}e^{-j2\pi q\frac{n\Delta y}{L_y}}$$
where $f_{mn}=f(m \Delta x, n\Delta y)$. The approach is then the following:

*

*Sample $f(x,y)$ at the points $(m \Delta x, n\Delta y)$;

*Perform the FFT of $f_{mn}$;

*Divide $\hat{f}_{pq}$ by $\left[\left(-j2\pi \frac{p}{L_x}\right)^2+\left(-j2\pi \frac{q}{L_y}\right)^2\right]$ to compute $\hat{\phi}_{pq}$;

*Perform the IFFT of $\hat{\phi}_{pq}$.

THE ISSUE
The above equation can be solved for any $(p,q)$, except for $p=q=0$, giving an ambiguity for the determination of the zero frequency Fourier coefficient $\hat{\phi}_{00}$.
THE SOLUTION
Seeing the issue from a different perspective, the boundary conditions define the solution of the above Poisson equation only
up to an additive constant: if $\phi(x,y)$ solves this equation and obeys the above periodic boundary conditions, then so does $\phi(x,y)+c$ for any real constant $c$. Accordingly, we need another condition to pin the constant $c$ down. The standard stipulation is that
$$\int_{0}^{L_x}\int_{0}^{L_y}\phi(x,y)dxdy = 0$$
Note that, on using the periodic boundary conditions and on integrating the Poisson equation, it can be easily shown that the forcing term $f(x,y)$ must also obey the condition
$$\int_{0}^{L_x}\int_{0}^{L_y}f(x,y)dxdy = 0$$
REFERENCE
Arieh Iserles, A first course in the numerical analysis of differential equations, Section 10.4, Cambridge Texts in Applied Mathematics.
A: This is conceptually similar the integrating constants that show up when you are solving a differential equation by other methods. Usually these integrating constants are decided by your boundary conditions.
Instead of doing any division you can simply rewrite it as a least-norm problem: $$\|-k^2\hat \phi +\hat f\|_2^2+\text{any more terms you may want}$$
If your $\hat \phi,\hat f$ are stored in vectors the $k^2$ in the expression above will be  a diagonal weight matrix which multiplies the $\hat \phi$ vector.

Oops sorry I forgot where the FFTW comes into place, you can use it to calculate the matrix-vector product, $F$ below is the Fourier transform matrix.
$$\|-k^2 \hat \phi + F f\|_2^2$$
So $Ff$ means multiply $f$ vector by $F$ and it will be the output of calling your FFTW library with $f$ vector as the input.
And for other terms you have it might be $F^{-1}$ you want to multiply with but that is simply the corresponding IFFT routine in your library.
