When checking if a numerical ODE solving method (euler,runge-kutta etc) is stable, we use the problem $y'=\lambda y$, and find the step size $h$ such that our method converges to zero as the number of time steps $n\to\infty$.


Intuitively, why do we use this ODE to find the stable region? Why don't we use another ODE?


This answer is highly inspired by:

Desmond J. HighamLloyd N. Trefethen, Stiffness of ODEs

Stability analysis is performed to study whether the numerical solution of the ODE $$\mathbf{y}' = \mathbf{A}(t, \mathbf{y}) + \mathbf{g}(t),\quad \mathbf{y}(0) = \mathbf{y}_0\in\mathbb{R}^d$$ blows up or stays bounded. To perform such an analysis, the ODE is simplified according to the following steps:

  1. The ODE is linearized, that is, you consider instead a problem of the form $$\mathbf{u}' = \mathbf{B}(t) \mathbf{u},\quad \mathbf{u}(0) = \mathbf{u}_0\in\mathbb{R}^d\,,$$ where $\mathbf{B}$ is the Jacobian of $\mathbf{A}$ at $\mathbf{y}_0$.
  2. The Jacobian is frozen that is, you consider instead $$\mathbf{u}' = \tilde{\mathbf{B}} \mathbf{u},\quad \mathbf{u}(0) = \mathbf{u}_0\in\mathbb{R}^d\,,$$ where $\tilde{\mathbf{B}}=\mathbf{B}(0)$.
  3. You assume that $\tilde{\mathbf{B}}$ is diagonalizable, that is, that there is an invertible matrix $\mathbf{S}$ such that $\mathbf{S}\tilde{\mathbf{B}}\mathbf{S}^{-1} = \mathbf{D}$, where $\mathbf{D}$ is a diagonal matrix containig the eigenvalues of $\tilde{\mathbf{B}}$. In this case, the function $\mathbf{z} = \mathbf{S}\mathbf{u}$ satisfies $$\mathbf{ z}' = \mathbf{D} \mathbf{z},\quad \mathbf{z}(0) = \mathbf{S}\mathbf{u}_0\in\mathbb{R}^d\,,$$ which is a decoupled system of ODEs of the form $y' = \lambda y$.

If the method you're using is stable for $y' = \lambda y$ for all $\lambda$ that are on the diagonal of $\mathbf{D}$, then you can use the method to approximate $\mathbf{z}$, and if your method is affine invariant (usually it is), you can compute $\mathbf{u}(T)$ at a final time $T$ by computing $\mathbf{z}(T)$ and multiplying it by $\mathbf{S}^{-1}$. This is why you consider $y' = \lambda y$ as a test case.

The limits of steps 1. and 2. are discussed in the paper I mentioned above.


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.