Possible ways to do stability analysis of non-linear, three-dimensional Differential Equations For example Lorenz system,
$$
\frac{d}{dt}\begin{pmatrix}
x\\ 
y\\ 
z
\end{pmatrix}=\begin{pmatrix}
-\sigma & \sigma & 0\\ 
\rho & -1 & -x\\ 
y & 0 & -\beta
\end{pmatrix}\begin{pmatrix}
x\\ 
y\\ 
z
\end{pmatrix}
$$
The article on Wiki describes the system of differential equations to have stable equilibrium for a specific value of $\rho$.
Given some three dimensional system of differential equations, I would like to know how to do the stability analysis? or explain how stability analysis of Lorenz system is done ...
 A: Consider the Lorenz equations:
$$x' = f(x,y,z) = \sigma(y-x)$$
$$y' = g(x,y,z) = \rho x - y -xz$$
$$z' = h(x,y,z)= -\beta z + xy$$
where $\sigma$ and $\beta$ are viewed as fixed positive constants and $\rho$ is a parameter.
Here is a guide to work through given that the Wiki describes what happens at each of the critical points. The outline of how to derive those is as follows:


*

*$(1)$ Find the critical points of the system, that is, where you simultaneously have $x' = y' = z' = 0$.


Clearly, a critical point is given by $(x, y, z) = (0,0,0)$.


*

*$(2)$ Find the Jacobian matrix of the system of equations.


The Jacobian matrix (using partial derivatives) is given by:
$$\tag 2 \displaystyle J = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y}& \frac{\partial f}{\partial z} \\\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}& \frac{\partial g}{\partial z}\\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y}& \frac{\partial h}{\partial z}\end{bmatrix} = \begin{bmatrix}-\sigma & \sigma & 0\\\rho - z & -1 & -x\\ y & x & -\beta\end{bmatrix}$$


*

*$(3)$ Lastly, evaluate the eigenvalues of the Jacobian evaluated at each of the critical points for stability analysis. So, evaluating $(2)$ at the CP $(0,0,0)$, yields the matrix:


$$\displaystyle J_{0,0,0} = \begin{bmatrix}-\sigma & \sigma & 0\\\rho & -1 & 0\\ 0 & 0 & -\beta\end{bmatrix}$$
To find the eigenvalues, we set-up and solve for the characteristic polynomial using:
$$| A - \lambda I| = 0 \rightarrow -\beta \lambda^2-\beta \lambda+\beta r \sigma-\beta \lambda \sigma-\beta \sigma-\lambda^3-\lambda^2+\lambda \rho \sigma-\lambda^2 \sigma-\lambda \sigma = 0$$
This leads to the three eigenvalues (note, we could have written these straight off given the special form of the diagonal matrix above) and their corresponding eigenvectors as:
$\displaystyle \lambda_1 = -\beta, v_1 = (0, 0, 1)$
$\displaystyle \lambda_2 = \frac{1}{2} \left(-\sqrt{4 \rho \sigma+\sigma^2-2 \sigma+1)}-\sigma-1\right), v_2 = -\frac{-1+\sigma+\sqrt{1-2 \sigma+4 \rho \sigma+\sigma^2}}{2 \rho}, 1, 0)$ 
$\displaystyle \lambda_3 = \frac{1}{2}\left(\sqrt{4 \rho \sigma+\sigma^2-2 \sigma+1)}-\sigma-1\right), v_3 = -\frac{-1+\sigma-\sqrt{1-2 \sigma+4 \rho \sigma+\sigma^2}}{2 \rho}, 1, 0)$
Of course those $\lambda$ are the eigenvalues we use to do a stability analysis.
To study the stability of this, we need some complex center manifold theory, so I am going to stay away from that. It leads to a pitchfork bifurcation at the origin.
It is worth noting that there are two other pairs of fixed points at:


*

*$x = y = \pm \sqrt{\beta(\rho - 1)}, z = \rho -1$


To find this, notice that the first equation gives us $y=x$, substituting this into the second equation gives $x(\rho - 1) - xz = 0 \rightarrow z = \rho-1$, and we substitute these previous two values into the third and have $x^2 = \beta z \rightarrow x =  \pm \sqrt{\beta(\rho - 1)}$. 
To simplify matters, lets linearize the system by removing the nonlinear terms and analyze the simpler system.
The linearization (eliminate the $xy$ and $xz$ nonlinear terms from the system) is given by:
$$x' = \sigma(y-x)$$
$$y' = \rho x - y$$
$$z' = -\beta z$$
If you look at the equation for $z'$, you notice that it is decoupled from the other two equations and it can be solved outright as $z(t) = c_1 e^{-\beta t}$ and $z(t) \rightarrow 0$ exponentially fast. The other two directions are governed by the system:
$$\begin{bmatrix}x'\\y'\end{bmatrix} = \begin{bmatrix}-\sigma & \sigma\\\rho & -1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}$$
You should be able to analyze this system from all the data provided now.
As you can see, these systems are very complex and include a rich set of mathematics to deal with, so it takes time to get there.
Here is another guide with some of the details for you to work out.
