The Stochastic Asymptotic Stability Analysis in Two Species Lotka-Volterra Model

Yuqin Li¹, Yuehua He²

¹Institute of Mathematics, Zunyi Vocational and Technical College, Zunyi, China.
²School of Science, Xuchang University, Xuchang, China.

DOI: 10.4236/am.2023.147028                

Abstract

The asymptotic stability of two species stochastic Lotka-Volterra model is explored in this paper. Firstly, the Lotka-Volterra model with random parameter is built and reduced into the equivalent deterministic system by orthogonal polynomial approximation. Then, the linear stability theory and Routh-Hurwitz criterion for nonlinear deterministic systems are applied to the equivalent one. At last, at the aid of Lyapunov second method, we obtain that as the random intensity or statistical parameter of random variable is changed, the stability about stochastic Lotka-Volterra model is different from the deterministic system.

Keywords

Asymptotic Stability, Stochastic Lotka-Volterra Model, Lyapunov Method

1. Introduction

Stability is a hot topic in the ecological literature. It is not only associated with the structure, function and evolution of ecosystem, but also closely related to strength and characteristics of external disturbance. It is well known that there exist various kinds of stabilities from the different aspects of research both internal and external, such as resistance, recovery, persistence, variability and so on. For example, the recovery refers to the ability to return to the original state after suffering the external disturbance. Although the stability of ecosystem has been studied comprehensively and systematically in theory, research on the stability of ecological system is and will be an eternal topic due to the ecosystem is constantly changing. Hence, study on stability has been vigorously done from mathematical and applied perspectives whether present or in the future.

Lotka-Volterra system as a model of undamped oscillations in autocatalytic chemical reactions, was firstly developed by Alfred J. Lotka independently in 1920 [1] and was later applied by Vito Volterra in 1926 [2] to treat predator-prey interactions in ecology. From then on, the application of Lotka-Volterra system has a wide range in various fields, such as population dynamics, epidemiology, physics, economics and chemistry [3] [4] [5] . Other applications occur in neural networks, game theory, plasma physics and so on [6] [7] . Owing to its unrealistic stability, Lotka-Volterra model, as a starting point is a more advanced model in the analysis of population dynamics, and receives much attention for a long time including the addition of small random perturbations, polynomial interactions, time delayed [8] [9] and diffusion effected [10] . Liu and Luo [11] gave the discriminant conditions of local stability, global stability and proved that the local stability and global stability are equivalent at the equilibrium position. Reference [12] has obtained the sufficient conditions of stability in a nonautonomous Lotka-Volterra system based on the differential and integral method. The condition for the existence of a globally stable equilibrium of n-dimensional Lotka-Volterra systems is presented in references [13] and [14] . In reference [11] , based on the large deviation principle, the influence of some small random perturbations on stability and extinction to predator-prey systems is explored. Up to now, some researchers have already studied some other dynamical behavior such as bifurcation, chaos and so on [15] [16] [17] . To our best knowledge, as a mathematical model, the system with random physical parameters is more close to the actual, so the research about stability under the influence of random internal parameters in these systems is tremendously done form a practical point of view.

Motivated by the above discussion, the statistical characteristic of random variable, a Lotka-Volterra system with random parameter is investigated by the orthogonal polynomial approximation in this paper.

This paper is organized as follows. Transformation of the stochastic Lotka-Volterra system into its equivalent deterministic one by orthogonal polynomial approximation is shown in Section 2. Section 3 investigates the locally asymptotic stability and globally asymptotic stability of stochastic Lotka-Volterra system. And in Section 4, the locally and globally asymptotic stability of one of the prey-predator models is discussed. Finally, conclusions are drawn in Section 5.

2. Orthogonal Polynomial Approximation for Stochastic Lotka-Volterra Model

The classic two species Lotka-Volterra mathematical model can be described as

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1=x_1\left({\bar{b}}_1+a_{11}{x}_1+a_{12}{x}_2\right)\\ {\stackrel{\dot{}}{x}}_2=x_2({\bar{b}}_2+a_{21}{x}_1+a_{22}{x}_2)\end{array}$ $\left({\stackrel{\dot{}}{x}}_i\triangleq \text{d}{x}_i/\text{d}t,i=1,2\right)$ (1)

where $x_1,x_2$ are all nonnegative variables and ${\bar{b}}_1,{\bar{b}}_2$ account for the self-regulation of each species. Biologically, we can interpret this system as follows, ${\bar{b}}_2<0$ stands for the natural mortality of heterotrophic organisms such as predators (or parasites). If ${\bar{b}}_1<0$ , ${\bar{b}}_2<0$ that means the mortality of two symbiotic species separated. The case ${\bar{b}}_1>0$ , ${\bar{b}}_2<0$ means that the species $x_1$ interacts strongly with species $x_2$ and weakly among themselves. Similarly, ${\bar{b}}_1<0$ , ${\bar{b}}_2>0$ , species $x_2$ interact more strongly with $x_1$ than they do with themselves.

The parameters $a_{ij}\left(i=1,2;j=1,2\right)$ describe how species $x_1$ and $x_2$ interact each other as follows:

1) $a_{12}<0,a_{21}>0$ , represent prey-predator models;

2) $a_{12}<0,a_{21}<0$ , represent competition model;

3) $a_{12}>0,a_{21}>0$ , represent cooperation model.

What’s more, $a_{11}<0,a_{22}<0$ refer to that the two species are controlled by the density and $a_{12}=0,a_{21}>0$ (or $a_{12}>0,a_{21}=0$ ) means dominance effect; $a_{12}=0,a_{21}<0$ (or $a_{12}<0,a_{21}=0$ ) means partial damage effect; $a_{12}=0,a_{21}=0$ represents neutral relationship. With these interpretations, only solutions of (1) with $x_1$ and $x_2$ nonnegative are physically interest. The equilibrium points of

(1) are $A_1\left(0,0\right)$ , $A_2\left(0,-\frac{{\bar{b}}_2}{a_{22}}\right)$ , $A_3\left(-\frac{{\bar{b}}_1}{a_{11}},0\right)$ , $A_4\left(\frac{a_{12}{\bar{b}}_2-a_{21}{\bar{b}}_1}{a_{21}{a}_{12}-a_{22}{a}_{11}},\frac{a_{11}{\bar{b}}_2-a_{21}{\bar{b}}_1}{a_{21}{a}_{12}-a_{22}{a}_{11}}\right)$ .

Obviously $A_1,A_2,A_3$ are unstable, and the stability of $A_4$ is all depends on the parameter. This article focuses on the nontrivial equilibrium position $A_4$ of system (1) which is the solution of the equation: $b_i+{\displaystyle \underset{j=1}{\overset{2}{\sum }}{a}_{ij}{x}_i^{*}}=0,i=1,2$ , where $x_1^{*}=\frac{a_{12}{\bar{b}}_2-a_{21}{\bar{b}}_1}{a_{21}{a}_{12}-a_{22}{a}_{11}}$ , $x_2^{*}=\frac{a_{11}{\bar{b}}_2-a_{21}{\bar{b}}_1}{a_{21}{a}_{12}-a_{22}{a}_{11}}$ .

Utilizing the coordinate transformation

$\{\begin{array}{l}{x}_1=x+x_1^{*},\\ x_2=y+x_2^{*}.\end{array}$ (2)

The nontrivial equilibrium point $A_4$ is converted to origin (0, 0). Then we can obtain the following Lotka-Volterra model

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{x}^2+a_{12}xy+{\bar{b}}_{1}x+2a_{11}{x}_1^{*}x+a_{12}{x}_2^{*}x+a_{12}{x}_2^{*}y+{\bar{b}}_1{x}_1^{*}+a_{11}{x}_1^{*}{}^2+a_{12}{x}_1^{*}{x}_2^{*},\\ \frac{\text{d}y}{\text{d}t}=a_{22}{y}^2+a_{21}xy+{\bar{b}}_{2}y+2a_{22}{x}_2^{*}y+a_{21}{x}_1^{*}y+a_{21}{x}_2^{*}x+{\bar{b}}_2{x}_2^{*}+a_{22}{x}_2^{*}{}^2+a_{21}{x}_1^{*}{x}_2^{*}.\end{array}$ (3)

The necessary and sufficient condition of locally asymptotic stability of deterministic dynamics Equation (3) at the zero solution is that $x_1^{*}{a}_{11}+x_2^{*}{a}_{22}<0$ , $a_{21}{a}_{12}-a_{22}{a}_{11}<0$ [1] . The sufficient condition of locally asymptotic stability of this deterministic dynamics equation at the zero solution is that the zero solution is locally asymptotic stable and two species are controlled by the density. Namely, $\{\begin{array}{l}{x}_1^{*}{a}_{11}+x_2^{*}{a}_{22}<0,a_{21}{a}_{12}-a_{22}{a}_{11}<0\\ a_{11}<0\\ a_{22}<0\end{array}$ .

For natural disaster or other reasons for each species, the parameters $b_i,\left(i=1,2\right)$ should be affected by uncertain factors. So suppose that the parameter ${\bar{b}}_i,\left(i=1,2\right)$ can be described as a random parameter, ${\bar{b}}_i=b_i+\delta u,\left(i=1,2\right)$ where $b_i,\left(i=1,2\right)$ is deterministic parameters of ${\bar{b}}_i,\left(i=1,2\right)$ , $\delta$ is regarded as intensity of random disturbance and u is random parameter which defined on [−1, 1] subject to some probability density function.

Corresponding to this random variable, the orthogonal ploynomial is choosen as,

$G_n^{\eta }\left(u\right)={\displaystyle \underset{k=0}{\overset{n}{\sum }}\frac{1}{k!\left(n-k\right)!}}\frac{\left(2\eta \right){\left(2\eta +n\right)}_k}{{\left(\eta +\frac{1}{2}\right)}_k}{\left(\frac{u-1}{2}\right)}^k$ , (4)

where

${\eta }_k=\eta \left(\eta +1\right)\cdots \left(\eta +k-1\right)=\frac{\Gamma \left(\eta +k\right)}{\Gamma \left(\eta \right)},{\eta }_0\equiv 1$

And the Orthogonal of Gegenbauer polynomial can be expressed as

${\displaystyle {\int }_{-1}^{1}p\left(u\right)G_m\left(u\right)G_n\left(u\right)\text{d}u}=\{\begin{array}{l}0,m\ne n\\ b_n^{\eta },m=n\end{array}$ , (5)

where $b_n^{\eta }=\frac{1}{2^{2\eta -1}n!}\cdot \frac{\pi \eta \Gamma \left(2\eta +n\right)}{\left(\eta +n\right)\Gamma \left(\eta \right)\Gamma \left(\eta +\frac{1}{2}\right)\Gamma (\; 1\; 2\; )}$

Substituting Equation (4) and Equation (5) into Equation (3), we get

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{x}^2+a_{12}xy+{\bar{b}}_{1}x+2a_{11}{x}_1^{*}x+a_{12}{x}_2^{*}x+a_{12}{x}_2^{*}y+{\bar{b}}_1{x}_1^{*}+a_{11}{x}_1^{*}{}^2+a_{12}{x}_1^{*}{x}_2^{*}+\delta ux,\\ \frac{\text{d}y}{\text{d}t}=a_{22}{y}^2+a_{21}xy+{\bar{b}}_{2}y+2a_{22}{x}_2^{*}y+a_{21}{x}_1^{*}y+a_{21}{x}_2^{*}x+{\bar{b}}_2{x}_2^{*}+a_{22}{x}_2^{*}{}^2+a_{21}{x}_1^{*}{x}_2^{*}+\delta uy.\end{array}$ (6)

By the following replacement,

$A=2a_{11}{x}_1^{*}+a_{12}{x}_2^{*}+b_1,B=a_{12}{x}_2^{*},C=a_{21}{x}_1^{*},D=a_{21}{x}_1^{*}+2a_{22}{x}_2^{*}+b_2$ ,

$P={\bar{b}}_1{x}_1^{*}+a_{11}{x}_1^{*}{}^2+a_{12}{x}_1^{*}{x}_2^{*},Q={\bar{b}}_2{x}_2^{*}+a_{22}{x}_2^{*}{}^2+a_{21}{x}_1^{*}{x}_2^{*}$ ,

we simplify the Lotka-Volterra model as

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{x}^2+a_{12}xy+Ax+By+\delta ux+P,\\ \frac{\text{d}y}{\text{d}t}=a_{22}{y}^2+a_{21}xy+Dy+Cx+\delta uy+Q.\end{array}$ (7)

The only one equilibrium point of Equation (7) is (0, 0). Furthermore, the response of Lotka-Volterra system with random parameter can be approximately expressed by the following Fourier series under the condition of the convergence in mean square

$\{\begin{array}{l}x\left(t,u\right)={\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i}\left(u\right),\\ y\left(t,u\right)={\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i}\left(u\right).\end{array}$ (8)

where $G_i\left(u\right)$ is the ith Gegenbauer orthogonal polynomial, M represents the largest order of the polynomial we have taken.

Substituting Equation (8) into Equation (7), we obtain

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}\right)}^2+a_{12}\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}\right)\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}\right)\\ +A{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}+B{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}+\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)u{G}_i\left(u\right)}+P,\\ \frac{\text{d}y}{\text{d}t}=a_{22}{\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}\right)}^2+a_{21}\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}\right)\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}\right)\\ +D{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}+C{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}+\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)u{G}_i\left(u\right)}+Q.\end{array}$ (9)

The coefficients of $G_i^{\eta }\left(\mu \right)G_j^{\eta }\left(\mu \right)$ can be donated as

$\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i^{\eta }\left(u\right)}\right)\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i^{\eta }\left(u\right)}\right)={\displaystyle \underset{i=0}{\overset{2M}{\sum }}{H}_i\left(t\right)G_i^{\eta }}\left(u\right),$

${\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i^{\eta }\left(u\right)}\right)}^2={\displaystyle \underset{i=0}{\overset{2M}{\sum }}{L}_i\left(t\right)G_i^{\eta }}\left(u\right),$ (10)

${\left({\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i^{\eta }\left(u\right)}\right)}^2={\displaystyle \underset{i=0}{\overset{2M}{\sum }}{M}_i\left(t\right)G_i^{\eta }}\left(u\right).$

where $H_i\left(t\right),L_i\left(t\right),M_i\left(t\right)\left(i=0,1,2,\cdots ,2\times M\right)$ which stands for the linear combination of non-linearity can be calculated by Maple.

Substituting Equation (10) into Equation (9), we obtain

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{M}_i\left(t\right)G_i\left(u\right)}+a_{12}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{H}_i\left(t\right)G_i\left(u\right)}+A{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}\\ +B{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}+\delta u{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}+P,\\ \frac{\text{d}y}{\text{d}t}=a_{22}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{L}_i\left(t\right)G_i\left(u\right)}+a_{21}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{H}_i\left(t\right)G_i\left(u\right)}+D{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}\\ +C{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}+\delta u{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}+Q.\end{array}$ (11)

With the help of the cycle recurrence formula of Gegenbauer polynomial

$u{G}_n^{\eta }\left(u\right)={\alpha }_n{G}_{n-1}^{\eta }\left(u\right)+{\gamma }_n{G}_{n+1}^{\eta }\left(u\right).$ (12)

where ${\alpha }_n=\frac{2\eta +n-1}{2\left(\eta +n\right)}$ , ${\gamma }_n=\frac{n+1}{2\left(\eta +n\right)}$ .

The stochastic term and the non-linearity in the right equation of system (11) can be written as $\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)u{G}_i^{\eta }\left(u\right)}=\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}\left({\alpha }_i{x}_{i-1}\left(t\right)G_{i-1}^{\eta }\left(u\right)+{\gamma }_i{x}_{i+1}\left(t\right)G_{i+1}^{\eta }\left(u\right)\right)},$ (13)

$\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)u{G}_i^{\eta }\left(u\right)}=\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}\left({\alpha }_i{y}_{i-1}\left(t\right)G_{i-1}^{\eta }\left(u\right)+{\gamma }_i{y}_{i+1}\left(t\right)G_{i+1}^{\eta }\left(u\right)\right)}.$ (14)

Substituting Equation (13) and Equation (14) into Equation (11), we get

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{11}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{M}_i\left(t\right)G_i\left(u\right)}+a_{12}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{H}_i\left(t\right)G_i\left(u\right)}+A{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}\\ +B{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}+\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}\left({\alpha }_i{x}_{i-1}\left(t\right)G_{i-1}^{\eta }\left(u\right)+{\gamma }_i{x}_{i+1}\left(t\right)G_{i+1}^{\eta }\left(u\right)\right)}+P,\\ \frac{\text{d}y}{\text{d}t}=a_{22}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{L}_i\left(t\right)G_i\left(u\right)}+a_{21}{\displaystyle \underset{i=0}{\overset{2M}{\sum }}{H}_i\left(t\right)G_i\left(u\right)}+D{\displaystyle \underset{i=0}{\overset{M}{\sum }}{y}_i\left(t\right)G_i\left(u\right)}\\ +C{\displaystyle \underset{i=0}{\overset{M}{\sum }}{x}_i\left(t\right)G_i\left(u\right)}+\delta {\displaystyle \underset{i=0}{\overset{M}{\sum }}\left({\alpha }_i{y}_{i-1}\left(t\right)G_{i-1}^{\eta }\left(u\right)+{\gamma }_i{y}_{i+1}\left(t\right)G_{i+1}^{\eta }\left(u\right)\right)}+Q.\end{array}$ (15)

Multiply both sides of Equation (14) by $G_i^{\eta }\left(u\right),\left(i=1,2,\cdots ,M\right)$ in sequence and take expectation. Based on the orthogonal of the polynomial approximation of convergence random function in the Hilbert spaces and the orthogonality of Gegenbauer orthogonal polynomials, we can finally obtain the equivalent deterministic equation. As $M\to \infty$ , the Lotka-Volterra system with random parameter is strictly equivalent to the system (7) under condition of the convergence in mean square. We denote the coefficient of $H_i,M_i$ , and $L_i\left(i=0,1,\cdots ,2M\right)$ in the linear combination as $K_i,S_i\left(i=0,1,\cdots ,M\right)$ respectively. According to the principle of approximation $x_{-1},y_{-1}$ are zero. The nonlinear term of

Equation (11) can be expanded into

$\{\begin{array}{l}\{\begin{array}{l}\frac{\text{d}{x}_0}{\text{d}t}=A{x}_0+B{y}_0+{\gamma }_0\delta x_1+K_0,\\ \frac{\text{d}{y}_0}{\text{d}t}=D{y}_0+C{x}_0+{\gamma }_0\delta y_1+S_0.\end{array}\\ \{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=A{x}_1+B{y}_1+\delta \left({\alpha }_1{x}_0+{\gamma }_1{x}_2\right)+K_1,\\ \frac{\text{d}{y}_1}{\text{d}t}=D{y}_1+C{x}_1+\delta \left({\alpha }_1{y}_0+{\gamma }_1{y}_2\right)+S_1.\end{array}\\ \cdots ,\\ \{\begin{array}{l}\frac{\text{d}{x}_M}{\text{d}t}=A{x}_M+B{y}_M+\delta \left({\alpha }_M{x}_{M-1}+{\gamma }_M{x}_{M+1}\right)+K_M,\\ \frac{\text{d}{y}_M}{\text{d}t}=D{y}_M+C{x}_M+\delta \left({\alpha }_M{y}_{M-1}+{\gamma }_M{y}_{M+1}\right)+S_M.\end{array}\end{array}$ (16)

3. The Stability Analysis of Zero Solution

In this section we take zero solution as the example and discuss the asymptotic stability in stochastic Lotka-Volterra system (7). The equivalent deterministic system (16) can be rewritten as

$\stackrel{\dot{}}{Z}\left(t\right)=JZ\left(t\right)+f\left(x_i\left(t\right),y_i\left(t\right)\right),$

where $Z\left(t\right)={\left\{x_0\left(t\right),y_0\left(n\right),x_1\left(t\right),y_1\left(n\right),\cdots ,x_M\left(t\right),y_M\left(n\right)\right\}}^{\text{T}}$ , J is the coefficient matrix and $f\left(x_i\left(t\right),y_i\left(t\right)\right)=\circ \left(Z\left(t\right)\right)$ .

In order to analyze the stability of Equation (16) at the equilibrium point. We first need two lemmas.

Lemma 1. (Routh-Hurwitz criterion) All of the eigenvalues of the Equation (16) have real parts strictly less than zero if and only if all elements in the first column of the Routh table are nonzero and have the same sign.

An elementary proof of the Routh-Hurwitz criterion can be found in references [18] and [19] .

Lemma 2. Suppose all of the eigenvalues of Equation (16) have negative real parts. Then the equilibrium solution is asymptotically stable [20] .

It is obvious that the zero solution is the equilibrium point of Equation (12). In order to facilitate the numerical analysis of this paper, we select $M=1$ , $\eta =1$ . According to Equation (12), the coefficients ${\alpha }_1,{\gamma }_0$ are $\frac{1}{2},\frac{1}{2}$ respectively.

Then the Equation (16) can be rewritten as

$\{\begin{array}{l}\{\begin{array}{l}\frac{\text{d}{x}_0}{\text{d}t}=A{x}_0+B{y}_0+1/2\delta x_1,\\ \frac{\text{d}{y}_0}{\text{d}t}=D{y}_0+C{x}_0+1/2\delta y_1.\end{array}\\ \{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=A{x}_1+B{y}_1+1/2\delta x_0,\\ \frac{\text{d}{y}_1}{\text{d}t}=D{y}_1+C{x}_1+1/2\delta y_0.\end{array}\end{array}$ (17)

The Jacobian matrix J at the equilibrium point of Equation (16) is

$J=\left(\begin{array}{cccc}A& B& 1/2\delta & 0\\ D& C& 0& 1/2\delta \\ 1/2\delta & 0& A& B\\ 0& 1/2\delta & C& D\end{array}\right).$ (18)

With aid of Maple, we obtain the characteristic polynomial of Jacobian matrix can be obtained as

$f\left(\lambda \right)=a_0{\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0,$ (19)

where $a_i\left(i=0,\cdots ,4\right)$ are coefficients of characteristic equation, which are shown as follows: $\begin{array}{l}{a}_0=1,\\ a_1=-2A-D-C,\\ a_2=-1/2{\delta }^2-BC+2DA+DC+2CA+A^2-BD,\\ a_3=1/4{\delta }^2\left(2A+C+D\right)-\left(A+C+D\right)\left(BC-CA-AD\right),\\ a_4=\left(1/4{\delta }^2-BC+AD\right)\left(1/4{\delta }^2-BD+AC\right).\end{array}$

The determinants which construct by the Routh-Hurwitz criterion are shown as follows:

$\begin{array}{l}1){\Delta }_1=a_1,\\ 2){\Delta }_2=\left|\begin{array}{cc}{a}_1& a_0\\ a_3& a_2\end{array}\right|,\\ 3){\Delta }_3=\left|\begin{array}{ccc}{a}_1& a_0& 0\\ a_3& a_2& a_1\\ a_5& a_4& a_3\end{array}\right|.\end{array}$

where $a_5$ is zero.

4. To Take a Prey-Predator Model as Example

In this section, the asymptotical stability of a prey-predator model is considered. And the complex dynamic behavior is discussed and numerical simulation is performed. We select the coefficients $a_{11},a_{12},b_1,a_{21},a_{22},b_2$ are $-1,-2,1,2,-1,-1$ respectively. Obviously, with all these coefficients, the stochastic Lotka-Volterra system (7) is a prey-predator model accurately. Then the Jacobian matrix (18) can be rewritten as

$J=\left(\begin{array}{cccc}-1/5& -2/5& 1/2\delta & 0\\ -1& 2/5& 0& 1/2\delta \\ 1/2\delta & 0& -1/5& -2/5\\ 0& 1/2\delta & 2/5& -1\end{array}\right).$ (20)

By the Mathematical software, the coefficients of characteristic equation of the Jacobian matrix J (20) are as follows.

$a_0=1$ , $a_1=12/5$ , $a_2=54/25-1/2\delta$ , $a_3=108/125-3/5{\delta }^2$ , $a_4=-9/50{\delta }^2+81/625+1/16{\delta }^4$ .

And all eigenvalues of Equation (19) are

${\lambda }_{1,2}=-\frac{\delta }{2}-\frac{3}{5},{\lambda }_{3,4}=\frac{\delta }{2}-\frac{3}{5}$ .

We can see ${\Delta }_1=a_1>0$ easily, So only if ${\Delta }_{2,3}>0$ , we can say all of the eigenvalues of Jacobian matrix (18) have negative real parts. The functions are constructed by Routh-Hurwitz criterion

$\{\begin{array}{l}3/5{\delta }^2-6/5\delta +108/35>0,\\ 12/25{\delta }^4+18/25{\delta }^3-3456/1250{\delta }^2+648/625\delta +47502/15125>0.\end{array}$ (21)

i.e. $-2/5<\delta <6/5$ .

Theorem 1. The necessary and sufficient condition of locally asymptotically stable at the zero equilibrium solution is $-2/5<\delta <6/5$ .

Prof. Depend on the lemma2, we know if all of the eigenvalues of Equation (19) have negative real parts. Then the equilibrium solution is asymptotically stable. Obviously, ${\lambda }_{1,2,3,4}$ have negative real parts when $-2/5<\delta <6/5$ , so we get it.

Theorem 2. The sufficient condition of globally asymptotically stable at the zero equilibrium solution is $-2/5<\delta <6/5$ .

Prof. (Using the Lyapunov function)

We choose the Lyapunov function $V\left(x,y\right)$ as $V\left(x,y\right)={\left(Dx-By\right)}^2+\left(AD-BC\right)x^2$ , Obviously $AD-BC>0$ , So Lyapunov function $V\left(x,y\right)$ is always positive. Then the total derivative of Equation (19) is ${\frac{\text{d}V}{\text{d}t}|}_{(19)}=\frac{\text{d}V}{\text{d}x}\ast \frac{\text{d}x}{\text{d}t}+\frac{\text{d}V}{\text{d}y}\ast \frac{\text{d}y}{\text{d}t}=4x^2\left(AD-BC\right)\left(D^{2}x-BDy\right)\left(Ax-By+1/2\delta x\right)$ ,

At the aid of the image of quadratic function, it’s not difficult to get at the condition of $-2/5<\delta <6/5$ , ${\frac{\text{d}V}{\text{d}t}|}_{(19)}$ less than zero constantly.

5. Conclusion

Orthogonal polynomial approximation is applied to study the stability in a stochastic Lotka-Volterra system with random parameter. Analysis shows that orthogonal polynomial approximation is effective to reduce the stochastic Lotka-Volterra system with random parameter into its equivalent deterministic system. The linear stability theory and Routh-Hurwitz criterion for nonlinear deterministic systems are applied to the equivalent one. By the mathematics analysis method, we have discovered that as the random intensity or statistical parameters of random variable are increased, results are different from the deterministic system which are characterized the stability of realistic models accurately.

Founding

This work is supported by National Natural Science Foundation of China (12001471).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lotka, A.J. (1920) Undamped Oscillations Derived From the Law of Mass Action. Journal of American Chemical Society, 42, 1595-1599. https://doi.org/10.1021/ja01453a010
[2]	Volterra, V. (1931) Variations and Fluctuations of the Number of Individuals in Animal Species Living Together. In: Chapman, R.N., Ed., Animal Ecology, McGraw-Hill, New York, 409-448.
[3]	Hou, X.J., Wang, B. and Zhang, Z.C. (2014) Mutual Inclusion in a Nonlocal Competitive Lotka Volterra System. Japan Journal of Industrial and Applied Mathematics, 31, 87-110. https://doi.org/10.1007/s13160-013-0126-0
[4]	Maharaj, A. and Leach, P.G.L. (2014) Application of Symmetry and Singularity Analyses to Mathematical Models of Biological Systems. Mathematics and Computers in Simulation, 96, 104-123. https://doi.org/10.1016/j.matcom.2013.06.005
[5]	Wang, M.-X. and Ma, Y.-Q. (2014) Population Evolution in Mutualistic Lotka-Volterra System with Spatial Diffusion. Physica A: Statistical Mechanics and Its Applications, 395, 228-235. https://doi.org/10.1016/j.physa.2013.10.019
[6]	Luo, Z.G., Luo, L.P. and Zeng, Y.H. (2014) Positive Periodic Solution for the Generalized Neutral Differential Equation with Multiple Delays and Impulse. Journal of Applied Mathematics, 2014, Article ID: 592513. https://doi.org/10.1155/2014/592513
[7]	Zhu, H., Campbell, S.A. and Wolkowicz, G.S.K. (2002) Bifurcation Analysis of a Predator-Prey System with Nonmonotonic Functional Response. SIAM Journal on Applied Mathematics, 63, 636-682. https://doi.org/10.1137/S0036139901397285
[8]	Zhang, H., Li, Y.Q. and Jing, B. (2014) Global Stability of Almost Periodic Solution of Multispecies Mutualism System with Time Delays and Impulsive Effects. Applied Mathematics and Computation, 232, 1138-1150. https://doi.org/10.1016/j.amc.2014.01.131
[9]	Gu, X., Wang, H.C., Wong, P.J.Y. and Xia, Y.H. (2014) Existence and Stability of Periodic Solution to Delayed Nonlinear Differential Equations. Abstract and Applied Analysis, 2014, Article ID: 156948. https://doi.org/10.1155/2014/156948
[10]	Shi, C.L. and Chen, X.Y. (2013) Global Asymptotic Stability in Non-Autonomous Lotka-Volterra Competitive System with Infinite Delay and Feedback Controls. Annals of Differential Equations, 29, 312-323.
[11]	Bobieński, M. (2011) Invariant Torus in 3D Lotka-Volterra System Appearing after Perturbation of Hopf Center. Qualitative Theory of Dynamical Systems, 10, 51-64. https://doi.org/10.1007/s12346-011-0037-x
[12]	Liu, H.Y., Luo, J.G. and Zhao, X.H. (2003) On Criteria for Global Stability of N-Dimensional Lotka-Volterra Systems. Annals of Differential Equations, 19, 343-351.
[13]	Khaminskii, R.Z., Klebaner, F.C. and Liptser, R. (2003) Some Results on the Lotka-Volterra Model and its Small Random Perturbations. Acta Applicandae Mathematics, 78, 201-206. https://doi.org/10.1023/A:1025737516835
[14]	Sun, C., Lin, Y. and Han, M. (2006) Stability and Hopf Bifurcation for an Epidemic Disease Model with Delay. Chaos, Solitons & Fractals, 30, 204-216. https://doi.org/10.1016/j.chaos.2005.08.167
[15]	Xu, C.J. and Wu, Y.S. (2013) Dynamics in a Lotka-Volterra Predator-Prey Model with Time-Varying Delays. Abstract and Applied Analysis, 2013, Article ID: 956703. https://doi.org/10.1155/2013/956703
[16]	Jun, Z. and Kim, C.-G. (2014) Positive Solutions for a Lotka-Volterra Prey-Predator Model with Cross-Diffusion of Fractional Type. Results in Mathematics, 65, 293-320. https://doi.org/10.1007/s00025-013-0346-2
[17]	Meng, X.Z., Liu, R. and Zhang, T.H. (2014) Adaptive Dynamics for a Non-Autonomous Lotka-Volterra Model with Size-Selective Disturbance. Nonlinear Analysis: Real World Applications, 16, 202-213. https://doi.org/10.1016/j.nonrwa.2013.09.019
[18]	Meinsma, G. (1995) Elementary Proof of the Routh-Hurwitz Test. Systems & Control Letters, 25, 237-242. https://doi.org/10.1016/0167-6911(94)00089-E
[19]	Gantmacher, F.R. (2000) Theory of Matrices, Volume 2. Chelsea Publishing, New York.
[20]	Marsden, J.E., Sirovich, S. and Antman, S.S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, New York.