# FFT for solving Poisson equation in 3 dimensions with periodic boundary conditions

In electrostatics, the laplacian of the electrostatic potential $\Delta V(\mathbf{r})$ arising from a charge distribution $\rho(\mathbf{r})$ is

$$\Delta V(\mathbf{r})=-\frac{\rho(\mathbf{r})}{\epsilon_0},$$

where $\epsilon_0$ is the vacuum permitivity.

I would like to use Fast Fourier Transforms (FFT) to solve this PDE using the FFTW library. FFTW computes $Y$, the Discrete Fourier Transform (DFT) of a 1D complex array $X$ of size $N$ with this definition for the forward transform: $$Y_{k}\triangleq\sum_{n=0}^{N-1}X_{n}e^{-2i\pi nk/N},$$ and the following definition for the backward (or inverse) transform: $$Y_{k}=\sum_{n=0}^{N-1}X_{n}e^{2i\pi nk/N}.$$

As far as I understand, the sampling interval in these definitions is equal to 1. Also, $FFT^{-1}(FFT(X))=N\cdot X$, which is not usual. The manual says this definition is not unitary and that normalization should be taken into consideration by the user.

My application is in 3 dimensions, with regular sampling in each direction, but the sampling frequency is not the same in all directions ($L_1$, $L_2$, $L_3 \neq 1$).

My question : could you help me get an expression of $\Delta V$ in k-space ? With a simplest definition of FFTs, it would be for instance $-\mathbf{k}^2 V(\mathbf{k})$. Here ... I can't get it.

Sorry if the question is a bit messy. I'm not a mathematician and this is quite new for me.

-
I think $\Delta V$ should be $-k^2 V$ in $k$-space. – Fabian Jan 11 '12 at 16:35
@Fabian corrected. Thanks. – max Jan 11 '12 at 17:36

You define $L_1, L_2, L_3$ to be the sampling interval in the 3 directions. Thus a point $\mathbf{r}$ in real space is given by $$\mathbf{r} = (L_1 n_1, L_2 n_2, L_3 n_3)$$ where $n_j \in \{0,\dots, N-1\}$.

Let us understand how $\partial V(\mathbf r)/\partial r_1$ is represented in Fourier space. As you sample your space, the derivative has to be replaced by a finite difference. So we have $$\frac{\partial V(\mathbf r)}{\partial r_1} \approx \frac{V(\mathbf{r} + L_1 \mathbf{e}_1) - V(\mathbf{r})}{L_1}= \frac{V_{n_1 +1, n_2, n_3} - V_{n_1-1,n_2,n_3}}{2L_1}$$ where I introduced the notation $V_{n_1,n_2,n_3} = V(\mathbf{r}=(L_1 n_1, L_2 n_2, L_3 n_3))$.

Given the definition of the Fourier transform in the question, you have $$\hat{V}_{k_1,k_2,k_3} = \sum_{n_1,n_2,n_3} V_{n_1,n_2,n_3} \exp\left(-2\pi i \sum_{j=1}^3 n_j k_j/N_j\right). \qquad\qquad(1)$$ Then $$V_{n_1,n_2,n_3} = \frac{1}{N_1 N_2 N_3} \sum_{k_1,k_2,k_3} \hat V_{k_1,k_2,k_3} \exp\left(2\pi i \sum_{j=1}^3 n_j k_j/N_j\right); \qquad\qquad(2)$$ Eq. (1) is the definition of $\hat{V}$ and $\hat{V}$ will be what the library will hopefully calculate for you. Eq. (2) then just follows from (1) --- and sadly does not coincide with the backtransform which the library offers.

Plugging (2) into the expression for $\partial V(\mathbf r)/\partial r_1$ yields $$\frac{\partial V(\mathbf r)}{\partial r_1} = \frac{1}{N_1 N_2 N_3} \sum_{k_1,k_2,k_3} \frac{i\sin(2\pi k_1)}{L_1} \hat V_{k_1,k_2,k_3} \exp\left(2\pi i \sum_{j=1}^3 n_j k_j/N_j\right).$$ Thus $\partial/\partial k_1$ is represented in the Fourier domain as $i\sin(2\pi k_1)/L_1$. Analogue, you can easily show that $$\frac\partial{\partial k_j} \mapsto \frac{i\sin(2\pi k_j)}{L_j}.$$

Using these results, we obtain $$\Delta \mapsto - \sum_{j=1}^3\frac{\sin(2\pi k_j)^2}{L_j^2}.$$

-
That's very nice I'll now be able to demonstrate the Froueri transform (whatever the definition) of $\Delta$. Thank you @Fabian. What is unclear to me is considering my practical problem and the actual definition of backward FFT in FFTW, what should I divide the forward FFT of $-\rho(\mathbf{r})/\epsilon_0$ by in order to calculate $V(\mathbf{r})$? – max Jan 11 '12 at 17:48
@MaximilienLevesque: keep the forward transform as defined in the package, so just solve $[\sum_j \sin^2(2\pi k_j)/L_j^2] \hat{V} = \hat{\rho}/\epsilon_0$. Just make sure that you include the factor $1/N_1 N_2 N_3$ when (at the end) transforming back to $V$. – Fabian Jan 11 '12 at 18:07