Discrete sine transform Poisson solver (source code)

= Discrete sine transform Poisson solver
{title2=$\widehat u_{k\ell}=\widehat f_{k\ell}/\lambda_{k\ell}$}

For $-\Delta u=f$ on a square with homogeneous <Dirichlet boundary conditions>, the five-point <finite difference> operator on a mesh $h=1/N$ is diagonalized by $\sin(\pi ki/N)\sin(\pi\ell j/N)$. Its eigenvalues are $\lambda_{k\ell}=4h^{-2}[\sin^2(\pi k/(2N))+\sin^2(\pi\ell/(2N))]$. A two-dimensional <discrete sine transform>, division by these nonzero eigenvalues and the inverse transform give the solution. Odd extension realizes each one-dimensional transform by a <Fast Fourier transform> of length $2N$, making the total cost $O(N^2\log N)$ and storage $O(N^2)$. These finite-difference eigenvalues should not be confused with the continuous spectral eigenvalues $\pi^2(k^2+\ell^2)$.