# Stability of Euler's Method for non-linear ODE

Consider the ODE $$y'(t) = \lambda y(t), \quad y(t_0) =y_0.$$ Euler's method $y_{i+1}=y_i+h\lambda y_i$ is stable (meaning that the solution decays or stays constant as $i \to \infty$) provided that $|1 + h\lambda| \leq 1$. This idea can be extended to systems of $n$ dimensions by placing a similar condition on the maximum coefficient instead of just $\lambda$.

This question is about non-linear, $n$ dimensional systems. Under what conditions is Euler's method stable for the system

$$\mathbf{y}'(t)=\mathbf{F}(\mathbf{y}(t)), \quad \mathbf{y}(t_0)= \mathbf{y}_0$$?

I suspect this will involve the Jacobian of the non-linear term. The only resource I've found is on slide 18 of these lecture slides, but the notation isn't explained (what is $l_{k+1}$?). It was probably introduced in a previous lecture. The integral of the Jacobian that appears here is interesting, it looks like some kind of weighted average between the numerical and the true solution. Where does this come from?

Given the linear problem $\vec{y}'=A\vec{y} + \vec{b}$, where $A$ is a diagonalizable n-dimensional matrix, you can find matrices $V,D$ s.t $AV=VD$, with $D=diag${$\lambda_1,...,\lambda_n$}.

Then $\vec{y}'=A\vec{y} + \vec{b}$ iff $\vec{z}'=D\vec{z} + V^{-1}\vec{b}$. Since we're using Euler's method, the restriction on the timestep is given by the smallest eignevalue of $A$.

If we consider now $\mathbf{y'}=\mathbf{F(y(t))}$,we can linearize $\mathbf{y'(t)})=\mathbf{F(y(t))} + J(\mathbf{y_n})(\mathbf{y}(t) - \mathbf{y_n})$, where $J$ is the jacobian matrix. If $J$ is diagonalizable, then we proceed as done in the linear case (if $J$ is not diagonalizable, then we use Jordan decomposition.

Problems arise because the eigenvalues of $J$ aren't usually enough to describe the behaviour of the exact solution.

We can have several cases: for example $J_f$ could has $\Re(\lambda_i)<0$, but the solution doesn't doesn't decay monotous for every t.

That's what happen if we consider the system $y'=Ay$, with $A=[ ( -\frac{1}{2t} , \frac{2}{t^3} ),( -\frac{t}{2} , -\frac{1}{2t} ) ]$

For more references, see Lambert J. (1991) Numerical methods for ODE systems

• Thanks for the answer! I just wanted to check some possible minor errors: (1) In the linear system, pre-multiplying by $V^{-1}$ and post-multiplying by $V$, should it be $\vec{z} = D\vec{z} + \vec{b}$? (2) When linearizing the non-linear system, should it be $\mathbf{y}'(t) = \mathbf{F}(\mathbf{y_n}) + J(\mathbf{y_n})(\mathbf{y}(t)-\mathbf{y_n})$? Feb 26, 2018 at 0:07
• You're right in (2). Check your calculation for (1): the term $V^{-1}b$ is right.
– VoB
Feb 26, 2018 at 10:42
• Gotcha: $V^{-1} \mathbf{y}'=V^{-1} A\mathbf{y}+V^{-1} \mathbf{b}$. Thanks! Feb 27, 2018 at 4:59
• I know it has been a while since you commented this, but I wanted to check something... you say step size is set by smallest eigenvalue, but I think that it ought to be largest. See these notes on linear stability analysis (slides 16-18 in particular) cs.princeton.edu/courses/archive/fall11/cos323/notes/… Does this change seem reasonable to you? May 20, 2021 at 0:58