# 5D Nonlinear Dynamic Evolutionary System in Real Estate Market

Authors: Jingyuan Zhang
Journal: Complexity
 (2021)

---

## Abstract

In this paper, we propose a new predator-prey nonlinear dynamic evolutionary model of real estate enterprises considering the large, medium, and small real estate enterprises for three different prey teams. A 5D predator-prey nonlinear dynamic evolutionary system in the real estate market is established, where the large, medium, and small real estate enterprises correspond to three differential equations, provincial and local officials, and the central government correspond to the other two differential equations. Nonlinear dynamic analysis on a 5D predator-prey evolutionary system in the real estate market, containing the analysis of equilibrium points and stabilities, is made. The corresponding discrete system is simulated, and the simulation results about Lyapunov spectrum, bifurcation diagram, sequence diagram, and phase diagram are given. Compared with the work of Yang et al. in which all real estate enterprises corresponded to one differential equation, in our proposed model, the large, medium, and small real estate enterprises correspond to three differential equations which is more accordant with the specific circumstance of real estate companies.


---

## Body

## 1. Introduction

Since the reform in 1978, China began to implement the system of market economy which brings the upward trend in real estate investment. According to the statistics of the National Administration of Industry and Commerce in 2018, there are altogether 97,000 real estate developers registered in the National Administration of Industry and Commerce. The total number of permanent residents in China, consisting of more than 500 million urban residents and 300 million migrant workers, is over 8.8 million. There are more than 90,000 real estate enterprises with a permanent urban population of more than 800 million. According to data released by the National Bureau of Statistics, housing prices in major cities are still rising mostly, especially for new commercial residential buildings. So in China, real estate enterprises are very large groups and are also very important.Nonlinear characteristics exist in systems with various research directions, such as chaotic circuit [1–7], neural network [8–15], and image encryption [16–18]. Competition of enterprises also has nonlinear characteristics. Puu [19] studied the adjustment process of three oligopolists, under Coumot and Stackelberg action. It is demonstrated that, with an isoelastic demand function and constant marginal costs, the system can result in periodic or in chaotic behavior. Based on the analysis, several useful issues are investigated either analytically or numerically. Agiza et al. [20] investigate the dynamics of a nonlinear discrete-time duopoly game, where the players have heterogeneous expectations. Two players with different expectations are considered: one is boundedly rational and the other thinks with adaptive expectations. The stability conditions of the equilibria were discussed. How the dynamics of the game depend on the model parameters was found. They demonstrate that as some parameters of the game are varied, the stability of Nash equilibrium is lost by period-doubling bifurcation. The chaotic features were justified numerically via computing Lyapunov exponents, sensitive dependence on initial conditions, and the fractal dimension. Ma and Pu [21] modeled a dynamic triopoly game characterized by firms with different expectations by three-dimensional nonlinear difference equations, where the market has quadratic inverse demand function and the firm possesses cubic total cost function. The local stability of Nash equilibrium was studied. Numerical simulations were presented to show that the triopoly game model behaved chaotically with the variation of the parameters. Elsadany [22, 23] used three nonlinear difference equations to study the dynamic Cournot game with the characteristics of three bounded rational actors and analyzed the stability of the system. The global complexity analysis is helpful for behavior taking some effective measures, avoiding the collapse of the output dynamic competition game, and obtaining some practical and theoretical significance in the practice. Ahmed and Agiza [24] derived the dynamical system of n competitors in a Cournot game and studied the stability of its fixed point “Nash equilibrium.” The effect of a modification of the price demand relation was pointed out. Zhao and Lv [25] studied a three-species food chain model with a Beddington-De Angelis functional response. The equilibrium states of the system were identified and their stability was analyzed analytically. The simulation results showed chaotic long-term behavior over a broad range of parameters. Elettreby and Hassan [26] proposed two different versions of the multiteam model where a team of two firms compete with another team. The firms in each team help each other. The equilibrium solutions and the conditions of their local asymptotic stability were studied. Elettreby and Mansour [27] studied an incomplete information dynamical system. Then, they suggested a modification of this system which was applied to the standard Cournot game. The equilibrium solutions and the conditions of their locally asymptotic stability for the static and the dynamic in monopoly and duopoly cases were studied. Elettreby [28] proposed a new multiteam prey-predator model, in which the prey teams help each other. Its local stability was studied. In the absence of a predator, there was no help between the prey teams. The global stability and persistence of the model without help were studied. Upadhyay et al. [29] studied the influence of top predator interference on the dynamics of food chain models including intermediate predators and top predators and found that there were different types of attraction sets including chaos. Elettreby and El-Metwally [30] applied the multiteam concept to the predator model and studied the global stability and persistence of the model without help. Liu et al. [31] investigated the appointed-time consensus problem under directed and periodical switching topologies. From a motion-planning perspective, a novel distributed appointed-time algorithm was developed for a multiagent system with double-integrator dynamics. Cotter and Roll [32] studied a comparative anatomy of residential REITs and private real estate markets: returns, risks, and distributional characteristics. REITs have somewhat less market risk than equity. Residential REIT characteristics differ from those of the S&P/Case–Shiller (SCS) private real estate markets. This is partly attributable to methods of index construction. They suggested that investment in residential real estate is far more risky than what might be inferred from the widely followed SCS series. The evolution of the real estate market is very complex and changeable, similar to the biological evolutionary process, falling within the scope of complexity. Motivated by this, the study uses the biological evolutionary process to explain the economic phenomenon in the real estate market. However, little research has quantitatively analyzed this process in the real estate market from the biological view. Yang and Tang [33] established a model on a three-dimension predator-prey evolutionary system in the real estate market. The model involved the relationship among private enterprises, provincial and local officials, and the central government in the real estate market using the population ecology theory of mutual relations. The complex dynamical behaviors of such a predator-prey model are investigated by means of numerical simulation. However, [33] considered all private enterprises as a prey team, and just one equation corresponds to the private enterprises. In fact, according to a 2015 real estate enterprise sales ranking, the private enterprises were divided into three levels [33]: first class: sales ≥ 100 billion; the second class: 20 billion ≤ sales < 100 billion; and the third class: sales < 20 billion. The first kind of private enterprises has great strength, can assist the third type of enterprises, or join cooperation evolution and development with the second category of private enterprises. Private enterprises in the second category compete or cooperate with other companies, having their own characteristics. The third-class private enterprises have relatively weak strength, requiring the help of other companies, in order to maintain survival in the real estate market. Thus, in [33], it is relatively unspecific for considering three types of real estate companies as a prey team, which does not fit the specific circumstance of real estate companies. In this paper, we propose a new model considering the large, medium, and small real estate enterprises for three different prey teams. A 5D predator-prey evolutionary system in the real estate market is established, where the large, medium, and small real estate enterprises correspond to three differential equations, provincial and local officials, and the central government correspond to the other two differential equations. Compared with literature [33] in which all real estate enterprises corresponded to a differential equation, in our proposed model, the large, medium, and small real estate enterprises correspond to three differential equations.


## 2. Model and Nonlinear Dynamic Analysis

### 2.1. Model

Letx be the density of large real estate companies, y the density of medium real estate companies, z the density of small real estate companies, u the density of provincial and local officials as predators, and v the center government, which can be growth-oriented central leaders who are intelligent designers of institutions that moderate the predator-prey relationship. The ecosystem on the top of the predator-prey interactions is the institutional framework subject to adjustment by the central government and a, b, c, d, e, f, h are parameters. Large enterprises can help medium-sized and small enterprises, and large enterprises, medium-sized enterprises, and small enterprises are affected by each other; and α, β, and γ are their influence coefficients. The nonlinear differential equation can be established as follows:(1)dxdt=ax1−x−xu+αxyzu,dydt=by1−y−yu+βxyzu,dzdt=cz1−z−zu+γxyzu,dudt=−du2+exu+fyu+gzu,dvdt=hv−v2u,where in the absence of any predation u, each team of preys grows logistically; this is a x (1 − x), by (1 − y), cz (1 − z); and the terms −xu, −yu, and −zu denote the reduction in the prey growth rate because of the effect of the predation. The teams of preys help each other against the predator, that is, a xyzu term exists. In the absence of any prey for sustenance, the predator’s death rate results in inverse decay, that is, the term −du2. The prey’s contribution to the predator growth rate is exu, fyu, and gzu, that is proportional to the available prey as well as the size of the predator population. hv measures the value of self-reproduction of top predator v. In this paper, it is central government −v2/u denotes the effect of predator u on v.


### 2.2. Nonlinear Dynamic Analysis for the Model

#### 2.2.1. Solution of Equilibrium Points

In order to find the equilibrium points, the right side of equation (1) is set as zero, and we can obtain equation (2) as follows:(2)ax1−x−xu+αxyzu=0,by1−y−yu+βxyzu=0,cz1−z−zu+γxyzu=0,−du2+exu+fyu+gzu=0,hv−v2u=0.Solving equation (2), we can get the equilibrium points Ej=x,y,z,u,v (j = 1, 2, …, 12) as follows:(3)E10,0,dcdc+g,gcdc+g,0,E20,0,dcdc+g,gcdc+g,hgcdc+g,E30,dbdb+f,0,fbdb+f,0,E40,dbdb+f,0,fbdb+f,hfbdb+f,E5dada+e,0,0,eada+e,0,E6dada+e,0,0,eada+e,heada+e,E70,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,0,E80,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,hbcf+hbcgbcd+cf+bg,E9acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,0,E10acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,hace+hacgacd+ce+ag,E11abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,0,E12abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,habe+habfabd+be+af.


#### 2.2.2. Stability Analysis of Equilibrium Point

Jacobian matrix of the system described by (1) can be calculated easily as shown in(4)J=a1−2x−u1−αyzαxzuαxyu−x1−αyz0βyzub1−2y−u1−βxzβxyu−y1−βxz0γyzuγxzuc1−2z−u1−γxy−z1−γxy0eufugu−2du+ex+fy+gz0000v2u2hu−2vu.Proposition 1.
Whenv is equal to 0, the equilibrium points (E1, E3, E5, E7, E9, and E11) are unstable.Proof.
From the 5th expression of equation (2), we can get vt = 0 or vt = hu (t). Substituting v=0, v=hu into equation (4), we can get (5)J=±ha1−2x−u1−αyzαxzuαxyu−x1−αyzβyzub1−2y−u1−βxzβxyu−y1−βxzγyzuγxzuc1−2z−u1−γxy−z1−γxyeufugu−2du+ex+fy+gz,where “±” means that “+” corresponds to v=0 and “−” corresponds to v=hu. It is known that when v=0, one of the eigenvalues of the Jacobian matrix is h. Because h > 0, the equilibrium point is not stable.
That is the end of the proof.
For the equilibrium pointE20,0,dc/dc+g,gc/dc+g,hgc/dc+g, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−gc/dc+g,b−gc/dc+g,−c,−dgc/dc+g. Obviously, all eigenvalues are negative and the equilibrium point E20,0,dc/dc+g,gc/dc+g,hgc/dc+g is locally asymptotically stable when the following conditions are met:(6)a<gcdc+g,b<gcdc+g.
For the equilibrium pointE40,db/db+f,0,fb/db+f,hfb/db+f, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−fb/db+f,c−fb/db+f,−b,−dfb/db+f. Obviously, all eigenvalues are negative and the equilibrium point E40,db/db+f,0,fb/db+f,hfb/db+f is locally asymptotically stable when the following conditions are satisfied:(7)a<fbdb+f,b<fbdb+f.
For the equilibrium pointE6da/da+e,0,0,ea/da+e,hea/da+e, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ==−h,−a,b−ea/da+e,c−ea/da+e,−dea/da+e. Obviously, all eigenvalues are negative and the equilibrium point E6da/da+e,0,0,ea/da+e,hea/da+e is locally asymptotically stable when the following conditions are satisfied:(8)b<eada+e,c<eada+e.
For the equilibrium pointE80,bcd−cg+bg/bcd+cf+bg,bcd+cf−bf/bcd+cf+bg,bcf+bcg/bcd+cf+bg,hbcf+hbcg/bcd+cf+bg, from the Jacobian matrix (4), the following equation can be obtained:(9)J=−ha−u1−αyz000yzub1−2y−u0−yyzu0c1−2z−u−zeufugu−2du+fy+gz=−ha−u1−αyzb1−2y−u0−y0c1−2z−u−zfugu−du.
The Jacobian determinant is as follows:(10)J−λE==−h−λa−u1−αyz−λb1−2y−u−λ0−y0c1−2z−u−λ−zfugu−du−λ=−−h−λA−λλ3+du−B−Cλ2+CB+guz+fuy−duB−duCλ+duBC−guzB+fuyC,where A=a−u1−αyz;B=b1−2y−u;C=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part:(11)A<0,du−B−C>0,CB+guz+fuy−duB−duC>0,duBC−guzB+fuyC>0,du−B−CCB+guz+fuy−duB−duC>duBC−guzB+fuyC.
For the equilibrium pointE10acd−cg+ag/acd+ce+ag,0,acd+ce−ae/acd+ce+ag,ace+acg/acd+ce+ag,hace+hacg/acd+ce+ag, from the Jacobian matrix (4), equation (12) can be obtained:(12)J=−ha1−2x−uxzu0−x0b−u1−βxz000xzuc1−2z−u−zeufugu−du=−hb−u1−βxza1−2x−u0−x0c1−2z−u−zeugu−du.
The Jacobian determinant is as follows:(13)J−λE=−h−λb−u1−βxz−λa1−2x−u−λ0−x0c1−2z−u−λ−zeugu−du−λ=−−h−λB1−λλ3+du−A1−C1λ2+A1C1+zgu+eux−duA1−duC1λ+duA1C1−zguA1+euxC1,where A1=a1−2x−u;B1=b−u1−βxz;C1=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part: (14)B1<0,du−A1−C1>0,A1C1+zgu+eux−duA1−duC1>0,duA1C1−zguA1+euxC1>0,du−A1−C1A1C1+zgu+eux−duA1−duC1>duA1C1−zguA1+euxC1.
For the equilibrium pointE12abd−bf+af/abd+be+af,abd+be−ae/abd+be+af,0,abe+abf/abd+be+af,habe+habf/abd+be+af, from the Jacobian matrix (4), equation (15) can be obtained:(15)J=−ha1−2x−u0xyu−x0b1−2y−uxyu−y00c−u1−γxy0eufugu−du=−hc−u1−γxya1−2x−u0−x0b1−2y−u−yeufu−du.
The Jacobian determinant is as follows:(16)J−λE=c−u1−γxy−λa1−2x−u−λ0−x0b1−2y−u−λ−yeufu−du−λ=−−h−λC2−λλ3+du−A2−B2λ2+A2B2+eux+fuy−duB2−duA2λ+duA2B2−euxB2+fuyA2=0,where A2=a1−2x−u;B2=b1−2y−u;C2=c−u1−γxy. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part: (17)C2<0,du−A2−B2>0,A2B2+eux+fuy−duB2−duA2>0,duA2B2−euxB2+fuyA2>0,du−B2−A2A2B2+eux+fuy−duB2−duA2>duA2B2−euxB2+fuyA2.
We have discussed the stability of all the equilibria of the system. Now, we validate the above conclusions, and the validation method is as follows: at first, we suppose the coefficients area, b, c, d, e, f, g, h, α, β, γ; then, by the Jacobi matrix, we find its eigenvalues. According to the characteristics of the positive and negative polarity of the eigenvalues, we can judge the stability of equilibrium points and then substitute into the supposed coefficient of a, b, c, d, e, f, g, h, α, β, γ into stability judgment inequalities (6), (8), (11), (14), and (17) and verify the inequality of stability judgment conditions. For example, if the characteristic roots of the Jacobian matrix are all negative corresponding to the preassumed coefficients, then the equilibrium point is stable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is true, the theoretical analysis process is proved to be correct. On the contrary, if the characteristic roots of the Jacobian matrix are not all negative corresponding to the preassumed coefficients, then the equilibrium point is unstable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is not true, the theoretical analysis process is proved also to be correct. Here, Table 1 shows the equilibrium points (viz, the solution of equation (2)), the corresponding eigenvalue, the stability, and the conditional inequality systems true or no when the preassumed coefficients are a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, γ = 0.8.Table 1
The equilibrium points, the corresponding eigenvalue, the stability, and the conditional inequality systems true or no.
Equilibrium pointsEigenvaluesStabilityIs the inequality set true?(0, 0, 0.9158, 0.2466, 0)−2.6495; −0.3545; 0.7534; 0.8934; 2UnstableNo(0, 0, 0.9158, 0.2466, 0.4931)−2;−0.3545;−2.6495; 0.7534; 0.8934UnstableNo(0, 0.9368, 0, 0.0721, 0)−0.1007; −1.061; 0.9279; 2.8579; 2UnstableNo(0, 0.9368, 0, 0.0721, 0.1441)−2;−1.061; 0.1007; 0.9279; 2.8579UnstableNo(0.8609, 0, 0, 0.1391, 0)−0.8217;−0.22; 1.0009; 2.7909; 2UnstableNo(0.8609, 0, 0, 0.1391, 0.2781)−2; −0.22; −0.8217; 1.0009; 2.7909UnstableNo(0, 0.7381, 0.8981, 0.2986, 0)−2.589; −0.4959;−0.7761; 0.7113; 2UnstableNo(0, 0.7381, 0.8981, 0.2986, 0.5971)−2;−−0.7761;−0.4959;−2.589; 0.7113UnstableNo(0.6563, 0, 0.8827, 0.3437, 0)−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945; 2UnstableNo(0.6563, 0, 0.8827, 0.3437, 0.6873)−2;−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945UnstableNo(0.806, 0.8298, 0, 0.194, 0)−0.3516; −0.7237; −0.9289; 2.8398; 2UnstableNo(0.806, 0.8298, 0, 0.194, 0.3881)−2; −0.3516; −0.9289;−0.7237; 2.8398UnstableNo


## 2.1. Model

Letx be the density of large real estate companies, y the density of medium real estate companies, z the density of small real estate companies, u the density of provincial and local officials as predators, and v the center government, which can be growth-oriented central leaders who are intelligent designers of institutions that moderate the predator-prey relationship. The ecosystem on the top of the predator-prey interactions is the institutional framework subject to adjustment by the central government and a, b, c, d, e, f, h are parameters. Large enterprises can help medium-sized and small enterprises, and large enterprises, medium-sized enterprises, and small enterprises are affected by each other; and α, β, and γ are their influence coefficients. The nonlinear differential equation can be established as follows:(1)dxdt=ax1−x−xu+αxyzu,dydt=by1−y−yu+βxyzu,dzdt=cz1−z−zu+γxyzu,dudt=−du2+exu+fyu+gzu,dvdt=hv−v2u,where in the absence of any predation u, each team of preys grows logistically; this is a x (1 − x), by (1 − y), cz (1 − z); and the terms −xu, −yu, and −zu denote the reduction in the prey growth rate because of the effect of the predation. The teams of preys help each other against the predator, that is, a xyzu term exists. In the absence of any prey for sustenance, the predator’s death rate results in inverse decay, that is, the term −du2. The prey’s contribution to the predator growth rate is exu, fyu, and gzu, that is proportional to the available prey as well as the size of the predator population. hv measures the value of self-reproduction of top predator v. In this paper, it is central government −v2/u denotes the effect of predator u on v.


## 2.2. Nonlinear Dynamic Analysis for the Model

### 2.2.1. Solution of Equilibrium Points

In order to find the equilibrium points, the right side of equation (1) is set as zero, and we can obtain equation (2) as follows:(2)ax1−x−xu+αxyzu=0,by1−y−yu+βxyzu=0,cz1−z−zu+γxyzu=0,−du2+exu+fyu+gzu=0,hv−v2u=0.Solving equation (2), we can get the equilibrium points Ej=x,y,z,u,v (j = 1, 2, …, 12) as follows:(3)E10,0,dcdc+g,gcdc+g,0,E20,0,dcdc+g,gcdc+g,hgcdc+g,E30,dbdb+f,0,fbdb+f,0,E40,dbdb+f,0,fbdb+f,hfbdb+f,E5dada+e,0,0,eada+e,0,E6dada+e,0,0,eada+e,heada+e,E70,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,0,E80,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,hbcf+hbcgbcd+cf+bg,E9acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,0,E10acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,hace+hacgacd+ce+ag,E11abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,0,E12abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,habe+habfabd+be+af.


### 2.2.2. Stability Analysis of Equilibrium Point

Jacobian matrix of the system described by (1) can be calculated easily as shown in(4)J=a1−2x−u1−αyzαxzuαxyu−x1−αyz0βyzub1−2y−u1−βxzβxyu−y1−βxz0γyzuγxzuc1−2z−u1−γxy−z1−γxy0eufugu−2du+ex+fy+gz0000v2u2hu−2vu.Proposition 1.
Whenv is equal to 0, the equilibrium points (E1, E3, E5, E7, E9, and E11) are unstable.Proof.
From the 5th expression of equation (2), we can get vt = 0 or vt = hu (t). Substituting v=0, v=hu into equation (4), we can get (5)J=±ha1−2x−u1−αyzαxzuαxyu−x1−αyzβyzub1−2y−u1−βxzβxyu−y1−βxzγyzuγxzuc1−2z−u1−γxy−z1−γxyeufugu−2du+ex+fy+gz,where “±” means that “+” corresponds to v=0 and “−” corresponds to v=hu. It is known that when v=0, one of the eigenvalues of the Jacobian matrix is h. Because h > 0, the equilibrium point is not stable.
That is the end of the proof.
For the equilibrium pointE20,0,dc/dc+g,gc/dc+g,hgc/dc+g, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−gc/dc+g,b−gc/dc+g,−c,−dgc/dc+g. Obviously, all eigenvalues are negative and the equilibrium point E20,0,dc/dc+g,gc/dc+g,hgc/dc+g is locally asymptotically stable when the following conditions are met:(6)a<gcdc+g,b<gcdc+g.
For the equilibrium pointE40,db/db+f,0,fb/db+f,hfb/db+f, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−fb/db+f,c−fb/db+f,−b,−dfb/db+f. Obviously, all eigenvalues are negative and the equilibrium point E40,db/db+f,0,fb/db+f,hfb/db+f is locally asymptotically stable when the following conditions are satisfied:(7)a<fbdb+f,b<fbdb+f.
For the equilibrium pointE6da/da+e,0,0,ea/da+e,hea/da+e, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ==−h,−a,b−ea/da+e,c−ea/da+e,−dea/da+e. Obviously, all eigenvalues are negative and the equilibrium point E6da/da+e,0,0,ea/da+e,hea/da+e is locally asymptotically stable when the following conditions are satisfied:(8)b<eada+e,c<eada+e.
For the equilibrium pointE80,bcd−cg+bg/bcd+cf+bg,bcd+cf−bf/bcd+cf+bg,bcf+bcg/bcd+cf+bg,hbcf+hbcg/bcd+cf+bg, from the Jacobian matrix (4), the following equation can be obtained:(9)J=−ha−u1−αyz000yzub1−2y−u0−yyzu0c1−2z−u−zeufugu−2du+fy+gz=−ha−u1−αyzb1−2y−u0−y0c1−2z−u−zfugu−du.
The Jacobian determinant is as follows:(10)J−λE==−h−λa−u1−αyz−λb1−2y−u−λ0−y0c1−2z−u−λ−zfugu−du−λ=−−h−λA−λλ3+du−B−Cλ2+CB+guz+fuy−duB−duCλ+duBC−guzB+fuyC,where A=a−u1−αyz;B=b1−2y−u;C=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part:(11)A<0,du−B−C>0,CB+guz+fuy−duB−duC>0,duBC−guzB+fuyC>0,du−B−CCB+guz+fuy−duB−duC>duBC−guzB+fuyC.
For the equilibrium pointE10acd−cg+ag/acd+ce+ag,0,acd+ce−ae/acd+ce+ag,ace+acg/acd+ce+ag,hace+hacg/acd+ce+ag, from the Jacobian matrix (4), equation (12) can be obtained:(12)J=−ha1−2x−uxzu0−x0b−u1−βxz000xzuc1−2z−u−zeufugu−du=−hb−u1−βxza1−2x−u0−x0c1−2z−u−zeugu−du.
The Jacobian determinant is as follows:(13)J−λE=−h−λb−u1−βxz−λa1−2x−u−λ0−x0c1−2z−u−λ−zeugu−du−λ=−−h−λB1−λλ3+du−A1−C1λ2+A1C1+zgu+eux−duA1−duC1λ+duA1C1−zguA1+euxC1,where A1=a1−2x−u;B1=b−u1−βxz;C1=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part: (14)B1<0,du−A1−C1>0,A1C1+zgu+eux−duA1−duC1>0,duA1C1−zguA1+euxC1>0,du−A1−C1A1C1+zgu+eux−duA1−duC1>duA1C1−zguA1+euxC1.
For the equilibrium pointE12abd−bf+af/abd+be+af,abd+be−ae/abd+be+af,0,abe+abf/abd+be+af,habe+habf/abd+be+af, from the Jacobian matrix (4), equation (15) can be obtained:(15)J=−ha1−2x−u0xyu−x0b1−2y−uxyu−y00c−u1−γxy0eufugu−du=−hc−u1−γxya1−2x−u0−x0b1−2y−u−yeufu−du.
The Jacobian determinant is as follows:(16)J−λE=c−u1−γxy−λa1−2x−u−λ0−x0b1−2y−u−λ−yeufu−du−λ=−−h−λC2−λλ3+du−A2−B2λ2+A2B2+eux+fuy−duB2−duA2λ+duA2B2−euxB2+fuyA2=0,where A2=a1−2x−u;B2=b1−2y−u;C2=c−u1−γxy. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part: (17)C2<0,du−A2−B2>0,A2B2+eux+fuy−duB2−duA2>0,duA2B2−euxB2+fuyA2>0,du−B2−A2A2B2+eux+fuy−duB2−duA2>duA2B2−euxB2+fuyA2.
We have discussed the stability of all the equilibria of the system. Now, we validate the above conclusions, and the validation method is as follows: at first, we suppose the coefficients area, b, c, d, e, f, g, h, α, β, γ; then, by the Jacobi matrix, we find its eigenvalues. According to the characteristics of the positive and negative polarity of the eigenvalues, we can judge the stability of equilibrium points and then substitute into the supposed coefficient of a, b, c, d, e, f, g, h, α, β, γ into stability judgment inequalities (6), (8), (11), (14), and (17) and verify the inequality of stability judgment conditions. For example, if the characteristic roots of the Jacobian matrix are all negative corresponding to the preassumed coefficients, then the equilibrium point is stable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is true, the theoretical analysis process is proved to be correct. On the contrary, if the characteristic roots of the Jacobian matrix are not all negative corresponding to the preassumed coefficients, then the equilibrium point is unstable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is not true, the theoretical analysis process is proved also to be correct. Here, Table 1 shows the equilibrium points (viz, the solution of equation (2)), the corresponding eigenvalue, the stability, and the conditional inequality systems true or no when the preassumed coefficients are a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, γ = 0.8.Table 1
The equilibrium points, the corresponding eigenvalue, the stability, and the conditional inequality systems true or no.
Equilibrium pointsEigenvaluesStabilityIs the inequality set true?(0, 0, 0.9158, 0.2466, 0)−2.6495; −0.3545; 0.7534; 0.8934; 2UnstableNo(0, 0, 0.9158, 0.2466, 0.4931)−2;−0.3545;−2.6495; 0.7534; 0.8934UnstableNo(0, 0.9368, 0, 0.0721, 0)−0.1007; −1.061; 0.9279; 2.8579; 2UnstableNo(0, 0.9368, 0, 0.0721, 0.1441)−2;−1.061; 0.1007; 0.9279; 2.8579UnstableNo(0.8609, 0, 0, 0.1391, 0)−0.8217;−0.22; 1.0009; 2.7909; 2UnstableNo(0.8609, 0, 0, 0.1391, 0.2781)−2; −0.22; −0.8217; 1.0009; 2.7909UnstableNo(0, 0.7381, 0.8981, 0.2986, 0)−2.589; −0.4959;−0.7761; 0.7113; 2UnstableNo(0, 0.7381, 0.8981, 0.2986, 0.5971)−2;−−0.7761;−0.4959;−2.589; 0.7113UnstableNo(0.6563, 0, 0.8827, 0.3437, 0)−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945; 2UnstableNo(0.6563, 0, 0.8827, 0.3437, 0.6873)−2;−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945UnstableNo(0.806, 0.8298, 0, 0.194, 0)−0.3516; −0.7237; −0.9289; 2.8398; 2UnstableNo(0.806, 0.8298, 0, 0.194, 0.3881)−2; −0.3516; −0.9289;−0.7237; 2.8398UnstableNo


## 2.2.1. Solution of Equilibrium Points

In order to find the equilibrium points, the right side of equation (1) is set as zero, and we can obtain equation (2) as follows:(2)ax1−x−xu+αxyzu=0,by1−y−yu+βxyzu=0,cz1−z−zu+γxyzu=0,−du2+exu+fyu+gzu=0,hv−v2u=0.Solving equation (2), we can get the equilibrium points Ej=x,y,z,u,v (j = 1, 2, …, 12) as follows:(3)E10,0,dcdc+g,gcdc+g,0,E20,0,dcdc+g,gcdc+g,hgcdc+g,E30,dbdb+f,0,fbdb+f,0,E40,dbdb+f,0,fbdb+f,hfbdb+f,E5dada+e,0,0,eada+e,0,E6dada+e,0,0,eada+e,heada+e,E70,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,0,E80,bcd−cg+bgbcd+cf+bg,bcd+cf−bfbcd+cf+bg,bcf+bcgbcd+cf+bg,hbcf+hbcgbcd+cf+bg,E9acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,0,E10acd−cg+agacd+ce+ag,0,acd+ce−aeacd+ce+ag,ace+acgacd+ce+ag,hace+hacgacd+ce+ag,E11abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,0,E12abd−bf+afabd+be+af,abd+be−aeabd+be+af,0,abe+abfabd+be+af,habe+habfabd+be+af.


## 2.2.2. Stability Analysis of Equilibrium Point

Jacobian matrix of the system described by (1) can be calculated easily as shown in(4)J=a1−2x−u1−αyzαxzuαxyu−x1−αyz0βyzub1−2y−u1−βxzβxyu−y1−βxz0γyzuγxzuc1−2z−u1−γxy−z1−γxy0eufugu−2du+ex+fy+gz0000v2u2hu−2vu.Proposition 1.
Whenv is equal to 0, the equilibrium points (E1, E3, E5, E7, E9, and E11) are unstable.Proof.
From the 5th expression of equation (2), we can get vt = 0 or vt = hu (t). Substituting v=0, v=hu into equation (4), we can get (5)J=±ha1−2x−u1−αyzαxzuαxyu−x1−αyzβyzub1−2y−u1−βxzβxyu−y1−βxzγyzuγxzuc1−2z−u1−γxy−z1−γxyeufugu−2du+ex+fy+gz,where “±” means that “+” corresponds to v=0 and “−” corresponds to v=hu. It is known that when v=0, one of the eigenvalues of the Jacobian matrix is h. Because h > 0, the equilibrium point is not stable.
That is the end of the proof.
For the equilibrium pointE20,0,dc/dc+g,gc/dc+g,hgc/dc+g, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−gc/dc+g,b−gc/dc+g,−c,−dgc/dc+g. Obviously, all eigenvalues are negative and the equilibrium point E20,0,dc/dc+g,gc/dc+g,hgc/dc+g is locally asymptotically stable when the following conditions are met:(6)a<gcdc+g,b<gcdc+g.
For the equilibrium pointE40,db/db+f,0,fb/db+f,hfb/db+f, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ=−h,a−fb/db+f,c−fb/db+f,−b,−dfb/db+f. Obviously, all eigenvalues are negative and the equilibrium point E40,db/db+f,0,fb/db+f,hfb/db+f is locally asymptotically stable when the following conditions are satisfied:(7)a<fbdb+f,b<fbdb+f.
For the equilibrium pointE6da/da+e,0,0,ea/da+e,hea/da+e, the eigenvalues of the corresponding Jacobian matrix can be obtained as λ==−h,−a,b−ea/da+e,c−ea/da+e,−dea/da+e. Obviously, all eigenvalues are negative and the equilibrium point E6da/da+e,0,0,ea/da+e,hea/da+e is locally asymptotically stable when the following conditions are satisfied:(8)b<eada+e,c<eada+e.
For the equilibrium pointE80,bcd−cg+bg/bcd+cf+bg,bcd+cf−bf/bcd+cf+bg,bcf+bcg/bcd+cf+bg,hbcf+hbcg/bcd+cf+bg, from the Jacobian matrix (4), the following equation can be obtained:(9)J=−ha−u1−αyz000yzub1−2y−u0−yyzu0c1−2z−u−zeufugu−2du+fy+gz=−ha−u1−αyzb1−2y−u0−y0c1−2z−u−zfugu−du.
The Jacobian determinant is as follows:(10)J−λE==−h−λa−u1−αyz−λb1−2y−u−λ0−y0c1−2z−u−λ−zfugu−du−λ=−−h−λA−λλ3+du−B−Cλ2+CB+guz+fuy−duB−duCλ+duBC−guzB+fuyC,where A=a−u1−αyz;B=b1−2y−u;C=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part:(11)A<0,du−B−C>0,CB+guz+fuy−duB−duC>0,duBC−guzB+fuyC>0,du−B−CCB+guz+fuy−duB−duC>duBC−guzB+fuyC.
For the equilibrium pointE10acd−cg+ag/acd+ce+ag,0,acd+ce−ae/acd+ce+ag,ace+acg/acd+ce+ag,hace+hacg/acd+ce+ag, from the Jacobian matrix (4), equation (12) can be obtained:(12)J=−ha1−2x−uxzu0−x0b−u1−βxz000xzuc1−2z−u−zeufugu−du=−hb−u1−βxza1−2x−u0−x0c1−2z−u−zeugu−du.
The Jacobian determinant is as follows:(13)J−λE=−h−λb−u1−βxz−λa1−2x−u−λ0−x0c1−2z−u−λ−zeugu−du−λ=−−h−λB1−λλ3+du−A1−C1λ2+A1C1+zgu+eux−duA1−duC1λ+duA1C1−zguA1+euxC1,where A1=a1−2x−u;B1=b−u1−βxz;C1=c1−2z−u. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of the system (1) at the equilibrium point have a negative real part: (14)B1<0,du−A1−C1>0,A1C1+zgu+eux−duA1−duC1>0,duA1C1−zguA1+euxC1>0,du−A1−C1A1C1+zgu+eux−duA1−duC1>duA1C1−zguA1+euxC1.
For the equilibrium pointE12abd−bf+af/abd+be+af,abd+be−ae/abd+be+af,0,abe+abf/abd+be+af,habe+habf/abd+be+af, from the Jacobian matrix (4), equation (15) can be obtained:(15)J=−ha1−2x−u0xyu−x0b1−2y−uxyu−y00c−u1−γxy0eufugu−du=−hc−u1−γxya1−2x−u0−x0b1−2y−u−yeufu−du.
The Jacobian determinant is as follows:(16)J−λE=c−u1−γxy−λa1−2x−u−λ0−x0b1−2y−u−λ−yeufu−du−λ=−−h−λC2−λλ3+du−A2−B2λ2+A2B2+eux+fuy−duB2−duA2λ+duA2B2−euxB2+fuyA2=0,where A2=a1−2x−u;B2=b1−2y−u;C2=c−u1−γxy. Using the Routh Hurwitz condition, the following condition is necessary and sufficient when all roots of the characteristic equation of system (1) at the equilibrium point have a negative real part: (17)C2<0,du−A2−B2>0,A2B2+eux+fuy−duB2−duA2>0,duA2B2−euxB2+fuyA2>0,du−B2−A2A2B2+eux+fuy−duB2−duA2>duA2B2−euxB2+fuyA2.
We have discussed the stability of all the equilibria of the system. Now, we validate the above conclusions, and the validation method is as follows: at first, we suppose the coefficients area, b, c, d, e, f, g, h, α, β, γ; then, by the Jacobi matrix, we find its eigenvalues. According to the characteristics of the positive and negative polarity of the eigenvalues, we can judge the stability of equilibrium points and then substitute into the supposed coefficient of a, b, c, d, e, f, g, h, α, β, γ into stability judgment inequalities (6), (8), (11), (14), and (17) and verify the inequality of stability judgment conditions. For example, if the characteristic roots of the Jacobian matrix are all negative corresponding to the preassumed coefficients, then the equilibrium point is stable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is true, the theoretical analysis process is proved to be correct. On the contrary, if the characteristic roots of the Jacobian matrix are not all negative corresponding to the preassumed coefficients, then the equilibrium point is unstable. Subsequently, the preassumed coefficients are replaced into the inequality system for stability judgment condition, and if the inequality set is not true, the theoretical analysis process is proved also to be correct. Here, Table 1 shows the equilibrium points (viz, the solution of equation (2)), the corresponding eigenvalue, the stability, and the conditional inequality systems true or no when the preassumed coefficients are a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, γ = 0.8.Table 1
The equilibrium points, the corresponding eigenvalue, the stability, and the conditional inequality systems true or no.
Equilibrium pointsEigenvaluesStabilityIs the inequality set true?(0, 0, 0.9158, 0.2466, 0)−2.6495; −0.3545; 0.7534; 0.8934; 2UnstableNo(0, 0, 0.9158, 0.2466, 0.4931)−2;−0.3545;−2.6495; 0.7534; 0.8934UnstableNo(0, 0.9368, 0, 0.0721, 0)−0.1007; −1.061; 0.9279; 2.8579; 2UnstableNo(0, 0.9368, 0, 0.0721, 0.1441)−2;−1.061; 0.1007; 0.9279; 2.8579UnstableNo(0.8609, 0, 0, 0.1391, 0)−0.8217;−0.22; 1.0009; 2.7909; 2UnstableNo(0.8609, 0, 0, 0.1391, 0.2781)−2; −0.22; −0.8217; 1.0009; 2.7909UnstableNo(0, 0.7381, 0.8981, 0.2986, 0)−2.589; −0.4959;−0.7761; 0.7113; 2UnstableNo(0, 0.7381, 0.8981, 0.2986, 0.5971)−2;−−0.7761;−0.4959;−2.589; 0.7113UnstableNo(0.6563, 0, 0.8827, 0.3437, 0)−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945; 2UnstableNo(0.6563, 0, 0.8827, 0.3437, 0.6873)−2;−0.5767 + 0.2056i;−0.5767 − 0.2056i;−2.5361; 1.1945UnstableNo(0.806, 0.8298, 0, 0.194, 0)−0.3516; −0.7237; −0.9289; 2.8398; 2UnstableNo(0.806, 0.8298, 0, 0.194, 0.3881)−2; −0.3516; −0.9289;−0.7237; 2.8398UnstableNo


## 3. System Simulation Results

To facilitate the simulation of the system, we discretize equation (1). The system equation after discretization is shown in (18)xn+1=xn+axn1−xn−xnun+αxnynznun,yn+1=yn+byn1−yn−ynun+βxnynznun,zn+1=zn+czn1−zn−znun+γxnynznun,un+1=un−du2n+exnun+fynun+gznun,vn+1=vn+hvn−v2nun.Using MATLAB simulation, we can obtain a Lyapunov spectrum and bifurcation diagram about parametersα, β, γ and output sequence diagrams as follows:


### 3.1.α Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters α whose results are shown in Figures 1 and 2. From Figure 1, the following conclusions can be obtained: when α < 0.8, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the chaos is generated in the system; when 0.8 < α < 1, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the periodic window is generated in the system; when α > 1, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the periodic window turns into chaos in the system. In short, as α increases, the system experiences chaos-period-chaos. Figure 2 is the flip bifurcation diagram versus α. According to Figure 2, when α is used as 0.05, 0.95, and 2, we can obtain sequence diagrams shown in Figure 3–5, respectively, and phase diagrams shown in Figures 6–9, respectively.Figure 1
Lyapunov spectrum versusα when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5.Figure 2
Flip bifurcation diagram versus parametersα when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5.Figure 3
x–n sequence diagram which denotes that the system is chaotic when α = 0.05.Figure 4
x–n sequence diagram which denotes that the system is periodic when α = 0.95.Figure 5
x–n sequence diagram which denotes that the system is chaotic when α = 2.Figure 6
x–y phase diagram which denotes that the system is chaotic when α = 0.05.Figure 7
y-z phase diagram which denotes that the system is chaotic when α = 0.05.Figure 8
x-y phase diagram which denotes that the system is periodic when α = 0.95.Figure 9
x-y phase diagram which denotes that the system is chaotic when α = 2.


### 3.2.β Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters β whose results are shown in Figures 10 and 11. From Figure 10, the following conclusions can be obtained: when β < 1.12, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the chaos is generated in the system; when 1.12 < β < 1.24, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the periodic window is generated in the system; when 1.24 < β < 3.5, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the periodic window turns into chaos in the system; when β > 3.5, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the chaos turns into a periodic window in the system.Figure 10
Lyapunov spectrum versusβ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5.Figure 11
Flip bifurcation diagram versus parametersβ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5.In short, asβ increases, the system experiences chaos-period-chaos-period. Figure 11 is flip bifurcation diagram versus β. According to Figure 11, when β is used as 0.7, 1.2, and 3, we can obtain sequence diagrams shown in Figures 12–14, respectively, and phase diagrams shown in Figures 15–18, respectively.Figure 12
x–n sequence diagram which denotes that the system is chaotic when β = 0.7.Figure 13
x–n sequence diagram which denotes that the system is periodic when β = 1.2.Figure 14
x–n sequence diagram which denotes that the system is chaotic when β = 3.Figure 15
x-y phase diagram which denotes that the system is chaotic when β = 0.7.Figure 16
x–y phase diagram which denotes that the system is periodic when β = 1.2.Figure 17
x–y phase diagram which denotes that the system is chaotic when β = 3.Figure 18
y-z phase diagram which denotes that the system is chaotic when β = 3.


### 3.3.γ Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters β whose results are shown in Figures 19 and 20. From Figure 19, when γ = 0.096 and γ = 0.184, we can know LE1 = 0, LE2 < 0, and LE3 < 0. From Figure 20, it is clear that when 0 < γ < 0.096, the system has a 4-fold periodic; when 0.096 < γ < 0.2, the system has a 8-fold periodic. From Figure 19, when 0.2 < γ < 0.376, we know that LE1 > 0, LE2 < 0, and LE3 < 0, and the system is chaotic; when γ = 0.376 and γ = 0.392, we know that LE1 = 0, LE2 < 0, and LE3 < 0. Hence, from Figure 20, we know that when 0.376 < γ < 0.39, the system is periodic. When γ > 0.392, we know that LE1 > 0, LE2 < 0, and LE3 < 0, and the system is chaotic. In short, as γ increases, the system experiences a route of 4-fold periodic, 8-fold periodic, and chaos-period-chaos. Figure 20 is the flip bifurcation diagram versus β. According to Figure 20, when γ is used as 0.08, 0.15, 0.39, and 0.8, we can obtain sequence diagrams shown in Figure 21–24, respectively, and phase diagrams shown in Figure 25–29, respectively.Figure 19
Lyapunov spectrum versusγ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2.Figure 20
Flip bifurcation diagram versus parametersγ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2.Figure 21
x-n sequence diagram which denotes that the system is 4-fold periodic when γ = 0.08.Figure 22
x-n sequence diagram which denotes that the system is 8-fold periodic when γ = 0.15.Figure 23
x-n sequence diagram which denotes that the system is periodic when γ = 0.39.Figure 24
x-n sequence diagram which denotes that the system is chaotic when γ = 0.8.Figure 25
x-y phase diagram which denotes that the system is 4-fold periodic when γ = 0.08.Figure 26
x-y phase diagram which denotes that the system is 8-fold periodic when γ = 0.15.Figure 27
x-y phase diagram which denotes that the system is periodic when γ = 0.39.Figure 28
x-y phase diagram which denotes that the system is chaotic when γ = 0.8.Figure 29
y-z phase diagram which denotes that the system is chaotic when γ = 0.8.


### 3.4. Multistability of the System

Multistability means that when parameters are fixed and the initial values are different, the system presents different attractors. When usinga = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.08, we simulate the multistability of the system, and simulation results are shown in Figure 30 where the blue points denote 4-fold periodic states corresponding to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).Figure 30
Whena = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2, γ = 0.08, multistability results of the system where the blue points denote 4-fold periodic states correspond to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).When usinga = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.15, we simulate the multistability of the system, and simulation results are shown in Figure 31 where the blue points denote 8-fold periodic states corresponding to initial value (0.1, 0.1, 0.1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).Figure 31
Whena = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.15, multistability results of the system where the blue points denote 8-fold periodic states corresponding to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).From the above results, we can know that, under different initial value conditions, three kinds of attractors coexist, so the system is multistable.In this paper, to facilitate the simulation of the system, we discretize equation (1) just as what literature [33] did. Compared with the corresponding continuous-time system, the discrete system has the following different points: (i) the system states are discrete and sequence output is discrete and (ii) the curves in the phase diagrams are formed by discrete points.


## 3.1.α Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters α whose results are shown in Figures 1 and 2. From Figure 1, the following conclusions can be obtained: when α < 0.8, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the chaos is generated in the system; when 0.8 < α < 1, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the periodic window is generated in the system; when α > 1, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the periodic window turns into chaos in the system. In short, as α increases, the system experiences chaos-period-chaos. Figure 2 is the flip bifurcation diagram versus α. According to Figure 2, when α is used as 0.05, 0.95, and 2, we can obtain sequence diagrams shown in Figure 3–5, respectively, and phase diagrams shown in Figures 6–9, respectively.Figure 1
Lyapunov spectrum versusα when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5.Figure 2
Flip bifurcation diagram versus parametersα when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, β = 0.2, γ = 0.5.Figure 3
x–n sequence diagram which denotes that the system is chaotic when α = 0.05.Figure 4
x–n sequence diagram which denotes that the system is periodic when α = 0.95.Figure 5
x–n sequence diagram which denotes that the system is chaotic when α = 2.Figure 6
x–y phase diagram which denotes that the system is chaotic when α = 0.05.Figure 7
y-z phase diagram which denotes that the system is chaotic when α = 0.05.Figure 8
x-y phase diagram which denotes that the system is periodic when α = 0.95.Figure 9
x-y phase diagram which denotes that the system is chaotic when α = 2.


## 3.2.β Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters β whose results are shown in Figures 10 and 11. From Figure 10, the following conclusions can be obtained: when β < 1.12, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the chaos is generated in the system; when 1.12 < β < 1.24, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the periodic window is generated in the system; when 1.24 < β < 3.5, we can know LE1 > 0, LE2 < 0, and LE3 < 0, and the periodic window turns into chaos in the system; when β > 3.5, we can know LE1 < 0, LE2 < 0, and LE3 < 0, and the chaos turns into a periodic window in the system.Figure 10
Lyapunov spectrum versusβ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5.Figure 11
Flip bifurcation diagram versus parametersβ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, γ = 0.5.In short, asβ increases, the system experiences chaos-period-chaos-period. Figure 11 is flip bifurcation diagram versus β. According to Figure 11, when β is used as 0.7, 1.2, and 3, we can obtain sequence diagrams shown in Figures 12–14, respectively, and phase diagrams shown in Figures 15–18, respectively.Figure 12
x–n sequence diagram which denotes that the system is chaotic when β = 0.7.Figure 13
x–n sequence diagram which denotes that the system is periodic when β = 1.2.Figure 14
x–n sequence diagram which denotes that the system is chaotic when β = 3.Figure 15
x-y phase diagram which denotes that the system is chaotic when β = 0.7.Figure 16
x–y phase diagram which denotes that the system is periodic when β = 1.2.Figure 17
x–y phase diagram which denotes that the system is chaotic when β = 3.Figure 18
y-z phase diagram which denotes that the system is chaotic when β = 3.


## 3.3.γ Is Changed and Other Parameters of System (18) Are Fixed

We usea = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2 and simulate the Lyapunov spectrum and bifurcation diagram versus parameters β whose results are shown in Figures 19 and 20. From Figure 19, when γ = 0.096 and γ = 0.184, we can know LE1 = 0, LE2 < 0, and LE3 < 0. From Figure 20, it is clear that when 0 < γ < 0.096, the system has a 4-fold periodic; when 0.096 < γ < 0.2, the system has a 8-fold periodic. From Figure 19, when 0.2 < γ < 0.376, we know that LE1 > 0, LE2 < 0, and LE3 < 0, and the system is chaotic; when γ = 0.376 and γ = 0.392, we know that LE1 = 0, LE2 < 0, and LE3 < 0. Hence, from Figure 20, we know that when 0.376 < γ < 0.39, the system is periodic. When γ > 0.392, we know that LE1 > 0, LE2 < 0, and LE3 < 0, and the system is chaotic. In short, as γ increases, the system experiences a route of 4-fold periodic, 8-fold periodic, and chaos-period-chaos. Figure 20 is the flip bifurcation diagram versus β. According to Figure 20, when γ is used as 0.08, 0.15, 0.39, and 0.8, we can obtain sequence diagrams shown in Figure 21–24, respectively, and phase diagrams shown in Figure 25–29, respectively.Figure 19
Lyapunov spectrum versusγ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2.Figure 20
Flip bifurcation diagram versus parametersγ when a = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2.Figure 21
x-n sequence diagram which denotes that the system is 4-fold periodic when γ = 0.08.Figure 22
x-n sequence diagram which denotes that the system is 8-fold periodic when γ = 0.15.Figure 23
x-n sequence diagram which denotes that the system is periodic when γ = 0.39.Figure 24
x-n sequence diagram which denotes that the system is chaotic when γ = 0.8.Figure 25
x-y phase diagram which denotes that the system is 4-fold periodic when γ = 0.08.Figure 26
x-y phase diagram which denotes that the system is 8-fold periodic when γ = 0.15.Figure 27
x-y phase diagram which denotes that the system is periodic when γ = 0.39.Figure 28
x-y phase diagram which denotes that the system is chaotic when γ = 0.8.Figure 29
y-z phase diagram which denotes that the system is chaotic when γ = 0.8.


## 3.4. Multistability of the System

Multistability means that when parameters are fixed and the initial values are different, the system presents different attractors. When usinga = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.08, we simulate the multistability of the system, and simulation results are shown in Figure 30 where the blue points denote 4-fold periodic states corresponding to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).Figure 30
Whena = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, and β = 2, γ = 0.08, multistability results of the system where the blue points denote 4-fold periodic states correspond to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).When usinga = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.15, we simulate the multistability of the system, and simulation results are shown in Figure 31 where the blue points denote 8-fold periodic states corresponding to initial value (0.1, 0.1, 0.1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).Figure 31
Whena = 1, b = 1.14, c = 2.93, d = 1.3, e = 0.21, f = 0.1, g = 0.35, h = 2, α = 0.05, β = 2, and γ = 0.15, multistability results of the system where the blue points denote 8-fold periodic states corresponding to initial value (0.1, 0.1, 0 .1, 0.1, 0.1), the black point denotes a stable state of equilibrium point corresponding to initial value (0.1, 0.1, 0, 0.1, 0.1), and the rose points denote chaotic state corresponding to initial value (0.1, 0, 0.1, 0.1, 0.1).From the above results, we can know that, under different initial value conditions, three kinds of attractors coexist, so the system is multistable.In this paper, to facilitate the simulation of the system, we discretize equation (1) just as what literature [33] did. Compared with the corresponding continuous-time system, the discrete system has the following different points: (i) the system states are discrete and sequence output is discrete and (ii) the curves in the phase diagrams are formed by discrete points.


## 4. Comparison with Literature [33]

In [33], it is relatively unspecific for considering all real estate companies as a prey team, which does not fit the specific circumstance of real estate companies. Our model considers the large, medium, and small real estate enterprises for three different prey teams. Compared with the 3D model in literature [33], our model is more accordant with the specific circumstance of real estate companies and emerges richer dynamic characteristics. In [33], the range of the positive Lyapunov exponent is narrower, and the corresponding chaos range is also narrower. On the contrary, in our paper, the range of the positive Lyapunov exponent is wider, and the corresponding chaos range is also wider. In addition, in our model, there are “period-chaos-period” or “chaos-period-chaos” bifurcation characteristics.


## 5. Conclusion

In this paper, a 5D predator-prey evolutionary model of real estate enterprises was proposed. This model considers the large, medium, and small real estate enterprises for three different prey teams. We made a nonlinear dynamic analysis on a 5D predator-prey evolutionary system in the real estate market and simulated the corresponding discrete system. The simulation results about the Lyapunov spectrum, bifurcation diagram, sequence diagram, and phase diagram were given. In this paper, the range of the positive Lyapunov exponent is wide, and the corresponding chaos range is also wide. The system has rich dynamic characteristics showing “period-chaos-period” or “chaos-period-chaos” bifurcation performance.



---

*Source: 6670222-2021-01-11.xml*