Heun's method:
$$w_{n+1} = w_n + \frac{h}{4}[f(t_n,w_n) + 3f(t_n+2/3*h, w_n +2/3*hf(t_n,w_n))]$$
Question:
Is Heun's method numerically stable for any choice of step size $h$? Why or why not?
I'm assuming we need to calculate the stability region for which after applying it to the test equation $y' = ky$ I got:
$$y_{n+1} = (1+kh+0.5(kh)^2)y_n$$
and so the stability region is for $z = kh$
$$|(1+z+0.5z^2)| \leq1$$
I'm stuck from here. How can I find the values $h$ for which this inequality is satisfied? I don't think I understand A stability too well. I think we need to take Re($k$) < 0 and solve right? I have that the method is stable for $0<z<2$ but I feel like this is wrong. Can anyone shed some light?
Best Answer
EDIT: We can rewrite $1+z+0.5z^2$ to $0.5(z+1)^2+0.5$ by completing the square. Since the function is always positive we can drop the absolute value, leading $0.5(z+1)^2+0.5\leq 1$. This is equivalent to $\lvert z+1\rvert\leq1$. This gives $z\leq0$ and $z\geq -2$.
Then, the stability interval is $$\mathrm{SI}=[-2,0].$$ The values of $h$ the method is stable for depends on the ODE. You have to calculate the eigenvalues $\lambda_i$ of the Jacobi matrix $f_x(t,x)$ of the right side of the ODE. Then the method is stable if $\lambda_i h\in \mathrm{SR}$ where $\mathrm{SR}$ is the stability region. Note that this is only an approximation. It does not need to be true for every ODE but works most of the time (at least if $f_x(t,x)$ is diagonalisable).
A method is called A-stable iff the complete left half plane $\{z\in\mathbb C\setminus\{0\} : \Re(z)\leq0\}$ is contained in the stability region $\mathrm{SR}$. This makes sure that our numerical method falls off to zero as our solution to the test equation. A-stable methods as the trapezoidal rule are used for stiff problems where other methods would need very small step sizes or explode for $t\to\infty$.