Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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$?

share|cite|improve this question
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.

share|cite|improve this answer
"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... – J. M. 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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.