Poisson equation, spectral method and Hermitian symmetry of DFT Let's say I have 2D rectangular grid $N \times M$ with the real values values $f_{nm}$.
I would like to solve 2D Poisson equation
$$\Delta u = f \,\,\, in \,\,\, [0,A] \times [0,B].$$
with periodic boundary conditions
\begin{equation}
\begin{aligned}
u(x,0) &= u(x,B), \quad x \in [0,A]\\
u(0,y) &= u(A,y), \quad y \in [0,B]\\
u_y(x,0) &= u_y(x,B), \quad x \in [0,A]\\
u_x(0,y) &= u_x(A,y), \quad y \in [0,B].\\
\end{aligned}
\end{equation}
using spectral method along the lines of
https://math.stackexchange.com/a/2583914/224869
and
https://stackoverflow.com/questions/42878802/solve-the-poisson-equation-using-fft-with-cuda
There are still issues I don't understand fully:

*

*There is a claim that this is equivalent to solving linear system with the same boundary conditions using five point stencil finite difference scheme for Laplacian. Is there a proof that this is true? I am talking about equivalence up to numerical error for rectangular domain not just about the fact that it solves the same equation.

*In the documentation of cufft, chapter multidimensional transforms, you can see there is some "group symmetry" so that when I know the transform
$$X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i k n/N}$$
for $k \in \{0, \cdots, N-1\}$ I know it any given $k\in \mathbb{N}$. How it comes that when doing transition to the continuous case, I have to use only those $k$ from the sequence, which are closest to $0$? What would fail if I used $k \in \{0, \cdots, N-1\}$? Moreover how to treat case of even $N$, where you have this "unbalanced harmonic"? Can you suggest some source where this transition between DFT and FFT is well covered?

*I would like to use real to complex transform, which is more memory efficient as it uses roughly half of the data. Namely for
$$X_{kl} = \sum_{n=0}^{N-1} \sum_{m=0}^{M-1} x_{nm} e^{-2\pi i (k n/N + l m/M)}$$
it stores only $k \in \{0, \cdots, N-1\}$ and $l \in \{0, \cdots, \lfloor M/2\rfloor\}$. This relies on the Hermitan symmetry $X_{kl} = \overline{X}_{N-k M-l}$. Will I break the symmetry when I start dividing by the terms related to the second derivative in Poisson equation? I would like to operate on that array of $k \in \{0, \cdots, N-1\}$ and $l \in \{0, \cdots, \lfloor M/2\rfloor\}$ and later use so called complex to real transform, which use that hermitian symmetry of input data. I think this is related to the treatment of unbalanced harmonics.

*Is it possible to change the method to compute Poisson equation with different type of boundary conditions? How the spectral method shall be modified when I need to prescribe zero Dirichlet conditions everywhere on the boundary? Is there a chance for generalization for any combination of Dirichlet and Neumann?

 A: I have asked this question some time ago and in the meantime I gained some insights into the problem. Thus I will try to address the questions I had and explain the relationships between continuous problem, discrete problem and FFT algorithm. I would also like to note that I have implemented this in this CUDA program.
First of all, because of the fact that we have just $N \times M$ values of the function $f$ at the gridpoints and we would like to solve continuous equation
$$ \Delta u = f,$$
we would need the mean we understand the objects $f$ and $u$. How $f$ will be defined outside the gridpoints and in which points of the rectangle $[0,A] \times [0,B]$ we actually place gridpoints. I will be using notation, where I identify $N$, $n$ with axis $x$ and $M$, $m$ with axis $y$.
We aim to construct continuous object $f$, which coincides with $f_{nm}$ at the grid points
$$f(n\frac{A}{N}, m \frac{B}{M}) = f_{nm}, \quad n \in \{0...N-1\}, m\in \{0 ... M-1\}.$$
Let's do FFT of $f_{nm}$ to obtain $F_{kl}$, now we have that
$$f_{mn} = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}$$
and thus
$$f(n\frac{A}{N}, m \frac{B}{M}) = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}$$
can be extended to the $f$ defined on the whole domain by the following formula
$$f(x, y) = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k x/A + l y/B)}.$$
Partial answer 2:
This definition of continuous object $f$ is possible but probably not the best one. If we subscribe to this definition, the function $f$ will be real at the gridpoints, but outside of them it will be complex function in general. So we use the trick based on the following identity that holds for the discrete IFFT
$$f_{mn} = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}=\sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{kl} e^{2\pi i (k n/N + l m/M)},$$
where we identify $F_{kl} = F_{k+M l+N}$. The identity holds because $e^{2\pi i} = e^{0} = 1$. For us however is more convenient to use the second sum in the identity to define continuous object $f$. When $N$ and $M$ will be odd, then the function
$$f(x, y) = \sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{kl} e^{2\pi i (k x/A + l y/B)}$$
will also be purely real and can also be represented as a expansion by means of trigonometric functions. This interpolation of our grid by the sum of trigonometric functions is actually the object we would like to use to solve our differential equation. Analogously we define the function $u$
$$u(x, y) = \sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor}  U_{kl} e^{2\pi i (k x/A + l y/B)},$$
and therefore its Laplace will have the form
$$\Delta u(x, y) = -(2 \pi)^2\sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor}  \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl} e^{2\pi i (k x/A + l y/B)}.$$
From the theory about Poisson's equation with periodic boundary conditions we know that $f$ and $u$ have to have the zero mean. This can be achieved only by first ensuring that the input data $f_{nm}$ have zero mean, second by fixing $U_{00}=F_{00}=0$.
Partial answer 1: Note what is now our discrete problem. It is to find a values of $u_{mn}$ on the gridpoints such that the Laplace operator on the $u$ defined as the "Fourier extension" the gridpoint values will coincide with the "Fourier extension" of the object $f_{mn}$. We don't discretize the differential operator on the grid. We rather construct Fourier or trigonometric extension of the functions involved and understand the derivatives in the point sense at the gridpoints as the derivatives of interpolating functions $f$ and $u$. From the restriction of the Poisson equation to the gridpoints we know, that we can find $u_{nm}$ by IFFT of $U_{kl}$, where
$$F_{kl} =  -(2 \pi)^2  \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl}.$$
From how we constructed $f$ and $u$ it is evident that the derivatives are understood as a point derivatives at the gridpoints of the interpolations by trigonometric functions. Thus they will not coincide with any finite difference scheme. The referenced article rather states that "when the problem is discretized by the finite difference (FD) method with a five-point stencil on an (n + 1) × (n + 1) equispaced grid, there is a FFT-based algorithm that computes the values of the solution... " I am not aware of such FFT-based algorithm, the referenced book has to be paid for, but for the algorithm I just have outlined the five-point stencil won't coincide with aforementioned definition of Laplacian on trigonometric extension. It does not exclude that it probably could be shown that these two definitions are in some sense close, e.g. by Taylor expansion.
Partial answer 2 - continue: There is also one more argument, that $k$ and $l$ up to $N-1$ and $M-1$ is not a good idea other that $f$ will be complex. The higher values of $|k|$ are there, the higher frequencies will be present in my interpolation of $f$. And the higher differences might be between point derivatives of $f$ and five-point stencil or whatever reasonable discretization of Laplace operator on my grid.
Partial answer 3: We have now the method how to solve this discrete problem by finding functions $f$, $u$ with the support on $[0,A] \times [0,B] $ with given values of $f$ on the gridpoints so that these two functions satisfy Poisson equation. It can be done for any $o_k \in \mathbb{N}$, $o_l \in \mathbb{N}$ and frequencies $k \in \{ o_k ... N - 1 - o_k\}$,  $l \in \{ o_l ... M - 1 - o_l\}$. We have seen that by setting $o_k=-\lfloor \frac{M-1}{2} \rfloor$ and $o_l=-\lfloor \frac{N-1}{2} \rfloor$ we obtain real $f$ in odd case. In even case however, there will still be the term correspondig to $k=N/2$
$$ \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(N/2)l} e^{2\pi i (N/2 x/A + l y/B)}$$
which needs to be symmetrized to obtain real function $f$ in the way
$$  \frac{1}{2} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(-N/2)l} e^{2\pi i (-N/2 x/A + l y/B)} + \frac{1}{2} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(N/2)l} e^{2\pi i (N/2 x/A + l y/B)}$$
In fact it would need to be symmetrized in $l$ as well. By doing so we have found real $f$ even in unbalanced case. And if moreover we use $k$ and $l$ that are centered around $0$ from the formula
$$F_{kl} =  -(2 \pi)^2  \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl}$$
we also see that $U_{-kl}=\bar{U}_{kl}$ because of $F_{-kl}=\bar{F}_{kl}$ and because of the choice of frequencies I does not break hermitan symmetry.
Partial answer 4:
There is a theory regarding coercive operators, that for the problem with Dirichlet boundary conditions on the part of the boundary there is unique weak solution to the Poisson equation, for sources see this answer. It is natural to ask if we can use spectral methods to solve Poisson equation on square with Neumann or Dirichlet boundary conditions instead of periodic boundary conditions. I am not sure about general answer yet. It seems however that the simplicity of the spectral method stems from the fact, that functions in the Fourier decompositions are eigenfunctions of Laplace operator with periodic boundary conditions. Whenever we lose this, the method gets more ugly.
It is apparently possible to use something called odd or even extension to solve problem with homogeneous Dirichlet or homogeneous Neumann. Key observation here is that any sine function in FFT on $[0,2A]$ is zero at $0$, $A$ and $2A$ and any cosine function in FFT on $[0,2A]$ will have zero derivative at $0$, $A$ and $2A$. And if we construct odd function, it will have zero even components in Fourier transform and vice versa.
So let's say we would like to have homogeneous Dirichlet and we have
$$f(x,y), x \in [0,A], y\in [0, B]$$ and we extend it to $[0,2A] \times [0, 2B]$ such that
$$f(x,y) = \begin{cases}
f(x,y), & x \in [0,A], y \in [0,B],\\
-f(2A-x,y), & x \in [A,2A], y \in [0,B],\\
-f(x,2A-y), & x \in [0,A], y \in [B,2B],\\
f(2A-x,2A-y) & x \in [A,2A], y \in [B,2B].
\end{cases}$$
Now we solve the original problem $\Delta u = f$ with periodic boundary conditions on bigger domain $[0,2A] \times [0, 2B]$. The restriction of the solution to the $[0,A] \times [0, B]$ will be solution to the original problem with Dirichlet boundary conditions. See also Figure 3 of this paper. By solving four times bigger problem, we managed to use this type of spectral method to solve the following boundary value problems
\begin{equation}
\begin{aligned}
u(x,0) &= u(x,B) = 0, \quad x \in [0,A],\\
u(0,y) &= u(A,y)=0,  \quad y \in [0,B].
\end{aligned}
\end{equation}
or
\begin{equation}
\begin{aligned}
u(x,0) &= u(x,B) = 0, \quad x \in [0,A],\\
u_x(0,y) &= u_x(A,y) = 0, \quad y \in [0,B].\\
\end{aligned}
\end{equation}
or
\begin{equation}
\begin{aligned}
u_y(x,0) &= u_y(x,B)=0, \quad x \in [0,A],\\
u_x(0,y) &= u_x(A,y)=0, \quad y \in [0,B].\\
\end{aligned}
\end{equation}
Of course that raises further question if there is a way to extend the set of problems even further. And if if it is necessary to solve four times bigger problem to have solution with some kind of normal boundary conditions or if we can computationally exploit the symmetry of the problem.
Apparently when homogeneous Dirichlet conditions are imposed and I need other type of Dirichlet conditions, it is sufficient to compute Laplace equation solution on the square with those Dirichlet conditions and sum the two solutions. This way the method can be generalized for any type of Dirichlet boundary conditions.
Note also when doing the odd or even extension we have to pay attention how FFT and IFFT of sequence [0,n-1] is defined with respect to the domain [0,A]. The last point $n-1$ is mapped not to A but to $A-A/N$ so if we do naïve extension by simple mirroring we obtain lot of aliasing between FFT and FT frequencies. From the perspective of keeping radix-2 algorithm it might be much faster than using different than $2N$ grid.
