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

| cite | improve this question | | | | |
  • $\begingroup$ 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. $\endgroup$ – 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.

| cite | improve this answer | | | | |
  • $\begingroup$ "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... $\endgroup$ – J. M. isn't a mathematician 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.

| cite | improve this answer | | | | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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