constructing a symmetric matrix for finite difference

I come across the following operator in a paper

$$\mathcal{I}\psi = \psi_{xxxx} + (r~\psi_x)_x$$,

where $$\psi=\psi(x)$$ and $$r=r(x)$$. Periodic boundary condition is employed. It claims that the operator $$\mathcal{I}$$ is symmetric. Suppose I use finite difference to discretize it, then the first part $$\psi_{xxxx}$$ is obviously symmetric. As an example, for a 6x6 finite difference matrix with 2nd order accuracy and periodic BC, the finite difference matrix (assume step size dx=1 for simplicity) is just

$$D_4 = \begin{pmatrix} 6 & -4 & 1 & 0 & 1 & -4\\ -4 & 6 & -4 & 1 & 0 & 1\\ 1 & -4 & 6 & -4 & 1 & 0\\ 0 & 1 & -4 & 6 & -4 & 1\\ 1 & 0 & 1 & -4 & 6 & -4\\ -4 & 1 & 0 & 1 & -4 & 6 \end{pmatrix}$$

The problem comes from the second part $$(r\psi_x)_x$$. This operator is self-adjoint (for periodic BC) as it can be shown that $$\int \phi(r\psi_x)_x dx = \int \psi(r\phi_x)_x dx$$ via integration by part twice, and from what I understand (correct me if I'm wrong) this implies the matrix representing the operator is symmetric (or hermitian if it is complex, just recalling what I learn from my physics QM class). But how do I construct the symmetric matrix representing this part? I thought about expanding it to give $$(r\psi_x)_x = r_x\psi_x + r\psi_xx$$, so that the matrix is given by diagonal(D1*r)*D1 + diagonal(r)*D2. D1 and D2 are the finite difference matrix for first and second derivative with periodic BC analogous to the one above, and are given below for 2nd-order accuracy and size 5x5 as examples:

$$D_1 = \begin{pmatrix} 0 & 1/2 & 0 & 0 & -1/2\\ -1/2 & 0 & 1/2 & 0 & 0\\ 0 & -1/2 & 0 & 1/2 & 0\\ 0 & 0 & -1/2 & 0 & 1/2\\ 1/2 & 0 & 0 & -1/2 & 0\\ \end{pmatrix}$$

$$D_2 = \begin{pmatrix} -2 & 1 & 0 & 0 & 1\\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 1 & 0 & 0 & 1 & -2\\ \end{pmatrix}$$

diagonal(r) gives a diagonal matrix where the diagonal is $$r(x)$$ represented by the vector r. * is just matrix multiplication. D2 is symmetric, but D1 is not. Multiplying D2 with the variable coefficient r also makes it asymmetric. What have I done wrong? Or are there something I misunderstand?

• $$r_x\cdot\psi_x = D_1diagonal(r)D_1diagonal(\psi) = diagonal(r)D_1^TD_1diagonal(\psi) = -diagonal(r)*D_1^2diagonal(\psi)$$ because $D_1$ is skew symmetric it squared is symmetric therefore the operator is symmetric
– EMP
May 2, 2022 at 20:53
• @EMP: wait, a diagonal matrix (which is not constant) destroys the symmetry in general. So even if the second derivative part is symmetric, multiplying it by diag(r) yields a non-symmetric matrix in general. May 2, 2022 at 23:00
• Periodic boundary conditions should directly trigger a Fourier transformation reflex. The derivatives become diagonal then, and the product will become a convolution. May 2, 2022 at 23:25
• The diagonal matrix will ruin the symmetry, but he wants the operator on the variables to be symmetric, yes? so it doesn't matter if $D_1 diag(r)$ is going to ruin the symmetry as long as the matrix $D_1^T D1$ is.
– EMP
May 3, 2022 at 0:21
• @EMP: i don't get that. In Quantum Mechanics, you often want to pass the system matrix to an eigenvalue solver to get the eigenenergies, for example. There it matters whether your whole matrix is symmetric (and not just the indivual factors). May 3, 2022 at 17:50

Here is an approach to achieve a symmetric finite difference operator for a space-dependent diffusion term $$(r \psi_x)_x$$. The trick is to apply central difference discretizations on first derivatives, but using grid spacings of $$\Delta x/2$$:

\begin{align*} (r \psi_x)_x(x_i) &= \frac{r(x_{i+1/2}) \psi_x(x_{i+1/2}) - r(x_{i-1/2}) \psi_x(x_{i-1/2})}{\Delta x} + \mathcal{O}(\Delta x^2) \\ &= \frac{r(x_{i+1/2}) \frac{\psi(x_{i+1}) - \psi(x_i)}{\Delta x} - r(x_{i-1/2}) \frac{\psi(x_{i}) - \psi(x_{i-1})}{\Delta x}}{\Delta x} + \mathcal{O}(\Delta x^2) \\ &= \frac{r(x_{i-1/2}) \psi(x_{i-1}) - (r(x_{i-1/2}) + r(x_{i+1/2})) \psi(x_i) + r(x_{i+1/2}) \psi(x_{i+1})}{\Delta x^2} + \mathcal{O}(\Delta x^2) \end{align*}

Despite using the grid spacing of $$\Delta x/2$$, the cell-centered values of $$\psi$$ cancel, and we only need its value at the standard stencil points $$[x_{i-1}, x_i, x_{i+1}]$$. In this form, $$r$$ is evaluated at cell-centered values. If it is not convenient to get cell-centered values, they can be approximated in a number of ways (see this, for example).

Here is a $$5\times 5$$ example of the finite difference operation

$$\frac{1}{\Delta x^2} \begin{bmatrix} -r_{9/2} - r_{1/2} & r_{1/2} & 0 & 0 & r_{9/2} \\ r_{1/2} & -r_{1/2}-r_{3/2} & r_{3/2} & 0 & 0 \\ 0 & r_{3/2} & -r_{3/2}-r_{5/2} & r_{5/2} & 0 \\ 0 & 0 & r_{5/2} & -r_{5/2}-r_{7/2} & -r_{9/2} \\ r_{9/2} & 0 & 0 & r_{7/2} & r_{7/2}-r_{9/2} \end{bmatrix} \begin{bmatrix} \psi_0 \\ \psi_1 \\ \psi_2 \\ \psi_3 \\ \psi_4 \end{bmatrix}$$

• are there any reasons for using grid spacing of $\Delta x/2$ instead of $\Delta x$? I think the latter also gives a symmetric matrix, except the three nonzero main diagonals are separated by one diagonal of zeroes and the denominator has a factor of 4? May 2, 2022 at 21:32
• The grid spacing of $\Delta x/2$ ensures the resulting finite difference is based on the standard grid points $[x_{i-1}, x_i, x_{i+1}]$. You could follow the same procedures with a spacing of $\Delta x$, but this would have a somewhat non-standard stencil and a larger error. May 3, 2022 at 2:52
• In the second step, using central difference to replace $\psi_x$, the error is $O(\Delta x^2)$, while there is a $\Delta x$ in the denominator, the error approximation should be $O(\Delta x)$. Feb 14 at 1:29
• Yes but the central differences applied in the second step have leading error terms which cancel. Feb 14 at 4:47