Hopf Bifurcation Analysis in a Modified Price Differential Equation Model with Two Delays

Abstract

The paper investigates the behavior of price differential equation model based on economic theory with two delays. The primary aim of this thesis is to provide a research method to explore the undeveloped areas of the price model with two delays. Firstly, we modify the traditional price model by considering demand function as a downward opening quadratic function, and supply and demand functions both depending on the price of the past and the present. Then the price model with two delays is established. Secondly, by considering the price model with one delay, we get the stable interval. Regarding another delay as a parameter, we studied the linear stability and local Hopf bifurcation. In addition, we pay attention to the direction and stability of the bifurcating periodic solutions which are derived by using the normal form theory and center manifold method. Afterwards, the study turns to simulate the results through numerical analysis, which shows that the provided method is valid.

1. Introduction

Recent developments in mathematical economics and in problems of business administration have led to extensive use of differential equation model. Bifurcations and chaos always show in the contemporary literature of economics as basic concepts. The paper of Shuhe [1] was the pioneering work in studying price differential equation model which provides a dynamic system to investigate sectoral dynamics of an economy phenomenon. Since demand depends on the past price, further study on price model with delay can be found in [2]. Reference [3] focused on the phenomenon of bifurcation which forms an integral part of qualitative approach to study dynamical systems. Similarly, [4–6], which investigated the local Hopf bifurcation and the existence of periodic solutions of price model, had important consequences for theoretical and empirical model building in economics. Reference [7] provided a brief survey of the literature of bifurcation detections in economic models. The existence of different bifurcate parameters leads to chaos and the causes of complicated phenomenon were argued in [8]. In short, few people studied price model with delay; what is more, no results involving price model with two delays have occurred.

In order to illustrate the economic phenomena with price varying accurately, a reasonable mathematical model of price is needed. Thus, we introduce a traditional price differential equation model in [1]: 𝑑²𝑃(𝑡)𝑑𝑡²(𝑃(𝑡(𝑡)=𝜇𝛿))𝑑𝑃𝑑𝑡−𝜇𝑏₀𝑃(𝑡)+𝜇𝑑₀−𝜇𝑔₀. (1) The meaning of parameters refers to [1]. We modify the traditional price differential equation model by considering the following factors.

Firstly, we denote the correlation coefficient between demand and price rising rate by 𝛿(𝑃(𝑡)). We consider demand function as a downward opening quadratic function and supply function as a linear function: 𝛿(𝑃(𝑡))=𝑏𝑃(𝑡)+𝛽²+𝐶₀,𝑏<0,𝛽<0,𝐶₀<0. (2)

Secondly, according to the cobweb theory, since manufacturers need a production cycle time from obtaining market information to adjust the production line, the role of price adjustment lags on the supply function of time. The purchase of the consumers also depends on the price change and decision whether to buy or not. Then, supply and demand functions both depend on the price of the past and the present. We introduce two delays 𝜏₁, 𝜏₂ to denote supply and demand functions, which depend on the price of the past, respectively. A differential system with two delays for price differential equation model is transformed into the following form: 𝑑²𝑃(𝑡)𝑑𝑡²=𝜇𝑏𝑃+𝛽²+𝐶₀𝑑𝑃𝑡−𝜏₂𝑑𝑡−𝜇𝑏₀𝑃(𝑡)−𝜇𝑎𝑃𝑡−𝜏₁+𝜇𝑑₀−𝑔₀. (3)

Let 𝑃(𝑡)=𝑥(𝑡), 𝑑𝑃(𝑡)/𝑑𝑡=𝑦(𝑡), and notice that supply and demand functions both depend on the price of the past and the present; then the model is described by the following autonomous system: 𝑥̇(𝑡)=𝑑𝑦(𝑡)+𝑘𝑦𝑡−𝜏₂,𝑦(̇𝑡)=𝜇𝑏𝑥(𝑡)+𝛽²+𝐶₀𝑦𝑡−𝜏₂−𝜇𝑏₀𝑥(𝑡)−𝜇𝑎𝑥𝑡−𝜏₁+𝜇𝑑₀−𝑔₀, (4) where 𝑑, 𝑘 refers to [2]. 𝑥(𝑡) is the price at time 𝑡, and 𝑦(𝑡) is the amount of supply at time 𝑡; 𝜇>0, 𝑏₀>0, 𝑎>0, 𝑏<0, 𝛽<0, 𝐶₀<0.

Different from the previous work in [1, 2, 4], the purpose of this paper is to investigate the stability of the Hopf bifurcation and the direction of bifurcation periodic solution of a price differential equation model with two delays. The structure of the paper is as follows. In Section 2, linear stability and local Hopf bifurcations are studied by using qualitative methods. In Section 3, we regard 𝜏₂ as bifurcation parameter and consider (4) with 𝜏₁ in its stable interval; the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions are derived by using the normal form theory and center manifold method. Afterwards, the presented numerical simulations in Section 4 have demonstrated the theoretical analysis.

2. Local Stability and Hopf Bifurcation

Obviously, system (4) always has an equilibrium 𝐸^∗=(𝑥^∗,𝑦^∗)=((𝑑₀−𝑔₀)/(𝑏₀+𝑎),0). In the following, we will investigate the effect of the delay 𝜏₁, 𝜏₂ on the dynamics of system (4).

Let 𝑢(𝑡)=𝑥(𝑡)−𝑥^∗, 𝑣(𝑡)=𝑦(𝑡)−𝑦^∗; the linearization of system (4) at zero steady state is 𝑢̇(𝑡)=𝑑𝑣(𝑡)+𝑘𝑣𝑡−𝜏₂,𝑣(̇𝑡)=𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑣𝑡−𝜏₂−𝜇𝑏₀𝑢(𝑡)−𝜇𝑎𝑢𝑡−𝜏₁. (5)

For convenience, as in literature of [5], we denote 𝐴=−𝜇[𝑏(𝑃₀+𝛽)²+𝐶₀)]>0, 𝐵=𝜇𝑏₀>0, 𝐶=𝜇𝑎>0. Then (5) become 𝑢̇(𝑡)=𝑑𝑣(𝑡)+𝑘𝑣𝑡−𝜏₂,𝑣(̇𝑡)=−𝐴𝑣𝑡−𝜏₂−𝐵𝑢(𝑡)−𝐶𝑢𝑡−𝜏₁. (6)

The corresponding characteristic equation is 𝜆²+𝐴𝑒^−𝜆𝜏2𝜆+𝐵𝑘𝑒^−𝜆𝜏2+𝐶𝑑𝑒^−𝜆𝜏1+𝐶𝑘𝑒^{−𝜆(𝜏1+𝜏2)}=0. (7)

To study the stability of equilibrium 𝐸^∗ of (4) and Hopf bifurcation, it is sufficient to analyze the distribution of the roots of (7). It is stable if all roots of (7) have negative real parts and unstable if one root has positive real part. In order to study the characteristic (7) with two delays, we first consider (7) with one delay. Without loss of generality, we choose 𝜏₁ as a parameter and employ Rouche's theorem in Cooke and Grossman [9]; we will find the stable interval for 𝜏₁. Then we consider (7) with 𝜏₁ in its stable intervals. Using Rouche's theorem again, regard 𝜏₂ as a parameter; we will find the stable interval (depends on 𝜏₁) for 𝜏₂. Then we obtain the stable interval for system (4).

Now we analyse the case when 𝜏₂=0; (7) becomes 𝜆²+𝐴𝜆+𝐵𝑘+𝐶𝑑𝑒^−𝜆𝜏1+𝐶𝑘𝑒^−𝜆𝜏1=0. (8) Using a procedure similar to [10], we make some hypotheses as follows: (H1)𝐴²≤4[𝐵𝑘+𝐶𝑑+𝐶𝑘],(H2)𝐴²>4[𝐵𝑘+𝐶𝑑+𝐶𝑘],(H3)𝐵²𝑘²<𝐶²(𝑑+𝑘)²,(H4)𝐵²𝑘²>𝐶²(𝑑+𝑘)², 2𝐵𝑘>𝐴², and (−2𝐵𝑘+𝐴²)²>4[𝐵²𝑘²−𝐶²(𝑑+𝑘)²],(H5)neither (H3) nor (H4).Obviously, when 𝜏₁=𝜏₂=0, (4) becomes a system of ODE. Under the hypothesis (H1), all roots of (7) have negative real parts if and only if 𝐴>0; under the hypothesis (H2), all roots of (8) have negative real parts if and only if 𝐴>0 and 𝐵𝑘+𝐶𝑑+𝐶𝑘>0. Above all, either (H1) or (H2) holds for 𝜏₁=𝜏₂=0, and all roots of (7) have negative real parts.

Applying the lemma in [9] again, we obtain the following results.

Lemma 1. For (8), one has the following: (i)if (H3) holds and 𝜏₁=𝜏⁽¹⁾_1,𝑛, then (8) has a pair of purely imaginary roots ±𝑖𝜔₊;(ii)if (H4) holds and 𝜏₁=𝜏⁽¹⁾_1,𝑛(𝑟𝑒𝑠.𝜏₁=𝜏⁽²⁾_1,𝑛) then (8) has a pair of imaginary roots of ±𝑖𝜔₊(𝑟𝑒𝑠.±𝑖𝜔₋);(iii)if (H5) holds and 𝜏₂>0, then (8) has no purely imaginary root, where 𝜔²_±=2𝐵𝑘−𝐴²2±14−2𝐵𝑘+𝐴²²−𝐵²𝑘²+𝐶²(𝑑+𝑘)^21/2,𝜏⁽¹⁾_1,𝑛=1𝜔₊𝜔arccos²⁺−𝐵𝑘+𝐶𝑑+𝐶𝑘2𝑛𝜋𝜔₊,𝜏⁽²⁾_1,𝑛=1𝜔₋𝜔arccos²⁻−𝐵𝑘+𝐶𝑑+𝐶𝑘2𝑛𝜋𝜔₋,(𝑛=0,1,…). (9)(10)

Proof. Let 𝜆=𝑖𝜔 be the root of (8); we thus have −𝜔²+𝐴𝜔𝑖+𝐵𝑘+𝐶(𝑑+𝑘)𝑒^−𝜆𝜏1=0. (11) Separating the real and imaginary parts −𝜔²+𝐵𝑘+𝐶(𝑑+𝑘)cos𝜔𝜏₁()=0,𝐴𝜔−𝐶𝑑+𝑘sin𝜔𝜏₁=0. (12) Eliminating the harmonic terms gives 𝜔⁴+−2𝐵𝑘+𝐴²𝜔²+𝐵²𝑘²−𝐶²(𝑑+𝑘)²=0. (13) Obviously 𝜔²₁⋅𝜔²₂=𝐵²𝑘²−𝐶²(𝑑+𝑘)²,𝜔²₁+𝜔²₂=2𝐵𝑘−𝐴². (14) From (12) and by direct computation, we obtain 𝜔²_±, 𝜏⁽¹⁾_1,𝑛, 𝜏⁽²⁾_1,𝑛. The result is straightforward.

Denote the minimum value of 𝜏_1,𝑛 by 𝜏⁰₁; that is, min(𝜏_1,𝑛)=𝜏⁰₁, and 𝜆_𝑘,𝑛=𝛼_𝑘,𝑛𝜏₁+𝑖𝜔_𝑘,𝑛𝜏₁,𝛼_1,𝑛𝜏⁽¹⁾_1,𝑛=0,𝜔_1,𝑛𝜏⁽¹⁾_1,𝑛=𝜔₊,𝛼_2,𝑛𝜏⁽²⁾_1,𝑛=0,𝜔_2,𝑛𝜏⁽²⁾_1,𝑛=𝜔₋. (15)

To see if 𝜏⁽¹⁾_1,𝑛 and 𝜏⁽²⁾_1,𝑛 are bifurcation values, we need to verify if the transversality conditions hold.

Lemma 2. The following transversality conditions 𝑑Re𝜆_1,𝑛𝜏⁽¹⁾_1,𝑛𝑑𝜏₁>0,𝑑Re𝜆_1,𝑛𝜏⁽²⁾_1,𝑛𝑑𝜏₁>0, (16) hold.

Proof. Differentiating (8) with respect to 𝜏₁ yields 2𝜆𝑑𝜆𝑑𝜏₁+𝐴𝑑𝜆𝑑𝜏₁+𝐶(𝑘+𝑑)𝑒^−𝜆𝜏1−𝜏₁𝑑𝜆𝑑𝜏₁−𝜆=0,𝑑𝜆𝑑𝜏₁=𝐶(𝑘+𝑑)𝜆𝑒^−𝜆𝜏12𝜆+𝐴−𝐶(𝑘+𝑑)𝜏₁𝑒^−𝜆𝜏1; (17) then we have Re𝑑𝜆𝑑𝜏₁=Re2𝜆+𝐴𝐶(𝑘+𝑑)𝜆𝑒^−𝜆𝜏1−𝜏₁𝜆=2𝜆+𝐴𝐶(𝑘+𝑑)𝜆𝑒^−𝜆𝜏1=1𝐶()⋅𝑑+𝑘2𝑖𝜔+𝐴𝑖𝜔cos𝜔𝜏₁−𝑖sin𝜔𝜏₁=1𝐶()⋅𝑑+𝑘𝐴𝜔sin𝜔𝜏₁+2𝜔²cos𝜔𝜏₁𝜔², (18) which satisfied 𝐶>0 and 𝑑+𝑘>0; then for 𝜆_1,𝑛(𝜏⁽¹⁾_1,𝑛), 𝑑Re𝜆_1,𝑛(𝜏⁽¹⁾_1,𝑛)/𝑑𝜏₁>0, and for 𝜆_1,𝑛(𝜏⁽²⁾_1,𝑛), 𝑑Re𝜆_1,𝑛(𝜏⁽²⁾_1,𝑛)/𝑑𝜏₁>0. We complete the proof.

Thus, we get the distribution of the characteristic roots of (8).

Lemma 3. For (8), one has the following: (i)if (H3) and either (1) (H1) 𝐴>0 or (2) (H2) 𝐴>0 and 𝐵𝑘+𝐶𝑑+𝐶𝑘>0 hold, then when 𝜏₁∈[0,𝜏⁽¹⁾_1,0), all roots of (8) have negative real parts, and when 𝜏₁>𝜏⁽¹⁾_1,0, (8) has at least one root with positive real part;(ii)if (H4) and either (H1) or (H2) hold, then there are k switches from stability to instability; that is, when 𝜏₁∈(𝜏⁽²⁾_1,𝑛,𝜏⁽¹⁾_1,𝑛+1), 𝑛=−1,0,1,…,𝑘−1, all roots of (8) have negative real parts, where 𝜏⁽²⁾_2,−1=0, 𝜏₁∈[𝜏⁽¹⁾_1,𝑛,𝜏⁽²⁾_1,𝑛+1) and 𝜏₁>𝜏⁽¹⁾_1,𝑘, 𝑛=0,1,…,𝑘−1, and (8) has at least one root with positive real part.That is to say, under those conditions when 𝜏₂=0, 𝜏₁∈[0,𝜏⁰₁), system (4) is asymptotically stable, and system (4) undergoes a Hopf bifurcation when 𝜏₁=𝜏⁰₁.

Then, we consider stable interval for 𝜏₁ in which all roots of (8) have negative real parts, regarding 𝜏₂>0 as a parameter.

Lemma 4. If all roots of (8) have negative real parts, then there exists a 𝜏⁰₂(𝜏₁)>0, such that when 𝜏₂∈[0,𝜏⁰₂(𝜏₁)) all roots of (7) have negative real parts.

Proof. All roots of (8) have negative real parts which means that system (4) is stable when 𝜏₂=0. In what follows, we consider (4) with fixed 𝜏₁ in its stable interval, regarding 𝜏₂ as a parameter. Let 𝑖𝑣(𝑣>0) be a root of (4); then we obtain −𝑣²+𝐴cos𝑣𝜏₂−𝑖sin𝑣𝜏₂𝑣𝑖+𝐵𝑘cos𝑣𝜏₂−𝑖sin𝑣𝜏₂+𝐶𝑑cos𝑣𝜏₁−𝑖sin𝑣𝜏₁+𝐶𝑘cos𝑣𝜏₁+𝜏₂−𝑖sin𝑣𝜏₁+𝜏₂=0. (19) Suppose that 𝐹(𝑣)=𝑣²−𝐴𝑣𝑖cos𝑣𝜏₂−𝐴𝑣sin𝑣𝜏₂−𝐵𝑘cos𝑣𝜏₂+𝐵𝑘𝑖sin𝑣𝜏₂−𝐶𝑑cos𝑣𝜏₂+𝐶𝑑𝑖sin𝑣𝜏₂−𝐶𝑘cos𝑣(𝜏₁+𝜏₂)+𝐶𝑘𝑖sin𝑣(𝜏₁+𝜏₂). Since 𝐹(0)=−(𝐵𝑘+𝐶𝑑+𝐶𝑘)<0 and 𝐹(+∞)=+∞, then (19) has at least one positive root. Without loss of generality, the roots of (19) are defined by 𝑣₁,𝑣₂,…,𝑣_𝑘. For every 𝑣_𝑖(𝑖=1,2,…,𝑘), there exists a sequence {𝜏^(𝑗)_2𝑖∣𝑗=1,2,…}, such that (19) holds. The expression of 𝜏^(𝑗)_2𝑖 and 𝑣_𝑖 can be derived by (19) for fixed 𝜏₁; we will calculate them directly by the use of Mathematica software in Section 4; here we omit them. Let 𝜏⁰₂={min𝜏^𝑗_2𝑖∣𝑖=1,2,…,𝑘,𝑗=1,2,…} and let 𝑣₀ be the positive and simple root of (19) when 𝜏₂=𝜏⁰₂. When 𝜏₂=𝜏⁰₂, (19) has a pair of purely imaginary roots ±𝑖𝑣^∗ for 𝜏₁∈[0,𝜏⁰₁). Then, as 𝜏₂ varies, the sum of the multiplicities of zeros in the open right half-plane can change only if a zero appears on or crosses the imaginary axis. In what follows, we assume that (H6) [𝑑Re𝜆(𝜏₂)/𝑑𝜏₂]_𝜏2=𝜏^𝑗_2𝑖≠0. Therefore, by the general Hopf bifurcation theorem for FDEs in Wei and Ruan [10], we get that when 𝜏₂∈[0,𝜏⁰₂) all roots of (7) have negative real parts.
Applying the discussion above and noticing that all roots of (8) have negative real parts, we know that there exist 𝜏⁰₂>0 such that all roots of (7) with 𝜏₂∈[0,𝜏⁰₂) have negative real parts. The proof is complete.

Summarizing the above lemmas and literature of Hale [11], we obtain the following sufficient conditions for all characteristic roots of (8) to have negative real parts.

Theorem 5. Suppose that (H6) holds and either (H1) or (H2) is satisfied. (i)If (H3) holds, then for any 𝜏₁∈[0,𝜏⁽¹⁾_1,0) there exists a 𝜏⁰₂(𝜏₁)>0 such that when 𝜏⁰₂∈[0,𝜏₂(𝜏₁)), all roots of (7) have negative real parts.(ii)If (H4) holds, then for any 𝜏₁∈∪^𝑘−1_𝑛=−1(𝜏⁽²⁾_1,𝑛,𝜏⁽¹⁾_1,𝑛+1) there exists a 𝜏₂(𝜏₁)>0, such that when 𝜏₁∈[0,𝜏⁽¹⁾_1,0) all roots of (7) have negative real parts, where 𝜏⁽¹⁾_1,𝑗 and 𝜏⁽²⁾_1,𝑗 are defined by (10).(iii)If (H5) holds, then for any 𝜏₁≥0, there exists a 𝜏₂(𝜏₁)>0, such that when 𝜏₂∈[0,𝜏₂(𝜏₁)) all roots of (7) have negative real parts.That is to say, under the conditions that 𝜏₁ is stable interval, there exist a 𝜏₂ (depend on 𝜏₁) such that when 𝜏₂∈[0,𝜏⁰₂), system (4) is asymptotically stable, and system (4) undergoes a Hopf bifurcation when 𝜏₂=𝜏^(𝑗)_2𝑖, 𝑖=1,2,…,𝑘; 𝑗=1,2….

3. Direction and Stability of the Bifurcating Periodic Solutions

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcate from the steady state and the equilibrium loses its stability when 𝜏₂=𝜏^(𝑗)_2𝑖, 𝑖=1,2,…,𝑘; 𝑗=1,2,… for fixed 𝜏₁, and the relationship between 𝜏₁ and 𝜏₂ can be derived by (19). Throughout this section, by using techniques of the normal form and center manifold theory due to Hale [11], we derive the algorithm for determining the direction of the Hopf bifurcations and the stability of the bifurcating periodic solutions at critical values on the center manifold.

Since the analysis is local, we regard 𝜏₂=𝜏⁰₂+𝛾,  𝛾∈𝑅 as bifurcation parameter. Choosing the space as 𝐶=𝐶([−𝜏⁰₂,0],𝑅²) and 𝑢_𝑡=𝑢(𝑡+𝜃)∈𝐶 for 𝜃∈[−𝜏⁰₂,0], system (3) is transformed into FDE as 𝑢̇(𝑡)=𝐿_𝛾𝜙+𝐹𝜙,𝛾, (20) with 𝐿_𝛾𝜙=𝐵₁𝜙(0)+𝐵₂𝜙−𝜏₁+𝐵₃𝜙−𝜏₂,𝐹=0𝜙,𝛾𝜇𝑏𝜙²₁(0)𝜙₂−𝜏₂+2𝜇𝑃₀+𝛽𝜙₁(0)𝜙₂(0), (21) where 𝐵₁=^0𝑑_{−𝜇𝑏00}, 𝐵₂=⁰⁰_−𝜇𝑎0, and 𝐵₃=^0𝑘_{02𝜇(𝑃0+𝛽)}²_+𝐶0. Obviously, 𝐿_𝛾(𝜙) is continuous linear function mapping 𝐶([−𝜏⁰₂,0],𝑅²) into 𝑅². By the Riesz representation theorem, there exists a matrix whose elements are bounded variation functions 𝜂(𝜃,𝛾) in 𝜃∈[−𝜏⁰₂,0] such that 𝐿_𝛾𝜙=⁰_−𝜏02𝑑𝜂𝜃,𝛾𝜙(𝜃),for𝜙∈𝐶. (22) In fact, we choose 𝜂𝜃,𝛾=𝐵₁𝛿(𝜃)+𝐵₂𝛿𝜃+𝜏₁+𝐵₃𝛿𝜃+𝜏₂, (23) where 𝛿(𝜃) is a delta function.

For 𝜙∈𝐶^′([−𝜏⁰₂,0],𝑅²) the operators 𝐴 and 𝑅 are defined as 𝐴𝜇𝜙(𝜃)=𝑑𝜙(𝜃)𝑑𝜃,𝜃∈−𝜏⁰₂,0_0−𝜏02,𝑑𝜂𝜙(𝑡𝑅𝜇𝜙(𝜃)=𝑡,𝜇),𝜃=0,0,𝜃∈−𝜏⁰₂,𝑓,0𝛾,𝜃,𝜃=0. (24) Let 𝜓∈𝐶^′[0,𝜏⁰₂]; the adjoint operator 𝐴^∗(0) corresponding to 𝐴(0) is defined as follows: 𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)𝑑𝑠,𝑠∈0,𝜏⁰₂,_0−𝜏02𝑑𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0. (25) Then system (20) can be written in the following form: 𝑢̇_𝑡=𝐴(𝛼)𝑢_𝑡+𝑅(𝛼)𝑢_𝑡, where 𝑢_𝑡=𝑢(𝑡+𝜃) for 𝜃∈[−1,0).

For 𝜙∈𝐶^′([−𝜏⁰₂,0],𝑅²) and 𝜓∈𝐶^′[0,𝜏⁰₂], define the adjoint bilinear: 𝜓,𝜙=𝜓(0)𝜙(0)−⁰_−𝜏02^𝜃_𝜀=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (26) where 𝜂(𝜃)=𝜂(𝜃,0).

Proposition 6. Let 𝑞(𝜃) and 𝑞^∗(𝑠) be eigenvectors of 𝐴 and 𝐴^∗ corresponding to 𝑖𝜔₀ and −𝑖𝜔₀, respectively, satisfying ⟨𝑞^∗,𝑞⟩=1 and ⟨𝑞^∗,𝑞⟩=0. Then 𝑞(𝜃)=𝑞₁,𝑞₂^𝑇𝑒^{𝑖𝜔0𝜃}=𝑑+𝑘𝑒^{−𝑖𝜔0𝜏2},𝑖𝜔^𝑇𝑒^{𝑖𝜔0𝜃},𝑞^∗(𝑠)=𝐷𝑞^∗₁,𝑞^∗₂𝑒^{−𝑖𝜔0𝑠}=𝜇𝑏₀+𝜇𝑎𝑒^{𝑖𝜔0𝜏1},𝑖𝜔𝑒^{−𝑖𝜔0𝑠}, (27) where 𝐷=𝑞₁𝑞^∗₁+𝑞₂𝑞^∗₂+𝜏⁰₁−𝜇𝑎𝑒^{−𝑖𝜔0𝜏01}𝑞₁𝑞^∗₂+𝑘𝑞₂𝑞^∗₁+2𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑞^∗₂𝑞₂𝜏⁰₂𝑒^{−𝑖𝜔0𝜏02−1}. (28)

Proof. We assume that 𝑞(𝜃)=(𝑞₁,𝑞₂)^𝑇𝑒^{𝑖𝜔0𝜃} is the eigenvector of 𝐴(0) corresponding to 𝑖𝜔 and 𝑞^∗(𝑠)=𝐷(𝑞^∗₁,𝑞^∗₂)^𝑇𝑒^{−𝑖𝜔0𝑠} is the eigenvector of 𝐴^∗(0) corresponding to −𝑖𝜔. It follows from the definition of 𝐴(0), 𝐴^∗(0), (22), and (23) that we have 𝐴𝑞(0)=∫⁰_−𝜏02𝑑(𝜂(𝑡,𝜇)𝜙(𝑡))=𝑖𝜔₀𝑞(0) and 𝐴^∗𝑞(0)=∫⁰_−𝜏02𝑑(𝜂^𝑇(𝑡,𝜇)𝜙(−𝑡))=−𝑖𝜔₀𝑞(0); we have 𝑖𝜔₀𝐼−𝐵₁+𝐵₂𝑒^{−𝑖𝜔0𝜏1}+𝐵₃𝑒^{−𝑖𝜔0𝜏2}𝑞(0)=0,−𝑖𝜔₀𝐼−𝐵₁+𝐵₂𝑒^{𝑖𝜔0𝜏1}+𝐵₃𝑒^{𝑖𝜔0𝜏2}𝑞^∗(0)=0, (29) where 𝐼 is identity matrix; that is, 𝑖𝜔₀−𝑑+𝑘𝑒^{−𝑖𝜔0𝜏2}𝜇𝑏₀+𝑎𝑒^{−𝑖𝜔0𝜏1}𝑖𝜔₀−𝜇𝑏𝑃₀+𝛽+𝐶₀𝑒^{−𝑖𝜔0𝜏2}×𝑞₁𝑞₂=0,−𝑖𝜔₀𝜇𝑏₀+𝑎𝑒^{𝑖𝜔0𝜏1}−𝑑+𝑘𝑒^{𝑖𝜔0𝜏2}−𝑖𝜔₀−𝜇𝑏𝑃₀+𝛽+𝐶₀𝑒^{𝑖𝜔0𝜏2}×𝑞^∗₁𝑞^∗₂=0. (30) By direct computation and considering 𝑞(𝜃)=𝑞(0)𝑒^𝑖𝜔𝜃, 𝑞^∗(𝑠)=𝑞^∗(0)𝑒^{−𝑖𝜔𝑠}, we obtain 𝑞(𝜃) and 𝑞^∗(𝑠). Now, we calculate ⟨𝑞^∗,𝑞⟩ as follows: 𝑞(𝑠)^∗(𝜃=,𝑞)𝐷𝑞^∗₁,𝑞^∗₂𝑞₁𝑞₂−⁰_−𝜏02^𝜃_𝜀=0𝑞^∗₁𝑞^∗₂^𝑇𝑒^{−𝑖𝜔0(𝜀−𝜃)}(𝜃)𝑞𝑑𝜂¹𝑞₂𝑒^{𝑖𝜔0𝜀}=𝑑𝜀𝐷𝑞₁𝑞^∗₁+𝑞₂𝑞^∗₂−⁰_−𝜏02𝑞^∗₁𝑞^∗₂^𝑇𝜃𝑒^{𝑖𝜔0𝜃}(𝜃)𝑞𝑑𝜂¹𝑞₂=𝐷𝑞₁𝑞^∗₁+𝑞₂𝑞^∗₂+𝜏⁰₁−𝜇𝑎𝑒^{−𝑖𝜔0𝜏01}𝑞₁𝑞^∗₂+𝑘𝑞₂𝑞^∗₁+2𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑞^∗₂𝑞₂𝜏⁰₂𝑒^{−𝑖𝜔0𝜏02}=1. (31) Since ⟨𝜓,𝐴𝜙⟩=⟨𝐴^∗𝜓,𝜙⟩, we get 𝑖𝜔𝑞^∗,𝑞=𝐴^∗𝑞^∗,𝑞=−𝑖𝜔₀𝑞^∗,𝑞=𝑖𝜔₀𝑞^∗,𝑞. (32) Therefore, ⟨𝑞^∗,𝑞⟩=0. This completes the proof.

Then, we construct the coordinates of the center manifold 𝐶₀ at 𝛾=0. Let 𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃). (33) On the center manifold 𝐶₀, we have 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), (34) where 𝑊(𝑧,𝑧,𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧+𝑊₀₂𝑧²2+𝑊₃₀𝑧³6+⋯, (35) and 𝑧 and 𝑧 are local coordinates for the center manifold 𝐶₀ in the direction of 𝑞 and 𝑞^∗, respectively. Since 𝛾=0, we have 𝑧^′(𝑡)=𝑖𝜔₀𝑧(𝑡)+𝑞^∗(𝜃),𝑓𝑊+2Re𝑧(𝑡)𝑞(𝜃)=𝑖𝜔₀𝑧(𝑡)+𝑞^∗(0)𝑓𝑊(𝑧,𝑧,0)+2Re𝑧(𝑡)𝑞(0)≜𝑖𝜔₀𝑧(𝑡)+𝑞^∗(0)𝑓₀(𝑧,𝑧), (36) where 𝑓₀(𝑧,𝑧)=𝑓_𝑧2𝑧²2+𝑓_𝑧²𝑧²2+𝑓_𝑧𝑧𝑧𝑧+⋯. (37) We rewrite in abbreviated form as 𝑧^′(𝑡)=𝑖𝜔₀𝑧+𝑔(𝑧,𝑧), (38) where 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (39)

By (20) and (38), we obtain 𝑊̇=𝑢̇_𝑡−𝑧𝑞̇−𝑧̇𝑞=𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝜃),𝜃∈−𝜏⁰₂,0,𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝜃)+𝑓₀(,𝜃=0≜𝐴𝑊+𝐻𝑧,𝑧,𝜃), (40) where 𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯. (41) Substituting (26) and (38) into 𝑊̇=𝑊_𝑧𝑧̇+𝑊_𝑧𝑧̇‌ on the center manifold 𝐶₀ and comparing the coefficients we get 𝐴−2𝑖𝜔₀𝐼𝑊₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃),𝐴+2𝑖𝜔₀𝐼𝑊₀₂(𝜃)=−𝐻₀₂(𝜃). (42) Comparing the coefficients with (41) gives that 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (43) From (42), (43), and the definition of 𝐴, we can derive the following equation: 𝑊̇₂₀(𝜃)=2𝑖𝜔₀𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₂₀𝑞(𝜃),𝑊̇₁₁(𝜃)=𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃). (44) Solving for 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃), we get 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝑞(0)𝑒^{𝑖𝜔0𝜃}−𝑔₀₂3𝑖𝜏₀𝜔₀𝑞(0)𝑒^{−𝑖𝜔0𝜃}+𝐸₁𝑒^{2𝑖𝜔0𝜃},𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝑞(0)𝑒^{𝑖𝜔0𝜃}+𝑖𝑔₁₁𝜔₀𝑞(0)𝑒^{−𝑖𝜔0𝜃}+𝐸₂, (45) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸²₁)∈𝑅², 𝐸₂=(𝐸₂⁽¹⁾,𝐸²₂)∈𝑅² are two constant vectors, which can be determined by setting 𝜃=0 in 𝐻.

By (33), 𝜙_𝑡(𝜃)=[𝑢₁(𝑡−𝜏₁),𝑢₂(𝑡−𝜏₂)]=𝑊(𝑡,𝜃)+𝑧𝑞(𝜃)+𝑧𝑞(𝜃), and noticing 𝑞(𝜃)=(𝑞₁,𝑞₂)^𝑇𝑒^{𝑖𝜔0𝜃} we have 𝜙_𝑡(𝜃)=𝑧𝑑+𝑘𝑒^{−𝑖𝜔0𝜏02}𝑒𝑖𝜔^{−𝑖𝜔0𝜃}+𝑧𝑑+𝑘𝑒^{𝑖𝜔0𝜏02}𝑒−𝑖𝜔^{−𝑖𝜔0𝜃}+𝑊(𝑡,𝜃). (46) Then it is easy to obtain 𝜙₁(0)=𝑧𝑑+𝑘𝑒^{−𝑖𝜔0𝜏02}+𝑧𝑑+𝑘𝑒^{𝑖𝜔0𝜏𝑜2}+𝑊⁽¹⁾(0),𝜙₂(0)=𝑧𝜔𝑖+𝑧(−𝑖𝜔)+𝑊⁽²⁾(0),𝜙₂𝑡−𝜏⁰₂=𝑧𝜔𝑖𝑒^{−𝑖𝜔0𝜏02}−𝑧𝑒^{𝑖𝜔0𝜏02}+𝑊²𝑡−𝜏⁰₂, (47) where 𝑊⁽¹⁾(0)=𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧2+𝑜|𝑧,𝑧|³,𝑊⁽²⁾(0)=𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧2+𝑜|𝑧,𝑧|³,𝑊⁽²⁾𝑡−𝜏₂=𝑊⁽²⁾₂₀−𝜏⁰₂𝑧²2+𝑊⁽²⁾₁₁−𝜏⁰₂𝑧𝑧+𝑊⁽²⁾₀₂−𝜏⁰₂𝑧2+𝑜|𝑧,𝑧|³,𝑓₀(𝑧,𝑧)=0𝜇𝑏𝜙²₁(0)𝜙₂−𝜏⁰₂+2𝜇𝑃₀+𝛽𝜙₁(0)𝜙₂(0). (48) Thus, from (39), it follows that 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝐷𝑞^∗₁𝑞^∗₂×0𝜇𝑏𝜙²₁(0)𝜙₂−𝜏⁰₂+2𝜇𝑃₀+𝛽𝜙₁(0)𝜙₂(0)=𝐷𝑞^∗₂𝜇𝜙²₁(0)𝜙₂−𝜏⁰₂+2𝜇𝑃₀+𝛽𝜙₁(0)𝜙₂(0=)𝐷2𝜇𝑃₀+𝛽𝑞^∗₂𝑞₁𝑞₂𝑧²+2𝜇𝑃₀×+𝛽𝑞^∗₂𝑞₁𝑞₂+𝑞₁𝑞₂𝑧𝑧+2𝜇𝑃₀+𝛽𝑞^∗₂𝑞₁𝑞₂𝑧²+𝜇𝑏𝑞^∗₂𝑞²₁𝑞₂𝑒^{𝑖𝜔0𝜏02}+2𝑞₁𝑞₁𝑞₂𝑒^{−𝑖𝜔0𝜏02}+2𝜇𝑃₀+𝛽𝑞^∗₂×𝑞₂𝑊⁽¹⁾₁₁(0)+𝑞₁𝑊⁽²⁾₁₁(0)+𝑞₂𝑊⁽¹⁾₂₀(0)+𝑞₁𝑊⁽²⁾₂₀(0)2𝑧²𝑧+𝑜𝜙⁴. (49) Comparing the coefficient with (39), we have 𝑔₂₀=4𝐷𝜇𝑃₀+𝛽𝑞^∗₂𝑞₁𝑞₂,𝑔₁₁=4𝐷𝜇𝑃₀+𝛽𝑞^∗₂Re𝑞₁𝑞₂,𝑔₀₂=4𝐷𝜇𝑃₀+𝛽𝑞^∗₂𝑞₁𝑞₂,𝑔₂₁=2𝐷𝜇𝑏𝑞^∗₂𝑞²₁𝑞₂𝑒^{𝑖𝜔0𝜏02}+2𝑞₁𝑞₁𝑞₂𝑒^{−𝑖𝜔0𝜏02}+2𝜇𝑃₀+𝛽𝑞^∗₂×𝑞₂𝑊⁽¹⁾₁₁(0)+𝑞₁𝑊⁽²⁾₁₁(0)+𝑞₂𝑊⁽¹⁾₂₀(0)+𝑞₁𝑊⁽²⁾₂₀(0)2. (50) From (40), we get that 𝐻(𝑧,𝑧,0)=−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀=−2Re(𝑔𝑞(0))+𝑓₀=−𝑔𝑞(0)−𝑔𝑞(0)+𝑓₀; that is, 𝐻₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₂₀𝑞(0)+04𝜇𝑃₀+𝛽𝑞₁𝑞₂,𝐻₁₁(0)=−𝑔₁₁𝑞(0)−𝑔₁₁𝑞(0)+02𝑃₀+𝛽𝑞₁𝑞₂+𝑞₁𝑞₂. (51) By the definition of 𝐴 and (42) we have 𝑑𝜂(⁰_−𝜏02𝜃)𝑊₂₀(𝜃)=2𝑖𝜔₀𝑊₂₀−𝐻₂₀(0),_0−𝜏02𝑑𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0). (52) Notice that 𝑖𝜔₀𝐼−⁰_−𝜏02𝑑𝜂(𝜃)𝑒^{𝑖𝜔0𝜃}𝑞(0)=0,−𝑖𝜔₀𝐼−⁰_−𝜏02𝑑𝜂(𝜃)𝑒^{𝑖𝜔0𝜃}𝑞(0)=0. (53) Substituting (45) and (52) into (51), we obtain 2𝑖𝜔₀−𝑑−𝑘𝑒^{2𝑖𝜔0𝜏02}𝜇𝑏₀+𝜇𝑎𝑒^{−2𝑖𝜔0𝜏01}2𝑖𝜔₀−𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑒^{2𝑖𝜔0𝜏02}𝐸₁=04𝜇𝑃₀+𝛽𝑞₁𝑞₂,𝜇𝑏0−𝑑−𝑘₀+𝑎−𝜇𝑏𝑃₀+𝛽²+𝐶₀𝐸₂=02𝑃₀+𝛽𝑞₁𝑞₂+𝑞₁𝑞₂. (54) By direct computation, we obtain𝐸₁⁽¹⁾=4𝜇𝑃₀+𝛽𝑞₁𝑞₂𝑑+𝑘𝑒^{−2𝑖𝜔0𝜏02}𝜇𝑏₀+𝜇𝑎𝑒^{−2𝑖𝜔0𝜏01}𝑑+𝑘𝑒^{−2𝑖𝜔0𝜏02}−2𝑖𝜔₀𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑒^{−2𝑖𝜔0𝜏02}−2𝑖𝜔₀,𝐸₁⁽²⁾=−6𝑖𝜔₀𝜇𝑃₀+𝛽𝑞₁𝑞₂−𝜇𝑏₀+𝜇𝑎𝑒^{−2𝑖𝜔0𝜏01}𝑑+𝑘𝑒^{−2𝑖𝜔0𝜏02}+2𝑖𝜔₀𝜇𝑏𝑃₀+𝛽²+𝐶₀𝑒^{−2𝑖𝜔0𝜏02}−2𝑖𝜔₀,𝐸₂⁽¹⁾=2𝑃₀+𝛽𝑞₁𝑞₂+𝑞₁𝑞₂𝜇𝑏₀−𝑎,𝐸₂⁽²⁾=0. (55)

Thus, we can compute the following values which determine the properties of bifurcating periodic solutions at the critical value 𝜏⁰₂: 𝐶₁(0)=𝑖2𝜏⁰₂𝜔₀𝑔₂₀𝑔₁₁−2𝑔₁₁−13𝑔₀₂²+𝑔₂₁2,𝜇₂𝐶=−Re₁(0)Re𝜆^′𝜏⁰₂,𝛽₂=2Re𝐶₁(0,𝑇₂)𝐶=−Im₁(0)+𝜇₂Im𝜆^′𝜏⁰₂𝜔₀. (56) More specifically (see Hassard et al. [12]), 𝜇₂ determines the direction of the Hopf bifurcation: if 𝜇₂>0 (𝜇₂<0), then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for 𝜏>𝜏⁰₂ (𝜏<𝜏⁰₂). 𝛽₂ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if 𝛽₂<0 (𝛽₂>0). 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0 (𝑇₂<0).

4. Numerical Simulation

In order to validate the theoretical analysis, we will present some numerical simulations. We take the following coefficients as an example: 𝜇=0.13, 𝑏₀=10, 𝑎=12, 𝑑₀=60, 𝑔₀=20, 𝛽=−2, 𝑏=1, 𝑘=1, and 𝑑=1, it is easy to obtain 𝐸^∗=(1.82,0), 𝜔₀≈0.89, and 𝜏⁰₁≈0.72. We choose 𝜏⁰₁=0.7; then from (19) and by means of Mathematica software, we get 𝑣₀≈2 and 𝜏⁰₂≈0.86. Taking (𝑥(0),𝑦(0)=(2,2)) as the initial conditions. By Theorem 5 and the above results, we know equilibrium 𝐸^∗ is locally asymptotically stable when 𝜏₂=0.8<𝜏⁰₂ as is illustrated in Figure 1. When 𝜏₂>𝜏⁰₂, 𝐸^∗ is unstable and periodic solutions occur from 𝐸^∗; we take 𝜏₂=0.9, the corresponding phase plots are shown in Figure 2.

Figure 1 

Equilibrium 𝐸∗ is locally asymptotically stable when 𝜏1=0.7, 𝜏2=0.8<𝜏02.

Figure 2 

Equilibrium 𝐸∗ is unstable when 𝜏1=0.7, 𝜏2=0.9>𝜏02.

Finally, the numerical simulation shows that it is a complex transformation process for the system changes from stable equilibrium to chaos.

5. Conclusion

Different from the previous work in [1, 2, 4], the main contribution of this paper lies in the following aspects. Firstly, we modify the traditional price differential equation model by considering demand function that is settled as a downward opening quadratic function and considering supply and demand functions that are both depending on the price of the past and the present. Then the price differential equation model with two delays is established. Secondly, to study the stability and Hopf bifurcation of system (4), we consider (4) with one delay 𝜏₁ and find the stable interval for 𝜏₁. In the following, regarding 𝜏₂ as a parameter, we obtain the stable interval for 𝜏₂. Then we get the stable interval for (4). In addition, we derive the algorithm for determining the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions by using techniques of normal form theory and center manifold method. Lastly, a numerical analysis confirms the effectiveness of our research results. The paper provided the preparative work for further discussion. For instance, we consider that supply cannot increase with the price without limit; the production capacity of enterprises and social resources are limited; then we can modify the supply function as a fractional linear function of price. The results in the paper enrich the toolbox for the qualitative analysis of mathematical economics and business administration.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors are grateful to the referees for their helpful comments and constructive suggestions.