1 Introduction
The Nicholson blowflies model,
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
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
Remark 1.1. Let
In this case, $K=\sum \nolimits ^m_{i=1}k_i$ and the equation (1.3) takes the form
If, in particular, $m=1$ , $\lambda _1=\lambda>0$ , then $K=k_1$ and
If $m=2$ , then the equilibrium $x^*$ is defined by the equation
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],
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
where $y_t$ is the trajectory of $y(s)$ , $s\le t$ , and
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
and (1.3), we transform $I(y_t)$ as follows.
where
Substituting (2.2) into (2.1) and neglecting $o(y)$ , we obtain the linear equation
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),
and equation (2.4) takes the form
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,
3.1 Delay-independent stability condition
Theorem 3.4. If
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,
where $F(s,x^*)$ is as defined in (2.3) and
Now, consider the additional functional $V_2(z_t)$ in the form
Then,
From (3.3), (3.4), (3.5) and (3.2) for the Lyapunov functional
we obtain
where
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
For equation (2.5), the stability condition (3.6) takes the form
from which it also follows that
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]
where
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
Note that
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],
where
Really, from (3.10) and (3.9),
which coincides with (2.4).
Theorem 3.7. Let $2\gamma>\sigma ^2$ and $\alpha _1$ be as defined in (3.1). If
and there exist positive numbers $p_1$ and $p_2$ such that the linear matrix inequality (LMI)
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)$ ,
Using the general method for construction of Lyapunov functionals [Reference Shaikhet15, Reference Shaikhet16], we now consider an additional functional $V_2$ in the form
Then,
Using Lemma A.1 (see Appendix A),
and
So,
As a result, for the functional $V(z_t)=V_1(z_t)+V_2(z_t)$ , we obtain
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
Example 3.9. Via Remark 1.1, for $m=2$ , equation (1.6) takes the form
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
So, condition (3.11) holds, that is,
In addition, via (3.10),
By virtue of MATLAB for the LMI (3.12), we found $p_1=0.7256$ and $p_2=4.7962$ , for which the matrix
is negative definite, that is,
(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}$ .
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.
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,
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
Here, in particular,
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.