# Von Newman stability analysis for 2D acoustic wave equation explicit

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

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

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

$$\Omega = r^2 \sum_{a=-N}^N w_a \left( \phi_{j+a k} + \phi_{j k+a} \right) \label{eq:w3}$$(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) ?

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

• Not an explicit answer per se, but if you're unaware you'll find it extremely easy to simply compute a contour plot of your growth factor $\epsilon$ for various values of $r$ to determine where in the complex plane a given method is stable. As you obviously know, finding a closed-form relation for the stability limits for anything but the most trivial methods quickly becomes an algebraic nightmare, and I don't believe it's a very informative exercise in the end anyway. – Aurelius Aug 25 '13 at 19:47
• Thanks @aurelius that's indeed a good idea I wasn't aware. That's indeed a algebraic nightmare, in the near future for more complex stuff maybe I will try the contour plots. thanks – eusoubrasileiro Aug 27 '13 at 0:28
• No problem, if you have Matlab handy here's a sample for a couple simple time integration methods: spitfire.princeton.edu/stability.m – Aurelius Aug 28 '13 at 14:13