# A-stability of Heun method for ODEs

I'm trying to determine the stability region of the Heun method for ODEs by using the equation $y' = ky$, where $k$ is a complex number, based on the method described here.

If the Heun method is:

$$y_{n+1} = y_n + 0.5\cdot h\bigl(f(t_n, y_n) + f(t_{n+1},y_n + 0.5\cdot h\cdot f(t_n, y_n)\bigr)$$

then when I insert $y' = zy$ for $f(t,y)$, my result simplifies to

$$y_{n+1} = (0.25\cdot h^2 \cdot z^2 + hz + 1)y_n$$

to judge from the wiki article, the stability region is then the area described by

$$\\{z \in \mathbb C \mid 0.25h^2z^2 + hz + 1 < 1\\}$$

Am I doing this right? What would such a region look like? Can someone help me get the intuition for this? And then I guess the method is A-stable if that region includes wherever $\Re < 0$?

-
Please, may I ask you?I want to know the stabilities such asbsolue stability, asymptitical stabiliy, A-stability and B - stability for numerical integration method.Thanks for this. –  user51346 Dec 1 '12 at 2:01

For the method that you wrote down, you indeed have $$y_{n+1} = (0.25\cdot h^2 \cdot z^2 + hz + 1)y_n.$$ The expression between parentheses is a function of $w = hz$, and the stability region consists of the numbers $w$ for which the modulus of the expression between parentheses is at most one: $$\text{stability region} = \{ w \in \mathbb{C} : |0.25w^2 + w + 1| < 1 \}.$$ So the stability region does not depend on the step size $h$.

In the particular case you are asking about, you can use that $$0.25w^2 + w + 1 = (0.5w+1)^2.$$ You should be able to find the region from this.

You are right about how to determine whether the method is A-stable. You should find that the method is not A-stable, because explicit methods are never A-stable (see also the Wikipedia page that you link to).

In general, to get a feeling for what the stability region looks like, one may start by restricting to the real axis. If $w$ is real, then $0.25w^2+w+1$ is also real, so the condition $|0.25w^2+w+1| < 1$ simplifies to $-1 < 0.25w^2+w+1 < 1$. But what people often do in practice is to plot the region using the computer.

A final note: Are you sure you copied the method correctly? The method $$y_{n+1} = y_n + 0.5\cdot h\bigl(f(t_n, y_n) + f(t_{n+1},y_n + h\cdot f(t_n, y_n)\bigr)$$ with no factor 0.5 in front of the last $h$, is more popular.

-
"explicit methods are never A-stable" - I was waiting for this to be mentioned. People would have never been interested in implicit methods if explicit methods worked nicely on stiff sets of equations... –  Guess who it is. Apr 27 '12 at 10:55

I think you mean to type(?) $$\{z \in \mathbb C : |0.25h^2z^2 + hz + 1| < 1\}$$ Why to impose this condition? Consider the following at time step $n$: $$y_{n} = \Big(\frac{h^2 z^2}{4} + hz + 1\Big)^n \,y_0$$ therefore having the first stability condition is the sufficient condition for a small disturbance in the initial value wouldn't get magnified marching in the time step: Consider $y_0^{\epsilon} = y_0 + \epsilon$, then $$|y^{\epsilon}_{n} - y_n| = \Big\vert\frac{h^2 z^2}{4} + hz + 1\Big\vert^n \,\epsilon < \epsilon$$

and $\epsilon$ can be happen at any time-step, or even could be the numerical error itself, this is why sometimes for a relatively not-so-small step-size $h$, after a few iterations, the numerical error gets magnified and the numerical solution blows up.

As for question about A-stableness, based on wiki, I am guessing the answer is Yes.

-