Absolute stability of numerical methods for ODEs

numerical methodsordinary differential equationsrunge-kutta-methodsstability-in-odes

I've troubles understanding the meaning of region of absolute stability for numerical methods for ODEs. I know that we can restric the study of stability of a certain method to the case of the test equation:
$$ \begin{cases}
y'(x) &= \lambda y(x) \\
y(0) &= 1 \\
\end{cases} $$
and that for Runge-Kutta methods the approximations $ \{ y_{n} \} $ will tend to $0$ for $n \rightarrow \infty$ if and only if the roots of the characteristic equation (with h standing for the lenght of the step) $$ \sum_{i}^{k}(\alpha_{i} -\beta_{i}h\lambda)t^{k-i} $$ have norm $<1$. The region of stability for the trapezoidal rule method is $\{ h\lambda \in \mathbb{C} : \mathbf{Re}(h\lambda) <0 \}$ and so I would expect the method to be unstable for any choice of $\lambda >0$, but in fact I've tried it on mathlab with the test equation using $\lambda = 1$ and the approximations $\{ y_{n}\} \rightarrow 0 $ for $n \rightarrow \infty$.
Thanks so much for everyone's help. Much appreciated.

Best Answer

Stability relates only to the fact that ODE with exact solutions that converge (in some sense) or get trapped in some bounded region see this behavior replicated in the numerical solutions. ODE solutions that do not show this behavior are of no interest for the question of stability.

A first quantitative measure for this general idea can be obtained for linear test systems, where the stability can be characterized by the product of eigenvalue and step size falling into some stability region in the complex plane. This leads to the A-stability.

Note that there is a difference between stability and convergence.

Convergence is concerned with what happens if the step size goes to zero. For multi-step methods consistency and zero-stability guarantees that.
Stability is concerned with how (moderately) large one can make the step size and still get a result that is qualitatively and somewhat quantitatively correct.

You get some interplay of these concepts in numerical solvers with adaptive step size. If the system approaches an equilibrium, then at some point the error tolerances are satisfied if the solution were to stay constant. This would mean that the step size could be arbitrarily large. But in fact it is restricted by the stability region of the method, and the largest eigenvalue of the Jacobian at the stationary point. This has practically the consequence that for a large range of error tolerances the method produces the same final step sizes.

Related Solutions

[Math] Stability and convergence for Heun’s method

Stability of explicit Runge-Kutta methods around a fixed point or a slow moving equilibrium is obtained when $Lh\sim 1$ with $L$ the $y$-Lipschitz constant and $h$ the step size. The exact constant on the right side depends on the order of the method, usually it is between 1 and 5. For first and second order methods that is only a guideline as the stability region does not contain a semi-circle around the origin in the negative complex half-plane.

To get a practical feeling for this, just solve your problem for various step sizes over the integration interval from $0$ to $1/3$, as $y(1/3)=0.1e^{-50}$ is a really small number, relative to the initial value and floating point precision.

def Heun(f,t0,y0,tf,h):
    def HeunStep(t,y,h): 
        k1 = h*f(t,y); 
        k2 = h*f(t+h,y+k1); 
        return y+0.5*(k1+k2);
    t, y = t0, y0
    while t+1.1*h<tf: t, y = t+h, HeunStep(t,y,h)
    return HeunStep(t,y,tf-t)

def f(t,y): return 49-150*(t+y)
t0,y0 = 0.0,1./3+0.1
for h in np.linspace(3.5,0.5,16)/150:
    print "h=%.4e 150*h=%.4e, y(1/3)=%.12g"%(h,150*h, Heun(f,t0,y0,1./3,h))

with the results

h=2.3333e-02 150*h=3.5000e+00, y(1/3)=3382845.45109
h=2.2000e-02 150*h=3.3000e+00, y(1/3)=1820496.14682
h=2.0667e-02 150*h=3.1000e+00, y(1/3)=558708.429462
h=1.9333e-02 150*h=2.9000e+00, y(1/3)=79763.1519923
h=1.8000e-02 150*h=2.7000e+00, y(1/3)=9204.09082254
h=1.6667e-02 150*h=2.5000e+00, y(1/3)=1648.41784102
h=1.5333e-02 150*h=2.3000e+00, y(1/3)=37.6118108918
h=1.4000e-02 150*h=2.1000e+00, y(1/3)=0.740437945724
h=1.2667e-02 150*h=1.9000e+00, y(1/3)=0.00432803787737
h=1.1333e-02 150*h=1.7000e+00, y(1/3)=1.0688786734e-05
h=1.0000e-02 150*h=1.5000e+00, y(1/3)=1.14794370277e-08
h=8.6667e-03 150*h=1.3000e+00, y(1/3)=5.57823318786e-12
h=7.3333e-03 150*h=1.1000e+00, y(1/3)=2.8015784137e-15
h=6.0000e-03 150*h=9.0000e-01, y(1/3)=2.34187669257e-17
h=4.6667e-03 150*h=7.0000e-01, y(1/3)=-1.38777878078e-17
h=3.3333e-03 150*h=5.0000e-01, y(1/3)=2.60208521397e-17

As the stability region suggests, $0<h<\frac2{150}$ are valid choices in the table above, giving results close to zero. Sensible results that are practially zero as floating point numbers (and relative to the initial value) are only obtained for about $h<\frac{1.5}{150}$.

Runge-Kutta Method and Stability Function

Yes, in implicit methods you need to solve the implicit non-linear system. Preferably using a super-linear method.
Yes, the stability function covers the behavior for linear problems, with a focus on non-expanding dynamics. You get to solve $\newcommand{\ones}{{\bf 1}}$ $$k=λ(y\ones+hBk)\implies k=λ(I-λhB)^{-1}\ones y$$ and $$y_{+1}=y+hc^Tλ(I-λhB)^{-1}y=(1+(λh)c^T(I-λhB)^{-1}\ones)y$$.

Then apply the determinant identity $\det(I+uv^T)=1+v^Tu$ in reverse $$1+(λh)c^T(I-λhB)^{-1}\ones=\det(I+(λh)\ones c^T(I-λhB)^{-1})=\frac{\det(I-λhB+(λh)\ones c^T)}{\det(I-λhB)}$$

Best Answer

Related Solutions

[Math] Stability and convergence for Heun’s method

Runge-Kutta Method and Stability Function

Related Question