Journal of Computational and Applied Mathematics 218 (2008) 247–258
www.elsevier.com/locate/cam
Stable cycles in a Cournot duopoly model of Kopel
W. Govaertsa,∗
, R. Khoshsiar Ghaziania,b
aDepartment of Applied Mathematics and Computer Science, Ghent University, Krijgslaan 281-S9, B-9000 Gent, Belgium
bDepartment of Mathematics, Faculty of Science, Shahrekord University, PO Box 115, Shahrekord, Iran
Received 17 August 2006; received in revised form 13 November 2006
Abstract
We consider a discrete map proposed by M. Kopel that models a nonlinear Cournot duopoly consisting of a market structure
between the two opposite cases of monopoly and competition. The stability of the fixed points of the discrete dynamical system
is analyzed. Synchronization of two dynamics parameters of the Cournot duopoly is considered in the computation of stability
boundaries formed by parts of codim-1 bifurcation curves. We discover more on the dynamics of the map by computing numerically
the critical normal form coefficients of all codim-1 and codim-2 bifurcation points and computing the associated two-parameter
codim-1 curves rooted in some codim-2 points. It enables us to compute the stability domains of the low-order iterates of the map.
We concentrate in particular on the second, third and fourth iterates and their relation to the period doubling, 1:3 and 1:4 resonant
Neimark–Sacker points.
© 2007 Elsevier B.V. All rights reserved.
MSC: 39A10
Keywords: Bifurcation of fixed points; Cycles; Normal form coefficients; Bistability
1. Introduction
The first well-known model which gives a mathematical description of competition in a duopoly market dates back
to the French economist Cournot [13] with the highlighted characteristics:
• Competing firms produce goods that are perfect substitutes.
• Both firms must consider the actions and reactions of the competitor.
• Each firm forms expectations of the other firm’s output in order to determine a profit maximizing quantity to produce
in the next time period (this situation is called strategic interdependence).
The model that he presented has been much studied for its ability to generate complex dynamics and also because of
its more general foreshadowing of game theory. It has often been noted that the Cournot equilibrium is but a special
case of the Nash Equilibrium [21], the more general formulation used by modern industrial organization economists
in studying oligopoly theory.
∗ Corresponding author.
E-mail address: Willy.Govaerts@UGent.be (W. Govaerts).
0377-0427/$ - see front matter © 2007 Elsevier B.V. All rights reserved.
doi:10.1016/j.cam.2007.01.012
248 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
Recently, several works have shown that the Cournot model may lead to complex behaviors such as cyclic and
chaotic, see, for example [1,17,22–24]. Among the first to do this was Puu [22,23] who found a variety of complex
dynamics arising in the Cournot duopoly case including the appearance of attractors with fractal dimension. Dynamics
of a Cournot game by players with bounded rationality has been studied in [5]. Local stability of a duopoly game with
heterogeneous expectation has been studied in [4]. Multistability, cyclic and chaotic behavior of a Cournot game have
been studied in [9], where in the model the reaction functions have the form of the logistic map. Some preliminary
results on the local bifurcations of a Kopel map were obtained in [17]. Explicit boundaries of local stability of the
fixed point of a Kopel map have been derived in [2]. Basins of attraction in a Kopel map have been studied in [8].
Other studies on the dynamics of oligopoly models with more firms and other modifications include Ahmed and
Agiza [6], Agiza [1] and Agiza et al. [3]. The development of complex oligopoly dynamics theory has been reviewed
in [25].
In this paper we consider the general case of a duopoly model, see [2], introduced in [17] with positive adjustment
coefficient . The main aim is to investigate the overall dynamic behavior of the model when > 0 and to compute
stability domains of the first, second, third and fourth iterates of the map. In Section 2 we introduce the model and discuss
the general stability and branching of the fixed points. In particular, we compute analytically the critical normal form
coefficients in the case of period doubling bifurcations to reveal sub- or supercriticality. In Section 3 we concentrate on
the economically relevant case 1 and numerically compute curves of codim-1 bifurcations and the critical normal
form coefficients of codim-2 bifurcation points, using the new software CL_MATCONTM [14,16]. These tools enable us
to compute stability boundaries of 2, 3 and 4-cycles. Furthermore, by considering the critical normal form coefficients
of the R4 resonance point, we determine the bifurcation scenario of the map near this point.
In Section 4 we briefly describe R3 and R2 bifurcation points in the region > 1 which is of no interest for the
economic model. In Section 5 we summarize our results and draw some conclusions.
2. The map and the local stability analysis of its fixed points
The model that we use is a two-dimensional map described in [17,2]. Two firms are homogeneous with regard
to their expectation formation and the action effect on each other. The duopoly Kopel model assumes that at each
discrete time t the two firms produce the quantities x1(t) and x2(t), respectively, and decide their productions for the
next period x1(t + 1) and x2(t + 1). The time evolution of the model is determined by the two-dimensional map
T : (x1(t), x2(t)) → (x1(t + 1), x2(t + 1)) defined by
T :
x1(t + 1) = (1 − )x1(t) + x2(t)(1 − x2(t)),
x2(t + 1) = (1 − )x2(t) + x1(t)(1 − x1(t)),
(1)
where and are two model parameters. The positive parameter measures the intensity of the effect that one firm’s
actions has on the other firm. Firms do not change their productions according to the computed optimal productions
(i.e., the ‘logistic’ reaction functions) but they prefer to choose a weighted average between the previous production
and the computed one, with weights 1 − and , respectively; is called the adjustment coefficient. The meaning of
the model implies that the parameter ∈ [0, 1]. However, it is best to ignore this restriction in a first global study of
the properties of the model, cf. [2].
2.1. Fixed points of the map
Bifurcation of maps have been studied intensively in the literature, cf [20,12,11,10]. A comprehensive discussion is
given in [18]. We further use the recent results in [16,19].
The fixed points of (1) and their stability were studied analytically in [2, Section 2.1–4]. We summarize the obtained
results briefly. For = 0, the fixed points of (1) are the solutions to
x∗
1 = x∗
2 (1 − x∗
2 ), x∗
2 = x∗
1 (1 − x∗
1 ). (2)
Besides the trivial solution E1 : (x∗
1 , x∗
2 ) = (0, 0), a positive symmetric fixed point exists for > 1, given by E2 :
(x∗
1 , x∗
2 ) = (( − 1)/ , ( − 1)/ ).
W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258 249
0 1 2 3 4 5 6 7 8
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
μ
ρ
ΩS
(E1) ΩS
(E2) ΩS
(E3,4)
Fig. 1. Stable regions of Ei, i = 1, 2, 3, 4.
Two further nonsymmetric Nash Equilibria, given by
E3 : (x∗
1 , x∗
2 ) =
+ 1 +
√
( + 1)( − 3)
2
,
+ 1 −
√
( + 1)( − 3)
2
, (3)
and its (x1, x2) → (x2, x1) reflection E4, exist for 3.
The study of the local stability of fixed points is based on the linearization of (1). In an equilibrium point the Jacobian
J(x1, x2) of (1) has the eigenvalues:
1,2 = (1 − ) ± (1 − 2x1)(1 − 2x2). (4)
Depending on the values of x1 and x2, these may be real or form a conjugate complex pair. A fixed point of (1) is
stable if
| j | < 1, j = 1, 2. (5)
Proposition 2.1. The equilibrium solution E1 is asymptotically stable for ( , ) ∈ S
(E1) where
S
(E1) = ( , )|0 < < 1, 0 < < 1( ) =
2
1 +
.
It loses stability via a flip bifurcation when crossing the threshold 1( ), 0 < < 1 and via branching along = 1.
Proof. The stability boundaries of E1 can be derived by imposing the stability conditions (5). These boundaries were
computed in [2] and are presented in Fig. 1 ( S
(E1)).
What remains to be proved is that E1 loses its stability and bifurcates to a new branch of fixed points at = 1. To
do this we show that the discriminant of the algebric branching equation (ABE), see [16], is positive. We consider the
Jacobian matrix FX = [Tx − I|T ], evaluated in E1 that is
− 0
− 0
. (6)
250 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
This matrix is clearly rank deficient along = 1. We first compute vectors 
1, 
2 and which form a basis for the
null space of N([Tx − I|T ]) and N([Tx − I|T ]∗), respectively. Possible choices are

1 =
1
√
2
,
1
√
2
, 0
T
, 
2 = (0, 0, 1)T
, =
1
√
2
,
1
√
2
T
.
Now we consider the ABE:
c11
2
+ 2c12 + c22
2
= 0, (7)
where cjk = , F0
YY 
j 
k , for j, k = 1, 2. Here the 2 × 3 × 3 tensor F0
YY is given by
F0
YY (:, :, 1) =
0 0 0
−2 0 (1 − x1) − x1
, (8)
F0
YY (:, :, 2) =
0 −2 (1 − x2) − x2
0 0 0
, (9)
F0
YY (:, :, 3) =
0 (1 − x2) − x2 0
(1 − x1) − x1 0 0
. (10)
We now obtain c11 = −
√
2 , c12 = , c22 = 0. So the discriminant of (7), c2
12 − c11c22 = 2 > 0 is clearly positive.
This shows that we have a branch point when = 1.
Proposition 2.2. E2 is asymptotically stable for ( , ) ∈ S
(E2) = S
(E21) ∪ S
(E22) where
S
(E21) = ( , )|1 < < 2, 0 < < 21( ) =
2
3 −
,
and
S
(E22) = ( , )|2 < < 3, 0 < < 22( ) =
2
− 1
.
It loses stability via a flip bifurcation point on the boundaries of
(i) = 21( ), 1 < < 2.
(ii) = 22( ), 2 < < 3.
Furthermore, it loses stability via a branch point when = 1 and 3.
Proof. The stability domain of E2 is given in [2] and presented in Fig. 1 ( S
(E2)). By the same procedure as in
Proposition 2.1, we can show that E2 bifurcates to the branches of fixed points E1 and E3(E4) at = 1 and 3,
respectively.
Proposition 2.3. E3 (E4) is asymptotically stable for ( , ) ∈ S
(E3) = S
(E31) ∪ S
(E32) where
S
(E31) =
⎧
⎪⎨
⎪⎩
( , )|3 < < 1 +
√
5, 0 < < 31( ) =
2
1 + 5 − ( − 1)2
⎫
⎪⎬
⎪⎭
,
and
S
(E32) = ( , )| > 1 +
√
5, 0 < < 32( ) =
2
( − 1)2
− 4
.
W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258 251
It loses stability:
(i) via a flip point when = 31( ), 3 < < 1 +
√
5.
(ii) via a Neimark–Sacker bifurcation point when = 32( ), > 1 +
√
5.
Moreover, it loses stability and bifurcates to the branch of E2 fixed points along = 3.
Proof. The stability boundaries of E3 were computed in [2] and are sketched in Fig. 1 ( S
(E3,4)). It is easy to prove
that E3 bifurcates to a branch of fixed points E2 at = 3, by the same procedure as in Proposition 2.1.
Proposition 2.4. The flip bifurcation in Proposition 2.1 is subcritical.
Proof. We show that E1 undergoes a subcritical flip when = 2/(1 + ), 0 < < 1. It is sufficient to show that the
critical normal form coefficient b,
b = 1
6 p, C(q, q, q) + 3B(q, (I − A(1)
)−1
B(q, q)) , (11)
derived by the parameter-dependent center manifold reduction is negative, see [18, Chapter 8 and 16], where A(1) is
the Jacobian of (1) at E1, B(., .) and C(., ., .) are the second and third order multilinear forms, respectively, p and q
are the left and right eigenvectors of A(1), respectively. These vectors are normalized by p, q = 1, q, q = 1, where
., . is the standard scalar product in R2. We obtain
q =
q1
q2
= p =
p1
p2
=
1√
2
− 1√
2
, (12)
[B(q, q)]1 =
n
j,k=1
j2((1 − )x1 + x2(1 − x2))
jxj jxk
qj qk = −2 q2q2 = − , (13)
[B(q, q)]2 =
n
j,k=1
j2((1 − )x2 + x1(1 − x1))
jxj jxk
qj qk = −2 q1q1 = − . (14)
Let = (I − A(1))−1
B(q, q), then we have = −1
−1
and find
[B(q, )]1 = − q2 2 =
√
2
2
− 1
, [B(q, )]2 − q1 1 = −
√
2
2
− 1
. (15)
Since the third order multilinear form C(q, q, q) is zero, the critical normal form coefficient b is given by b= 2/( −1).
It is clear that b < 0, since 0 < < 1 and > 0 in S
(E1).
Proposition 2.5. The flip point in Proposition 2.2 is sub- or supercritical in cases (i) and (ii), respectively.
Proof. First we consider case (i) and show that the flip point is subcritical. It is sufficient to show that b < 0 where b
is defined in (11). The normalized left and right eigenvectors for A(2) are given by
q =
q1
q2
= p =
p1
p2
=
1√
2
− 1√
2
, (16)
where A(2) is the Jacobian of (1) at E2. B(q, q) is given by
[B(q, q)]1 = −2 q2q2 = − , [B(q, q)]2 = −2 q1q1 = − . (17)
252 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
We proceed with the computation of = (I − A(2))−1
B(q, q) and obtain
= 1−
1−
, b =
2
1 −
. (18)
So b < 0, since 1 < < 2 and > 0 in S
(E2).
In case (ii) we obtain b = 2/3( − 1). So b > 0, since 2 < < 3 and > 0 in S
(E2).
Proposition 2.6. The flip bifurcation in Proposition 2.3 is subcritical.
Proof. Similar to the previous cases we show that the critical normal form coefficient b < 0. The Jacobian matrix (1)
at E3 (E4) is
A(3)
=
1 − − (1 +
√
( + 1)( − 3))
− (1 −
√
( + 1)( − 3)) 1 −
, (19)
and has a multiplier −1 when = 2/(1 + 5 − ( − 1)2
), 3 < < 1 +
√
5. The left and right eigenvectors associated
to the eigenvalue −1 are given by
q =
− 4 − 2 + 2
−1 + −3 + 2 − 2
, p =
−1 + −3 + 2 − 2
− 4 − 2 + 2
. (20)
To avoid complicated computations we do not normalize p and q, since rescaling does not change the sign of b provided
p, q > 0 (it can be proved easily that this is the case). B(q, q) is computed as:
B(q, q) =
⎛
⎜
⎜
⎜
⎝
4 (−1 + (−3 + 2 − 2 )2
)
(1 + 4 − 2 + 2 )(−4 + 2 − 2 )
−
4
(1 + 4 − 2 + 2 )
⎞
⎟
⎟
⎟
⎠
. (21)
The vector = (I − A(3))−1
B(q, q), is given by
=
⎛
⎜
⎜
⎜
⎜
⎝
2 (−1 + (−3 + 2 − 2 )2
( + 1)( − 3)(−4 + 2 − 2 )
+
2(1 +
√
( + 1)( − 3))
( + 1)( − 3)
2(−1 +
√
( + 1)( − 3)) (−1 + (−3 + 2 − 2 )2
)
( + 1)( − 3)(−4 + 2 − 2 )
−
2
( + 1)( − 3)
⎞
⎟
⎟
⎟
⎟
⎠
. (22)
B(q, ) can be computed from (20) and (22):
B(q, )=
⎛
⎜
⎜
⎜
⎝
−
8(12−4 2+8 + (−3+ 2−2 ) 2−2 (−3+ 2−2 ) )(−1+ (−3+ 2−2 )) 2
(4− 2+2 )
3
2 ( −3)( +1)(1+ (4− 2+2 ))
−
8(−6−6 (−3+ 2−2 )+2 2−4 + (−3 + 2 − 2 ) 2 − 2 (−3 + 2 − 2 ) ) 2
(−4 + 2 − 2 )( − 3)( + 1)(1 + 4 − 2 + 2 )
⎞
⎟
⎟
⎟
⎠
. (23)
Finally, the normal form coefficient b can be computed:
b = −
576( 2 − 2 − 2 −3 + 2 − 2 − 2) 2
(1 + 4 − 2 + 2 )(−4 + 2 − 2 )2
. (24)
W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258 253
We will prove that the factor h1( ) = 2 − 2 − 2( −3 + 2 − 2 + 1) in (24) is positive when 3 < < 1 +
√
5.
Equivalently we have to prove that ( 2 −2 −2)2
−4( 2 −2 −3) 0. Since dh1/d ( )=4( −1)( 2 −2 −4) < 0,
h1(3) = 1 and h1(1 +
√
5) = 0, we conclude h1 0. So b < 0.
We remark that our numerical evidence indicates that the Neimark–Sacker bifurcation in Proposition 2.3 is supercritical.
This is based on the numerical computation of the normal form coefficient d (see [18, Chapter 8 , 16]). However,
we were not able to prove this analytically.
3. Numerical bifurcation analysis in the economically relevant region
In this section we concentrate on the region 1 which is economically relevant. Since a complete analytical
bifurcation study of the iterates of (1) is not feasible, we perform a numerical bifurcation analysis by using the
MATLAB package CL_MATCONTM, see [14,16]. The bifurcation analysis is based on continuation methods, tracing
out the solution manifolds of fixed points while some of the parameters of the map vary, see [7].
3.1. Numerical bifurcation of E2
By continuation of E2 with = 2.5 and free, we see that E2 loses stability via a supercritical PD point when
crossing the hyperbola = 22( ). A stable 2-cycle is given by C2 ={X2
1, X2
2} where X2
1 =(0.658212, 0.658212), X2
2 =
(0.527341, 0.527341), for = 1.366229. This 2-cycle loses stability at a supercritical PD point (of the second iterate)
for = 1.490763. A stable 4-cycle is given by C4 = {X4
1, X4
2, X4
3, X4
4} where X4
1 = (0.3851221, 0.479532), X4
2 =
(0.745563, 0.649279), X4
3 = (0.479532, 0.38512201), X4
4 = (0.649279, 0.745563).
The multipliers of the fixed point of the fourth iterate in X4
1 are 0.406438 and 0.129274. This 4-cycle with the
parameter values is depicted in Fig. 2. We note that the 4-cycle is invariant under the reflection (x1, x2) → (x2, x1).
This 4-cycle loses stability via a supercritical Neimark–Sacker bifurcation.
Stability regions of the 2-cycles ( S,i
2 , i = 1, 2, 3) and 4-cycles ( S,i
4 , i = 1, 2, 3) are given in Fig. 3. They stretch
into the economically relevant region 1. In this figure the regions S,2
2 and S,2
4 indicate bistability of 2- and 4-cycles
with E3 (E4), respectively. We note that the stability region of the 2-cycle is bounded by the PD2 curve and a curve of
branch points of the second iterate, when 3. This curve of branch points bifurcates from the LPPD point on the PD
0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75
0.4
0.45
0.5
0.55
0.6
0.65
0.7
0.75
x1
x2
X1
X3
X2
X4
Fig. 2. A stable 4-cycle for = 1.509191 and = 2.5.
254 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8
0.8
1
1.2
1.4
1.6
1.8
2
R4
R3
μ
ρ
PD
PD2
NS4
LPPD
R2
PD (ρ31(μ))
NS (ρ32(μ))
Ω4
S,3
Ω2
S,2
Ω4
S,2
Ω4
S,1
Ω2
S,1
Ω2
S,3
Fig. 3. The Neimark–Sacker curve of Run 2, the flip curve of Run 6 of Section 4 and the curve of branch points of the second iterate, the stability
regions of S
2 and S
4 in ( , ) space.
curve of the original map. This curve is shown by ∗ in Fig. 3 and is completely in the economically relevant region.
We note that the LPPD point is on the boundary of the economically relevant region.
3.2. Numerical bifurcation study of E3 (E4)
We now do a continuation of the fixed point E3 starting from = 4, = 0.1 in the stable region bounded by the
curve = 32( ). The parameter is free, we call this Run 1:
label = NSm, x = (0.904508 0.345492 0.400000)
normal form coefficient of NSm = -7.372800e+000
A supercritical Neimark–Sacker bifurcation point is detected for =0.4. Thus, the fixed point E3 (E4) is transformed
from stable to unstable through an NS point at which a closed invariant curve is created around the unstable fixed point
E3 (E4). We now compute the Neimark–Sacker curve, by starting from the NS point in Run 1, with free parameters
and , this is Run 2:
label = R4 , x = ( 0.849938 0.439960 1.000000 3.449490 0.0000 )
normal form coefficient of R4 : A = -3.000000e+000 + -2.019371e-017i
label = R3 , x = ( 0.825542 0.476627 1.500000 3.309401 -0.500)
normal form coefficient of R3 : Re(c_1) = -1.333333e+000
label = R2 , x = ( 0.809017 0.500000 2.000000 3.236068 -1.00 )
normal form coefficient of R2 : [c , d] = 1.340433e+003, -3.351046e+003
A picture of the Neimark–Sacker curve of Run 2 is given in Fig. 3.
Since the R2 and R3 points are not in the region 1 we postpone their study to Section 4.
We now consider the R4 point in Run 2. Since |A| > 1, two cycles of period 4 of the map are born. A stable
4-cycle for = 0.990844 and = 3.466353 is given by C4 = {X1, X2, X3, X4} where X1 = (0.841774, 0.407047),
X2 = (0.836685, 0.461186), X3 = (0.861140, 0.473539), X4 = (0.864133, 0.4150395). We present this cycle in
Fig. 4. We note that it is not invariant under the reflection (x1, x2) → (x2, x1). The multipliers in X1 are 1 =0.901140
and 2 = 0.675526, confirming the stability of the 4-cycle.
To determine the stability domain of the 4-cycle we compute in Run 3 two branches of fold curves of the fourth
iterate, emanating from the R4 point, by switching at the R4 point. These fold curves exist because |A| > 1, where A
W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258 255
0.83 0.835 0.84 0.845 0.85 0.855 0.86 0.865
0.41
0.42
0.43
0.44
0.45
0.46
0.47
0.48
0.49
x1
x2
x1
x2
x3
x4
Fig. 4. A stable 4-cycle for = 0.9999617 and = 3.449802.
is the normal form coefficient of the R4 point. The stable fixed points of the fourth iterate exist in the wedge between
the two fold curves. We note that there is no bistability with fixed points of the original map.
label = CP, x = (0.849982 0.439945 0.999745 3.449889)
normal form coefficient of CP s = 4.009280e+002
label = LPPD , x = (0.841586 0.354516 0.935299 3.566686)
normal form coefficient of LPPD :[a/e , be] =
2.574002e+000, -5.829597e+001
label = CP, x = ( 0.849982 0.439945 0.999745 3.449889 )
normal form coefficient of CP s = 4.009280e+002
label = LPPD , x = ( 0.836428 0.522216 1.071080 3.486079 )
normal form coefficient of LPPD :[a/e , be] =
1.733856e+000, -2.471512e+001
Two cusps, CP, and two LPPD bifurcation points are detected on the fold branches of the fourth iterate. The CP
points are merely the R4 point from which we started. We can further compute the stability boundaries of the 4-cycle
by computing the flip curve of the fourth iterate rooted at the detected LPPD points. The stable region S
4 of C4 is
bounded by two fold curves and a flip curve of the fourth iterate, see Fig. 5. Furthermore, if we continue the fixed point
of the fourth iterate starting from X1, it loses stability via a supercritical PD point where = 3.545530. It means that
a stable 8-cycle is born when > 3.545530.
We note that we have bistability of three different 4-cycles in a region bounded by the curves of the PD of the second
iterate, a fold and the PD curve of the fourth iterate. This region is shown as S
4,4 in Fig. 6. Furthermore, we have a
small bistability region of two 4-cycles and a 2-cycle. This bistability region is shown as S
2,4 in Fig. 6.
4. Bifurcations of E3 (E4) in the region > 1
Now we consider the R2 point computed in Run 3 of Section 3.2. Since the first component of the normal form
coefficient c =1.340433e +003 is positive, the bifurcation scenario near the R2 point is analogous to [18, Fig. 9] (case
s = 1). For the parameter values in the wedge between the PD ( 31) and NS ( 32) curves, the map has an unstable
2-cycle that coexists with a stable fixed point.
256 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
3.4 3.45 3.5 3.55 3.6
0.9
0.95
1
1.05
1.1
1.15
μ
ρ
PD4
NS
LPPD
LP4
Ω4
SR4
Fig. 5. Two fold bifurcation curves of the fourth iterate emanate from the R4 point in ( , ) space.
2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8
0.8
0.9
1
1.1
1.2
1.3
μ
ρ
R4
LPPD
PD (ρ31)
PD2
NS4
NS (ρ32)PD (ρ22)
LP4
PD4
LPPD
ΩS,2
4
ΩS,2
2
ΩS
4,4
ΩS
2,4
Fig. 6. Bistability regions of 4-cycles, 2-cycles and fixed points.
Next we consider the resonance 1:3 point in Run 3 of Section 3.2. Since its normal form coefficient is negative, the
bifurcation picture near the R3 point is qualitatively the same as presented in [18, Fig. 9]. In particular, there is a region
near the R3 point where a stable invariant closed curve coexists with an unstable fixed point. For parameter values
close to the R3 point, the map has a saddle cycle of period three.
Furthermore, a Neutral Saddle bifurcation curve of fixed points of the third iterate emanates. We compute this curve
by branch switching at the R3 point of Run 5. This curve is presented in Fig. 7. Further, a stable 3-cycle exists not far
from the R3 point (this is not guaranteed by the theory but it is found in many examples, e.g [15]). The stability region
of this cycle is bounded by fold and NS bifurcation curves of the third iterate of the map ( S
3). These boundary curves
W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258 257
3.28 3.3 3.32 3.34 3.36 3.38
1.25
1.3
1.35
1.4
1.45
1.5
1.55
1.6
μ
ρ
R3
NS3
LP3
NS
ΩS
3
Fig. 7. Two stability boundary curves (LP and NS) for the stable 3-cycle, together with the NS curve of Run 3 and the curve of Neutral Saddle
bifurcation points of the third iterate.
are given in Fig. 7. We have bistability of the fixed point E3 (E4) with the fixed point of the third iterate of the map in
the region that is bounded by the fold and NS curves and the hyperbola = 32( ).
5. Conclusions
We studied the dynamical behavior of the Kopel model and computed the stability boundaries of 2-, 3- and 4-cycles.
We showed analytically that the trivial fixed point E1 undergoes a subcritical flip bifurcation when = 1( ). The
nontrivial symmetric fixed point E2 loses its stability via a subcritical flip point when crossing the curve = 21( ),
whereas it loses stability along the curve = 22( ) via a supercritical flip point.
For E3 (E4) the transition of stability to unstability is possible via a subcritical flip point when crossing the curve
= 31( ) and via a supercritical Neimark–Sacker point when crossing the curve = 32( ). E1 branches to E2 along
= 1 and E2 branches to E3 (E4) when = 3.
In the case of E2, the stability domain of the 2-cycle is bounded by two flip curves of the first and second iterates
when 1 < 3 and by flip and branch point curves of the second iterate when 3. These stability domains are shown
as S,i
2 , i = 1, 2, 3 in Fig. 3. We have bistability of E3 (E4) and a 2-cycle in S,2
2 . The stability domain of the 4-cycle
is bounded by curves of flip and Neimark–Sacker bifurcations of the second and fourth iterates, respectively. These
stability domains are shown as S,i
4 , i = 1, 2, 3 in Fig. 3. Moreover, we have bistability of E3 (E4) and a 4-cycle in
S,2
4 .
In the case of E3 (E4), the stability domain of the 4-cycle is bounded by two branches of fold curves and a flip curve
of the fourth iterate, and shown as S
4 in Fig. 5.
The stability domain of the 3-cycle is bounded by curves of fold and Neimark–Sacker bifurcations of the third iterate;
we have bistability of the fixed point of the map and the 3-cycle.
We have bistability of three different 4-cycles in S
4,4 in Fig. 6 and bistability of two 4-cycles and a 2-cycle in S
2,4
in Fig. 6.
We note that for ∈ [0, 1], the stability regions of cycles are economically interesting. This applies to the stability
of 2- and 4-cycles in Fig. 3, stability of the 4-cycle in Fig. 5, bistability of two 4-cycles as well as bistability of 2- and
4-cycles in Fig. 6.
258 W. Govaerts, R. Khoshsiar Ghaziani / Journal of Computational and Applied Mathematics 218 (2008) 247–258
Acknowledgment
The authors thank the referee for several suggestions which have substantially improved the paper, in particular with
respect to the economic background of the paper and the relevance of the obtained results.
References
[1] H.N. Agiza, Explicit stability zones for Cournot games with 3 and 4 competitors, Chaos Solitons Fractals 9 (1998) 1955–1966.
[2] H.N. Agiza, On the analysis of stability, bifurcation, chaos and chaos control of Kopel map, Chaos Solitons Fractals 10 (1999) 1909–1916.
[3] H.N. Agiza, G.I. Bischi, M. Kopel, Multistability in a dynamic Cournot game with three oligopolists, Math. Comput. Simulation 51 (1999)
63–90.
[4] H.N. Agiza, A.A. Elsadany, M. Kopel, Nonlinear dynamics in the Cournot duopoly game with heterogeneous players, Physica A 320 (2003)
512–524.
[5] H.N. Agiza, A.S. Hegazi, A.A. Elsadany, Complex dynamics and synchronization of a duopoly game with bounded rationality, Math. Comput.
Simulation 58 (2002) 133–146.
[6] E. Ahmed, H.N. Agiza, Dynamics of a Cournot game with n-competitors, Chaos Solitons Fractals 9 (1998) 1513–1517.
[7] E.L. Allgower, K. Georg, Numerical Continuation Methods: An Introduction, Springer, Berlin, 1990.
[8] D.R. Anderson, N.G. Myran, D.L. White, Basin of attraction in a Cournot duopoly model of Kopel, J. Differential Equations Appl. 11 (10)
(2005) 879–887.
[9] G.I. Bischi, A. Mammana, L. Gardini, Multistability and cyclic attractors in duopoly games, Chaos Solitons Fractals 11 (2000) 543–564.
[10] J.P. Carcassès, Determination of different configurations of fold and flip bifurcation curves of one or two-dimensional map, Internat. J. Bifur
Chaos 3 (4) (1993) 869–902.
[11] J.P. Carcassès, Singularities of the parametric plane of an n-dimensional map, Determination of different configurations of fold and flip
bifurcation curves, Int. J. Bifurcation Chaos 5 (2) (1995) 419–447.
[12] J.P. Carcassès, C. Mira, C. Simó, M. Bosch, J.C. Tatjer, “Cross-road area–spring area” transition (I) Parameter plane representation, Internat.
J. Bifurcation Chaos 1 (2) (1991) 339–348.
[13] A. Cournot, Researches into the principles of wealth, 1963 (English Translation), Irwin Paperback Classics in Economics (original: Recherches
sur les principes mathématiques de la théorie des richesses, 1838).
[14] A. Dhooge, W. Govaerts, Yu.A. Kuznetsov, W. Mestrom, A. Riet, CL_MATCONT: a continuation toolbox in Matlab, 2004, http://
www.matcont.UGent.be .
[15] W. Govaerts, R. Khoshsiar Ghaziani, Numerical bifurcation analysis of a nonlinear stage structured cannibalism population model, J. Differential
Equations Appl 12 (10) (2006) 1069–1085.
[16] W. Govaerts, R. Khoshsiar Ghaziani, Yu. A. Kuznetsov, H. Meijer, Bifurcation analysis of periodic orbits of maps in Matlab, 2006, preprint.
[17] M. Kopel, Simple and complex adjustment dynamics in Cournot duopoly models, Chaos Solitons Fractals 12 (1996) 2031–2048.
[18] Yu.A. Kuznetsov, Elements of Applied Bifurcation Theory, third ed., Springer, New York, 2004.
[19] Yu.A. Kuznetsov, H.G.E. Meijer, Numerical normal forms for codim 2 bifurcations of maps with at most two critical eigenvalues, SIAM J. Sci.
Comput. 26 (2005) 1932–1954.
[20] C. Mira, J.P. Carcassès, M. Bosch, C. Simó, J.C. Tatjer, “Cross-road area–spring area” transition (II) foliated parametric representation, Internat.
J. Bifurcation Chaos 1 (2) (1991) 339–348.
[21] J. Nash, Non-cooperative games, Ann. Math. 54 (1951) 286–295.
[22] T. Puu, Chaos in duopoly pricing, Chaos Solitons Fractals 1 (1991) 573–581.
[23] T. Puu, The Chaotic duopolists revisited, J. Econom. Behav. Organ. 33 (1998) 385–394.
[24] D. Rand, Exotic phenomena in games duopoly models, J. Math. Econom 5 (1978) 173–184.
[25] J.B. Rosser, The development of complex oligopoly dynamics theory, in: T. Puu, I. Sushko (Eds.), Text Book Oligopoly Dynamics: Models
and Tools, Springer, Berlin, 2002, pp. 31–83.