Role of non-hermitian coefficient matrices in the discretization of self-adjoint operators

Question

What is the role of non-symmetric coefficient matrices in the solution of partial differential equations with self-adjoint operators?

Particularly, I'm thinking about time-propagation of a linear equation

$$i\partial_t v =M v$$ where $M$ is self-adjoint. (Or equivalently, without the complex unit $i$ and then $M$ is skew-hermitian, or in the real case, symmetric or skew-symmetric),

Discretization of this equation on an equidistant grid, together with the usual finite-difference formulas for derivatives, leads to a hermitian coefficient matrix. The same is true for a spectral expansion. The arising matrices have nice properties: one can diagonalize them and get real eigenvalues and orthogonal eigenvectors, one can apply linear algebra routines tailored to symmetric matrices, and so on.

On the other hand, there are several successful methods (finite differences on non-equally spaced grids, one-sided finite differences, Chebyshev and other pseudospectral methods, etc.) which allow for a more accurate expansion (e.g. finer grid in important regions), but usually lead to non-hermitian coefficient matrices.

What is preferable? Should I keep the symmetry as long as possible? Or can one rely as well on non-symmetric matrices?

This is a question which is puzzling me for a long time and where I have not found a concrete analysis yet. However, almost any reference seems to give statements such as "symmetry in the coefficient matrix is preferable" due to "stability" or whatever.

Questions (more detailed):

1. What is the theoretical impact of this symmetry breaking? Problems I see: eigenvalues can become complex, so there is no guarantee that variational optimization still holds. This also implies that conservation laws do not hold anymore. Is this a problem in practice, or does a good choice of non-symmetric discretization (numerically) avoids these problems?

2. What is the drawback in using non-symmetric coefficient matrices? (Besides the obvious one that no Linear Algebra routines specifically designed to symmetric matrices can be used).

3. If there should be a disadvantage by using non-symmetric matrices: Is it reasonable to introduce symmetry artificially? E.g. by symmetrizing $[\tilde M_{ij} = 1/2(M_{ij} + M_{ji})]$, setting up the second derivative matrix (which is often the source of asymmetry) in terms of the first $(D_2 = D_1^T D_1)$ or by more sophisticated constructions like variable changes, etc.?

4. Summarizing the previous points: are there any references which consider this topic (under whatever aspect)?

Answer 1

I'll answer why skew-symmetric matrices are sometimes preferable.

The first thing to note is that this is usually problem dependent. But for a general 1D PDE like

$\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0$

skew-symmetry has a distinct advantage. Suppose we discretize this equation, then we'll get

$u_t + Du = 0$

Now suppose that D is skew-symmetric, then we have

$D^T = -D$

Let's now find out what the derivative of $u^2$ is

$(u^Tu)_t = (u^T)_tu + u^Tu_t $

Substituting from our expression for $u_t$ we get,

$(u^Tu)_t = -(Du)^Tu - u^TDu $

$(u^Tu)_t = -u^TD^Tu - u^TDu $

$(u^Tu)_t = -u^T(D^T + D)u $

Recall the skew-symmetry matrix property, $D^T = -D$

So,

$(u^Tu)_t = 0$

There are a couple of things to note. First $(u^Tu)_t $ is what we usually call kinetic energy (with a factor of 2).

Second, the original equation had a divergence term,

$ \frac{\partial f}{\partial x} $

that means that if we do a volume integral, any energy can only come in through the boundaries. So, $(u^Tu)_t = 0$ should be satisfied in general unless there's a dissipation term in the PDE. And that is exactly what the skew-symmetry operator is doing.

I'll link to a presentation that I found useful.

http://calcul.math.cnrs.fr/Documents/Ecoles/CEMRACS2012/Julius_Reiss.pdf

Again, all this depends on individual problems. A different PDE may not like this method. A different operator matrix may have some other property that is preferable.

For example, in

https://arxiv.org/pdf/1605.09763v1.pdf

The convection operator that is,

$u.\nabla u $

is transformed to the skew-symmetric form

$ \frac{1}{2}(u.\nabla u + (\nabla.u)u) $

There is in fact a lot of literature in which a similar transformation is done. The reason for doing it is because often a uniform mesh is not available so there is less control available over the discretized operator, so instead the PDE term is itself transformed to enforce skew-symmetry.

Another common technique is to make the finite difference operator into a summation by parts operator which is essentially converting it into skew-symmetric. This is a paper that shows the technique.

http://link.springer.com/article/10.1007/s10543-015-0599-0

But obviously this is only useful when finite difference is an option.

So, in essence, more and more studies are being done to create schemes so that desired properties be satisfied.