5

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

Question:

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

MeowBlingBling
  • 540

1 Answers1

6

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.

Albe
  • 465