NOVEL STABILITY CONDITIONS FOR SOME GENERALIZATION OF NICHOLSON’S BLOWFLIES MODEL WITH STOCHASTIC PERTURBATIONS

Abstract

We consider a generalization of the well-known nonlinear Nicholson blowflies model with stochastic perturbations. Stability in probability of the positive equilibrium of the considered equation is studied. Two types of stability conditions: delay-dependent and delay-independent conditions are obtained, using the method of Lyapunov functionals and the method of linear matrix inequalities. The obtained results are illustrated by numerical simulations by means of some examples. The results are new, and complement the existing ones.

1 Introduction

The Nicholson blowflies model,

(1.1) $$ \begin{align} \dot x(t)=-ax(t)+kx(t-h)e^{-\lambda x(t-h)}, \end{align} $$

is primarily developed to describe the periodic oscillation in Nicholson’s classic experiments [12], where k denotes the maximum daily egg production rate per capita. The size at which the blowfly population reproduces at its maximum rate is denoted by $1/\lambda $ , and the per capita daily adult death rate and the generation time are denoted by a and h, respectively. This model and its different variations have received a lot of attention and are very popular in research, especially in the qualitative theory of differential equations, due to their applicability. For more details, we refer the reader to [2, 3, 5, 7–10, 15–17] and references therein.

In this work, we consider a very general framework for a model which includes discrete and distributed delays as particular cases. Also, since the randomness is quite evident in the real world, we consider the model with stochastic perturbations. These incorporations make the model very realistic and genuine.

Consider the integro-differential equation

(1.2) $$ \begin{align}\begin{aligned} \dot x(t)&=-ax(t)-bx(t-\tau)+\int^\infty_0x(t-s)e^{-\lambda(s)x(t-s)}\, dK(s),\quad t\ge0, \\ x(s)&=\phi(s), \quad s\le0,\end{aligned} \end{align} $$

where $\lambda (s)>0$ , $1/\lambda (s)$ is the size at which the blowfly population reproduces at its maximum rate at the time moment s, the kernel $K(s)$ is a nondecreasing function of bounded variation on $[0,\infty )$ , that is, $K=\int ^\infty _0dK(s)<\infty $ , and the integral is understood in the Stieltjes sense [14]. In particular, this means that both distributed and discrete delays can be used depending on the concrete choice of the kernel $K(s)$ . For instance, in the case when $b=0$ , $\lambda (s)=\lambda $ and $dK(s)=k\delta (s-h)ds$ , where $k>0$ and $\delta (s)$ is Dirac’s function [16], we obtain the classical Nicholson’s blowflies model (1.1).

It is clear that equation (1.2) has the zero equilibrium, and its positive equilibrium $x^*$ is a solution of the equation

(1.3) $$ \begin{align} \int^\infty_0e^{-\lambda(s)x^*}dK(s)=a+b<K. \end{align} $$

Remark 1.1. Let

$$ \begin{align*} dK(s)=\sum^m_{i=1}k_i\delta(s-h_i)\,ds, \quad \lambda(h_i)=\lambda_i. \end{align*} $$

In this case, $K=\sum \nolimits ^m_{i=1}k_i$ and the equation (1.3) takes the form

$$ \begin{align*} \sum^m_{i=1}k_ie^{-\lambda_ix^*}=a+b<\sum^m_{i=1}k_i. \end{align*} $$

If, in particular, $m=1$ , $\lambda _1=\lambda>0$ , then $K=k_1$ and

(1.4) $$ \begin{align} x^*=\frac{1}{\lambda}\ln\frac{K}{a+b}>0. \end{align} $$

If $m=2$ , then the equilibrium $x^*$ is defined by the equation

(1.5) $$ \begin{align} k_1e^{-\lambda_1x^*}+k_2e^{-\lambda_2x^*}=a+b<k_1+k_2, \end{align} $$

which, generally speaking, can be solved numerically but sometimes analytically too.

Example 1.2. Let $m=2$ , $k_1=\lambda _1=a+b=1$ , $k_2=\lambda _2=2$ . Then $x^*=\ln 2$ .

Suppose that equation (1.2) is exposed to stochastic perturbations, which are of the type of white noise and are directly proportional to the deviation of the system state $x(t)$ from the equilibrium $x^*$ and influence $\dot x(t)$ immediately. So, equation (1.2) is transformed into the Ito stochastic differential equation [6],

(1.6) $$ \begin{align} dx(t)=(-ax(t)-bx(t-\tau)+\int^\infty_0x(t-s)e^{-\lambda(s)x(t-s)}\,dK(s))\,dt+\sigma(x(t)-x^*)\,dw(t), \end{align} $$

where $\sigma $ is a constant and $w(t)$ is the standard Wiener process on a complete probability space $\{\Omega , \mathfrak F, \mathbf P\}$ . For more details on this topic, we refer the reader to [6, 16].

Note that the equilibrium $x^*$ of equation (1.2) is also a solution of equation (1.6).

2 Centralization and linearization

To centralize equation (1.6) around the equilibrium $x^*$ , put $x(t)=y(t)+x^*$ . Substituting this into (1.6) gives

(2.1) $$ \begin{align} dy(t)=(-ay(t)-by(t-\tau)+I(y_t))\,dt+\sigma y(t)\,dw(t), \end{align} $$

where $y_t$ is the trajectory of $y(s)$ , $s\le t$ , and

$$ \begin{align*} I(y_t)=&\int^\infty_0(y(t-s)+x^*)e^{-\lambda(s)(y(t-s)+x^*)}\,dK(s)-(a+b)x^*. \end{align*} $$

Note that, via (1.3), $I(0)=0$ and stability of the equilibrium $x^*$ of equation (1.6) is equivalent to stability of the zero solution of equation (2.1). Using the representation

$$ \begin{align*}e^{-\lambda y}=1-\lambda y+o(y),\quad \text{where}\ o(y)\text{ means that}\quad \lim_{y\to0}\dfrac{o(y)}{y}=0, \end{align*} $$

and (1.3), we transform $I(y_t)$ as follows.

(2.2) $$ \begin{align} I(y_t)&=\int^\infty_0(y(t-s)+x^*)e^{-\lambda(s)y(t-s)}e^{-\lambda(s)x^*}dK(s)-(a+b)x^* \nonumber \\ &=\int^\infty_0y(t-s)e^{-\lambda(s)x^*}dK(s)+x^*\int^\infty_0(1-\lambda(s)y(t-s))e^{-\lambda(s)x^*}dK(s) \nonumber\\ &\quad +o(y)-x^*\int^\infty_0e^{-\lambda(s)x^*}dK(s) \nonumber \\ &=\int^\infty_0y(t-s)F(s,x^*)\,dK(s)+o(y), \end{align} $$

where

(2.3) $$ \begin{align} F(s,x^*)=(1-\lambda(s)x^*)e^{-\lambda(s)x^*}. \end{align} $$

Substituting (2.2) into (2.1) and neglecting $o(y)$ , we obtain the linear equation

(2.4) $$ \begin{align} dz(t)=\bigg(-az(t)-bz(t-\tau)+\int^\infty_0z(t-s)F(s,x^*)\,dK(s)\bigg)+\sigma z(t)\,dw(t), \end{align} $$

which is the linear part of equation (2.1).

Remark 2.1. Note that, for equation (1.1) $b=0$ , via (2.3) and (1.4),

$$ \begin{align*} F(x^*)=(1-\lambda x^*)e^{-\lambda x^*} =\bigg(1-\ln\frac{k}{a}\bigg)\frac{a}{k}, \end{align*} $$

and equation (2.4) takes the form

(2.5) $$ \begin{align} dz(t)=\bigg(-az(t)-a\bigg(\ln\frac{k}{a}-1\bigg)z(t-h)\bigg)\,dt+\sigma z(t)\,dw(t). \end{align} $$

3 Stability

Definition 3.1. The zero solution of equation (2.1) is called stable in probability if, for any $\varepsilon _1>0$ and $\varepsilon _2 \in (0,1)$ , there exists $\delta>0$ such that the solution $y(t,\phi )$ of equation (2.1) satisfies the condition $\mathbf P\{\sup \nolimits _{t\ge 0}|y(t,\phi )|>\varepsilon _1\}<\varepsilon _2$ for any initial function $\phi $ such that $\mathbf P\{\sup \nolimits _{s\in [-\tau ,0]}|\phi (s)|<\delta \}=1$ .

Definition 3.2. Let $\mathbf E$ be the expectation. The zero solution of equation (2.4) is called:

• • mean square stable if, for each $\varepsilon>0$ , there exists a $\delta>0$ such that $\mathbf E|z(t,\phi )|^2<\varepsilon $ , $t\ge 0$ , provided that $\sup \nolimits _{s\in [-\tau ,0]}\mathbf E|\phi (s)|^2<\delta $ ; and

• • asymptotically mean square stable if it is mean square stable and $\lim \nolimits _{t\to \infty } \mathbf E|z(t,\phi )|^2=0$ for each initial function $\phi $ .

Remark 3.3. It is known [16] that a condition of asymptotic mean square stability of the zero solution of the linear equation (2.4) at the same time is a condition for stability in probability of the zero solution of the nonlinear equation (2.1) and, therefore, is a condition for stability in probability of the equilibrium $x^*$ of the initial nonlinear equation (1.6).

The following notation is used below, that is,

(3.1) $$ \begin{align} \alpha_i=\int^\infty_0s^i|F(s,x^*)|\,dK(s), \quad i=0,1. \end{align} $$

3.1 Delay-independent stability condition

Theorem 3.4. If

(3.2) $$ \begin{align} \rho=a-\alpha_0-\tfrac{1}{2}\sigma^2>|b|, \end{align} $$

then the equilibrium $x^*$ of equation (1.6) is stable in probability.

Proof. Via Remark 3.3, it is enough to prove that, by condition (3.2), the zero solution of equation (2.4) is asymptotically mean square stable. Using the general method for construction of Lyapunov functionals [15, 16], we construct a Lyapunov functional $V(z_t)$ in the form $V(z_t)=V_1(z(t))+V_2(z_t)$ , where $V_1(z(t))=z^2(t)$ and the additional functional $V_2(z_t)$ is chosen below.

Let L be the generator [6, 16] of equation (2.4). Then,

(3.3) $$ \begin{align} LV_1(z(t))&=2z(t)\bigg(-az(t)-bz(t-\tau)+\int^\infty_0z(t-s)F(s,x^*)\,dK(s)\bigg)+\sigma^2z^2(t) \nonumber\\ &=(-2a+\sigma^2)z^2(t)-2bz(t)z(t-\tau)+J(z_t), \end{align} $$

where $F(s,x^*)$ is as defined in (2.3) and

(3.4) $$ \begin{align} J(z_t)&=2\int^\infty_0z(t)z(t-s)F(s,x^*)\,dK(s) \nonumber\\ &\le \int^\infty_0(z^2(t)+z^2(t-s))|F(s,x^*)|\,dK(s) \nonumber\\ &=\alpha_0z^2(t)+\int^\infty_0z^2(t-s)|F(s,x^*)|\,dK(s). \end{align} $$

Now, consider the additional functional $V_2(z_t)$ in the form

$$ \begin{align*} V_2(z_t)=\rho\int^t_{t-\tau}z^2(s)\,ds+\int^\infty_0\int^t_{t-s}z^2(\theta)\,d\theta|F(s,x^*)|\,dK(s). \end{align*} $$

Then,

(3.5) $$ \begin{align} LV_2(z_t)=\rho(z^2(t)-z^2(t-\tau))+\alpha_0z^2(t)-\int^\infty_0z^2(t-s)|F(s,x^*)|\,dK(s). \end{align} $$

From (3.3), (3.4), (3.5) and (3.2) for the Lyapunov functional

$$ \begin{align*}V(z_t)=V_1(z(t))+V_2(z_t)\end{align*} $$

we obtain

$$ \begin{align*} LV(z_t)&\le(-2a+\sigma^2)z^2(t)-2bz(t)z(t-\tau)+2\alpha_0z^2(t)+\rho(z^2(t)-z^2(t-\tau)) \nonumber\\ &=(-2a+2\alpha_0+\rho+\sigma^2)z^2(t)-2bz(t)z(t-\tau)-\rho z^2(t-\tau) \nonumber\\ &=\eta'D\eta, \end{align*} $$

where

$$ \begin{align*} D = \left[\begin{matrix} -\rho & -b \\ - b & -\rho \end{matrix}\right], \quad \eta=\left[\begin{matrix} z(t) \\ z(t-\tau) \end{matrix}\right]. \end{align*} $$

Via condition (3.2), the matrix D is negative definite, that is, $LV(z_t)\le -cz^2(t)$ for some $c>0$ . From the Lyapunov-type theorem [16], it follows that the zero solution of the equation (2.4) is asymptotically mean square stable. This completes the proof.

Remark 3.5. If $dK(s)=0$ , then, from (3.2), the known [16] sufficient condition of asymptotic mean square stability follows for the zero solution of the linear delay differential equation

(3.6) $$ \begin{align} a>|b|+\tfrac{1}{2}\sigma^2. \end{align} $$

For equation (2.5), the stability condition (3.6) takes the form

$$ \begin{align*} a>a\bigg|\ln\frac{k}{a}-1 \bigg|+\frac{1}{2}\sigma^2, \end{align*} $$

from which it also follows that

$$ \begin{align*} a\exp\bigg(\frac{\sigma^2}{2a}\bigg)<k<a\exp\bigg(2-\frac{\sigma^2}{2a}\bigg), \quad \sigma^2<2a. \end{align*} $$

Remark 3.6. The necessary and sufficient condition for asymptotic mean square stability of the zero solution of the equation (2.5) is [15, 16]

(3.7) $$ \begin{align} p\;G<1, \end{align} $$

where

(3.8) $$ \begin{align} p=\dfrac{1}{2}\sigma^2,\quad \!G=\left\{\begin{matrix} \dfrac{1+\dfrac{a}{q}\bigg(\ln\dfrac{k}{a}-1\bigg)\sin(qh)}{a[1+\bigg(\ln\dfrac{k}{a}-1\bigg)\cos(qh)]}, & k>a{e}^2, & q=a\sqrt{\ln\dfrac{k}{a}\bigg(\ln\dfrac{k}{a}-2\bigg)}, \\ \\ \dfrac{1+ah}{2a}, & k=a{e}^2, \\ \\ \dfrac{1+\dfrac{a}{q}\bigg(\ln\dfrac{k}{a}-1\bigg)\sinh(qh)}{a[1+\bigg(\ln\dfrac{k}{a}-1\bigg)\cosh(qh)]}, & a<k<a{e}^2, & q=a\sqrt {\ln\dfrac{k}{a}\bigg(2-\ln\dfrac{k}{a}\bigg)}. \end{matrix} \right. \end{align} $$

In particular, if $p>0$ , $h=0$ , then the stability condition (3.7), (3.8) takes the form $a\ln (k/a)>p$ ; if $p=0$ , $h>0$ , then the region of stability is bounded by the lines $a=0$ , $a=k$ and

$$ \begin{align*}1+\bigg(\ln\frac{k}{a}-1\bigg)\cos(qh)=0 \quad \text{for }k>a{e}^2.\end{align*} $$

Note that

$$ \begin{align*}\sinh (x)=\tfrac{1}{2}(e^x-e^{-x})\quad \mathrm{and} \quad \cosh (x)=\tfrac{1}{2}(e^x+e^{-x})\end{align*} $$

are hyperbolic sine and hyperbolic cosine functions, respectively.

3.2 Delay-dependent stability condition

Note that the linear equation (2.4) can be presented in the form of a stochastic differential equation of a neutral type [16],

(3.9) $$ \begin{align} dZ(t)=-\gamma z(t)\,dt+\sigma z(t)\,dw(t), \end{align} $$

where

(3.10) $$ \begin{align}\begin{aligned} Z(t)&=z(t)-bI_1(t)+I_2(t),\quad \gamma=a+b-\beta,\quad \beta=\int^\infty_0F(s,x^*)\,dK(s), \\ I_1(t)&=\int^t_{t-\tau}z(s)\,ds,\quad I_2(t)=\int^\infty_0\int^t_{t-s}z(\theta)\,d\theta F(s,x^*)\,dK(s).\end{aligned} \end{align} $$

Really, from (3.10) and (3.9),

$$ \begin{align*} dZ(t)&=dz(t)+\bigg(-bz(t)+bz(t-\tau)+\beta z(t)-\int^\infty_0z(t-s)F(s,x^*)\,dK(s)\bigg)\,dt , \nonumber\\ &=(-az(t)-bz(t)+\beta z(t))\, dt+\sigma z(t)\,dw(t), \end{align*} $$

which coincides with (2.4).

Theorem 3.7. Let $2\gamma>\sigma ^2$ and $\alpha _1$ be as defined in (3.1). If

(3.11) $$ \begin{align} |b|\tau+\alpha_1<1 \end{align} $$

and there exist positive numbers $p_1$ and $p_2$ such that the linear matrix inequality (LMI)

(3.12) $$ \begin{align} Q=\left[\begin{matrix} -(2\gamma-p_1\tau^2-p_2\alpha_1^2-\sigma^2) & b\gamma & -\gamma\\[.2cm] b\gamma & -p_1 & 0 \\[.2cm] -\gamma & 0 & -p_2 \end{matrix}\right]<0 \end{align} $$

holds, then the positive equilibrium $x^*$ of equation (1.6) is stable in probability.

Proof. Via Remark 3.3, it is enough to prove that, by conditions (3.11) and (3.12), the zero solution of equation (2.4) is asymptotically mean square stable. Let L be the generator [6, 16] of equation (3.9). Via (3.9) and (3.10), for the functional $V_1(z_t)=Z^2(t)$ ,

$$ \begin{align*} LV_1(z_t)&=-2\gamma Z(t)z(t)+\sigma^2z^2(t)\\ &=-2\gamma(z(t)-bI_1(t)+I_2(t))z(t)+\sigma^2z^2(t)\\ &=-(2\gamma-\sigma^2)z^2(t)+2b\gamma z(t)I_1(t)-2\gamma z(t)I_2(t). \end{align*} $$

Using the general method for construction of Lyapunov functionals [15, 16], we now consider an additional functional $V_2$ in the form

$$ \begin{align*} V_2(z_t)=p_1\tau\int^t_{t-\tau}(s-t+\tau)z^2(s)\,ds+p_2\alpha_1\int^\infty_0\int^t_{t-s}(\theta-t+s)z^2(\theta)\, d\theta|F(s,x^*)|\,dK(s). \end{align*} $$

Then,

$$ \begin{align*} LV_2(z_t)&=p_1\tau\bigg(\tau z^2(t)-\int^t_{t-\tau}z^2(s)\,ds\bigg) +p_2\alpha_1\int^\infty_0\bigg(sz^2(t)-\int^t_{t-s}z^2(\theta)\,d\theta\bigg)|F(s,x^*)|\,dK(s)\\ &=(p_1\tau^2+p_2\alpha^2_1)z^2(t)-p_1\tau\int^t_{t-\tau}z^2(s)\,ds -p_2\alpha_1\int^\infty_0\!\int^t_{t-s}z^2(\theta)\,d\theta|F(s,x^*)|\,dK(s). \end{align*} $$

Using Lemma A.1 (see Appendix A),

$$ \begin{align*} I_1^2(t)\le\tau\int^t_{t-\tau}z^2(s)\,ds \end{align*} $$

and

$$ \begin{align*} I_2^2(t)&\le \bigg(\int^\infty_0\int^t_{t-s}|z(\theta)|\,d\theta|F(s,x^*)|\,dK(s)\bigg)^2\\ &\le\alpha_1\int^\infty_0\int^t_{t-s}z^2(\theta)\,d\theta|F(s,x^*)|\,dK(s). \end{align*} $$

So,

$$ \begin{align*} LV_2(z_t)\le(p_1\tau^2+p_2\alpha_1^2)z^2(t)-p_1I_1^2(t)-p_2I_2^2(t). \end{align*} $$

As a result, for the functional $V(z_t)=V_1(z_t)+V_2(z_t)$ , we obtain

$$ \begin{align*} LV(z_t)&\le-(2\gamma-p_1\tau^2-p_2\alpha_1^2-\sigma^2)z^2(t)+2b\gamma z(t)I_1(z_t)-2\gamma z(t)I_2(t)\\ &\quad -p_1I_1^2(t)-p_2I_2^2(t)\\ &=\eta'(t)Q\eta(t), \end{align*} $$

where $\eta (t)=( z(t), I_1(t), I_2(t))'$ and matrix Q is as defined in (3.12). From the LMI $Q<0$ , it follows that $LV(z_t)\le -cz^2(t)$ for some $c>0$ . From this and (3.11), it follows that the zero solution of the neutral-type equation is asymptotically mean square stable [16]. The proof is now complete.

Remark 3.8. For equation (1.1), we have $b=0$ and $dK(s)=k\delta (s-h)ds$ , where $\delta (s)$ is Dirac’s function [16]. Then condition (3.11) takes the form $\alpha _1=ah|\ln ({k}/{a})-1|<1$ and can be presented as

$$ \begin{align*} a\exp(1-(ah)^{-1})<k<a\exp(1+(ah)^{-1}). \end{align*} $$

Example 3.9. Via Remark 1.1, for $m=2$ , equation (1.6) takes the form

(3.13) $$ \begin{align} \hspace{-16pt}dx(t)&=(-ax(t)-bx(t-\tau)+k_1x(t-h_1)e^{-\lambda_1x(t-h_1)}+k_2x(t-h_2)e^{-\lambda_2x(t-h_2)})\,dt \nonumber \\ \hspace{-16pt}&\quad +\sigma(x(t)-x^*)\,dw(t). \end{align} $$

Let $m=2$ , $a=b=0.5$ , $k_1=\lambda _1=1$ , $k_2=\lambda _2=2$ , $\tau =0.8$ , $h_1=0.4$ , $h_2=0.6$ and ${\sigma =1}$ . Note that, for the given values of the parameters, condition (3.2) does not hold. We now check conditions (3.11) and (3.12). Via Example 1.2, the equilibrium $x^*=\ln 2=0.6931$ . From (3.1), it follows that

$$ \begin{align*} \alpha_1=\tfrac{1}{2}\{(1-\ln2)h_1+(2\ln2-1)h_2\}=0.1773. \end{align*} $$

So, condition (3.11) holds, that is,

$$ \begin{align*}|b|\tau+\alpha_1=0.4+0.1773=0.5773<1.\end{align*} $$

In addition, via (3.10),

$$ \begin{align*} \beta=\tfrac{1}{2}(1-\ln2)+1-2\ln2=1.5-2.5\ln2=-0.2329,\quad \gamma=1.2329. \end{align*} $$

By virtue of MATLAB for the LMI (3.12), we found $p_1=0.7256$ and $p_2=4.7962$ , for which the matrix

$$ \begin{align*} Q=\left[\begin{matrix} -0.8506 & 0.6164 & -1.2329\\[.2cm] 0.6164 & -0.7256 & 0 \\[.2cm] -1.2329 & 0 & -4.7963 \end{matrix}\right] , \end{align*} $$

is negative definite, that is,

$$ \begin{align*}S_1=-6.3725<0, \quad S_2=6.2775>0,\quad S_1S_2=-40.0037<S_3=-0.0350<0\end{align*} $$

(see Lemma A.3, Appendix A). So, the equilibrium $x^*=\ln 2$ of the equation (3.13) is stable in probability.

Using the Euler–Maruyama scheme [11, 13, 16] for a difference analogue of equation (3.13), the following simulations were obtained.

In Figure 1, 100 trajectories of the solution to equation (3.13) are shown for the values of the parameters given above and the initial function $x(s)=1.18\cos (s)$ , ${s\in [-\tau ,0]}$ . One can see that all trajectories converge to the stable equilibrium ${x^*=\ln 2}$ .

Figure 1 Stable equilibrium. 100 trajectories of the solution to equation (3.13): $a=0.5$ , $b=0.5$ , $x(s)=1.18\cos (s)$ , $s\in [-\tau ,0]$ .

Now, put $b=0.3$ . From (1.5), it follows that $x^*=0.43$ ; via MATLAB, it is shown that, in this case, the LMI (3.12) is impossible.

In Figure 2, 100 trajectories of the solution to equation (3.13) are shown for $b=0.3$ and the same values of all other parameters. The equilibrium $x^*$ is unstable; therefore, trajectories do not converge to $x^*$ and fill the whole space.

Figure 2 Unstable equilibrium. 100 trajectories of the solution to equation (3.13): $a=0.5$ , $b=0.3$ , $x(s)=x^*\cos (s)$ , $s\in [-\tau ,0]$ .

4 Conclusion

Stability is one of the important concepts in the study of mathematical models. We consider a very general kind of stochastic model that is motivated by Nicholson’s model. We obtain stability conditions, using the general method for construction of Lyapunov functionals and the LMI method. The obtained stability conditions can be easily verified, and are less restrictive. Via numerical simulation of the solutions to the considered equation, some graphs are obtained to demonstrate stable and unstable equilibria. The research method used here can be extended to many more general and complicated models in different applications.