# A Dynamic Duopoly Cournot Model with R&D Efforts and Its Dynamic Behavior Analysis

Authors: Wei Zhou; Jie Zhou; Tong Chu; Hui Li
Journal: Complexity
 (2020)

---

## Abstract

In this paper, a dynamic two-stage Cournot duopoly game with R&D efforts is built. Then, the local stability of the equilibrium points are discussed, and the stability condition of the Nash equilibrium point is also deduced through Jury criterion. The complex dynamical behaviors of the built model are investigated by numerical simulations. We found that the unique route to chaos is flip bifurcation, and the increase of adjusting speed will cause the system to lose stability and produce more complex dynamic behavior. In addition, we also found the phenomenon of multistability in the given model. Several kinds of coexistence of attractors are shown. In particular, we found that boundary attractors can coexist with internal attractors, which also aggravates the complexity of the system. At last, the chaotic state in the built system has been successfully controlled.


---

## Body

## 1. Introduction

In recent years, the market competition mode among firms are no longer subjected to price or quantity, but more towards soft power, including technical level, talent resources, and corporate culture. With more and more systematic and mature development of the organization management for a great quantity of large firms, R&D (the abbreviation of Research and Development) is becoming more and more important for them. The visionary entrepreneurs will spend a lot of money to build their own exclusive R&D team. It may be risky for the capitalists who pursue maximum profits to do such things which may be harmful to their revenues. However, the reasons for these firms that lay much stress on R&D activity can be listed as follows: on the one hand, R&D activities can improve production efficiency, reduce production costs, and then the entrepreneurs can gain surplus value proliferation; on the other hand, R&D activities can also improve quality of products so as to increase advantage competitiveness of firms in the market. Hence, R&D activity gradually attracts a lot of attentions of economists.In the process of carrying out R&D activities, due to the flow of technical talents, information diffusion, and other factors, R&D spillover effect may emerge among firms [1, 2]. Numerous scholars studied and advanced this problem. Zhou et al. [3] studied the types of bifurcation, intermittent chaos, and coexistences in a dynamical two-stage Cournot duopoly game with R&D spillover effect. Zhou and Wang [4] studied the stability and multistability of a duopoly game model with R&D spillover. While in the case that the firms take highly confidential measures for its R&D achievements, the positive effects of R&D achievements only promote the firms that carry out R&D activities, but not the competitors. In other words, there is no spillover effect. Based on this assumption, D’Aspremont and Jacquemin [5] built the foundation of the performance of R&D cooperation in the product differentiation market. Luckraz [6] established a Cournot duopoly game model with isoelastic demand function and also considered the case without R&D spillovers. The work of Chatterjee [7] showed that the firms in the market undertake R&D activity and studied that R&D incentive may depend on the information structures they had. In addition, Ferreira et al. [8] constructed myopic optimal discrete and continuous R&D dynamics.The most extensive research on spillover effect is technology spillover. AJ Model [5] and KMZ model [9] are the most representative ones. In both cases, they are both two-stage game models. Also, the two-stage game presents good adaptability in dealing with problems about uncertainty of R&D market and mutuality of competitive strategies. That is to say, on the first stage, both these two firms conduct R&D activities to reduce their production costs. And on the second stage, they compete in market in term of the price or the quantity. The investigations on the two-stage model have gradually become the focus of economists. While at the beginning of the study, many scholars were limited in studying the properties of the static model [10–12]. All the investigations about the static model just illustrated the equilibrium solutions under different R&D strategies. No one can deny that the static model does have some limitations. One limitation is that it can only analyze the single supply and demand balance between firms. Once the supply and demand factors change, the corresponding supply relationship will shift.In this research, we introduce nonlinear dynamics theory to discuss evolution process among firms’ games. Many scholars had carried out some studies about this issue. Bischi and Lamantia [13, 14] revealed the local and global dynamical properties of the two-stage oligopoly game model with adaptive dynamic mechanism and described its complex evolutionary behaviors. Askar et al. [15] investigated the properties in a dynamic Cournot duopoly game model with nonlinear demand function. A homogeneous Cournot duopoly market possessing two different mechanisms was proposed by Elettreby [16]. Most of the similar works can be seen in [17–19]. Additionally, the dynamical two-stage game has raised much interest in the economic dynamical system. For instance, Zhang et al. [20] proposed a two-stage Cournot duopoly game of semicollusion in production and analyzed its complex dynamics, such as local stabilities and existence. Ma et al. [21, 22] have focused the analysis of complex dynamical behaviors of the two-stage dynamic games.Furthermore, although in the electronic information world, the corresponding channels for firms to obtain market information have proliferated. However, it cannot be ignored that if the firms want to get complete information, they still need to pay its costs. That is to say, few firms are willing to pay high cost to obtain all the information about market demand and consumer preferences. So, most firms have incomplete information about the market, so they could be regarded to be bounded rational [23] when making decisions. In fact, compared with those having complete information of demand function, it is much easier to obtain the existing information, such as local estimation of marginal profit. That is, gradient adjustment mechanism which has attracted many scholars to carry out deep research studies on this topic (see e.g., [23–26]). Cavalli and Naimzada [27] took the factors such as rationality and information set into consideration and adopted a traditional dynamical adjustment mechanism based on classical best reply function. Elsadany [28] interpreted a Cournot duopoly dynamics model with different strategies including bounded rationality. Fanti and Gori [29] put forward a Bertrand duopoly game model with heterogeneous product expectations and further studied the properties of the model. And Naimzada and Sbragia [30] dealt with an oligopoly game model with nonlinear demand function and bounded rationality.Additionally, other topics worth studying in nonlinear dynamics are global dynamics, synchronization, and multistability. In recent years, great deals of scholars have deepened the study of these topics. Askar [31] dealt with the dynamical characteristics with log-concave demand of equilibrium points, such as stability, bifurcation, and chaos. Andaluz and Jarne [32] analyzed the dynamics of not only quantity-setting competition but also price-setting competition. Agliari et al. [19] observed that some bifurcations can lead to changes in the structure of attractors and their basins. And a higher degree of production differentiation tends to reduce competition. Tu and Wang [33] built a dynamical R&D two-stage input competition triopoly game model with bounded rationality and further analyzed the dynamics properties of this model. Bischi et al. [34, 35] had studied the phenomena of dynamic properties, such as synchronization, intermittency, coexisting attractors, and complex basins. Fanti et al. [36] provided local and global dynamic characteristics in a duopoly game model with price-setting competition and market share delegation.Basing on the above analyses, the behaviors including dynamic evolution, intermittency, and multistability are exploited to interpret the complexity of competition. We take use of a Cournot duopoly model with differentiated products to consider the impact of adjustment speeds. And the influence of firm’s adjustment speed on the development of firms is studied through combining a two-stage dynamic game with the method of numerical simulation and the properties of nonlinear dynamics. It is worth to focus on whether the firm will gain more advantages at a lower adjustment speed.Furthermore, another noticed feature is chaos control. Since chaos could increase the disorder and unpredictability of the established model. In the sense of economics, if actual economic market is in a chaotic state, then the resources of firms will not be controlled by market supply, demand relationship, and the government macrocontrol to some extent. This may disrupt the development of market economy. However, as for firms in the market, they prefer to pursue a long-term and stable development, which can raise their accumulation of capitals to expand themselves industrial scale. Hence, it is necessary to control chaos. A large number of scholars have carried out relevant research studies and achieved some results. For instance, feedback control was given by Elsadany et al. [37, 38]. Ma and Ji [39] studied the dynamical repeated Cournot model by using the chaos control method with feedback control. Parameters adjustment was provided to control chaos by Askar and Al-khedhairi [40].The main purpose of this research is to investigate the influence of continuous change of system parameters on the eventual dynamical scenarios of the dynamic duopoly model with R&D efforts. The local stability of boundary equilibrium and Nash equilibrium is provided by Jury stability criterion and numerical simulation to obtain the internal complexity of the model. We show that the successive change of system parameters not only causes the system to lose its stability but also causes the high complexity of eventual behaviors of the system, such as the number of coexisting attractors, the critical bifurcation of attractors, and the global bifurcation of attracting domain.This research is organized as follows. In Section2, a two-stage dynamical game model without R&D efforts is established, where these two firms’ products are of the same kind but different products. In addition, its stability and the conditions of stability of equilibrium points are discussed. In Section 3, the complex evolution law based on speed of adjustment is discussed in details. Through numerical simulation, the complexity of bifurcations, evolution of strange attractors, intermittent chaos, multistability, and control of chaos are studied. At last, the impacts of adjustment speeds on built model are concluded in Section 4.


## 2. The model

Assume that there are two technology-based firms operating in the market, labeled byi (i = 1, 2), and they supply homogeneous products to the same market. Research and Development (R&D) is a kind of creative behavior for the firms to increase production efficiency and decrease costs so as to gain more market share. While in the actual economics environment, R&D activity is always carried out in a stepwise manner: research phase and development phase. The process of competition between two R&D firms can be simulated by a two-stage game. On the first stage, these two firms conduct R&D activities to reduce their production costs. And on the second stage, they compete in the market in terms of the price or the quantity. The difference between the R&D levels of the firms can result the difference of goods in quality, function, and externalization. Hence, the inverse demand function can be written as(1)piqi,q3−i=a−qi−bq3−1,i=1,2,where pi is the price of firm i (i = 1, 2), qi and q3−i represents the outputs of firm i and its rival, respectively, a>0 denotes the maximum scale of this market, b is defined as the differentiating degree of the products of two firms, and b∈0,1. Moreover, if b=0, it means that the products of firm 1 and firm 2 are mutually independent, and the oligopoly market will be degenerated into two completely independent monopoly markets. If 0<b<1, it indicates that the products of these two firms can replace each other. And if b=1, it means that there is no difference between the products of these two firms. In other words, the products of firm 1 and firm 2 are totally identical.Consider now the R&D efforts of firmi (i = 1, 2) labeled by xi. Then, the effective marginal cost of firm i can be represented as(2)Cixi=c−xi,i=1,2,where c represents the marginal cost of two firms without presence of R&D activities and a>c>0. Equation (2) states that technology is totally confidential and there is no R&D spillover. The investments of firm i (i = 1, 2) can be seen as a quadratic function with respect to its R&D efforts, which can be expressed as Ixi=γxi2/2 (i = 1, 2), where γ is a parameter about the technical innovation cost of firm i. The innovation ability of firm i is stronger when the technical innovation cost parameter γ is tinny. Therefore, by following the above assumptions, the profits of firm 1 and firm 2 can be represented as follows:(3)π1q1,q2,x1,x2=p1q1,q2−C1x1,x2q1−γx122,π2q1,q2,x1,x2=p2q1,q2−C2x1,x2q2−γx222.Substituting equations (1) and (2) into equation (3), we can get the specific expression of profit function, which could be given as(4)π1q1,q2,x1,x2=a−q1−bq2−c−x1q1−γx122,π2q1,q2,x1,x2=a−q2−bq1−c−x2q2−γx222.The marginal profits of these two firms captured by taking the derivative of (4) with respect to qi (i = 1, 2) are given as(5)∂π1∂q1=a−2q1−bq2−c+x1,∂π2∂q2=a−2q2−bq1−c+x2.Letting∂πi/∂qi=0 (i = 1, 2), the subgame perfect Nash equilibrium (SPNE) of firm i (i = 1, 2) in the second stage can be obtained as(6)q1∗=a−c2−b+2x1−bx24−b2,q2∗=a−c2−b+2x2−bx14−b2.The profit function about R&D effortsxi captured by taking equation (6) into equation (4) in reverse order can be obtained as(7)π1x1,x2=a−c2−b+2x1−bx224−b22−γx122,π2x1,x2=a−c2−b+2x2−bx124−b22−γx222.Computing the derivative of equation (7) with respect to xi (i = 1, 2), the local marginal profits about the R&D efforts of two firms are given as(8)∂π1∂x1=4a−c2−b+2x1−bx24−b22−γx1,∂π1∂π1=4a−c2−b+2x2−bx14−b22−γx2.Actually, it is impossible to make perfect rational decision with incomplete market information, such as incomplete demand information, the investments of their competitors, and the level of technical innovation. These incomplete elements make the firms “myopic” in some aspects. These “myopic” firms can only adjust their strategies through continuous “trial and error.” That is, the output at the next period is decided according to their own experience, so these firms can be seen as bounded rational. The firm shall have a decision on its R&D effortxi at t + 1 period by using a local estimate of marginal profit ∂πi/∂xi at t period. More in detail, if ∂πi/∂xi>0, the firm i will increase its R&D investment in the next period. And the firm will keep invariant outputs in next period if the local estimate of marginal profit meets ∂πi/∂xi=0. Then, if ∂πi/∂xi<0, the R&D investment of the firm i in the next period will be reduced. On the basis of the above thoughts, a dynamical gradient adjustment mechanism is introduced [20]; then, the dynamic model can be constructed as(9)x1t+1=x1t+α1x1t∂π1∂x1,x2t+1=x2t+α2x2t∂π2∂x2.Taking equation (8) into equation (9), the two-dimensional nonlinear discrete-time evolution of dynamic system (9) is as follows:(10)x1t+1=x1t+α1x1t4a−c2−b+2x1t−bx2t4−b22−γx1t,x2t+1=x2t+α2x2t4a−c2−b+2x2t−bx1t4−b22−γx2t,where αi>0 denotes the adjustment speed, which determines the variation ratio of R&D efforts xi of firm i in the next period.


### 2.1. Stability Analysis of the Equilibrium Points

System (10) is a two-dimensional noninvertible map whose iteration decides the trajectory of R&D effort. In order to discuss local qualitative behaviors of the solutions of system (10), we need to find equilibrium points of the game between these two firms. Let xit+1=xiti=1,2; then, we can get the nonlinear equations as follows:(11)α1x1t4a−c2−b+2x1t−bx2t4−b22−γx1t=0,α2x2t4a−c2−b+2x2t−bx1t4−b22−γx2t=0.By solving these nonlinear equations, four equilibrium points of system (10) can be obtained. They are(12)E10,0,E20,−42−ba−c8−γ4−b22,E3−42−ba−c8−γ4−b22,0,E4x∗,x∗,where x∗=−42−ba−c/8−4b−γ4−b22. The fixed points E1, E2, and E3 are boundary equilibrium points. And the fixed point E4 is an unique Nash equilibrium point. In the sense of economy, the inequalities a>c and 0≤b≤1 should be satisfied. And if the other parameters satisfy the condition(13)γ>84−b22,all the four equilibrium points of system (10) will be nonnegative. The corresponding Jacobian matrix of system (10) at any point x1,x2 has the form(14)J=1+α142−ba−c+28−γ4−b22x1−4bx24−b22−4α1b4−b22x1−4α2b4−b22x21+α242−ba−c+28−γ4−b22x2−4bx14−b22.In addition, the local stability of equilibrium points of system (10) depends on the eigenvalues of Jacobian matrix at the given equilibrium point.Proposition 1.
The boundary equilibrium pointE1 is an unstable node.Proof.
The Jacobian matrix atE1=0,0 is given as follows:(15)JE1=1+4α12−ba−c4−b22001+4α22−ba−c4−b22.
Obviously, the Jacobian matrixJE1 is a diagonal matrix; thus, its corresponding distinct eigenvalues are(16)λ1=1+4α12−ba−c4−b22,λ2=1+4α22−ba−c4−b22.
With the constrictsa>c and 0≤b≤1, λ1,2>1 is satisfied, then we can get that the equilibrium point E10,0 is an unstable node.Proposition 2.
The types of local stability of the equilibrium pointE2 can be classified as follows:(1)
If the condition8+4b−γ4−b22<0 meets, the equilibrium point E2 is an unstable node(2)
If the conditions8+4b−γ4−b22>0, 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, and α2>4−b22/22−ba−c meet, then the equilibrium point E2 is an unstable saddle(3)
If the conditions8+4b−γ4−b22>0, 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, and α2<4−b22/22−ba−c meet, then the equilibrium point E2 is a stable node(4)
Else if8+4b−γ4−b22>0 and α1>4−b228−γ4−b22/22−ba−c8+4b−γ4−b22 are satisfied, then the equilibrium point E2 is an unstable pointProof.
The Jacobian matrix atE2 is(17)JE2=1+4α12−ba−c8+4b−γ4−b224−b228−γ4−b22016α2b2−ba−c4−b228−γ4−b221−4α22−ba−c4−b22,where JE2 is a lower triangular matrix; thus, its eigenvalues are main diagonal elements, that is,(18)λ1=1+4α12−ba−c8+4b−γ4−b224−b228−γ4−b22,λ2=1−4α22−ba−c4−b22.
The type of local stability of this boundary equilibrium pointE2 is determined by different values of the parameters:(1)
If8+4b−γ4−b22<0, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22>1; thus, the equilibrium point E2 is an unstable node.(2)
If8+4b−γ4−b22>0 and 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22 are satisfied, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22<1 meets. And if α2>4−b22/22−ba−c, then λ2=1+4α22−ba−c/4−b22>1 meets; hence, the equilibrium point E2 is an unstable saddle.(3)
If8+4b−γ4−b22>0 and 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22<1. And if α2<4−b22/22−ba−c, then λ2=1−4α22−ba−c/4−b22<1, so the equilibrium point E2 is a stable node.(4)
If the parameters meet with8+4b−γ4−b22>0 and α1>4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22>1; thus, the equilibrium point E2 is an unstable point.Proposition 3.
The types of equilibrium pointE3 are similar to the boundary equilibrium point E2.Proof.
The Jacobian matrix atE3 can be written as follows:(19)JE3=1−4α12−ba−c4−b2216α1b2−ba−c4−b228−γ4−b2201+4α22−ba−c8+4b−γ4−b224−b228−γ4−b22.
Matrix (19) is an upper triangular matrix and the corresponding eigenvalues are(20)λ1=1−4α12−ba−c4−b22,λ2=1+4α22−ba−c8+4b−γ4−b224−b228−γ4−b22.
In addition, the boundary equilibrium pointsE2 and E3 are symmetrically located around the diagonal line in the plane R2. And the local stability conditions of E3 are symmetrical to the conditions of E2 listed in Proposition 2, so the conclusions of E3 are no longer repeated here.
In sense of economics, the boundary equilibrium points interpret a situation that one of these two firms has withdrawn from the market. The equilibrium pointsE2 and E3 can be regarded as the situation that one of these two firms takes the lead in the market. That is to say, the oligopoly market becomes the monopoly market, and the local stability of equilibrium points implies the stability of economic market in a short term. Actually, neither of the above two market situations is expected, only when these two companies restrict each other can they be conducive to the stable development of the market and the country. That is, the game between firms reaches an equilibrium, named “Nash equilibrium.” Nash equilibrium [41] is on the premise of maximizing its own profits. Meanwhile, it also ensures the stable development of the market.
The Jacobian matrix at the Nash equilibriumE4 of system (10) is formulated as(21)JE4=1+4α12−ba−c+28−2b−γ4−b22x∗4−b224α1bx∗4−b224α2bx∗4−b221+4α22−ba−c+28−4b−γ4−b22x∗4−b22,where x∗=−42−ba−c/8−4b−γ4−b22. To simplify the process of analyses, the auxiliary variables A=42−ba−c/4−b22>0, B=8−γ4−b22/4−b22<0, and C=4b/4−b22>0 are introduced. Thus, the matrix JE4 can be rewritten as follows:(22)JE4=1+α1A+2B−Cx∗α1x∗Cα2x∗C1+α2A+2B−Cx∗.
The characteristic polynomial of matrix (22) is(23)Pλ=λ2−TrJ+DetJ=0,where TrJ and DetJ denote the trace and the determinant of matrix (22) at Nash equilibrium point E4. And the specific formulation of TrJ and DetJ are given as(24)TrJ=2+α1+α2A+2B−Cx∗,DetJ=1+α1A+2B−Cx∗1+α2A+2B−Cx∗−α1α2C2x∗2=1+α1+α2A+2B−Cx∗+α1α2A2+2A2B−Cx∗+4BB−Cx∗2.Proposition 4.
Neimark–Sacker bifurcation will not occur at the Nash equilibrium pointE4.Proof.
The corresponding discriminant of roots of characteristic polynomial (23) is given as(25)Δ=α1−α22A+2B−Cx∗2+4α1α2C2x∗2>0,which means that the dynamical system (10) will not generate complex eigenvalues at the Nash equilibrium point E4; hence, the Neimark–Sacker bifurcation will not occur at the Nash equilibrium point E4.
On the basis of canonical stability analysis, the sufficient condition of local asymptotical stability of Nash equilibrium pointE4 is that the eigenvalues of corresponding Jacobian matrix of E4 is within the unit cycle. According to Jury criterion [15], the condition of local asymptotical stability of Nash equilibrium point E4 can be expressed in detail as(i)
1+TrJ+DetJ=4+α1+α2A+2B−Cx∗+α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0(ii)
1−TrJ+DetJ=α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0(iii)
1−DetJ=−α1+α2A+2B−Cx∗−α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0Proposition 5.
WhenB+C<0 is satisfied, the second condition of Jury criterion holds. That is,(26)1−TrJ+DetJ=α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0.Proof.
First, lety=A2+2B2B−Cx+4BB−Cx2, and then its roots can be computed as x1=−A/2B−C and x2=−A/2B. We recall that A=42−ba−c/4−b22>0, B=8−γ4−b22/4−b22<0, and C=4b/4−b22>0; then, we can obtain(27)B−C=8−4b−γ4−b224−b22<0.
The numerical relationship betweenx1 and x2 is that x1=−A/2B−C<x2=−A/2B. If the condition A2+2B2B−Cx∗+4BB−Cx∗2>0 holds, then x∗satisfies the condition x∗<x1 or x∗>x2. Following the constriction x∗=−42−ba−c/8−γ4−b22=−A/2B−C, x∗=−A/B−C>x1=−A/2B−C can be captured. Hence, x∗<x1 does not hold. Only when the condition x∗>x2 meets, the second condition of Jury criterion holds. The corner case is(28)x∗−x2=−AB−C+A2B=−AB+C2BB−C.
Hence, whenB+C<0 and x∗>x2 meet, the second condition of Jury criterion holds.Proposition 6.
If the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold.Proof.
Letα1=0 of first condition of Jury criterion, and we can get 4+α1+α2A+2B−Cx∗>0. As for A>0, B<0, C>0, B−C<0, and x∗=−A/B−C, we have(29)A+2B−Cx∗=A+−A2B−CB−C=−ABB−C<0.
After this, we can deduce thatα2<−2/A+2B−Cx∗.
LetΔ=A2+2B2B−Cx∗+4BB−Cx∗2>0, if 0<α2≤−A+2B−Cx∗/Δ, the third condition of Jury criterion will hold. That is,(30)4+α1+α2A+2B−Cx∗+α1α2Δ=α2A+2B−Cx∗+α1A+2B−Cx∗+α2Δ<0,
In the same time, it can be easily proved that(31)A+2B−Cx∗Δα22+4Δα2+4A+2B−Cx∗<0.which is satisfied for any value of α2. Following the above formulations, we can get −2/A+2B−Cx∗<−2A+2B−Cx∗/Δ. The inequalities −A+2B−Cx∗/Δ<α2<−2/A+2B−Cx∗ and(32)α1<−4−2α2A+2B−Cx∗2A+2B−Cx∗+α2Δ<−α2A+2B−Cx∗2A+2B−Cx∗+α2Δ,can be captured through jury criterion and inequality (31). Then, we can deduce the conclusion that if the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold. In summary, system (10) has an unique Nash equilibrium point E4. When B+C<0 and α1<−4−2α2A+2B−Cx∗/2A+2B−Cx∗+α2Δ, system (10) is local asymptotical stable at Nash equilibrium point E4. Under the above constrictions of the system parameters, these two firms will reach equilibrium state and get a maximum value for total profit after several games. Then, the market will keep stable development.


## 2.1. Stability Analysis of the Equilibrium Points

System (10) is a two-dimensional noninvertible map whose iteration decides the trajectory of R&D effort. In order to discuss local qualitative behaviors of the solutions of system (10), we need to find equilibrium points of the game between these two firms. Let xit+1=xiti=1,2; then, we can get the nonlinear equations as follows:(11)α1x1t4a−c2−b+2x1t−bx2t4−b22−γx1t=0,α2x2t4a−c2−b+2x2t−bx1t4−b22−γx2t=0.By solving these nonlinear equations, four equilibrium points of system (10) can be obtained. They are(12)E10,0,E20,−42−ba−c8−γ4−b22,E3−42−ba−c8−γ4−b22,0,E4x∗,x∗,where x∗=−42−ba−c/8−4b−γ4−b22. The fixed points E1, E2, and E3 are boundary equilibrium points. And the fixed point E4 is an unique Nash equilibrium point. In the sense of economy, the inequalities a>c and 0≤b≤1 should be satisfied. And if the other parameters satisfy the condition(13)γ>84−b22,all the four equilibrium points of system (10) will be nonnegative. The corresponding Jacobian matrix of system (10) at any point x1,x2 has the form(14)J=1+α142−ba−c+28−γ4−b22x1−4bx24−b22−4α1b4−b22x1−4α2b4−b22x21+α242−ba−c+28−γ4−b22x2−4bx14−b22.In addition, the local stability of equilibrium points of system (10) depends on the eigenvalues of Jacobian matrix at the given equilibrium point.Proposition 1.
The boundary equilibrium pointE1 is an unstable node.Proof.
The Jacobian matrix atE1=0,0 is given as follows:(15)JE1=1+4α12−ba−c4−b22001+4α22−ba−c4−b22.
Obviously, the Jacobian matrixJE1 is a diagonal matrix; thus, its corresponding distinct eigenvalues are(16)λ1=1+4α12−ba−c4−b22,λ2=1+4α22−ba−c4−b22.
With the constrictsa>c and 0≤b≤1, λ1,2>1 is satisfied, then we can get that the equilibrium point E10,0 is an unstable node.Proposition 2.
The types of local stability of the equilibrium pointE2 can be classified as follows:(1)
If the condition8+4b−γ4−b22<0 meets, the equilibrium point E2 is an unstable node(2)
If the conditions8+4b−γ4−b22>0, 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, and α2>4−b22/22−ba−c meet, then the equilibrium point E2 is an unstable saddle(3)
If the conditions8+4b−γ4−b22>0, 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, and α2<4−b22/22−ba−c meet, then the equilibrium point E2 is a stable node(4)
Else if8+4b−γ4−b22>0 and α1>4−b228−γ4−b22/22−ba−c8+4b−γ4−b22 are satisfied, then the equilibrium point E2 is an unstable pointProof.
The Jacobian matrix atE2 is(17)JE2=1+4α12−ba−c8+4b−γ4−b224−b228−γ4−b22016α2b2−ba−c4−b228−γ4−b221−4α22−ba−c4−b22,where JE2 is a lower triangular matrix; thus, its eigenvalues are main diagonal elements, that is,(18)λ1=1+4α12−ba−c8+4b−γ4−b224−b228−γ4−b22,λ2=1−4α22−ba−c4−b22.
The type of local stability of this boundary equilibrium pointE2 is determined by different values of the parameters:(1)
If8+4b−γ4−b22<0, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22>1; thus, the equilibrium point E2 is an unstable node.(2)
If8+4b−γ4−b22>0 and 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22 are satisfied, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22<1 meets. And if α2>4−b22/22−ba−c, then λ2=1+4α22−ba−c/4−b22>1 meets; hence, the equilibrium point E2 is an unstable saddle.(3)
If8+4b−γ4−b22>0 and 0<α1<4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22<1. And if α2<4−b22/22−ba−c, then λ2=1−4α22−ba−c/4−b22<1, so the equilibrium point E2 is a stable node.(4)
If the parameters meet with8+4b−γ4−b22>0 and α1>4−b228−γ4−b22/22−ba−c8+4b−γ4−b22, then λ1=1+4α12−ba−c8+4b−γ4−b22/4−b228−γ4−b22>1; thus, the equilibrium point E2 is an unstable point.Proposition 3.
The types of equilibrium pointE3 are similar to the boundary equilibrium point E2.Proof.
The Jacobian matrix atE3 can be written as follows:(19)JE3=1−4α12−ba−c4−b2216α1b2−ba−c4−b228−γ4−b2201+4α22−ba−c8+4b−γ4−b224−b228−γ4−b22.
Matrix (19) is an upper triangular matrix and the corresponding eigenvalues are(20)λ1=1−4α12−ba−c4−b22,λ2=1+4α22−ba−c8+4b−γ4−b224−b228−γ4−b22.
In addition, the boundary equilibrium pointsE2 and E3 are symmetrically located around the diagonal line in the plane R2. And the local stability conditions of E3 are symmetrical to the conditions of E2 listed in Proposition 2, so the conclusions of E3 are no longer repeated here.
In sense of economics, the boundary equilibrium points interpret a situation that one of these two firms has withdrawn from the market. The equilibrium pointsE2 and E3 can be regarded as the situation that one of these two firms takes the lead in the market. That is to say, the oligopoly market becomes the monopoly market, and the local stability of equilibrium points implies the stability of economic market in a short term. Actually, neither of the above two market situations is expected, only when these two companies restrict each other can they be conducive to the stable development of the market and the country. That is, the game between firms reaches an equilibrium, named “Nash equilibrium.” Nash equilibrium [41] is on the premise of maximizing its own profits. Meanwhile, it also ensures the stable development of the market.
The Jacobian matrix at the Nash equilibriumE4 of system (10) is formulated as(21)JE4=1+4α12−ba−c+28−2b−γ4−b22x∗4−b224α1bx∗4−b224α2bx∗4−b221+4α22−ba−c+28−4b−γ4−b22x∗4−b22,where x∗=−42−ba−c/8−4b−γ4−b22. To simplify the process of analyses, the auxiliary variables A=42−ba−c/4−b22>0, B=8−γ4−b22/4−b22<0, and C=4b/4−b22>0 are introduced. Thus, the matrix JE4 can be rewritten as follows:(22)JE4=1+α1A+2B−Cx∗α1x∗Cα2x∗C1+α2A+2B−Cx∗.
The characteristic polynomial of matrix (22) is(23)Pλ=λ2−TrJ+DetJ=0,where TrJ and DetJ denote the trace and the determinant of matrix (22) at Nash equilibrium point E4. And the specific formulation of TrJ and DetJ are given as(24)TrJ=2+α1+α2A+2B−Cx∗,DetJ=1+α1A+2B−Cx∗1+α2A+2B−Cx∗−α1α2C2x∗2=1+α1+α2A+2B−Cx∗+α1α2A2+2A2B−Cx∗+4BB−Cx∗2.Proposition 4.
Neimark–Sacker bifurcation will not occur at the Nash equilibrium pointE4.Proof.
The corresponding discriminant of roots of characteristic polynomial (23) is given as(25)Δ=α1−α22A+2B−Cx∗2+4α1α2C2x∗2>0,which means that the dynamical system (10) will not generate complex eigenvalues at the Nash equilibrium point E4; hence, the Neimark–Sacker bifurcation will not occur at the Nash equilibrium point E4.
On the basis of canonical stability analysis, the sufficient condition of local asymptotical stability of Nash equilibrium pointE4 is that the eigenvalues of corresponding Jacobian matrix of E4 is within the unit cycle. According to Jury criterion [15], the condition of local asymptotical stability of Nash equilibrium point E4 can be expressed in detail as(i)
1+TrJ+DetJ=4+α1+α2A+2B−Cx∗+α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0(ii)
1−TrJ+DetJ=α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0(iii)
1−DetJ=−α1+α2A+2B−Cx∗−α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0Proposition 5.
WhenB+C<0 is satisfied, the second condition of Jury criterion holds. That is,(26)1−TrJ+DetJ=α1α2A2+2A2B−Cx∗+4BB−Cx∗2>0.Proof.
First, lety=A2+2B2B−Cx+4BB−Cx2, and then its roots can be computed as x1=−A/2B−C and x2=−A/2B. We recall that A=42−ba−c/4−b22>0, B=8−γ4−b22/4−b22<0, and C=4b/4−b22>0; then, we can obtain(27)B−C=8−4b−γ4−b224−b22<0.
The numerical relationship betweenx1 and x2 is that x1=−A/2B−C<x2=−A/2B. If the condition A2+2B2B−Cx∗+4BB−Cx∗2>0 holds, then x∗satisfies the condition x∗<x1 or x∗>x2. Following the constriction x∗=−42−ba−c/8−γ4−b22=−A/2B−C, x∗=−A/B−C>x1=−A/2B−C can be captured. Hence, x∗<x1 does not hold. Only when the condition x∗>x2 meets, the second condition of Jury criterion holds. The corner case is(28)x∗−x2=−AB−C+A2B=−AB+C2BB−C.
Hence, whenB+C<0 and x∗>x2 meet, the second condition of Jury criterion holds.Proposition 6.
If the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold.Proof.
Letα1=0 of first condition of Jury criterion, and we can get 4+α1+α2A+2B−Cx∗>0. As for A>0, B<0, C>0, B−C<0, and x∗=−A/B−C, we have(29)A+2B−Cx∗=A+−A2B−CB−C=−ABB−C<0.
After this, we can deduce thatα2<−2/A+2B−Cx∗.
LetΔ=A2+2B2B−Cx∗+4BB−Cx∗2>0, if 0<α2≤−A+2B−Cx∗/Δ, the third condition of Jury criterion will hold. That is,(30)4+α1+α2A+2B−Cx∗+α1α2Δ=α2A+2B−Cx∗+α1A+2B−Cx∗+α2Δ<0,
In the same time, it can be easily proved that(31)A+2B−Cx∗Δα22+4Δα2+4A+2B−Cx∗<0.which is satisfied for any value of α2. Following the above formulations, we can get −2/A+2B−Cx∗<−2A+2B−Cx∗/Δ. The inequalities −A+2B−Cx∗/Δ<α2<−2/A+2B−Cx∗ and(32)α1<−4−2α2A+2B−Cx∗2A+2B−Cx∗+α2Δ<−α2A+2B−Cx∗2A+2B−Cx∗+α2Δ,can be captured through jury criterion and inequality (31). Then, we can deduce the conclusion that if the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold. In summary, system (10) has an unique Nash equilibrium point E4. When B+C<0 and α1<−4−2α2A+2B−Cx∗/2A+2B−Cx∗+α2Δ, system (10) is local asymptotical stable at Nash equilibrium point E4. Under the above constrictions of the system parameters, these two firms will reach equilibrium state and get a maximum value for total profit after several games. Then, the market will keep stable development.


## 3. Analysis of Complex Evolution Law Based on Speed of Adjustment

The local stability of boundary equilibrium points and Nash equilibrium point are discussed by the way of analytics. In this part, in order to investigate the intrinsic laws of evolution of system (10) and to mainly analyze its local and global properties, the method of numerical simulation through two-dimensional bifurcation diagram, one-dimensional bifurcation diagram, two-dimensional largest Lyapunov exponent, and basin of attraction is used.


### 3.1. The Routes to Chaos

Firstly, the parametersa,c,b,andγ are fixed as a=61, c=51.5, b=0.85, and γ=2, and then the speeds of adjustment α1 and α2 are chosen as the bifurcation parameters. The global bifurcation route with respect to speed of adjustment α1,α2 is illustrated in Figure 1, where Figure 1(a) is a two-dimensional bifurcation and Figure 1(b) is a two-dimensional largest Lyapunov exponent corresponding to Figure 1(a).Figure 1
(a) The two-dimensional bifurcation with respect toα1,α2 of system (10). (b) The two-dimensional largest Lyapunov exponent corresponding to Figure 1(a).
(a)(b)In the right of Figure1(a) is a color bar with different colors, that is to say, these different colors represent different periods. Specifically speaking, the brown area denotes the stable area of Nash equilibrium point, i.e., period-1. The green area denotes period-2, the orange area denotes period-3, the yellow area denotes period-4, and so on; the color of deep black denotes period more than 30. As a result of the constraints of algorithm, the states of quasi-period, chaotic, and escaping of system (10) cannot be classified clearly, so all three states are marked by dark black in Figure 1(a). In addition, when the largest Lyapunov exponent is equal to zero then the corresponding state of system (10) is quasi-period, and when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (10) is chaotic. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (10) is escaping area. So, the three states mentioned above can be distinguished by this method displayed in Figure 1(b), where the color bar in the right of this figure represents the largest Lyapunov exponent.The Nash equilibrium point corresponding to Figure2(a) is E42.59,2.59. With the increase of α2, the corresponding eigenvalues at the Nash equilibrium point will increase to λ1=0.97376 and λ2=−0.69956. Then, we can deduce that this Nash equilibrium point E4 is a stable node. And more in detail, when the speed of adjustment α2<0.608, then the state of system (10) is at the Nash equilibrium and will keep stable. While if the speed of adjustment satisfy that α2=0.608, then the stability of this Nash equilibrium will change and a flip bifurcation will occur. Once the speed of adjustment α2 meets the condition that α2>0.791, system (10) will enter into a chaotic state. And in the corresponding largest Lyapunov exponent to Figure 2(a), when the largest Lyapunov exponent is less than zero, the state of system (10) is stable, and system (10) enters into the chaotic state when the largest Lyapunov exponent is greater than zero. And when the largest Lyapunov exponent is equal to zero, a bifurcation takes place in system (10).Figure 2
(a) The one-dimensional bifurcation with respect toα2 of system (10). (b) The corresponding largest Lyapunov exponent of Figure 2(a).
(a)(b)As it is known to all, the flip bifurcation is a classical trajectory for economic system to change the state from stable into chaotic. A new way of development for firms will generate when the economy develops to a certain stage, and there is only one way for firms to further develop the economy and they will face a dilemma in the choice of the new and old way. As a result, the route of period-doubling is generated in constant choice between firms in dilemma. When the speeds of adjustment of these two firms are relatively small, then the two firms will be in the Nash equilibrium state and maintain its profit maximizing, respectively. In the same time, these two firms mutually hold each other back. Once the speed of adjustment exceeds the threshold value, the stability of Nash equilibrium point will be destroyed and then the period-doubling will occur. Finally, the constant choice of the firms leads to the disorder state of the market.


### 3.2. The Evolution of Strange Attractor

In the case of a nonlinear dynamic system, its final state inflects the eventual behavior of system because it is a key point to focus on. Furthermore, the final state of the nonlinear system can be indicated by the attractors representing asymptotical behaviors of the nonlinear system if the number of iteration times tends to infinite. In this research, we discuss the impact of adjustment speed of firm 2α2 on eventual behavior of system (10) by investigating the evolution of attractors, where the parameters are chosen as the same as in Figure 1(a), and the adjustment speed of firm 1 is fixed as α1=0.63. The attractor is a period-4 focus with rough selvedge when the parameter equals to 0.625, which can be seen in Figure 3(a). As shown in Figures 3(b) and 3(c), as the adjustment speed of firm 2 increases, this period-4 focus is transformed to four invariant cycles with selvedge via a Neimark–Sacker bifurcation, and these four invariant cycles become more and more larger and the rough selvedge disappears simultaneously. While when the parameter α2 further grows up to 0.657, these four invariant cycles break to form the attractor shown in Figure 3(d); system (10) enters into “periodic window.” And the broken attractor continues to spread then recombine to form four chaotic attractors, which can be seen in Figure 3(e). When the adjustment speed of firm 2 α2 continues to increase until α2=0.7, the four chaotic attractors mentioned above conduct the dynamic behaviors that spread and merge and finally form a piece of chaotic attractor, see Figure 3(f).Figure 3
The process of evolution of attractors in system (10). (a) α2=0.625; (b) α2=0.626; (c) α2=0.656; (d) α2=0.658; (e) α2=0.679; (f) α2=0.7.
(a)(b)(c)(d)(e)(f)The increase of adjustment speed leads to the complexity of eventual behavior of system (10). This process of evolution presents the properties of “certainty” and “irregular development,” where “certainty” means that the way of bifurcation, that is to say, period-4 Neimark–Sacker bifurcation, and “irregular development” means that the interior structure of the attractors, as well as complexity of interior structure of the chaotic attractor, as shown in Figures 3(e) and 3(f). No matter if from phenomenon development or interior mechanism, system (10) presents a high of certainty and complexity. That means the economic system makes mistakes in predicting, and due to the influence of external condition the development of economy cannot be forecasted anymore.


### 3.3. Multistability

In the process of studying the nonlinear dynamic system, we can observe that the stability of equilibrium points may vary with the change of parameters and initial conditions by numerical simulation. This variety can cause single-level or multilevel bifurcation so as to lead to the system displaying complex dynamical behaviors, such as coexistence of attractor, fractal, and chaos. The coexistence of the attractor implies the motion of multistability of system. That is to say, the bifurcation of the nonlinear system varies with the change of the parameters and initial conditions, which will lead to the change of number of solution. And the number of bifurcation solution varies raised a phenomenon of coexistence of multiattractors, so there is a close relationship between the motion of multistability and bifurcation. In addition, the phenomenon of coexistence of the attractor also indicates that the global bifurcation behavior of the system and the discussion of global bifurcation behavior are mainly reflected in two aspects. One is about the variation of the structure, shape, and amount of attractors with parameters, and another is about the variation of the structure and shape of basin of attraction with parameters.Then, we fixed the values of parameters asa=50, c=38, b=0.45, γ=2.3, and α1=0.4925. This set of parameters has a case that the Nash equilibrium point E42.759,2.759 is an unstable node, and the corresponding eigenvalues, respectively, are λ1=1.526 and λ2=1.1877. As shown in Figure 4, with the increase of the adjustment speed α2 of system (10) can lead to the phenomenon about the global bifurcation of structure, shape, and amount of attractors.Figure 4
The evolution of attracting basin of system (10) with respect to the parameter α2. (a) α2=0.4865; (b) α2=0.489; (c) α2=0.496; (d) α2=0.497.
(a)(b)(c)(d)When the adjustment speedα2 is equal to 0.4865, the chaotic attractor coexists with a period-2 attractor, see Figure 4(a), where the red area denotes the escaping area. The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.06358486 and λLE2=−∞; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.04076098 and λLE2=−0.10188996; the corresponding attractor is the chaotic attractor and its attracting area is the white area.When the adjustment speedα2 increases to α2=0.489, the period-4 attractor suddenly adds into the basin of attraction, which can be seen in Figure 4(b). The increase of parameter raises the change of the structure of the chaotic attractor and the number of the attractor, at the same time, also gives rise to the change of corresponding largest Lyapunov exponent of attractor. The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.18013771 and λLE2=−∞; the corresponding attractor is a period-2 attractor. The initial value x1,x2=2,2 corresponds to two Lyapunov exponents, λLE1=−0.04197454 and λLE2=−0.06106757; the corresponding attractor is period-4 attractor. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.06520552 and λLE2=0.00617542; the corresponding chaotic attractor becomes a hyperchaotic attractor. From Figure 4(b), it is easy to know that the attracting area of period-2 attractors is the light blue area, the attracting area of period-4 attractors is the blue area, and the attracting area of super-chaotic attractor is the white area.And if we further increase the value of parameterα2 to α2=0.496, then we can get an interesting phenomenon that the period-4 attractors become the period-8 attractors via the period-doubling bifurcation, and the chaotic attractors constantly enlarge so the scenario of coexistence, as shown in Figure 4(c). The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.26726652 and λLE1=−0.26726652; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value x1,x2=2,2 corresponds to two Lyapunov exponents, λLE1=−0.19091335 and λLE2=−0.29042750; the corresponding attractor is period-4 attractor and its attracting area is the blue area. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.11730683 and λLE2=0.05142154; the corresponding attractor is still hyperchaotic attractor and its attracting area is the white area.And then if we continue to increase the value of parameterα2, the chaotic attractor will contact with the boundary of the blue area and a “contact bifurcation” will occur, which can be seen in Figure 4(c). As shown in Figure 4(d), this “contact bifurcation” can lead to the chaotic attractor burst and then form a “ghost” existing in the former attracting area. Based on above analysis, it is obvious that, with the increase of the parameter α2, the largest Lyapunov exponent of period-2 decreases but the corresponding attracting area is the light blue area keeping invariantly, while the largest Lyapunov exponent of the chaotic attractor increases, where the attractor becomes hyperchaotic from chaotic, and its attracting area is the white area. This white area is constantly occupied by the blue attracting area so as to shrink with the presence of period-4 attractors. In addition, the number of attractors has increased from two to three and the structure of the attractor changes. The last bifurcation form of the chaotic attractor is the “critical bifurcation,” when this chaotic attractor contacts the boundary of its basin of attraction.The parameters are given asa=48.344, c=43.965, b=0.936, γ=1.202, and α1=1.557. A series of special phenomena about the coexistence are discussed in Figure 5, which indicates the one-dimensional bifurcation of different initial values of the adjustment speed ... The initial value corresponds to Figure 5(a) which is x1,x2=1,4, and the initial values corresponds to Figure 5(c) which is x1,x2=3,0.5, where the blue curve represents the R&D effort x1 of the firm 1 and the red curve represents the R&D effort x2 of the firm 2. The diagrams in Figures 5(b) and 5(d) are the largest Lyapunov exponent corresponding to Figures 5(a) and 5(c), respectively. Let us compare the one-dimensional bifurcation Figure 5(a) with Figure 5(c), and we can draw a result that these two bifurcating diagrams have an obvious difference between them when the adjustment speed α2 satisfies this condition as α2∈0,0.5∪1.536,1.6, which means that due to the different choices of initial values, a huge difference of the eventual state of the system will arise. And if we constantly increase the value of the parameter, α2 tends to a relative large number; then, firm 1 and firm 2 will withdraw from the market in turn. Especially, the way of bifurcation of system (10) is a nonnormal bifurcation in Figure 5(c) when parameter α2 belongs to the interval 0.4,0.5. As shown in the local enlargement of Figure 5(d), the corresponding largest Lyapunov will present a form of oscillation in this interval and the detailed discussion is shown in Figure 6.Figure 5
The bifurcation of system (10) with respect to parameter α2. (a) The initial value is given as x1,x2=1,4. (b) The corresponding largest Lyapunov exponent to Figure 5(a). (c) The initial value is given as x1,x2=3,0.5. (d) The corresponding largest Lyapunov exponent to Figure 5(c).
(a)(b)(c)(d)Figure 6
The evolution of basin of attraction of system (10) with respect to parameter α2. (a) α2=0.43; (b) α2=0.468; (c) α2=0.472; (d) α2=0.48; (e) α2=1.56; (f) α2=1.6.
(a)(b)(c)(d)(e)(f)In the case of this group of parameters, the boundary equilibrium pointE2=0,4.99 is a stable node and the corresponding eigenvalues are λ1=0.98888 and λ2=0.1788, respectively. When the parameter α2 is α2=0.43, the phenomenon of period-2 attractors coexisting with the stable boundary equilibrium E2 is shown in Figure 6(a), where the red area represents the escaping area. And the attractor corresponding to the initial values x1,x2=1,4 is this boundary equilibrium point E2, and its attracting area is the white area. The attractor corresponding to the initial values x1,x2=3,0.5 is the period-2 attractor, and its attracting area is the light blue area, where the diagram in right upper corner of Figure 6(a) is the enlargement of period-2 attractors.When the parameterα2 is α2=0.468, the period-2 attractors become the chaotic attractors through a nonnormal bifurcation, see Figure 6(b), where the corresponding largest Lyapunov exponents are λLE1=0.0251229 and λLE2=−0.3369075. From the enlargement diagram in the right upper corner of Figure 6(b), we can clearly see these chaotic attractors have an obvious level of structure, by using the colors of green, red, and black separately represent the inner, middle, and outer layers of the attractor. Each layer is surrounded by numerous points, and the outer layer is not connected with the inner layer and the middle layer, but the boundary of the inner layer and the inner layer is nested. In addition, with increase of the adjustment speed α2, the white attracting area corresponding to the boundary equilibrium point E2 shrinks much more smaller, while the light blue attracting area expands more and more larger. When the value of adjustment speed α2 grows up to 0.472, the three layers of chaotic attractors split from each other and then the outer expands outward continuously to form the attractor, as shown in Figure 5(c). In this case, the attractors are periodic attractors, and the corresponding largest Lyapunov exponent are λLE1=−0.00826758 and λLE2=−0.05828646. And let the value of parameter α2 further increase, then the inner and middle layers of the periodic attractors vanish. Finally, the outer layer constantly expands outward and mutually merge into the periodic attractors, as shown in Figure 6(c), where the corresponding largest Lyapunov exponent are λLE1=−0.00731851 and λLE2=−0.0565693.With the increase of adjustment speedα2, the periodic attractor will contact its attracting area so as that “critical bifurcation” occurs. And until the value of adjustment speed α2 increases to 1.56, the boundary equilibrium point E2 becomes the boundary chaotic attractor represented by the blue straight line via flip bifurcation, as shown in Figure 6(e). In the same time, the periodic attractor undergoes the “critical bifurcation” and then reforms the boundary chaotic attractor represented by the black straight line, where the attracting area corresponding to the chaotic attractor represented by the blue straight line is the white area, and the attracting area of another boundary chaotic attractor is the light blue area. In addition, these two different attracting areas locate symmetric along the diagonal. If the value of parameter α2 continues to increase, then the basin of attraction will go through a bifurcation, named “global bifurcation,” and the chaotic attractor represented by blue line disappears, forming a “hole,” as shown in Figure 6(f) in the light blue attracting area.From the study above, we can deduce that there is no relationship between the bifurcation of changing qualitative property of basin of attraction and the bifurcation of changing qualitative property of the attractor. After the local bifurcation, the game between the firms will not converge to Nash equilibrium point representing the global optimal strategy anymore. Even though this game starting near Nash equilibrium point, system (10) will also tend to different attractors, and these attractors may be periodic or aperiodic. Furthermore, these bifurcations may also lead to the lack of predictabilities of system (10). The global bifurcation with respect to the boundary of basin of attraction increases the uncertainty of the game fate of the firms as given in the initial strategy. Hence, a slight variation in initial condition of the firms or a small external disturbance in the process of adjustment, will cause a huge influence on the evolution behaviors of the firms in the long run.


### 3.4. The Control of Chaos

With the increase of the speed of adjustment, the oligopoly market will enter into the disorder and chaotic state, and the nonlinear dynamic system is at the chaotic state. It is harmful to the development of the firms when this kind of fierce economic fluctuation serves in the market. For the sake of maintaining the longrun and stable development of two firms in this market, we can take control of system (10), aiming at expanding the stable range at Nash equilibrium of system (10) by introducing a control factor. Based on the opinions above, we exploit the method of state feedback and parameter variable control [40] to take control of system (10). And then the form of system (10) after conducting the control factor can be written as follows:(33)x1t+1=1−μ1x1t+α1x1t4a−c2−b+2x1t−bx2t4−b22−γx1t+μ1x1t,x2t+1=1−μ1x2t+α2x2t4a−c2−b+2x2t−bx1t4−b22−γx2t+μ2x2t,where the parameters μ1 and μ2 are control factors and the restrictive condition is μ1,μ2∈0,1.Firstly, the parameters are given asa=61, c=51.5, b=0.85, γ=2, α1=0.87, and α2=0.87, and we choose the control factors μ1 and μ2 as the parameters of bifurcation. The route of global bifurcation with respect to the control parameters μ1,μ2 is indicated in Figure 7, where Figure 7(a) is a two-dimensional bifurcation and Figure 7(b) is the largest Lyapunov exponent corresponding to Figure 7(a).Figure 7
(a) The two-dimensional bifurcation of system (33) with respect to a pair of parameters μ1,μ2. (b) The corresponding two-dimensional largest Lyapunov exponent to Figure 7(a).
(a)(b)It can be observed from Figure7(a) that there is only one route for system (33) which leads to the Nash equilibrium from chaotic, that is, inverse flip bifurcation. The control parameters μ1 and μ2 successively pass through the dark black area, yellow area, green area, and eventually enter into brown area, which implies that system (33) changes from chaotic to period-4, period-2, and other period-halving and finally reaches the Nash equilibrium. Figure 7(b) is a diagram of the largest Lyapunov exponent corresponding to Figure 7(a), where the color bar in the right of this figure represents the largest Lyapunov exponent. When the largest Lyapunov exponent is equal to zero then the corresponding state of system (33) is quasi-period, that is, the dark grey area. And when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (33) is the chaotic state, that is, grey area. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (33) is escaping trajectories, that is, the white area. Fix the control factor as μ1=0.4 and μ1=0.65; then, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent are displayed in Figures 8(a) and 8(b), respectively, where the red curve represents the R&D effort x1 of firm 1, the blue curve represents the R&D effort x2 of firm 2, and the black curve represents the largest Lyapunov exponent (abbreviated as Lyp).Figure 8
The bifurcation and the largest Lyapunov exponent of system (33) with respect to the parameter μ2. (a) μ1=0.4 and (b) μ1=0.65.
(a)(b)From Figure8, we can observe that when control factor μ1 is fixed as μ1=0.4, and if another control factor satisfies the condition that 0<μ2<0.168, then system (33) is at the escaping state, that is, the market disappears. If the value of μ2 meets 0.168<μ2<0.232, then system (33) is at the chaotic state. If the value of μ2 meets 0.232<μ2<0.478, then system (33) conducts an inverse flip bifurcation. And if μ2>0.478 is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of μ2 meets 0.168<μ2<0.232, then Lyp>0. The corresponding largest Lyapunov exponent is less than zero when μ2>0.232. When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place a bifurcation. Another case is that the control factor μ1 increases to 0.65. In this case, when the condition 0<μ2<0.139 is met, then system (33) is at the escaping state. When the condition 0.139<μ2<0.141is met, then system (33) is at the chaotic state. When the condition 0.141<μ2<0.347 is satisfied, then system (33) conducts an inverse flip bifurcation. And if μ2>0.347 is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of μ2 meets 0.168<μ2<0.141 then Lyp>0. The corresponding largest Lyapunov exponent is less than zero when μ2>0.141. When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place at bifurcation. In addition, we can also know that the larger the control factor μ2 is, the faster system (33) enters into Nash equilibrium.For sake of eliminating or slowing down the speed of system (10) entering into chaos, the control factors are introduced and play a stabilizing role in system (33). We will illustrate how they affect the behavior of system (10) entering into chaos through the control factors and the influence on the adjustment speeds α1 and α2. As shown in Figure 9, a two-dimensional bifurcation with respect to adjustment speeds α1,α2 when the parameters are fixed as a=61, c=51.5, b=0.85, and γ=2.Figure 9
The two-dimensional bifurcation of system (33) with respect to a pair of parameters α1,α2. (a) μ1=μ2=0; (b) μ1=μ2=0.3; (c) μ1=μ2=0.5; (d) μ1=μ2=0.95.
(a)(b)(c)(d)In this group of parameters, the Nash equilibrium point isE42.59,2.59. The increase of control factors and adjustment speeds have no impacts on Nash equilibrium point E4. It can be seen from Figure 9 that the stability region with respect to α1,α2 gradually expands at the increase of control factors. However, it does not affect the global shape of two-dimensional bifurcation with respect to α1,α2, but increases the stability region of the Nash equilibrium. The purpose of introducing the control factor is to expand the region of the Nash equilibrium point attracting area, so as to further delay the speed of system (10) to enter into chaos with increasing of parameters, so as to keep the market stable as long as possible.Furthermore, we will illustrate the influence of control factors on local bifurcation of system (33). Under the action of control factors, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent about adjustment speed α2 are shown in Figure 10, where the parameters are the same as in Figure 9, and the adjustment speed of firm 1 is fixed as α1=0.8. In Figure 10, the blue curve represents R&D effort x1 of firm 1, the red curve represents R&D effort x2 of firm 2, and the local enlargement is the one-dimensional bifurcation of R&D effort x1 with respect to adjustment speed α2. When the control factors are chosen as μ1=μ2=0, system (33) is equivalent to system (10). At present, if the adjustment speed of firm 2 is less than 0.465, that is, α2<0.465, then system (33) will enter into the chaotic state and the market enters into a disorder state. And if α2>0.465, then system is at escaping state, and the corresponding largest Lyapunov exponent is greater than zero, that is, Lyp>0. If the largest Lyapunov exponent is less than zero, that is, Lyp<0, then system (33) is at the period state, see Figure 10(a) for more details. When the control factors are chosen as μ1=μ2=0.3, and the adjustment speed of firm 2 is chosen in the interval that 0<α2<0.534, then system (33) is at the Nash equilibrium state, see Figure 10(b). System (33) will take place a bifurcation, i.e., flip bifurcation, when the adjustment speed of firm 2 is chosen in the interval that 0.534<α2<0.934. When the condition 0.934<α2<1.05 is satisfied, system (33) enters into chaos. Lastly, when the condition α2>1.05 is met, system (33) is at the escaping state. The corresponding largest Lyapunov exponent is negative when 0<α2<0.934. And when α2>0.934, the corresponding largest Lyapunov exponent is positive, that is, Lyp>0. In addition, if we further increase the value of control factors to 0.5 (that is, μ1=μ2=0.5), then we can get the results as follows. System (33) is at the Nash equilibrium when 0<α2<1.1. When 1.1<α2<1.453, system (33) takes place at bifurcation. When 1.453<α2<1.49, then system (33) will enter into chaos. And when α2>1.49, then system (33) is at the escaping state. If 0<α2<1.453, then the corresponding largest Lyapunov exponent is negative, that is, Lyp<0. If α2>1.453, then the corresponding largest Lyapunov exponent is positive, that is,Lyp>0. The detailed dynamical behaviors under μ1=μ2=0.5are shown in Figure 10(c). While when 0<α2<12.34 and μ1=μ2=0.95, system (33) is in a stable period-1 state, see Figure 10(d). From Figure 10(d), we find that system (33) takes place at a series of bifurcations when 12.24<α2<15.88. And when α2>15.88, then system (33) is at the escaping state. If 0<α2<15.88, then the corresponding largest Lyapunov exponent is negative, that is, Lyp<0. The control factors are increased from 0 to 0.95, the stability of system (33) gradually gets stronger and the local bifurcation point is increased from 0.534 to 12.34. It means that the increase of control factors are more beneficial for system (33) to maintain stability, delay the bifurcation point, and keep local stability of the Nash equilibrium point in the long run. However, as for the firms, any slight change of parameters may lead to a huge fluctuation in the market, which makes it difficult for the firms to keep pace with or catch the changing trend of the market. While this fluctuation can be eliminated or slowed down by introducing a control strategy, it is a benefit to a long-term stable development for the market.Figure 10
The bifurcation and the largest Lyapunov exponent of system (33) with respect to parameter α2. (a) μ1=μ2=0; (b) μ1=μ2=0.3; (c) μ1=μ2=0.5; (d) μ1=μ2=0.95.
(a)(b)(c)(d)


## 3.1. The Routes to Chaos

Firstly, the parametersa,c,b,andγ are fixed as a=61, c=51.5, b=0.85, and γ=2, and then the speeds of adjustment α1 and α2 are chosen as the bifurcation parameters. The global bifurcation route with respect to speed of adjustment α1,α2 is illustrated in Figure 1, where Figure 1(a) is a two-dimensional bifurcation and Figure 1(b) is a two-dimensional largest Lyapunov exponent corresponding to Figure 1(a).Figure 1
(a) The two-dimensional bifurcation with respect toα1,α2 of system (10). (b) The two-dimensional largest Lyapunov exponent corresponding to Figure 1(a).
(a)(b)In the right of Figure1(a) is a color bar with different colors, that is to say, these different colors represent different periods. Specifically speaking, the brown area denotes the stable area of Nash equilibrium point, i.e., period-1. The green area denotes period-2, the orange area denotes period-3, the yellow area denotes period-4, and so on; the color of deep black denotes period more than 30. As a result of the constraints of algorithm, the states of quasi-period, chaotic, and escaping of system (10) cannot be classified clearly, so all three states are marked by dark black in Figure 1(a). In addition, when the largest Lyapunov exponent is equal to zero then the corresponding state of system (10) is quasi-period, and when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (10) is chaotic. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (10) is escaping area. So, the three states mentioned above can be distinguished by this method displayed in Figure 1(b), where the color bar in the right of this figure represents the largest Lyapunov exponent.The Nash equilibrium point corresponding to Figure2(a) is E42.59,2.59. With the increase of α2, the corresponding eigenvalues at the Nash equilibrium point will increase to λ1=0.97376 and λ2=−0.69956. Then, we can deduce that this Nash equilibrium point E4 is a stable node. And more in detail, when the speed of adjustment α2<0.608, then the state of system (10) is at the Nash equilibrium and will keep stable. While if the speed of adjustment satisfy that α2=0.608, then the stability of this Nash equilibrium will change and a flip bifurcation will occur. Once the speed of adjustment α2 meets the condition that α2>0.791, system (10) will enter into a chaotic state. And in the corresponding largest Lyapunov exponent to Figure 2(a), when the largest Lyapunov exponent is less than zero, the state of system (10) is stable, and system (10) enters into the chaotic state when the largest Lyapunov exponent is greater than zero. And when the largest Lyapunov exponent is equal to zero, a bifurcation takes place in system (10).Figure 2
(a) The one-dimensional bifurcation with respect toα2 of system (10). (b) The corresponding largest Lyapunov exponent of Figure 2(a).
(a)(b)As it is known to all, the flip bifurcation is a classical trajectory for economic system to change the state from stable into chaotic. A new way of development for firms will generate when the economy develops to a certain stage, and there is only one way for firms to further develop the economy and they will face a dilemma in the choice of the new and old way. As a result, the route of period-doubling is generated in constant choice between firms in dilemma. When the speeds of adjustment of these two firms are relatively small, then the two firms will be in the Nash equilibrium state and maintain its profit maximizing, respectively. In the same time, these two firms mutually hold each other back. Once the speed of adjustment exceeds the threshold value, the stability of Nash equilibrium point will be destroyed and then the period-doubling will occur. Finally, the constant choice of the firms leads to the disorder state of the market.


## 3.2. The Evolution of Strange Attractor

In the case of a nonlinear dynamic system, its final state inflects the eventual behavior of system because it is a key point to focus on. Furthermore, the final state of the nonlinear system can be indicated by the attractors representing asymptotical behaviors of the nonlinear system if the number of iteration times tends to infinite. In this research, we discuss the impact of adjustment speed of firm 2α2 on eventual behavior of system (10) by investigating the evolution of attractors, where the parameters are chosen as the same as in Figure 1(a), and the adjustment speed of firm 1 is fixed as α1=0.63. The attractor is a period-4 focus with rough selvedge when the parameter equals to 0.625, which can be seen in Figure 3(a). As shown in Figures 3(b) and 3(c), as the adjustment speed of firm 2 increases, this period-4 focus is transformed to four invariant cycles with selvedge via a Neimark–Sacker bifurcation, and these four invariant cycles become more and more larger and the rough selvedge disappears simultaneously. While when the parameter α2 further grows up to 0.657, these four invariant cycles break to form the attractor shown in Figure 3(d); system (10) enters into “periodic window.” And the broken attractor continues to spread then recombine to form four chaotic attractors, which can be seen in Figure 3(e). When the adjustment speed of firm 2 α2 continues to increase until α2=0.7, the four chaotic attractors mentioned above conduct the dynamic behaviors that spread and merge and finally form a piece of chaotic attractor, see Figure 3(f).Figure 3
The process of evolution of attractors in system (10). (a) α2=0.625; (b) α2=0.626; (c) α2=0.656; (d) α2=0.658; (e) α2=0.679; (f) α2=0.7.
(a)(b)(c)(d)(e)(f)The increase of adjustment speed leads to the complexity of eventual behavior of system (10). This process of evolution presents the properties of “certainty” and “irregular development,” where “certainty” means that the way of bifurcation, that is to say, period-4 Neimark–Sacker bifurcation, and “irregular development” means that the interior structure of the attractors, as well as complexity of interior structure of the chaotic attractor, as shown in Figures 3(e) and 3(f). No matter if from phenomenon development or interior mechanism, system (10) presents a high of certainty and complexity. That means the economic system makes mistakes in predicting, and due to the influence of external condition the development of economy cannot be forecasted anymore.


## 3.3. Multistability

In the process of studying the nonlinear dynamic system, we can observe that the stability of equilibrium points may vary with the change of parameters and initial conditions by numerical simulation. This variety can cause single-level or multilevel bifurcation so as to lead to the system displaying complex dynamical behaviors, such as coexistence of attractor, fractal, and chaos. The coexistence of the attractor implies the motion of multistability of system. That is to say, the bifurcation of the nonlinear system varies with the change of the parameters and initial conditions, which will lead to the change of number of solution. And the number of bifurcation solution varies raised a phenomenon of coexistence of multiattractors, so there is a close relationship between the motion of multistability and bifurcation. In addition, the phenomenon of coexistence of the attractor also indicates that the global bifurcation behavior of the system and the discussion of global bifurcation behavior are mainly reflected in two aspects. One is about the variation of the structure, shape, and amount of attractors with parameters, and another is about the variation of the structure and shape of basin of attraction with parameters.Then, we fixed the values of parameters asa=50, c=38, b=0.45, γ=2.3, and α1=0.4925. This set of parameters has a case that the Nash equilibrium point E42.759,2.759 is an unstable node, and the corresponding eigenvalues, respectively, are λ1=1.526 and λ2=1.1877. As shown in Figure 4, with the increase of the adjustment speed α2 of system (10) can lead to the phenomenon about the global bifurcation of structure, shape, and amount of attractors.Figure 4
The evolution of attracting basin of system (10) with respect to the parameter α2. (a) α2=0.4865; (b) α2=0.489; (c) α2=0.496; (d) α2=0.497.
(a)(b)(c)(d)When the adjustment speedα2 is equal to 0.4865, the chaotic attractor coexists with a period-2 attractor, see Figure 4(a), where the red area denotes the escaping area. The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.06358486 and λLE2=−∞; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.04076098 and λLE2=−0.10188996; the corresponding attractor is the chaotic attractor and its attracting area is the white area.When the adjustment speedα2 increases to α2=0.489, the period-4 attractor suddenly adds into the basin of attraction, which can be seen in Figure 4(b). The increase of parameter raises the change of the structure of the chaotic attractor and the number of the attractor, at the same time, also gives rise to the change of corresponding largest Lyapunov exponent of attractor. The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.18013771 and λLE2=−∞; the corresponding attractor is a period-2 attractor. The initial value x1,x2=2,2 corresponds to two Lyapunov exponents, λLE1=−0.04197454 and λLE2=−0.06106757; the corresponding attractor is period-4 attractor. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.06520552 and λLE2=0.00617542; the corresponding chaotic attractor becomes a hyperchaotic attractor. From Figure 4(b), it is easy to know that the attracting area of period-2 attractors is the light blue area, the attracting area of period-4 attractors is the blue area, and the attracting area of super-chaotic attractor is the white area.And if we further increase the value of parameterα2 to α2=0.496, then we can get an interesting phenomenon that the period-4 attractors become the period-8 attractors via the period-doubling bifurcation, and the chaotic attractors constantly enlarge so the scenario of coexistence, as shown in Figure 4(c). The initial value x1,x2=2,3.3 corresponds to two Lyapunov exponents, λLE1=−0.26726652 and λLE1=−0.26726652; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value x1,x2=2,2 corresponds to two Lyapunov exponents, λLE1=−0.19091335 and λLE2=−0.29042750; the corresponding attractor is period-4 attractor and its attracting area is the blue area. The initial value x1,x2=2,2.4 corresponds to two Lyapunov exponents, λLE1=0.11730683 and λLE2=0.05142154; the corresponding attractor is still hyperchaotic attractor and its attracting area is the white area.And then if we continue to increase the value of parameterα2, the chaotic attractor will contact with the boundary of the blue area and a “contact bifurcation” will occur, which can be seen in Figure 4(c). As shown in Figure 4(d), this “contact bifurcation” can lead to the chaotic attractor burst and then form a “ghost” existing in the former attracting area. Based on above analysis, it is obvious that, with the increase of the parameter α2, the largest Lyapunov exponent of period-2 decreases but the corresponding attracting area is the light blue area keeping invariantly, while the largest Lyapunov exponent of the chaotic attractor increases, where the attractor becomes hyperchaotic from chaotic, and its attracting area is the white area. This white area is constantly occupied by the blue attracting area so as to shrink with the presence of period-4 attractors. In addition, the number of attractors has increased from two to three and the structure of the attractor changes. The last bifurcation form of the chaotic attractor is the “critical bifurcation,” when this chaotic attractor contacts the boundary of its basin of attraction.The parameters are given asa=48.344, c=43.965, b=0.936, γ=1.202, and α1=1.557. A series of special phenomena about the coexistence are discussed in Figure 5, which indicates the one-dimensional bifurcation of different initial values of the adjustment speed ... The initial value corresponds to Figure 5(a) which is x1,x2=1,4, and the initial values corresponds to Figure 5(c) which is x1,x2=3,0.5, where the blue curve represents the R&D effort x1 of the firm 1 and the red curve represents the R&D effort x2 of the firm 2. The diagrams in Figures 5(b) and 5(d) are the largest Lyapunov exponent corresponding to Figures 5(a) and 5(c), respectively. Let us compare the one-dimensional bifurcation Figure 5(a) with Figure 5(c), and we can draw a result that these two bifurcating diagrams have an obvious difference between them when the adjustment speed α2 satisfies this condition as α2∈0,0.5∪1.536,1.6, which means that due to the different choices of initial values, a huge difference of the eventual state of the system will arise. And if we constantly increase the value of the parameter, α2 tends to a relative large number; then, firm 1 and firm 2 will withdraw from the market in turn. Especially, the way of bifurcation of system (10) is a nonnormal bifurcation in Figure 5(c) when parameter α2 belongs to the interval 0.4,0.5. As shown in the local enlargement of Figure 5(d), the corresponding largest Lyapunov will present a form of oscillation in this interval and the detailed discussion is shown in Figure 6.Figure 5
The bifurcation of system (10) with respect to parameter α2. (a) The initial value is given as x1,x2=1,4. (b) The corresponding largest Lyapunov exponent to Figure 5(a). (c) The initial value is given as x1,x2=3,0.5. (d) The corresponding largest Lyapunov exponent to Figure 5(c).
(a)(b)(c)(d)Figure 6
The evolution of basin of attraction of system (10) with respect to parameter α2. (a) α2=0.43; (b) α2=0.468; (c) α2=0.472; (d) α2=0.48; (e) α2=1.56; (f) α2=1.6.
(a)(b)(c)(d)(e)(f)In the case of this group of parameters, the boundary equilibrium pointE2=0,4.99 is a stable node and the corresponding eigenvalues are λ1=0.98888 and λ2=0.1788, respectively. When the parameter α2 is α2=0.43, the phenomenon of period-2 attractors coexisting with the stable boundary equilibrium E2 is shown in Figure 6(a), where the red area represents the escaping area. And the attractor corresponding to the initial values x1,x2=1,4 is this boundary equilibrium point E2, and its attracting area is the white area. The attractor corresponding to the initial values x1,x2=3,0.5 is the period-2 attractor, and its attracting area is the light blue area, where the diagram in right upper corner of Figure 6(a) is the enlargement of period-2 attractors.When the parameterα2 is α2=0.468, the period-2 attractors become the chaotic attractors through a nonnormal bifurcation, see Figure 6(b), where the corresponding largest Lyapunov exponents are λLE1=0.0251229 and λLE2=−0.3369075. From the enlargement diagram in the right upper corner of Figure 6(b), we can clearly see these chaotic attractors have an obvious level of structure, by using the colors of green, red, and black separately represent the inner, middle, and outer layers of the attractor. Each layer is surrounded by numerous points, and the outer layer is not connected with the inner layer and the middle layer, but the boundary of the inner layer and the inner layer is nested. In addition, with increase of the adjustment speed α2, the white attracting area corresponding to the boundary equilibrium point E2 shrinks much more smaller, while the light blue attracting area expands more and more larger. When the value of adjustment speed α2 grows up to 0.472, the three layers of chaotic attractors split from each other and then the outer expands outward continuously to form the attractor, as shown in Figure 5(c). In this case, the attractors are periodic attractors, and the corresponding largest Lyapunov exponent are λLE1=−0.00826758 and λLE2=−0.05828646. And let the value of parameter α2 further increase, then the inner and middle layers of the periodic attractors vanish. Finally, the outer layer constantly expands outward and mutually merge into the periodic attractors, as shown in Figure 6(c), where the corresponding largest Lyapunov exponent are λLE1=−0.00731851 and λLE2=−0.0565693.With the increase of adjustment speedα2, the periodic attractor will contact its attracting area so as that “critical bifurcation” occurs. And until the value of adjustment speed α2 increases to 1.56, the boundary equilibrium point E2 becomes the boundary chaotic attractor represented by the blue straight line via flip bifurcation, as shown in Figure 6(e). In the same time, the periodic attractor undergoes the “critical bifurcation” and then reforms the boundary chaotic attractor represented by the black straight line, where the attracting area corresponding to the chaotic attractor represented by the blue straight line is the white area, and the attracting area of another boundary chaotic attractor is the light blue area. In addition, these two different attracting areas locate symmetric along the diagonal. If the value of parameter α2 continues to increase, then the basin of attraction will go through a bifurcation, named “global bifurcation,” and the chaotic attractor represented by blue line disappears, forming a “hole,” as shown in Figure 6(f) in the light blue attracting area.From the study above, we can deduce that there is no relationship between the bifurcation of changing qualitative property of basin of attraction and the bifurcation of changing qualitative property of the attractor. After the local bifurcation, the game between the firms will not converge to Nash equilibrium point representing the global optimal strategy anymore. Even though this game starting near Nash equilibrium point, system (10) will also tend to different attractors, and these attractors may be periodic or aperiodic. Furthermore, these bifurcations may also lead to the lack of predictabilities of system (10). The global bifurcation with respect to the boundary of basin of attraction increases the uncertainty of the game fate of the firms as given in the initial strategy. Hence, a slight variation in initial condition of the firms or a small external disturbance in the process of adjustment, will cause a huge influence on the evolution behaviors of the firms in the long run.


## 3.4. The Control of Chaos

With the increase of the speed of adjustment, the oligopoly market will enter into the disorder and chaotic state, and the nonlinear dynamic system is at the chaotic state. It is harmful to the development of the firms when this kind of fierce economic fluctuation serves in the market. For the sake of maintaining the longrun and stable development of two firms in this market, we can take control of system (10), aiming at expanding the stable range at Nash equilibrium of system (10) by introducing a control factor. Based on the opinions above, we exploit the method of state feedback and parameter variable control [40] to take control of system (10). And then the form of system (10) after conducting the control factor can be written as follows:(33)x1t+1=1−μ1x1t+α1x1t4a−c2−b+2x1t−bx2t4−b22−γx1t+μ1x1t,x2t+1=1−μ1x2t+α2x2t4a−c2−b+2x2t−bx1t4−b22−γx2t+μ2x2t,where the parameters μ1 and μ2 are control factors and the restrictive condition is μ1,μ2∈0,1.Firstly, the parameters are given asa=61, c=51.5, b=0.85, γ=2, α1=0.87, and α2=0.87, and we choose the control factors μ1 and μ2 as the parameters of bifurcation. The route of global bifurcation with respect to the control parameters μ1,μ2 is indicated in Figure 7, where Figure 7(a) is a two-dimensional bifurcation and Figure 7(b) is the largest Lyapunov exponent corresponding to Figure 7(a).Figure 7
(a) The two-dimensional bifurcation of system (33) with respect to a pair of parameters μ1,μ2. (b) The corresponding two-dimensional largest Lyapunov exponent to Figure 7(a).
(a)(b)It can be observed from Figure7(a) that there is only one route for system (33) which leads to the Nash equilibrium from chaotic, that is, inverse flip bifurcation. The control parameters μ1 and μ2 successively pass through the dark black area, yellow area, green area, and eventually enter into brown area, which implies that system (33) changes from chaotic to period-4, period-2, and other period-halving and finally reaches the Nash equilibrium. Figure 7(b) is a diagram of the largest Lyapunov exponent corresponding to Figure 7(a), where the color bar in the right of this figure represents the largest Lyapunov exponent. When the largest Lyapunov exponent is equal to zero then the corresponding state of system (33) is quasi-period, that is, the dark grey area. And when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (33) is the chaotic state, that is, grey area. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (33) is escaping trajectories, that is, the white area. Fix the control factor as μ1=0.4 and μ1=0.65; then, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent are displayed in Figures 8(a) and 8(b), respectively, where the red curve represents the R&D effort x1 of firm 1, the blue curve represents the R&D effort x2 of firm 2, and the black curve represents the largest Lyapunov exponent (abbreviated as Lyp).Figure 8
The bifurcation and the largest Lyapunov exponent of system (33) with respect to the parameter μ2. (a) μ1=0.4 and (b) μ1=0.65.
(a)(b)From Figure8, we can observe that when control factor μ1 is fixed as μ1=0.4, and if another control factor satisfies the condition that 0<μ2<0.168, then system (33) is at the escaping state, that is, the market disappears. If the value of μ2 meets 0.168<μ2<0.232, then system (33) is at the chaotic state. If the value of μ2 meets 0.232<μ2<0.478, then system (33) conducts an inverse flip bifurcation. And if μ2>0.478 is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of μ2 meets 0.168<μ2<0.232, then Lyp>0. The corresponding largest Lyapunov exponent is less than zero when μ2>0.232. When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place a bifurcation. Another case is that the control factor μ1 increases to 0.65. In this case, when the condition 0<μ2<0.139 is met, then system (33) is at the escaping state. When the condition 0.139<μ2<0.141is met, then system (33) is at the chaotic state. When the condition 0.141<μ2<0.347 is satisfied, then system (33) conducts an inverse flip bifurcation. And if μ2>0.347 is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of μ2 meets 0.168<μ2<0.141 then Lyp>0. The corresponding largest Lyapunov exponent is less than zero when μ2>0.141. When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place at bifurcation. In addition, we can also know that the larger the control factor μ2 is, the faster system (33) enters into Nash equilibrium.For sake of eliminating or slowing down the speed of system (10) entering into chaos, the control factors are introduced and play a stabilizing role in system (33). We will illustrate how they affect the behavior of system (10) entering into chaos through the control factors and the influence on the adjustment speeds α1 and α2. As shown in Figure 9, a two-dimensional bifurcation with respect to adjustment speeds α1,α2 when the parameters are fixed as a=61, c=51.5, b=0.85, and γ=2.Figure 9
The two-dimensional bifurcation of system (33) with respect to a pair of parameters α1,α2. (a) μ1=μ2=0; (b) μ1=μ2=0.3; (c) μ1=μ2=0.5; (d) μ1=μ2=0.95.
(a)(b)(c)(d)In this group of parameters, the Nash equilibrium point isE42.59,2.59. The increase of control factors and adjustment speeds have no impacts on Nash equilibrium point E4. It can be seen from Figure 9 that the stability region with respect to α1,α2 gradually expands at the increase of control factors. However, it does not affect the global shape of two-dimensional bifurcation with respect to α1,α2, but increases the stability region of the Nash equilibrium. The purpose of introducing the control factor is to expand the region of the Nash equilibrium point attracting area, so as to further delay the speed of system (10) to enter into chaos with increasing of parameters, so as to keep the market stable as long as possible.Furthermore, we will illustrate the influence of control factors on local bifurcation of system (33). Under the action of control factors, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent about adjustment speed α2 are shown in Figure 10, where the parameters are the same as in Figure 9, and the adjustment speed of firm 1 is fixed as α1=0.8. In Figure 10, the blue curve represents R&D effort x1 of firm 1, the red curve represents R&D effort x2 of firm 2, and the local enlargement is the one-dimensional bifurcation of R&D effort x1 with respect to adjustment speed α2. When the control factors are chosen as μ1=μ2=0, system (33) is equivalent to system (10). At present, if the adjustment speed of firm 2 is less than 0.465, that is, α2<0.465, then system (33) will enter into the chaotic state and the market enters into a disorder state. And if α2>0.465, then system is at escaping state, and the corresponding largest Lyapunov exponent is greater than zero, that is, Lyp>0. If the largest Lyapunov exponent is less than zero, that is, Lyp<0, then system (33) is at the period state, see Figure 10(a) for more details. When the control factors are chosen as μ1=μ2=0.3, and the adjustment speed of firm 2 is chosen in the interval that 0<α2<0.534, then system (33) is at the Nash equilibrium state, see Figure 10(b). System (33) will take place a bifurcation, i.e., flip bifurcation, when the adjustment speed of firm 2 is chosen in the interval that 0.534<α2<0.934. When the condition 0.934<α2<1.05 is satisfied, system (33) enters into chaos. Lastly, when the condition α2>1.05 is met, system (33) is at the escaping state. The corresponding largest Lyapunov exponent is negative when 0<α2<0.934. And when α2>0.934, the corresponding largest Lyapunov exponent is positive, that is, Lyp>0. In addition, if we further increase the value of control factors to 0.5 (that is, μ1=μ2=0.5), then we can get the results as follows. System (33) is at the Nash equilibrium when 0<α2<1.1. When 1.1<α2<1.453, system (33) takes place at bifurcation. When 1.453<α2<1.49, then system (33) will enter into chaos. And when α2>1.49, then system (33) is at the escaping state. If 0<α2<1.453, then the corresponding largest Lyapunov exponent is negative, that is, Lyp<0. If α2>1.453, then the corresponding largest Lyapunov exponent is positive, that is,Lyp>0. The detailed dynamical behaviors under μ1=μ2=0.5are shown in Figure 10(c). While when 0<α2<12.34 and μ1=μ2=0.95, system (33) is in a stable period-1 state, see Figure 10(d). From Figure 10(d), we find that system (33) takes place at a series of bifurcations when 12.24<α2<15.88. And when α2>15.88, then system (33) is at the escaping state. If 0<α2<15.88, then the corresponding largest Lyapunov exponent is negative, that is, Lyp<0. The control factors are increased from 0 to 0.95, the stability of system (33) gradually gets stronger and the local bifurcation point is increased from 0.534 to 12.34. It means that the increase of control factors are more beneficial for system (33) to maintain stability, delay the bifurcation point, and keep local stability of the Nash equilibrium point in the long run. However, as for the firms, any slight change of parameters may lead to a huge fluctuation in the market, which makes it difficult for the firms to keep pace with or catch the changing trend of the market. While this fluctuation can be eliminated or slowed down by introducing a control strategy, it is a benefit to a long-term stable development for the market.Figure 10
The bifurcation and the largest Lyapunov exponent of system (33) with respect to parameter α2. (a) μ1=μ2=0; (b) μ1=μ2=0.3; (c) μ1=μ2=0.5; (d) μ1=μ2=0.95.
(a)(b)(c)(d)


## 4. Conclusion

In this research, a nonlinear game model of the R&D competition between duopoly firms is studied in detail. The game between these two firms is mainly described by the way of a two-stage game. In the first stage, all of firms determine the level of their R&D efforts to reduce the production cost. And in the second stage, supposing that all firms compete in the market with a Cournot form and determine the outputs to maximize themselves profits, respectively, the local and global properties are discussed by theoretical analysis and numerical simulation. Theoretically, the local stability condition of the equilibrium points and the kind of local bifurcation are analyzed through stability theory and the Jury criterion.In numerical simulation, we choose the adjusting speeds of these two firms as bifurcation parameters. The evolution law of internal complexity of the nonlinear economic system can be explained by the tools as one-dimensional bifurcation, two-dimensional bifurcation, the one-dimensional largest Lyapunov exponent, the two-dimensional largest Lyapunov exponent, and the basin of attraction. The system will transform from the stable Nash equilibrium state to the chaotic state via a flip bifurcation, which can be verified by the diagrams of one-dimensional bifurcation and the corresponding largest Lyapunov exponent. As it is known to all, the flip bifurcation is a classical route for the economic system to change the state from stable into chaotic. A new way of development for firms will generate when the economy develops to a certain stage, and there is only one way for firms to further develop the economy and they will face a dilemma in the choice of the new and old way. As a result, the route of period-doubling is generated in constant choice between firms in dilemma.The increase of adjusting speed leads to the complexity of eventual behavior of system (10). This process of evolution presents the properties of “certainty” and “irregular development,” where “certainty” means that the way of bifurcation, that is to say, period-4 Neimark–Sacker bifurcation, and “irregular development” means that the interior structure of strange attractors. No matter from phenomenon development or interior mechanism, system (10) presents a high of certainty and complexity. That means the economic system makes mistakes in predicting, due to the influence of external condition, the development of economy cannot be forecasted anymore. Both the variation of parameters and the choice of initial condition will cause system (10) to generate multistability. With the increase of adjustment speed of firms, the properties of multistability are mainly reflected in the number and structure of coexisting attractors and global bifurcation of the system. There are four distinct coexisting phenomena displayed here, that is, periodic attractors coexist with periodic attractors, boundary equilibrium points coexist with strange attractors, periodic attractors coexist with chaotic attractors, and multiboundary chaotic attractors coexist. In the case of the boundary attractors coexisting with strange attractors with the layered structure, the largest Lyapunov exponent corresponds to strange attractors displaying a fluctuating and oscillating phenomenon as adjusting speed increases. This is a special phenomenon of coexistence and seldom reported by other researchers. As we all know, the boundary attractor corresponds to the fact that the duopoly market has degenerated into a monopoly market. That is, one of these two firms was forced out of the market. Therefore, the coexistence of the boundary attractor and internal attractor implies that not only the parameters can affect the market structure but also the initial state of market can affect the market structure.The increase of the adjustment speed will lead to change of multistability, and at the same time, it will also cause “critical bifurcation” so as to form “ghost,” as well as the “global bifurcation” of attracting basin so as to form “hole.” The increasing of control factors are more beneficial for system (33) to maintain stability, delay the bifurcation point, and keep local stability of the Nash equilibrium point in the long run. However, as for the firms, any slight change of parameters may lead to a huge fluctuation in the market, which makes it difficult for the firms to keep pace with or catch the changing trend of the market. While this fluctuation can be eliminated or slowed down by introducing a control strategy, it is a benefit to a long-term stable development for the market.



---

*Source: 9634878-2020-08-11.xml*