# On the Stability and Multistability in a Duopoly Game with R&D Spillover and Price Competition

Authors: Wei Zhou; Xiao-Xue Wang
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

In this paper, a dynamical two-stage game with R&D competition and joint profit maximization is built. The stability of all the equilibrium points is discussed through Jury condition, and the stability region of the Nash equilibrium point is then given. The influence of the parameters on the system is discussed, and we find that the firm can even benefit from chaos, when it has higher innovation efficiency and higher adjusting speed. And then the coexistence of multiple attractors is studied using basin of attraction. Our research result shows that the coexisting attractors can be observed in the two-parameter bifurcation diagram. At last, the boundary of feasible region, global bifurcations, and formation mechanism of fractal structure of attracting basin are analyzed through critical curves and noninvertible map theory.


---

## Body

## 1. Introduction

With the rapid development of economics, the competition among firms has become more and more fierce. The firms, which are the most important part of the market, must form a core competitiveness that is different from other firms, if they want to get a good survival and remain invincible in the market. Research and development (or R&D for short) has become the main driving force for the development of enterprises and it is also an important way for enterprises to obtain the core competitiveness. In fact, the R&D activities of enterprises can reduce production cost, improve the quality of products, expand the market share, and further improve its market competitiveness. However, in the real economic environment, R&D spillovers often arise when the R&D activities are carried out. In fact, the R&D spillovers are inevitable due to the information exchange of R&D among firms and the flow of human resources. In recent years, the problem of competition and cooperation in R&D investment has attracted more and more attention of entrepreneurs and economists. The classic research begins in the 1980s; a famous duopoly model with technology spillover was built by D’Aspremont and Jacquemin [1]. In this model, the competition between firms was divided into two stages. At the first stage, both firms decided their own R&D level, while, at the second stage, all the firms determined their outputs to the market. This famous model is also called AJ model by followers. Then, Kamien et al. [2] extended the model of D’Aspremont and Jacquemin to n firms and got the KMZ model. In KMZ model, the forms of R&D are also extended to total competition, R&D cartel, research joint ventures (RJVs), and RJVs cartel.Both AJ model and KMZ model are two-stage game models. Actually, the two-stage game is better than one-stage game for solving the uncertainty of the market and the interaction of competitors’ strategies when the R&D activity is considered. Following the AJ model and the KMZ model, Ziss [3] built a two-stage duopoly game model with R&D spillovers. Amir [4] studied the effects of R&D spillover and pointed out that the R&D cost of firms would increase with an increasing amount of technology outflow. Lambertinti et al. [5] extended the research results of Amir and proposed that the cooperative innovation will increase information sharing of market. Katsoulacos and Ulph [6] have led to a study of endogenous spillovers in R&D cooperation. Yin [7] studied the asymmetric R&D cooperation structure and pointed out that one of the important motivations for providing cooperative innovation is to increase the market share. Alallah [8] analyzed the influence of spillover rate on cooperative motivation and joint profit under asymmetrical R&D spillover. Bischi et al. [9, 10] established a two-stage dynamic game according to knowledge share and R&D competition in the market and analyzed the stability of the built model. Zhang et al. [11] considered a two-stage duopoly model with semi-collusion in production, and the complex chaotic behavior of the given model was researched in their works.The influence of chaos theory on modern science is limited not only to nature science, but also to other fields, such as philosophy [12], social science [13], and economics [14]. The economic system is essentially nonlinear, so the chaotic phenomenon must exist in the economic system. In recent years, it has attracted the attention of a growing number of scholars to explore the evolution of the economic system and explain the complicated economic phenomena by using chaos theory. Gangopadhyay [15] built a dynamic model of enterprise merger, and the bifurcation behavior and the coexistence of multiple attractors in the built model were discussed. Canovas et al. [16] set up a piecewise smooth model for market competition, and the asymptotic dynamic behaviors of this model were discussed. Bischi et al. [17] proposed an inertial model of market share based on the optimal response and naïve expectation, the influence of the heterogeneity of the enterprise was also analyzed, and it was found that there are very rich local and global bifurcation phenomena in the built system. Li and Ma [18] considered a bounded rational dual-channel game, and the complex dynamical behavior of their model was simulated in their works. Cavalli et al. [19] built a nonlinear dynamical duopoly game by using bounded rationality and “LMA” (the abbreviation of “local monopolistic approximation”) adjustment strategies, and then they analyzed the coexistence of periodic points and the “lobes” basin of attraction was formed by periodic points. Based on the research results of [19], Cavalli and Naimzada [20] built a heterogeneous duopoly Cournot game with reduced rationality, isoelastic demand function, and linear total cost. And they found that the Nash equilibrium point may lose its stability through a flip bifurcation or a Neimark-Sacker bifurcation. In addition, a lot of researchers have discussed the complex dynamical behaviors of nonlinear oligopolies from different aspects, such as differentiated goods [21–25], heterogeneous firms [26–31], and delayed decisions [32, 33].Another important issue for dynamical economical model is the coexistence of attractors. During the long run of dynamical economical system, the initial conditions play a very important role in determining the asymptotic behavior of the system. And thus, it is especially important to discuss the dependence of system on the initial conditions using basin of attraction and to analyze the evolution of attracting basin. In recent years, this topic has attracted a lot of attention of many scholars. Bischi et al. [34] discussed the coexistence of multi-attractors and cyclic attractors in a dynamic Cournot game. Agiza et al. [35] analyzed the multi-stability of a three-dimensional dynamical oligopoly model. Both the routes to complex attractors and the creation of basins with complex structures were studied in their paper. Cavalli and Naimzada [19] built an oligopoly model with different rational players, and then they explored the complex dynamical behaviors of the given model induced by coexistence of multiple attractors. Gori and Sodini [36] presented a nonlinear duopoly game with price competition, and the stability and coexistence of attractors were studied using the basin of attraction.Through analyzing the above related research results, it is not difficult to find out that the R&D competition between duopoly firms that produce complementary products is seldom studied. However, the R&D competition between firms that produce complementary products is very common in the real market. A typical example is the cooperation and competition between “small program” and “APP” software. One firm developing “small program” and another firm developing “APP” always compete in R&D level as the development of them has a lot of commonalities. But they collaborate in the production level because they are complementary in functionality. They choose to cooperate in producing and selling their products, so as to maximize their joint profits. Therefore, a two-stage duopoly game under the background of R&D spillover, complementary products, and joint profits maximization is proposed in this present work.In addition, the stability and the complex dynamical behavior in the built two-stage duopoly game is another topic in this research. As the stability of equilibrium point is often influenced by various parameters, the stability condition and stability region with various sets of parameters will be discussed in this article. And the complicated dynamical phenomena maybe arise when the equilibrium points lose their stability. The coexistence of many different attractors is the most important dynamical behaviors, which can be found in our model. The coexistence of attractors (i.e., multistability) indicates the long-run dependence of market dynamics on the initial conditions. Different initial conditions not only will result in a vast difference for firms after several games, but also may cause the coexistence of different situations. So by modeling the practical problem and analyzing the stability and multistability of the built model, it is beneficial for firms to make their decisions and predictions.The plan of this article is given as follows. In the second section, a duopoly two-stage game based on R&D spillover, complementary goods, and maximization of joint profits is built. At the first stage of this game, two firms compete in the R&D level, but the R&D spillovers are permitted in this model. However, at the second stage, these two firms produce complementary goods and choose cooperation at the production level. In Section3, the equilibrium points of this model are solved and the stability conditions are discussed. In Section 4, the stability of the Nash equilibrium point and the bifurcation of this system are analyzed using numerical simulation. And then, in Section 5, the coexistence among many attractors is investigated. The evolution and formulation mechanism of attracting basins are studied. Results are summarized in the last section.


## 2. The Model

Suppose there are two firms, labeled byi (i=1,2), in a market, which conduct R&D and produce complementary goods. In order to better analyze the R&D investment behavior of the two firms, we consider the following two-stage game. In the first stage of the game, both firms compete in the R&D level and choose the R&D efforts as their decision variables to reduce the production cost and improve the quality of their products, while, at the second stage of the game, both firms choose the market prices as their decision variables. And they choose cooperation, share the market information, and pursue the maximization of their joint profit at the production level. Obviously, the R&D investments of the two firms in the first stage will affect their profits in the second stage.Letpi>0 and qi>0 be the price and the quantity of product of firm i (i=1,2), respectively. And the output of firm i is determined by the prices of the goods produced by firm i and its rival through an inverse demand function, which has the form [27](1)qi=ai-bipi-bjpji,j=1,2,i≠jwhere ai>0 represent the market sizes and bi>0, bj>0 are the price sensitivity of consumers.As we know, the spillover of knowledge is always inevitable in the process of R&D. That is, the R&D results of firmi not only can reduce the production cost of its own, but also can reduce the production cost of its rival. So if we suppose xi and xj are the R&D efforts of firms i and j, respectively, then the effective marginal cost of firm i with R&D efforts can be given as(2)Cixi,xj=ci-xi-βxj,i,j=1,2,i≠jwhere ci>0 represents the unit cost of produced goods without R&D efforts. If we assume that the two firms have equal knowledge levels when they are not engaged in R&D, that is, the productive efficiency of these two firms is the same, then we can let ci=c. And β∈0,1 is related to the R&D spillover. If β=0, it means that the technology between firms is completely confidential and the R&D spillover is ignored. And if β=1, it means that R&D technology is fully shared among competing firms. In this paper, we only consider the intermediate condition 0<β<1.However, the firms have to pay certain R&D expenditures when they conduct R&D activities. And the R&D expenditure can be assumed to be a quadratic function of R&D effortxi, which reflects the decreasing returns to R&D effort. And the R&D expenditure of firm i can be given as (3)gxi=12γxi2,i=1,2where γ>0 is the cost parameter of firm’s technological innovation, which indicates the efficiency of using or producing the unique technology or knowledge resources for an enterprise. The smaller the parameter γ, the stronger the innovation ability of firm i.The individual profit of firmi can be given as πi=pi-Ciqi-gxi,(i=1,2). Through substituting (2) and (3) into πi, we can get (4)π1=p1-c-x1-βx2q1-12γx12(5)π2=p2-c-x2-βx1q2-12γx22For convenience, the parametersa1 and a2 are set as the same in the following discussion. That is, a1=a2=a. And the similar hypothesis is also imposed on parameters b1 and b2; it means that b1=b2=b. The marginal profit of firm i can be obtained by taking the first derivative of profits πi with respect to prices pi(i=1,2). That is,(6)∂π1∂p1=a-bp1+p2-bp1-c+x1+βx2(7)∂π2∂p2=a-bp1+p2-bp2-c+x2+βx1Through solving the first order conditions∂πi/∂pi=0(i=1,2), we can get the optimal prices of products produced by firm 1 and firm 2, respectively. And they are given as (8)p1∗=a+bc+b1-2βx2+bβ-2x13b(9)p2∗=a+bc+b1-2βx1+bβ-2x23bBy substituting (8) and (9) into (1), the optimal quantity of firm i can be obtained, (10)qi∗=a-2bc+bβ+1xi+xj3,i,j=1,2,i≠jThen the profit of firmi can be obtained by substituting (8), (9), and (10) into (4) and (5), respectively. (11)π1=a-2bc+bβ+1x1+x29b2-12γx12(12)π2=a-2bc+bβ+1x2+x19b2-12γx22Since, as we have assumed, the two firms choose collaboration in the production level, we consider the joint profit, namely,Π=π1+π2. By combining equation (11) and (12), the specific forms of the joint profit of these two firms can be given.(13)Π=2a-2bc+bβ+1x1+x29b2-12γx12-12γx22The marginal joint marginal profits of each firm can be obtained by taking the first derivative of joint profitΠ with respect to the R&D efforts xi(i=1,2). That is,(14)∂Π∂x1=4b9β+12x1+x2+49a-2bcβ+1-γx1(15)∂Π∂x2=4b9β+12x1+x2+49a-2bcβ+1-γx2As we have presumed above, the two firms choose to cooperate in the production level. Therefore, they will share the market information in order to improve their joint profit. However, the two firms will not share the information in the R&D level, as they are competitive in this process. And the decision-making of them cannot be completely rational, due to the lack of market information. So, the firms are not able to optimize their joint profit with respect to the R&D effort in one shot. They can only adjust their R&D strategies over time towards an optimal strategy along a direction of local estimation of the expected profits gradient, according to gradient adjustment mechanism. For a detailed description of the economic rationale of the gradient adjustment mechanism we refer to the description in [37].According to the gradient adjustment mechanism, we know that if firmi has positive marginal profit in period t, it will increase its R&D effort in periodt+1. On the contrary, if its marginal profit is negative, it will reduce its R&D effort. In the actual market, the firm dynamically adjusts its strategy at each discrete time, and the resulting dynamical system can be described by the following couple of equations:(16)x1t+1=x1t+α1x1t∂Π∂x1x2t+1=x2t+α2x2t∂Π∂x2where xi=xit(i=1,2) is the R&D effort of firm i at discrete time periods t=0,1,2⋯ and αi>0(i=1,2) is the speed of adjustment for firm i(i=1,2). The dynamic game model can be obtained by substituting equations (14) and (15) into (16).(17)x1t+1=x1t+α1x1t4b9β+12x1t+x2t+49a-2bcβ+1-γx1tx2t+1=x2t+α2x2t4b9β+12x1t+x2t+49a-2bcβ+1-γx2t


## 3. Stability Analysis of Equilibrium Point

In this section, the local stability around the equilibrium point will be analyzed through Jury condition [38]. The equilibrium points of system (17) can be solved by setting xit+1=xit=xi, that is, (18)x1+α1x14b9β+12x1+x2+49a-2bcβ+1-γx1=x1x2+α2x24b9β+12x1+x2+49a-2bcβ+1-γx2=x2Through solving the above nonlinear algebraic equations, four equilibrium points of system (17) can be obtained, (19)E0=0,0,E1=0,4A9γ-4Bb,E2=4A9γ-4Bb,0,E3=4A9γ-8Bb,4A9γ-8Bbwhere A=a-2bc1+β, B=β+12. The equilibrium points E0, E1, and E2 are called boundary equilibrium points as they locate at the coordinate axes. And the equilibrium point E3 is called interior equilibrium point, which is also the unique Nash equilibrium. Obviously, if we assume the parameters meet the conditions as a>2bc and 9γ>8Bb, then the economic meaningfulness of these equilibria can be guaranteed.It is obvious that all the coordinate axesxi(i=1,2) are invariant sub-manifold of system (17), since xit=0 implies xit+1=0. It means that any initial conditions starting from coordinate axis will stay on the coordinate axis forever after any time of iteration. And this situation can be regarded as the market competitors have degenerated from duopoly to monopoly. Or we can say that one of the firms has gone out of the market. The eventuality in which one firm exits the market and how the dynamics behave have been discussed by Gori et al. [38]. At this moment, the system (17) can be rewritten as a one-dimensional map,(20)xi′=1+49a-2bcβ+1αixi-αi99γ-4β+12bxi2where xi′ represents the evolution operator of xi. The above equation is obtained by letting xi′=xi(t+1) and xj(t)=xj(t+1)=0 in system (17), where i,j=1,2 and i≠j. Obviously, the above equation is conjugate to the famous Logistic map l′=τl1-l topologically through a proper linear transformation, which is given as(21)xi=9+4a-2bcβ+1αi9γ-4β+12bland the corresponding value of τ is(22)τ=1+49a-2bcβ+1αiAs all the boundary equilibrium points locate on the coordinate axes ofx1,x2, the stability of all the boundary equilibrium points of system (17) is equivalent to the corresponding points of the Logistic map. According to the stability theory, the local stability of an equilibrium point can be determined by analyzing the eigenvalues of Jacobian matrix of system (17) evaluated at the equilibrium point on the complex plane. The Jacobian matrix of system (17) at any point x1,x2 can be given as(23)Jx1,x2=1+α149Bb2x1+x2+4A9-2γx149Bbα1x149Bbα2x21+α249Bbx1+2x2+4A9-2γx2Proposition 1.
The trivial equilibrium pointE0 is an unstable node.Proof.
At the trivial equilibrium pointE0=0,0, the Jacobian matrix (23) becomes a diagonal matrix, which can be given as(24)JE0=1+α14A9001+α24A9Therefore the eigenvalues of this diagonal matrix are the diagonal entries, which are λ1=1+α14A/9 and λ2=1+α24A/9. According to the nonnegativity condition a>2bc, we know that A>0. And thus, both of the eigenvalues are larger than the unity, that is, λ1>1 and λ2>1. And hence the trivial equilibrium point E0 is an unstable node along the direction of the coordinate axes.In the actual market,E0=0,0 is on behalf of the zero R&D efforts of these two firms. It is obvious that this boundary equilibrium means that there is no investment in R&D for both firms. And if the initial condition is chosen from the coordinate axes of the plane x1,x2 (except the origin 0,0), the system will eventually evolve to the equilibrium point E1 or E2 as E0=0,0 is a repelling node. However, the system will stay at E0 forever, if the initial condition is chosen as 0,0, even though the trivial equilibrium E0 is unstable.Proposition 2.
The boundary equilibriums pointE1 is a saddle point when α2<9/2A. But E1 is an unstable node when α2>9/2A.Proof.
At the boundary equilibrium pointE1=0,4A/9γ-4Bb, the Jacobian matrix (23) becomes a lower triangular matrix, which is given as(25)JE1=1+α14A9γ-4Bbγ016A99γ-4BbBbα21-α24A9The eigenvalues of this lower triangular matrix are given by the diagonal entries λ1=1+α14A/9γ-4Bbγ and λ2=1-α24A/9. Through the assumption given above, we know that 9γ>8Bb, so λ1>1. Again with the nonnegative condition we know that A>0, so λ2<1 will also be satisfied. However, when α2>9/2A, the eigenvalue λ2 will be less than -1, so the boundary equilibrium point is an unstable node. But when α2<9/2A, the eigenvalue λ2 will be located in the unit circle, so the boundary equilibrium point E1 is an unstable saddle.It is clear that the boundary equilibrium pointE1=0,4A/9γ-4Bb presents that only firm 2 conducts R&D, while firm 1 is eliminated or left out of the competitive market and therefore it can be viewed as a single monopoly market. As the equilibrium point E2 is symmetric with E1 in the rectangular coordinate system, the stability analysis of E2 is very similar to that of E1 and is not repeated. In order to analyze the stability of the Nash equilibrium point, the following proposition is firstly given.Proposition 3.
The Neimark-Sacker bifurcation cannot occur at the Nash equilibrium pointE3.Proof.
At the Nash equilibrium pointE3, the Jacobian matrix becomes the following form:(26)JE3=1-α1CD4CBbα14CBbα21-α2CDwhere C=4A/99γ-8Bb and D=9γ-4Bb. Let TrE3 and DetE3 be the trace and the determinant of the Jacobian matrix JE3, respectively. The following equation gives the characteristic equation of matrix JE3:(27)λ2-Tr⁡E3λ+DetE3=0where the trace and the determinant of the Jacobian matrix are given as (28)TrE3=2-α1+α2CD(29)DetE3=1-α1+α2CD+4ACγα1α2According to Eq. (28) and Eq. (29), the discriminant of the characteristic polynomial is given as (30)Tr2-4Det=4A99γ-8Bb2α1-α229γ-4Bb2+8Bb2α1α2As the parameters α1 and α2 meet α1>0 and α2>0, the discriminant of the characteristic polynomial will be always nonnegative. It means that there is no complex root for this quadratic equation. And hence, the Neimark-Sacker bifurcation cannot occur.Proposition 4.
The unique Nash equilibriumE3 is stable when the system parameters meet the following conditions: (31)A2Aα1α2γ-9γ-4Bbα1+α2>4.59γ-8BbProof.
According to the Jury condition [38], the Nash equilibrium point E3 is local asymptotic stable if and only if all the following conditions are satisfied.(32)i1-Tr+Det=4ACγα1α2>0ii1+Tr+Det=4-2α1+α2CD+4ACγα1α2>0iii1-Det=α1+α2CD-4ACγα1α2>0
The conditioni apparently holds as all the parameters A, C, γ, α1, and α2 are all positive. And the condition iii is unnecessary as a consequence of Proposition 3. In order to solve the inequality ii, we let 1+tr+Det=0. And then the following threshold of α1 can be gained through solving the given equation. (33)α1∗=CDα2-2C2Aγα2-D=2A9γ-4Bbα2-99γ-8Bb2A2Aγα2-9γ-4Bb
And the inequalityii gives α1>α1∗. Through the discussion given above, the stability region of Nash equilibrium point can be given as(34)A2Aα1α2γ-9γ-4Bbα1+α2>4.59γ-8BbIn system (17), the equilibrium E3 is stable when the parameters satisfy the stability condition, while, when one or more parameters violate the stability condition, the Nash equilibrium E3 will lose its stability through some types of bifurcation. In the next section, the stability of the equilibrium point is studied using numerical simulation. Then the bifurcation with varying parameters and the complex dynamical behaviors of this system when the equilibrium point loses its stability are also analyzed.


## 4. Bifurcation Analysis Using Numerical Simulation

In this section, some numerical simulations will be conducted in order to analyze the stability of the equilibrium point and the complex dynamical behaviors of system (17).Firstly, the system parameters are fixed asa=7.12, b=0.024, c=0.98, β=0.38, and γ=0.44, and the initial condition is chosen as x1,x2=2.33,1.3. The stability region of the Nash equilibrium point E3 is given in Figure 1. In Figure 1, the area with dark color represents the stability region of E3, while the white area represents the instability region of E3. We can also see from Figure 1 that the stability region is situated on the lower left corner of this figure. The shape of the stability region just likes a curved edge rectangle in the parameter plane of α1,α2. And the boundary between stability region and instability region corresponds to a bifurcation curve. As there is no possibility for Neimark-Sacker bifurcation, this bifurcation curve refers to a period doubling bifurcation curve. It means that if the parameter α1 or/and α2 cross this curve into the white area from the dark area, then a period doubling bifurcation occurs. The endpoints of this bifurcation curve also can be solved through the equation 1+Tr+Det=0, if we let α1=0 and α2=0 in this equation, respectively. And the coordinate axis of these two endpoints can be given as(35)A1=99γ-8Bb2A9γ-4Bb,0,A2=0,99γ-8Bb2A9γ-4BbFigure 1
The stability region of the Nash equilibrium pointE3 on the parameter plane α1,α2.In order to further analyze the influence of the system parameters on the stability region of the Nash equilibrium pointE3 on the plane α1,α2, the stability curves with different values of γ are plotted in the plane α1,α2. For example, if we fix all the other parameters as a=7.12, b=0.024, c=0.98, and β=0.38, the value of parameter γ is chosen as 0.1, 0.44, and 1.0, respectively. The stability curves in the parameter plane α1,α2 with different values of γ are shown in Figure 2. In Figure 2, the stability curves with black color, blue color, and red color represent γ=0.1, γ=0.44, and γ=1.0, respectively. And the area formed by these stability curves and coordinate axes are the stability regions with different values of γ. From Figure 2 we can see that the stability region increases very quickly as the value of γ increases from 0.1 to 0.44. However, the change of stability region is very small when γ increases from 0.44 to 1.0. In fact, there is hardly any change in the stability region as γ continues to increase further. And the shape of stability curve is also changing as γ varies; see Figure 2 for details.Figure 2
The stability curves of the Nash equilibrium pointE3 on the parameter plane α1,α2 with different values of γ.The above analysis is mainly about the stability region of the Nash equilibrium point in the parameter planeα1,α2. However, the complex dynamical behaviors outside the stability region have not been revealed. In order to analyze the complex dynamical behaviors of system (17) when the Nash equilibrium point loses its stability, the bifurcation diagram is employed in this article. However, the tool we used here is not only the one-parameter bifurcation diagram but also two-parameter bifurcation diagram, in order to analyze the complex dynamical behavior of system (17) in depth.Figures3(a) and 3(b) show the two-parameter bifurcation diagram with fixed parameter as a=7.12, b=0.024, c=0.98, β=0.38, and γ=0.44 and initial conditions as x1,x2=2.33,1.3 in the parameter plane of α1,α2. The different colors in these figures represent different number of points of the numerically evaluated attractor, which can be seen from the color bar in the right of these figures. In Figure 3(a) deep yellow represents the period-1 attractor, which is the stability region of the Nash equilibrium point. Furthermore, the deep yellow region in Figure 3(a) also corresponds to the stability region in Figure 1. The dark green color represents the region of period-2 attractor. And the yellow color represents the region of period-4 attractor, and so on. From Figures 3(a) and 3(b), we can see that the system will lose stability through period doubling bifurcation as one/both of the parameters αi,i=1,2, increase. In particular, in the black region the trajectories either converge toward a large period attractor, a chaotic or quasi-periodic one or they diverge. It is worth pointing out that the quasi-periodic attractor is a special attractor, which will arise after the Neimark-Sacker bifurcation. And the largest Lyapunov exponent equals zeros, if the system is in a quasi-periodic state.Figure 3
(a) The two-parameter bifurcation diagram in the parameter planeα1,α2 with fixed parameters as a=7.12, b=0.024, c=0.98, β=0.38, and γ=0.44. (b) Partial enlargement of (a). (c) The two-parameter bifurcation diagram in the parameter plane γ,α2 with fixed parameters as a=2.679, b=0.061, c=0.105, α1=0.666, and β=0.375. (d) Partial enlargement of (c).
(a) (b) (c) (d)Figure3(b) is a partial enlargement of Figure 3(a). From Figure 3(b), we can find some fractal structures. And symmetry of these figures can also be detected. It also means that the symmetry of parameters α1 and α2 also exists in system (17). Therefore, only the changes of dynamical behavior with parameter α1 are studied in this article. There is no need for us to study the effects of parameter α2 on system (17) due to the symmetry.If we fix the value ofα2 as α2=0.4 in Figure 3(a), the detailed dynamical behavior of system (17) can be found in Figures 4(a) and 4(b). In Figures 4(a) and 4(b), the one-parameter bifurcation diagram and the corresponding largest Lyapunov exponent are given, where the bifurcation parameter is chosen as the speed of adjustment of firm 1. And the parameters chosen in Figure 4 are same with those of Figure 3(a). From Figure 4(a), we can see that the Nash equilibrium point is stable when α1∈0,0.41. At about α1=0.41, the Nash equilibrium point loses its stability through a period doubling bifurcation, and a stable state of period-2 appears. As the value of α1 further increases, the stable period-2 orbit loses its stability again, and a period-4 cycle arises, and so on. At last, the system enters into chaos. However, some periodic windows can also be found in the chaotic region, which can be seen from Figure 4(a). Through the above analysis, a common conclusion can be gained that the speed of adjustment plays a vital role in system (17). In other words, smaller speed of adjustment means more stable R&D investments and more stable market environment. In fact, the R&D investments for firms can exhibit a more stable behavior if the reaction speeds for them to the market information are slow enough. Another useful tool to analyze the complex dynamical behavior is the largest Lyapunov exponent. The largest Lyapunov exponent diagram can be used to show the degree of separation of nonlinear system variables running over time, and it can also be used to prove whether the system is chaotic. When the system is in a stable state, the largest Lyapunov exponent is less than zero. While the system enters chaos, its largest Lyapunov exponent is greater than zero. In Figure 4(b), the largest Lyapunov exponent diagram corresponding to Figure 4(a) is plotted. The evolution process of system (17) is confirmed again.Figure 4
(a) The one-parameter bifurcation diagram witha=7.12, b=0.024, c=0.98, β=0.38, γ=0.44, and α2=0.4, where the bifurcation parameter is chosen as α1. (b) The corresponding largest Lyapunov exponent of (a). (c) The one-parameter bifurcation diagram with a=2.679, b=0.061, c=0.105, α1=0.666, and α2=0.9, where the bifurcation parameter is chosen as γ. (d) The corresponding largest Lyapunov exponent of (c).
(a) (b) (c) (d)Moreover, the role of parameterγ is also very important for the system (17). If we fix the parameter sets as a=2.679, b=0.061, c=0.105, α1=0.666, and β=0.375, the two-parameter bifurcation diagrams with parameters α2 and γ are shown in Figures 3(c) and 3(d). In Figures 3(c) and 3(d), the initial conditions are chosen as x1=6.15 and x2=5.75. From Figure 3(c), we can see that the dynamical behaviors of system (17) are very complex, when the value of parameter α2 meets α2>1.4 (to be honest, this value is not very accurate), while, if α2<1.4, we can see that the system is in a chaotic state when the value of parameter γ is relatively small. And then the stable periodic states will arise. However, the system will lose stability again, when the value of γ is large enough. And the escaping trajectories will appear. Figure 3(d) is the partial enlargement of Figure 3(c). Although Figure 3(c) looks ugly, Figure 3(d) does show a good fractal structure. Furthermore, if we fix the value of α2 as α2=0.9 in Figures 3(c) and 3(d), the one-parameter bifurcation diagram is given in Figure 4(c). From Figure 4(c), we can see that the system is in a chaotic state when γ<0.19225. And then a multi-periodic window appears. The chaotic window reappears, subsequently. But when the value of γ meets γ≈0.1994, the system enters a multi-periodic state again. And then a secondary Neimark-Sacker bifurcation occurs at γ≈0.211. The Nash equilibrium point becomes stable at γ≈0.3. Finally, the Nash equilibrium point loses its stability again at γ≈0.58, and the escaping trajectories arise. The escaping trajectory cannot be detected, as the value(s) of x1 or/and x2 will go to infinity. Such changes can also be seen from the corresponding largest Lyapunov exponent diagram, which is given in Figure 4(d).In order to further explain the influence of parameterγ on system (17), the two-dimensional bifurcation diagram with different parameter sets is given in Figure 5. The parameters are chosen as a=7.12, b=0.024, c=0.98, and β=0.38, which are same with those of Figure 3(a). But the value of parameter γ is different from Figure 3(a). In Figures 5(a) and 5(b), the value of γ is chosen as γ=0.1, and Figure 5(b) is the partial enlargement of Figure 5(a). However, the value of γ is chosen as γ=1.0 in Figures 5(c) and 5(d), and Figure 5(d) is the corresponding partial enlargement of Figure 5(c). By comparing Figures 3(a) and 5, we can see that the Nash equilibrium point always loses its stability through a period doubling bifurcation, and then a period-2 cycle appears. After that the system enters to chaotic state through a series of period doubling bifurcation. The detailed bifurcation process can be found in Figure 6(a), where the parameters are same with Figure 5(a) except that the value of α2 is fixed as α2=0.345. From Figure 6(a), a clear secondary Neimark-Sacker bifurcation can be found. In addition, we can also find the existence of periodic windows on the way to chaos.Figure 5
The two-parameter bifurcation diagram in the parameter planeα1,α2 with fixed parameters as a=7.12, b=0.024, c=0.98, and β=0.38. (a) The two-parameter bifurcation diagram with γ=0.1. (b) Partial enlargement of (a). (c) The two-parameter bifurcation diagram with γ=1.0. (d) Partial enlargement of (c).
(a) (b) (c) (d)Figure 6
The one-parameter bifurcation diagrams, corresponding equilibrium profits, and average profits of firm 1 and firm 2 with fixed parameters asa=7.12, b=0.024, c=0.98, β=0.38, and α2=0.345, where the blue color corresponds to firm 1, while the red color represents that of firm 2. (a) The one-parameter bifurcation diagram when γ=0.1. (b) The equilibrium profits and average profits corresponding to (a). (c) The one-parameter bifurcation diagram when γ=0.44. (d) The equilibrium profits and average profits corresponding to (c). (e) The one-parameter bifurcation diagram when γ=1.0. (f) The equilibrium profits and average profits corresponding to (e).
(a) (b) (c) (d) (e) (f)From Figures3(a) and 5, we can also see that the parameter window into quasi-periodic state becomes very narrow when γ increases to some certain value. When γ=1.0, the parameter window of Neimark-Sacker bifurcation is even hard to observe; see Figure 5(c). But this phenomenon can be revealed from the partial enlargement of Figure 5(c); see Figure 5(d).The symmetry structure and fractal structure can also be observed from Figures3(a) and 5. As the value of γ increases, the fractal structure of these figures has also changed. But the symmetry structure of all the figures will never change. This phenomenon also directly reflects the symmetry of parameters α1 and α2. In addition, we can also find that there are a lot of sparse black points, which have covered the light blue area and the yellow area in the upper right corner of Figure 5(d). Here, these sparse black points are called noisy points, as they have affected the beauty of Figure 5(d). However, the presence of the noisy points often implies the coexistence of multiple attractors. And the problem of coexistence of attractors will be discussed in the following section.The parameterγ plays a very important role in the analysis of high-tech enterprises and it can reflect the R&D capability of enterprises. The smaller value of parameter γ means the stronger innovation ability of enterprises. The following analysis is mainly focused on the influence of parameter γ on the dynamical behaviors and profits of firms in the given model. To be specific, the equilibrium profits and the average profits of firm 1 and firm 2 will be discussed in order to enclose the relation between parameter γ and the behaviors of firms. The equilibrium profits of firm i can be gained by substituting the value of E3 into (11) and (12), respectively. The specific mathematical expressions of equilibrium profits will be given as(36)π1Eq=π2Eq=9γa-2bc2-8γba-2bcβ+12b9γ-8β+12b2And the average profits of firm 1 and firm 2 can be calculated through the mean value of the profit of firmi at different period. The detailed expressions of average profits of firm 1 and firm 2 can be written as (37)π1Av=1n∑t=1na-2bc+bβ+1x1t+x2t9b2-12γx12tπ2Av=1n∑t=1na-2bc+bβ+1x1t+x2t9b2-12γx22twhere n is the number of iteration. In the following numerical simulation, the value of n is chosen as n=5000. Figure 6 gives the one-parameter bifurcation diagrams, the corresponding equilibrium profits, and average profits of firm 1 and firm 2, where the blue color corresponds to firm 1, while the red color represents firm 2. In Figure 6, the fixed parameters are given as a=7.12, b=0.024, c=0.98, β=0.38, and α2=0.345. However, the value of γ is fixed as γ=0.1, γ=0.44, and γ=1.0, respectively. From the value of E3, we can see that the R&D efforts of firm 1 and firm 2 are equal at the Nash equilibrium point. And thus the equilibrium profits are same as well for firm 1 and firm 2. This phenomenon can also be seen from Figure 6. However, the R&D efforts and average profits of firm 1 and firm 2 are no longer equal any more, once the Nash equilibrium point loses its stability, as α1 increases. Through comparing Figures 6(b), 6(d), and 6(f), we can find that the average profits will fall when the value of γ increases. It is easy to understand that larger value of γ represents a lower R&D capability, so the profits will fall with increasing of γ. And the amplitude of fluctuation of R&D efforts is smaller for firm with slower speed of adjustment; see Figures 6(a), 6(c), and 6(e).From Figure6, we can see that the average profits would be smaller than the equilibrium profits when the Nash equilibrium point loses stability. However, the average profits are close to the equilibrium profits, if the firm’s speed of adjustment is large enough; see Figures 6(d) and 6(f). It is worth pointing out that the average profit is even higher than the equilibrium profit for the firm with large speed of adjustment, once the system enters into chaos; see the blue curves in Figures 6(d) and 6(f), while the average profit of another firm fluctuates very fierce and is far below its equilibrium profit; see the red curves in Figures 6(d) and 6(f). This phenomenon can be concluded as chaos is usually harmful to the entire market, but firms with large speed of adjustment sometimes can even benefit from chaos.


## 5. Multistability

As we have noticed, there are some noisy points in the upper right corner of Figure5(d). These noisy points are probably caused by the coexistence of multiple attractors as we have mentioned above. To give a more explicit picture about this problem, the one-parameter bifurcation diagrams with different initial conditions are plotted if we fix the parameters as a=7.12, b=0.024, c=0.98, β=0.38, and γ=0.44 and choose the parameter α1=α2=α as the bifurcation parameter. The one-parameter bifurcation diagrams with initial conditions x1,x2=2.33,1.3 and x1,x2=2.33,2.3 are given in Figures 7(a) and 7(b), respectively. The differences between these two figures are very clear. In order to highlight these differences, the bifurcation diagrams are plotted with different colors. In Figure 7(a), a period doubling bifurcation firstly happens at α≈0.4175, and then a Neimark-Sacker bifurcation arises at α≈0.51175. Later, a period-2 point appears again. After that the system enters into chaos through a period doubling bifurcation. In Figure 7(b), the dynamical behavior of system (17) is relatively simple, but a transcritical bifurcation can be observed.Figure 7
The one-parameter bifurcation diagram witha=7.12, b=0.024, c=0.98, β=0.38, γ=0.44, and α1=α2=α. (a) The one-parameter bifurcation diagram in relation to α with initial conditions x1,x2=2.33,1.3; (b) The one-parameter bifurcation diagram in relation to α with initial conditions x1,x2=2.33,2.3.
(a) (b)In Figure7, different colors represent different meanings. From left to right, there is no coexistence attractor in the red region. Blue corresponds to the coexistence of attractors. Three different forms of coexistence of multiple attractors, which are the coexistence of two period-2 attractors, the coexistence of a period-2 attractor with a 2-cyclic attractor, and the coexistence of a period-2 attractor with a chaotic attractor, can be found in the blue region. However, the coexistence of multiple attractors will disappear again in the green area.The multi-stability in economic system has a very important research meaning. In fact, the choices of initial conditions are very crucial for firms in the real market. As the system we built here can be regarded as a repeated process, the state of firms at timet can be considered as the initial conditions for firms at timet+1. Furthermore, once the initial conditions are given for system (17), the long and asymptotic behaviors of this system after many iterations can be predicted. All the time, a recognized tool to analyze the coexistence of attractors is the basin of attraction. And thus, the basin of attraction is employed in our article. Firstly, the definition of basin of attraction is given here.Suppose there is an attractor of a nonlinear system, which can be expressed byA. Then all the basins of attractor A can be written as BA=⋃n=0∞T-nUA, where UA is the immediate basin of A. Obviously, UA⊆BA, and the points mapped to U(A) after a finite number of iterations also belong to BA. T-nx represents the set ofrank-n preimages of x. The basin of attraction can be divided into the attractor, domain of attraction, and domain of escape. And the domain of attraction refers to the set of initial conditions that converge to the same attractor after a series of iterations. The domain of attraction can be divided into connected sets and nonconnected sets. The domain of escape refers to the set of initial conditions that diverge to infinity after a series of iterations. If we fix the parameters as a=7.12, b=0.024, c=0.98, β=0.38, γ=0.44, α1=0.41844, and α2=0.4184, a typical basin of attraction of system (17) is given in Figure 8. In Figure 8, the domain of attraction and the escaping domain are represented by yellow area and dark blue area, respectively.Figure 8
Basin of attraction of system (17). The yellow region is the basin of attraction of the Nash equilibrium point E3, and the dark blue region represents domain of escaping. O-11, O-12, and O-13 are therank-1 preimages of the origin O=0,0; ω1-1 and ω2-1 are therank-1 preimages of ω1 and ω2. The boundaries of attracting domain of E3 are constituted by these lines.As it can be seen from Figure8, the shape of the basin of Nash equilibrium point is a quadrilateral OO-11O-13O-12, where O=0,0 is the origin of coordinate axis and the rest vertices of this quadrilateral O-11, O-12, and O-13 are therank-1 preimages of the origin O=0,0. It is clear that the points of O-11 and O-12 are located on the horizontal axis x1 and the vertical axis x2, respectively. And their coordinates can be gained through solving the algebraic equation xi′=0 in formula (20), which are given as (38)O-11=9+4Aα19γ-4Bbα1,0,O-12=0,9+4Aα29γ-4Bbα2For the convenience of analysis, the line-segmentsOO-11 and OO-12 are labeled as ω1 and ω2, i.e., OO-1i=ωi,i=1,2. And therank-1 preimages of ωi,i=1,2, are labeled as ω1-1 and ω2-1, respectively. Therank-1 preimages of ωi,i=1,2, can be understood as a set of all the preimages of the points on ωi,i=1,2. In order to get the coordinates of points located on ω1-1, it is supposed that the coordinate of any point on ω1 can be represented as x,0; then the preimages of x,0 should satisfy the following equations:(39)x11+α14b9β+12x1+x2+49a-2bcβ+1-γx1=xx21+α24b9β+12x1+x2+49a-2bcβ+1-γx2=0From the second equation of (39) we can find that therank-1 preimage of ω1 should meet x2=0 or the linear function given below:(40)x2=4Bb9γ-4Bbx1+9+4Aα29γ-4Bbα2The function x2=0 represents that ω1 is one of therank-1 preimages of itself. However, the linear function (40) gives anotherrank-1 preimage of ω1, that is, ω1-1. Through similar method, the coordinates of points located on ω2-1 also can be solved. And the following linear function can determine the straight line ω2-1:(41)x1=4Bb9γ-4Bbx2+9+4Aα19γ-4Bbα1It is obvious that the line determined by (40) and the line determined by (41) have a unique intersection O-13, which is one of therank-1 preimages of the origin O=0,0. The coordinate of point O-13 can be gained through solving the combination of (40) and (41), which can be represented as(42)O-13=4Bbα1-α2+9γα2+4Aγα1α29γ-8Bbγα1α2,4Bbα2-α1+9γα1+4Aγα1α29γ-8Bbγα1α2A more intuitive explanation of the pointsO, O-11, O-12, O-13 and these straight lines ω1, ω2, ω1-1 and ω2-1 can be found in Figure 8. It is clear that the boundaries of the attracting domain of the Nash equilibrium point E3 are constituted by these lines. And the points O, O-11, O-12, and O-13 are just the intersections of these straight lines. From the detailed expressions of (40) and (41), we find that the structure of domain of attraction is related to the system parameters. Any change in these parameters will lead to changes of structure of basins.For example, we fix all the parameters asa=7.12, b=0.024, c=0.98, β=0.38, and γ=0.44, which are same with the parameters given in Figure 7. And the parameter α=α1=α2∈0.53325,0.60275 varies. Firstly, the parameter α is fixed as 0.53325, and the corresponding basins can be found in Figure 9(a). At this very moment, a chaotic attractor coexists with a period-2 cycle. The attracting domain of the chaotic attractor is represented by yellow color, and the attracting domain of the period-2 cycle is represented by light green, while the dark blue area means the domain of escaping. And the shape of Figure 9(a) is just like a quadrilateral “sieve” with fractal structure. As parameter α further increases to 0.55675, the chaotic attractor enlarges, but the fractal structure of the attracting basin has hardly changed. The chaotic attractor has almost contacted with its attracting domain; see Figure 9(b). It means that a global bifurcation called “contact bifurcation” will arise. In Figure 9(c), the parameter α is fixed as 0.557, we can find that both the chaotic attractor and its corresponding domain of attraction have been destroyed, and the destroyed attractor will fill the whole domain of attraction. This phenomenon can also be observed from the time series plot, which is given in Figure 10. In Figure 10, the values of parameters are same with those in Figure 9(c) and the initial conditions are chosen as 2.33,1.3. We can find that the evolution is disorderly in Figure 10. Comparing Figures 9(b) and 9(c), it can be found that the truth of the matter is that the attracting basin of the “dead” chaotic attractor is also occupied by the period-2 attractor. We refer to the long-run transition as the “ghost” of the destroyed attractor; see Bischi, etc. [17, 39] for more detail. Furthermore, if we choose α as 0.58117, it can be found that the “ghost” of the destroyed attractor still exists, and the coexistence of a period-4 point (red points in Figure 9(d)) and an 8-piece quasi-periodic attractor (black cycle in Figure 9(d)) can be found. But the attracting domain of the period-4 points is confined to the diagonal. As α further increases to 0.58225, the 8-piece quasi-periodic attractor has evolved into a 4-piece quasi-periodic attractor. The coexisting period-4 attractor is hard to distinguish from the quasi-periodic attractor; see Figure 9(e). As α continues to increase, the 4-piece quasi-periodic attractor becomes a 4-piece chaotic attractor. And some disconnected sets arise in the attracting basins of these two attractors; see Figure 9(f). The 4-piece chaotic attractor starts to contact with its attracting domain. This means that global bifurcation occurs, and the structure of attracting domain will change once again; see Figure 9(g). In Figure 9(g), we fix the value of α as 0.6025, and the coexistence of a 4-piece chaotic attractor and a multi-periodic cycle in the diagonal can be found. The attracting basin of the multi-periodic cycle is only the diagonal line. With a further increase of the parameter α, the 4-piece chaotic attractor starts to merge into a one-piece chaotic attractor at α=0.60275; see Figure 9(h). However, a weak chaotic attractor is beginning to appear on the diagonal, too. And its attracting domain is just the diagonal line.Figure 9
Basins of attraction and coexisting attractors of system (17) with a=7.12, b=0.024, c=0.98, β=0.38, γ=0.44, and α=α1=α2∈0.53325,0.60275, where different colors represent attracting domains of different attractors. (a) α=0.53325; (b) α=0.55675; (c) α=0.557; (d) α=0.58117; (e) α=0.58225; (f) α=0.58525; (g) α=0.6025; (h) α=0.60275.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 10
The time series plot with parametersa=7.12, b=0.024, c=0.98, β=0.38, γ=0.44, and α1=α2=0.557.In Figure9, the multistability of system (17) with α1=α2 has been discussed. Then, the cases with α1≠α2 will be discussed. Firstly, the parameter sets are fixed as a=4.19, b=0.04, c=0.13, β=0.32, and γ=0.098, and the parameters α1 and α2 are chosen as the bifurcation parameters. In Figure 11(a), a two-parameter bifurcation diagram is plotted. And the local enlargement of Figure 11(a) is given in Figure 11(b). The noisy points can be found in the lower right corner of Figure 11(b). In the previous section, we have mentioned that these noisy points are related with the coexistence of attractors. To throw further light on this problem, the values of α1 and α2 are chosen as 0.524 and 0.2037, respectively. This set of parameter corresponds to the overlapping region of fluorescein and turquoise in Figure 11(b). The basins of attraction at this set of parameters are given in Figure 11(c). From Figure 11(c), we can find that there are two attractors coexisting: one is a period-2 cycle (labeled by A2), and the other is a period-6 cycle (labeled by B6). The immediate basin of A2 is the area where A2 locates, while the total basin of A2, which can be written as BA2=⋃n=0∞UA2, is the entire region with light blue. In fact, the light blue region occupies the vast majority of the feasible region; see Figure 11(c). Accordingly, the total basin of B6 is very sparse, which can be represented as BB6=⋃n=0∞UB6 and is expressed with yellow color in Figure 11(c). The feasible region can be represented as F=BA2∪BB6=⋃n=0∞UA2∪⋃n=0∞UB6.Figure 11
(a) The two-parameter bifurcation diagram in the parameter planeα1,α2 with fixed parameters as a=4.19, b=0.04, c=0.13, β=0.32, and γ=0.098. (b) Partial enlargement of (a). (c) The basin of attraction, where a period-2 cycle coexists with a period-6 cycle when α1=0.524 and α2=0.2037. (d) The basin of attraction, where a period-4 cycle coexists with a 7-piece chaotic attractor when α1=0.5665 and α2=0.2037.
(a) (b) (c) (d)If the values ofα1 and α2 are changed to 0.5665 and 0.2037, respectively, other coexisting attractors can be found, see Figure 11(d), where a period-4 cycle (labeled by A4) coexists with a chaotic attractor constituted by 7-pieces (labeled by B7). In Figure 11(d), the total basin of A4 is the area of light blue, and the total basin of B7 is the area of yellow. In Figure 11(d), the yellow region is further increased, while the light blue region is reduced compared with Figure 11(c). The feasible regions F are the total basins of these two attractors. The shape of Figure 11(d) is similar to that of Figure 11(c), and it is just like a knife coin in ancient China.It is generally known that the critical set is the primary tool to analyze the topological structure of attracting basin in noninvertible maps. If we letxit+1=xi′, xit=xi in system (17), and solve the algebraic equations about x1 and x2, which are given as (43)x11+α14b9β+12x1+x2+49a-2bcβ+1-γx1=x1′x21+α24b9β+12x1+x2+49a-2bcβ+1-γx2=x2′Then, therank-1 preimages of point x1′,x2′ can be obtained. Actually, the preimages of point x1′,x2′ may be more than one, as the nonlinear equations (43) are constituted by two quadratic equations. And therefore, the system can be regarded as a two-dimensional noninvertible map, which can be abbreviated as x1′,x2′=Tx1,x2. The noninvertible map represents many-to-one. It means that the image of two distinct points on the two-dimensional plane may be same. That is, even if x1,x2≠y1,y2, Tx1,x2 may be equal to Ty1,y2, too. In fact, if we solve the nonlinear algebraic equations (43), there are maybe four, maybe two, or even no solutions for equations (43) in the real field. And the critical curves can be used to divide the phase space into the regions called Z4, Z2, and Z0. The symbol Zi represents the region in phase space, where the points have i preimages. The boundary of these regions, also called critical curves, can be denoted as LC. Suppose that LC-1 are the preimages sets of points on LC, that is, LC=TLC-1. And LC-1 can be defined as a locus of the points where the determination of Jacobian matrix vanishes, that is, (44)LC-1⊆J0=x1,x2∈R2∣det⁡Jx1,x2=0As the Jacobian matrix of system (17) has given in (23), thus the points on LC-1 can be determined by the following equation:(45)Ex12+x22+Fx1x2+Gx1+x2+H1x1+H2x2+I=0where E=8Bb4Bb-9γα1α2/81, F=44Bb-9γ2α1α2/81, G=4A12Bb-18γα1α2/81, H1=24Bb-9γα1+4Bbα2/9, H2=24Bb-9γα2+4Bbα1/9, and I=9+4Aα19+4Aα2/81.It can be found that there are two branches inLC-1, i.e., LC-1(a) and LC-1(b). Therefore, there are also two branches in LC, that is, LCa=TLC-1a and LCb=TLC-1b. A more intuitive image of LC and LC-1 can be found in Figure 12. The sketch of LC-1 has been given in Figure 12(a), while the sketch of LC has been given in Figure 12(b). In Figure 12(b), we find that the branch LCa separates Z2 from Z0, while the other branch LCb separates Z2 from Z4. For example, the origin point O=0,0 has four preimages as O-1i,i=1,2,3 and itself, so the origin O∈Z4.Figure 12
(a) The sketch ofLC-1 with a=1.95, b=0.044, c=1.58, β=0.9, γ=0.77, α1=1.57, and α2=1.625. (b) Critical curves LC=TLC-1. The curve LCa separates Z2 from Z0, and the other curve LCb separates Z2 from Z4.
(a) (b)In the above discussion, the coexistence of attractors in system (17) has been shown. However, another interesting phenomenon in the basin of attraction called “hole” by Bischi, etc., [40] can be found in system (17). The formation mechanism of these holes and their evolutions will be discussed in the following through the critical curves, which we have given above.First of all, the two-parameter bifurcation diagram is given in Figure13, where the parameters are fixed as a=1.95, b=0.044, c=1.58, β=0.9, and γ=0.77. The specific meaning of the color in Figure 13 can be referred to the color bar in the right side of these figures. Figure 13(b) is the local enlargement of Figure 13(a). The obvious fractal structure can be found in these two figures. But there are also a lot of noisy points in these figures, which have already affected the beauty of these figures. Similarly, these noisy points also give us the information of the coexistence of multiple attractors.Figure 13
(a) The two-parameter bifurcation diagram in the parameter planeα1,α2 with fixed parameters as a=1.95, b=0.044, c=1.58, β=0.9, and γ=0.77. (b) Partial enlargement of (a).
(a) (b)As a matter of fact, if we fixα1 and α2 as 1.572 and 1.4893, the basin of attraction is given in Figure 14(a). The rough boundaries of feasible region can be determined using Eqs. (40) and (41), but there is also a small gap in the feasible region; see the topmost of Figure 14(a) for detail. And the boundary of basin is tangent to the critical curve LCa at the gap. No coexisting of attractors can be found at such a situation. Only one chaotic attractor exists in the feasible region. That is, the feasible region is just the attracting basin of that chaotic attractor. The critical curve LCa is also the boundary of that chaotic attractor. It means that the chaotic attractor has contacted with the boundary of its attracting domain; then a contact bifurcation will happen soon.Figure 14
The corresponding basin of attraction for different values ofα2, and the parameters are fixed as a=1.95, b=0.044, c=1.58, β=0.9, γ=0.77, and α1=1.572. (a) α2=1.4893: a chaotic attractor and its attracting domain. (b) When the parameter α2 is equal to 1.49, a few holes start to arise in the feasible region. (c) At α2=1.52, more holes arise. (d) At α2=1.62, the fractal structure of the attracting domain is very obvious, and many small holes are interconnected to larger holes. (e) A chaotic attractor and its attracting domain when we fix parameter α2 as 1.679, and the chaotic attractor starts to contact with the boundary of its basin; a global bifurcation will happen. (f) At α2=1.6815, the chaotic attractor disappears and the ghost of its basin is still surviving.
(a) (b) (c) (d) (e) (f)Whenα2 increases to 1.49, the contact bifurcation happens and the chaotic attractor disappears. And the feasible region is full of the “ghost” of the disappeared attractor; see Figure 14(b). At the same time, the critical curve LCa starts intersecting with the gap, and a very tiny hole H0 in the infeasible region arises. As the critical curve LCa separates the phase space into Z0 and Z2, so the closed hole H0 should have two preimages. And the two preimages of H0, which are represented as H-11 and H-12, form another two holes in the feasible region; see Figure 14(b) for detail. The holes H-11 and H-12 lie in the Z4 region, which is different from the research result of other scholars [40]. Therefore, there should be four preimages for the hole H-11 and H-12, respectively. One of the preimages of hole H-11 and H-12 is given as H-21 and H-22, respectively. And other preimages are too small to see.Asα2 further increases to 1.52, the main hole H0 becomes larger than that of Figure 14(b), and more holes arise. These holes can be regarded as therank-k(k=0,1,2,…) preimages of H0; see Figure 14(c). By comparison with Figure 14(b), the “ghost” of chaotic attractor has disappeared and a new period-2 cycle emerges in Figure 14(c). With a further increase of α2 to 1.62, the main hole H0 is getting larger, and the preimages of H0 are going to get bigger, too. In the meanwhile, the attractor has evolved from a period-2 cycle to a period-4 cycle; see Figure 14(d). It is very clear that the holes in the feasible region will grow bigger and bigger and burst into larger holes, as the main hole H0 is no longer closed; see Figure 14(e). In Figure 14(e), we can see that the number of holes has been reduced, but the scale of these holes has increased. Furthermore, a chaotic attractor arises when α2=1.679. However, the newborn chaotic attractor is going to die, as there is no place for it to live. And then a contact bifurcation happens at α2=1.6815, the newborn chaotic attractor has already died, and only its ghost still remains in the feasible region; see Figure 14(f).


## 6. Conclusion

In this paper, a two-stage dynamical duopoly game of R&D competition between two high-tech enterprises is studied. At the first stage, both firms that conduct R&D so as to reduce their production cost and improve the quality of products compete in the R&D level. And at the second stage, both firms choose cooperation and pursue the maximization of joint profit. And then, the existence and stability of all the equilibrium points are discussed. It has been proved that the boundary equilibrium points are always unstable, and the Nash equilibrium point is stable only when the parameters meet some conditions. The stability region of Nash equilibrium point is derived through Jury condition. We find that the Neimark-Sacker bifurcation cannot occur as both eigenvalues of the Jacobian matrix at Nash equilibrium point are always real numbers. And a numerical simulation is employed in order to check the theoretical results. The stability region is plotted numerically with different values of parameterγ. We find that the stability region is easy to be influenced by the parameter γ. The stability region is bigger when the value of γ is larger. But the stability region will not change when the parameter γ is large enough. We also find that the parameter window into quasi-periodic state becomes very narrow when γ increases to some certain value. In addition, the influence of adjusting speed on the system is also discussed in this paper.Another topic of this paper is the multi-stability. The multistability always implies path-dependence, which means the long run behavior of firms is deeply influenced by historical accidents. That is, a small perturbation of initial conditions will affect the evolution of the system. And hence, the choice of initial conditions plays a very important role in the long run of the built system. For the sake of analyzing the effect of initial conditions on the long run evolution of the system, the basin of attraction is employed. In this research, the two-parameter bifurcation diagram is used to find the coexistence of multiple attractors. We found that there exist noisy points in the two-parameter bifurcation diagram, when the phenomenon of coexistence of multiple attractors arises. And then the coexistence of multi-attractors is studied with two different cases:α1=α2 and α1≠α2, respectively. Our research result shows that the topological structures of these two cases are totally different. The global bifurcation, such as “contact” bifurcation, and the phenomenon of “holes” in attracting basins have been found through numerical simulation. And at last, the formation mechanisms of these phenomena have been analyzed through critical curves and noninvertible maps.



---

*Source: 2369898-2019-08-07.xml*