Self-adjoint and nonpositive differential operators

Question

I recently tumbled over a statement in a geophysics paper (PDF here). They have a wave equation which they formulate as

$$ \frac{1}{v_0}\frac{\partial^2}{\partial t^2} \begin{pmatrix}p \\ r\end{pmatrix} = \begin{pmatrix} 1+2\epsilon&\sqrt{1+2\delta}\\\sqrt{1+2\delta} &1\end{pmatrix} \begin{pmatrix}G_{\bar x \bar x}+G_{\bar y \bar y}&0\\0&G_{\bar z \bar z} \end{pmatrix} \begin{pmatrix} p\\r \end{pmatrix}\tag{20} $$

and they claim that

To achieve stability, the rotated differential operators $G_{\bar x \bar x}$, $G_{\bar y \bar y}$, and $G_{\bar z \bar z}$ should be self-adjoint and nonpositive definite as are the second-order derivative operators ($\tfrac{\partial^2}{\partial x^2}$, $\tfrac{\partial^2}{\partial y^2}$ and $\tfrac{\partial^2}{\partial z^2}$).

(see eq. 14 and statement under eq. 20). They also claim that self-adjointness and nonpositiveness of the differential operators of the wave equation are necessary to conserve the energy in this system, and that if they were not self-adjoint, numerical instabilities occur.

We have solved the problem by introducing the self-adjointness to the operator matrices in equation 20 to make sure that energy is conserved during the wave propagation to avoid amplitude blowup in the modeling.

In this case the wave equation consists of two coupled elliptical PDEs. What happens in general, when some operators are not bounded and linear?

Unfortunately I don't have the mathematical background to understand this statement.

Answer 1

This is much to do with the possible eigenvalues of the operators. Normal operators on a Hilbert space are closely analogous to complex numbers, with the adjoint taking the role of the conjugate; these relations are typically inherited directly to the operator's eigenvalues. Thus, if a linear operator $L$ has an eigenfunction $f$ with eigenvalue $\lambda$, $$Lf=\lambda f,$$ then saying "$L$ is self-adjoint" means that $L^\dagger=L$ which translates to $\lambda^\ast=\lambda$, i.e. that $\lambda$ be real. Similarly, $L$ being nonpositive implies that $\lambda\leq0$.

In your case, the behaviour can be reduced to an equation of the form $$ \frac{\partial^2 p}{\partial t^2}(x,t)=\hat Lp(x,t), $$ where $\hat L$ is some differential operator. In general the solution will not be of this form, but you can take a first stab of the problem by inserting in an eigenfunction of the differential operator for the spatial dependence. That is, you use the trial solution $p(x,t)=p_0(x)T(t)$, where $\hat Lp_0=\lambda p_0$. This hugely simplifies the time-propagation equation, which reduces to the solvable form $$ \frac{\partial^2 }{\partial t^2}T=\lambda T. $$

While the solutions of this equation are formally all the same (i.e. $T(t)=T_+e^{\sqrt{\lambda}t}+T_-e^{-\sqrt{\lambda}t}$) regardless of what $\lambda$ is, the behaviour will be very different and depend, sometimes sensitively, on $\lambda$:

• If $\lambda>0$, then at least one of the exponentials $e^{\pm\sqrt{\lambda}t}$ will have a blow-up.
• If $\lambda$ has an imaginary part, however small, then one of the two square roots $\pm\sqrt{\lambda}$ will have a positive real part, and the corresponding contribution to $T(t)$ will oscillate at a blowing-up amplitude.
• If $\lambda$ is negative or zero, then both roots $\pm\sqrt{\lambda}$ will be imaginary or zero, and both exponentials will be completely oscillatory and have bounded amplitude for all time.

It is clear that only the third case is consistent with conservation of energy. In terms of the differential operator, it corresponds to a condition of self-adjointness (i.e. $\lambda\in\mathbb R$) and non-negativity of the operator.