# On Complex Dynamics of Differentiated Products: Cournot Duopoly Model under Average Profit Maximization

Authors: S.S. Askar
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

In this paper, some important dynamic characteristics such as multistability and synchronization phenomena are investigated for a game of an economic Cournot duopoly whose time evolution is received by the iteration of a noninvertible map in the plane. In the asymmetric case, the equilibrium points of game’s map are calculated, and their stability conditions are obtained. The obtained results show that the Nash equilibrium point loses its stability through flip bifurcation. Under some restrictions, the map’s coordinate axes form an invariant manifold, and hence their dynamics are studied based on a one-dimensional discrete dynamic map. In the symmetric case where both firms are identical, the map has the property of symmetry, and this implies that the diagonalq1=q2 forms an invariant manifold and therefore synchronization phenomena occur. Global analysis of the behavior of the noninvertible map is carried out through studying critical manifolds of the map that categorize it as Z4−Z2−Z0 type. Furthermore, global bifurcation of the basins of attraction is confirmed through contact between the critical curves and the boundaries of escaping domain.


---

## Body

## 1. Introduction

Dynamic economic games whose evolutions has been depended on discrete time maps have attracted many researchers due to the interesting dynamic behaviors and bifurcations types occurring in such maps. Analyzing the dynamic characteristics of such maps has been dated at least to [1] by Puu who passed away in 2020 after enriching this direction with many useful books and papers. In the current paper, we focus on a special type of these maps that are known as duopoly maps. In duopoly game, the market possesses only two competing firms (or players) whose target is seeking the optimal quantities in the case of Cournot by maximizing their profits, or seeking the optimal prices in the case of Bertrand. Cournot, the French economist who started the core of such kind of games in 1838, introduced a duopoly model adopting quantities as decision variables. Later in 1880s and based on Cournot model, Bertrand suggested his model known as Bertrand model and considered prices as competitive means. In 1950, the Nash equilibrium point detected by the famous American mathematician John Nash has been raised in the equilibrium theory of noncooperative games. In general, this theory provided research with significant tools in studying oligopoly models. As a result, Nash point became important aspect in studying such models and was found in several Cournot and duopoly games. Indeed, one has to highlight game theory as an important tool of applied mathematics that was used and adopted to study the interaction among players in those economic games. It has been used in studying the strategic interactions among rational competitors and decision-makers. Game theory was extensively developed in the 1950s by many scholars and was explicitly applied to all fields of social science, as well as logic, biology, economics systems, and computer science. Recently, the theory of game has extensively been applied in economics for studying competitions among competing firms. Even though static equilibrium analysis of such games is limited and has not provided valuable information on games’ evolutions, it has been investigated by many researchers (see, for example, [2]). Studying the game evolution through dynamics gives more information on the stability/instability of Nash equilibrium, the types of bifurcations by which it loses its stability, and the topological properties of basin of attraction in the phase plane. Rand [3] was the first to introduce and analyze some of these dynamic characteristics.Firm in duopoly competition simultaneously chooses its decision or strategy taking into account two important factors: its own actions and actions by its competitor. In literature, several works have been reported for studying the dynamic characteristics of duopoly games. These dynamic characteristics included bifurcation and its types, either flip or Neimark–Sacker; peculiar shapes of attractive basins for attracting sets and chaotic attractors; and multistability phenomena. Furthermore, studying such games requires adopting specific kinds of utility consumer preferences. The most popular one that is well-known as Cobb–Douglas utility was adopted in many studies in literature. The reason is that Cobb–Douglas utility is a particular functional form of the production function. It has been widely used to represent the technological relationship between the inputs and outputs that can be produced by those inputs. Moreover, literature reported many other utility functions such as constant elasticity of substitution (CES) and Singh and Vives utility. One can get more information and properties of those utilities from literature ([4–10]).The dynamic characteristics of duopoly models are complex and reveal interesting information on their evolutions. Therefore, many interesting results on these complex dynamics have reported many interesting results on the stability of Nash equilibrium point of these games and opened the routes to attract many researchers. Some related works and results for these models are reported here. For example, a game of two competing firms whose decision variables are different (Cournot duopoly game) has been analyzed in [11]. Other studies on duopoly and triopoly in which firms seek the optimality of production have been introduced and investigated in ([12–14]). In [15], a Cournot duopoly game whose players are rational, adopting bounded rationality rule and seeking maximization of their relative profits, has been introduced and analyzed. Using a theoretical framework, a Cournot–Bertrand game has been formalized in [8]. In [16], Cournot–Bertrand model has been studied by Naimzada et al. using a two-dimensional discrete linear map. They focused on investigating the dynamic characteristics of Nash equilibrium points and its instability through two types of bifurcations. Additionally, analyzing the game’s evolution requires representing it by a discrete time map so that firms can update of their productions. This requires some kinds of adjustment rules. Literature has presented and reported many of these rules. The first rule among those rules which has been extensively adopted is known as bounded rationality. It has been deeply adopted in many studies for process of modeling discrete time maps representing the evolutions of these games. It is described as a gradient rule because it depends on firms’ marginal profits and requires firms to carry out an estimation of their marginal profits to determine whether they increase or decrease, so they can update their outputs next period of time. There are also other mechanisms which have been used in the modeling process such as tit-for-tat rule and local monopolistic approximation mechanism ([17–27]).Based on a linear inverse demand functions and adopting the bounded rationality mechanism, a Cournot duopoly game with players seeking the maximization of their average profits is introduced and analyzed in this manuscript. Most firms do not rank profits as the major goal. The work of modem firms is so complex that they do not think merely about profit maximization. The basis of the difference between the objectives of the neoclassical firm and the modem corporation arises from the fact that the profit maximization objective relates to the entrepreneurial behavior while modem corporations are motivated by different objectives because of the separate roles of shareholders and managers. Due to the above discussion, we have introduced a general case with different objectives. The game is described by a nonlinear discrete dynamic map that possesses four equilibrium points, three of which are boundary points, and the fourth is an interior point and coincides with Nash equilibrium point. We use in this manuscript local and global analysis to investigate the dynamics of these equilibrium points. The analysis shows that Nash point loses its local stability through flip bifurcation. The structure of the game’s map shows that it is being attracted to the origin, and hence its coordinate axes form an invariant manifold. Consequently, the map’s dynamics are studied by a one-dimensional map conjugating the famous logistic map. It is also shown that under some restrictions, the synchronization phenomena are studied and show that the diagonal in map’s phase plane forms an invariant manifold. Furthermore, our findings show that the game’s map is noninvertible and its phase plane is divided into three zones that make it belong toZ4−Z2−Z0 type.The sections in the current paper are summarized as follows. In Section2, the discrete form of the game’s map is introduced. The equilibrium point and its stability are investigated both analytically and numerically through local analysis in Section 3. In Section 4, the invariant manifold is studied. The basin of attraction and critical curves are analyzed in Section 5. In Section 6, the global analysis and synchronization phenomena are investigated based on some numerical experiments.


## 2. The Model

Suppose an economic market consisting of two competing firms (or players) whose decision variables are quantities,q1 and q2. Both firms adopt linear inverse demand functions given by(1)p1=a−q1−bq2,p2=a−q2−bq1.The horizontal differentiated products mainly refer to difference in extrinsic features between the products, like color and shape. Such differentiation mostly adapts to differences in consumer preferences. In (1), the parameter b refers to the degree of horizontal product differentiation that is taken in the interval 0,1. For b=0, each firm behaves as a monopoly player due to the independent varieties between firms’ products. b=1 means that both products are perfect substitutes or are homogenous, and hence p1=p2. In this paper, we focus on the case of b∈0,1. The parameter a denotes maximum price in market. Suppose that the competing firms adopt different nonnegative marginal costs and hence their total costs are taken as Ciqi=ciqi,i=1,2 (c1≠c2). Now, the total profit for each firm takes the following form.(2)π1=a−c1−q1−bq2q1,π2=a−c2−q2−bq1q2.Unlike many studies in literature, the current game assumes that both firms focus on maximization of their average profits rather than profits only. The average profits then become(3)φ1=π1+r1π2,φ2=π2+r2π1.The parameterri,i=1,2, where 0<ri<1, is directly used to measure altruism. That is to say, it explains the firms’ care of average performance. When ri=0,i=1,2, it means firms give care to their own profits only. Substituting (4) into (5), one gets(4)φ1=a−c1q1+r1a−c2q2−q12−r1q22−b1+r1q1q2,φ2=a−c2q2+r2a−c1q1−q22−r2q12−b1+r2q1q2,and then the marginal average profits ∂φi/∂qi,i=1,2 become(5)∂φ1∂q1=a−c1−2q1−b1+r1q2,∂φ2∂q2=a−c2−2q2−b1+r2q1.Maximizing the average profits requires firms to know some characteristics of market and the nature of competition. Knowing the characteristics of market needs complete information on the market and the strategies of firms’ rivals, and these become impossible. Therefore, firms believe that updating output productions requires watching the changes taking place in the marginal average profits. That is to say, if∂φi/∂qi>0, this means firms will increase their production next period of time; otherwise, they decrease productions or become naive. Accordingly, firms estimate production based on bounded rationality rule that has the following form.(6)qit+1=qit+viqi∂φi∂qi;i=1,2.Adopting bounded rationality rule makes us assume thatviqi=νiqi,i=1,2, where νi>0,i=1,2 is known as the speed of adjustment parameter. This gives qit+1−qit/qit∝∂φi/∂qi which makes relative production directly proportional to the firms marginal average profits. Equation (6) possesses some important characteristics that can be found in literature ([25–28]). Substituting (5) into (6), one simply has a two-dimensional discrete map that will be used throughout the manuscript to discuss the game’s repetition (or evolution).(7)Tq1,q2:q1t+1=q1t+ν1q1ta−c1−2q1t−b1+r1q2t,q2t+1=q2t+ν2q2ta−c2−2q2t−b1+r2q1t.The parametert=0,1,2,… denotes time steps. Setting the fixed points’ condition, Tq1,q2=q1,q2 in (7), one has four fixed points given by(8)eo=0,0,e1=a−c12,0,e2=0,a−c22,e∗=2a−c1−ba−c21+r14−b21+r11+r2,2a−c2−ba−c11+r24−b21+r11+r2.


## 3. Local Stability

In order to study the stability of the fixed points (8), we recall the Jacobian matrix at these points, and hence the following propositions are raised (their proofs are given in Appendix).Proposition 1.
The pointeo is an unstable repelling node.Proposition 2.
The pointe1 is a stable point provided that 0<ν1<2/a−c1 and a−c2<b/21+r2a−c1. It is an unstable saddle if 0<ν1<2/a−c1 and a−c2>b/21+r2a−c1.Proposition 3.
The pointe2 is a stable point provided that 0<ν2<2/a−c2 and a−c1<b/21+r1a−c2. It is an unstable saddle if 0<ν2<2/a−c2 and a−c1>b/21+r1a−c2.Proposition 4.
The pointe∗ that is called Nash point is locally stable provided that(9)0<2Av1+2Bv2−4ABv1v2+b21+r11+r2AB<4,A=2a−c1−ba−c21+r14−b21+r11+r2,B=2a−c2−ba−c11+r24−b21+r11+r2.Proposition 5.
The pointe∗ loses its stability due to flip bifurcation only.The obtained results show that the stability condition of Nash point depends on the parametersa,b,c1,c2,r1,r2 and the speed of adjustment parameters v1 and v2. Some numerical experiments are carried out here to validate the theoretical results given above. As in many studies in literature, we keep our focus on studying the influences of speed parameters on the stability of Nash while fixing the other parameters. Let us suppose the set of parameters a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5, and v1=1 (or v2=3.852) on varying v2 (or on varying v1). Figures 1(a) and 1(b) show the one-dimensional period-doubling bifurcation diagram on varying v1 with respect to the two decisional variables q1 and q2. As one can see, both figures present flip bifurcation diagrams. In Figures 1(c) and 1(d), one can see the flip bifurcation on varying the speed parameter, v2. In Figure 1(e), the phase plane of period-8 cycle with Nash equilibrium point is depicted. Figure 1(f) shows a one-piece chaotic attractor. In Figures 1(g) and 1(h), a two-piece disconnecting chaotic attractor is given. In these figures the marginal costs are fixed.Figure 1
(a)–(d) 1D bifurcation diagram on varying bothv1 and v2 at a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5, and v1=1 (or v2=3.852) on varying v2 (or on varying v1). (e) The phase plane of period-8 cycle at a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=3.33, and v2=3.852. (f) The phase plane of one-piece chaotic attractor at a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=4.5696, and v2=3.852. (g) Two-piece chaotic attractor at a=1,b=0.55,c1=0.4,c2=0.3,r1=0.3,r2=0.25,v1=4.6689, and v2=3.6. (h) Two-piece chaotic attractor at a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=4.03, and v2=3.852.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4. The Invariant Manifold

Looking at the map (7), one can see that the point 0,0 is trapping the map. That is to say, if q1t=0 or q2t=0, one gets q1t+1=0 or q2t+1=0, and hence both coordinate axes become invariant axes for the map. Therefore, the invariant set made of the axes q1 and q2 forms an invariant manifold for T. Accordingly, the dynamics of the map (7) on the invariant axes can be described by the one-dimensional map given by(10)qit+1=1+a−ciνi1−2νi1+a−ciνiqitqit,i=1,2.It is easy to see that the map (10) is topologically equivalent to(11)uit+1=μiuit1−uit,i=1,2,through the transformation given below.(12)qi=1+a−ciνi2ν1ui,i=1,2,and μi=1+a−ciνi, i=1,2.


### 4.1. Dynamic Analysis

For equation (11), suppose the following function.(13)σu=μu1−u.Function (13) shows the logistic equation, and μ is the parameter responsible for its dynamics. It has a first derivative given by σ´u=μ1−2u. Its maximum value occurs at μ/4 at the point u=1/2. In addition, σ0=0 and σ1=0, which means that σu∈0,1 on the open interval 0<μ<4 and for all u∈0,1. Moreover, it is easy to see that (11) has only two fixed points, that is, u¯=0 and u¯=1−1/μ, where μ>1. Consequently, σ´0=μ, and then u¯=0 is stable if μ∈0,1 and unstable otherwise. The other fixed point u¯=1−1/μ is stable if μ∈1,3 and unstable if μ>3. In addition, one can see that the map (11) has two fixed points, u¯=0 and u¯=1−1/μ, where μ>1. Hence, σ´0=μ, and then u¯=0 is stable if μ∈0,1 and unstable otherwise. The other fixed point u¯=1−1/μ is stable if μ∈1,3 and unstable if μ>3.The simulation in Figure2(a) presents a 1D bifurcation diagram on varying μ that confirms that the point u¯=1−1/μ is stable in μ∈1,3. In Figure 2(b) a cobweb diagram is given at μ=2.99 in the stability interval and shows a stable fixed point, while in Figure 2(c), an unstable behavior is plotted at μ=3.75, and hence the instability of the fixed point is confirmed.Figure 2
(a) 1D bifurcation diagram of mapσu with respect to μ1. Iterations of σμ at (b) μ = 2.99 and (c) μ = 3.75. (d) 1D bifurcation diagram for the map (14) at a = 1, b = 0.35, c = 0.4, and r = 0.7.
(a)(b)(c)(d)Proposition 6.
Once the critical valueνi=2/a−ci,i=1,2, is reached, the trajectories of T starting on the invariant axes diverge when νi∈2/a−ci,+∞,i=1,2.Proof.
The proof is straightforward asμi=1+a−cνi, and our discussion given above confirms the topological equivalence between the map (7) and (11) on which the trajectories of (11) diverge when μi>3,i=1,2.


### 4.2. The Symmetric Case

In this subsection, the symmetric casev1=v2=v,c1=c2=c, and r1=r2=r is studied. As one can see, the map T defined in (7) is symmetric. That is to say, its structure does not change if the decision variables q1 and q2 are swapped; i.e., T∘S=S∘T, where S:q1,q2⟶q2,q1. This means that the diagonal defined by Δ=q1,q2:q1=q2 forms an invariant manifold. It implies that the dynamics of any trajectories starting on the diagonal (q10=q20) will get back to the diagonal for every t. In this case, the dynamic behavior of map (7) under the restriction imposed on Δ will be described by the one-dimensional map given below.(14)TΔ:q′=gq≔q+νqa−c−2+b1+r1q.That is to say, any synchronized trajectories (i.e.,q1t=q2t for every time step t) will be governed by TΔ:Δ⟶Δ, and the following proposition is raised.Proposition 7.
The map (14) is unimodal and concave, limq⟶0+gq=0, and limq⟶+∞gq=−∞. This means that there is a threshold q¯=1+va−c/v2+b1+r at gq¯=0. When the model becomes unimodal, one gets the following critical point.(15)qcr=1+va−c2v2+b1+r.Proof.
The proof is straightforward by looking at the functiongq that has a global maximum without any local maximum, and hence it is unimodal. The sign of the second derivative of gq proves that the map is concave and finally g′qcr=0.
Now, the mapTΔ has only one fixed point, that is, q¯, and this point is locally stable if ν<2/a−c while flip bifurcation takes place at v=vflip:=2/a−c. Figure 2(d) shows a 1D bifurcation diagram where flip bifurcation takes place at v≈3.33.


## 4.1. Dynamic Analysis

For equation (11), suppose the following function.(13)σu=μu1−u.Function (13) shows the logistic equation, and μ is the parameter responsible for its dynamics. It has a first derivative given by σ´u=μ1−2u. Its maximum value occurs at μ/4 at the point u=1/2. In addition, σ0=0 and σ1=0, which means that σu∈0,1 on the open interval 0<μ<4 and for all u∈0,1. Moreover, it is easy to see that (11) has only two fixed points, that is, u¯=0 and u¯=1−1/μ, where μ>1. Consequently, σ´0=μ, and then u¯=0 is stable if μ∈0,1 and unstable otherwise. The other fixed point u¯=1−1/μ is stable if μ∈1,3 and unstable if μ>3. In addition, one can see that the map (11) has two fixed points, u¯=0 and u¯=1−1/μ, where μ>1. Hence, σ´0=μ, and then u¯=0 is stable if μ∈0,1 and unstable otherwise. The other fixed point u¯=1−1/μ is stable if μ∈1,3 and unstable if μ>3.The simulation in Figure2(a) presents a 1D bifurcation diagram on varying μ that confirms that the point u¯=1−1/μ is stable in μ∈1,3. In Figure 2(b) a cobweb diagram is given at μ=2.99 in the stability interval and shows a stable fixed point, while in Figure 2(c), an unstable behavior is plotted at μ=3.75, and hence the instability of the fixed point is confirmed.Figure 2
(a) 1D bifurcation diagram of mapσu with respect to μ1. Iterations of σμ at (b) μ = 2.99 and (c) μ = 3.75. (d) 1D bifurcation diagram for the map (14) at a = 1, b = 0.35, c = 0.4, and r = 0.7.
(a)(b)(c)(d)Proposition 6.
Once the critical valueνi=2/a−ci,i=1,2, is reached, the trajectories of T starting on the invariant axes diverge when νi∈2/a−ci,+∞,i=1,2.Proof.
The proof is straightforward asμi=1+a−cνi, and our discussion given above confirms the topological equivalence between the map (7) and (11) on which the trajectories of (11) diverge when μi>3,i=1,2.


## 4.2. The Symmetric Case

In this subsection, the symmetric casev1=v2=v,c1=c2=c, and r1=r2=r is studied. As one can see, the map T defined in (7) is symmetric. That is to say, its structure does not change if the decision variables q1 and q2 are swapped; i.e., T∘S=S∘T, where S:q1,q2⟶q2,q1. This means that the diagonal defined by Δ=q1,q2:q1=q2 forms an invariant manifold. It implies that the dynamics of any trajectories starting on the diagonal (q10=q20) will get back to the diagonal for every t. In this case, the dynamic behavior of map (7) under the restriction imposed on Δ will be described by the one-dimensional map given below.(14)TΔ:q′=gq≔q+νqa−c−2+b1+r1q.That is to say, any synchronized trajectories (i.e.,q1t=q2t for every time step t) will be governed by TΔ:Δ⟶Δ, and the following proposition is raised.Proposition 7.
The map (14) is unimodal and concave, limq⟶0+gq=0, and limq⟶+∞gq=−∞. This means that there is a threshold q¯=1+va−c/v2+b1+r at gq¯=0. When the model becomes unimodal, one gets the following critical point.(15)qcr=1+va−c2v2+b1+r.Proof.
The proof is straightforward by looking at the functiongq that has a global maximum without any local maximum, and hence it is unimodal. The sign of the second derivative of gq proves that the map is concave and finally g′qcr=0.
Now, the mapTΔ has only one fixed point, that is, q¯, and this point is locally stable if ν<2/a−c while flip bifurcation takes place at v=vflip:=2/a−c. Figure 2(d) shows a 1D bifurcation diagram where flip bifurcation takes place at v≈3.33.


## 5. Basin of Attraction

Here, some properties of the map (7) regarding its attractive basin are given. According to [29], it is possible to see that the coordinate axes and their preimages of any rank construct boundaries for initial conditions of any nondiverging trajectories. It can be seen that the dynamics of map (7) can be restricted to the invariant axes and can be rewritten in the one-dimensional map as follows.(16)q´i=giqi≔qi1+νia−ci−2νiqi,i=1,2,which possesses a positive fixed point, qi∗=1+νia−ci/2νi,i=1,2. It is easy to see that giqi is concave and increasing function provided that qi<1+νia−ci/4vi. In addition, for any points in the form q1,0,q1>0, or 0,q2,q2>0, one gets s or q´1=0,q´2>0. Conversely, if qi>1+νia−ci/4vi, then giqi is concave and unimodal. In the last case, one can get the nonnegative points q˜ibound=1+νia−ci/2vi and q˜icr=1+νia−ci/4vi at q´i=0 and g´iq˜icr=0, respectively. If g´iq˜icr<q˜ibound, i.e., if νia−ci<3, this means that trajectories along the invariant axes are bounded if the initial states lie on the line segment ωi=O,Oi−1,i=1,2, where Oi−1=q˜ibound≠0 is the rank-1 preimage of the point of origin. In contrast, any other trajectories along the invariant axes starting from an initial state outside ωi become unfeasible. Furthermore, calculating the eigenvalues at a generic point q1,q2 on the invariant axis qi gives λ˜i=1+νia−ci−4viqi,i=1,2. This means that trajectories commencing close to these invariant axes with positive initial states are repelled by those axes. In Figure 3(b), the line segments ωi,i=1,2, and their inverses ωi−1,i=1,2, are depicted at the set of parameter values a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=4.5696, and v2=3.852 where the point O3−1 is the fourth preimage point of the origin and is given by(17)O3−1=q∘1,q∘2,q∘1=2v2−b1−α1v1+ν1v22a−c1−ba−c21−α14−b21−r11−r2,q∘1=2v1−b1−α2v2+ν1v22a−c2−ba−c11−α24−b21−r11−r2.Figure 3
(a) The shape of mapTΔ at the values a=1,b=0.35,c=0.4,r=0.7, and ν=4.9. At a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=4.5696, and v2=3.852v2=3.852. (b) The segments ωi,i=1,2, and their inverses ωi−1. (c)-(d) The locus LC−1 and the critical curve LC with the corresponding zones.
(a)(b)(c)(d)It is also noted that this point is the intersection of both inversesω1−1 and ω2−1.


### 5.1. Critical Curves

Besides the fact that the map (7) is nonlinear, it is also noninvertible map, which means that the rank-1 preimages of a point q´1,q´2 may not exist or may be multivalued. As one can see, the origin point possesses four rank-1 preimages, and hence the phase plane of map (7) includes Z4 zone, the zone whose points have four rank-1 preimages. It may also include Z2, the zone whose points possess two rank-1 preimages, and Z0, the zone with no rank-1 preimages. To investigate these zones, let us calculate the rank-1 preimages for the points along the invariant axes and take the forms p,0 and 0,q. Substituting these points into the map (7), one gets the following algebraic systems.(18a)p=q1+ν1q1a−c1−2q1−b1+r1q2,0=q2+ν2q2a−c2−2q2−b1+r2q1,(18b)0=q1+ν1q1a−c1−2q1−b1+r1q2,q=q2+ν2q2a−c2−2q2−b1+r2q1.Solving the second equation of (18a) gives the following: q2=0 or 0=1+ν2a−c2−2q2−b1+r2q1, and then one gets two case.Case 1.
Substitutingq2=0 in the first equation of (18a) gives the following quadratic equation.(19a)2v1q12−1+v1a−c1q1+p=0.
This means that the point in the formp,0 may have no preimages or two distinct real-1 preimages.Case 2.
At0=1+ν2a−c2−2q2−b1+r2q1, one gets the algebraic system(19b)p=q1+ν1q1a−c1−2q1−b1+r1q2,0=1+ν2a−c2−2q2−b1+r2q1.
Solving (19b), we get the following quadratic equation.(19c)v14−b21−r11−r2q12−2a−c1v1−b21−r11+a−c2v1q1+p=0.
This means that in this case, the pointp,0 may have no preimages or two distinct real rank-1 preimages. Combining the two cases, one can conclude that such points in the invariant axes may have four, two, or no preimages. For the interior points in the form p,q≠0,0, it is complicated to categorize such points as belonging to one of the zones Z0,Z2, or Z4 due to the difficulty of solving the algebraic system (18a) at this point. For this reason, numerical simulation is used to calculate critical curves. The critical curves are responsible for dividing the map’s phase plane into these zones. The critical curve of rank-1 is denoted by LC and forms a locus of coincident points of rank-1 preimages located on a set denoted by LC−1. Since the map (7) is continuously differentiable, LC−1 can be defined as the locus of points where the associated Jacobian determinant is vanished; i.e., LC−1⊆q1,q2∈ℝ2:detJq1,q2=0, where Jq1,q2 is defined by (A.1) in the appendix. Therefore, LC represents the rank-1 image of LC−1 under the map T given in (7); i.e., LC=TLC−1. The expression for detJq1,q2=0 is cumbersome and is omitted here, and we instead do numerical simulation for presenting both LC−1 and LC. Figure 3(c) shows that LC−1 consists of two branches LC−1a and LC−1b while Figure 3(d) presents the critical curve that is also formed by two branches, LCa and LCb. It is also clear that the map (7) is not invertible as its phase plane is divided into three zones, Z0,Z2, and Z4.


## 5.1. Critical Curves

Besides the fact that the map (7) is nonlinear, it is also noninvertible map, which means that the rank-1 preimages of a point q´1,q´2 may not exist or may be multivalued. As one can see, the origin point possesses four rank-1 preimages, and hence the phase plane of map (7) includes Z4 zone, the zone whose points have four rank-1 preimages. It may also include Z2, the zone whose points possess two rank-1 preimages, and Z0, the zone with no rank-1 preimages. To investigate these zones, let us calculate the rank-1 preimages for the points along the invariant axes and take the forms p,0 and 0,q. Substituting these points into the map (7), one gets the following algebraic systems.(18a)p=q1+ν1q1a−c1−2q1−b1+r1q2,0=q2+ν2q2a−c2−2q2−b1+r2q1,(18b)0=q1+ν1q1a−c1−2q1−b1+r1q2,q=q2+ν2q2a−c2−2q2−b1+r2q1.Solving the second equation of (18a) gives the following: q2=0 or 0=1+ν2a−c2−2q2−b1+r2q1, and then one gets two case.Case 1.
Substitutingq2=0 in the first equation of (18a) gives the following quadratic equation.(19a)2v1q12−1+v1a−c1q1+p=0.
This means that the point in the formp,0 may have no preimages or two distinct real-1 preimages.Case 2.
At0=1+ν2a−c2−2q2−b1+r2q1, one gets the algebraic system(19b)p=q1+ν1q1a−c1−2q1−b1+r1q2,0=1+ν2a−c2−2q2−b1+r2q1.
Solving (19b), we get the following quadratic equation.(19c)v14−b21−r11−r2q12−2a−c1v1−b21−r11+a−c2v1q1+p=0.
This means that in this case, the pointp,0 may have no preimages or two distinct real rank-1 preimages. Combining the two cases, one can conclude that such points in the invariant axes may have four, two, or no preimages. For the interior points in the form p,q≠0,0, it is complicated to categorize such points as belonging to one of the zones Z0,Z2, or Z4 due to the difficulty of solving the algebraic system (18a) at this point. For this reason, numerical simulation is used to calculate critical curves. The critical curves are responsible for dividing the map’s phase plane into these zones. The critical curve of rank-1 is denoted by LC and forms a locus of coincident points of rank-1 preimages located on a set denoted by LC−1. Since the map (7) is continuously differentiable, LC−1 can be defined as the locus of points where the associated Jacobian determinant is vanished; i.e., LC−1⊆q1,q2∈ℝ2:detJq1,q2=0, where Jq1,q2 is defined by (A.1) in the appendix. Therefore, LC represents the rank-1 image of LC−1 under the map T given in (7); i.e., LC=TLC−1. The expression for detJq1,q2=0 is cumbersome and is omitted here, and we instead do numerical simulation for presenting both LC−1 and LC. Figure 3(c) shows that LC−1 consists of two branches LC−1a and LC−1b while Figure 3(d) presents the critical curve that is also formed by two branches, LCa and LCb. It is also clear that the map (7) is not invertible as its phase plane is divided into three zones, Z0,Z2, and Z4.


## 6. Global Analysis

In economic models, global analysis is required to investigate the topological structure of basin of attraction that is not recognized by local analysis and whose qualitative changes occur due to change in parameter values. The long-term behavior of variables of these economic models given initial conditions is indeed necessary in order to show the influences that may result if these initials are taken far away from a fixed point of an attracting set. The numerical analysis shows two important aspects for the map (7), the asymmetric and symmetric cases. Let us first begin analyzing the asymmetric case by adopting the following set of parameter values: a=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5,v1=3.3285, and v2=3.852. This set is given in Figure 4(a) which represents a period-8 cycle with Nash equilibrium point and its attractive basins. These basins are colored orange and green while the escaping (or infeasible points) domain is plotted in gray color. It is also clear that the escaping domain forms a connecting set. By increasing the speed parameter of the first firm to v1=4.03 while keeping the other values fixed, the dynamic of map (7) turns to a chaotic attractor that consists of two disconnecting areas as shown in Figure 4(b). It can be observed that the branch LCb of the critical curve LC is moving upward. With further increase in this parameter to v1=4.4675, a one-piece chaotic attractor emerges in Figure 4(c). Its attracting domain is displayed in white color and is separated from the infeasible domain by the segments ω1−1 and ω2−1. At v1=4.647, with the other parameter values being fixed, the chaotic attractor continues to appear, the critical curve LCb becomes tangent to the boundary ω2−1, and hence contact bifurcation is raised as given in Figure 4(d). This contact takes place only when one branch LCb of the critical curve which separates Z2 zone from Z0 zone becomes tangent to ω2−1. This contact with the boundary of infeasible domain occurs from one side due to asymmetry between the two branches of map (7). This contact bifurcation becomes more visible in Figure 4(e) at the parameters set ν1=4.63 and ν2=3.95 while the other parameters are fixed. As one can see, the contact bifurcation causes disconnection in the infeasible domain because an area h0 from Z0 zone moves into Z2 zone; i.e., LCb∩ω2−1=h0, which is out of the quadrilateral OO1−1O3−1O2−1. This kind of contact bifurcation is well-known (see [29]). Just after it, the area h0 from the infeasible domain enters inside the region Z2. As a result, the points belonging to this region have now two distinct preimages. Since this region belongs to Z0, their preimages will belong to the same region and form holes lying in the attraction domain that is bounded by the quadrilateral shape OO1−1O3−1O2−1. Indeed, only at νi=2/a−ci,i=1,bothh0 and its preimage h0−1 (that is formed by two areas connected by the branch LC−1b) have a contact with coordinate axes. It is also clear that both h0 and h0−1 are entirely located in the region Z2. Therefore, each of them has further preimages that form smaller holes of the infeasible domain located in the attracting domain and are bounded by rank-2 preimages of the segments ω1 and ω2.Figure 4
Ata=1,b=0.35,c1=0.4,c2=0.3,r1=0.7,r2=0.5, the basin of attraction for (a) period-8 cycle, ν1=3.3285,v2=3.852. (b) Two-piece chaotic attractor, ν1=4.03,v2=3.852. (c) One-piece chaotic attractor, ν1=4.4675,v2=3.852. (d) One-piece chaotic attractor, ν1=4.647,v2=3.852. (e) One-piece chaotic attractor, ν1=4.63,v2=3.3.95. At a=1,b=0.35,c=0.4,r=0.25, the basin of attraction for (f) period-4 cycle, ν=4.1772185. (g) Four closed rings, ν=4.2199558. (h) Period-20 cycle, ν=4.2818322.
(a)(b)(c)(d)(e)(f)(g)(h)In contrast, the symmetric case is now investigated. This case happens atv1=v2=v,c1=c2=c, and r1=r2=r. This means the map (7) can be rewritten in the following form.(20)T∗q1,q2:q1t+1=q1t+νq1ta−c−2q1t−b1+rq2t,q2t+1=q2t+νq2ta−c−2q2t−b1+rq1t.The map (20) is symmetric on the diagonal Δ, and studying its dynamics through local analysis around Δ is related to synchronization (intermittency). In this regard, one can say that a dynamic behavior of the map is synonymized if there exists a step time t∘ such that q1t=q2t for any t>t∘. Therefore, in this case, the dynamics can be studied by recalling the Jacobian of map (20) given by(21)Jq1,q2=ℓqmqmqℓq,where, ℓq=1+va−c−4+1+rbq and mq=−b1+rvq. Then, the eigenvalues can be given as follows.(22)λ‖=ℓq+mq=1+va−c−22+1+rbq,λ⊥=ℓq−mq=1+va−c−4q,and their corresponding eigenvectors are 1,1 and 1,−1, respectively. In this case, e∗=a−c/2+b1+r,a−c/2+b1+r, and then the following proposition is raised.Proposition 8.
The pointe∗ is locally stable provided that v<2/a−c. At v=vf=2/a−c, the point e∗ undergoes flip bifurcation.Proof.
Substitutinge∗ into (22), we get(23)λ‖=1−a−c2−1+rb2+b1+rv,λ⊥=1+va−c.
Simple calculations makeλ⊥<1 if v<2/a−c, and λ‖<1 if v<22+1+rb/a−c2−1+rb. Combining the two conditions on ν gives v<2/a−c, and this completes the proof.
From the above proposition, one can deduce that the parameterv plays a key role in the stability of e∗. This means that high reactivity of each firm with respect to a marginal change in its average profit should exist. Furthermore, the importance of the above proposition comes from its generalization when an attracting cycle for the map T∗ on Δ exists. To see that, suppose that a k-cycle given by q1,q1,…,qm,qm of the map T∗ is embedded in the invariant diagonal Δ which corresponds to the cycle q1,…,qm. Then, the multipliers for this k-cycle are given by(24)λ‖k=∏i=1kℓqi+mqi,λ⊥k=∏i=1kℓqi−mqi,whose eigenvectors are 1,1 and 1,−1. As shown, the stability of e∗ is governed by λ⊥<1 and is confirmed by λ‖<1. In contrast, a complex chaotic attractor H for the map T∗ says that H is asymptotically stable if and only if all trajectories belonging to H are transversely attracting. This asymptotic stability condition of H is given in terms of transverse Lyapunov exponent as follows.(25)Λ⊥=limn⟶∞∑i=1nlnλ⊥qi,where q0∈H and qi is the trajectories generated by the map(26)Ts:qt+1=qt+νqta−c−2+b1+rqt.
The following definitions are given [29].Definition 1.
A chaotic attractorH is asymptotically stable if it is Lyapunov stable; i.e., for every neighborhood U of H, there exists a neighborhood V of H such that TnV⊂U∀n≥0 and the attractive basin BH possesses a neighborhood of H.
A spectrum of Lyapunov exponents can be defined according to the initial conditions as follows [30].(27)Λ⊥min<⋯<Λ⊥nat<⋯<Λ⊥max,where Λ⊥nat denotes a Lyapunov exponent evaluated at a generic trajectory in the chaotic attractor H. If Λ⊥max is nonpositive, then a set is called Lyapunov attractor. While Λ⊥max>0 and Λ⊥nat<0, the set becomes no longer Lyapunov stable and becomes Milnor attractor. It is defined as follows.Definition 2.
A closed invariant setH is called a weak attractor in Milnor sense if its stable set (the basin of attraction) BH has positive Lebesgue measure.
Now, we give some experiments through numerical analysis to confirm the above findings. Let us assume the following set:a=1,b=0.35,c=0.4,v=4.1772185, and r=0.25. At this set, the dynamic of map (20) is represented by period-4 cycle given in Figure 4(f). Its attractive basin and the basin of Nash equilibrium point are colored green and orange, respectively, and are bounded by ω1 and ω2 and their inverses. The infeasible domain of this cycle is given in gray color. By keeping the same set of values and increasing v to 4.2199558, four closed invariant rings emerge around the diagonal as depicted in Figure 4(g). All these rings are symmetric around the invariant diagonal Δ and lie within the region Z2 only when two of them are tangent to the branch of critical curve LCb. By further increasing ν to 4.2818322, the four rings convert into a period-20 cycle as shown in Figure 4(h) whose basin of attraction is colored cyan while the basin of Nash point is denoted by orange color. At ν=4.414216, while the other parameter values are fixed, a Milnor attractor is raised on the invariant diagonal as presented in Figure 5(a). This means the dynamics of map (20) are captured by a two-piece chaotic attractor on the invariant diagonal Δ that is in Milnor sense. This indicates that any cycle embedded in the diagonal Δ is transversely unstable, and hence all trajectories approaching that diagonal are exploded away by transversely unstable cycles. By further increasing ν to 4.763885, a one-piece chaotic attractor emerges as shown in Figure 5(b). This chaotic attractor is symmetric on the diagonal because synchronization still occurs. The symmetry of this chaotic attractor is given in Figure 5(c) where its boundary at the same set of parameter values is given by the critical curves LCa,LCb,LC1a,LC1b,LC2a, and LC2b. In Figure 5(d), the transverse Lyapunov exponent is given. It gives better understanding on the map’s dynamics in the long term. As one can see, the bifurcation diagram corresponding to the transverse Lyapunov exponent confirms the existence of chaotic attractor on the invariant diagonal in the interval ν∈4.4,4.7 except for windows of some periodic cycles. It is also noted that in this interval of the speed parameter ν, the Lyapunov exponent changes its sign many times until it eventually becomes positive for large values of the speed parameter ν. From economic perspective, this implies that the long-term dynamics of the map (20) can be very sharp due to small changes in this parameter. In other words, firms commencing with the same initial states may be different in the long term.Figure 5
Ata=1,b=0.35,c=0.4,r=0.25, the basin of attraction for (a) Milnor attractor, ν=4.414216. (b) One-piece chaotic attractor, ν=4.763885. (c) The boundaries of chaotic attractor in (b). (d) Bifurcation diagram on varying ν and the corresponding transverse Lyapunov exponent for ν∈4.4,4.7.
(a)(b)(c)(d)


## 7. Conclusion

In this paper, some dynamic characteristics such as multistability, synchronization phenomena, and complex structure of some attractive basins are investigated for a nonlinear two-dimensional noninvertible map describing a Cournot duopoly game with an invariant one-dimensional manifold where synchronization takes place. It has been proved that, with identical firms, the game’s map becomes symmetric and both coordinate axes and diagonalq1=q2 form invariant manifolds. The main results in this paper concern the global analysis in the asymmetric and symmetric case where synchronization has taken place. Such global analysis has focused on the properties of critical curves of the noninvertible map. The phase plane of game’s map has been divided into three zones of preimages, that is, Z0,Z2, and Z4. We have shown that for a set of parameter values, the one-dimensional Milnor attractor to which the synchronized trajectories arise has been obtained. After that, for increasing the value of the speed parameter, the Milnor attractor is transformed into a chaotic attractor. Indeed, weak chaotic attractor in the Milnor sense has been confirmed by transverse Lyapunov exponent.



---

*Source: 8677470-2022-03-09.xml*