constructing a symmetric matrix for finite difference

Question

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?

Answer 1

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