1 Introduction
The Nicholson blowflies model,
$$ \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 [Reference Nicholson12], 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 [Reference Berezansky, Braverman and Idels2, Reference Bradul and Shaikhet3, Reference Ding and Alzabut5, Reference Gurney, Blythe and Nisbet7–Reference Liu10, Reference Shaikhet15–Reference So and Yu17] 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
$$ \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 [Reference Rudin14]. 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 [Reference Shaikhet16], 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
$$ \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
$$ \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
$$ \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 [Reference Gikhman and Skorokhod6],
$$ \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 [Reference Gikhman and Skorokhod6, Reference Shaikhet16].
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
$$ \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.
$$ \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
$$ \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
$$ \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
$$ \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 [Reference Shaikhet16] 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,
$$ \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
$$ \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 [Reference Shaikhet15, Reference Shaikhet16], 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 [Reference Gikhman and Skorokhod6, Reference Shaikhet16] of equation (2.4). Then,
$$ \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
$$ \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,
$$ \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 [Reference Shaikhet16], 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 [Reference Shaikhet16] sufficient condition of asymptotic mean square stability follows for the zero solution of the linear delay differential equation
$$ \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 [Reference Shaikhet15, Reference Shaikhet16]
$$ \begin{align} p\;G<1, \end{align} $$
where
$$ \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 [Reference Shaikhet16],
$$ \begin{align} dZ(t)=-\gamma z(t)\,dt+\sigma z(t)\,dw(t), \end{align} $$
where
$$ \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
$$ \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)
$$ \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 [Reference Gikhman and Skorokhod6, Reference Shaikhet16] 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 [Reference Shaikhet15, Reference Shaikhet16], 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 [Reference Shaikhet16]. 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 [Reference Shaikhet16]. 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
$$ \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 [Reference Maruyama11, Reference Pang, Deng and Mao13, Reference Shaikhet16] 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.
Appendix A
Lemma A.1 [Reference Burgos, Cortes, Shaikhet and Villanueva4].
Let
$R\in \mathbf R^{n\times n}$
be a positive definite matrix,
$y=\int _Dx(s)\mu \, (ds)$
, where
$y, x(s)\in \mathbf R^n$
, and
$\mu (ds)$
is some measure on D such that
$\mu (D)<\infty $
. Then,
$$ \begin{align*} y'Ry\le\mu(D)\int_Dx'(s)Rx(s)\mu\,(ds). \end{align*} $$
Definition A.2 [Reference Abbas and Shaikhet1, Reference Shaikhet16].
The trace of the kth order of an
$n\times n$
-matrix
$A=\|a_{ij}\|$
is defined as
$$ \begin{align*}\begin{array}{ll} S_k=\displaystyle\sum\limits_{1\le i_1<\cdots<i_k\le n}\left|\begin{array}{lll} \; a_{i_1i_1} &\cdots& a_{i_1i_k} \\ \cdots & \cdots & \cdots \\ \; a_{i_ki_1} &\cdots& a_{i_ki_k}\end{array}\right|,\quad k=1,\ldots,n. \end{array}\end{align*} $$
Here, in particular,
$$ \begin{align*}S_1=Tr(A),\quad S_n=\det(A), \quad S_{n-1}=\sum\limits^n_{i=1}A_{ii},\end{align*} $$
where
$A_{ii}$
is the algebraic complement of the diagonal element
$a_{ii}$
of the matrix A.
Lemma A.3 [Reference Abbas and Shaikhet1, Reference Shaikhet16].
A
$3\times 3$
-matrix A is the Hurwitz matrix [Reference Rudin14] if and only if
$S_1<0$
,
$S_1S_2<S_3<0$
.
Acknowledgement
The authors are very grateful to an anonymous reviewer for their useful remarks.
