Modified Equation and Stability for Centred Finite Differences for Wave Equation

Question

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.,

Answer 1

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