Self-adjoint differential operators I'm having a hard time understanding the deal with self-adjoint differential opertors  used to solve a set of two coupled 2nd order PDEs.
The thing is, that the solution of the PDEs becomes numerically unstable and I've heared that this is due to the fact, that the used operators were not self-adjoint and the energy is not preserved in this case.
The two coupled 2nd order PDEs are:
$$\frac{\partial ^2p}{\partial t^2}=V_{px}^2 {H_2} p + \alpha V_{pz}^2 {H_1} q + V_{sz}^2{H_1}(p - \alpha q) + S\tag{1}$$
$$\frac{\partial ^2q}{\partial t^2}=\frac{V_{pn}^2}{\alpha}{H_2} p + V_{pz}^2 {H_1} q - V_{sz}^2{H_2} \left(\frac{1}{\alpha}p - q \right) + S\tag{2}$$
where p is the pressure wave field and q is an auxiliary wave field, $S$ is the Source term $V_{px}$ and $V_{sz}$ are seismic velocities into the x - or z-direction respectively, $\alpha = 1$ and $H_1$ and $H_2$ are the rotated differential operators:
\begin{eqnarray}
{H_1}& =& \sin ^2 \theta \cos ^2 \phi \frac{\partial ^2}{\partial x ^2}+\sin ^2 \theta \sin ^2 \phi \frac{\partial ^2}{\partial y ^2} + \cos ^2 \theta \frac{\partial ^2}{\partial z ^2}+\\ &&\sin ^2 \theta \sin 2 \phi \frac{\partial ^2}{\partial x \partial z} + \sin 2 \theta \sin \phi \frac{\partial ^2}{\partial x \partial y}+ \sin 2 \theta \cos \phi \frac{\partial ^2}{\partial x \partial y}\end{eqnarray}
$$
{H_2} = \frac{\partial ^2}{\partial x ^2}+\frac{\partial ^2}{\partial y ^2}+\frac{\partial ^2}{\partial z ^2} - {H_1}.
$$
where $\phi$ is an azimuth angle and $\theta$ is a tilt angle.
In this case I am solving for the solution of a pressure wavefield.
EDIT
Is there a physical explanation for self-adjoint operators?
The paper I am referring to can be found here were equation 14 and 15 resemble my postet equations.
 A: I do not know it this is an answer, since I am not sure to have understood your question. The structure of the equation is formally hyperbolic:
$$\frac{\partial^2 \psi}{\partial t^2} - A\psi = S\quad (1)$$
where $\psi =(p,q)^t$. 
If $A$ were self-adjoint and non-negative (or non positive, changing a sign and inserting a further $i$ in front of $\sqrt{-A}$ as I say below), one would construct another self-adjoint operator $\sqrt{A}$ using the spectral theory, and (1) would be re-written as:
$$\left(\frac{\partial}{\partial t} - \sqrt{A} \right)\left(\frac{\partial}{\partial t} + \sqrt{A} \right)\psi(t) = S(t)\:.\quad (2)$$
This equation can be solved interpreting the derivative as a derivative in the strong operator topology in the Hilbert space of the theory. The solution $\psi=\psi(t)$ is a map valued in the said Hilbert space. So for every fixed $t$, in your case, $\psi(t)= \psi(t|\vec{x})$ is an element $L^2(\mathbb R^3)\oplus L^2(\mathbb R^3)$ or some associated Sobolev space. Then one should prove that these solutions are also solutions in proper sense.
Equation (2) has a canonical solution, in the said sense, obtained iterating the solution of 
$$\left(\frac{\partial}{\partial t} \pm  \sqrt{A} \right) \Phi(t) = Z(t)$$
which is:
$$\Phi(t) = \Phi(0) + e^{\mp t\sqrt{A}}\int_0^t e^{\pm \tau\sqrt{A}} Z(\tau) d\tau\:.$$ 
Iteration introduces the first derivative $\partial_t \psi(0)$ as second initial datum, together with $\psi(0)$. The solution of (2) eventually depends on those initial data.
There are however problems with the domains of the involved operators, in general, especially because $e^{t\sqrt{A}}$ is not bounded for $t>0$.
If $A$ is self-adjoint and non-positive, $-A$ is non-negative and thus (2) can be re-written as:
$$\left(\frac{\partial}{\partial t} - i\sqrt{-A} \right)\left(\frac{\partial}{\partial t} + i\sqrt{-A} \right)\psi(t) = S(t)\:,\quad (3)$$
and the solution can be obtained referring to 
$$\left(\frac{\partial}{\partial t} \pm  i\sqrt{-A} \right) \Phi(t) = Z(t)$$
which solves as:
$$\Phi(t) = \Phi(0) + e^{\mp it\sqrt{-A}}\int_0^t e^{\pm i\tau\sqrt{-A}} Z(\tau) d\tau\:.$$ 
The situation with domains here improves because $e^{it\sqrt{-A}}$ is bounded (is unitary) for every $t \in \mathbb R$ and thus its domain is the whole Hilbert space.
If $A$ is not self-adjoint what I wrote above does not apply form scratch since $\sqrt{\pm A}$ is not well-defined (it could still be if $A$ were normal).
