I am trying to use the modified equation to derive the stability condition for the finite difference approximation

$$ \frac{u(x,t+\Delta t) - 2 u(x, t) + u(x, t -\Delta t)}{\Delta t^2} = c^2 \frac{u(x+h,t) - 2 u(x,t) + u(x-h,t)}{h^2} $$

to the wave equation

$$ u_{tt}(x,t) = c^2 u_{xx}(x,t). $$

I know that the resulting stability limit should be

$$ \frac{c^2 \Delta t^2}{h^2} \leq 1 $$

but this is not what I am getting. Please tell me where I go wrong! Here is my line of argument:

Truncation Error. If we insert the continuous solution $u(x,t)$ into the left hand side of the finite difference stencil and use Taylor expansions \begin{align*} u(x, t + \Delta t) &= u(x,t) + \Delta t u_t(x,t) + \frac{1}{2} \Delta t^2 u_{tt}(x,t) + \frac{1}{6} \Delta t^3 u_{ttt} + \frac{1}{24} \Delta t^4 u_{tttt} + \mathcal{O}(\Delta t^5) \\ u(x, t - \Delta t) &= u(x,t) - \Delta t u_t(x,t) + \frac{1}{2} \Delta t^2 u_{tt}(x,t) - \frac{1}{6} \Delta t^3 u_{ttt} + \frac{1}{24} \Delta t^4 u_{tttt} + \mathcal{O}(\Delta t^5) \\ \end{align*} we can show that $$ \frac{u(x,t+\Delta t) - 2 u(x, t) + u(x, t -\Delta t)}{\Delta t^2} = u_{tt}(x,t) + \frac{1}{12} \Delta t^2 u_{tttt} + \mathcal{O}(\Delta t^3) $$

Making an almost identical argument for the right hand side, we find that $$ c^2 \frac{u(x+h,t) - 2 u(x,t) + u(x-h,t)}{h^2} = c^2 u_{xx}(x,t) + \frac{1}{12} c^2 h^2 u_{xxxx}(x,t) + \mathcal{O}(h^3). $$

Taken together, we find that

$$ u_{tt}(x,t) - c^2 u_{xx}(x,t) + \frac{1}{12} \left( \Delta t^2 u_{tttt} - c^2 h^2 u_{xxxx} \right) = T(x,t) + \mathcal{O}(\Delta t^3) + \mathcal{O}(h^3). $$

where $T(x,t)$ is the local truncation error. This is a standard argument to show that the finite difference stencil is second order accurate.

Modified equation. Following ideas introduced by Warming and Hyett, I now derive the modified equation. Since $u(x,t)$ is assumed to be the continuous solution, if it is smooth enough we have $u_{tttt} = c^2 u_{xxtt} = c^2 u_{ttxx} = c^4 u_{xxxx}$. Therefore, the expression above becomes

$$ u_{tt} - c^2 u_{xx} + \frac{c^2}{12} \left( c^2 \Delta t^2 - h^2 \right) u_{xxxx} = T(x,t) + \mathcal{O}(\Delta t^3) + \mathcal{O}(h^3). $$

Thus the finite difference is a second order approximation to the wave equation, but a third order approximation to the modified equation

$$ u_{tt} - c^2 u_{xx} + a u_{xxxx} = 0 $$

with $a = \frac{c^2}{12} \left( c^2 \Delta t^2 - h^2 \right)$ and we can analyse the behaviour of this equation to understand how our finite difference behaves.

Stability. Finally, to assess stability, we insert a plane wave

$$ u(x,t) = e^{i (k x - \omega t)} = e^{i (k x - \mathbf{R}(\omega))} e^{\mathbf{I}(\omega) t} $$

into the modified equation and figure out the dispersion relation. A frequency $\omega$ with a positive imaginary part means a solution that grows exponentially in time, indicating instability. Inserting the plane wave into the modified equation yields the dispersion relation

$$ \omega = \pm \sqrt{ c^2 k^2 + a k^4 } = \pm c k \sqrt{1 + \frac{a}{c^2} k^2}. $$

Now for $a > 0$, the radicand is always positive and the root remains real. Therefore, $\omega$ does not have an imaginary part and the solution remains stable.

But: $a > 0$ corresponds to $c^2 \Delta t^2 - h^2 > 0$ or $\frac{c^2 \Delta t^2}{h^2} > 1$ which is clearly nonsense, given that the actual stability criterion is the other way round.

This looks like there should be a stupid, simple sign error somewhere but I can't seem to find it. Any help is much appreciated.

Just to clarify, my question is where my argument goes wrong. I am aware that there are other ways to derive the stability condition.

Warming, R. F.; Hyett, B. J., The modified equation approach to the stability and accuracy analysis of finite-difference methods, J. Comput. Phys. 14, 159-179 (1974). ZBL0291.65023.,

  • $\begingroup$ What you wrote seems ok, though you should use the scheme to do the elimination of higher derivatives rather than the pde, as pointed out by Warming/Hyett. I have not seen anyone apply this to second order wave equation. There seems to be many issues with this approach, see e.g., doi.org/10.1016/0021-9991(90)90093-G and doi.org/10.1137/0907067, and its probably better not to use this approach to decide on scheme stability. $\endgroup$
    – cfdlab
    Dec 5 '19 at 13:47
  • $\begingroup$ Interesting! Just reading through your references and the first on explicitly says that there are issues with multi-step methods like Leapfrog. $\endgroup$
    – Daniel
    Dec 5 '19 at 13:54

I'm rewriting my answer. In fact, you don't need Taylor expansion to find out why $\frac{c \Delta t}{\Delta x} < 1$. I define second order numerical time and spatial differential operators as respectively:

$$D_{tt} u = \frac{u(x,t+\Delta t) + u(x,t-\Delta t) - 2 u(x,t)}{\Delta t^{2}}$$

$$D_{xx} u = \frac{u(x+\Delta x,t)+u(x-\Delta x,t)-2u(x,t)}{\Delta x^{2}}$$

The discrete form of your original wave equation is:

$$D_{tt} u = c^{2} D_{xx} u$$

And your solution is: $u(x,t) = \exp(ikx-i\omega t)$

So by putting that solution into that discretized wave equation and use discretized differential operators, we have:

$$D_{tt} u = -\frac{4}{\Delta t^{2}} \sin^{2}\Big(\frac{\omega \Delta t}{2}\Big) \exp(ikx-i\omega t)$$

$$D_{xx} u = -\frac{4}{\Delta x^{2}} \sin^{2}\Big(\frac{k\Delta x}{2}\Big) \exp(ikx-i\omega t)$$

So, finally:

$$\Big |\sin\Big(\frac{\omega \Delta t}{2}\Big)\Big| = \frac{c\Delta t}{\Delta x} \Big | \sin\Big(\frac{k\Delta x}{2}\Big) \Big |$$

You know that always: $\Big | \sin \Big ( \frac{\omega \Delta t}{2} \Big ) \Big | < 1$, so if $\frac{c \Delta t}{\Delta x} > 1$, for some values of $k$, which obviously your numerical discretization will be unstable, you would have $\frac{c\Delta t}{\Delta x} \Big | \sin \Big ( \frac{k \Delta x}{2} \Big ) \Big | > 1$ that is obviously wrong. So, you should always have $\frac{c \Delta t}{\Delta x} < 1$ to maintain stability of your discretized scheme.

It's good to see that for small $\Delta t$ and $\Delta x$, you have:

$$\lim_{\Delta t \rightarrow 0} D_{tt} u = \lim_{\Delta t \rightarrow 0} -\frac{4}{\Delta t^{2}} \sin^{2} \Big ( \frac{\omega \Delta t}{2} \Big )\exp(ikx-i\omega t) = -\omega^{2} \exp(ikx-i\omega t) = \frac{\partial^{2} u}{\partial t^{2}}$$

$$\lim_{\Delta x \rightarrow 0} D_{xx} u = \lim_{\Delta x \rightarrow 0} -\frac{4}{\Delta x^{2}} \sin^{2} \Big ( \frac{k \Delta x}{2} \Big )\exp(ikx-i\omega t) = -k^{2} \exp(ikx-i\omega t) = \frac{\partial^{2} u}{\partial x^{2}}$$

Furthermore, you can call $\sqrt{\frac{4}{\Delta t^{2}} \sin^{2} \Big ( \frac{\omega \Delta t}{2} \Big )}$ and $\sqrt{\frac{4}{\Delta x^{2}} \sin^{2} \Big ( \frac{k \Delta x}{2} \Big )}$ numerical frequency ($\tilde{\omega}$) and numerical wave number ($\tilde{k}$) respectively. So, the numerical dispersion relation is:

$$\tilde{\omega}^{2} = c^{2} \tilde{k}^{2}$$

Where you can easily deduce $\frac{c\Delta t}{\Delta x} < 1$, due to the fact that $\tilde{\omega} < \frac{2}{\Delta t}$ or $\tilde{k} < \frac{2}{\Delta x}$.

  • $\begingroup$ Wait, I agree up to u_xx = -k^2! But now the next derivative adds another +ik, hence u_xxx = -i k^3. Then another ik for u_xxxx = -i k^3 * ik = -i * i * k^4 = + k^4. Am I missing something? $\endgroup$
    – Daniel
    Dec 4 '19 at 6:04
  • $\begingroup$ Although this feels embarrassing, I put it into Wolfram Alpha - it confirms that (i*k)^4 = + k^4 ... hence the fourth derivative should come up with a positive sign in the dispersion relation? $\endgroup$
    – Daniel
    Dec 4 '19 at 6:07
  • $\begingroup$ @Daniel Sorry, answer is updated. $\endgroup$ Dec 4 '19 at 20:52
  • $\begingroup$ Thanks, Alone Programmer! I appreciate the answer. But I still can't see where my original argument goes wrong. $\endgroup$
    – Daniel
    Dec 5 '19 at 7:06
  • $\begingroup$ For multi-level methods as found for second order wave equation, the correct procedure to do Fourier stability analysis involves finding the roots of a characteristic equation and checking that all of them are less than unity in magnitude. See e.g., Section (1.9) of Gustaffson, Kreiss, Oliger: Time dependent problems and difference methods. $\endgroup$
    – cfdlab
    Dec 5 '19 at 13:54

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.