If one has the Poisson problem (assume $\int_{\Omega} f = 0$ and $\int_{\Omega} u = 0$):
\begin{alignat}{3} \Delta u(x) &= f(x), &\quad&x\in\Omega \\ \partial_nu(x) &= 0, &\quad&x\in\partial\Omega, \end{alignat}
then a standard approach is to discretise it using finite differences. For example in 1D we have:
\begin{equation} (\Delta u)_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}. \end{equation}
Enforcing the Neumann boundary conditions one ends up with a system of the form:
\begin{equation} \begin{bmatrix} -1 & 1 & 0 & 0 & 0 & \ldots & 0\\ 1 & -2 & 1 & 0 & 0 & \ldots & 0\\ 0 & 1 & -2 & 1 & 0 & \ldots & 0\\ &&&\ldots&& \\ 0 & \ldots & 0 & 1 & -2 & 1 & 0\\ 0 & \ldots & 0 & 0 & 1 & -2 & 1\\ 0 & \ldots & 0 & 0 & 0 & 1 & -1\\ \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ \ldots \\ u_n \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \\ \ldots \\ f_n \end{bmatrix} \end{equation}
I have seen (e.g. https://elonen.iki.fi/code/misc-notes/neumann-cosine/) the DCT being suggested as a solution for this problem. If I term the above matrix $L$, then as I understand it, the DCT matrix $V^T$ is supposed to diagonalize $L = V\Lambda V^T:$ \begin{align} Lu&=f \\ V^TLu&=V^Tf \\ V^TV\Lambda V^T u &= V^Tf \\ V^T u &= \Lambda^{-1}V^T f \\ u &= V \Lambda^{-1} V^T f. \end{align} I know for instance that the DFT diagonalizes circulant matrices, however I wasn't aware that the DCT diagonalizes matrices of the above form. I found the following work on the structure of matrices that the discrete cosine and sine transforms diagonalize: "Matrices diagonalized by the discrete cosine and discrete sine transforms". I decided to verify that this is indeed the case. I took as an example: \begin{equation} L = \begin{bmatrix} -1 & 1 & 0 \\ 1 & -2 & 1 \\ 0 & 1 & -1 \end{bmatrix} \end{equation} and used the form described in the paper to try and match it: \begin{equation} \begin{bmatrix} t_0 & \sqrt{2}t_1 & t_2 \\ \sqrt{2}t_1 & t_0+t_2 & \sqrt{2}t_1 \\ t_2 & \sqrt{2}t_1 & t_0 \end{bmatrix} \, \textrm{or} \, \begin{bmatrix} t_0 & \sqrt{2}t_1 & \sqrt{2}t_2 \\ \sqrt{2}t_1 & t_0+t_2 & t_1+t_2 \\ \sqrt{2}t_2 & t_1+t_2 & t_0+t_1 \end{bmatrix}. \end{equation}
However neither of the described forms (corresponding to different types of DCT) can produce the above matrix, as $t_0=-1, t_2 = 0 \implies t_0+t_2 = -1 \ne -2$. Am I misunderstanding something?
In fact, since the DCT acts as if the signal had been mirrored, I would expect that it should be able to diagonalize any circulant matrix $C\in\mathbb{R}^{(n+2m)\times(n+2m)}, \, m\leq n$ that has been restricted as follows:
\begin{equation} \begin{bmatrix} \boldsymbol{0}_{n\times m} & \boldsymbol{I}_{n} & \boldsymbol{0}_{n\times m} \end{bmatrix} C\begin{bmatrix} \boldsymbol{A}_{m\times n} \\ \boldsymbol{I}_{n} \\ \boldsymbol{B}_{m\times n} \end{bmatrix}, \end{equation}
where $\boldsymbol{I}_n$ is the $n\times n$ identity matrix and $\boldsymbol{A}_{m \times n}$ is the exchange matrix (ones on the antidiagonal) with its top $(n-m)$ rows removed, while $\boldsymbol{B}_{m \times n}$ is the exchange matrix with it bottom $(n-m)$ rows removed. Essentially the above acts as if the vector to be multiplied on the right was reflected on the boundaries for $m$ elements, e.g. $v = (v_1, \ldots, v_n)$ results in: \begin{equation} u = (v_m, \ldots, v_1, v_1, \ldots, v_n, v_n, \ldots, v_{n-m+1}), \end{equation} then the full $C$ is applied, and finally the extra elements are thrown away.