3
$\begingroup$

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: \begin{equation} Y_{k}\triangleq\sum_{n=0}^{N-1}X_{n}e^{-2i\pi nk/N}, \end{equation} and the following definition for the backward (or inverse) transform: \begin{equation} Y_{k}=\sum_{n=0}^{N-1}X_{n}e^{2i\pi nk/N}. \end{equation}

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.

  • 0
    @Fabian corrected. Thanks.2012-01-11

1 Answers 1

4

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}. $

  • 0
    @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$.2012-01-11