# On the Dynamics of Cournot Duopoly Game with Governmental Taxes

Authors: S. S. Askar
Journal: Complexity
 (2022)

---

## Abstract

A quadratic utility function is introduced in this paper to study the dynamic characteristics of Cournot duopoly game. Based on the bounded rationality mechanism, a discrete dynamical map that describes the game’s dynamic is obtained. The map possesses only one equilibrium point which is Nash point. The stability conditions for this point are analyzed. These conditions show that the point becomes unstable due to two bifurcation types that are flip and Neimark–Sacker. The synchronization property for that map is studied. Through local and global analysis, some dynamics of attracting sets are investigated. This analysis gives some insights on the basis of those sets and the shape of the critical curves. It also shows some lobes found in those attracting sets which are constructed due to the origin focal point.


---

## Body

## 1. Introduction

The Cournot duopoly game as described in literature consists of two players (or two competing firms), each wants to seek the optimal quantity of his/her production. Many scholars in literature have studied the dynamic of such games which have been described by discrete maps ([1–5]). The Nash equilibrium point of such games and its stability conditions were the core of study by scholars. Different types of bifurcations have been reported in literature as causes for the equilibrium point to be unstable. In addition, many utility functions have been adopted to model such games. Of which are Cobb–Douglas, constant elasticity of substitution (CES) and others. Interested readers are advised to see some applications of those utility functions and their properties in literature ([1–10]).The present paper introduces a simple quadratic utility function. This function gives under Lagrangian function and its first-order condition price functions same as given by Cobb–Douglas. Our introduced utility function may be considered as a special case of Singh and Vives function [1]; however, the later does not give linear prices as reported in literature [6] if one apply Lagrangian function and its first-order conditions on it. Anyway, we highlight here some important works that have studied such kind of games. For instance, in [11], a Cournot–Bertrand duopoly game whose products are considered to be differentiated has been studied. Few useful investigations on such games and their optimality have been reported in the following studies ([12–14]). In [15], Tremblay et al. have analyzed the differentiated products in a Cournot–Bertrand duopoly game. They have focused on studying the static game and the stability/instability conditions of the game’s Nash equilibrium point. In [8], a theoretical framework of a Cournot–Bertrand game has been introduced. Based on a two-dimensional discrete linear map, Naimzada et al. [16] have analyzed a model of Cournot–Bertrand and studied its dynamic characteristics. Naimzada et al. have adopted different mechanisms such as the best response and adaptive adjustment to introduce the model which has been used to describe the dynamic of the game.When studying such games, different types of adjustment mechanisms have been used in the modelling process. The modelling process means introducing the discrete map that describes the dynamic of game at discrete time steps. The most popular mechanism that has been used in the modelling process of such games is the bounded rationality approach. It has been considered as a gradient-based approach. Other mechanisms that has been used in few studies in literature are the naive mechanism, the tit-for-tat approach, and the approximation of local monopolistic or LMA mechanism. Information about those mechanisms and their properties can be founded in literature ([17–21]). Other interesting works on the extensions of those mechanisms and their applications have been reported ([22–26]).The present paper belongs to the category of Cournot duopoly game on which both competing firms adopt the bounded rationality approach and want to seek the optimal quantities of production in order to achieve the profit maximization. The current game differs from those in literature [27], on which we introduce a new utility function that has not been adopted before. In this paper, local and global analyses are performed to investigate the game’s dynamics. This includes the investigation of multiple stable attractors and analyzes the attractive basins for some attracting sets. The main results in such study focus on analyzing the dynamics of the unique Nash equilibrium of the map’s game. This includes investigating the types of bifurcations by which the equilibrium point may be unstable. The obtained results show that there are two types of bifurcations by which the Nash point becomes unstable. These types are flip and Neimark–Sacker bifurcations. In addition, the form of the game’s map possesses a focal point that is the origin which gives rise to some lobes affecting on the shape of the attractive basins of some attracting sets. Furthermore, the synchronization property is studied. Moreover, the critical curves are calculated and show that the phase plane of the game’s map belongs to Z1−Z3 type.The structure of paper is given as follows. The quadratic utility function and the discrete dynamic map describing the game are introduced in Section2. This section includes also the stability investigation of the Nash equilibrium point and the route to chaos due to two different types of bifurcations. In addition, local and global analyses including the synchronization property are discussed in this section. Furthermore, the critical curves that divide the phase plane of the game’s map into Z1 and Z3 regions are calculated. Finally, Section 3 concludes our obtained results and suggests some future works.


## 2. The Model with Tax

Let us first assume the following quadratic utility function:(1)U=Q−Q2=q1+q2−q1+q22,where Q=q1+q2 and qi,i=1,2 denotes the quantity produced by the player (firm) i. Our suggested model in this manuscript consists of two competing firms whose decision variables are the quantities produced by them. It is simple to see that the utility given in (1) is convex which economically means that the good qi,i=1,2 has an increasing marginal utility. Furthermore, the marginal utility of the good q1 does not depend on the good q2 (as ∂2U/∂q1q2=−2≠0). In addition, this utility is homogeneous. Assuming the budget constraint p1q1+p2q2=m,m>0, then we get the following maximization problem:(2)MaxUq1,q2,s.tp1q1+p2q2=m,where pi>0,i=1,2 represents the price of good qi,i=1,2. Solving (2), we get (see Appendix A)(3)pi=1Q=1q1+q2.Now, we assume that both firms detect the optimum according to maximizing their profits as follows:(4)maxq1,q2π1q1,q2=p1q1−Cq1−Taxq1,π2q1,q2=p2q2−Cq2−Taxq2,where Ciqi represents the cost of the quantity qi,i=1,2 and is taken as a linear cost Ciqi=ciqi where ∂Ciqi/∂qi=ci denotes a constant marginal cost. We assume also that the government has imposed a tax on each quantity as Taxqi=riqi,i=1,2 and ri∈0,1. Now, (4) can be represented by(5)π1q1,q2=q1Q−c1+r1q1,π2q1,q2=q2Q−c2+r2q2.Recalling the mechanism of bounded rationality [8], both firms can update their outputs as follows:(6)q1t+1=q1t+k1q1t∂π1q1t,q2t∂q1t,q2t+1=q2t+k2q2t∂π2q1t,q2t∂q2,where t=0,1,2,…. Using (5) in (6), we get the two-dimensional discrete dynamic map that describes the current game.(7)Tq1,q2:q1t+1=q1t+k1q1tq2Q2−c1−r1,q2t+1=q2t+k2q2tq1Q2−c2−r2.Now, we discuss the dynamic characteristics of the equilibrium point of map (7) in the next section.


### 2.1. The Equilibrium Point and Its Stability

The equilibrium point of the above map is a unique Nash equilibrium point and is obtained by settingq1t+1=q1t=q¯1 and q2t+1=q2t=q¯2. It has the following form:(8)O=q¯1,q¯2=c2+r2c1+c2+r1+r22,c1+r1c1+c2+r1+r22,which is a positive point. We should highlight here that the equilibrium point (8) is the same as the one given in [28] but when the model has free taxes r1=r2=0. The following propositions are used in the sequel to classify the type of the point O.Proposition 1.
Let the Jacobian matrix the map (7) at any equilibrium point possesses two eigenvalues λ1 and λ2. Then, this equilibrium point can be classified according to the following:(i)
Ifλ1,2<1, then it becomes an attracting node and it is stable(ii)
Ifλ1,2>1, then it is an unstable repelling node(iii)
Ifλ1<1 and λ2>1, (or λ1>1 and λ2<1) then it is an unstable saddle point(iv)
It is a nonhyperbolic point ifλ1=1 and λ2≠1 (or λ1≠1 and λ2=1)
In case, the eigenvalues have a complicated analytical form, we use the following proposition.Proposition 2.
Letτ and δ be the trace and determinant of the Jacobian matrix of (7), and suppose the following(9)g1=1−τ+δ,g−1=1+τ+δ,Ω=1−δ.
Then, we have the following properties:(i)
The equilibrium point is locally asymptotically stable ifg1>0,g−1>0, and Ω>0.(ii)
The equilibrium point becomes unstable through a flip bifurcation ifg1>0,g−1=0, and Ω>0. In this case, the two eigenvalues are real and pass through −1.(iii)
The equilibrium point becomes unstable through a transcritical or fold bifurcation ifg1>0,g−1>0, and Ω>0. In this case, the two eigenvalues are real and pass through 1.(iv)
The equilibrium point becomes unstable through Neimark–Sacker ifg1>0,g−1>0, and Ω<0. In this case, the two eigenvalues are complex and their modulus passes through 1.
The Jacobian matrix for the map (7) at the equilibrium point O becomes(10)JO=1−2c1+r1c2+r2k1c1+r1+c2+r2−c2+r2c1+r1−c2−r2k1c1+r1+c2+r2c1+r1c1+r1−c2−r2k2c1+r1+c2+r21−2c1+r1c2+r2k2c1+r1+c2+r2.
And its eigenvalues are(11)λ1,2=1−c1+r1c2+r2k1+k2c1+r1+c2+r2±±2c1+r1c2+r2c1+r1+c2+r2k12+k22−c12+r12+c22+r22+2c1r1+2c2r2k1k2c1+r1c2+r2.Proposition 3.
The equilibrium pointO is locally asymptotically stable if the following conditions are satisfied:(12)c1+r1c2+r2k1k2−4c1+r1c2+r2c1+r1+c2+r2k1+k2+4>0,2c1+r1c2+r2c1+r1+c2+r2k1+k2−c1+r1c2+r2k1k2>0.Proof.
Simple calculations show thatτ and δ for the Jacobian given in (10) become(13)τ=21−c1+r1c2+r2k1+k2c1+r1+c2+r2,δ=1−2c1+r1c2+r2k1+k2c1+r1+c2+r2+c1+r1c2+r2k1k2.
Substituting (13) in (9), we get g1=c1+r1c2+r2k1k2 which is always positive since ci,ri and ki,i=1,2 are positive parameters. Simplifying the other two conditions completes the proof.Proposition 4.
The equilibrium pointO may be destabilized due to(i)
flip bifurcation ifk1k2<4/c1+r1c2+r2(ii)
Neimark–Sacker bifurcationk1+k2<4c1+r1+c2+r2/c1+r1c2+r2Proof.
Substituting (13) in g−1 and Ω then taking g−1=0, we get(14)c1+r1c2+r2k1k2=4c1+r1c2+r2k1+k2c1+r1+c2+r2.
Substituting (14) in Ω>0 completes the proof of (i). For the part (ii), we put Ω=0, and then we get(15)c1+r1c2+r2k1k2=8c1+r1c2+r2k1+k2c1+r1+c2+r2.
Then substituting (15) in g−1>0 completes the proof of (ii).


### 2.2. The Synchronization Property

Synchronized trajectories in such game is an important property that may give more information about the dynamic of the model’s game. Indeed, such property may be occurred and there may be a transversely stable orbit on the diagonalΔ=q1,q2:q1=q2. In case, such dynamics exist due to nonsynchronizing trajectories, and it would be important to detect the initial conditions leading to synchronization property. The possibility of such property arises when an invariant one-dimensional submanifold of ℝ2 exists. This means that the synchronized trajectories are described by(16)q1t,q2t=Tt:q10,q20|q1t=q2t∀t≥0.These trajectories given in (16) are regularized by the restriction of the map T on the invariant submanifold by which synchronized dynamics occur and are described by the one-dimensional map:(17)T|ΔΔ⟶Δ.For trajectories beginning outside the map given in (17) are said to be synchronized ones if q1t−q2t⟶0 as t⟶∞. Our map given in (7) possesses the six parameters, k1,k2,c1,c2,r1, and r2. Now, we assume the following:(18)k1=k2=k,c1=c2=candr1=r2=r.Under this assumption, the map (7) becomes(19)Tsq1,q2:q1t+1=q1t+kq1tq2Q2−c−r,q2t+1=q2t+kq2tq1Q2−c−r.The restrictionTs|ΔΔ⟶Δ conjugates the following one-dimensional map:(20)q′=q+kq14q−c−r.Studying the transverse stability of the synchronized attractors of the system (19) requires to calculate its Jacobian at the equilibrium point Os=1/4c+r,1/4c+r as follows:(21)JTsq1,q2=1−kc+r001−kc+r,Whose eigenvalues areλ1=λ2=1−kc+r. It is simple to see that λ1,2<1 if k<2/c+r, and hence Os is an attracting stable node.


### 2.3. Local Analysis

In this section, we give some numerical experiments to validate the above results. As in [6], we assume the following values, c1=2,c2=1.5,r1=0.13,r2=0.15,k1=0.2, and k2=0.1. The Jacobian (10) at these values becomes(22)JO≈0.6281−0.41910.02710.8141,whose eigenvalues are real and equal λ1≈0.6344 and λ2≈0.8077. One can see that λ1,2<1, and hence the equilibrium point O=0.2194787379,0.1508916323 is the local stable point. Any increase in the parameters r1,r2,k1, and k2 leads to unstable equilibrium point through flip bifurcation. In Figure 1(a), we take the parameter k1 as the bifurcation parameter and fix the other parameters’ values to c1=2,c2=1.5,r1=0.13,r2=0.15, and k2=0.1. As it can be seen that the point O is locally stable for all the values of k1 till this parameter reaches the point of cycle of period two and it becomes unstable. As k1 increases further, higher periodic cycles are born, then the map enters the chaos area, and then becomes chaotic (as confirmed in the Largest Lyapunov exponent (LLE) given in Figure 1(c)). We should highlight here that the governmental taxes are different with 13% imposed on the quantity produced by the first firm while the tax imposed on the second firm is 15%. Numerical experiments show that higher governmental taxes imposed lead to shrinking of the stability region as will be discussed in the global analysis later on. Figure 1(b) shows that the impact on the parameter k1 is slightly different on the quantities produced by both firms. This may be due to the cost of productions and taxes imposed that are slightly different. The same discussions are for the impact of the parameter k2 when it is taken as the bifurcation diagram. Figures 1(d)–1(f) show the bifurcation diagram with respect to k2 and the corresponding LLE at the parameters’ values, c1=2,c2=1.5,r1=0.13,r2=0.15, and k1=0.2. On the other hand, we give in Figure 1(g)–1(l) the taxation impact on the equilibrium point and the corresponding LLE. Figure 1(g) shows that at the parameters’ values c1=2,c2=1.5,k1=1.1,k=0.9, and r2=0.15, the equilibrium point is stable at taxation rate r1<37%; then, it becomes unstable through Neimark–Sacker bifurcation. As r1 increases further, a closed ring is born that is followed by a cycle of period 7 and then the map enters the chaotic region. The same discussion and observations are for the other taxation parameter r2 on where its effect is given in Figures 1(j) and 1(k). The Figures 1(j) and 1(k) show the impact of r2 on the quantity q2. Simulation shows that in the interval 0.256756,0.26667, the quantity q2 bifurcates into two chaotic bands not period 2-cycle, and this is observed in the corresponding Lyapunov exponent given in Figure 1(l).Figure 1
(a–l) Different types of bifurcation diagrams with respect to the parametersr1,r2,k1, and k2 when r1<r2 and their corresponding Lyapunov exponents. Each figure is plotted at some parameters values given in the text and within the figure itself.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Now, we assume the set of parameters’ valuesr1>r2, c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1, and k2=1. The Jacobian [10] at these values becomes(23)JO≈−0.93671−0.236710.31329−0.93671,Whose eigenvalues are complex and equalλ1,2≈−0.93671±0.27232i. One can see that λ1,2<1, and hence the equilibrium point O=0.1089568979,0.1442076590 is local stable point. Any increase in the parameters r1,r2,k1, and k2 leads to unstable equilibrium point through Neimark–Sacker (N-S) bifurcation. In Figures 2(a) and 2(b), we take the parameter k1 as the bifurcation parameter and fix the other parameters’ values to c1=2,c2=1.5,r1=0.25,r2=0.20, and k2=1. It is clear that the equilibrium point loses its stability through N-S bifurcation. The same observations are obtained when we take the parameters r1,r2, and k2 as the bifurcation’s parameter. The Figures 2(c)–2(h) give the N-S bifurcations for those parameters.Figure 2
(a–h) The N-S bifurcation diagrams with respect to the parametersr1,r2,k1, and k2 when r1>r2. Each figure is plotted at some parameters values given in the text and within the figure itself.
(a)(b)(c)(d)(e)(f)(g)(h)


### 2.4. Global Analysis

The above local analysis regarding the bifurcation diagrams does not permit further understanding about the future evolution of the map (7). For this reason, some complex behaviors and their basins of attraction are given here. As said before, the map (7) depends on the parameters k1,k2,c1,c2,r1, and r2. Let us first investigate the influences of the parameters k1 and k2 by assuming the following parameters values, c1=2,c2=1.5,r1=0.13,r2=0.15. This set of parameters values gives the two-dimensional bifurcation diagram in the k1,k2−plane presented by Figure 3(a). It is clear that this plane is divided into different periods by different colors besides the stability region. The figure shows that the system behaves chaotically when the speed of adjustment parameters k1 and k2 are all high. Any changing in those two parameters while keeping the other parameters fixed makes the system enter the chaotic region through Flip and Neimark Sacker bifurcation. The system’s attractors and multistability have been investigated by many scholars in literature ([4, 6, 16]). Such investigations have shown one of the important conditions of chaos which is the sensitivity to the initial conditions. In economic market, different strategies by decision makers may lead to different directions to the development of the firm’s direction. Consequently, adopting a good way of choosing initial conditions is very important which in turn gives an indication of whether the firm will behave well or not in the future. This includes the coexistence of multiple attractors and their influences on the stability of the equilibrium point. The basins of attraction can help in investigating such attractors. It can be used to analyze the convergence of system after a series of iterations (games) based on certain initial conditions. It helps the firm’s decision makers to choose a range of initial values by which the firm can develop better. For instance, let us assume the following initial set: c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422, and k2=0.66044. This set gives rise to a period of 3-cycle with its basins of attraction plotted in Figure 3(b). In Figure 3(b), the gray color characterizes the divergence (or the unfeasible trajectories), while the other colors represent the attractive basin of the equilibrium point O and the basin of attraction of period 3-cycle coexisting with O. The peculiar shape of the basin of attraction given in Figure 3(b) will be because of the origin point as will be discussed in the next section. Figure 3(c) gives a chaotic situation of the system around the equilibrium point as k1 increases to 1.26, and the other parameters values are fixed to c1=2,c2=1.5,r1=0.13,r2=0.15,k2=0.1. Increasing the taxes and the speed of adjustment parameters gives rise to instability situation of the equilibrium point due to Neimark–Sacker bifurcation. For instance, the parameter values c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.027, and k2=0.1 present a quasiperiodic dynamic possessing an attracting invariant closed ring around the equilibrium point. Other closed invariant curve is given when increasing the speed of adjustment parameter k1 to 1.1. As shown previously, the stability conditions depend on the tax parameters and any changes in those parameters give rise to instability of the equilibrium point. For this reason, the two-dimensional bifurcation diagram in the r1,r2−plane is depicted in Figure 3(f). We finish this section by giving in Figure 4, an example of period 2-cycle at fixing the parameters values: c1=2,c2=1.5,r1=0.027,r2=0.448,k1=1.1, and k2=0.9.Figure 3
(a) 2D-bifurcation diagram atc1=2,c2=1.5,r1=0.13,r2=0.15. (b) The period of 3-cycle at c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422, and k2=0.66044. (c) Chaotic dynamic behavior at c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.26,k2=0.1. (d) Closed invariant curve at c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.027,k2=0.1. (e) Closed ring at c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.1,k2=0.1. (f) 2D-bifurcation diagram at c1=2,c2=1.5,k1=1.1,k2=0.9.
(a)(b)(c)(d)(e)(f)Figure 4
The period of 2-cycle atc1=2,c2=1.5,r1=0.072,r2=0.448,k1=1.1, and k2=0.9.


### 2.5. Critical Curves and Preimage Zones

The phase plane of the map (7) at any sets of parameters’ values possesses some important characteristics. Looking at this map, one can observe that at q1t=0 or q2t=0, one gets q1t+1=0 or p2t+1=0, respectively. This means that the origin point 0,0 traps the map (7), and therefore, it is important to start with this point in order to calculate the boundaries of the attractive basins for any attracting set for that map. To do that, we set q1t+1=q´1 and q2t+1=q´2 in (7) where ‘indicates time evolution. This means that the map’s time evolution is attained by the iteration of T:q1,q2⟶q´1,q´2 as follows:(24)T:q´1=q1+k1q1q2Q2−c1−r1,q´2=q2+k2q2q1Q2−c2−r2.If the inverse mapT−1:q´1,q´2⟶q1,q2 is existed and gets unique values in each point in the range, then the map is called invertible. Hence, the point q´1,q´2∈ℝ2 is called a rank-1 image while q1,q2 is called rank-1 preimages. If for an image q´1,q´2, there are at least two rank-1 preimages and we call T a noninvertible map. The following proposition shows that the point 0,0 has two real rank-1 preimages, and hence, the map is noninvertible.Proposition 5.
The origin point0,0 has two real rank-1 preimages.Proof.
Settingq´1=0 and q´2=0 in (23), one gets a system of algebraic equations whose solutions are the two preimages as follows:(25)O−10=0,0,O−11=k12k2c2k2+r2k2−1k1k2c1+c2+r1+r2−k1−k22,k1k22c1k1+r1k1−1k1k2c1+c2+r1+r2−k1−k22.
This completes the proof.
Figure5 shows that phase plane of the map is divided by the two regions that are Z1 and Z3, and hence, the map (7) is of type Z1−Z3. These two regions are calculated by the critical curves LC and LC−1. The later can be calculated by vanishing the Jacobian of (7) which gives(26)ℓ1ℓ2q23+ℓ1ℓ2q13+k2ℓ1q12+k1ℓ2q22+3ℓ1ℓ2q1q22+3ℓ1ℓ2q12q2+ℓ3q1q2=0,ℓi=1−ciki−riki;i=1,2ℓ3=k1k2c1+c2+r1+r2−k1−k2.
BecauseLC is calculated from the relation LC=TLC−1, and (26) contains many parameters, it is plotted at the parameter values c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1,k2=0.9. It is important to highlight that the map (7) is defined at any point in the phase plane except the line q2=−q1 and its preimages of any order. Particularly, there is point in the phase plane where the map (7) takes the form 0/0. This point is called a focal point, and in our case, it is the origin point 0,0. It is responsible for constructing the lobes found in the basins of attraction of any attracting set. These lobes have been discussed in detail in [28]. Such lobes appear in Figure 3(b), and we highlight them in Figure 6.Figure 5
The critical curvesLC and LC−1 and the region Zi, i=1,3 at the set of parameters: c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1,k2=0.9.Figure 6
The Lobes constructed with period of 3-cycle at:c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422 and k2=0.66044.


## 2.1. The Equilibrium Point and Its Stability

The equilibrium point of the above map is a unique Nash equilibrium point and is obtained by settingq1t+1=q1t=q¯1 and q2t+1=q2t=q¯2. It has the following form:(8)O=q¯1,q¯2=c2+r2c1+c2+r1+r22,c1+r1c1+c2+r1+r22,which is a positive point. We should highlight here that the equilibrium point (8) is the same as the one given in [28] but when the model has free taxes r1=r2=0. The following propositions are used in the sequel to classify the type of the point O.Proposition 1.
Let the Jacobian matrix the map (7) at any equilibrium point possesses two eigenvalues λ1 and λ2. Then, this equilibrium point can be classified according to the following:(i)
Ifλ1,2<1, then it becomes an attracting node and it is stable(ii)
Ifλ1,2>1, then it is an unstable repelling node(iii)
Ifλ1<1 and λ2>1, (or λ1>1 and λ2<1) then it is an unstable saddle point(iv)
It is a nonhyperbolic point ifλ1=1 and λ2≠1 (or λ1≠1 and λ2=1)
In case, the eigenvalues have a complicated analytical form, we use the following proposition.Proposition 2.
Letτ and δ be the trace and determinant of the Jacobian matrix of (7), and suppose the following(9)g1=1−τ+δ,g−1=1+τ+δ,Ω=1−δ.
Then, we have the following properties:(i)
The equilibrium point is locally asymptotically stable ifg1>0,g−1>0, and Ω>0.(ii)
The equilibrium point becomes unstable through a flip bifurcation ifg1>0,g−1=0, and Ω>0. In this case, the two eigenvalues are real and pass through −1.(iii)
The equilibrium point becomes unstable through a transcritical or fold bifurcation ifg1>0,g−1>0, and Ω>0. In this case, the two eigenvalues are real and pass through 1.(iv)
The equilibrium point becomes unstable through Neimark–Sacker ifg1>0,g−1>0, and Ω<0. In this case, the two eigenvalues are complex and their modulus passes through 1.
The Jacobian matrix for the map (7) at the equilibrium point O becomes(10)JO=1−2c1+r1c2+r2k1c1+r1+c2+r2−c2+r2c1+r1−c2−r2k1c1+r1+c2+r2c1+r1c1+r1−c2−r2k2c1+r1+c2+r21−2c1+r1c2+r2k2c1+r1+c2+r2.
And its eigenvalues are(11)λ1,2=1−c1+r1c2+r2k1+k2c1+r1+c2+r2±±2c1+r1c2+r2c1+r1+c2+r2k12+k22−c12+r12+c22+r22+2c1r1+2c2r2k1k2c1+r1c2+r2.Proposition 3.
The equilibrium pointO is locally asymptotically stable if the following conditions are satisfied:(12)c1+r1c2+r2k1k2−4c1+r1c2+r2c1+r1+c2+r2k1+k2+4>0,2c1+r1c2+r2c1+r1+c2+r2k1+k2−c1+r1c2+r2k1k2>0.Proof.
Simple calculations show thatτ and δ for the Jacobian given in (10) become(13)τ=21−c1+r1c2+r2k1+k2c1+r1+c2+r2,δ=1−2c1+r1c2+r2k1+k2c1+r1+c2+r2+c1+r1c2+r2k1k2.
Substituting (13) in (9), we get g1=c1+r1c2+r2k1k2 which is always positive since ci,ri and ki,i=1,2 are positive parameters. Simplifying the other two conditions completes the proof.Proposition 4.
The equilibrium pointO may be destabilized due to(i)
flip bifurcation ifk1k2<4/c1+r1c2+r2(ii)
Neimark–Sacker bifurcationk1+k2<4c1+r1+c2+r2/c1+r1c2+r2Proof.
Substituting (13) in g−1 and Ω then taking g−1=0, we get(14)c1+r1c2+r2k1k2=4c1+r1c2+r2k1+k2c1+r1+c2+r2.
Substituting (14) in Ω>0 completes the proof of (i). For the part (ii), we put Ω=0, and then we get(15)c1+r1c2+r2k1k2=8c1+r1c2+r2k1+k2c1+r1+c2+r2.
Then substituting (15) in g−1>0 completes the proof of (ii).


## 2.2. The Synchronization Property

Synchronized trajectories in such game is an important property that may give more information about the dynamic of the model’s game. Indeed, such property may be occurred and there may be a transversely stable orbit on the diagonalΔ=q1,q2:q1=q2. In case, such dynamics exist due to nonsynchronizing trajectories, and it would be important to detect the initial conditions leading to synchronization property. The possibility of such property arises when an invariant one-dimensional submanifold of ℝ2 exists. This means that the synchronized trajectories are described by(16)q1t,q2t=Tt:q10,q20|q1t=q2t∀t≥0.These trajectories given in (16) are regularized by the restriction of the map T on the invariant submanifold by which synchronized dynamics occur and are described by the one-dimensional map:(17)T|ΔΔ⟶Δ.For trajectories beginning outside the map given in (17) are said to be synchronized ones if q1t−q2t⟶0 as t⟶∞. Our map given in (7) possesses the six parameters, k1,k2,c1,c2,r1, and r2. Now, we assume the following:(18)k1=k2=k,c1=c2=candr1=r2=r.Under this assumption, the map (7) becomes(19)Tsq1,q2:q1t+1=q1t+kq1tq2Q2−c−r,q2t+1=q2t+kq2tq1Q2−c−r.The restrictionTs|ΔΔ⟶Δ conjugates the following one-dimensional map:(20)q′=q+kq14q−c−r.Studying the transverse stability of the synchronized attractors of the system (19) requires to calculate its Jacobian at the equilibrium point Os=1/4c+r,1/4c+r as follows:(21)JTsq1,q2=1−kc+r001−kc+r,Whose eigenvalues areλ1=λ2=1−kc+r. It is simple to see that λ1,2<1 if k<2/c+r, and hence Os is an attracting stable node.


## 2.3. Local Analysis

In this section, we give some numerical experiments to validate the above results. As in [6], we assume the following values, c1=2,c2=1.5,r1=0.13,r2=0.15,k1=0.2, and k2=0.1. The Jacobian (10) at these values becomes(22)JO≈0.6281−0.41910.02710.8141,whose eigenvalues are real and equal λ1≈0.6344 and λ2≈0.8077. One can see that λ1,2<1, and hence the equilibrium point O=0.2194787379,0.1508916323 is the local stable point. Any increase in the parameters r1,r2,k1, and k2 leads to unstable equilibrium point through flip bifurcation. In Figure 1(a), we take the parameter k1 as the bifurcation parameter and fix the other parameters’ values to c1=2,c2=1.5,r1=0.13,r2=0.15, and k2=0.1. As it can be seen that the point O is locally stable for all the values of k1 till this parameter reaches the point of cycle of period two and it becomes unstable. As k1 increases further, higher periodic cycles are born, then the map enters the chaos area, and then becomes chaotic (as confirmed in the Largest Lyapunov exponent (LLE) given in Figure 1(c)). We should highlight here that the governmental taxes are different with 13% imposed on the quantity produced by the first firm while the tax imposed on the second firm is 15%. Numerical experiments show that higher governmental taxes imposed lead to shrinking of the stability region as will be discussed in the global analysis later on. Figure 1(b) shows that the impact on the parameter k1 is slightly different on the quantities produced by both firms. This may be due to the cost of productions and taxes imposed that are slightly different. The same discussions are for the impact of the parameter k2 when it is taken as the bifurcation diagram. Figures 1(d)–1(f) show the bifurcation diagram with respect to k2 and the corresponding LLE at the parameters’ values, c1=2,c2=1.5,r1=0.13,r2=0.15, and k1=0.2. On the other hand, we give in Figure 1(g)–1(l) the taxation impact on the equilibrium point and the corresponding LLE. Figure 1(g) shows that at the parameters’ values c1=2,c2=1.5,k1=1.1,k=0.9, and r2=0.15, the equilibrium point is stable at taxation rate r1<37%; then, it becomes unstable through Neimark–Sacker bifurcation. As r1 increases further, a closed ring is born that is followed by a cycle of period 7 and then the map enters the chaotic region. The same discussion and observations are for the other taxation parameter r2 on where its effect is given in Figures 1(j) and 1(k). The Figures 1(j) and 1(k) show the impact of r2 on the quantity q2. Simulation shows that in the interval 0.256756,0.26667, the quantity q2 bifurcates into two chaotic bands not period 2-cycle, and this is observed in the corresponding Lyapunov exponent given in Figure 1(l).Figure 1
(a–l) Different types of bifurcation diagrams with respect to the parametersr1,r2,k1, and k2 when r1<r2 and their corresponding Lyapunov exponents. Each figure is plotted at some parameters values given in the text and within the figure itself.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Now, we assume the set of parameters’ valuesr1>r2, c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1, and k2=1. The Jacobian [10] at these values becomes(23)JO≈−0.93671−0.236710.31329−0.93671,Whose eigenvalues are complex and equalλ1,2≈−0.93671±0.27232i. One can see that λ1,2<1, and hence the equilibrium point O=0.1089568979,0.1442076590 is local stable point. Any increase in the parameters r1,r2,k1, and k2 leads to unstable equilibrium point through Neimark–Sacker (N-S) bifurcation. In Figures 2(a) and 2(b), we take the parameter k1 as the bifurcation parameter and fix the other parameters’ values to c1=2,c2=1.5,r1=0.25,r2=0.20, and k2=1. It is clear that the equilibrium point loses its stability through N-S bifurcation. The same observations are obtained when we take the parameters r1,r2, and k2 as the bifurcation’s parameter. The Figures 2(c)–2(h) give the N-S bifurcations for those parameters.Figure 2
(a–h) The N-S bifurcation diagrams with respect to the parametersr1,r2,k1, and k2 when r1>r2. Each figure is plotted at some parameters values given in the text and within the figure itself.
(a)(b)(c)(d)(e)(f)(g)(h)


## 2.4. Global Analysis

The above local analysis regarding the bifurcation diagrams does not permit further understanding about the future evolution of the map (7). For this reason, some complex behaviors and their basins of attraction are given here. As said before, the map (7) depends on the parameters k1,k2,c1,c2,r1, and r2. Let us first investigate the influences of the parameters k1 and k2 by assuming the following parameters values, c1=2,c2=1.5,r1=0.13,r2=0.15. This set of parameters values gives the two-dimensional bifurcation diagram in the k1,k2−plane presented by Figure 3(a). It is clear that this plane is divided into different periods by different colors besides the stability region. The figure shows that the system behaves chaotically when the speed of adjustment parameters k1 and k2 are all high. Any changing in those two parameters while keeping the other parameters fixed makes the system enter the chaotic region through Flip and Neimark Sacker bifurcation. The system’s attractors and multistability have been investigated by many scholars in literature ([4, 6, 16]). Such investigations have shown one of the important conditions of chaos which is the sensitivity to the initial conditions. In economic market, different strategies by decision makers may lead to different directions to the development of the firm’s direction. Consequently, adopting a good way of choosing initial conditions is very important which in turn gives an indication of whether the firm will behave well or not in the future. This includes the coexistence of multiple attractors and their influences on the stability of the equilibrium point. The basins of attraction can help in investigating such attractors. It can be used to analyze the convergence of system after a series of iterations (games) based on certain initial conditions. It helps the firm’s decision makers to choose a range of initial values by which the firm can develop better. For instance, let us assume the following initial set: c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422, and k2=0.66044. This set gives rise to a period of 3-cycle with its basins of attraction plotted in Figure 3(b). In Figure 3(b), the gray color characterizes the divergence (or the unfeasible trajectories), while the other colors represent the attractive basin of the equilibrium point O and the basin of attraction of period 3-cycle coexisting with O. The peculiar shape of the basin of attraction given in Figure 3(b) will be because of the origin point as will be discussed in the next section. Figure 3(c) gives a chaotic situation of the system around the equilibrium point as k1 increases to 1.26, and the other parameters values are fixed to c1=2,c2=1.5,r1=0.13,r2=0.15,k2=0.1. Increasing the taxes and the speed of adjustment parameters gives rise to instability situation of the equilibrium point due to Neimark–Sacker bifurcation. For instance, the parameter values c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.027, and k2=0.1 present a quasiperiodic dynamic possessing an attracting invariant closed ring around the equilibrium point. Other closed invariant curve is given when increasing the speed of adjustment parameter k1 to 1.1. As shown previously, the stability conditions depend on the tax parameters and any changes in those parameters give rise to instability of the equilibrium point. For this reason, the two-dimensional bifurcation diagram in the r1,r2−plane is depicted in Figure 3(f). We finish this section by giving in Figure 4, an example of period 2-cycle at fixing the parameters values: c1=2,c2=1.5,r1=0.027,r2=0.448,k1=1.1, and k2=0.9.Figure 3
(a) 2D-bifurcation diagram atc1=2,c2=1.5,r1=0.13,r2=0.15. (b) The period of 3-cycle at c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422, and k2=0.66044. (c) Chaotic dynamic behavior at c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.26,k2=0.1. (d) Closed invariant curve at c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.027,k2=0.1. (e) Closed ring at c1=2,c2=1.5,r1=0.25,r2=0.20,k1=1.1,k2=0.1. (f) 2D-bifurcation diagram at c1=2,c2=1.5,k1=1.1,k2=0.9.
(a)(b)(c)(d)(e)(f)Figure 4
The period of 2-cycle atc1=2,c2=1.5,r1=0.072,r2=0.448,k1=1.1, and k2=0.9.


## 2.5. Critical Curves and Preimage Zones

The phase plane of the map (7) at any sets of parameters’ values possesses some important characteristics. Looking at this map, one can observe that at q1t=0 or q2t=0, one gets q1t+1=0 or p2t+1=0, respectively. This means that the origin point 0,0 traps the map (7), and therefore, it is important to start with this point in order to calculate the boundaries of the attractive basins for any attracting set for that map. To do that, we set q1t+1=q´1 and q2t+1=q´2 in (7) where ‘indicates time evolution. This means that the map’s time evolution is attained by the iteration of T:q1,q2⟶q´1,q´2 as follows:(24)T:q´1=q1+k1q1q2Q2−c1−r1,q´2=q2+k2q2q1Q2−c2−r2.If the inverse mapT−1:q´1,q´2⟶q1,q2 is existed and gets unique values in each point in the range, then the map is called invertible. Hence, the point q´1,q´2∈ℝ2 is called a rank-1 image while q1,q2 is called rank-1 preimages. If for an image q´1,q´2, there are at least two rank-1 preimages and we call T a noninvertible map. The following proposition shows that the point 0,0 has two real rank-1 preimages, and hence, the map is noninvertible.Proposition 5.
The origin point0,0 has two real rank-1 preimages.Proof.
Settingq´1=0 and q´2=0 in (23), one gets a system of algebraic equations whose solutions are the two preimages as follows:(25)O−10=0,0,O−11=k12k2c2k2+r2k2−1k1k2c1+c2+r1+r2−k1−k22,k1k22c1k1+r1k1−1k1k2c1+c2+r1+r2−k1−k22.
This completes the proof.
Figure5 shows that phase plane of the map is divided by the two regions that are Z1 and Z3, and hence, the map (7) is of type Z1−Z3. These two regions are calculated by the critical curves LC and LC−1. The later can be calculated by vanishing the Jacobian of (7) which gives(26)ℓ1ℓ2q23+ℓ1ℓ2q13+k2ℓ1q12+k1ℓ2q22+3ℓ1ℓ2q1q22+3ℓ1ℓ2q12q2+ℓ3q1q2=0,ℓi=1−ciki−riki;i=1,2ℓ3=k1k2c1+c2+r1+r2−k1−k2.
BecauseLC is calculated from the relation LC=TLC−1, and (26) contains many parameters, it is plotted at the parameter values c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1,k2=0.9. It is important to highlight that the map (7) is defined at any point in the phase plane except the line q2=−q1 and its preimages of any order. Particularly, there is point in the phase plane where the map (7) takes the form 0/0. This point is called a focal point, and in our case, it is the origin point 0,0. It is responsible for constructing the lobes found in the basins of attraction of any attracting set. These lobes have been discussed in detail in [28]. Such lobes appear in Figure 3(b), and we highlight them in Figure 6.Figure 5
The critical curvesLC and LC−1 and the region Zi, i=1,3 at the set of parameters: c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1,k2=0.9.Figure 6
The Lobes constructed with period of 3-cycle at:c1=2,c2=1.5,r1=0.13,r2=0.15,k1=1.39422 and k2=0.66044.


## 3. Conclusion

In the present paper, we have introduced a Cournot duopoly game whose players are rational and seek the optimality of quantity production in order to achieve profit maximization. The utility function introduced in this paper may be considered as a special case of the Singh and Vives function [1]. Our utility function has given inverse demand functions which are the same as those obtained from Cobb–Douglas utility. The map’s game possesses a unique equilibrium point which is Nash point which loses its stability through flip and Neimark–Sacker bifurcation. As shown in the obtained results, increasing the tax parameters gives rise to Neimark–Sacker bifurcation, then route to chaos is obtained, and then unpredictable behaviors for the dynamics of the game are arisen. It has been observed that the high values for the tax parameters affect the stability region. The region has decreased for high values of the taxes and the speed of adjustment parameters. Furthermore, the attractive basins for some attracting sets have been presented. The shapes of those basins are peculiar because they have contained some lobes due to the origin focal point. Future study directions will focus on the application of the introduced utility in economic models of public enterprises with consumer surplus.



---

*Source: 5195337-2022-04-14.xml*