Complex Dynamical Behavior of a Bounded Rational Duopoly Game with Consumer Surplus

In this chapter, we assume that two bounded rational firms not only pursue profit maximization but also take consumer surplus into account, so the objections of all the firms are combination of their profits and the consumer surplus. And then a dynamical duopoly Cournot model with bounded rationality is established. The existence and stability of the boundary equilibrium points and the Nash equilibrium of the model are discussed, respectively. And then the stability condition of the Nash equilibrium is given. The complex dynamical behavior of the system varies with parameters in the parameter space is studied by using the so-called 2D bifurcation diagram. The coexistence of multiple attractors is discussed through analyzing basins of attraction. It is found that not only two attractors can coexist, but also three or even four attractors may coexist in the established model. Then, the topological structure of basins of attraction and the global dynamics of the system are discussed through invertible map, critical curve and transverse Lyapunov exponent. At last, the synchronization phenomenon of the built model is studied.


Introduction
Oligopoly is a market between perfect monopoly and perfect competition [1]. With the application of chaos theory and nonlinear dynamic system into oligopoly models, the static game evolves into a dynamic game. Especially in recent years, with the rapid development of computer technology, a powerful tool has been provided for dealing with the complex nonlinear problems. And hence, the economists and the mathematicians can simulate the complex dynamical behavior of oligopoly market by using computer technology. Recently, a large number of scholars have improved the oligopoly models, and introduced bounded rationality (see [2,3]), incomplete information [4], time delay [5], market entering and entering barriers [6], differentiated products [7] and other factors [8,9] into the classical oligopoly models, and the bifurcation and chaos phenomenon were founded in the process of repeated game.
However, all of the above discussions are mainly based on private enterprises, which pursuit the maximization of their own profits. In fact, the public ownership enterprises, which always aim at maximizing social welfare, and mixed ownership enterprises, which always aim at maximizing the weighted average of the social welfare and their own profits, are also widespread in the real economic environment. De Fraja and Delbono [10] found that the social welfare might be higher when a public ownership enterprise is a profit-maximizer rather than a social-welfaremaximizer. Matsumura [11] proposed that the social welfare could be improved through partial privatization of public enterprises. The research of Fujiwara [12] suggested that partial privatized public enterprises are more efficient than private enterprises. Elsadany and Awad [13] explored the complex dynamical behavior of competition between two partial public enterprises under the assumption of bounded rationality. However, the global dynamical behavior and synchronization behavior of semi-public enterprises, which corporate social responsibility into their objectives, are rarely studied. In this chapter, the occurrence of synchronization, the coexistence of attractors and the global dynamic of a duopoly game corporate consumer surplus are mainly discussed.

The model
Considering a duopolistic market where two firms produce homogeneous goods. In order to study the long-term behaviors of the duopoly market with quantity competition, we briefly present the economic setup leading to the final model in this chapter. The price and quantity of product of firm i are given by p i and q i respectively, with i ¼ 1, 2. We also assume the existence of a continuum of identical consumers which have preferences toward q 1 and q 2 .
Following Dixit [14] and Singh and Vives [15], we suppose that the utility function used in this chapter is quadratic and can be given by, where q 1 , q 2 are the quantity of goods produced by firm 1 and firm 2, respectively. a > 0 represents the maximum price of a unit's commodity, b > 0 represents the amount of its price decreases when the price of the product increases by one unit.
Suppose that the budget constraint of consumer is, where p 1 and p 2 denote the prices of goods produced by firm 1 and firm 2, respectively. And M denotes the budget of the consumers on the product. The utility function of consumers is maximized under the budget constraint, and then the inverse demand function of the two firms can be obtained as, This chapter discusses homogenous products, so here it is assumed that all these two players have the identical marginal cost. Therefore, the cost function of firm 1 and firm 2 are same and can be given by, C q ð Þ ¼ cq, where c > 0 denotes the marginal cost of the goods and a > c always holds. Then, the profits of firm i, i ¼ 1, 2 can be obtained as follows, In the real market, there are a lot of firms, who not only pursue their own profits but also take corporate social responsibility into account. A large number of empirical studies have shown how the introduced corporate social responsibility affects firm's performance, where we interpret corporate social responsibility as either consumer surplus (for short CS) or social welfare (for short SW). In this chapter we take CS into account to analyze which firms have an incentive to exhibit corporate social responsibility as a means of maximizing their profits in a Cournot competition. Based on the above assumptions and the definition of consumer surplus, CS can be written as, where p ∈ p; a ð Þ is the price variable. According to the above assumptions, the objective function of the firm i can be given as, where α i represents the weight of the consumer surplus in the objective function of firm i, and 0 ≤ α i ≤ 1 always holds. By substituting (4) and (5) into (6), the objection function of firm i can be given by, And the first-order condition of the objection function (7) is given as, It is now significant to specify the information set of both players regarding the objection functions, to determine the behaviors of the players with the change of time. We assume a discrete time t ∈ Z þ ð Þdynamic setting, where two firms with bounded rationality make decisions at the same time. That is, all firms do not have complete knowledge of their competitors' decisions and the market demands. So it can only use the local estimation of the steepest slope of the objection function at period t to determine the output at period t þ 1. By following Bischi et al. [16] and Fanti et al. [17], the adjustment mechanism of quantities with the change of time of firm i can be obtained as, where v i > 0, i ¼ 1, 2 is an adjustment parameter of firm i. The firm i will increase its output at period t þ 1, if ∂O i t ð Þ ∂q i t ð Þ > 0. But the firm i will reduce its output at period t þ 1, if ∂O i t ð Þ ∂q i t ð Þ > 0. By substituting (8) into (9), we can get a two-dimensional map as, Since the output of a firm cannot be negative, the initial conditions of map T belong to (10), the fixed points of the system are obtained. Besides the trivial equilibrium E 0 ¼ 0; 0 ð Þ, system (10) admits the following non-trivial fixed points (boundary equilibrium points), ; 0 and the only Nash equilibrium is

Stability properties
The local stability analyses of system (10) near the fixed points are too difficult to carry on. For the sake of analyzing the local stability of the system, we firstly let α 1 ¼ α 2 ¼ α in system (10). And the Jacobian matrix of map T at any fixed point q 1 ; q 2 À Á can be given as, Then all the equilibrium points are substituted into the Jacobian matrix (12). According to the characteristic values of the Jacobian matrix evaluated at each equilibrium, the type and stability of the equilibrium can be analyzed and the following results can be obtained. Proposition 1. The equilibrium point E 0 is always an unstable node. Proof. It is clear that the Jacobian matrix of map T, evaluated at the boundary equilibrium point E 0 can be written as, The eigenvalues of J E 0 ð Þ are given by Proof. By substituting the equilibrium E 1 into (11), the Jacobian matrix of map T evaluated at the boundary equilibrium point E 1 can be written as, The eigenvalues of J E 1 ð Þ are given by Since (11) holds and v 1 > 0, then λ 1 > 1 always meets. According to (11) and v 2 > 0, we can deduce that λ 2 j j< 1 if and only if v 2 < . For the purpose of research of the local stability near the Nash equilibrium, we should compute the Jacobian matrix evaluated at the Nash equilibrium E * as, It can be seen that the form of the Jacobian matrix is so complex. In order to simplify the calculation, let Then the trace and the determinant of the Jacobian matrix evaluated at the Nash equilibrium E * can be given as, According to Jury condition, if we substitute the specific mathematical expressions of q * 1 , q * 2 into the above two equations, then the following set of inequalities can be gotten through a complex calculation, Since all the equilibrium points are non-negative when the parameters meet 0 ≤ α < 3 5 , a > c and v i > 0, i ¼ 1; 2 ð Þ. So we can get A > 0, B < 0, B þ D < 0 and B À D < 0, then the set of inequalities (14) are equivalent to Then the stability region of the Nash equilibrium can be obtained by substituting A, B and D into inequalities (15), which are given as, The stability condition of the Nash equilibrium gives a parameters region, in which the Nash equilibrium is always stable. For the sake of better analysis of the stability of the Nash equilibrium under different set of parameters, a useful tool called "two-dimensional bifurcation diagram" (also called 2-D bifurcation diagram) is employed. From (16), we can find that the stability region of the Nash equilibrium is related to the difference of parameters a and c, that is a À c. So we only discuss the values of a À c, rather than the values of a and c in the rest of this chapter. Figure 1 is a two-dimensional bifurcation diagram of system (10) with a set of fixed parameter a À c, b and α. Figure 1a is a two-dimensional bifurcation diagram when a À c ¼ 2:1, α ¼ 0:50608821 and b ¼ 2:00232504. We can observe that there are two different routes to chaos, when the parameters are chosen as this set of parameters. The system enters chaos through flip bifurcation when the parameter v 1 ; v 2 ð Þpasses through green, yellow and light green from the brown region to the black region. It means that when the firms change their speed of adjustment according to the path, a periodic fluctuation of system (10) will happen. That is, the period motion will increase exponentially until it enters chaos. And the chaotic behavior of this system can be understood as the confusion of the market competition, and one of the two firms may be even out of the market with an increasing speed of adjustment. But if the parameter passes through the green region from the brown region to the black region, the system will first undergo a flip bifurcation, and then enters quasi-period motion through a Neimark-Sacker bifurcation. The system enters quasi-period from period-2 when firms determine their speed of adjustment along this path. Figure 1b is a partial amplification of Figure 1a, and we can observe that there are many scattered points of different colors, which is caused by the coexistence of multi-attractors with different period. Figure 2a shows the coexistence of attractors with the parameters chosen as v 1 ¼ 2:44, and v 2 ¼ 2:45, where the scatter points are shown in Figure 1b. we can observe a period-6 cycle coexisting with a period-4 cycle. Figure 1c is the twodimensional bifurcation diagram, where the fixed parameters are given as a À c ¼ 0:88017028, α ¼ 0:27445462 and b ¼ 0:52714274. At this set of parameters, the system enters chaos through a flip bifurcation. Figure 1d is a partial enlargement of Figure 1c. Similarly, the parameter space of Figure 1d is also chosen according to the area with scattered points of Figure 1c. Figure 2b shows the coexisted attractors and their basins of attraction at this set of parameters. We can observe that there are three attractors coexisting. Figure 3 shows a series of two-dimensional bifurcation diagram under different parameters. It shows a very beautiful gallery, from which we can enjoy the system (10) with full complex dynamics phenomenon. We can observe from Figure 3 that the difference between the maximum price of the unit product a and the marginal cost c affects the size of the stable region. The weight of the consumer surplus α affects the shape of the two-dimensional bifurcation diagram, while the parameter b almost has hardly effect on the two-dimensional bifurcation of the system. Therefore, the game can be balanced more quickly by reducing the difference a À c. In Figure 3a, it can be observed that the chaotic area surrounded by the period-4 area is like a "hand" and the 8-period area is like a small "bottle" raised by the beautiful hand. It is observed that the shape of Figure 3b is similar to Figure 3a due to a tiny adjustment of parameter α. Since the difference between the parameters a and c is reduced, the stability area of the Nash equilibrium becomes larger. Figure 3c is like a "volcanic eruption." It can be observed that there is an inward cave in the diagonal that is like a "crater." As the parameters vary, the hole in Figure 3c continues to sink inward and become larger and larger. From Figure 3d we can see that the period-4 arrives at a quasi-periodic motion directly. Therefore, we should change the weight of consumer surplus α slightly in order to maintain the market fluctuations not fierceness.

Global dynamics and synchronization
The type and stability of the equilibrium points have been analyzed as above. And the boundary equilibrium E 1 and the boundary equilibrium E 2 are in symmetrical positions with respect to the main diagonal line It is also clear that the unique Nash equilibrium E * of system (10) is located on the main diagonal Δ. So we mainly study the dynamical behavior of the system on the diagonal. We choose the initial conditions near the diagonal, and the phenomenon via finite iteration back to the diagonal is called synchronization. The synchronization of chaotic systems was quite interesting and unexpected. In fact, due to the nonlinear system usually has sensitive dependence on initial conditions, a property which implies that the slightly change of initial conditions will lead to an exponential difference between the trajectories of two identical systems, making it impossible for two separated and even identical systems to synchronize. Therefore, the small coupling between two chaotic oscillators makes the system asymptotically converge to same trajectory, which is worth studying.
Subsequently, we assume that both firms have the same speed of adjustment. It means that the latter discussion is based on v 1 ¼ v 2 ¼ v and α 1 ¼ α 2 ¼ α. In this case the two players are identical; the system T can then be rewritten as follows, It can be proved that the map T 0 has symmetry property, i.e., there exists a map S : q 1 ; q 2 À Á ! q 2 ; q 1 À Á , which makes T 0 ∘S ¼ S∘T 0 . The symmetry property of the map T 0 implies that the diagonal Δ is a one-dimensional sub-manifold of system (17), i.e., T 0 Δ ð Þ ⊆ Δ. However, the phenomenon of synchronization occurs when the diagonal Δ is invariant one-dimensional submanifold of the system (17). Therefore, the phenomenon of synchronization of the system can be analyzed by studying the invariant set. We can also use critical curve and noninvertible map to describe the global dynamical behaviors of a 2-dimensional map.

Critical curve and noninvertible map
We divide the discrete dynamical system into invertible and noninvertible. The invertible discrete dynamical system refers that an image q 1 0 ; q 2 0 À Á of the map T 0 is correspond to the only preimage q 1 ; q 2 À Á . The noninvertible discrete dynamical system implies that the map T 0 is multi-valued, i.e., the image of T 0 has more than one preimages. In a noninvertible discrete dynamical system, the curve that divides the phase space into regions with a different number of rank-1 preimages is called critical curve, denoted by LC. And the regions can be represented by Z i , i ∈ N ð Þ. For example, a point belonging to area Z 0 has no preimage and a point belonging to area Z 2 has two preimages. Let us denote the rank-1 preimages of critical curve LC under map T 0 as LC À1 . The set LC is the 2-dimensional generalization of the critical value or local extremes of 1-dimensional noninvertible map. Its preimages LC À1 are corresponding to local extreme point (critical point) in the one-dimensional noninvertible map. Since the map (17) is a continuously differentiable map, LC À1 belongs to the locus of points where the Jacobian determinant of T 0 vanishes, i.e., LC À1 ¼ q 1 ; q 2 À Á ∈ R 2 detDT 0 ¼ 0 j g È . In this case, curve LC À1 can be determined by the following equation, The noninvertible properties play a significant role in analyzing the global behavior of a nonlinear discrete dynamical model. So the critical curve is a powerful tool for us to study these complex structures. Using the segment of critical curve as well as their preimages of any rank, and we will get the boundary of the basins of attraction as shown in Figure 10f.

Invariant sets
The dynamics of the system on the diagonal is studied by analyzing the invariant sets. Firstly, we can prove that the coordinates are invariant sets of map T 0 . Let q 2 t ð Þ ¼ 0, then we can obtain q 2 t þ 1 ð Þ¼0, and the first equation of (17) can be rewritten as, It is easy to verify that the dynamics on the axis q 2 is also controlled by the map (19). It means that the system (17) can be regarded as a 1-dimensional map at the coordinate axes. The map (19) is topologically conjugate to the standard logistic Þthrough a linear transformation, which is given as, and the parameter ω can be presented as ω Thus the nonlinear dynamics of system (17) on the invariant axes can be analyzed through the standard logistic map.
It can also be proved that the diagonal Δ is an invariant set of system (17), i.e., the trajectory starting from the diagonal Δ will stay forever on it. Therefore, the dynamical behavior of system (17) can be analyzed through the map T 0 which is restricted to the diagonal. If we let q 1 ¼ q 2 ¼ q, then the dynamics generated by T 0 on the diagonal Δ can be analyzed through the following map, Similarly, through the following linear transformation we can also prove that the map (21) is topologically conjugate to the standard logistic map y t þ 1 ð Þ¼μy t ð Þ 1 À y t ð Þ ð Þ , where Through the standard logistic map, we can easily analyze the dynamical behavior of the two-dimensional map T 0 on the diagonal Δ. Under this situation, the Nash equilibrium E * of the system (10) is identical with the fixed point of map T Δ 0 . Since μ ¼ 1 þ v 1 À α ð Þ a À c ð Þ, we take different values of the bifurcation parameter μ of the logistic map, and Figure 4a gives different bifurcation curves of the system on the parameter plane α; v ð Þ. The flip bifurcation occurs when the system parameter v equals v ¼ 2 1Àα ð Þ aÀc ð Þ , and the Nash equilibrium E * loses its stability and forms a period-2 cycle around E * . At v ¼ Þ aÀc ð Þ , the period-2 cycle generates a period-4 cycle after a flip bifurcation. When μ≈3:5699, the standard period doubling cascade ends and the system enters chaos. When v > 3 1Àα ð Þ aÀc ð Þ , the general trajectory of the map T Δ 0 is divergent. As shown in Figure 4b, which the parameters is the same as Figure 4a, a twodimensional bifurcation diagram of the system with v and α is obtained. Since it has been proved that the map T Δ 0 is topologically conjugate to the logistic map, we can find that the bifurcation curves of the two graphs are the same. In Figure 4a, the curve C 1 is correspond to the equation v ¼ 2 1Àα ð Þ aÀc ð Þ , and the region below it represents the set of points of v and α at 1 < μ < 3. In this region, the fix point is stable. That is, the synchronization trajectory converges to the Nash equilibrium point. In Figure 4b, it corresponds to the period-1 region below the green region. When the point above the curve C 2 passes through the curve C 0 in Figure 4a, the system goes into chaos through a period doubling cascade. In Figure 4a, the curve C ∞ is correspond to the equation v ¼  Figure 4b. Through the above analysis, the following proposition can be derived, Proposition 3. If we let v 1 ¼ v 2 ¼ v, the parameters a, b and c are fixed for system (17). Then, a threshold α 0 ¼ 1 À 3 v aÀc ð Þ of the weight of consumer surplus α or a threshold v 0 ¼ 3 1Àα ð Þ aÀc ð Þ of the speed of adjustment v does exist such that synchronized trajectories of the system (17) are divergent when ∀α ∈ 0; α 0 ½ Þor ∀v ∈ v 0 ; þ∞ ð Þ. In order to analyze the effect of any slight perturbation of one parameter on the system, we study the transverse stability of an attractor A of map T 0 . And the Jacobian matrix of map T 0 on the diagonal can be obtained as follow, Then, the characteristic values of the Jacobian matrix J q; q ð Þ evaluated at any point on the diagonal are given by, where the corresponding eigenvectors are 1; 1 ð Þ and 1; À1 ð Þ, respectively. And the eigenvalue λ k is related to the invariant manifolds on the diagonal.
It is assumed that a period-k cycle q 1 ð Þ; q 1 ð Þ ð Þ ; q 2 ð Þ; q 2 ð Þ ð Þ ; ⋯; q k ð Þ; q k ð Þ ð Þ f g embedded into the invariant set Δ of the map T 0 is correspond to the cycle q 1 ð Þ; q 2 ð Þ; ⋯; q k ð Þ f g of the map T Δ 0 when the synchronized phenomenon occurs, the two multipliers are given as, Since the stability conditions of the period-k cycle on the diagonal Δ of system (17) is same with the one-dimensional map T Δ 0 , here we only study the transverse stability of the one-dimensional map T Δ 0 . Under this situation, the transverse eigenvalue evaluated at the Nash equilibrium point E * is given by Through Eq. (27), we can draw the following conclusions directly. That is, when all the parameters satisfy 0 < v a À c ð Þ 1 À α ð Þ 2 þ 5α < 3, the Nash equilibrium E * is transversely attractive.
As we know that an attractor A of T 0 is asymptotically stable if and only if all the trajectories that belong to attractor A are transversely attractive. To study the stability of the attractor, we can calculate its transverse Lyapunov exponent as, where q 0 ð Þ ∈ A and q i ð Þ is a generic trajectory generated by the map T Δ . If the initial condition q 0 ð Þ belongs to a period-k cycle, then Λ ⊥ ¼ ln λ k ⊥ . In this case, if Λ ⊥ < 0, then the period-k cycle is transversely stable. When the initial condition q 0 ð Þ belongs to a generic aperiodic trajectory embedded in the chaotic attractors, then the transverse Lyapunov exponent Λ ⊥ is the natural transverse Lyapunov exponent Λ nat ⊥ . Since many unstable cycles along the diagonal are embedded in the chaotic attractor A, a spectrum of transverse Lyapunov exponents can be determined by the inequality If all cycles embedded in A are transversely stable (Λ max ⊥ < 0) then A is asymptotically stable in the Lyapunov sense. If some cycles embedded in the chaotic attractor A are transversely unstable (Λ max ⊥ > 0 and Λ nat ⊥ < 0) then A is not stable in the Lyapunov sense, but it is a stable Milnor attractor. So we can look for the Milnor attractors by transverse Lyapunov exponents. Figure 5 gives the natural transverse Lyapunov exponent and the bifurcation diagram with the fixed parameter v when a À c ¼ 5:15, b ¼ 0:1911895 and α ¼ 0:4. Later, we will exhibit the attractors and their basins of attraction corresponding to different values v under this set of parameters, and analyze the changes of attractors and their basins of attraction when the parameter v varies.

Global bifurcation and basins of attraction
A closed invariant set A is a attractor which means that it is asymptotically stable, i.e., a neighborhood U of A does exist such that T 0 U ð Þ⊆U and We also define a asymptotically stable invariant set as attractor. A basin of attraction may contain one or more attractors that may coexist with a set of repel points that produce either intermittent chaos or a blurry boundary. The basin of attraction of attractor A is the set of those initial conditions that cause the trajectory to converge to A, i.e., For the sake of analyzing the topological structure of the basin of attraction B A ð Þ, we study the boundary of B A ð Þ firstly. Suppose that the map T 0 has a unique attractor A at finite distance, let ∂B A ð Þ be the boundary of the basin B A ð Þ, then it is also the boundary of the basin of infinity B ∞ ð Þ generated by unbounded trajectories. Firstly, we take the dynamics of system (17) into account and restrict it to the invariant axis. (17), we can obtain the bounded trajectories along the invariant axes, where 0 i À1 , i ¼ 1; 2 ð Þis the rank-1 preimage of the origin. It has been obtained previously that the dynamical behavior of system (17) on the coordinate axis is governed by the map (19), so that 0 i À1 , i ¼ 1; 2 ð Þcan be computed by the following algebraic system The result is given as, Since ε 1 and ε 2 are the segments of the boundary ∂B A ð Þ, and ∂B A ð Þ is also the boundary of the basin of infinity B ∞ ð Þ, their rank-k preimages T 0 Àk ð Þ ε i ð Þ, i ¼ 1; 2 ð Þ also belong to ∂B A ð Þ. We can compute the rank-1 preimages of a point according to the algebraic system as follows We can easily obtain the rank-1 preimages of the origin, which are , there are itself and O 3 À1 besides, i.e., Through the discussion above, we can get the following propositions, Proposition 4. Let 1 < v 1 À α ð Þ a À c ð Þ< 3 and ε i ¼ 0; 0 i À1 Â Ã , i ¼ 1, 2 be the segments of the coordinate axes q i , i ¼ 1, 2, then we can obtain the boundary of B A ð Þ as follow, Basins of attraction may be connected or not. The connected basins of attraction are divided into simple connected and complex one, and the complex connected basins of attraction means the existence of holes. If A is a connected attractor, the direct basin of attraction D 0 of A is the largest connected area of the entire attractor domain D containing A. The system (17) has the coexistence of attractors in a set of given parameters, the basin of attraction D refers to the union of the domain of attraction of all attractors in such a situation. Figure 6 shows the coexistence of attractors and their basins of attraction for given parameters a À c ¼ 5:3, b ¼ 0:234 and v ¼ 0:85. In Figure 6a, the parameter α is chosen as α ¼ 0:4, there are two attractors coexisting, one is a Milnor attractor A located on the diagonal and the other consisting of 4-piece chaos attractor is in symmetrical positions with respect to the diagonal, i.e., F ¼ ∪ 4 i¼1 F i . The basin of attraction is composed of the union of the attractive domain of two attractors. The attractive domain of the Milnor attractor A is the complex connected set, and the attraction domain of the attractor F is non-connected set. And the boundary of the 4-piece chaos attractor is just contact with the critical curve, and it is because of this contact that the system undergoes a global bifurcation. Figure 6b is the attractor and the attractive basin at α ¼ 0:387275 after the global bifurcation occurred in the Figure 6a. We can find 4 attractors coexisted in this figure the attractor F in Figure 6a undergoes a global bifurcation and turns into 3 period-4 cycles, and the attractor A is also a Milnor attractor. The attraction domain of the period-4 cycle is composed of some complex connected sets being in symmetrical positions with respect to the diagonal. As is shown in Figure 6b, there are many holes in the attracting domain, but it is a non-connected set. The basin of attraction of the Milnor attractor lying on the diagonal is still a complex connected set.
We have analyzed the global bifurcations that occur when the attractor's boundary contact to the critical curve, and we discuss the global bifurcation when the attractor contacts to the boundary of its basin of attraction. We also denote global bifurcation as "boundary crisis," the attractor is destroyed when it contacts to its basin of attraction. Figure 7 shows the coexistence of attractors and their basins of attraction corresponding to the parameter a À c, when the parameters are chosen as α ¼ 0:5023335, b ¼ 0:4 and v 1 ¼ v 2 ¼ 0:85. We can see that as the difference between the maximum price a of a unit commodity and the marginal cost c increases, a period-4 cycle turns into 4-piece chaos attractor being in symmetrical positions with respect to the diagonal and finally merges into 2-piece chaos attractor. However, the 2-piece attractor being in symmetrical positions with respect to the diagonal grows larger as the parameter a À c increasing, until it contacts to its basin's boundary, and eventually occurs a global bifurcation, causing itself and its basin are destroyed until it disappears. We can also see its "ghost" in Figure 7d. This means that trajectories of the initial conditions that belong to the basin of attraction spend a long number of steps in the region occupied by the former attractor before converging to the other attractor. Figure 8 is the bifurcation diagram of the system at this set of parameters, and the bifurcation parameter is chosen as a À c. Figures 6 and 7 give two different global bifurcations, such bifurcations which can be restored clearly by numerical simulation method only. With the set of parameters in Figure 9 being identical to Figure 4, we select different speed of adjustment to analyze the change of attractors and their basins of attraction. We can observe that as the speed of adjustment changes from 0.79 to 0.9, the period-4 cycle being in symmetrical positions with respect to the diagonal of the system generates smooth limit cycle via a Neimark-Sacker bifurcation, as shown in Figure 9b, and the limit cycle becomes non-smooth gradually, and finally forms four-piece chaotic attractor, as shown in Figure 9d. The basin of attraction shrinks as the speed of adjustment v increasing. It is implied that when both firms choose a lower speed of Figure 7. Basin of attractions for parameters b ¼ 0.4, α ¼ 0.502335 and v ¼ 0.85 (a) a À c ¼ 6, a two-piece chaotic attractor coexists with a period-4 cycle, (b) a À c ¼ 6.1, the 4-piece chaotic attractor coexist with a two-piece chaotic attractor, (c) a À c ¼ 6.1896, the 2-piece chaotic attractor formed by the a 4-piece chaotic attractor coexist with a two-piece chaotic attractor, and the many holes created by the global bifurcation and (d) a À c ¼ 6.26, a two-piece chaotic has a contact with its basin's boundary, and it is destroyed. adjustment, they can reach the balance easily in the game. However, the period-2 cycle embedded in the diagonal becomes a period-4 cycle, period-8 cycle, etc. That is, a flip-bifurcation happens. And finally a Milnor attractor forms with the increasing speed of adjustment. Its basin of attraction increases with the increasing speed of adjustment gradually. In the bifurcation diagram of Figure 5b, we can observe the process of entire bifurcation process.

Synchronization
In this section we study the formation mechanism of the synchronization trajectories. The trajectories starting from different initial conditions return to the diagonal eventually, i.e., q 1 0 ð Þ 6 ¼ q 2 0 ð Þ. A t * does exist such that q 1 t ð Þ ¼ q 2 t ð Þ when t > t * , and we define such trajectories as synchronization. However, when the diagonal Δ is an invariant sub-manifold, synchronized dynamics occur. We have proved that the map T 0 can be obtained by two identical one-dimensional coupling maps, and the synchronization trajectory can be controlled by a map T Δ 0 which is topologically conjugate to the standard logistic map. When we choose the Figure 10.
Parameter values are chosen as v ¼ 1, a À c ¼ 5 and b ¼ 0.234. (a) Four-piece Milnor attractor of system T belonging to the diagonal for α ¼ 0.48365, (b) the displacement q 1 À q 2 versus time for the same parameters as in (a), (c) α ¼ 0.485092, a 16-cyclic chaotic attractor is in symmetrical positions with respect to the diagonal, (d) α ¼ 0.437609, a trajectory in the phase space (q 1 , q 2 ) whose transient part is out of diagonal that synchronizes along the Milnor attractor in the long run, (e) α ¼ 0.476955, a two-cyclic chaotic attractor coexists with a period-2 cycle and (f) boundary of the chaotic area obtained by ∂A ¼ ∪ 6 k¼1 T 0 k ð Þ Γ ð Þ.
parameters as v ¼ 1, a À c ¼ 5 and b ¼ 0:234, the weight α varies, and we can observe that the dynamic behavior of system is controlled by the attractor on the diagonal in Figure 10. When α ¼ 0:4837015, we can observe a Milnor attractor in Figure 10a. This means that cycle embedded in the diagonal are transversely unstable and blowout phenomenon occurs when the trajectory is near diagonal. The trajectory converges to the unique Milnor attractor embedded in the diagonal after experiencing a long transient. Figure 10b shows that the evolution of q 1 À q 2 versus time and synchronization is observed after a long transient. This is a typical on-off intermittency phenomenon. We can observe 16-piece chaos attractor being in symmetrical positions with respect to the diagonal in Figure 10c, when α increases to 0:485092. Figure 10d shows a chaotic attractor when α decreases to 0:437609, then synchronization occurs. As shown in Figure 10f, we adopt the trial-and-error method, with suitable part of LC À1 taken as the starting part of Γ ¼ A∩LC À1 to obtain the boundary of the chaotic attractor A and the entire basin of attraction in Figure 10e, i.e., the boundary of the chaotic attractor A is ∂A ¼ ∪ 6 k¼1 T 0 k ð Þ Γ ð Þ.

Conclusion
In this chapter, the nonlinear dynamics of a Cournot duopoly game with bounded rationality is investigated. Unlike the existing literature, we suppose that the two firms not only pursue profit maximization but also take consumer surplus into account. Meanwhile, the objection of firms is supposed as the weighted sum of profit and consumer surplus. Based on the theory of gradient adjustment, all the firms adjust the output of next period according to the estimation of "marginal goal." The existence and stability of fixed points are analyzed. It is found that the boundary equilibrium point is always unstable, no matter what the parameters of the system are satisfied. At the same time, with the two-dimensional bifurcation diagram as the tool, the stability of the Nash equilibrium is analyzed. We found that the Nash equilibrium will lose its stability when the speed of adjustment of firms is too large, which maybe lead the market into chaos. The stability region of the Nash equilibrium will be only affected by the weight of consumer surplus. And the parameters a À c and b have hardly effect on the stability region of the Nash equilibrium. Meanwhile, we found that the two-dimensional bifurcation diagram have a beautiful fractal structure, but there are also many scattered points which is due to the coexistence of multiple attractors of the system through numerical simulation. By selecting corresponding parameters in the two-dimensional bifurcation diagram with scattered points, we draw the corresponding basin of attraction, and found the model not only has two attractors coexistence phenomenon, but also has 3, or even 4 attractors coexistence phenomenon.
Moreover, with the theory of invertible mapping and the critical curves of the system, the topological structure of basin of attraction is analyzed. By calculating the transverse Lyapunov exponent, the weak chaotic attractor of the system in the sense of Milnor is found, and the synchronization of the system is further studied. If we fix the other parameters of the system, and only change the weight of the firm to the consumer surplus, we can find on-off intermittency phenomenon and synchronization phenomenon. With the increasing of α, the synchronization phenomenon is vanished and a 16-piece chaotic attractor being in symmetrical position with respect to the main diagonal is produced. Under another set of parameters, and the parameter α is chosen as the bifurcation parameter. Through numerical simulation, it can be found that when the critical curve contact with the boundary of the basin, a global bifurcation is obtained. The global bifurcation makes the basin of attraction of the attractor non connected. In addition, if we fixed parameters of the system, and change the values of the parameters a À c only, we find another global bifurcation called "boundary crisis," i.e., when the attractor contact with its boundary of the basin of attraction, one of the attractors and its basin of attraction will be destroyed.