4
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Aurelius Aug 25 '13 at 19:47
  • $\begingroup$ 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 $\endgroup$ – eusoubrasileiro Aug 27 '13 at 0:28
  • 1
    $\begingroup$ No problem, if you have Matlab handy here's a sample for a couple simple time integration methods: spitfire.princeton.edu/stability.m $\endgroup$ – Aurelius Aug 28 '13 at 14:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.