I'm working on a prey-predator model. I'm using the following system of differential equations for it:
\begin{align} x'&=-a_1x+a_2xy+a_3xz\\
y'&=b_1y-b_2xy\\
z'&=c_1z-c_2xz
\end{align}
Where $a_i, b_i, c_i >0$.
One of the stationary points is $P=(\frac{b_1}{b_2},\frac{a_1}{a_2},0)$.
Question:
How can I determine the stability of this point $P$?
Attempt:
First I wrote the equation as:
\begin{align} \frac{\mathrm d \underline{v}}{\mathrm d t}=\underline{F}(\underline{v}), \hspace{10pt} \underline{v}=(x,y,z)\end{align}
I looked at the linearized ODE, and I found:
\begin{align} \frac{\mathrm d \underline{\hat{v}}}{\mathrm t}=
\begin{pmatrix}
0 & \dfrac{a_2b_1}{b_2} & \dfrac{a_3b_1}{b_2}\\
-\dfrac{a_1b_2}{a_2} & 0 & 0 \\
0 & 0 &c_1-\dfrac{c_2b_1}{b_2}\\
\end{pmatrix}
\underline{\hat{v}}
\end{align}
The eigenvalues are $c_1-\dfrac{c_2b_1}{b_2}, \pm i \sqrt[]{a_1b_1}$. My problem is the pair with zero real part. I have learned a theorem that only says something when all eigenvalues have negative real part (then it's stable) or at least one has positive eigenvalue (then it's unstable). My other approach was Lyapunov's theorem. I found (with a little bit puzzling) the following Lyapunov function:
\begin{align}V(x,y,z) = a_2y+ b_2x -a_1\left(1+\ln (y)-\ln\left(\frac{a_1}{a_2}\right)\right)-b_1\left(1+\ln(x)-\ln\left(\frac{b_1}{b_2}\right)\right) \end{align}
Differentiating it brought me eventually to this:
\begin{align} V'=a_3b_2xz-a_3b_1z\end{align}
Now the only thing that I can say is that $V'$ is both positive and negative in every neighborhood of $P$, so the theorem doesn't say something about that case.
What I also thought about is to Taylor approximate the function $F$ till the second order term, like it can be done with one-dimensional ODE if the first derivative is equal to zero. But then here I get a matrix in a matrix (a tensor?) which is pretty vague. I can not visualise what is going on there, but if someone can explain how to do it that way, then you are welcome. I have read many articles about this problem on internet. I couldn't find something useful/understandable. Some are talking about manifolds, but I have not learned that yet. What I understand is that it is something like a solution curve.
I really couldn't sleep well because of this problem. Your help is appreciated! Thanks in advance.
Update
I have realised that the function $V$ that I have used missed an imortant requirment, it was $0$ in all points $(\frac{b_1}{b_2},\frac{a_1}{a_2},z)$, so that was not a Lyapunov function. My new approach was to pick:
\begin{align} V(x,y,z)=\left(x-\frac{b_1}{b_2}\right)^2+\left( y-\frac{a_1}{a_2}\right)^2+z^2\end{align}
Okay, this one is Lyapunov for sure. I differentiate w.r.t. $t$ and get:
\begin{align} V'=(-a_1x+a_2xy+a_3xz)\left(x-\frac{b_1}{b_2}\right) + (b_1y-b_2xy)\left( y-\frac{a_1}{a_2}\right) +z(c_1z-c_2xz)
\end{align}
This is zero in $(x,y,z)=P$ but that is also a saddle point, which means that is both negative and positive in every neighborhood of $P$. I also tried some Lyapunov functions of the form:
\begin{align} V(x,y,z)=f(x-\frac{b_1}{b_2})+g(y-\frac{a_1}{a_2})+h(z)\end{align}
where $f,g,z$ are all even functions that has minimum in zero. Examples that I have tried are: $-\cos(x), -e^{-x^2}$ and they also made $V'$ have a saddle point in $P$.
I also transformed everything to cilindar coordinates, but the equation that I have got was not so beautiful. I believe there is a way to prove the stability of this point.
Can you guys help me? I really don't know what to do after this.
Best Answer
Heh, it's actually simpler than expected. The system is stable but not asymptotically stable. Here is a rough sketch of a proof:
First, for simplicity, set \begin{align*} \beta &= \frac{b_1}{b_2}\\ \alpha&= \frac{a_1}{a_2}\\ \end{align*} and then center your system at the equilibrium point $P = \left(\frac{b_1}{b_2}, \frac{a_1}{a_2}, 0 \right) = \left(\beta, \alpha, 0 \right)$ which means shift the variables as follows \begin{align*} x &\mapsto x + \frac{b_1}{b_2}\\ y &\mapsto y + \frac{a_1}{a_2}\\ z &\mapsto z \end{align*} The system becomes \begin{align*} \dot{x} &= \beta a_2 \, y + \beta a_3 \, z + a_2 \, xy + a_3 \, xz\\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy\\ \dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz \end{align*} which could also be reorganized as \begin{align*} \dot{x} &= \beta a_2 \, y + a_2 \, xy + a_3(\beta + x) z\\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy\\ \dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz \end{align*} Now observe that the $x,y$-coordinate plane, also given by the equation $z=0$, is invariant under the flow of the system, meaning that the vector field is tangent to $z=0$. You can even check directly that if you take a solution $\big(x(t), y(t) \big)$ of the 2D system \begin{align*} \dot{x} &= \beta a_2 \, y + a_2 \, xy \\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy\\ \end{align*} then $\big(x(t), \, y(t), \, 0 \big)$ is a solution of the original system \begin{align*} \dot{x} &= \beta a_2 \, y + a_2 \, xy + a_3(\beta + x) z\\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy\\ \dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz \end{align*} All of this yields that the plane $z=0$ is a center manifold of your system and the restriction of the system to that invariant manifold is \begin{align*} \dot{x} &= \beta a_2 \, y + a_2 \, xy \\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy\\ \end{align*} As you already know, the latter is a conservative system (i.e. integrable, has a first integral) Thus, all trajectories on $z=0$ are ovals (closed orbits, cycles) and consequently, all orbits close enough to the invariant plane $z=0$ and in a neighborhood of the equilibrium point $P$, will tend to approach it more and more, as long as, the restriction $c_1 < \frac{b_1}{b_2} c_2$ is satisfied. Thus, they will follow the cyclic closed orbits on the invariant manifold. Hence, the equilibrium point is stable, but not asymptotically.
Observe that whenever $c_1 > \frac{b_1}{b_2} c_2$ the equilibrium $P$ is automatically unstable.
Here is the rigorous proof. The following theorem is a direct corollary from the theory of center manifolds (see for example these slides). One can also prove it, although it requires some techniques from analysis. Notice that in this case, the center manifold is known explicitly and is $z=0$.
Theorem. There exist two small open neighborhood $U, \, V \, \in \, \mathbb{R}^2$ of the point $0 = (0,0,0)$ and a smooth change of variables $$H \, : \, U \, \to \, V$$ $$H \, : \, (u,v,w) \, \mapsto \, (x,y,z)$$ of the form \begin{align*} u &= x + h_1(x,y,z)\\ v &= y + h_2(x,y,z)\\ w &= z \end{align*} where the functions $h_1$ and $h_2$ are smooth of order at least two, so that $H$ transforms the system of ODEs \begin{align} \dot{u} &= \beta a_2 \, v + a_2 \, uv\\ \dot{v} &= - \alpha b_2\, u - b_2 \, uv \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ (1)\\ \dot{w} &= \left(c_1 - \beta c_2\right) w - g(u,v)w \end{align} into the original system of ODEs \begin{align} \dot{x} &= \beta a_2 \, y + a_2 \, xy + a_3(\beta + x) z \\ \dot{y} &= - \alpha b_2\, x - b_2 \, xy \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ (2)\\ \dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz \end{align}
Since system (1) is definitely stable, but not asymptotically stable, system (2) is also stable but not asymptotically because (2) is obtained from (1) by the change of variables $H$ and any change of variables preserves stability (which is a topological property).
System (1) is stable because the first two equations are separate from the third. If you pick their first integral, let's call it $J(u,v)$, one can check directly that the function $J(u,v)$ is also a first integral of (1). Then, it is easy to prove the stability of (1). On each level surface $J(u,v) = const$, a general solution $(u(t), v(t), w(t))$ close to $w=0$, which starts from a point on the surface $J(u,v) = const$, stays on the surface $J(u,v) = const$ and because of the coefficient $c_1 - \beta c_2 < 0$ it converges to the solution $(u(t), v(t), 0)$, simply because $\lim_{t \to \infty} \, w(t) = 0$. However,
$(u(t), v(t), 0)$ is actually a closed cycle (a periodic solution). Therefore system (1) is stable, but not asymptotically, because $(u(t), v(t), 0)$ does not converge to the origin.