Von Newman stability analysis for 2D acoustic wave equation explicit

Question

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

\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 \right) + 2 U_{jk}^{n} - U_{jk}^{n-1} \nonumber \\ U_{jk}^{n+1} = \left( \frac{\Delta t V_{jk}}{\Delta s} \right) ^2 \sum_{a=-N}^N w_a \left( U_{j+a k}^n + U_{j k+a}^n \right) + 2 U_{jk}^{n} - U_{jk}^{n-1} \label{eq:w1} \end{eqnarray}(1)

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

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

\begin{equation} U_{jk}^{n+1} = \left( \frac{\Delta t V_{jk}}{\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} - U_{jk}^{n-1} \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{eq:w2} \end{eqnarray}(2)

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

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

$$ 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{eq:w3} \end{equation}(3)

we get:

\begin{eqnarray} 1 = \left( \Omega + 2 \right) \epsilon^{-1} -\epsilon^{-2} \nonumber \\ \quad \text{making} \ \ \epsilon^{-1} = \mu \nonumber \\ \mu^2 - \left( \Omega + 2 \right) \mu + 1 = 0 \nonumber \\ \mu = \frac{(\Omega+2) \pm \sqrt{\Omega^2 + 4\Omega}}{2} \label{eq:w4} \end{eqnarray}(4)

back to expand $ \Omega $ defined in (3):

\begin{eqnarray} \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 \end{eqnarray}

\begin{eqnarray} &=& 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 \end{eqnarray}

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

\begin{eqnarray} &=& 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{eqnarray}

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

\begin{eqnarray} \Omega &=& r^2 \left( 2\cos(\Delta s p) + 2\cos(\Delta s q) - 4\right) \nonumber \\ &=& -4r^2 \left( \sin^2(\frac{\Delta s p}{2}) + \sin^2(\frac{\Delta s q}{2}) \right) \label{eq:w5} \end{eqnarray}(5)

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

We can also write (5) using $ \beta = \left( \sin^2(\frac{\Delta s p}{2}) + \sin^2(\frac{\Delta s q}{2}) \right) $ as :

$$ \Omega = -4r^2\beta $$

Replacing back to (4) :

\begin{eqnarray} \mu &=& \frac{(\Omega+2) \pm \sqrt{\Omega^2 + 4\Omega}}{2} \nonumber \\ \mu &=& -2r^2\beta+1 \pm 2\sqrt{r^2\beta(r^2\beta-1)} \nonumber \end{eqnarray}

I am a little lost how to find if $ | \mu | >= 1 $ or what limitations I have in $r$ for this requirement, that is the same as needing $ | \epsilon | <= 1 $.

Is there any easier alternative to Von Newman that also could be applied to the general explicit form in (1) ?

Answer 1

After more than 2 months and no answer. I post my own answer this is as far as I could get (not final answer though).

I found the general formula for stability criteria (in a paper[1]). That is given by:

$$ r \leq \frac{2}{\sqrt{\sum_{a=-N}^{N} (|w_a^1| + |w_a^2|)}} $$

With $ r = \frac{V \Delta t}{\Delta s} $ and $w_a$ is the centered finite differences weights and the indexes 1 e 2 refer to the x and y dimensions.

But I couldn't get to this general formula, I could just get to the criteria to the 2nd order that was the post N=1.

Not certain if this a proof by contradiction. (Also forgive my bad math I am really eager to learn)

Suppose $ \Delta > 0 $ condition holds for $ | \epsilon | \leq 1 $ that using $ \epsilon^{-1} = \mu $ means $ |\mu| \geq 1 $. Thus this requires for $ r^2 \beta - 1 > 0 $ to be $ r > \frac{1}{ \sqrt{ \beta} } $ that can be satisfied by using $ r = \frac{1}{ \sqrt{ \beta} } + \psi $ with $ \psi > 0 $ positive, real.

Going back to for the first root $ \mu{'} $, we have:

\begin{eqnarray} \mu^{'} &=& -2r^2\beta+1 + 2\sqrt{r^2\beta(r^2\beta-1)} \nonumber \\ &=& -2 \left(1+\frac{2\psi\beta}{\sqrt{\beta}} + \psi^2\beta\right)+1+ 2\sqrt{\left(1+\frac{2\psi\beta}{\sqrt{\beta}} + \psi^2\beta\right)\left[\left(1+\frac{2\psi\beta}{\sqrt{\beta}} + \psi^2\beta\right)-1 \right]} \nonumber \\ &=& -2 \left(1+2\psi\sqrt{\beta}+\psi^2\beta\right) +1 + 2 \sqrt{\left(1+2\psi\sqrt{\beta}+\psi^2\beta\right)\left[\left(1+2\psi\sqrt{\beta}+\psi^2\beta\right)-1\right]} \nonumber \\ &=& -2 A + 1 + 2\sqrt{A^2-A} \nonumber \end{eqnarray}

With:

$$ A = \left(1+2\psi\sqrt{\beta}+\psi^2\beta\right) $$

Note that since $ \psi > 1 $ then $ A > 1 $ always. Using the requirement for stability: \begin{eqnarray} | \mu^{'} | &\geq& 1 \nonumber \\ \left|-2 A + 1 + 2\sqrt{A^2-A}\right| &\geq& 1 \nonumber \\ \end{eqnarray} To satisfy the inequality, two possibilities \begin{eqnarray} -2 A + 1 + 2\sqrt{A^2-A} &\leq& -1 \nonumber \\ -2 A + 1 + 2\sqrt{A^2-A} &\geq& 1 \nonumber \\ -2 A + 2\sqrt{A^2-A} &\leq& -2 \text{ (2)} \nonumber \\ -2 A + 2\sqrt{A^2-A} &\geq& 0 \text{ (3)} \end{eqnarray}

At (2) for $ A > 1 $ left hand side cannot hold, always $ > -2 $. At (3) for $ A > 1 $ also cannot hold, $ -1 > left hand side > -2 $ We don't even need to look at the second root.

This implies that $ \Delta > 0 $ doesn't satisfy the stability criteria.

Now suppose $ \Delta = 0 $ condition holds} for $ | \epsilon | \leq 1 $ that using $ \epsilon^{-1} = \mu $ means $ |\mu| \geq 1 $. Thus this requires that $ r = \frac{1}{ \sqrt{ \beta} } $ \

For booth roots $ \mu $, we have:

\begin{eqnarray} \mu &=& -2r^2\beta+1 + 2\sqrt{r^2\beta(r^2\beta-1)} \nonumber \\ &=& -2+1 \\ &=& -1 \nonumber \end{eqnarray}

That clearly holds.

Finally suppose $ \Delta < 0 $ condition holds for $ | \epsilon | \leq 1 $ Thus this requires $ r^2 \beta - 1 < 0 $ that can be satisfied by using $ r = \frac{1}{ \sqrt{ \beta} } - \psi $ with $ \psi > 0 $ positive, real.

Again going back for the first root $ \mu{'} $, we have:

\begin{eqnarray} \mu^{'} &=& -2r^2\beta+1 + 2\sqrt{r^2\beta(r^2\beta-1)} \nonumber \\ &=& -2 \left(1-2\psi\sqrt{\beta}+\psi^2\beta\right) +1 + 2 \sqrt{\left(1-2\psi\sqrt{\beta}+\psi^2\beta\right)\left[\left(1-2\psi\sqrt{\beta}+\psi^2\beta\right)-1\right]} \nonumber \\ \end{eqnarray} Rearranging due the imaginary part \begin{eqnarray} &=& -2 \left(1-2\psi\sqrt{\beta}+\psi^2\beta\right) +1 + 2 i\sqrt{\left[1-\left(1-2\psi\sqrt{\beta}+\psi^2\beta\right)\right]\left(1-2\psi\sqrt{\beta}+\psi^2\beta\right)} \nonumber \\ &=& -2 A + 1 + 2i\sqrt{A^2-A} \end{eqnarray}

With $i=\sqrt{-1}$ imaginary unit and:

$$ A = \left(1-2\psi\sqrt{\beta}+\psi^2\beta\right) $$

Note that since $ \psi > 1 $ then $ A < 1 $ always. Then using the requirement for stability and complex number modulus (the other root is just conjugate of this so same modulus):

\begin{eqnarray} | \mu | &\geq& 1 \nonumber \\ \left|-2 A + 1 + 2i\sqrt{A^2-A}\right| &\geq& 1 \nonumber \\ \sqrt{(-2 A + 1)^2 + 4(A^2-A) } &\geq& 1 \nonumber \\ \sqrt{ 1 } &\geq& 1 \nonumber \end{eqnarray}

So this condition also holds.

Thus the solution is $ r \leq \frac{1}{ \sqrt{ \beta} } $ that is maximum given $\beta = 2 $ and then $ r \leq \frac{1}{ \sqrt{ 2 } } $ That agrees with general formula presented first for N=1.

[1] A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation - pg T38 - Geophysics Vol. 76 No. 2 2011