# Stability, Multistability, and Complexity Analysis in a Dynamic Duopoly Game with Exponential Demand Function

Authors: Hui Li; Wei Zhou; Tong Chu
Journal: Complexity
 (2021)

---

## Abstract

In this paper, a discrete-time dynamic duopoly model, with nonlinear demand and cost functions, is established. The properties of existence and local stability of equilibrium points have been verified and analyzed. The stability conditions are also given with the help of the Jury criterion. With changing of the values of parameters, the system shows some new and interesting phenomena in terms to stability and multistability, such as V-shaped stable structures (also called Isoperiodic Stable Structures) and different shape basins of attraction of coexisting attractors. The eye-shaped structures appear where the period-doubling and period-halving bifurcations occur. Finally, by utilizing critical curves, the changes in the topological structure of basin of attraction and the reason of “holes” formation are analyzed. As a result, the generation of global bifurcation, such as contact bifurcation or final bifurcation, is usually accompanied by the contact of critical curves and boundary.


---

## Body

## 1. Introduction

Nonlinear dynamical system can describe many complicated and curious phenomena, and the theory of nonlinear dynamics is widely applied in many fields of scientific research studies. As an important branch of modern mathematics, nonlinear theory makes up the defect of the linear function and can be well used to characterize the regularities of ecosystem [1], economical model [2, 3], and so on. A function in the exponential form [4] has a good property of nonnegativity, which can make sure that the solutions of the exponential function are all greater than zero. According to this practical economic significance, it is necessary to adopt this nonlinear type of demand function with exponential. In other words, the exponential demand function guarantees that the price of the product is always positive.For the process of dynamical system evolution, most scholars usually pay attention to or are interested in whether the dynamical system will keep in the steady state for a long time rather than experience a short transient state. Hence, the game evolution plays an important role in describing decision-making process, showing the collision of strategic sparks and the wisdom of strategies. In a dynamical system, nonlinear dynamics theory can clearly explain evolution regularities and complex dynamic behaviors. Especially in the economic field [5, 6], with the deepening of the research level, many scholars have gradually taken many practical factors into account in modeling process and combining nonlinear dynamics with game theory. In the real market, the cost function may not be monotonous or has a linear relationship with the total supply, that is to say, the nonlinear cost form is closer to the reality. Brianzoni et al. [7] changed the form of the cost function from linear to nonlinear, and they analyzed dynamic characteristics of a Bertrand competition system. Pu and Ma [8] investigated a dynamic game model with four oligopolistic players, in which a cubic total-cost demand function was considered. It is worth noticing that a group of scholars, Peng et al. [9] considered a quadratic form of the cost function to establish a remanufacturing duopoly model and to study its complex dynamic characteristics. In addition, an economic duopoly dynamical system with taking the nonlinear cost function into consideration has been built in this paper, and a lot of very peculiar and complex dynamic characteristics will arise, such as bifurcations, chaos, and coexisting attractors, which we will discuss in later sections.In real market competitions, the more market information the firm has, the more rational it is. However, in order to completely master whole market information, the firms must need to pay an extremely high amount of information cost. It is rational that no one wants to pay such high information cost. As a result, most firms in real market do grasp some market information but not all, which leads the firms to form “myopic” in some respects, so the firms will constantly adjust their decision in order to maximize their own interests. Such a decision-making process is also known as “trial and error.” The classical adjustment rule, i.e., “rule of thumb” (which can be also called as the gradient adjustment mechanism), does not require the firms to have a high degree of rationalities and has been widely used in many investigations; for instance, Zhou et al. [10], Guo and Ma [11], Elsadany [12], and Cao et al. [13]. Furthermore, the boundedly rational firms adjust and update their competitive strategies in the next period by utilizing local estimation of their marginal profits in the previous period. Therefore, the firms, who implement such local adjustment rule, do not need to have whole knowledge on demand and cost of the market.At the end of the last century, Bischi et al. [14] built a dynamic Cournot duopoly model and found that an interesting phenomenon of synchronization may occur, where the global dynamics had been analyzed by critical curves. At the beginning of this century, Bischi et al. [15, 16], Agliari et al. [17], Agiza and Elsadany [18], and other scholars had done much work on studying global bifurcation of a dynamical system. Their remarkable achievements laid the foundation for the analysis of complex dynamic phenomena by using noninvertible mapping and critical curves. Global bifurcations usually happen when the contacts occur between critical curves and boundaries of basins of attraction or between borders of attractors and boundaries of basins of attraction. After global bifurcation, many interesting and complex phenomena will appear. Generally speaking, complex dynamic phenomena usually include multistability, such as coexistence of attractors between periodic points or between chaotic attractors and periodic points, the occurrence of “holes” [15, 16], and so on. Due to these complex phenomena, the structure and stability of basin of attraction may dramatically change. Zhou and Wang [19] analyzed the properties in terms of bifurcation and multistability in a two-stage duopoly game with differentiated product and R&D spillovers. Cavalli and Naimzada [20] also had described complex dynamic behaviors and the phenomenon of multistability. The critical curve is one of the noticeable features for the noninvertible map, and we can use critical curves to analyze and investigate complex dynamic behaviors appearing in the dynamical system.In addition, the method of critical curves (see in Ref. [21] and references therein) is an efficient approach to further investigate global dynamics. Bischi and Lamantia [22] studied complex structures of basins of attraction and the phenomenon of coexisting attractors in a discrete dynamic economy game. Gao [23] also had carried out some works on complex properties about a noninvertible iteration map. Besides this, the paper also used center manifold and bifurcation to give the conditions when classical bifurcations arise, that is, pitch-fork bifurcation, flip bifurcation, and Neimark–Sacker bifurcation. Most recently, the scholars Bischi and Baiardi [24, 25], Cao et al. [13], and Zhang et al. [26] all utilized critical curves to study the complex dynamics of the built models.The primary purpose of this paper is to establish a discrete-time dynamic duopoly game model with bounded rationality, where both demand and cost functions are considered in nonlinear forms. Furthermore, the dynamic properties of the built system can be understood more clearly by numerical simulation, so we mainly use the single-parameter (1D) bifurcation diagram, the double-parameter (2D) bifurcation diagram, the largest Lyapunov exponent diagram, and the basin of attraction. By changing the values of the parameters, there are much complex phenomena that may arise, and we can further study their influences on the stability of the system or the structure of the basin of attraction through utilizing the theory of critical curves in subsequent sections.The rest parts of this paper are organized as follows. In Section2, a dynamic duopoly Cournot game under the exponential demand function and quadratic cost function is established, where these two firms are both considered as having bounded rationality. The existence of equilibrium points has been verified. The locally stable conditions of equilibrium points have been discussed in Section 3. In Section 4, we have analyzed the distinct routes of the system to chaotic states and found that there are some interesting and curious phenomena of stability and multistability, including the basin of attraction, coexisting attractors, and so on, by numerical simulation. Furthermore, we also use critical curves to study the steady state of the system and the change process of the structure of basin of attraction. Finally, some conclusions have been drawn in Section 5.


## 2. Modeling and the Existence of Equilibria

It is supposed that there are two firmsi (labeled by i=1and2) in the economic market, producing homogeneous products, and both these two kinds of homogeneous products are supplied to a same market. The prices of products may have some relations with the total quantity or the production costs or some other elements. In this paper, it is assumed that the price of firm ii=1and2 is an exponential inverse demand function defined by the outputs qi≥0 of two firms, that is, PQ=ae−q1+q2, in which a>0 is the maximum selling price of the products in the market and Q=q1+q2 is the number of total quantity of these two kinds of homogeneous products. The exponential inverse demand function has good nonnegativity, so we can be sure that the solutions of the exponential function are all greater than zero, even tending to infinity, or even tending quite near to zero. In addition, in order to be closer to reality, we consider that the cost function of the firm i has a quadratic form, that is, Ciqi=ciqi2, where ci>0 is the marginal cost, and the relation between a and ci satisfies a>ci. Hence, the profits of these two firms can be, respectively, written as follows:(1)∏1=ae−q1+q2q1−c1q12,∏2=ae−q1+q2q2−c2q22.Let∏i, respectively, deduce the derivative with regard to the output qi, i=1and2. Then, we can get the corresponding marginal profit functions as follows:(2)∂∏1∂q1=a1−q1e−q1+q2−2c1q1,∂∏2∂q2=a1−q2e−q1+q2−2c2q2.In the fierce market competition, for the sake of becoming completely rational and having whole information of the market, the firms must need to pay great number of information costs, whereas most firms are not willing to do this activity. In fact, in addition to that, it is really difficult to master some aspects of market information, let alone all market information. Hence, in this situation, the two firms we considered before can be seen as having incomplete market information. In other words, the firmii=1and2 is boundedly rational. The firms adjust their production strategies in next period t+1 in terms of their marginal profit in the previous period t. It means that if the marginal profit in t period is greater than zero (or less than zero), then the firm i will increase (or decrease) its outputs qi. If the marginal profit in t period is equal to zero, then the firm i will keep its volume of production unchanged. The “gradient adjustment mechanism” is a useful method to simulate how the firms adjust their decisions so as to make sure to get maximum interests. Thus, we introduce this classical adjustment rule into our model, and then, the dynamic adjustment mechanism of the firm i with respect to its outputs can be described in the following form:(3)qit+1=qit+viqit∂∏it∂qit,t∈Z,in which vi>0 represents the adjustment speed of the firm i. Through substituting equation (2) into equation (3), the discrete-time dynamic adjustment mechanism of the two firms can be specifically expressed as follows:(4)q1t+1=q1t+v1q1ta1−q1teq1t+q2t−2c1q1t,q2t+1=q2t+v2q2ta1−q2teq1t+q2t−2c2q2t.For simplicity, the auxiliary variables have been introduced in system (4), that is, αi=civi and βi=avii=1and2, and we have that αiandβi∈0,+∞. We introduce an evolutionary operator T and auxiliary variables as αi=civi and βi=avi, i=1and2, where αiandβi∈0,+∞. Hence, a discrete-time dynamic duopoly Cournot game can be further simplified as the following:(5)T:q1t+1=q1t+β1q1t1−q1te−q1t+q2t−2α1q12t,q2t+1=q2t+β2q2t1−q2te−q1t+q2t−2α2q22t,in which T is an evolutionary operator.According to the knowledge of nonlinear dynamics, letqit+1=qit, i=1and2, so we can know that the equilibrium point in system (5) must satisfy the conditions listed as follows:(i)
q1=0 or(ii)
β11−q1e−q1+q2=2α1q1(iii)
q2=0 or(iv)
β21−q2e−q1+q2=2α2q2On the basis of the above conditions, we can deduce that there are four equilibrium points existing in system (5). For instance, when condition (i) and condition (iii) simultaneously hold, the equilibrium point E1=0,0 of system (5) can be obtained. For the rest of the boundary equilibrium points, we can draw the following conclusions:(1)
The boundary equilibrium pointE2=x∗,0 can be obtained by condition (ii) and condition (iii), where x∗ satisfies the equation β11−x∗e−x∗=2α1x∗. If we rewrite this relationship as a function with regard to a variable z, that is, gz=2α1z−β11−z/ez, then we have g0=−β1<0 and g1=2α1>0. So we can know that there is at least one point x∗∈0,1 that satisfies β11−x∗e−x∗=2α1x∗. Finally, the existence of E2=x∗,0 is verified.(2)
Analogously, the boundary equilibrium pointE3=0,y∗ can be obtained by condition (i) and condition (iv), where y∗ satisfies the equation β21−y∗e−y∗=2α2y∗. If we rewrite this equation as a function with respect to a variable z, that is, wz=2α2z−β21−z/ez, then we will get two inequalities that w0=−β2<0 and w1=2α2>0. So we can know that there is at least one point y∗∈0,1 that satisfies β21−y∗e−y∗=2α2y∗. Finally, the existence of E3=0,y∗ is also verified.Hence, we do verify the existence of these three boundary equilibrium points. Then, the analysis of the Nash equilibrium point and its existence in system (5) is given by the proposition as the following.Proposition 1.
The unique Nash equilibrium pointE4=m∗,n∗ can be obtained through combining of condition (ii) with condition (iv), and m∗ satisfies the equation 2α1m∗=β11−m∗e−σm∗−σm∗2+m∗2+m∗/σ−σm∗+m∗, and n∗ satisfies the equation 2α2n∗=β21−n∗e−σn∗+σn∗2−n∗2+n∗/1+σn∗−n∗, where σ=α2β1/α1β2.Proof.
First, the equationem∗+n∗=β11−m∗/2α1m∗=β21−n∗/2α2n∗ can be obtained from condition (ii) and condition (iv), where n∗=m∗/σ−σm∗+m∗ and m∗=σn∗/1+σn∗−n∗. Then, let h1u=2α1u−β11−ue−σu−σu2+u2+u/σ−σu+u, we can easily get two inequalities that h10=−β1<0 and h11=2α1>0. Hence, we can know that there is at least one point m∗ in 0,1 that satisfies the equation 2α1m∗=β11−m∗e−σm∗−σm∗2+m∗2+m∗/σ−σm∗+m∗. Analogously, let h2u=2α2u−β21−ue−σu+σu2−u2+u/1+σu−u, we have two inequalities that h20=−β2<0 and h21=2α2>0. So there is at least one point n∗∈0,1 that satisfies the equation 2α2n∗=β21−n∗e−σn∗+σn∗2−n∗2+n∗/1+σn∗−n∗, where σ=α2β1/α1β2. Finally, the existence of the Nash equilibrium point E4=m∗,n∗ is verified.
In a real economic system, the equilibrium pointE1=0,0 represents the situation that these two firms both go bankrupt or get out of the market. The boundary equilibrium points E2=x∗,0 and E3=0,y∗ can be regarded as the situation that one of these two firms gets out of the market, while another one becomes a monopoly.


## 3. Local Stability Analyses of Equilibria

From the knowledge of nonlinear dynamics, for the sake of judging the local stability of equilibrium points, we only need to judge the relationship between 1 and the eigenvalueλii=1and2 of the Jacobian matrix. That is to say, if λi>1i=1and2, then the equilibrium point is an unstable node. If λi<1i=1and2, then the equilibrium point is a stable node. While if one of λii=1,2 is great than 1 and the other is smaller than 1, then the equilibrium point is a saddle point. The Jacobian matrix corresponding to system (5) can be expressed as follows:(6)J=1−4α1q1+β11−3q1+q12e−q1+q2β1q1q1−1e−q1+q2β2q2q2−1e−q1+q21−4α2q2+β21−3q2+q22e−q1+q2.The local stability of boundary equilibrium points of system (5) has been discussed, and the following propositions can be drawn.Proposition 2.
The boundary equilibrium pointE1=0,0 is an unstable node.Proof.
In system (5), the Jacobian matrix corresponding to the point E1=0,0 can be written as(7)JE1=1+β1001+β2.
It is clear that Jacobian matrix (7) is a diagonal matrix, and the corresponding eigenvalues are λ1=1+β1 and λ2=1+β2, respectively. Since the nonnegativity of β1 and β2, we have that λ1>1 and λ2>1. Therefore, the equilibrium point E1=0,0 is an unstable node.Proposition 3.
The boundary equilibrium pointE2=x∗,0 is unstable.(1)
Whenβ1>−2ex∗/x∗2−x∗−1, then the boundary equilibrium point E2=x∗,0 is a saddle point(2)
Whenβ1<−2ex∗/x∗2−x∗−1, then the boundary equilibrium point E2=x∗,0 is an unstable focusProof.
Apparently, the Jacobian matrix at equilibrium pointE2=x∗,0 is an upper triangular matrix, that is,(8)JE2=1+β1e−x∗x∗2−x∗−1β1x∗x∗−1e−x∗01+β2e−x∗.
The eigenvalues corresponding to this Jacobian matrix areλ1=1+β1e−x∗x∗2−x∗−1 and λ2=1+β2e−x∗, respectively. Since x∗∈0,1, we have that x∗2−x∗−1<0 and β1>0. Hence, λ1<1 and λ2>1. However, the relationship between λ1 and 1 needs to be discussed by using classification as follows.(1)
Whenβ1>−2ex∗/x∗2−x∗−1, then we can deduce that λ1>−1. Since λ1<1, then we can know that the moduli of λ1 and λ2 satisfy λ1<1 and λ2>1. As a result, E2=x∗,0 is a saddle point.(2)
Conversely, whenβ1<−2ex∗/x∗2−x∗−1, then we have a conclusion that λ1<−1. Furthermore, combining the conclusions λ1>1 with λ2>1, we can conclude that E2=x∗,0 is an unstable focus.Proposition 4.
The boundary equilibrium pointE3=0,y∗ is unstable.(1)
Whenβ2>−2ey∗/y∗2−y∗−1, the boundary equilibrium point E3=0,y∗ is a saddle point(2)
Whenβ2<−2ey∗/y∗2−y∗−1, the boundary equilibrium point E3=0,y∗ is an unstable focusProof.
Obviously, the Jacobian matrix atE3=0,y∗ is a lower triangular matrix as the following:(9)JE3=1+β1e−y∗0β2y∗y∗−1e−y∗1+β2e−y∗y∗2−y∗−1.
The eigenvalues corresponding to this Jacobian matrix areλ1=1+β1e−y∗ and λ2=1+β2e−y∗y∗2−y∗−1, respectively. This is a similar situation like in Proposition 3, so we do same steps as in Proposition 3. Since y∗ belongs to interval 0,1, then we can get two inequalities y∗2−y∗−1<0 and β2>0, and then, we already know that λ1>1 and λ2<1. However, the relation between λ2 and 1 needs to be discussed by using classification as follows.(1)
Whenβ2>−2ey∗/y∗2−y∗−1, then we can deduce that λ2>−1. Then, we get the moduli of λ1 and λ2 satisfy λ1>1 and λ2<1. As a result, E3=0,y∗ is a saddle point.(2)
Whenβ2<−2ey∗/y∗2−y∗−1, then we can deduce that λ2<−1. Then, we get the moduli of λ1 and λ2 satisfy λ1>1 and λ2>1. As a result, E3=0,y∗ is an unstable focus.
As for the unique Nash equilibrium of system (5), the corresponding Jacobian matrix at E4=m∗,n∗ can be expressed as(10)JE4=1+β1e−m∗+n∗m∗2−m∗−1β1e−m∗+n∗m∗m∗−1β2e−m∗+n∗n∗n∗−11+β2e−m∗+n∗n∗2−n∗−1,in which m∗ satisfies the equation 2α1m∗=β11−m∗e−m∗+n∗ and n∗ satisfies the equation 2α2n∗=β21−n∗e−m∗+n∗. Therefore, the trace and the determinant of the Jacobian matrix JE4 are denoted as Tr and Det, respectively. The specific expressions corresponding to Tr and Det can be, respectively, written as follows:(11)Tr=2+β1e−m∗+n∗m∗2−m∗−1+β2e−m∗+n∗n∗2−n∗−1,Det=1+e−m∗+n∗β1m∗2−m∗−1+β2n∗2−n∗−1+β1β2e−2m∗+n∗1−m∗2−n∗2+m∗+n∗.
As we all know, the Nash equilibrium pointE4 is stable when the moduli of its eigenvalues corresponding to the Jacobian matrix JE4 are all less than 1, while the forms of eigenvalues under the context of the exponential term are quite sophisticated. Therefore, we need to study the local stability around the Nash equilibrium E4 with the help of the Jury criterion.Proposition 5.
The Nash equilibrium pointE4=m∗,n∗ is asymptotically stable if and only if one of the following six conditions holds:(1)
When2eA/−C<eA−B/D, the values of the parameters β1 and β2 satisfy 0<β2<2eA/−C and β1<4e2A+2β2eAC/2eA−B+β2−D(2)
When2eA/−C<eA−B/D, the values of the parameters β1 and β2 satisfy 2eA/−C<β2<eA−B/D and β1>0(3)
When2eA/−C<eA−B/D, the values of the parameters β1 and β2 satisfy eA−B/D<β2<2eA−B/D and 0<β1<β2eA−C/eAB+β2D(4)
When2eA/−C<eA−B/D, the values of the parameters β1 and β2 satisfy β2>2eA−B/D and β1>4e2A+2β2eAC/2eA−B+β2D(5)
When2eA/−C>eA−B/D, the values of the parameters β1 and β2 satisfy 0<β2<eA−B/D and β1>0 or 2eA/−C<β2<2eA−B/D and β1>0(6)
When2eA/−C>eA−B/D, the values of the parameters β1 and β2 satisfy eA−B/D<β2<2eA/−C and β1>β2eA−C/eAB+β2D, or β2>2eA−B/D and β1>β2eA−C/eAB+β2D
whereA=m∗+n∗, B=m∗2−m∗−1, C=n∗2−n∗−1, and D=1−m∗2−n∗2+m∗+n∗.Proof.
According to the Jury criterion, we can verify the local stability ofE4. In other words, if all of conditions of the Jury criterion hold, then we can know that the Nash equilibrium E4 is asymptotically stable. The specific expression of the Jury criterion is given as follows:(12)i1−Tr+Det>0,ii1+Tr+Det>0,iii1−Det>0.
Substituting the specific formulations ofTr and Det into equation (12), then we will get equivalent expressions as follows:(13)iβ1β2e−2m∗+n∗1−m∗2−n∗2+m∗+n∗>0,ii4+2β1e−m∗+n∗m∗2−m∗−1+2β2e−m∗+n∗n∗2−n∗−1+β1β2e−2m∗+n∗1−m∗2−n∗2+m∗+n∗>0,iii−β1e−m∗+n∗m∗2−m∗−1−β2e−m∗+n∗n∗2−n∗−1−β1β2e−2m∗+n∗1−m∗2−n∗2+m∗+n∗>0.
For simplicity, letA=m∗+n∗, B=m∗2−m∗−1, C=n∗2−n∗−1, and D=1−m∗2−n∗2+m∗+n∗. Since the conditions m∗∈0,1 and n∗∈0,1, then we can know that A>0, B<0, C<0, and D>0. Then, the three conditions in equation (13) are equivalent to(14)i′β1β2D>0,ii′4e2A+2eAβ1B+2eAβ2C+β1β2D>0,iii′−β1eAB−β2eAC−β1β2D>0.
According toβi>0 and D>0, it can be judged that the condition i′ in equation (14) (that is, the first condition in the Jury criterion) always holds. Therefore, system (5) will not cross the eigenvalue 1 of the characteristic equation, namely, system (5) will not experience a transcritical bifurcation. Next, we just need to judge whether the last two conditions in the Jury criterion hold, and then, we can define local stability conditions of the Nash equilibrium point E4. Obviously, the condition ii′ is equivalent to the form of β12eAB+β2D+4e2A+2eAβ2C>0. In addition, in order to further analyze the local stability of E4 and get the value ranges of the parameters β1 and β2, the inequality β12eAB+β2D+4e2A+2eAβ2C>0 can be divided into two parts for discussion. One part is equivalent to β12eAB+β2D=0, that is, β2=2eA−B/D. Another part is equivalent to 4e2A+2eAβ2C=0, that is, β2=2eA/−C. Additionally, since 2eA/−C−2eA−B/D<0, that is, 2eA/−C<2eA−B/D, the auxiliary parameter βi involves the adjustment speed of the firm i, and its value range can be divided into three cases as the following.(i′)
When0<β2<2eA/−C, β1<4e2A+2β2eAC/2eA−B+β2−D.(ii′)
Whenβ2>2eA−B/D, β1>4e2A+2β2eAC/2eA−B+β2D.(iii′)
When2eA/−C<β2<2eA−B/D, β1>4e2A+2β2eAC/2eA−B+β2D, while 4e2A+2β2eAC/2eA−B+β2D<0, since the nonnegativity of the parameter β1, β1>0 in this case.The above are the value ranges of the auxiliary parameterβi that need to be satisfied, when the equivalent condition ii′ in the Jury criterion holds.The situation is discussed as following, when the equivalent conditioniii′ in the Jury criterion holds. Analogously, we will get an inequality β2eAC+β1eAB+β2D>0, and then, it will be discussed by divided into two cases, that is, β2eAC=0 and β1eAB+β2D=0. In other words,(15)β2eAC=0⇔β2=0,β1eAB+β2D=0⇔β2=eA−BD.The discussion in terms ofβ1 has two situations.(iv′)
whenβ2>eA−B/D, β1<β2eA−C/eAB+β2D.(v′)
when0<β2<2eA−B/D, β1<β2eA−C/eAB+β2D; we already know that β2eA−C/eAB+β2D<0, and since the nonnegativity of the parameter β1, so we consider that β1>0 in this situation.
In order to make the Jury criterion hold, we only need to take the intersection of two of the abovementioned five conditions, that is,i′, ii′, iii′, iv′, and v′. That is to say, if one of the above six conditions in Proposition 4 holds, then E4=m∗,n∗ is locally asymptotical stable.


## 4. Numerical Simulation

In the nonlinear dynamical system, the phase diagram or bifurcation diagram may include a lot of information in terms of evolution process of the dynamical system. Furthermore, due to the exponential term that appears in the equation, it is very difficult for us to get an analytical solution. Therefore, it is of great necessity to exploit the method of numerical simulation to intuitively and deeply explore and display complex dynamics of a discrete-time dynamic duopoly Cournot game proposed above. The main investigation tools that we have used are 1D bifurcation diagram, 2D bifurcation diagram, the largest Lyapunov exponent diagram, the basin of attraction, etc. Hence, in this section, the dynamic behaviors of system (5) are investigated through those numerical methods mentioned above.


### 4.1. Complex Dynamic Behaviors

First, a set of parameters is fixed asa=11.1542, c1=1.8044, and c2=1.3949, and we choose adjustment speeds v1 and v2 as the bifurcation parameters. As shown in Figure 1, the color bar is on the right side of figures, in which different colors denote different numbers of the period. The steel blue represents the stable area of the Nash equilibrium point, and the orange represents that the points in this area are at a period-2 state, and analogously, the grass green represents the period-4. The divergent area is denoted as the color white, which is also marked as “ET” in the color bar. However, the meaning of the dark black area is more complex, where the dynamic states included in this area, such as the quasi-period state and chaotic state. From Figure 1(a), abundant dynamic phenomena can be observed. As the values of v1 and v2 gradually increase, the bifurcation route of system (5) will start from the period-1 steel blue stable area, then pass through the period-2 orange area and period-4 grass green region, and finally enter into the dark black region, which means that system (5) will change from a stable state to an unstable state. Corresponding to the real economic market, this phenomenon means that if two firms blindly improve their adjustment speeds without strategy, the market will become unstable and eventually may enter into chaos. Furthermore, Figure 1(b) is a local enlargement of Figure 1(a). The structures of shrimp-like are called Isoperiodic Stable Structures (the abbreviation is ISSs), seen in Ref. [27], which can be clearly shown in Figure 1(b). The location where the shrimp-like structures are at is denoted by a white straight line, and the expression of this line is v2=−5.6v1+1.904.Figure 1
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. (a) The double-parameters’ bifurcation diagram of system (5). (b) The local enlargement of Figure 1(a).
(a)(b)In Figure2, the red curve denotes the outputs q1 of firm 1, and the black curve denotes the outputs q2 of firm 2. First, choose the adjustment speed v1 as the bifurcating parameter, and the other parameters are fixed as a=11.1542, c1=1.8044, c2=1.3949, and v2=0.1. Under this group of parameters, system (5) will go through a flip bifurcation at first and finally drop into a chaotic state, as shown in Figure 2(a). The bifurcation diagram of q2 with respect to v1, marked by a black box, is displayed in the enlargement. From Figure 2(a), when the adjustment speed of firm 1 satisfies v1<0.4197, system (5) is at a period-1 stable state. In addition, v1=0.4197 is a bifurcation point, which means that the stability of system (5) at this point will completely change, and then, a period-doubling bifurcation subsequently occurs. Once the value of v1 is greater than 0.5371, system (5) will uncontrollably reach a disordered chaotic state.Figure 2
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. The red curve denotes q1 of firm 1, and the black curve denotes q2 of firm 2. (a) The 1D bifurcation diagram when v2=0.1. (b) Corresponding to Figure 1(b), the 1D bifurcation diagram on the white straight line is displayed with its largest Lyapunov exponent diagram. Fixed the parameters as a=11.1542, v1=0.3625, and v2=0.4517. (c) The 1D bifurcation diagram when c2=1.8. (d) The 1D bifurcation diagram when c1=1.8.
(a)(b)(c)(d)Under same parameter conditions, it can be seen from Figure2(b) that the dynamic behaviors of system (5) are very complicated when shrimp-like structures appear. At this time, the blue curve represents the largest Lyapunov exponent curve corresponding to the red bifurcation curve. It is well known that we can easily determine the state of a dynamical system through calculating its Lyapunov exponents, and we usually need to calculate the largest one. Therefore, the dynamical system corresponds to a stable state, if the largest Lyapunov exponent is less than zero. Instead, if the largest Lyapunov exponent is greater than zero, it means that the dynamical system is at a chaotic state. In addition, it is worth noting that the largest Lyapunov exponents fluctuate up and down around zero, and then, the dynamical system is at a quasi-period state. There are many phenomena of “periodic window” that exist, and system (5) constantly transforms between the chaotic state and the periodic stable state. With the parameter v1 further increasing, system (5) will enter a stable state from a chaotic state for a short time and finally return to the chaotic state.When the parameters are changed toa=11.1542, v1=0.3625, and v2=0.4517, then we can get two bifurcation diagrams with respect to costs of two firms, respectively. As shown in Figure 2(c), where c2=1.8, when c1<0.7755, system (5) is at the period-1 stable state. At c1=0.7755, a typical flip bifurcation (also called period-doubling bifurcation) occurs. In addition, it is worth to notice that an obvious phenomenon called “jump” appears near c1=1.8723, which means that there may be a phenomenon of coexistence. Similarly, we choose c1=1.8 in Figure 2(d),then we can know that c2=0.6330 is a bifurcation point, and once system (5) crosses this point, then a flip bifurcation will occur. The dynamic phenomena between 1.3500<c2<1.4180 are very complicated. There is also a “jump” phenomenon that exists near c2=1.6083. With constant increase of costs of two firms, dynamical system (5) will eventually enter into chaos. However, through comparing Figures 2(c) and 2(d), we can draw a result that, under the premise of same costs, the speed of the system entering into chaos with c2 is higher than that with c1.Another group of parameters are fixed asa=11.6054, c1=2.9154, and c2=2.9639, and we choose adjustment speeds v1 and v2 as the bifurcating parameters. A totally different 2D bifurcation diagram is shown in Figure 3(a), in which many colorful and symmetrical structures exist. Similarly, the dark black area includes complex phenomena such as chaos and quasiperiodic states, and the white area denotes the divergent area. It is worth noticing that the dynamic behaviors, marked in a red box, have significances to investigate. Hence, we enlarge it in Figure 3(b) and increase the number of periods to 200. As shown in Figure 3(b), abundant and complex dynamic phenomena are displayed. There are many fractal structures that can be seen, and the eye-shaped structure is symmetrically arranged on both sides of the ribbon. The dynamical behaviors, located in a white box of Figure 3(b), are enlarged in Figure 3(c). Obviously, there are apparent V-shaped islands called ISSs structures, which are marked as a, b, and c in the yellow box. It is known that this kind of structures are Lyapunov stable structures or Lyapunov stable islands, that is, if parameters are chosen in these structures, the state of the dynamical system is regular and stable, seen in Ref. [27]. Since the 2D Lyapunov exponent diagram is a useful tool to clearly distinguish stable state, quasi-periodic state, chaotic state, and divergent state, we can clearly see the states of the dynamical system. That is to say, if the largest Lyapunov exponent is less than 1 corresponding to the black area in Figure 3(d), the dynamical system corresponds to the periodic state or stable state. When the largest Lyapunov exponent is equal to zero corresponding to the yellow area in Figure 3(d), then the dynamical system is at the quasi-period state. If the largest Lyapunov exponent is greater than zero and less than 2 corresponding to the red area, then the dynamical system corresponds to the chaotic state. Finally, when the largest Lyapunov exponent is equal to 2 corresponding to the white area, the dynamical system is at the divergent state. Therefore, by exploiting the 2D Lyapunov exponent diagram, we can easily distinguish the states of system (5) under this group of parameters and verify that ISSs’ structures are indeed Lyapunov stable structures.Figure 3
Fixed the parameters asa=11.6054, c1=2.9154, and c2=2.9639. (a) The double-parameter bifurcation diagram with two bifurcating parameters. (b) The enlargement corresponds to the red box in Figure 3(a), while we increase the number of periods to 200 and abundant dynamic phenomena appear. (c) The enlargement corresponding to the white box in Figure 3(b), where we can clearly see that ISSs (shrimp-like) structures exist, which are labeled by a, b, and c (d). The 2D Lyapunov exponent diagram corresponding to Figure 3(b).
(a)(b)(c)(d)The single-parameter diagrams of the eye-shaped structure with respect to adjustment speeds are displayed in Figure4, where the red curve represents q1 of firm 1 and the blue curve represents q2 of firm 2. As shown in Figure 4(a), this bifurcation diagram corresponds to the red line of the yellow box in Figure 3(b), when we fix the adjustment speed of firm 2 as 0.3193. The dynamic behaviors of this eye-shaped structure are quite complex. We can clearly observe that there is a period-doubling bifurcation at first, and then, system (5) will transiently enter into chaos; finally, a period-halving bifurcation occurs, when v1∈0.345,0.364. This way of bifurcation is also called “hub of periodicity.” Similarly, we can also see this kind of bifurcation in Figure 4(b) when we fix the adjustment speed of firm 2 as 0.3193, which corresponds to the black line of the yellow box in Figure 3(b). When v2∈0.31,0.322, system (5) has a same bifurcation way as that in Figure 4(a). That is to say, the system will experience a period-doubling bifurcation at first and then enter into chaos, and another period-halving bifurcation arises.Figure 4
Fixed the parameters asa=11.6054, c1=2.9154, and c2=2.9639. The red curve represents q1 of firm 1, and the blue curve represents q2 of firm 2. (a) The single-parameter bifurcation diagram corresponds to the red line of the yellow box in Figure 3(b), when v2=0.3193. (b) The single-parameter bifurcation diagram corresponds to the black line of the yellow box in Figure 3(b), when v1=0.3539.
(a)(b)In a nutshell, we can also know that both firms, considered in this paper, enter into chaos through a flip bifurcation, and the flip bifurcation is a typical route for a firm to enter into chaos from a steady state after the institutional or technological reform under an ever-changing era. Furthermore, both the two firms are in the Nash equilibrium state when the adjustment speeds are relatively small. At this time, we can adjust strategy within the stable threshold to maximize their respective interests. However, once the adjustment speed exceeds the stable threshold, the stability of the Nash equilibrium point will be destroyed, and the economic market will enter into chaos after the firms ceaselessly do choice games.


### 4.2. Multistability

The basin of attraction is used to further investigate complex dynamic behaviors of system (5) in long run. It can be seen from the area near v1≈0.45 and v2≈0.35 in Figure 1(a) that the stable period-4 grass green area is doped with a small part of period-12 lake blue area and period-24 dark blue area, which means that there may be the phenomenon of coexisting attractors.We keep fixing the parameters asa=11.1542, c1=1.8044, and c2=1.3949 and choose the adjustment speeds v1 and v2 as bifurcation parameters; then, we can get a group of basins of attraction, as shown in Figure 5. The blue points represent attractors, and their stable region is the yellow area. The invariant cycles are denoted as red, and their stable region is marked as the dark sky blue area. The divergent area and chaotic area are denoted by dark blue. In Figure 5(a), the adjustment speeds of two firms are chosen as v1=0.3605 and v2=0.4517, respectively. There is a coexistence of period-12 attractors and period-4 invariant cycles, and there are many “holes” on boundaries of the basin of attraction. Furthermore, both the values of v1 and v2 increase in Figure 5(b), where v1=0.3625 and v2=0.4527. Apparently, the “holes” at the boundaries of the basin of attraction have not changed significantly. At this time, the complex dynamic behavior of system (5) is depicted as the coexistence of period-24 attractors and period-4 invariant cycles. However, we can clearly see that the dark sky blue stable area is constantly occupied by the yellow attracting area. We can find that the period-4 invariant cycles will keep on becoming bigger until their borders contact the boundary of the yellow-attracting area of periodic attractors; then, a “contact bifurcation” [15] will occur, and the invariant cycles will break and disappear.Figure 5
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. A group of basins of attraction corresponding to system (5) when varying the adjustment speeds v1 and v2. The yellow area is the stable region of blue attractors, and the dark sky blue area is the stable region of red invariant cycles. In addition, the divergent area and chaotic area are denoted by dark blue. (a) v1=0.3605 and v2=0.4517. (b) v1=0.3625 and v2=0.4527.
(a)(b)Furthermore, we assume from the 2D bifurcation diagram in Figure1(a) that there is a phenomenon of coexistence of attractors after choosing appropriate parameters in the parameter space. With the adjustment speeds of these two firms gradually increasing, the number of coexisting attractors increases from 12 to 24, and it means that system (5) may get unstable as v1 and v2 increase. Again, it shows that, in the economic market, system (5) will remain relatively stable when both firms have small adjustment speeds.As shown in Figure6, the values of parameters are chosen as a=10.9788, c1=0.0806, and c2=0.2799, and we similarly choose v1 and v2 as bifurcation parameters. The shapes of these basins of attraction are very peculiar. The green area represents the divergent region, and the points in this region will tend to infinity and finally trap into the chaotic state. The yellow area is the flexible area corresponding to black periodic points or black invariant cycles. In other words, the stable region of these black attractors is the yellow region. The flexible area of other pink attractors is the gray region. In addition, the phenomenon of multistability occurs when v1=1.1940 and v2=0.6810, and there are period-4 pink attractors coexisting with period-12 black attractors, which can be clearly seen in Figure 6(a). The boundary shape of basin of attraction is irregular, and there are many large “holes” in the upper and lower boundaries of basin of attraction. Then, we keep parameter v1 unchanged and increase the value of v2 to 0.6840, and there still is a phenomenon of coexistence between distinct numbers of periodic attractors in Figure 6(b). That is, the period-12 attractors bifurcate to the period-24 attractors via a period-doubling bifurcation. Moreover, viewing this situation from a whole, yellow region occupies a large number of spaces of the gray region, and the whole shape of basin does not significantly change. There is a pretty fractal structure on the head of basin of attraction in Figure 6(b), which is enlarged in Figure 6(c).Figure 6
Fixed the parameters asa=10.9788, c1=0.0806, and c2=0.2799. A group of basins of attraction corresponding to system (5) when varying the adjustment speeds v1 and v2. (a) v1=1.1940 and v2=0.6810. (b)v1=1.1940 and v2=0.6840. (c) The local enlargement corresponding to Figure 6(b) shows a manifest fractal structure. (d) v1=1.1671 and v2=0.7110. (e) v1=1.1671 and v2=0.7264. (f) v1=1.1671 and v2=0.7307.
(a)(b)(c)(d)(e)(f)However, when the bifurcation parameters are varied asv1=1.1671 and v2=0.7110, another different coexisting situation will arise. Four-piece black invariant cycles coexist with period-16 pink points, and there is the phenomenon of “holes” appearing at all boundaries of basin of attraction, as shown in Figure 6(d). The boundaries of invariant cycles almost contact with the borders of their flexible area, and a “contact bifurcation” may happen. In the following steps, we just change the value of adjustment speed v2. It can be clearly seen in Figure 6(e) that the black invariant cycles have completely disappeared due to a contact bifurcation. Instead, only two pieces of chaotic attractors exist inside the basin of attraction, and the “holes” on the boundaries of basin of attraction slightly become much bigger and bigger. Interestingly, these two pieces chaotic attractors gradually become larger and finally connect each other to form one piece chaotic attractor, where there may be a homoclinic bifurcation that happens in Figure 6(f). Consequently, the firms should not substantially increase their adjustment speeds in order to improve their own interests because this will not only not generate considerable profits but also will lead to a disordered competition among firms, then enter into a chaotic state, and never able to return to a stable state. Even if we make small adjustments to the adjustment speed, it will have an irreversible effect on the whole economic system in the market. Therefore, the system will be more stable if the firm maintains a relatively small adjustment speed.


### 4.3. Global Dynamics and Critical Curves

Through the above simple investigations of basin of attraction, we will exploit the method of critical curves to further study some global properties of system (5) and the reason of “holes” formation and its evolution process. It is well known that the term LC (which is a brief expression of the term critical curve) originates from notion of critical points in one-dimensional map. The global bifurcation is the main reason, which causes qualitative change of the basin of attraction. Firstly, the rank-1 preimage LC−1 of critical curves LC is formed by a trajectory of points where the determinant of the Jacobian matrix disappears. In other words, LC and LC−1 satisfy a relation that LC=TLC−1, that is, the points on LC−1 will be mapped to LC through mapping T. In this situation, the points on the LC−1 need to satisfy a condition, that is, the determinant of Jacobian matrix (6) is equal to zero. Those points satisfy the equation as follows:(16)−4α2q2−14q12−3q1+1β1−4α1q1−14q22−3q2+1β2e−q1−q2−2β1β2q2−12q22+q2−4q2+32q1−q222+3q22−12e−2q1−2q2+16α1q1−14α2q2−14=0.Obviously, the expression of equationDetJ=0 is very complex so that we cannot directly obtain analytic solutions with respect to q1 and q2. However, we can graphically display the shape of LC−1 by numerical simulation. Two magenta branches of a hyperbola consist the LC−1, which is denoted as LC−1a and LC−1b, respectively. The critical curve LC is formed by two curves, marked by LCa and LCb, where LCa=TLC−1a and LCb=TLC−1b. It is apparent that these two critical curves separate the plane q1,q2 into three zones Z0, Z2, and Z4. That is to say, LCb separates zone Z0 from zone Z2 and LCa separates zone Z2 from Z4. The points in zone Z0 have no preimages, while the points in zone Z2 have two different rank-1 preimages, and the points in zone Z4 have four distinct rank-1 preimages. It is usual that the origin point O=0,0 always locates in the zone Z4, which means that the origin point O=0,0 has four different rank-1 preimages. One of them is itself, O=0,0. Two of them are located on coordinate axes, O−11=q1∗,0 and O−12=0,q2∗, where q1∗ satisfies the equation β11−q1e−q1=2α1q1 and q2∗ satisfies the equation β21−q2e−q2=2α2q2. The last one is O−13=q1∗∗,q2∗∗, which is the intersection of preimages of two lines OO−11 and OO−12. That is to say, line O−12O−13 is the preimage of line OO−11, and line O−11O−13 is the preimage of line OO−12. q1∗∗ and q2∗∗ simultaneously satisfy the formula as follows:(17)0=β21−q2e−q1+q2−2α2q2,0=β11−q1e−q1+q2−2α1q1.These four preimages of origin pointO=0,0 form the size and vertexes of the basin of attraction.As shown in Figure7, the parameters are fixed as a=11.1542, c1=1.8044, and c2=1.3949. According to above-proposed Propositions 2–4, we can easily judge out that origin point O=0,0 is always an unstable node and both O−11 and O−12 are saddle points. Moreover, O−11 has a stable manifold along line OO−11 on q1-axis, and O−12 also has a stable manifold along line OO−12 on q2-axis. In Figure 7(a), when bifurcation parameters are chosen as v1=0.3625 and v2=0.4517, we can obviously see that only period-2 blue points exist inside the basin of attraction and its stable area is marked as yellow. There are no “holes” on the boundary of basin of attraction because the critical curves, LCa and LCb, are totally inside the basin of attraction. The stable manifolds on the boundaries of basin, OO−11 and OO−12, have not been destroyed. However, one of critical curves, LCb, almost contacts with the boundary line O−12O−13 and already contacts with one of two period-2 points. With the adjustment speed v2 increasing to 0.4257, the period-2 points bifurcate to the period-4 points via a flip bifurcation. Since the critical curve LCb once has crossed rank-1 preimage point O−12, then the stability of stable manifold of O−12 will not exist anymore. It can be clearly seen from Figure 7(b) that the phenomenon “holes” occurs on the line OO−12 and its preimage O−11O−13. This kind of bifurcation which has been mentioned in the previous section can be called a “contact bifurcation.” As v2 constantly increases, a new phenomenon of coexisting attractors arises, that is, period-12 points coexist with period-4 invariant cycles, as shown in Figure 7(c). The size of “holes” apparently becomes bigger than before. Moreover, it is worth noticing that period-4 red invariant cycles gradually get bigger till they contact with a branch of critical curves, LCb. Then, there will be another “contact bifurcation” occurrence. As shown in Figure 7(d), four invariant cycles simultaneously break, and period-12 newly forms into four pieces of chaotic attractors. There are some phenomena of “lakes” that appear inside the basin of attraction. These four pieces of chaotic attractors get bigger as the parameters v1 and v2 increase till they contact with critical curves LCa and LCb. Then, there will be a “final bifurcation” [15], which causes the destruction of four pieces of chaotic attractors and leads the adjustment mechanism to lose its predictability. Although these four pieces of chaotic attractors disappear after a “final bifurcation,” its skeleton still exists inside the original basin of attraction, which is formed by uncountable dense points. Finally, the structure of basin of attraction has been completely destroyed with dense repelling points nested inside. As mentioned above, this phenomenon can be called as “ghost.” Hence, a contact between the boundary of basin of attraction and critical curves or a contact between the boundary of chaotic attractors and critical curves both do have the negative influence on the stability of the structure of the basin of attraction. They may destroy the stability property of system (5) and may cause nonpredictability of tendency of game process.Figure 7
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. The boundaries of basin of attraction are denoted as the black lines. The magenta curves represent preimages of critical curves. The critical curves of system (5) are denoted as green curves, and these two critical curves separate the plane into three zones, labeled by Z0, Z2, and Z4. (a) v1=0.3625 and v2=0.4057. (b) v1=0.3625 and v2=0.4257. (c) v1=0.3625 and v2=0.4517. (d) v1=0.3875 and v2=0.4517.
(a)(b)(c)(d)A group of basins of attraction with interesting shapes of “holes” located on the boundary can be observed, when we fix the parameters asa=11.1542, v1=0.3625, and v2=0.4517 and choose parameters c1 and c2 as bifurcation parameters. The meanings of magenta curves, green curves, and two black lines are same as in Figure 7. As shown in Figure 8, we change the values of production cost of two firms to investigate their influence on the shape and stability property of the basins of attraction. The origin point O=0,0 is still an unstable node, and its manifolds along both q1-axis and q2-axis which are unstable, while both O−11 and O−12 are saddle points so that O−11 has a stable manifolds along line OO−11 on q1-axis and O−12 also has a stable manifold along line OO−12 on q2-axis. As shown in Figure 8(a), when we fix the bifurcation parameters as c1=2.8860 and c2=0.9090, there are period-12 blue periodic points coexisting with period-4 black periodic points, where the stable areas corresponding to period-4 blue points and period-4 black points are marked as pink and orange, respectively. The blue area represents the divergent area. Another significant phenomenon is that a branch of critical curves, LCb, have already crossed both O−11 and O−12, which means that two stable manifolds along line OO−11 on q1-axis and line OO−12 on q2-axis both have been completely destroyed. As a consequence, the phenomena of “holes” in the shape of aliens, but of different sizes, appear on all boundaries of the basin of attraction.Because the critical curve LCb has contacted with period-12 blue points, and then, it is highly possible that a “contact bifurcation” may happen. Then, we fix the value of c1 unchanged and increase the value of c2 little by little to 0.9590, so we can clearly see in Figure 8(b) that a “contact bifurcation” do occurs, and there are phenomena of “lakes” that also appear besides the phenomena of “holes” on boundaries of basin. In addition, period-4 black points coexist with period-12 chaotic attractors, where the stable attracting area of period-4 black points is still the orange area and the stable attracting area of chaotic attractors is the pink area. In this situation, the critical curve LCb contacts with chaotic attractors. Moreover, after this kind of contact, the chaotic attractors will first become bigger for a very short time and then break but still exist in original basin of attraction, as shown in Figure 8(c). However, this phenomenon is not so called “ghost,” because the period-4 black points still exist inside the basin of attraction and have their own domain of attraction in plane q1,q2. We also can clearly see that the “lakes” notably grow bigger.Figure 8
Fixed the parameters asa=11.1542, v1=0.3625, and v2=0.4517. The boundaries of basin of attraction are denoted as the black lines. The magenta curves represent preimages of critical curves. The critical curves of system (5) are denoted as green curves, and these two critical curves separate the plane into three zones, labeled by Z0, Z2, and Z4. (a) c1=2.8860 and c2=0.9090. (b) c1=2.8860 and c2=0.9590. (c) c1=2.8860 and c2=1.0000. (d) c1=2.8860 and c2=1.0750.
(a)(b)(c)(d)As the value ofc2 constantly increases to 1.0750, period-4 points bifurcate to two pieces of chaotic attractors, and its stable area is still the orange area, while the pink area disappears, which can be clearly displayed in Figure 8(d). The “lakes” inside the basin of attraction have been connected with the “holes” on the boundary of basin of attraction and newly form interesting shapes of “holes” located on all sides of basin of attraction. It is apparent that the boundary of two pieces of black chaotic attractors almost contacts again with both two branches of critical curves. Finally, there will be a “final bifurcation” occurrence. Consequently, these two pieces of chaotic attractors simultaneously break out but their skeletons of the basin of attraction still exist, which are formed by infinitely many unstable points. Hence, increasing the value of a pair of parameters c1 and c2 similarly has significant influence on the stable property and structure of the basin of system (5). The firms in the market better consider to take some actions to decrease its production cost, as one of them has high cost that may irreversibly lead the game to a chaotic state so that the future destiny of firms cannot be forecasted.


## 4.1. Complex Dynamic Behaviors

First, a set of parameters is fixed asa=11.1542, c1=1.8044, and c2=1.3949, and we choose adjustment speeds v1 and v2 as the bifurcation parameters. As shown in Figure 1, the color bar is on the right side of figures, in which different colors denote different numbers of the period. The steel blue represents the stable area of the Nash equilibrium point, and the orange represents that the points in this area are at a period-2 state, and analogously, the grass green represents the period-4. The divergent area is denoted as the color white, which is also marked as “ET” in the color bar. However, the meaning of the dark black area is more complex, where the dynamic states included in this area, such as the quasi-period state and chaotic state. From Figure 1(a), abundant dynamic phenomena can be observed. As the values of v1 and v2 gradually increase, the bifurcation route of system (5) will start from the period-1 steel blue stable area, then pass through the period-2 orange area and period-4 grass green region, and finally enter into the dark black region, which means that system (5) will change from a stable state to an unstable state. Corresponding to the real economic market, this phenomenon means that if two firms blindly improve their adjustment speeds without strategy, the market will become unstable and eventually may enter into chaos. Furthermore, Figure 1(b) is a local enlargement of Figure 1(a). The structures of shrimp-like are called Isoperiodic Stable Structures (the abbreviation is ISSs), seen in Ref. [27], which can be clearly shown in Figure 1(b). The location where the shrimp-like structures are at is denoted by a white straight line, and the expression of this line is v2=−5.6v1+1.904.Figure 1
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. (a) The double-parameters’ bifurcation diagram of system (5). (b) The local enlargement of Figure 1(a).
(a)(b)In Figure2, the red curve denotes the outputs q1 of firm 1, and the black curve denotes the outputs q2 of firm 2. First, choose the adjustment speed v1 as the bifurcating parameter, and the other parameters are fixed as a=11.1542, c1=1.8044, c2=1.3949, and v2=0.1. Under this group of parameters, system (5) will go through a flip bifurcation at first and finally drop into a chaotic state, as shown in Figure 2(a). The bifurcation diagram of q2 with respect to v1, marked by a black box, is displayed in the enlargement. From Figure 2(a), when the adjustment speed of firm 1 satisfies v1<0.4197, system (5) is at a period-1 stable state. In addition, v1=0.4197 is a bifurcation point, which means that the stability of system (5) at this point will completely change, and then, a period-doubling bifurcation subsequently occurs. Once the value of v1 is greater than 0.5371, system (5) will uncontrollably reach a disordered chaotic state.Figure 2
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. The red curve denotes q1 of firm 1, and the black curve denotes q2 of firm 2. (a) The 1D bifurcation diagram when v2=0.1. (b) Corresponding to Figure 1(b), the 1D bifurcation diagram on the white straight line is displayed with its largest Lyapunov exponent diagram. Fixed the parameters as a=11.1542, v1=0.3625, and v2=0.4517. (c) The 1D bifurcation diagram when c2=1.8. (d) The 1D bifurcation diagram when c1=1.8.
(a)(b)(c)(d)Under same parameter conditions, it can be seen from Figure2(b) that the dynamic behaviors of system (5) are very complicated when shrimp-like structures appear. At this time, the blue curve represents the largest Lyapunov exponent curve corresponding to the red bifurcation curve. It is well known that we can easily determine the state of a dynamical system through calculating its Lyapunov exponents, and we usually need to calculate the largest one. Therefore, the dynamical system corresponds to a stable state, if the largest Lyapunov exponent is less than zero. Instead, if the largest Lyapunov exponent is greater than zero, it means that the dynamical system is at a chaotic state. In addition, it is worth noting that the largest Lyapunov exponents fluctuate up and down around zero, and then, the dynamical system is at a quasi-period state. There are many phenomena of “periodic window” that exist, and system (5) constantly transforms between the chaotic state and the periodic stable state. With the parameter v1 further increasing, system (5) will enter a stable state from a chaotic state for a short time and finally return to the chaotic state.When the parameters are changed toa=11.1542, v1=0.3625, and v2=0.4517, then we can get two bifurcation diagrams with respect to costs of two firms, respectively. As shown in Figure 2(c), where c2=1.8, when c1<0.7755, system (5) is at the period-1 stable state. At c1=0.7755, a typical flip bifurcation (also called period-doubling bifurcation) occurs. In addition, it is worth to notice that an obvious phenomenon called “jump” appears near c1=1.8723, which means that there may be a phenomenon of coexistence. Similarly, we choose c1=1.8 in Figure 2(d),then we can know that c2=0.6330 is a bifurcation point, and once system (5) crosses this point, then a flip bifurcation will occur. The dynamic phenomena between 1.3500<c2<1.4180 are very complicated. There is also a “jump” phenomenon that exists near c2=1.6083. With constant increase of costs of two firms, dynamical system (5) will eventually enter into chaos. However, through comparing Figures 2(c) and 2(d), we can draw a result that, under the premise of same costs, the speed of the system entering into chaos with c2 is higher than that with c1.Another group of parameters are fixed asa=11.6054, c1=2.9154, and c2=2.9639, and we choose adjustment speeds v1 and v2 as the bifurcating parameters. A totally different 2D bifurcation diagram is shown in Figure 3(a), in which many colorful and symmetrical structures exist. Similarly, the dark black area includes complex phenomena such as chaos and quasiperiodic states, and the white area denotes the divergent area. It is worth noticing that the dynamic behaviors, marked in a red box, have significances to investigate. Hence, we enlarge it in Figure 3(b) and increase the number of periods to 200. As shown in Figure 3(b), abundant and complex dynamic phenomena are displayed. There are many fractal structures that can be seen, and the eye-shaped structure is symmetrically arranged on both sides of the ribbon. The dynamical behaviors, located in a white box of Figure 3(b), are enlarged in Figure 3(c). Obviously, there are apparent V-shaped islands called ISSs structures, which are marked as a, b, and c in the yellow box. It is known that this kind of structures are Lyapunov stable structures or Lyapunov stable islands, that is, if parameters are chosen in these structures, the state of the dynamical system is regular and stable, seen in Ref. [27]. Since the 2D Lyapunov exponent diagram is a useful tool to clearly distinguish stable state, quasi-periodic state, chaotic state, and divergent state, we can clearly see the states of the dynamical system. That is to say, if the largest Lyapunov exponent is less than 1 corresponding to the black area in Figure 3(d), the dynamical system corresponds to the periodic state or stable state. When the largest Lyapunov exponent is equal to zero corresponding to the yellow area in Figure 3(d), then the dynamical system is at the quasi-period state. If the largest Lyapunov exponent is greater than zero and less than 2 corresponding to the red area, then the dynamical system corresponds to the chaotic state. Finally, when the largest Lyapunov exponent is equal to 2 corresponding to the white area, the dynamical system is at the divergent state. Therefore, by exploiting the 2D Lyapunov exponent diagram, we can easily distinguish the states of system (5) under this group of parameters and verify that ISSs’ structures are indeed Lyapunov stable structures.Figure 3
Fixed the parameters asa=11.6054, c1=2.9154, and c2=2.9639. (a) The double-parameter bifurcation diagram with two bifurcating parameters. (b) The enlargement corresponds to the red box in Figure 3(a), while we increase the number of periods to 200 and abundant dynamic phenomena appear. (c) The enlargement corresponding to the white box in Figure 3(b), where we can clearly see that ISSs (shrimp-like) structures exist, which are labeled by a, b, and c (d). The 2D Lyapunov exponent diagram corresponding to Figure 3(b).
(a)(b)(c)(d)The single-parameter diagrams of the eye-shaped structure with respect to adjustment speeds are displayed in Figure4, where the red curve represents q1 of firm 1 and the blue curve represents q2 of firm 2. As shown in Figure 4(a), this bifurcation diagram corresponds to the red line of the yellow box in Figure 3(b), when we fix the adjustment speed of firm 2 as 0.3193. The dynamic behaviors of this eye-shaped structure are quite complex. We can clearly observe that there is a period-doubling bifurcation at first, and then, system (5) will transiently enter into chaos; finally, a period-halving bifurcation occurs, when v1∈0.345,0.364. This way of bifurcation is also called “hub of periodicity.” Similarly, we can also see this kind of bifurcation in Figure 4(b) when we fix the adjustment speed of firm 2 as 0.3193, which corresponds to the black line of the yellow box in Figure 3(b). When v2∈0.31,0.322, system (5) has a same bifurcation way as that in Figure 4(a). That is to say, the system will experience a period-doubling bifurcation at first and then enter into chaos, and another period-halving bifurcation arises.Figure 4
Fixed the parameters asa=11.6054, c1=2.9154, and c2=2.9639. The red curve represents q1 of firm 1, and the blue curve represents q2 of firm 2. (a) The single-parameter bifurcation diagram corresponds to the red line of the yellow box in Figure 3(b), when v2=0.3193. (b) The single-parameter bifurcation diagram corresponds to the black line of the yellow box in Figure 3(b), when v1=0.3539.
(a)(b)In a nutshell, we can also know that both firms, considered in this paper, enter into chaos through a flip bifurcation, and the flip bifurcation is a typical route for a firm to enter into chaos from a steady state after the institutional or technological reform under an ever-changing era. Furthermore, both the two firms are in the Nash equilibrium state when the adjustment speeds are relatively small. At this time, we can adjust strategy within the stable threshold to maximize their respective interests. However, once the adjustment speed exceeds the stable threshold, the stability of the Nash equilibrium point will be destroyed, and the economic market will enter into chaos after the firms ceaselessly do choice games.


## 4.2. Multistability

The basin of attraction is used to further investigate complex dynamic behaviors of system (5) in long run. It can be seen from the area near v1≈0.45 and v2≈0.35 in Figure 1(a) that the stable period-4 grass green area is doped with a small part of period-12 lake blue area and period-24 dark blue area, which means that there may be the phenomenon of coexisting attractors.We keep fixing the parameters asa=11.1542, c1=1.8044, and c2=1.3949 and choose the adjustment speeds v1 and v2 as bifurcation parameters; then, we can get a group of basins of attraction, as shown in Figure 5. The blue points represent attractors, and their stable region is the yellow area. The invariant cycles are denoted as red, and their stable region is marked as the dark sky blue area. The divergent area and chaotic area are denoted by dark blue. In Figure 5(a), the adjustment speeds of two firms are chosen as v1=0.3605 and v2=0.4517, respectively. There is a coexistence of period-12 attractors and period-4 invariant cycles, and there are many “holes” on boundaries of the basin of attraction. Furthermore, both the values of v1 and v2 increase in Figure 5(b), where v1=0.3625 and v2=0.4527. Apparently, the “holes” at the boundaries of the basin of attraction have not changed significantly. At this time, the complex dynamic behavior of system (5) is depicted as the coexistence of period-24 attractors and period-4 invariant cycles. However, we can clearly see that the dark sky blue stable area is constantly occupied by the yellow attracting area. We can find that the period-4 invariant cycles will keep on becoming bigger until their borders contact the boundary of the yellow-attracting area of periodic attractors; then, a “contact bifurcation” [15] will occur, and the invariant cycles will break and disappear.Figure 5
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. A group of basins of attraction corresponding to system (5) when varying the adjustment speeds v1 and v2. The yellow area is the stable region of blue attractors, and the dark sky blue area is the stable region of red invariant cycles. In addition, the divergent area and chaotic area are denoted by dark blue. (a) v1=0.3605 and v2=0.4517. (b) v1=0.3625 and v2=0.4527.
(a)(b)Furthermore, we assume from the 2D bifurcation diagram in Figure1(a) that there is a phenomenon of coexistence of attractors after choosing appropriate parameters in the parameter space. With the adjustment speeds of these two firms gradually increasing, the number of coexisting attractors increases from 12 to 24, and it means that system (5) may get unstable as v1 and v2 increase. Again, it shows that, in the economic market, system (5) will remain relatively stable when both firms have small adjustment speeds.As shown in Figure6, the values of parameters are chosen as a=10.9788, c1=0.0806, and c2=0.2799, and we similarly choose v1 and v2 as bifurcation parameters. The shapes of these basins of attraction are very peculiar. The green area represents the divergent region, and the points in this region will tend to infinity and finally trap into the chaotic state. The yellow area is the flexible area corresponding to black periodic points or black invariant cycles. In other words, the stable region of these black attractors is the yellow region. The flexible area of other pink attractors is the gray region. In addition, the phenomenon of multistability occurs when v1=1.1940 and v2=0.6810, and there are period-4 pink attractors coexisting with period-12 black attractors, which can be clearly seen in Figure 6(a). The boundary shape of basin of attraction is irregular, and there are many large “holes” in the upper and lower boundaries of basin of attraction. Then, we keep parameter v1 unchanged and increase the value of v2 to 0.6840, and there still is a phenomenon of coexistence between distinct numbers of periodic attractors in Figure 6(b). That is, the period-12 attractors bifurcate to the period-24 attractors via a period-doubling bifurcation. Moreover, viewing this situation from a whole, yellow region occupies a large number of spaces of the gray region, and the whole shape of basin does not significantly change. There is a pretty fractal structure on the head of basin of attraction in Figure 6(b), which is enlarged in Figure 6(c).Figure 6
Fixed the parameters asa=10.9788, c1=0.0806, and c2=0.2799. A group of basins of attraction corresponding to system (5) when varying the adjustment speeds v1 and v2. (a) v1=1.1940 and v2=0.6810. (b)v1=1.1940 and v2=0.6840. (c) The local enlargement corresponding to Figure 6(b) shows a manifest fractal structure. (d) v1=1.1671 and v2=0.7110. (e) v1=1.1671 and v2=0.7264. (f) v1=1.1671 and v2=0.7307.
(a)(b)(c)(d)(e)(f)However, when the bifurcation parameters are varied asv1=1.1671 and v2=0.7110, another different coexisting situation will arise. Four-piece black invariant cycles coexist with period-16 pink points, and there is the phenomenon of “holes” appearing at all boundaries of basin of attraction, as shown in Figure 6(d). The boundaries of invariant cycles almost contact with the borders of their flexible area, and a “contact bifurcation” may happen. In the following steps, we just change the value of adjustment speed v2. It can be clearly seen in Figure 6(e) that the black invariant cycles have completely disappeared due to a contact bifurcation. Instead, only two pieces of chaotic attractors exist inside the basin of attraction, and the “holes” on the boundaries of basin of attraction slightly become much bigger and bigger. Interestingly, these two pieces chaotic attractors gradually become larger and finally connect each other to form one piece chaotic attractor, where there may be a homoclinic bifurcation that happens in Figure 6(f). Consequently, the firms should not substantially increase their adjustment speeds in order to improve their own interests because this will not only not generate considerable profits but also will lead to a disordered competition among firms, then enter into a chaotic state, and never able to return to a stable state. Even if we make small adjustments to the adjustment speed, it will have an irreversible effect on the whole economic system in the market. Therefore, the system will be more stable if the firm maintains a relatively small adjustment speed.


## 4.3. Global Dynamics and Critical Curves

Through the above simple investigations of basin of attraction, we will exploit the method of critical curves to further study some global properties of system (5) and the reason of “holes” formation and its evolution process. It is well known that the term LC (which is a brief expression of the term critical curve) originates from notion of critical points in one-dimensional map. The global bifurcation is the main reason, which causes qualitative change of the basin of attraction. Firstly, the rank-1 preimage LC−1 of critical curves LC is formed by a trajectory of points where the determinant of the Jacobian matrix disappears. In other words, LC and LC−1 satisfy a relation that LC=TLC−1, that is, the points on LC−1 will be mapped to LC through mapping T. In this situation, the points on the LC−1 need to satisfy a condition, that is, the determinant of Jacobian matrix (6) is equal to zero. Those points satisfy the equation as follows:(16)−4α2q2−14q12−3q1+1β1−4α1q1−14q22−3q2+1β2e−q1−q2−2β1β2q2−12q22+q2−4q2+32q1−q222+3q22−12e−2q1−2q2+16α1q1−14α2q2−14=0.Obviously, the expression of equationDetJ=0 is very complex so that we cannot directly obtain analytic solutions with respect to q1 and q2. However, we can graphically display the shape of LC−1 by numerical simulation. Two magenta branches of a hyperbola consist the LC−1, which is denoted as LC−1a and LC−1b, respectively. The critical curve LC is formed by two curves, marked by LCa and LCb, where LCa=TLC−1a and LCb=TLC−1b. It is apparent that these two critical curves separate the plane q1,q2 into three zones Z0, Z2, and Z4. That is to say, LCb separates zone Z0 from zone Z2 and LCa separates zone Z2 from Z4. The points in zone Z0 have no preimages, while the points in zone Z2 have two different rank-1 preimages, and the points in zone Z4 have four distinct rank-1 preimages. It is usual that the origin point O=0,0 always locates in the zone Z4, which means that the origin point O=0,0 has four different rank-1 preimages. One of them is itself, O=0,0. Two of them are located on coordinate axes, O−11=q1∗,0 and O−12=0,q2∗, where q1∗ satisfies the equation β11−q1e−q1=2α1q1 and q2∗ satisfies the equation β21−q2e−q2=2α2q2. The last one is O−13=q1∗∗,q2∗∗, which is the intersection of preimages of two lines OO−11 and OO−12. That is to say, line O−12O−13 is the preimage of line OO−11, and line O−11O−13 is the preimage of line OO−12. q1∗∗ and q2∗∗ simultaneously satisfy the formula as follows:(17)0=β21−q2e−q1+q2−2α2q2,0=β11−q1e−q1+q2−2α1q1.These four preimages of origin pointO=0,0 form the size and vertexes of the basin of attraction.As shown in Figure7, the parameters are fixed as a=11.1542, c1=1.8044, and c2=1.3949. According to above-proposed Propositions 2–4, we can easily judge out that origin point O=0,0 is always an unstable node and both O−11 and O−12 are saddle points. Moreover, O−11 has a stable manifold along line OO−11 on q1-axis, and O−12 also has a stable manifold along line OO−12 on q2-axis. In Figure 7(a), when bifurcation parameters are chosen as v1=0.3625 and v2=0.4517, we can obviously see that only period-2 blue points exist inside the basin of attraction and its stable area is marked as yellow. There are no “holes” on the boundary of basin of attraction because the critical curves, LCa and LCb, are totally inside the basin of attraction. The stable manifolds on the boundaries of basin, OO−11 and OO−12, have not been destroyed. However, one of critical curves, LCb, almost contacts with the boundary line O−12O−13 and already contacts with one of two period-2 points. With the adjustment speed v2 increasing to 0.4257, the period-2 points bifurcate to the period-4 points via a flip bifurcation. Since the critical curve LCb once has crossed rank-1 preimage point O−12, then the stability of stable manifold of O−12 will not exist anymore. It can be clearly seen from Figure 7(b) that the phenomenon “holes” occurs on the line OO−12 and its preimage O−11O−13. This kind of bifurcation which has been mentioned in the previous section can be called a “contact bifurcation.” As v2 constantly increases, a new phenomenon of coexisting attractors arises, that is, period-12 points coexist with period-4 invariant cycles, as shown in Figure 7(c). The size of “holes” apparently becomes bigger than before. Moreover, it is worth noticing that period-4 red invariant cycles gradually get bigger till they contact with a branch of critical curves, LCb. Then, there will be another “contact bifurcation” occurrence. As shown in Figure 7(d), four invariant cycles simultaneously break, and period-12 newly forms into four pieces of chaotic attractors. There are some phenomena of “lakes” that appear inside the basin of attraction. These four pieces of chaotic attractors get bigger as the parameters v1 and v2 increase till they contact with critical curves LCa and LCb. Then, there will be a “final bifurcation” [15], which causes the destruction of four pieces of chaotic attractors and leads the adjustment mechanism to lose its predictability. Although these four pieces of chaotic attractors disappear after a “final bifurcation,” its skeleton still exists inside the original basin of attraction, which is formed by uncountable dense points. Finally, the structure of basin of attraction has been completely destroyed with dense repelling points nested inside. As mentioned above, this phenomenon can be called as “ghost.” Hence, a contact between the boundary of basin of attraction and critical curves or a contact between the boundary of chaotic attractors and critical curves both do have the negative influence on the stability of the structure of the basin of attraction. They may destroy the stability property of system (5) and may cause nonpredictability of tendency of game process.Figure 7
Fixed the parameters asa=11.1542, c1=1.8044, and c2=1.3949. The boundaries of basin of attraction are denoted as the black lines. The magenta curves represent preimages of critical curves. The critical curves of system (5) are denoted as green curves, and these two critical curves separate the plane into three zones, labeled by Z0, Z2, and Z4. (a) v1=0.3625 and v2=0.4057. (b) v1=0.3625 and v2=0.4257. (c) v1=0.3625 and v2=0.4517. (d) v1=0.3875 and v2=0.4517.
(a)(b)(c)(d)A group of basins of attraction with interesting shapes of “holes” located on the boundary can be observed, when we fix the parameters asa=11.1542, v1=0.3625, and v2=0.4517 and choose parameters c1 and c2 as bifurcation parameters. The meanings of magenta curves, green curves, and two black lines are same as in Figure 7. As shown in Figure 8, we change the values of production cost of two firms to investigate their influence on the shape and stability property of the basins of attraction. The origin point O=0,0 is still an unstable node, and its manifolds along both q1-axis and q2-axis which are unstable, while both O−11 and O−12 are saddle points so that O−11 has a stable manifolds along line OO−11 on q1-axis and O−12 also has a stable manifold along line OO−12 on q2-axis. As shown in Figure 8(a), when we fix the bifurcation parameters as c1=2.8860 and c2=0.9090, there are period-12 blue periodic points coexisting with period-4 black periodic points, where the stable areas corresponding to period-4 blue points and period-4 black points are marked as pink and orange, respectively. The blue area represents the divergent area. Another significant phenomenon is that a branch of critical curves, LCb, have already crossed both O−11 and O−12, which means that two stable manifolds along line OO−11 on q1-axis and line OO−12 on q2-axis both have been completely destroyed. As a consequence, the phenomena of “holes” in the shape of aliens, but of different sizes, appear on all boundaries of the basin of attraction.Because the critical curve LCb has contacted with period-12 blue points, and then, it is highly possible that a “contact bifurcation” may happen. Then, we fix the value of c1 unchanged and increase the value of c2 little by little to 0.9590, so we can clearly see in Figure 8(b) that a “contact bifurcation” do occurs, and there are phenomena of “lakes” that also appear besides the phenomena of “holes” on boundaries of basin. In addition, period-4 black points coexist with period-12 chaotic attractors, where the stable attracting area of period-4 black points is still the orange area and the stable attracting area of chaotic attractors is the pink area. In this situation, the critical curve LCb contacts with chaotic attractors. Moreover, after this kind of contact, the chaotic attractors will first become bigger for a very short time and then break but still exist in original basin of attraction, as shown in Figure 8(c). However, this phenomenon is not so called “ghost,” because the period-4 black points still exist inside the basin of attraction and have their own domain of attraction in plane q1,q2. We also can clearly see that the “lakes” notably grow bigger.Figure 8
Fixed the parameters asa=11.1542, v1=0.3625, and v2=0.4517. The boundaries of basin of attraction are denoted as the black lines. The magenta curves represent preimages of critical curves. The critical curves of system (5) are denoted as green curves, and these two critical curves separate the plane into three zones, labeled by Z0, Z2, and Z4. (a) c1=2.8860 and c2=0.9090. (b) c1=2.8860 and c2=0.9590. (c) c1=2.8860 and c2=1.0000. (d) c1=2.8860 and c2=1.0750.
(a)(b)(c)(d)As the value ofc2 constantly increases to 1.0750, period-4 points bifurcate to two pieces of chaotic attractors, and its stable area is still the orange area, while the pink area disappears, which can be clearly displayed in Figure 8(d). The “lakes” inside the basin of attraction have been connected with the “holes” on the boundary of basin of attraction and newly form interesting shapes of “holes” located on all sides of basin of attraction. It is apparent that the boundary of two pieces of black chaotic attractors almost contacts again with both two branches of critical curves. Finally, there will be a “final bifurcation” occurrence. Consequently, these two pieces of chaotic attractors simultaneously break out but their skeletons of the basin of attraction still exist, which are formed by infinitely many unstable points. Hence, increasing the value of a pair of parameters c1 and c2 similarly has significant influence on the stable property and structure of the basin of system (5). The firms in the market better consider to take some actions to decrease its production cost, as one of them has high cost that may irreversibly lead the game to a chaotic state so that the future destiny of firms cannot be forecasted.


## 5. Conclusions

In this paper, on the basis of considering the exponential demand function, the cost function is improved to quadratic cost. The firms considered in the market both have incomplete market information. Therefore, they are boundedly rational. Then, a discrete-time evolution model of the duopoly game with bounded rationality is built. The gradient adjustment mechanism is introduced to adjust the production strategy. Although we have analyzed the local stability property of both boundary equilibrium points and Nash equilibrium point, we cannot obtain a clear analytic solution. Hence, we utilize numerical simulation to further study more properties and the influence of parameters on the stability of the established model. With the parameters continually changing, there are complex dynamic phenomena that arise, such as V-shaped structures (ISSs), eye-shaped structures, and the coexistence between different periods periodic points, or between periodic points and invariant cycles, or between periodic points and chaotic attractors. Through the method of critical curves, we describe the change process of the structure of basin of attraction and then the causes of “holes” formation. As a result, once the parameters exceed the stable threshold, the global bifurcations such as “contact bifurcation” and “final bifurcation” may appear, and then, the market will irreversibly enter into a chaotic state after the firms ceaselessly do choice games.



---

*Source: 9978320-2021-05-30.xml*