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:

\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}

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:

\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}

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}$.

\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}

with

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}

and with the propagation velocity

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

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

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:

\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}

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

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

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)}$$

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

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:

=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

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.

Monday, April 29, 2013

1D Wave Equation Implicit and Von Neumann Stability Analysis

The wave equation 1D acoustic and variable velocity $V_j$ is defined as:

\begin{eqnarray}
\frac{\partial^2 U}{\partial t ^2} = V^2 \frac{\partial^2 U}{\partial x ^2}
\label{imp1}
\end{eqnarray}

And the implicit difference schema can be written as bellow for the unknown time $n+1$ (using the times $n$ and $n-1$ known):

\begin{eqnarray}
U_j^{n+1}  =  \left(\frac{\Delta t  V_j^n}{\Delta s}\right)^2 \left(  U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1}  \right) + 2 U_j^n - U_j^{n-1}
\label{imp2}
\end{eqnarray}
where $\Delta x = \Delta s$. \\
or as:

\left(U_j^{n-1} - 2 U_j^n \right) r_j^{2} =   U_{j+1}^{n+1} - (2 + r_j^2) U_j^{n+1} + U_{j-1}^{n+1}

With $r_j = \frac{\Delta s}{\Delta t V_j^n}$ \\

Applying the boundary condition $U_0^t=U_{N-1}^t=0$, the linear system of equations, in matrix form, for every unknown $U_j^{n+1}$ for $j \in [0, N]$ being $NDs$ the 1D domain for the solution in $x$ is given by:

\begin{eqnarray}
\begin{pmatrix}
-(2+r_1^2)   &      1      &             &            &            \\
1      & -(2+r_2^2)  &      1      &            &            \\
&      1      & -(2+r_3^2)  &     1      &            \\
&             &      1      & -(2+r_4^2) &    1       \\
&             &             &     1      & -(2+r_5^2) \\
\end{pmatrix}
\begin{pmatrix}
\ U_{1}^{n+1} \ \\
\ U_{2}^{n+1} \ \\
\ U_{3}^{n+1} \ \\
\ U_{4}^{n+1} \ \\
\ U_{5}^{n+1} \ \\
\end{pmatrix}
=
\begin{pmatrix}
\left(U_j^{n-1} - 2 U_j^n \right) r_1^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_2^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_3^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_4^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_5^{2} \\
\end{pmatrix}
\end{eqnarray}

Simplified to the case where $N=7$ and $U_{6} = U_{0} = 0$

Von Neumann Stability Analysis

Since $u(x,t) = e^{i(wt+kx)}$ is a solution for (\ref{imp1}) the discrete schema can be written as:

\begin{eqnarray}
u(x_j,t_n) &=& e^{i(wt_n+kx_j)} \nonumber \\
&=& e^{iwn\Delta t} e^{ikj\Delta s} \nonumber \\
&=& \epsilon^n e^{ikj\Delta s}
\label{imp3}
\end{eqnarray}

Where $\epsilon = e^{iw \Delta t}$ and $i = \sqrt{-1}$ .
To maintain stability we should make sure that $\epsilon \leq 1$ not growing exponentially with increasing time steps. So applying (\ref{imp3}) in (\ref{imp2}) we can analyse the growth factor $\epsilon$ .

U_j^{n+1}  =  \left( \frac{\Delta t  V_j^n}{\Delta s} \right) ^2 \left(  U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1}  \right) + 2 U_j^n - U_j^{n-1} \nonumber

\begin{split}
\epsilon^{n+1} e^{ikj\Delta s}  &=  r^2 \left(  \epsilon^{n+1} e^{ik(j+1)\Delta s}  - 2 \epsilon^{n+1} e^{ikj\Delta s} + \epsilon^{n+1} e^{ik(j-1)\Delta s}  \right) \\
&\quad + 2 \epsilon^n e^{ikj\Delta s}  - \epsilon^{n-1} e^{ikj\Delta s} \nonumber
\end{split}

\begin{eqnarray}
1 = r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 2 \epsilon^{-1} - \epsilon^{-2} \nonumber \\
\epsilon^{-2} - 2 \epsilon^{-1} -  r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 1 = 0  \nonumber \\
\epsilon^{-2} - 2 \epsilon^{-1} + 4 r^2 \sin^2 \left( k \Delta s / 2 \right) + 1 = 0
\label{imp4}
\end{eqnarray}

Where $r = \left( \frac{\Delta t V_j^n}{\Delta s} \right)$

At equation (\ref{imp4}) we can substitute $\phi = \epsilon^{-1}$ turning it to a second degree bellow:

\begin{eqnarray}
\phi^2 - 2 \phi +4 r^2 \sin^2 \left(k\Delta s/2 \right) +1 = 0  \nonumber \\
\phi^2 - 2 \phi + c = 0  \nonumber
\end{eqnarray}

With:
$$c = 1 +4 r^2 \sin^2 \left(k\Delta s/2 \right)$$

Has roots $\phi^{'}$ and $\phi^{''}$ as bellow. Replacing $c$ also.

\begin{eqnarray}
&=& \frac{2 \pm \sqrt{4 - 4c}}{2} = 1 \pm \sqrt{1-c} \nonumber \\
&=& 1 \pm \sqrt{-4 r^2 \sin^2 \left(k\Delta s/2 \right) } \nonumber \\
&=& 1 \pm i \ 2 r \sin \left(k\Delta s/2 \right)  \nonumber \\
\phi^{'} &=& 1 + i \ 2 r \sin \left(k\Delta s/2 \right)   \nonumber \\
\phi^{''} &=& 1 - i \ 2 r \sin \left(k\Delta s/2 \right)   \nonumber
\end{eqnarray}

Analyzing the modulus, since $\phi^{'}$ and $\phi^{''}$ are conjugate pairs of the same complex number the modulus is the same for both.
We get:

$$\| \phi \| = \| \phi^{'} \| = \| \phi^{''} \| = \sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}$$

Going back to $\epsilon$ and getting its modulus we get :

\begin{eqnarray}
\| \epsilon^{'} \| &=& \| \epsilon^{''} \| = \frac{1}{\| \phi  \|} \nonumber \\
&=& \frac{1}{\sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}}
\label{imp5}
\end{eqnarray}

Looking at (\ref{imp5}) we can see that always $r^2 > 0$ and $sin^2(x) \in \ [0,1]$ so:

$$\| \epsilon \| \leq 1 \ \ \forall \ r, \ k,\ \Delta s$$

It shows that FTCS in (\ref{imp2}) is unconditionally stable