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.

-
You should be much more precise both mathematically and physically. It is impossible to understand (and help you) otherwise. –  V. Moretti Feb 15 at 21:57
Cross-posted from math.stackexchange.com/q/677772/11127 and to mathoverflow.net/q/157685/13917 . In general, it is frown upon to cross-post simultaneously, because it may waste potential answerer's time. As a minimum OP should mention the cross-post (on both sites!). The preferred procedure is to not cross-post, and if the post hasn't received an acceptable answer after, say, a couple of days, then OP could flag for migration. –  Qmechanic Feb 15 at 22:00
yes it is a cross post since they mentioned, that this question would get answered more likely in this forum. I am going to delete the question which does not get the answer at the end –  MichaelScott Feb 15 at 22:06
Could you write what coupled PDE you're solving (you can use MathJax here)? –  Kyle Kanos Feb 15 at 22:08
Barring the mathematical meaning, self-adjointness could mean several things in physics. It depends on the context. I really do not understand your question. What are you asking? –  V. Moretti Feb 15 at 22:09

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).

-
Thanks for a detailed answer. I tried to undesrstand what you did: a) you constructed a self-adjoint operator, b) you gave the solution for $\psi$ which is not bounded for t > 0 if the operators are not bounded? unfortunately I did not get the last part with $e^{it\sqrt{-A}}$. were do you get this from or is it a typo? –  MichaelScott Feb 15 at 23:01
The solution, constructed as I showed could not exists if $S$ does not belong to the right domain. I am completing my answer including the case $i\sqrt{-A}$ –  V. Moretti Feb 15 at 23:04
I stay in Europe, here is getting quite late, I should sleep! Tomorrow I will be back if necessary. Bye –  V. Moretti Feb 15 at 23:08
Thanks a lot! I appreciate your help! unfortunately I am no mathematician and it is quite difficult for me transferring your answer to my problem. So the PDEs form a wave equation of acoustic waves in a medium. S is a source function that produces the initial pressure anomaly which travels through the medium. Now, with the operators given, the amplitudes grow exponentially over time. so I assumed, that the problem was inherent to the diff. operators. as people state, that this was because they were not self-adjoined, I came up with this question. is your answer transferable to my problem? –  MichaelScott Feb 15 at 23:14
Can you find normal modes of your operators $H$? Is there another way to write down $H_1$? To understand its physical origin. Why polar coordinates? –  V. Moretti Feb 16 at 7:48