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

In this paper, a discrete-time dynamic duopoly model, with nonlinear demand and cost functions, is established.,e properties of existence and local stability of equilibrium points have been verified and analyzed. ,e 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. ,e eye-shaped structures appear where the period-doubling and periodhalving 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.


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 totalcost 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." e 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. erefore, 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. eir 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 twostage duopoly game with differentiated product and R&D spillovers. Cavalli and Naimzada [20] also had described complex dynamic behaviors and the phenomenon of multistability. e 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. e 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 singleparameter (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. e rest parts of this paper are organized as follows. In Section 2, 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. e existence of equilibrium points has been verified. e 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.

Modeling and the Existence of Equilibria
It is supposed that there are two firms i (labeled by i � 1 and 2) in the economic market, producing homogeneous products, and both these two kinds of homogeneous products are supplied to a same market. e 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 i(i � 1 and 2) is an exponential inverse demand function defined by the outputs q i ≥ 0 of two firms, that is, P(Q) � ae − (q 1 +q 2 ) , in which a > 0 is the maximum selling price of the products in the market and Q � q 1 + q 2 is the number of total quantity of these two kinds of homogeneous products. e 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 2 Complexity form, that is, C i (q i ) � c i q 2 i , where c i > 0 is the marginal cost, and the relation between a and c i satisfies a > c i . Hence, the profits of these two firms can be, respectively, written as follows: (1) Let i , respectively, deduce the derivative with regard to the output q i , i � 1 and 2. en, we can get the corresponding marginal profit functions as follows: 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 firm i(i � 1 and 2) is boundedly rational. e 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 q i . If the marginal profit in t period is equal to zero, then the firm i will keep its volume of production unchanged. e "gradient adjustment mechanism" is a useful method to simulate how the firms adjust their decisions so as to make sure to get maximum interests.
us, 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: in which v i > 0 represents the adjustment speed of the firm i. rough substituting equation (2) into equation (3), the discrete-time dynamic adjustment mechanism of the two firms can be specifically expressed as follows: For simplicity, the auxiliary variables have been introduced in system (4), that is, α i � c i v i and β i � av i (i � 1 and 2), and we have that α i and β i ∈ (0, +∞). We introduce an evolutionary operator T and auxiliary variables as α i � c i v i and β i � av i , i � 1 and 2, where α i and β i ∈ (0, +∞). Hence, a discrete-time dynamic duopoly Cournot game can be further simplified as the following: in which T is an evolutionary operator. According to the knowledge of nonlinear dynamics, let q i (t + 1) � q i (t), i � 1 and 2, so we can know that the equilibrium point in system (5) must satisfy the conditions listed as follows: (i) q 1 � 0 or (ii) β 1 (1 − q 1 )e − (q 1 +q 2 ) � 2α 1 q 1 (iii) q 2 � 0 or (iv) β 2 (1 − q 2 )e − (q 1 +q 2 ) � 2α 2 q 2 On 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 E 1 � (0, 0) of system (5) can be obtained. For the rest of the boundary equilibrium points, we can draw the following conclusions: (1) e boundary equilibrium point E 2 � (x * , 0) can be obtained by condition (ii) and condition (iii), where x * satisfies the equation β 1 (1 − x * )e − x * � 2α 1 x * . If we rewrite this relationship as a function with regard Complexity to a variable z, that is, g(z) � 2α 1 z − (β 1 (1 − z)/e z ), then we have g(0) � − β 1 < 0 and g(1) � 2α 1 > 0. So we can know that there is at least one point Finally, the existence of E 2 � (x * , 0) is verified. (2) Analogously, the boundary equilibrium point E 3 � (0, y * ) can be obtained by condition (i) and condition (iv), where y * satisfies the equation β 2 (1 − y * )e − y * � 2α 2 y * . If we rewrite this equation as a function with respect to a variable z, that is, w(z) � 2α 2 z − (β 2 (1 − z)/e z ), then we will get two inequalities that w(0) � − β 2 < 0 and w(1) � 2α 2 > 0. So we can know that there is at least one point y * ∈ (0, 1) that satisfies β 2 (1 − y * )e − y * � 2α 2 y * . Finally, the existence of E 3 � (0, y * ) is also verified.
Hence, we do verify the existence of these three boundary equilibrium points. en, the analysis of the Nash equilibrium point and its existence in system (5) is given by the proposition as the following.
In a real economic system, the equilibrium point E 1 � (0, 0) represents the situation that these two firms both go bankrupt or get out of the market. e boundary equilibrium points E 2 � (x * , 0) and E 3 � (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.

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 |λ i |(i � 1 and 2) of the Jacobian matrix.
at is to say, if |λ i | > 1(i � 1 and 2), then the equilibrium point is an unstable node. If |λ i | < 1(i � 1 and 2), then the equilibrium point is a stable node. While if one of |λ i |(i � 1, 2) is great than 1 and the other is smaller than 1, then the equilibrium point is a saddle point. e Jacobian matrix corresponding to system (5) can be expressed as follows: e local stability of boundary equilibrium points of system (5) has been discussed, and the following propositions can be drawn.
As for the unique Nash equilibrium of system (5), the corresponding Jacobian matrix at E 4 � (m * , n * ) can be expressed as in which m * satisfies the equation 2α 1 m * � β 1 (1 − m * )e − (m * +n * ) and n * satisfies the equation 2α 2 n * � β 2 (1 − n * )e − (m * +n * ) . erefore, the trace and the determinant of the Jacobian matrix J(E 4 ) are denoted as Tr and Det, respectively. e specific expressions corresponding to Tr and Det can be, respectively, written as follows: As we all know, the Nash equilibrium point E 4 is stable when the moduli of its eigenvalues corresponding to the Jacobian matrix J(E 4 ) are all less than 1, while the forms of eigenvalues under the context of the exponential term are quite sophisticated. erefore, we need to study the local stability around the Nash equilibrium E 4 with the help of the Jury criterion.
Proof. According to the Jury criterion, we can verify the local stability of E 4 . In other words, if all of conditions of the Jury criterion hold, then we can know that the Nash equilibrium E 4 is asymptotically stable. e specific expression of the Jury criterion is given as follows: Substituting the specific formulations of Tr and Det into equation (12), then we will get equivalent expressions as follows: Since the conditions m * ∈ (0, 1) and n * ∈ (0, 1), then we can know that A > 0, B < 0, C < 0, and D > 0. en, the three conditions in equation (13) are equivalent to 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. erefore, 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 E 4 . Obviously, the condition (ii) ′ is equivalent to the form of In addition, in order to further analyze the local stability of E 4 and get the value ranges of the parameters β 1 and β 2 , the inequality β 1 (− 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.
e 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.
e situation is discussed as following, when the equivalent condition (iii) ′ in the Jury criterion holds. Analogously, we will get an inequality β 2 e A C + β 1 (e A B + β 2 D) > 0, and then, it will be discussed by divided into two cases, that is, β 2 e A C � 0 and β 1 (e A B + β 2 D) � 0. In other words, e discussion in terms of β 1 has two situations. (iv′) when /(e A B + β 2 D)); we already know that ((β 2 e A (− C))/(e A B + β 2 D)) < 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 ′ ). at is to say, if one of the above six conditions in Proposition 4 holds, then E 4 � (m * , n * ) is locally asymptotical stable.

Numerical Simulation
In the nonlinear dynamical system, the phase diagram or bifurcation diagram may include a lot of information in 6 Complexity 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. erefore, 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. e 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.

Complex Dynamic Behaviors.
First, a set of parameters is fixed as a � 11.1542, c 1 � 1.8044, and c 2 � 1.3949, and we choose adjustment speeds v 1 and v 2 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. e 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. e 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 v 1 and v 2 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). e 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). e location where the shrimp-like structures are at is denoted by a white straight line, and the expression of this line is v 2 � − 5.6v 1 + 1.904.
In Figure 2, the red curve denotes the outputs q 1 of firm 1, and the black curve denotes the outputs q 2 of firm 2. First, choose the adjustment speed v 1 as the bifurcating parameter, and the other parameters are fixed as a � 11.1542, c 1 � 1.8044, c 2 � 1.3949, and v 2 � 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). e bifurcation diagram of q 2 with respect to v 1 , marked by a black box, is displayed in the enlargement. From Figure 2(a), when the adjustment speed of firm 1 satisfies v 1 < 0.4197, system (5) is at a period-1 stable state. In addition, v 1 � 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 v 1 is greater than 0.5371, system (5) will uncontrollably reach a disordered chaotic state.
Under same parameter conditions, it can be seen from Figure 2(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. erefore, 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. ere 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 v 1 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 to a � 11.1542, v 1 � 0.3625, and v 2 � 0.4517, then we can get two bifurcation diagrams with respect to costs of two firms, respectively. As shown in Figure 2(c), where c 2 � 1.8, when c 1 < 0.7755, system (5) is at the period-1 stable state. At c 1 � 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 c 1 � 1.8723, which means that there may be a phenomenon of coexistence. Similarly, we choose c 1 � 1.8 in Figure 2(d),then we can know that c 2 � 0.6330 is a bifurcation point, and once system (5) crosses this point, then a flip bifurcation will occur. e dynamic phenomena between 1.3500 < c 2 < 1.4180 are very complicated.
ere is also a "jump" phenomenon that exists near c 2 � 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 c 2 is higher than that with c 1 .
Another group of parameters are fixed as a � 11.6054, c 1 � 2.9154, and c 2 � 2.9639, and we choose adjustment speeds v 1 and v 2 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. ere are many fractal structures that can be seen, and the eye-shaped structure is symmetrically arranged on both sides of the ribbon. e dynamical behaviors, located in a white box of Figure 3(b), are enlarged in Figure 3(c).

Complexity 7
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

Complexity
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. at 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. erefore, 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. e single-parameter diagrams of the eye-shaped structure with respect to adjustment speeds are displayed in Figure 4, where the red curve represents q 1 of firm 1 and the blue curve represents q 2 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. e 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 periodhalving bifurcation occurs, when v 1 ∈ (0.345, 0.364). is 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 v 2 ∈ (0.31, 0.322), system (5) has a same bifurcation way as that in Figure 4(a). at is to say, the system will experience a period-doubling bifurcation at first and then enter into chaos, and another period-halving bifurcation arises.
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. e 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 v 1 ≈ 0.45 and v 2 ≈ 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 as a � 11.1542, c 1 � 1.8044, and c 2 � 1.3949 and choose the adjustment speeds v 1 and v 2 as bifurcation parameters; then, we can get a group of basins of attraction, as shown in Figure 5. e blue points represent attractors, and their stable region is the yellow area. e invariant cycles are denoted as red, and their stable region is marked as the dark sky blue area. e divergent area and chaotic area are denoted by dark blue. In Figure 5(a), the adjustment speeds of two firms are chosen as v 1 � 0.3605 and v 2 � 0.4517, respectively. ere 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 v 1 and v 2 increase in Figure 5(b), where v 1 � 0.3625 and v 2 � 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.
Furthermore, we assume from the 2D bifurcation diagram in Figure 1(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 v 1 and v 2 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 Figure 6, the values of parameters are chosen as a � 10.9788, c 1 � 0.0806, and c 2 � 0.2799, and we similarly choose v 1 and v 2 as bifurcation parameters. e shapes of these basins of attraction are very peculiar. e green area represents the divergent region, and the points in this region will tend to infinity and finally trap into the chaotic state. e 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. e flexible area of other pink attractors is the gray region. In addition, the phenomenon of multistability occurs when v 1 � 1.1940 and v 2 � 0.6810, and there are period-4 pink attractors coexisting with period-12 black attractors, which can be clearly seen in Figure 6(a). e boundary shape of basin of attraction is irregular, and there are many large "holes" in the upper and lower boundaries of basin of attraction. en, we keep parameter v 1 unchanged and increase the value of v 2 to 0.6840, and there still is a phenomenon of coexistence between distinct numbers of periodic attractors in Figure 6 at is, the period-12 attractors bifurcate to the period-24 attractors via a perioddoubling 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. ere is a pretty fractal structure on the head of basin of attraction in Figure 6(b), which is enlarged in Figure 6(c).
However, when the bifurcation parameters are varied as v 1 � 1.1671 and v 2 � 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 "contact bifurcation" may happen. In the following steps, we just change the value of adjustment speed v 2 . 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. erefore, the system will be more stable if the firm maintains a relatively small adjustment speed.

Global Dynamics and Critical Curves.
rough 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. e 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 � T(LC − 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. ose points satisfy the equation as follows: Obviously, the expression of equation Det(J) � 0 is very complex so that we cannot directly obtain analytic solutions with respect to q 1 and q 2 . 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 . It is apparent that these two critical curves separate the plane (q 1 , q 2 ) into three zones Z 0 , Z 2 , and Z 4 . at is to say, LC (b) separates zone Z 0 from zone Z 2 and LC (a) separates zone Z 2 from Z 4 . e points in zone Z 0 have no preimages, while the , where q * 1 satisfies the equation β 1 (1 − q 1 )e − q 1 � 2α 1 q 1 and q * 2 satisfies the equation

12
Complexity ese four preimages of origin point O � (0, 0) form the size and vertexes of the basin of attraction.
As shown in Figure 7, the parameters are fixed as a � 11.1542, c 1 � 1.8044, and c 2 � 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 (1) − 1 and O (2) − 1 are saddle points. Moreover, O (1) − 1 has a stable manifold along line OO (1) − 1 on q 1 -axis, and O (2) − 1 also has a stable manifold along line OO (2) − 1 on q 2 -axis. In Figure 7(a), when bifurcation parameters are chosen as v 1 � 0.3625 and v 2 � 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. ere are no "holes" on the boundary of basin of attraction because the critical curves, LC (a) and LC (b) , are totally inside the basin of attraction. e stable manifolds on the boundaries of basin, OO (1) − 1 and OO (2) − 1 , have not been destroyed. However, one of critical curves, LC (b) , almost contacts with the boundary line O (2) − 1 O (3) − 1 and already contacts with one of two period-2 points. With the adjustment speed v 2 increasing to 0.4257, the period-2 points bifurcate to the period-4 points via a flip bifurcation. Since the critical curve LC (b) once has crossed rank-1 preimage point O (2) − 1 , then the stability of stable manifold of O (2) − 1 will not exist anymore. It can be clearly seen from Figure 7(b) that the phenomenon "holes" occurs on the line OO (2) − 1 and its preimage O (1) is kind of bifurcation which has been mentioned in the previous section can be called a "contact bifurcation." As v 2 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). e 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, LC (b) . en, 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. ere are some phenomena of "lakes" that appear inside the basin of attraction. ese four pieces of chaotic attractors get bigger as the parameters v 1 and v 2 increase till they contact with critical curves LC (a) and LC (b) .
en, 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. ey may destroy the stability property of system (5) and may cause nonpredictability of tendency of game process.
A group of basins of attraction with interesting shapes of "holes" located on the boundary can be observed, when we fix the parameters as a � 11.1542, v 1 � 0.3625, and v 2 � 0.4517 and choose parameters c 1 and c 2 as bifurcation parameters. e 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. e origin point O � (0, 0) is still an unstable node, and its manifolds along both q 1 -axis and q 2 -axis which are unstable, while both O (1) − 1 and O (2) − 1 are saddle points so that O (1) − 1 has a stable manifolds along line OO (1) − 1 on q 1 -axis and O (2) − 1 also has a stable manifold along line OO (2) − 1 on q 2 -axis. As shown in Figure 8(a), when we fix the bifurcation parameters as c 1 � 2.8860 and c 2 � 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.
e blue area represents the divergent area. Another significant phenomenon is that a branch of critical curves, LC (b) , have already crossed both O (1) − 1 and O (2) − 1 , which means that two stable manifolds along line OO (1) − 1 on q 1 -axis and line OO (2) − 1 on q 2 -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 LC (b) has contacted with period-12 blue points, and then, it is highly possible that a "contact bifurcation" may happen. en, we fix the value of c 1 unchanged and increase the value of c 2 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 LC (b) 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 (q 1 , q 2 ). We also can clearly see that the "lakes" notably grow bigger.
As the value of c 2 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).
e "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.   14 Complexity 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 c 1 and c 2 similarly has significant influence on the stable property and structure of the basin of system (5). e 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.

Conclusions
In this paper, on the basis of considering the exponential demand function, the cost function is improved to quadratic cost. e firms considered in the market both have incomplete market information. erefore, they are boundedly rational. en, a discrete-time evolution model of the duopoly game with bounded rationality is built. e 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. rough 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.

Data Availability
No data were used to support the findings of this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.