Linear Stability Analysis of Euler’s Method

nonlinear dynamicsnumerical methodsstability-in-odes

I am somewhat familiar with linearization of non-linear ODEs for local stability analysis.

It's been explained to me like so. Let $\eta(t) = x(t) – x_0$ represent a small change of the solution away from $x_0$.

Then we have:

$$x(t) = \eta(t) + x_0$$

and

$$\dot{x} = \dot{\eta} = f(\eta(t) + x_0)$$

Then if we expand $f(\eta(t) + x_0)$ about $x_0$:

$$f(\eta(t) + x_0) = f(x_0) + f'(x_0)\eta(t) + O(\eta^2)$$

Now here, in my dynamics class we were referring to $x_0$ as not just any point, but specifically one that is a fixed point or equilibrium. In that case the first term disappears and neglecting higher order terms we get:

$$f(x) = \dot{\eta} = \eta f'(x_0)$$

Which is an ODE in $\eta$. We see by solving it that we can determine stability by the sign of $f'(x_0)$. If it is negative we have a decaying exponential, and if it is positive, the solution tends to infinity.

Ok, so with that preamble out of the way, I'm reading some notes I found online from a numerical methods class, that seem to be applying this exact line of thinking to the stability of a numerical method. I'm trying to find the connection:

If a numerical method is stable in the above sense for a certain range of
values of λ, then it is possible to show that it will be stable for the ODE
$y' = f(t, y)$ as long as $\frac{\partial f}{\partial y}$ is in that range of λ (and f is smooth enough). We won’t prove this theorem here.

Let us consider a few examples.

Example 14. For the forward Euler method applied to $y_0 = λy$, we get
$y_n+1 = y_n + hλy_n = (1 + hλ)y_n$.
The iterates $y_n$ tend to zero provided $|1 + hλ| < 1$, where the | · | denote
the complex modulus. Write $hλ = x + iy$, so that the inequality becomes
$(1 + x)^2 + y^2 < 1$. This is the equation of a disc in the complex (hλ)-plane,
with center at −1, and radius 1. If hλ sits inside this disk, the method
is stable, and otherwise it isn’t. We say that the forward Euler method is
conditionally stable: typically, we require both Re(λ) < 0 and a small step
size h in order to guarantee stability.

Question 1: It appears to me that the logic behind $\lambda$ here is arrived at with the same thought process as in my example above. Is this true?

Question 2: When applying Euler's method to the ODE $y' = \lambda y$, how is it this statement true: "The iterates $y_n$ tend to zero provided $|1 + hλ| < 1$"

I'm afraid the answer to question 2 is something so easy and obvious, and I hate to ask, but I just can't see it.

Best Answer

Regarding question 1, the answer is yes, it is the same reasoning that motivates this stability definition. Regarding question 2, you just need to write down what would be Euler's method applied to the ode $y'=\lambda y, y(0)=1$:

$$ y_{n+1} = y_n + h f(t_n,y_n) = y_n + h \lambda y_n = (1+\lambda h) y_n. $$

In this case you can get the explicit formula $y_n= (1+\lambda h)^n$, and you conclude that $y_n$ converges to zero iif $|1+\lambda h| < 1$.

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}$.

[Math] Stability region of Heun’s method

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

Best Answer

Related Solutions

[Math] Stability and convergence for Heun’s method

[Math] Stability region of Heun’s method

Related Question