# Modified Equation and Stability for Centred Finite Differences for Wave Equation

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

• 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. – cfdlab Dec 5 '19 at 13:47
• Interesting! Just reading through your references and the first on explicitly says that there are issues with multi-step methods like Leapfrog. – 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}$$.

• 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? – Daniel Dec 4 '19 at 6:04
• 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? – Daniel Dec 4 '19 at 6:07
• @Daniel Sorry, answer is updated. – Alone Programmer Dec 4 '19 at 20:52
• Thanks, Alone Programmer! I appreciate the answer. But I still can't see where my original argument goes wrong. – Daniel Dec 5 '19 at 7:06
• 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. – cfdlab Dec 5 '19 at 13:54