#### 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.

#### 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 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. 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 firms (labeled by ) 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 is an exponential inverse demand function defined by the outputs of two firms, that is, , in which is the maximum selling price of the products in the market and 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 has a quadratic form, that is, , where is the marginal cost, and the relation between and satisfies . Hence, the profits of these two firms can be, respectively, written as follows:

Let , respectively, deduce the derivative with regard to the output , . Then, 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 is boundedly rational. The firms adjust their production strategies in next period in terms of their marginal profit in the previous period . It means that if the marginal profit in period is greater than zero (or less than zero), then the firm will increase (or decrease) its outputs . If the marginal profit in period is equal to zero, then the firm 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 with respect to its outputs can be described in the following form:in which represents the adjustment speed of the firm . Through 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, and , and we have that . We introduce an evolutionary operator and auxiliary variables as and , , where . Hence, a discrete-time dynamic duopoly Cournot game can be further simplified as the following:in which is an evolutionary operator.

According to the knowledge of nonlinear dynamics, let , , so we can know that the equilibrium point in system (5) must satisfy the conditions listed as follows:(i) or(ii)(iii) or(iv)

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 of system (5) can be obtained. For the rest of the boundary equilibrium points, we can draw the following conclusions:(1)The boundary equilibrium point can be obtained by condition (ii) and condition (iii), where satisfies the equation . If we rewrite this relationship as a function with regard to a variable , that is, , then we have and . So we can know that there is at least one point that satisfies . Finally, the existence of is verified.(2)Analogously, the boundary equilibrium point can be obtained by condition (i) and condition (iv), where satisfies the equation . If we rewrite this equation as a function with respect to a variable , that is, , then we will get two inequalities that and . So we can know that there is at least one point that satisfies . Finally, the existence of 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 point can be obtained through combining of condition (ii) with condition (iv), and satisfies the equation , and satisfies the equation , where .*

*Proof. *First, the equation can be obtained from condition (ii) and condition (iv), where and . Then, let , we can easily get two inequalities that and . Hence, we can know that there is at least one point in that satisfies the equation . Analogously, let , we have two inequalities that and . So there is at least one point that satisfies the equation , where . Finally, the existence of the Nash equilibrium point is verified.

In a real economic system, the equilibrium point represents the situation that these two firms both go bankrupt or get out of the market. The boundary equilibrium points and 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 of the Jacobian matrix. That is to say, if , then the equilibrium point is an unstable node. If , then the equilibrium point is a stable node. While if one of 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:

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 point is an unstable node.*

*Proof. *In system (5), the Jacobian matrix corresponding to the point can be written asIt is clear that Jacobian matrix (7) is a diagonal matrix, and the corresponding eigenvalues are and , respectively. Since the nonnegativity of and , we have that and . Therefore, the equilibrium point is an unstable node.

Proposition 3. *The boundary equilibrium point is unstable.*(1)*When , then the boundary equilibrium point is a saddle point*(2)*When , then the boundary equilibrium point is an unstable focus*

*Proof. *Apparently, the Jacobian matrix at equilibrium point is an upper triangular matrix, that is,The eigenvalues corresponding to this Jacobian matrix are and , respectively. Since , we have that and . Hence, and . However, the relationship between and 1 needs to be discussed by using classification as follows.(1)When , then we can deduce that . Since , then we can know that the moduli of and satisfy and . As a result, is a saddle point.(2)Conversely, when , then we have a conclusion that . Furthermore, combining the conclusions with , we can conclude that is an unstable focus.

Proposition 4. *The boundary equilibrium point is unstable.*(1)*When , the boundary equilibrium point is a saddle point*(2)*When , the boundary equilibrium point is an unstable focus*

*Proof. *Obviously, the Jacobian matrix at is a lower triangular matrix as the following:The eigenvalues corresponding to this Jacobian matrix are and , respectively. This is a similar situation like in Proposition 3, so we do same steps as in Proposition 3. Since belongs to interval , then we can get two inequalities and , and then, we already know that and . However, the relation between and 1 needs to be discussed by using classification as follows.(1)When , then we can deduce that . Then, we get the moduli of and satisfy and . As a result, is a saddle point.(2)When , then we can deduce that . Then, we get the moduli of and satisfy and . As a result, is an unstable focus.As for the unique Nash equilibrium of system (5), the corresponding Jacobian matrix at can be expressed asin which satisfies the equation and satisfies the equation . Therefore, the trace and the determinant of the Jacobian matrix are denoted as and , respectively. The specific expressions corresponding to and can be, respectively, written as follows:As we all know, the Nash equilibrium point is stable when the moduli of its eigenvalues corresponding to the Jacobian matrix 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 with the help of the Jury criterion.

Proposition 5. *The Nash equilibrium point is asymptotically stable if and only if one of the following six conditions holds:*(1)*When , the values of the parameters and satisfy and *(2)*When , the values of the parameters and satisfy and *(3)*When , the values of the parameters and satisfy and *(4)*When , the values of the parameters and satisfy and *(5)*When , the values of the parameters and satisfy and or and *(6)*When , the values of the parameters and satisfy and , or and **where , , , and .*

*Proof. *According to the Jury criterion, we can verify the local stability of . In other words, if all of conditions of the Jury criterion hold, then we can know that the Nash equilibrium is asymptotically stable. The specific expression of the Jury criterion is given as follows:Substituting the specific formulations of and into equation (12), then we will get equivalent expressions as follows:For simplicity, let , , , and . Since the conditions and , then we can know that , , , and . Then, the three conditions in equation (13) are equivalent toAccording to and , it can be judged that the condition 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 . Obviously, the condition is equivalent to the form of . In addition, in order to further analyze the local stability of and get the value ranges of the parameters and , the inequality can be divided into two parts for discussion. One part is equivalent to , that is, . Another part is equivalent to , that is, . Additionally, since , that is, , the auxiliary parameter involves the adjustment speed of the firm , and its value range can be divided into three cases as the following.(i^{′})When , .(ii^{′})When , .(iii′)When, , while , since the nonnegativity of the parameter , in this case. The above are the value ranges of the auxiliary parameter that need to be satisfied, when the equivalent condition in the Jury criterion holds. The situation is discussed as following, when the equivalent condition in the Jury criterion holds. Analogously, we will get an inequality , and then, it will be discussed by divided into two cases, that is, and . In other words, The discussion in terms of has two situations.(iv^{′})when , .(v^{′})when , ; we already know that , and since the nonnegativity of the parameter , so we consider that 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, , , , , and . That is to say, if one of the above six conditions in Proposition 4 holds, then 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 as , , and , and we choose adjustment speeds and 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 and 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 .

**(a)**

**(b)**

In Figure 2, the red curve denotes the outputs of firm 1, and the black curve denotes the outputs of firm 2. First, choose the adjustment speed as the bifurcating parameter, and the other parameters are fixed as , , , and . 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 with respect to , marked by a black box, is displayed in the enlargement. From Figure 2(a), when the adjustment speed of firm 1 satisfies , system (5) is at a period-1 stable state. In addition, 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 is greater than , system (5) will uncontrollably reach a disordered chaotic state.

**(a)**

**(b)**

**(c)**

**(d)**

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. 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 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 , , and , then we can get two bifurcation diagrams with respect to costs of two firms, respectively. As shown in Figure 2(c), where , when , system (5) is at the period-1 stable state. At , 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 , which means that there may be a phenomenon of coexistence. Similarly, we choose in Figure 2(d),then we can know that is a bifurcation point, and once system (5) crosses this point, then a flip bifurcation will occur. The dynamic phenomena between are very complicated. There is also a “jump” phenomenon that exists near . 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 is higher than that with .

Another group of parameters are fixed as , , and , and we choose adjustment speeds and 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**,

**, and**

*b***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.**

*c***(a)**

**(b)**

**(c)**

**(d)**

The single-parameter diagrams of the eye-shaped structure with respect to adjustment speeds are displayed in Figure 4, where the red curve represents of firm 1 and the blue curve represents 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 . 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 , 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.

**(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 and 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 , , and and choose the adjustment speeds and 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 and , 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 and increase in Figure 5(b), where and . 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.

**(a)**

**(b)**

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 and 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 , , and , and we similarly choose and 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 and , 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 unchanged and increase the value of 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).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

However, when the bifurcation parameters are varied as and , 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 . 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 (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 of critical curves is formed by a trajectory of points where the determinant of the Jacobian matrix disappears. In other words, and satisfy a relation that , that is, the points on will be mapped to through mapping . In this situation, the points on the need to satisfy a condition, that is, the determinant of Jacobian matrix (6) is equal to zero. Those points satisfy the equation as follows:

Obviously, the expression of equation is very complex so that we cannot directly obtain analytic solutions with respect to and . However, we can graphically display the shape of by numerical simulation. Two magenta branches of a hyperbola consist the , which is denoted as and , respectively. The critical curve is formed by two curves, marked by and , where and . It is apparent that these two critical curves separate the plane into three zones , , and . That is to say, separates zone from zone and separates zone from . The points in zone have no preimages, while the points in zone have two different rank-1 preimages, and the points in zone have four distinct rank-1 preimages. It is usual that the origin point always locates in the zone , which means that the origin point has four different rank-1 preimages. One of them is itself, . Two of them are located on coordinate axes, and , where satisfies the equation and satisfies the equation . The last one is , which is the intersection of preimages of two lines and . That is to say, line is the preimage of line , and line is the preimage of line . and simultaneously satisfy the formula as follows:

These four preimages of origin point form the size and vertexes of the basin of attraction.

As shown in Figure 7, the parameters are fixed as , , and . According to above-proposed Propositions 2–4, we can easily judge out that origin point is always an unstable node and both and are saddle points. Moreover, has a stable manifold along line on -axis, and also has a stable manifold along line on -axis. In Figure 7(a), when bifurcation parameters are chosen as and , 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, and , are totally inside the basin of attraction. The stable manifolds on the boundaries of basin, and , have not been destroyed. However, one of critical curves, , almost contacts with the boundary line and already contacts with one of two period-2 points. With the adjustment speed increasing to 0.4257, the period-2 points bifurcate to the period-4 points via a flip bifurcation. Since the critical curve once has crossed rank-1 preimage point , then the stability of stable manifold of will not exist anymore. It can be clearly seen from Figure 7(b) that the phenomenon “holes” occurs on the line and its preimage . This kind of bifurcation which has been mentioned in the previous section can be called a “contact bifurcation.” As 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, . 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 and increase till they contact with critical curves and . 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.

**(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 as , , and and choose parameters and 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 is still an unstable node, and its manifolds along both -axis and -axis which are unstable, while both and are saddle points so that has a stable manifolds along line on -axis and also has a stable manifold along line on -axis. As shown in Figure 8(a), when we fix the bifurcation parameters as and , 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, , have already crossed both and , which means that two stable manifolds along line on -axis and line on -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 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 unchanged and increase the value of 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 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 . We also can clearly see that the “lakes” notably grow bigger.

**(a)**

**(b)**

**(c)**

**(d)**

As the value of 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 and 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.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the Foundation of A Hundred Youth Training Program of Lanzhou Jiaotong University.