Solve Poisson Equation Using FFT

Question

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?

Answer 1

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:

1. Sample $f(x,y)$ at the points $(m \Delta x, n\Delta y)$;
2. Perform the FFT of $f_{mn}$;
3. 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}$;
4. 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.

Answer 2

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.