Stability Analysis

General Theory

A group of ordinary differential equations (ODE) reads: \[ \frac{d\boldsymbol{r}}{dt} = \boldsymbol{F}(\boldsymbol{r}) \] in the vector form. Suppose that \(\boldsymbol{r^{\star}}\) is an equilibrium point. \(\boldsymbol{F}(\boldsymbol{r^{\star}})=0\). Now, we expand ODEs near the equillibrium point in Taylor series: \[ \frac{d\boldsymbol{r}}{dt} = \boldsymbol{F}(\boldsymbol{r^{\star}}) + \boldsymbol{J}(r^{\star})(\boldsymbol{r}-\boldsymbol{r^{\star}}) + \cdots \] \[ =\boldsymbol{J}(r^{\star})(\boldsymbol{r}-\boldsymbol{r^{\star}}) + \cdots \] The partial derivative above is a Jacobian matrix. Now, ODEs are linearized.From the theory of linear differential equations, the solution can be written as a superposition of terms like \(e^{\lambda_jt}\) where \(\lambda_j\) are the set of Jacobian's eigenvalues. Standard analysis result in the following theorem:

An equilibrium point \(\boldsymbol{r^{\star}}\) of ODEs is stable if all the eigenvalues of \(\boldsymbol{J}(\boldsymbol{r^{\star}})\), the Jacobian evaluated at \(\boldsymbol{r^{\star}}\), have negative real parts. The equilibrium point is unstable if at least one of the eigenvalues has a positive real part.

Definitely, if we want to see reaction oscillation, the equilibrium point must be unstable.

Local instability of FKN model

In FKN model \[ \boldsymbol{F} = \begin{bmatrix} \frac{qy - xy + x - x^2}{\epsilon} \\ \frac{-qy-xy+2fz}{\delta} \\ x-z \end{bmatrix} \quad \boldsymbol{r} = \begin{bmatrix}x \\ y \\ z\end{bmatrix} \] Generally speaking, all parameters must be positive. And Jacobian reads: \[ \boldsymbol{J}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{r}}=\begin{bmatrix} \frac{- 2 x - y + 1}{\epsilon} & \frac{q - x}{\epsilon} & 0 \\ -\frac{y}{\delta} & \frac{- q - x}{\delta} & \frac{2 f}{\delta} \\ 1 & 0 & -1 \end{bmatrix} \] There are two positive solutions of \(\boldsymbol{F}(\boldsymbol{r})=0\) The first one is trivial: \(\boldsymbol{r}^{\star} = \begin{bmatrix} 0 & 0&0\end{bmatrix}^T\). The determinant equation around this solution reads:

\[ - \lambda^{3} + \lambda^{2} \left(-1 + \frac{1}{\epsilon} - \frac{q}{\delta}\right) + \lambda \left(\frac{1}{\epsilon} - \frac{q}{\delta} + \frac{q}{\delta \epsilon}\right) + \frac{2 f q}{\delta \epsilon} + \frac{q}{\delta \epsilon} = 0 \]

This is a third order polynominal with positive zero order term. It is easy to find there must exist a positive real solution. Therefore, the trivial solution is always unstable. The other one is more interesting. The equilibrium point reads:

\[ x^{\star}=z^{\star}=- f - \frac{q}{2} + \frac{\sqrt{4 f^{2} + 12 f q - 4 f + q^{2} + 2 q + 1}}{2} + \frac{1}{2} \] \[ y^{\star} = \frac{2fx^{\star}}{q+x^{\star}} \] Instead of tuning all parameters, only \(f\) change and others are set as suggested values in [1]. To see oscillations, \(f\) must fall between 0.25 and 1.2 approximately which is close to the analytical result in [1].

Numerical experiments

To verify linear instability analysis, several numerical experiments were carried out. It is clear to see that if \(f\) is out estimated interval, the solution is approaching the equilibrium point. Otherwise, oscillation happens.

1: Murray, J.D. (2002). Mathematical Biology I: An Introduction. New York, NY: Springer Science+Business Media, Inc.