Stability analysis of Heun's method

Question

I am using Heun's method with a third order upwind spatial scheme, which is suggested by Shao (2008) to be used for solving the horizontal advection part of the advection-diffusion equation.

This is what I got:

\begin{align} &C^{\ast} = C^{n} - A^{n} \delta t\\ &C^{n+1} = C^{n} - \dfrac{1}{2}\left( A^{n} + A^{\ast} \right) \delta t \end{align}

Assuming $u=cte>0$, we have,

\begin{equation} A = \dfrac{u}{\delta x} ( \dfrac{1}{6} C_{i-2} - C_{i-1} + \dfrac{1}{2} C_{i} + \dfrac{1}{3} C_{i+1} ) \end{equation}

\begin{equation} \hat{C}_{\ast} = \hat{C}_{n} - \dfrac{u \delta t}{\delta x} ( \dfrac{1}{6} e^{-2imh} - e^{-imh} + \dfrac{1}{2} + \dfrac{1}{3} e^{imh} ) \hat{C}_{n} \end{equation}

\begin{equation} \hat{C}_{\ast} = \left[ 1 - r F(h) \right] \hat{C}_{n} \end{equation}

In which,

\begin{equation} r = \dfrac{u \delta t}{\delta x} \end{equation} \begin{equation} F(h) = \dfrac{1}{2}+\dfrac{1}{6} e^{-2imh} - e^{-imh} + \dfrac{1}{3} e^{imh} \end{equation}

\begin{equation} \hat{C}_{n+1} = \hat{C}_{n} - \dfrac{r}{2} F(h) \hat{C}_{n} - \dfrac{r}{2} F(h) \hat{C}_{\ast} \end{equation}

\begin{equation} g = \dfrac{\hat{C}_{n+1}}{\hat{C}_{n}} = 1 - r F(h) - \dfrac{r^2}{2} F^2(h) \end{equation}

Enforcing \begin{equation} |g| \leq 1 \end{equation}

should give the stability conditions. I have drawn $g$'s graph

and it seems the scheme is always unstable. But I know it is not. What am I doing wrong? Do you have any suggestions?

---

Shao, Yaping. Physics and modelling of wind erosion. Springer, 2008.

Answer 1

Notice that $\hat C_*=(1-r F(h))\hat C_n$, but the sign in front of $r$ is lost when you use $\hat C_*$ inside $\hat C_{n+1}$.

I'm looking at p.149 of Numerical Solution of Time-Dependent Advection-Diffusion-Reaction equations by Hundsdorfer and Verwer on Google books. (I'm going to use their signs.)

The stability region of Heun's method is $$|g(z)| = \big|1+z+\frac12z^2\big|< 1, $$ and if I substitute $$ z = \nu h(\omega) = \nu\big(-\tfrac16 e^{-2i\omega} + e^{-i\omega}-\tfrac12-\tfrac13e^{i\omega} \big), $$ and ask Mathematica to plot the values of $g(\nu h(\omega))$ for $|\omega|<\pi$, $0<\nu<1$, it gives me

which I think is what you should get instead of your plot.

Also, to find the maximal value of $\nu$, I solved $\partial_\omega^4\big |_{\omega=0} |g(\nu h(\omega))| = 0$ (which isn't quite rigorous as this uses a Taylor series approximation), and got $$ 0.87358046473629887, $$ matching Hundsdorfer and Verwer's number of 0.87 in Table 1.2.

Since flipping the signs in your expressions appears to solve the issue and matches Hundsdorfer & Verwer, I think that's the error here.

[Edit to explain 0.87] For any fixed $\nu$, it you plot the function $\omega\mapsto |g(\nu h(\omega))|$, it will have its maximum either at the origin (value 1), or it will be convex at the origin and the maximum will be slightly to the side (plot it and see). The first three derivatives of $|g(\nu h(\omega))|$ w.r.t. $\omega$ at $\omega=0$ vanish identically for all $\nu$, so the function can be approximated as $1+\frac1{24}\omega^4 \partial_\omega^4$. The condition that its maximum value be $1$ can then be approximated (although I think without proof that it's also exactly correct) by the condition that the first non-constant non-zero Taylor term $\partial_\omega^4\big|_{\omega=0}$ has a maximum at $\omega=0$, so to find the maximal value of $\nu$ it is enough to solve $0 = \partial_\omega^4\big|_{\omega=0}$ because it's negative for smaller values of $\nu$.

To solve it:

sage: from mpmath import *
sage: h = lambda w: -(mpf("1/6")*expj(-2*w)-expj(-w)+mpf("1/2")+expj(w)/3)
sage: g = lambda z: 1+z+z**2/2
sage: findroot(lambda t: diff(lambda w: abs(g(t * h(w))), 0, 4), 0.87)
mpf('0.87358046473629887')