Linear Stability Analysis
The system of differential equations is linearized using perturbations. For simplicity, we consider the following two-dimensional system,
(1)\[\begin{split}\frac{x_1(t)}{dt} &= f_1(x_1, x_2) \\
\frac{x_2(t)}{dt} &= f_2(x_1, x_2).\end{split}\]
The equilibrium point \(x_{1,0}, x_{2,0}\) is determined by
\[\begin{split}f_1(x_{1,0}, x_{2,0}) &= 0 \\
f_2(x_{1,0}, x_{2,0}) &= 0.\end{split}\]
We perturb the equations around its equilibrium point,
\[\begin{split}x_1(t) &\to x_{1,0} + \delta x_1(t) \\
x_2(t) &\to x_{2,0} + \delta x_2(t),\end{split}\]
where \(\delta x_1(t)\) and \(\delta x_2(t)\) are small quantities.
The equation (1) becomes
(2)\[\begin{split}\frac{d\delta x_1(t)}{dt} &= f_1(x_{1,0} + \delta x_1(t), x_{2,0} + \delta x_2(t)) \\
\frac{d\delta x_2(t)}{dt} &= f_2(x_{1,0} + \delta x_1(t), x_{2,0} + \delta x_2(t)).\end{split}\]
We apply Taylor series expansion to \(f_1(x_{1,0} + \delta x_1(t), x_{2,0} + \delta x_2(t))\),
\[\begin{split}&f_1(x_{1,0} + \delta x_1(t), x_{2,0} + \delta x_2(t)) \\
=& \left.\frac{d f_1(x_1, x_2)}{d x_1} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_1 + \left.\frac{d f_1(x_1, x_2)}{d x_2} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_2 + \mathscr{O},\end{split}\]
where we have used \(f_1(x_{1,0}, x_{2,0}) = 0\). The perturbed equation (2) is linearized
(3)\[\begin{split}\frac{d\delta x_1(t)}{dt} &= \left.\frac{d f_1(x_1, x_2)}{d x_1} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_1 + \left.\frac{d f_1(x_1, x_2)}{d x_2} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_2 \\
\frac{d\delta x_2(t)}{dt} &= \left.\frac{d f_2(x_1, x_2)}{d x_1} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_1 + \left.\frac{d f_2(x_1, x_2)}{d x_2} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}} \delta x_2.\end{split}\]
The equation (3) can be rewritten into a matrix form
(4)\[\begin{split}\frac{d}{dt} \begin{pmatrix}
\delta x_1(t) \\
\delta x_2(t)
\end{pmatrix} = \left. \begin{pmatrix}
\frac{d f_1(x_1, x_2)}{d x_1} & \frac{d f_1(x_1, x_2)}{d x_2} \\
\frac{d f_2(x_1, x_2)}{d x_1} & \frac{d f_2(x_1, x_2)}{d x_2}
\end{pmatrix} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}}\begin{pmatrix}
\delta x_1(t) \\
\delta x_2(t)
\end{pmatrix}.\end{split}\]
For simplicity, we define the Jacobian matrix
\[\begin{split}\mathbf J = \left. \begin{pmatrix}
\frac{d f_1(x_1, x_2)}{d x_1} & \frac{d f_1(x_1, x_2)}{d x_2} \\
\frac{d f_2(x_1, x_2)}{d x_1} & \frac{d f_2(x_1, x_2)}{d x_2}
\end{pmatrix} \right\vert_{x_1=x_{1,0}, x_2=x_{2,0}},\end{split}\]
so that the matrix form of the linearized equation (4) becomes,
\[\begin{split}\frac{d}{dt} \begin{pmatrix}
\delta x_1(t) \\
\delta x_2(t)
\end{pmatrix} = \mathbf J \begin{pmatrix}
\delta x_1(t) \\
\delta x_2(t)
\end{pmatrix}.\end{split}\]
To investigate the stability of the differential system, we assume that
\[\begin{split}\begin{pmatrix}
\delta x_1(t) \\
\delta x_2(t)
\end{pmatrix} = \begin{pmatrix}
\delta x_1(t_0) \\
\delta x_2(t_0)
\end{pmatrix} e^{\lambda t},\end{split}\]
which leads to linear equations
\[\begin{split}\begin{pmatrix}
\delta x_1(t_0) \\
\delta x_2(t_0)
\end{pmatrix} \lambda = \mathbf J \begin{pmatrix}
\delta x_1(t_0) \\
\delta x_2(t_0)
\end{pmatrix}.\end{split}\]
For non-trivial solutions, we require
\[\operatorname{Det}(\mathbf J - \lambda \mathbf I) = 0,\]
which is also the eigen value problem of the Jacobian matrix. We expand the determinant
(5)\[\lambda^2 - (J_{11} + J_{22}) \lambda + (J_{11}J_{22} - J_{12}J_{21}) = 0.\]
For real positive solutions \(\lambda>0\), we get an exponential growing result for the linearized equation. Any deviation from the equilibrium point leads to a run-away process and the system moves further away from the equilibrium point. For real negative solutions \(\lambda < 0\), the system will move back to the equilibrium point given any deviations from the equilibrium. Imaginary solutions of \(\lambda\) leads to oscillations.