Monday, November 7, 2016

Regex for files patterns

The `rename` (Perl utility not in all Linux distros), `find` and `grep` are the tools for all renaming of files

lists all files (f) in the folder "folder" and its subfolders
find folder/ -type f

shows the matches of file names by the grep regex Perl expression
find folder/ -type f | grep -P '(\d-\w+-(?:\d+|\d+\w)-ES)' --color

remove duplicate well names in file names appended with '-' (like 1-RSM-1A-ES-1-RSM-1A-ES.lis)
find arranging/ -type f -exec rename -n 's/[-](\d-\w+-(?:\d+|\d+\w)-ES)//g' {} ";"


REGEX summary for BDEP well files ANP

this is the more generic regex pattern for ANP well names
including deviation, repeated etc. 1A, 2D etc. separator is '-'

First a number from 1 to 9. Second a group of characters minimum 2. Third the number of the well. Forth the state where the well was drilled, two characters. A separator between these fields may exist many times or may not even exist.

([1-9])(?:-|_)*([A-Z]+)(?:-|_)*(\d+[A-Z]{0,2})(?:-|_)*([A-Z]{2})

?: turns it in a passive group

http://www.petrobras.com.br/fatos-e-dados/entenda-as-siglas-que-dao-nome-a-um-poco-exploratorio-pela-nomenclatura-da-anp.htm

Tuesday, October 22, 2013

Equivalent Staggered Grid Acoustic Wave Equation

Equivalent Staggered grid  (Di Bartolo et al., 2012) comes from the first order system of equations that come from the continuity equation and Newton's Second Law. First the vectorial linear equation of the motion for the particle velocity field $\mathbf{v}$ is given by:

\begin{equation}
\rho( \mathbf{r}) \frac{\partial \mathbf{v}(\mathbf{r},t)}{\partial t} + \mathbf{\nabla} p(\mathbf{r},t) = \mathbf{f}(\mathbf{r},t)
\label{s1}
\end{equation}

where $\mathbf{r}$ is the position vector, $t$ is the time, $\rho$ os the density of the medium, $\mathbf{f}$ is the density of external body force applied and $\mathbf{\nabla} $ is the gradient operator.
The equivalent scalar equation for pressure $p$ is given by:

\begin{equation}
\frac{1}{\kappa( \mathbf{r})} \frac{\partial p(\mathbf{r},t)}{\partial t} + \mathbf{\nabla} \cdot \mathbf{v}(\mathbf{r},t) = \frac{\partial i_v(\mathbf{r},t)}{\partial t}
\label{s2}
\end{equation}

where $\kappa$ is the adiabatic compression modulus of the medium and $i_v$ , which represents a source, is a volume density injection (an air gun, for example). The $\mathbf{\nabla}$ applies the divergence to the vectorial velocity field $\mathbf{v} $

The Staggered Grid schemes have the great advantage of being able to model regions of water-rock interface, allowing application in the presence of materials with any Poisson's ration (including water), even when density contrasts are considered. In summary are able to model coupling between elastic and acoustic media.

Di Bartolo et. al. 2011 gives an overview of the Staggered Grid methods and its advantages against non-Staggered Methods. Finally its own formulation defined as Equivalent Staggered Grid is based on the following versions of equations (1) and (2), manipulated to remove the vectorial field $\mathbf{v} $.

\begin{equation}
\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) - \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} = -s(\mathbf{r},t)
\label{s3}
\end{equation}

with

\begin{equation}
s(\mathbf{r},t) = \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t}-\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left(\frac{1}{\rho(\mathbf{r})} \mathbf{f}(\mathbf{r},t)\right)
\label{s4}
\end{equation}

and with the propagation velocity

\begin{equation}
c(\mathbf{r}) = \sqrt{\frac{\kappa(\mathbf{r})}{\rho(\mathbf{r})}}
\label{s5}
\end{equation}

For practical implementation the source term $s(\mathbf{r},t)$ is simplified to $ \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} $ using a null density force $ \mathbf{f}(\mathbf{r},t) $

So the equation (3) can be written as :

\begin{eqnarray}
\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + s(\mathbf{r},t) &=& \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} \nonumber \\
c^2(\mathbf{r}) \rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + c^2(\mathbf{r}) \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} &=& \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} \\
\kappa(\mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + \kappa(\mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} &=& \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t}
\end{eqnarray}

From (6) the finite differences marching scheme is defined using the grid properties are $ c(\mathbf{r}) $ and $ \rho( \mathbf{r}) $ and $i_v(\mathbf{r},t)$ can be any combination of punctual sources. The schema is of 2nd order in time using central differences for the second derivative of pressure related to time

For each dimension the divergence term in (6) is written as:

$$ \mathbf{\nabla} \cdot \left( b (\mathbf{r}) \mathbf{\nabla} p(\mathbf{r},t) \right) = \frac{\partial}{\partial x}\left( b (\mathbf{r}) \frac{\partial p}{\partial x}\right)$$

Translated to finite differences and equivalent staggered grid takes the form:

$$ \frac{\partial}{\partial x}\left( b (\mathbf{r}) \frac{\partial p}{\partial x}\right)^n_{i} $$
with $b$ called buoyancy the inverse of density.

$$
=\frac{9}{8} \left[ \frac{b_{i+\frac{1}{2}}}{\Delta x} \left( \frac{9}{8} \frac{ p^n_{i+1}-p^n_i }{\Delta x} - \frac{1}{8} \frac{ p^n_{i+2}-p^n_{i-1} }{3\Delta x}\right) -\frac{b_{i-\frac{1}{2}}}{\Delta x} \left( \frac{9}{8} \frac{ p^n_{i}-p^n_{i-1}}{\Delta x} - \frac{1}{8} \frac{ p^n_{i+1}-p^n_{i-2} }{3\Delta x}\right) \right]
$$
$$
-\frac{1}{8} \left[ \frac{b_{i+\frac{3}{2}}}{3\Delta x} \left( \frac{9}{8} \frac{ p^n_{i+2}-p^n_{i+1} }{\Delta x} - \frac{1}{8} \frac{ p^n_{i+3}-p^n_{i} }{3\Delta x}\right) -\frac{b_{i-\frac{3}{2}}}{3\Delta x} \left( \frac{9}{8} \frac{ p^n_{i-1}-p^n_{i-2}}{\Delta x} - \frac{1}{8} \frac{ p^n_{i}-p^n_{i-3} }{3\Delta x}\right) \right]
$$

better form by Andre

$$
=\frac{9}{64 \Delta x ^2} \left[ b_{i+\frac{1}{2}} \left( 27 \frac{ p^n_{i+1}-p^n_i }{3} - \frac{ p^n_{i+2}-p^n_{i-1} }{3}\right) -b_{i-\frac{1}{2}} \left( 27 \frac{ p^n_{i}-p^n_{i-1}}{3} - \frac{ p^n_{i+1}-p^n_{i-2} }{3}\right) \right]
$$
$$
-\frac{1}{64 \Delta x ^2} \left[ \frac{b_{i+\frac{3}{2}}}{3} \left( 27 \frac{ p^n_{i+2}-p^n_{i+1} }{3} - \frac{ p^n_{i+3}-p^n_{i} }{3}\right) -\frac{b_{i-\frac{3}{2}}}{3} \left( 27 \frac{ p^n_{i-1}-p^n_{i-2}}{3} - \frac{ p^n_{i}-p^n_{i-3} }{3}\right) \right]
$$

$$
=\frac{3}{64 \Delta x ^2} \left[ b_{i+\frac{1}{2}} \left( 27 ( p^n_{i+1}-p^n_i ) - p^n_{i+2}+p^n_{i-1} \right) -b_{i-\frac{1}{2}} \left( 27 (p^n_{i}-p^n_{i-1}) - p^n_{i+1}+p^n_{i-2} \right) \right]
$$
$$
-\frac{1}{576\Delta x ^2} \left[ b_{i+\frac{3}{2}} \left( 27 (p^n_{i+2}-p^n_{i+1} ) - p^n_{i+3}+p^n_{i} \right) -b_{i-\frac{3}{2}} \left( 27 (p^n_{i-1}-p^n_{i-2}) -p^n_{i}+p^n_{i-3} \right) \right]
$$

Then opening $b$ using linear interpolation

$$
=\frac{3}{128 \Delta x ^2} \left[ (b_{i+1}+b_i) \left( 27 ( p^n_{i+1}-p^n_i ) - p^n_{i+2}+p^n_{i-1} \right) -(b_{i-1}+b_i) \left( 27 (p^n_{i}-p^n_{i-1}) - p^n_{i+1}+p^n_{i-2} \right) \right]
$$
$$
-\frac{1}{1152\Delta x ^2} \left[ (b_{i+2}+b_{i+1}) \left( 27 (p^n_{i+2}-p^n_{i+1} ) - p^n_{i+3}+p^n_{i} \right) -(b_{i-2}+b_{i-1}) \left( 27 (p^n_{i-1}-p^n_{i-2}) -p^n_{i}+p^n_{i-3} \right) \right]
$$

 Di Bartolo, L. A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation - Geophysics 2012

Friday, August 9, 2013

Git and Git Hub reminds

Some short reminds:

Going back to a previous push :

(local):
git reset --hard ce82c8108042e555b57013 
where ce82c8108042e555b57013 is the SHA key 

(remote):
git push -f origin

Getting all (remote) branches to (local)

git fetch origin

Cloning a (remote) branch to (local)

git checkout -b  ''  origin/''


...

Thursday, August 1, 2013

Explicit 2D acoustic wave equation

Easy as a cake.

Explicit 2D acoustic wave equation

The scalar wave equation, a fluid medium with variable velocity, is defined by (\ref{ExpA}) Where $U$ and $V$ represents the pressure (or scalar displacement) and velocity fields and S(t) the source, and can be simplified in (\ref{ExpB}) for 2D case:

\begin{eqnarray}
\frac{\partial^2 U}{\partial t^2} \frac{1}{V^2} &=& \nabla^2 U + S(t) \label{ExpA} \\
\frac{\partial^2 U}{\partial t^2} &=& V^2 \left( \frac{\partial^2 U}{\partial x^2}+ \frac{\partial^2 U}{\partial z^2} +S(t) \right)
\label{ExpB}
\end{eqnarray}

From equation (\ref{ExpB}) we can approximate derivatives by Taylor series, manipulating we can reach backward differences in time and centered in space N'th order; using $Dx=Dz=Ds$, remember that $j=jDs$, $k=kDs$ and $n=nDt$.

\begin{eqnarray}
U_{jk}^{n+1} &=& \left( \frac{\Delta t V_{jk} }{\Delta s} \right) ^2 \left( \sum_{a=-N}^N w_a U_{j+a k}^n + \sum_{a=-N}^N w_a U_{j k+a}^n +S_{jk}^n {\Delta s}^2 \right) + 2 U_{jk}^{n} - U_{jk}^{n-1} \\
U_{jk}^{n+1} &=& \left( \frac{\Delta t V_{jk}}{\Delta s} \right) ^2 \left[ \sum_{a=-N}^N w_a \left( U_{j+a k}^n + U_{j k+a}^n \right) \right] + S_{jk}^n {\Delta t}^2 V_{jk}^2 + 2 U_{jk}^{n} - U_{jk}^{n-1} \\
U_{jk}^{n+1} &=& \left( \frac{\Delta t V_{jk}}{\Delta s} \right) ^2 \left[ \sum_{a=-N}^N w_a \left( U_{j+a k}^n + U_{j k+a}^n \right) + S_{jk}^n {\Delta s}^2\right] + 2 U_{jk}^{n} - U_{jk}^{n-1} \label{ExpC}
\end{eqnarray}

For forth order space, we have $N=2$ and $w$ is:
$$ w = \frac{1}{12} [-1, 16, -30, 16, -1] $$

Where, $ R $ must always be bounded (the stability criteria):

\[ R = \frac{\Delta t V_{max}}{\Delta s} \]

Where $V_{max}$ is the maximum velocity.
Simplified from Chen [2] the forth order schema requires :

$$ \Delta t \leq \frac{2 \Delta s}{ V_{max} \sqrt{2/12(1+16+30+16+1)}} \implies \Delta t \leq \frac{ \Delta s \sqrt{3}}{ V_{max} \sqrt{8}} = \frac{ \Delta s \sqrt{3}}{ V_{max} \sqrt{8}}$$

Notes regarding implementation:


* Numerical approximations turn the $R = \frac{\Delta t V_{max}}{\Delta s}$ in higher than accepted depending on the way the sums and products are done
* The Laplacian should not be contaminated, should just multiply by R.
* The last two lines in the equation including equation (\ref{ExpC}) are acceptable implementations.


References

[1] Alford R.M., Kelly K.R., Boore D.M. (1974) Accuracy of finite-difference modeling of the acoustic wave equation Geophysics, 39 (6), P. 834-842

[2] Chen, Jing-Bo (2011) A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation Geophysics, v. 76, p. T37-T42

Simple code example from Madagascar Seismic Package here

Wednesday, May 22, 2013

Digging out end up in Javascript

Ok, I accept this has nothing to do with the previous post ... but anyway since I found this code, I did in college years ago about linear momentum, and translated it to JavaScript using Html5 canvas and found it cool, here I post it just for fun...



This text is displayed if your browser does not support HTML5 Canvas.





Tuesday, May 14, 2013

Miss behavior MathJax and Blogger

If you face any problem seeing equations (numbering), try to visualize the post you want individually.
Unfortunately the numbering of equations gets messed up when viewing multiple posts at once.

Sorry :S For now no workaround!

Implicit 2D acoustic wave equation and Von Newman stability analysis


Implicit 2D acoustic wave equation


The acoustic wave equation, for a medium of constant density, is defined by:

\[ \frac{\partial^2 U}{\partial t^2} = V^2 \nabla^2 U \]

Where $U$ and $V$ represents the pressure (displacement) and velocity fields. For 2D case it gets simplified to:

\begin{equation}
\frac{\partial^2 U}{\partial t^2} = V^2 \left( \frac{\partial^2 U}{\partial x^2}+ \frac{\partial^2 U}{\partial z^2} \right)
\label{1}
\end{equation}

From equation (\ref{1}) we can approximate derivatives by Taylor series, manipulating we can reach backward differences in time and centered in space first and second order; using $Dx=Dz=Ds$, remember that $j=jDs$, $k=kDs$ and $n=nDt$.

\begin{eqnarray}
\frac{U_{jk}^n -2 U_{jk}^{n-1} + U_{jk}^{n-2}}{\Delta t^2} =  \left( \frac{V_{jk}^n}{\Delta s} \right) ^2 \left(  U_{j+1k}^n - 4 U_{jk}^n + U_{jk+1}^n + U_{j-1k}^n + U_{jk-1}^n  \right)
\label{2} \\
U_{jk}^n  =  \left( \frac{\Delta t  V_{jk}^n}{\Delta s} \right) ^2 \left(  U_{j+1k}^n - 4 U_{jk}^n + U_{jk+1}^n + U_{j-1k}^n + U_{jk-1}^n  \right) + 2 U_{jk}^{n-1} - U_{jk}^{n-2} \nonumber
\end{eqnarray}

Rearranging to solve in matrix form for the unknowns $ U^n $ (time n) around $U_{jk}$, we get

\begin{eqnarray}
 U_{jk}^{n-2} -2 U_{jk}^{n-1} &=&  \left( \frac{\Delta t  V_{jk}^n}{\Delta s} \right) ^2 \left[  U_{j+1k}^n - 4 U_{jk}^n - \left( \frac{\Delta s}{\Delta t  V_{jk}^n} \right) ^2  U_{jk}^n + U_{jk+1}^n + U_{j-1k}^n + U_{jk-1}^n  \right] \nonumber \\
\left( \frac{\Delta s}{ \Delta t  V_{jk}^n} \right) ^2 \left( U_{jk}^{n-2} -2 U_{jk}^{n-1} \right) &=& U_{j+1k}^n - U_{jk}^n \left[ 4+ \left( \frac{\Delta s}{\Delta t  V_{jk}^n} \right) ^2  \right]  + U_{jk+1}^n + U_{j-1k}^n + U_{jk-1}^n \label{3}
\end{eqnarray}

Note: Centered differences in  (\ref{2}) would reach a explicit method with stability restrictions. So although (\ref{2}) is first order in time the analyses in the end shows it's unconditionally stable.

Rewritting (\ref{3}) to a simplified form, introducing some new variables:

\begin{eqnarray}
\left( \frac{\Delta s}{ \Delta t  V_{jk}^n} \right) ^2 \left( U_{jk}^{n-2} -2 U_{jk}^{n-1} \right) =
  U_{j-1k}^n - \left[ 4+ \left( \frac{\Delta s}{\Delta t  V_{jk}^n} \right) ^2  \right] U_{jk}^n + U_{j+1k}^n + U_{jk+1}^n + U_{jk-1}^n  \nonumber \\
r^{-2} \left( U_{jk}^{n-2} -2 U_{jk}^{n-1} \right) =
  U_{j-1k}^n + \gamma \  U_{jk}^n + U_{j+1k}^n + U_{jk+1}^n + U_{jk-1}^n
\label{4}
\end{eqnarray}

Where:

$$ r = \frac{\Delta t  V_{jk}^n}{ \Delta s} $$
$$ \gamma = -4 -r^{-2}  $$

For a tiny grid 3x3 dimensions we can apply a Dirichlet boundary condition in (\ref{4}). That means $ U_{\Omega} = 0 \ \forall \ n $ and $ U_{jk}^n $ with $ j, \ k \in Z \ \left[ 0, 2  \right] $ , illustrated bellow:

\begin{array}{| c | c | c | c | c |}
\hline  U_{\Omega} & U_{\Omega} & U_{\Omega} & U_{\Omega} & U_{\Omega} \\\hline
U_{\Omega} & U_{00}^n &  U_{10}^n &  U_{20}^n & U_{\Omega} \\\hline
U_{\Omega} & U_{10}^n &  U_{11}^n & U_{21}^n & U_{\Omega}\\\hline
U_{\Omega} & U_{20}^n &  U_{21}^n & U_{22}^n & U_{\Omega} \\\hline
U_{\Omega} & U_{\Omega} & U_{\Omega} & U_{\Omega} & U_{\Omega} \\\hline
\end{array} (4.a)

\begin{array}{| c | c | c | c | c |}
\hline 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\hline
0.0 & P_{0} & P_{1} &  P_{2}  & 0.0 \\\hline
0.0 & P_{3} & P_{4} & P_{5} & 0.0 \\\hline
0.0 & P_{6} & P_{7} & P_{8} & 0.0 \\\hline
0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\hline
\end{array} (4.b)

We can assembly the linear system for this simple grid as a 9x9 matrix, using the equation (\ref{4}) for every unknown $ P_i $ . Not showing the independent term $  r^{-2} \left( U_{jk}^{n-2} -2 U_{jk}^{n-1} \right) $. Calculation of the coefficients can be made simple  using (4.b), empty position means 0 as coeficiente; also represented by $ \varnothing $.

\begin{array}{ l }
P_0 = \varnothing  0  1  \varnothing 3  \\
P_1 = 0 1 2 \varnothing 4 \\
P_2 = 1 2 \varnothing \varnothing 5 \\
P_3 = \varnothing 3 4 0 6 \\
P_4 = 3 4 5 1 7 \\
P_5 = 4 5 \varnothing 2 8 \\
P_6 = \varnothing 6 7 3 \varnothing \\
P_7 = 6 7 8 4 \varnothing \\
P_8 = 7 8 \varnothing 5 \varnothing \\
\end{array}


\begin{array}{| c | c | c | c | c | c | c | c | c | c |}
\hline
         &  P_0  &  P_1  &  P_2  &  P_3  &  P_4  &  P_5  &  P_6  &  P_7  &  P_8  \\ \hline
  P_0  & \gamma &    1    &         &    1    &         &         &         &         &         \\ \hline
  P_1  &    1    & \gamma &     1   &         &     1   &         &         &         &         \\ \hline
  P_2  &         &    1    & \gamma &         &         &    1    &         &         &         \\ \hline
 P_3  &    1    &         &         & \gamma &    1    &         &    1    &         &         \\ \hline
  P_4  &         &    1    &         &    1    & \gamma &    1    &         &    1    &         \\ \hline
  P_5  &         &         &    1    &         &    1    & \gamma &         &         &    1    \\ \hline
  P_6  &         &         &         &    1    &         &         & \gamma &    1    &         \\ \hline
  P_7  &         &         &         &         &    1    &         &    1    & \gamma &    1    \\ \hline
     P_8  &         &         &         &         &         &    1    &         &    1    & \gamma \\ \hline
\end{array}

Or the full human form, matrix notation using the independent term :

\begin{eqnarray}
\begin{pmatrix}
\gamma   &    1    &         &    1    &         &         &         &         &         \\
    1    & \gamma  &     1   &         &     1   &         &         &         &         \\
         &    1    & \gamma  &         &         &    1    &         &         &         \\
    1    &         &         & \gamma  &    1    &         &    1    &         &         \\
         &    1    &         &    1    & \gamma  &    1    &         &    1    &         \\
         &         &    1    &         &    1    & \gamma  &         &         &    1    \\
         &         &         &    1    &         &         & \gamma  &    1    &         \\
         &         &         &         &    1    &         &    1    & \gamma  &    1    \\
         &         &         &         &         &    1    &         &    1    & \gamma  \\
\end{pmatrix}
\begin{pmatrix}
 \ U_{00}^n \ \\
 \ U_{01}^n \ \\
 \ U_{02}^n \ \\
 \ U_{10}^n \ \\
 \ U_{11}^n \ \\
 \ U_{12}^n \ \\
 \ U_{20}^n \ \\
 \ U_{21}^n \ \\
 \ U_{22}^n \ \\
\end{pmatrix}
=
r^{-2}
\begin{pmatrix}
U_{00}^{n-2} \\
U_{01}^{n-2} \\
U_{02}^{n-2} \\
U_{10}^{n-2} \\
U_{11}^{n-2} \\
U_{12}^{n-2} \\
U_{20}^{n-2} \\
U_{21}^{n-2} \\
U_{22}^{n-2} \\
\end{pmatrix}
-2r^{-2}
\begin{pmatrix}
U_{00}^{n-1} \\
U_{01}^{n-1} \\
U_{02}^{n-1} \\
U_{10}^{n-1} \\
U_{11}^{n-1} \\
U_{12}^{n-1} \\
U_{20}^{n-1} \\
U_{21}^{n-1} \\
U_{22}^{n-1} \\
\end{pmatrix}
\end{eqnarray}

Analyzing above a recursive schema can be created. From a grid of dimensions $ N_x \times N_z $ we have a linear system $L$ with dimension $ M \times M $ ( $ M = N_x N_z $) . We have five rules for every line in $ L $, looking at $ L $ as a list of lines, list with dimension $ M $, these rules are based upon the principal diagonal, defined as $L[i][i]$.

Where $ i = \ \in [0, M] $

$$ j = i \% N_x \nonumber \ \mbox{integer divison remainder} $$
$$ k = i / N_x \nonumber \ \mbox{integer division} $$
$$ U_{j-1k}^n = 1 \ \ \forall \ j-1 \geq 0 \implies L[i][i-1] = 1 \nonumber $$
$$ U_{jk}^n =  \gamma \ \ \implies L[i][i] = \gamma \nonumber $$
$$ U_{j+1k}^n = 1 \ \ \forall \ j+1 < N_x \implies L[i][i+1] = 1 \nonumber $$
$$ U_{jk-1}^n = 1 \ \ \forall \ z-1 \geq 0 \implies L[i][i-N_x] = 1 \nonumber $$
$$ U_{jk+1}^n = 1 \ \ \forall \ z+1 < N_z \implies L[i][i+N_x] = 1 \nonumber $$

Can also be seeing, looking each diagonal as a list (reminded that we have 5 defined diagonals):

* Leading diagonal defined by $L[p] = \gamma_{jk}$ with $p \in [0, M] $

* Leading diagonal upper +1 $L[p] = 1$ with $p \in [0, M-1] $
   Leading diagonal lower -1  $L[p] = 1$ with $p \in [0, M-1] $
   Except where $ p \%Nx == 0 $

* Leading diagonal upper in $+Nx$ $L[p] = 1$ with $p \in [0, M-Nx] $
   Leading diagonal lower in $-Nx$  $L[p] = 1$ with $p \in [0, M-Nx] $

This matrix similar to the 2D poisson matrix, with $5N_xN_z-2N_x-2N_z$ non zero elements.

Von Newman stability analysis


Von Newman stability analysis for acoustic wave equation implicit centered differences: 1st order time and space (N 2)'th order:

\begin{eqnarray}
U_{jk}^{n+2}  =  \left( \frac{\Delta t  V_{jk} }{\Delta s} \right) ^2 \left(  \sum_{a=-N}^N w_a U_{j+a k}^{n+2} + \sum_{a=-N}^N w_a U_{j k+a}^{n+2} \right) + 2 U_{jk}^{n+1} - U_{jk}^n  \\
U_{jk}^{n+2}  =  \left( \frac{\Delta t  V_{jk}}{\Delta s} \right) ^2  \sum_{a=-N}^N  w_a \left( U_{j+a k}^{n+2} + U_{j k+a}^{n+2} \right) + 2 U_{jk}^{n+1} - U_{jk}^n \label{a6}
\end{eqnarray}

For forth order space, we have $N=2$ and $w$ is:
$$ w = \frac{1}{12} [-1, 16, -30, 16, -1] $$

Can be simplified to 1st order (N=1):

\begin{equation}
U_{jk}^{n+2}  =  \left( \frac{\Delta t  V_{jk}}{\Delta s} \right) ^2 \left(  U_{j+1k}^{n+2} - 4 U_{jk}^{n+2} + U_{jk+1}^{n+2} + U_{j-1k}^{n+2} + U_{jk-1}^{n+2}  \right) + 2 U_{jk}^{n+1} - U_{jk}^{n} \nonumber
\end{equation}

Using the discrete solution for 2D wave equation, where $ i = \sqrt{-1} $, $ n = n \Delta t $, $ j = j \Delta x $ and $ k = k \Delta z $. Last using $ \Delta x = \Delta z = \Delta s $, follows that the discrete solution can be written as:

\begin{eqnarray}
U_{jk}^n = e^{i \left( \omega t + px + qz \right)} \nonumber \\
U_{jk}^n = \epsilon^n e^{i \left( pj\Delta s + qk\Delta s \right)}  \nonumber \\
U_{jk}^n = \epsilon^n e^{i \Delta s \left( pj + qk \right)}  \label{a7}
\end{eqnarray}

Where $\epsilon $ is the growth factor, and should be $ |\epsilon| \leq 1$ for stability. \\

Replacing (\ref{a7}) in (\ref{a6}), using the identities bellow and simplifying dividing both sides by $ U_{jk}^{n+2} $

$$ r = \frac{\Delta t  V_{jk}}{\Delta s} $$
$$ \phi_{j+l\ k+m} = e^{i \Delta s \left( pl+qm \right)} $$

\begin{equation}
\Omega = r^2 \sum_{a=-N}^N  w_a \left( \phi_{j+a k} + \phi_{j k+a}  \right) \label{a8}
\end{equation}

we get:

\begin{eqnarray}
1  =  \Omega + 2 \epsilon^{-1} -\epsilon^{-2} \nonumber \\
\quad \text{making} \ \ \epsilon^{-1} = \mu \nonumber \\
\mu^2 - 2 \mu + 1 -\Omega = 0 \nonumber \\
\mu = \frac{ 2 \pm \sqrt{ 4 - 4 +4 \Omega}}{2} \nonumber \\
\mu = 1 \pm \sqrt{ \Omega} \label{a9}
\end{eqnarray}

back to expand $ \Omega $ defined in (\ref{a8}):

$$ \Omega = r^2 \sum_{a=-N}^N  w_a \left( \phi_{j+a k} + \phi_{j k+a}  \right) \nonumber  $$
$$ = r^2 \sum_{a=-N}^{N} w_a ( e^{i \Delta s \ p a} + e^{i \Delta s \ q a} ) \nonumber $$

$$ =r^2
\begin{pmatrix}
 \cdots & e^{-i \Delta s 2 p} + e^{-i \Delta s 2 q} & e^{-i \Delta s p} + e^{-i \Delta s q} & e^0+e^0 & e^{i \Delta s p} + e^{i \Delta s q} & e^{i \Delta s 2 p} + e^{i \Delta s 2 q} & \cdots \\
\end{pmatrix}
\begin{pmatrix}
\cdots \\
w_{-2} \\
w_{-1} \\
w_0 \\
w_1 \\
w_2 \\
\cdots
\end{pmatrix}
\nonumber $$

Since $w$ is even $ w_a = w_{-a} $ and $ e^{i\theta} + e^{-i\theta} = 2 \cos{\theta} $ we can rewrite as:

\begin{equation}
=r^2
\begin{pmatrix}
 \cdots & 2\cos( \Delta s 2 p) + 2\cos(\Delta s 2 q) & 2\cos(\Delta s p) + 2\cos(\Delta s q) & 2 \\
\end{pmatrix}
\begin{pmatrix}
\cdots \\
w_{2} \\
w_{1} \\
w_0 \\
\end{pmatrix}
\nonumber
\end{equation}

For the simplest case 2nd order $ N=1 $ we have $ (w_1, w_0) = (1, -2) $

$$ \Omega = r^2 \left( 2\cos(\Delta s p) + 2\cos(\Delta s q) - 4\right) $$
$$ = -4r^2 \left( \sin^2(\frac{\Delta s p}{2}) + \sin^2(\frac{\Delta s q}{2}) \right) \label{a10} $$

Note:  $ 2 \cos(\theta) - 2 = -4 \sin ^2 (\theta) $ .

Replacing back to (\ref{a9}) :

$$ \mu = 1 \pm \sqrt{\Omega} $$
$$ \mu = 1 \pm  i 2 \sqrt{r^2 \left( \sin^2(\frac{\Delta s p}{2}) + \sin^2(\frac{\Delta s q}{2}) \right)} $$

Back to $\epsilon$ :

\[ \frac{1}{a^2+b^2} = | \epsilon^{'} | = | \epsilon^{''} | = \frac{1}{1 + 4r^2 \left( \sin^2(\frac{\Delta s p}{2}) + \sin^2(\frac{\Delta s q}{2}) \right)} \]

That is for $ \forall \ r, \ \Delta \ s, \ p, \ q \ \implies | \epsilon | \leq 1 $

Follows (\ref{a6}) with $N=1$ is unconditionally stable.