# Dynamic Analysis on a Diffusive Two-Enterprise Interaction Model with Two Delays

Authors: Yanxia Zhang; Long Li
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The oscillation and instability of systems caused by time delays have been widely studied over the past several decades. In nature, the phenomenon of diffusion is universal. Therefore, it is necessary to investigate the dynamic behavior of reaction-diffusion systems with time delays. In this study, a two-enterprise interaction model with diffusion and delay effects is considered. By analyzing the distribution of the roots of the corresponding characteristic equation, some conditions for the stability of the unique positive equilibrium and the existence of Hopf bifurcation at the steady state are investigated. As the sum of the time delays changes, there are a series of periodic solutions at the trivial steady-state solution of the system. In addition, the direction of Hopf bifurcation and the stability of the periodic solutions are discussed by using the normal form theory and the center manifold reduction of partial functional differential equations. Finally, numerical simulation experiments are conducted to illustrate the validity of the theoretical conclusions.


---

## Body

## 1. Introduction

As we all know, it is of great significance for an enterprise to consider its developmental progress and growth trend, grasp its position, and comprehensively analyze its interactive relationships with other enterprises in the business ecosystem. In recent years, the dynamic relationships among enterprises described by corresponding mathematical models have been widely studied [1–6]. In [5], Xu proposed the following model for the dynamic competition and cooperation between two enterprises and applied coincidence degree theory to analyze the existence of periodic solutions:(1)x˙1t=r1x1t1−x1tK1−αx2t−c22K2,x˙2t=r2x2t1−x2tK2+βx1t−c12K1,x10≥0,x20≥0,where xiti=1,2 denotes the output of two enterprises, respectively; rii=1,2 represents the intrinsic growth rate of two enterprises, respectively; K1 and K2 are the carrying capacities of the natural market, i.e., load capacities of the two enterprises under unrestricted conditions; α is the consumption coefficient from the enterprise with output x2t to the one with output x1t; β is the transformation coefficient from the enterprise with output x1t to the one with output x2t; and cii=1,2 represents the initial output of the two enterprises, respectively. It was found that the outputs of the two enterprises reach an equilibrium state and at least have a positive periodic solution when the intrinsic growth rates, natural market loads, competition coefficients, and initial productions of the two enterprises meet certain conditions. Considering that the real dynamic ecological system always depends on information from the historical system, Liao et al. [3, 4] proposed two different competition and cooperation models for two enterprises with a delay (or delays) in succession. They investigated the effect of the delay (or delays) on the stability of the positive equilibrium by choosing τ (or τ1,τ2) as the bifurcation parameter (or parameters). On the basis of [3], Guerrini [1] further proposed that when the two delays are small, an approximate, simple method can be used to obtain the roots of the characteristic equation. An analysis of the dynamic behavior of a competition and cooperation model for two enterprises with two small delays was also provided. Considering the influence of periodic solutions on the development and operation of enterprises, Xu and Li [6] investigated almost periodic solutions for a competition and cooperation model of two enterprises with time-varying delays and feedback controls. The existence and global stability of the almost periodic solutions were obtained by establishing a suitable Lyapunov functional. In fact, the internal and external effects of time delays were found to be completely different for enterprises. Hence, Li et al. [2] proposed the following differential equations model with four different delays:(2)x˙1t=r1x1t1−x1t−τ1K1−αx2t−τ2−c22K2,x˙2t=r2x2t1−x2t−τ3K2+βx1t−τ4−c12K1,x1t=φt,x2t=ψt,t∈−max1≤i≤4τi,0,where φt,ψt∈C−max1≤i≤4τi,0,ℝ, and τi≥01≤i≤4 represent the delay effects within and between enterprises, respectively. Specifically, τ1 and τ3 denote the gestation periods of the outputs for the two enterprises. τ2 and τ4 represent the block delay and promoting delay, respectively, imposed by one enterprise on another. One case where τi=τ1≤i≤4 in (2) was considered in [4], another case where τ3=τ1, τ4=τ2 was discussed in [3], and some corresponding results were obtained. In [2], Li et al. researched the stability of the positive equilibrium and the existence of Hopf bifurcation, deduced the formula determining the direction of Hopf bifurcation, and achieved the stability of bifurcating periodic solutions by setting τ1=τ3=0 and choosing τ=τ2+τ4>0 as the parameter. It has shown that time delay has an important impact on the steady-state equilibrium of the interactive relationship between enterprises. The latest research results of some delayed systems can be found in [7–13], which focus on multiple time delays and feedback control, fractional-order delay model, and discrete model with multiple delays.However, most of these studies only considered the effect of time delays and lost sight of the effect of spatial diffusion. It is well known that diffusion factors are indispensable in the modeling of ecological and biological systems and often utilized to describe the spatial distributions of the densities of some certain substances, such as plants, animals, and other organisms. At present, there are a growing number of dynamic models described by reaction-diffusion equations with time delays in various application areas [14–27]. For example, some scholars [19, 23, 25–27] proposed several different diffusive predator-prey systems with delays to analyze the effect of diffusion on biological systems. Huang et al. [18] researched a plant-pollinator model with diffusion and time delay effects and discussed its bifurcation and temporal periodic patterns. Some scholars [16, 22] established delayed reaction-diffusion systems to investigate the spatiotemporal transmission of bacteria and viruses in epidemiology. As a result of the impact of the growth pole on the economic activities of industrial groups, Hu et al. [17] provided a Kaldor–Kalecki model of the business cycle with diffusion and time delay and described the influence of the diffusion coefficient on the spatial pattern of the system.Generally, in the development of the regional economy, the enterprises located in the center of the city have the functions of production, information, transportation, trade, service, and decision-making, which act like a “magnetic pole” that attracts production from the surrounding areas and spreads its own production factors to the surrounding areas. Especially when the industrial economic development in the central area is at the mature stage, the diffusion effect continues to increase. Therefore, based on the above argument, we further consider spatial diffusion with zero-flux boundary conditions. A delayed reaction-diffusion system for a two-enterprise interaction should be described by the following model:(3)∂x1∂t=D1Δx1+r1x1t,x1−x1t,xK1−αx2t−τ2,x−c22K2,x∈Ω,t>0,∂x2∂t=D2Δx2+r2x2t,x1−x2t,xK2+βx1t−τ1,x−c12K1,x∈Ω,t>0,∂x1∂v=∂x2∂v=0,t≥0,x∈∂Ω,x1t,x=φt,x≥0,x2t,x=ψt,x≥0,t,x∈−max1≤i≤2τi,0×Ω¯,where xit,xi=1,2 denotes the output of two enterprises at time t and spatial position x, respectively. xit−τi,xi=1,2 represents the output of two enterprises in the gestation period at spatial state x, respectively. 1−xit,x/Kii=1,2 reflects the prevention effect on the output of enterprise ii=1,2 caused by its own consumption of limited resources, respectively. As the two enterprises exist in the same environment, the impact of the delay, caused by enterprise 2 consuming the same limited resources in the process of development, on enterprise 1 is expressed by −αx2t−τ2,x−c22/K2, and the promoting delay by which enterprise 1 supports enterprise 2 is denoted as βx1t−τ1,x−c12/K1. In [3], it was suggested that the promoting delay and blocking delay in the gestation period of output for enterprises have important economic significance for their development through changing parameters α and β. Δ is the Laplacian operator, and D1,D2>0 denote the diffusion coefficients. Ω represents a bounded domain in ℝNN≥1 with a smooth boundary ∂Ω, and v is the outward unit normal vector on ∂Ω. Similar to [2–4], by introducing the terms αi=ri/Kii=1,2,b1=αr1/K2,b2=βr2/K1,di=ri−aici,yit,x=xit,x−cii=1,2, we obtain the following model:(4)∂y1∂t=D1Δy1+y1t,x+c1d1−a1y1t,x−b1y22t−τ2,x,x∈Ω,t>0,∂y2∂t=D2Δy2+y2t,x+c2d2−a2y2t,x+b2y12t−τ1,x,x∈Ω,t>0,∂y1∂v=∂y2∂v=0,t≥0,x∈∂Ω,y1t,x=φt,x≥0,y2t,x=ψt,x≥0,t,x∈−max1≤i≤2τi,0×Ω¯.For simplicity, we reduce the spatial domainΩ to one dimension 0,π. Although there is an abundant literature about the stability and Hopf bifurcation of delayed differential equation models with diffusion, stability analysis of the interactive behaviors among enterprises in the business ecosystem, especially reaction-diffusion financial systems with delays, is very rare. The main contributions of this study are as follows: first, the interaction between two enterprises is discussed from the perspective of ecology, and the factors such as multiple time delay and diffusion affecting the dynamic properties of the system are introduced into the model, which makes the given model more in line with the actual background in economics than the previous models. Second, the complex dynamic properties of the system are discussed, including the stability of the system and the qualitative analysis of the existence of Hopf bifurcation. The parameter calculation expressions of bifurcation direction, stability of bifurcation periodic solution, and period size are determined. Finally, we explain the results from the perspective of application and supplement the relevant conclusions of the existing literature.The rest of this study is organized as follows. In Section2, we consider the distribution of the roots of the characteristic equation associated linearized system and explore the stability of the positive steady state and the existence of Hopf bifurcation of system (4) to satisfy some requirements. In Section 3, according to the normal form method and the center manifold reduction developed by Hassard et al. [28], we determine the direction of Hopf bifurcation, give an explicit rule, and analyze the stability of bifurcating periodic solutions. In addition, some numerical simulations are conducted in Section 4 to verify our theoretical predictions. The conclusions are presented in Section 5.


## 2. Stability and Hopf Bifurcation

In this section, takingτ as the varying parameter, we analyze the distribution of the roots of the corresponding characteristic equation of (4) at the positive constant steady state and the stability of the positive equilibrium point. Some conditions of Hopf bifurcation of (4) are provided.According to [3, 4], there is a unique positive constant steady state Ey1∗,y2∗ of system (4), which provides that the following condition holds:(5)H1:a22d1>b1d22.Letu1t,x=y1t,x−y1∗,u2t,x=y2t,x−y2∗; then, (4) is equivalent to the following system:(6)∂u1t,x∂t=D1Δu1t,x+u1t,x+y1∗+c1d1−a1u1t,x+y1∗−b1u2t−τ,x+y2∗2,∂u2t,x∂t=D2Δu2t,x+u2t,x+y2∗+c2d2−a2u2t,x+y2∗+b2u1t,x+y1∗2,∂u1∂v=∂u2∂v=0,t≥0,x∈∂Ω,u1t,x=φt,x−y1∗,u2t,x=ψt,x−y2∗,t,x∈−max1≤i≤2τi,0×Ω¯,where τ=τ1+τ2. (6) can be written as follows by setting d1−a1y1∗−b1y2∗2=0 and d2−a2y2∗+b2y1∗2=0.(7)∂u1t,x∂t=D1Δu1t,x+u1t,x+y1∗+c1−a1u1t,x−b1u22t−τ,x−2b1y2∗u2t−τ,x,∂u2t,x∂t=D2Δu2t,x+u2t,x+y2∗+c2−a2u2t,x+b2u12t,x+2b2y1∗u1t,x.We assume thatφ,ψ∈ℓ=C−τ,0,X, and X is defined by(8)X=u1,u2:ui∈W2,20,π,∂uit,x∂v=0,i=1,2,x=0,π,with the inner product ⋅,⋅.Denoteuit=uit,⋅i=1,2 and U=u1,u2T. For ϕ=ϕ1,ϕ2T∈ℓ, system (7) can be rewritten as the following abstract differential equation in the phase space ℓ=C−τ,0,X:(9)U˙t=DΔUt+LUt+FUt,where D=diagD1,D2 and L:ℓ⟶X, F:ℓ⟶X are defined by the following:(10)Lϕ=A10A2B2ϕ10ϕ20+0B100ϕ1−τϕ2−τ,Fϕ=F1+F2F3+F4,where(11)A1=−a1y1∗+c1<0,A2=2b2y1∗y2∗+c2>0,B1=−2b1y2∗y1∗+c1<0,B2=−a2y2∗+c2<0,F1=−a1ϕ120−b1y1∗+c1ϕ22−τ,F2=−2b1y2∗ϕ10ϕ2−τ−b1ϕ10ϕ22−τ,F3=−a2ϕ220+b2y2∗+c2ϕ120,F4=2b2y1∗ϕ20ϕ10+b2ϕ20ϕ120,ϕθ=ϕ1θ,ϕ2θT∈ℓ.It is easy to see that the linearization of (9) is given by(12)U˙t=DΔUt+LUt,and its characteristic equation is(13)λy−DΔy−Leλy=0,where y∈domΔ\0,domΔ⊂X.According to the properties of the Laplacian operator on the bounded domain [14, 18], the operator Δ on X has the eigenvalues −k2 and the relative eigenfunctions βk1=coskx0, βk2=0coskx, for k∈ℕ0=0,1,2,⋯. βk1,βk2 is a basis of the phase space X. So, any element y∈X can be expanded to a Fourier series as follows:(14)y=∑k=0∞YkTβk1βk2,where YkT=y,βk1,y,βk2.By using some simple computations, we have(15)LϕTβk1βk2=LϕTβk1βk2,k∈ℕ0.Hence, (13) is equivalent to(16)∑k=0∞YkTλI2+Dk2−A1B1e−λτA2B2βk1βk2=0.Therefore, all the characteristic roots of (13) are equivalent to the roots of the following equation:(17)λ2+−A1+D1k2−B2+D2k2λ+−A1+D1k2−B2+D2k2−A2B1e−λτ=0.Letp=−A1+B2+D1k2+D2k2>0,q=−A2B1>0, and r=−A1+D1k2−B2+D2k2>0; then, (17) can be rewritten as(18)λ2+pλ+r+qe−λτ=0.It is clear thatλ=0 is not a root of (18) for ∀k∈ℕ0.Whenτ=0, (18) can be converted into(19)λ2+pλ+r+q=0.Supposeλ1 and λ2 are two roots of (19). For ∀k∈ℕ0, it is easy to see that(20)λ1+λ2=−−A1+B2+D1k2+D2k2<0,λ1λ2=−A1+D1k2−B2+D2k2−A2B1>0.Therefore, if (H1) holds, the unique positive equilibrium pointEy1∗,y2∗ of (4) is asymptotically stable for τ=0.Furthermore, suppose thatτ>0 and λ=iωω>0 is a solution of (18). Then, ω satisfies the following equation for some k∈ℕ0:(21)−ω2+iωp+r+qcosωτ−isinωτ=0.By separating the imaginary and real parts of (21), we have(22)ω2−r=qcosωτ,pω=qsinωτ.It follows from (22) that(23)ω4+p2−2rω2+r2−q2=0.Letz=ω2; then, (23) becomes(24)z2+p2−2rz+r2−q2=0.By performing some calculations, one can show that(25)p2−2r=−A1+D1k22+−B2+D2k22,r2−q2=−A1+D1k22−B2+D2k22−A22B12.Note that(26)r+q=−A1+D1k2−B2+D2k2−A2B1>0.Hence, if the following condition holds,(27)H2:r−q=−A1+D1k2−B2+D2k2+A2B1>0,then (24) has no positive real roots, and the characteristic equation has no purely imaginary roots for all τ>0. Therefore, we have the following conclusion if (H1) and (H2) hold.Theorem 1.
Assume that conditions (H1) and (H2) hold; then, the positive equilibrium pointEy1∗,y2∗ of (7) is locally asymptotically stable for all τ≥0.Now, we consider the existence of Hopf bifurcation at the positive equilibrium of (7).Whenk=0, we have(28)r2−q2=A1B2+A2B1A1B2−A2B1,where A1B2−A2B1>0.Hence, ifA1B2+A2B1<0, then r2−q2<0, and (23) has a unique positive root ω+0, where(29)ω+0=122r−p2+p2−2r2−4r2−p21/2,and z0=ω+02. Equation (18) with k=0 has a pair of purely imaginary roots ±iω+0. By (22), we can get the values for τ under this scenario:(30)τj0=1ω+0arctanpω+0ω+02−r+jπω+0,j∈ℕ0.Hence, all roots of (18), except ±iω+0, have no zero real parts. According to Theorem 2.7 by Li et al. [2], we know that if (H1) and A1B2+A2B1<0 hold, then the equilibrium Ey1∗,y2∗ remains stable for τ<τ00=minτj0,j∈ℕ0, and we verify that(31)dReλdττ=τ00>0,where λτ=ατ+iωτ is a root of (18), ατj0=0, and ωτj0=ω0.Apparently, the transversality condition holds in this case, and there is a sign of Hopf bifurcation atω=ω0 and τ=τ00.Furthermore, we discuss the case whenk≥1. Let λ=iω1ω1>0 be a root of (18) with k≥1. Similarly, we have(32)ω14+p2−2rω12+r2−q2=0.By performing some calculations, we can get(33)r2−q2≥−A1+D12−B2+D22−A22B12.It is easy to see that if the following condition holds,(34)H3:−A1+D1−B2+D2>A2B1,then (18) has no purely imaginary roots for k≥1. Thus, the following conclusions are true.Theorem 2.
Assume that condition (H1) holds.ω+0 and τj0j∈ℕ0 be given by (29) and (30), respectively. Then, for system (7), the following conclusions are valid:(i)
When (H2) holds, the positive equilibriumE is locally asymptotically stable for all τ≥0(ii)
WhenA1B2+A2B1<0 and (H3) holds, the positive equilibrium E is locally asymptotically stable if τ∈0,τ00 and unstable if τ>τ00(iii)
Whenτ=τj0j∈ℕ0, (7) occurs as Hopf bifurcation at EBy the same measure, we now make the following assumptions for a certaink0∈ℕ0\0:(35)H4−A1+D1k02−B2+D2k02+A2B1<0,H5−A1+D1k0−12−B2+D2k0−12>A2B1,H6−A1+D1k0+12−B2+D2k0+12>A2B1.Under condition (H4), equation (23) with k=k0 has only a positive solution:(36)ω+k0=222r−p2+2r−p22−4r2−q2k=k01/2.Hence, when (H5) and (H6) hold, fork∈1,2,…,k0−1, we can deduce that(37)−A1+D1k2−B2+D2k2+A2B1≥−A1+D1k0−12−B2+D2k0−12−A2B1>0,and for k∈k0+1,k0+2,⋯,(38)−A1+D1k2−B2+D2k2+A2B1≥−A1+D1k0+12−B2+D2k0+12−A2B1>0.Therefore, it is easy to show that (18) with k∈ℕ0\k0 has no purely imaginary roots, and it follows from (22) that(39)τjk0=1ω+k0arctanpw+k0ω+k02−r+jπω+k0,j∈ℕ0.Lemma 1.
Assume that conditions (H1) and (H4)–(H6) hold; then, (18) with k=k0 has a pair of purely imaginary roots ±iω+k0 with τ=τjk0j∈ℕ0, and all roots of (18), except ±iω+k0, have no zero real parts.Letλτ=ατ+iωτ be a root of (18) near τ=τjk0 satisfying ατjk0=0 and ωτjk0=ω+k0 for j∈ℕ0; then, the following transversality condition is obtained.Lemma 2.
Assume that conditions (H1) and (H4)–(H6) hold. Then,(40)signRedλdττ=τjk0>0.Proof.
Differentiating both sides of (18) with respect to τ yields(41)dλdτ−1=2λ+D1k02+D2k02−A1+B2eλτ−λA2B1−τλ.
It follows that(42)Redλdττ=τjk0−1=Re2λ+D1k02+D2k02−A1+B2eλτ−λA2B1−τλτ=τjk0=Re2iω+k0+D1k02+D2k02−A1+B2eiτjk0ω+k0−iω+k0A2B1=Re−2ω+k0+iD1k02+D2k02−A1+B2eiτjk0ω+k0ω+k0A2B1=−2ω+k0cosτjk0ω+k0−D1k02+D2k02−A1+B2sinτjk0ω+k0ω+k0A2B1.
According to the above analysis,(43)signRedλdττ=τjk0=signRedλdτ−1τ=τjk0>0.
From Lemmas1 and 2, and the Hopf bifurcation theorem [28], we can have the following result.Theorem 3.
Assume that condition (H1) holds.ω+k0 and τjk0 are given by (36) and (39), respectively. Then, for system (7), the following results are true:(i)
Whenτ∈0,τ0k0, the positive equilibrium point E is asymptotically stable(ii)
Whenτ>τ0k0, the positive equilibrium point E is unstable(iii)
Whenτ=τjk0j∈ℕ0, (7) occurs as Hopf bifurcation at E


## 3. Direction of Hopf Bifurcation and Stability of Bifurcating Periodic Solution

In the preceding section, we obtained a series of periodic solutions bifurcating from the positive constant steady stateE of (7) when τ increases past the critical values τjkk=0,k0;j∈ℕ0. In the present section, by using the normal form method and the center manifold theorem presented in [28–30], we compute the directions, stability, and periods of the bifurcating periodic solution. For brevity, we denote τjk by τ⌢ for some fixed j∈ℕ0.Letuit,x=uiτt,xi=1,2 and τ=τ⌢+μ. Then, system (7) can be transformed into(44)∂u1t,x∂t=τ⌢+μD1Δu1t,x+τ⌢+μu1t,x+y1∗+c1−a1u1t,x−b1u22t−1,x−2b1y2∗u2t−1,x,∂u2t,x∂t=τ⌢+μD2Δu2t,x+τ⌢+μu2t,x+y2∗+c2−a2u2t,x+b2u12t,x+2b2y1∗u1t,x.Letuit=uit,⋅i=1,2 and U=u1,u2T. Then, in the phase space ℓ=C−τ,0,X, system (44) can be rewritten as(45)dUtdt=τ⌢DΔUt+Lτ⌢Ut+FUt,μ,where D=diagD1,D2, L:ℓ⟶X, and F:ℓ×R⟶X are given as(46)Lμϕ=μA1ϕ10+B1ϕ2−1A2ϕ10+B2ϕ20,Fϕ,μ=μDΔϕ0+Lμϕ+fϕ,μ,where(47)fϕ,μ=τ⌢+μf1+f2f3+f4,f1=−a1ϕ120−b1y1∗+c1ϕ22−1,f2=−2b1y2∗ϕ10ϕ2−1−b1ϕ10ϕ22−1,f3=F3,f4=F4,for ϕθ=ϕ1θ,ϕ2θT∈ℓ. From the discussion in the last section, we know that (44) has an equilibrium (0,0) and Λ0=−iω+kτ⌢,iω+kτ⌢k=0,k0 are a pair purely imaginary eigenvalues of characteristic equation of the linear equation(48)U˙t=τ⌢DΔUt+Lτ⌢Ut,and the functional differential equation(49)z˙t=Lkτ⌢zt,where Lkτ⌢ is a continuous linear operator mapping C−1,0,ℝ2 into ℝ2.According to the Riesz representation theorem, for the ordinary functional differential equation,(50)X˙t=−τ⌢Dk2Xt+Lτ⌢Xt,there exists a 2×2 matrix function of the bounded variation ηkθ,τ⌢, θ∈−1,0, such that(51)−τ⌢Dk2ϕ0+Lτ⌢ϕ=∫−10dηkθ,τ⌢ϕθ,where ϕθ∈C−1,0,ℝ2. In fact, we can take(52)ηkθ,τ^=τ⌢−D1k2+A10A2−D2k2+B2,θ=0,0,θ∈−1,0,−τ⌢0B100,θ=−1.LetAτ^ be the infinitesimal generator corresponding to (50) and A∗ be the formal adjoint of Aτ^ under the bilinear inner product(53)ψs,ϕθ=ψ0ϕ0−∫−10∫ξ=0θψξ−θdηkθ,τ⌢ϕξdξ=ψ0ϕ0+τ⌢∫−10ψξ+10B100ϕξdξ,for ϕ∈C−1,0,ℝ2, and ψ∈C0,1,ℝ2. Then, Aτ⌢ and A∗ are a pair of adjoint operators [31]. Since ±iω+kτ⌢ are a pair of purely imaginary eigenvalues of Aτ⌢, they are also a pair of eigenvalues of A∗. Let P and P∗ be the center spaces of Aτ⌢ and A∗, respectively, associated with Λ0. Thus, P∗ is the adjoint space of P and dimP=dimP∗=2 [30].Lemma 3.
Let(54)ξ=iω+k+D1k2−A1B1e−iω+kτ⌢,η=−iω+k+D1k2+A1A2.Then,p1θ=eiω+kτ⌢θ1,ξT,p2θ=p1θ¯,θ∈−1,0 is a basis of P associated with Λ0, and q1s=1,ηe−iω+kτ⌢s,q2s=q1s¯,s∈0,1 is a basis of Q associated with Λ0.LetΦ=Φ1,Φ2 and Ψ∗=Ψ1∗,Ψ2∗T with Φ1θ=p1θ+p2θ/2,Φ2θ=p1θ−p2θ/2i, and Ψ1∗s=q1s+q2s/2,Ψ2∗s=q1s−q2s/2i.Through calculation, we can get the following concrete expressions:(55)Φ1θ=Reeiω+kτ⌢θReξeiω+kτ⌢θ=cosω+kτ⌢θM1cosω+kτ⌢θ−M2sinω+kτ⌢θB1,Φ2θ=Imeiω+kτ⌢θImξeiω+kτ⌢=sinω+kω+kτ⌢θM1sinω+kτ⌢θ+M2cosω+kτ⌢θB1,Ψ1∗s=Ree−iω+kτ⌢sReηe−iω+kτ⌢sT=cosω+kτ⌢s−M3A2T,Ψ2∗s=Ime−iω+kτ⌢sImηe−iω+kτ⌢sT=−sinω+kτ⌢s−M4A2T,for θ∈−1,0 and s∈0,1. Also,(56)M1=D1k2−A1cosω+kτ⌢−ω+ksinω+kτ⌢,M2=ω+kcosω+kτ⌢+D1k2−A1sinω+kτ⌢,M3=D1k2+A1cosω+kτ⌢s+ω+ksinω+kτ⌢s,M4=−D1k2+A1sinω+kτ⌢s+ω+kcosω+kτ⌢s.Then, by (53) and some calculations, we can further obtain(57)Ψ1∗,Φ1=1−D1k2+A1M1A2B1+τ⌢sinω+kτ⌢2ω+kτ⌢+cosω+kτ⌢2M1+sinω+kτ⌢2M2,Ψ1∗,Φ2=−D1k2+A1M2A2+τ⌢sinω+kτ⌢2ω+kτ⌢+cosω+kτ⌢2M2−sinω+kτ⌢2M1,Ψ2∗,Φ1=−ω+kM1A2B1−τ⌢−cosω+kτ⌢2−sinω+kτ⌢2ω+kτ⌢M2+sinω+k2M1,Ψ2∗,Φ2=−ω+kM2A2B1−τ^cosω+kτ2−sinω+kτ2ω+kτM1+sinω+kτ⌢2M2.Then, we define(58)Ψ∗,Φ=Ψj∗,Φk=Ψ1∗,Φ1Ψ1∗,Φ2Ψ2∗,Φ1Ψ2∗,Φ2,and structure a new basis Ψ for Q as(59)Ψ=Ψ1,Ψ2T=Ψ∗,Φ−1Ψ∗.Obviously,Ψ,Φ=I2 is a second-order unit matrix. Furthermore, we define fk=βk1,βk2 and c⋅fk=c1βk1+c2βk2, where βk1=coskx0, βk2=0coskx, and c=c1,c2T∈C−1,0,X. Thus, the center subspace of linear equation (48) is given as(60)PCNℓϕ=ΦΨ,ϕ,fk⋅fk,ϕ∈ℓ,where ℓ=PCNℓ⊕Psℓ and Psℓ is the complementary subspace of PCNℓ in ℓ.LetAτ^ be the infinitesimal generator corresponding to (48); then, (44) can be rewritten in the following form:(61)U˙t=Aτ^Ut+Rμ,Ut,where(62)Rμ,Ut=0,θ∈−1,0,FUt,μ,θ=0.Using the decompositionℓ=PCNℓ⊕Psℓ and (60), the solution of (45) on the center manifold can be rewritten as(63)Ut=Φx1tx2t⋅fk+hx1,x2,μ,where x1tx2t=Ψ,Ut,fk,hx1,x2,μ∈Psℓ, h0,0,0=0, and Dh0,0,0=0. In fact, on the center manifold, the solution of (45) is given by(64)Ut=Φx1tx2t⋅fk+hx1,x2,0.Letz=x1−ix2,Ψ0=Ψ10,Ψ20T. By the definitions of Φ1 and Φ2, it is easy to obtain p1=Φ1+iΦ2, and we further obtain(65)Φx1tx2t⋅fk=Φ1,Φ212z+z¯i2z−z¯⋅fk=12p1z+p1¯z¯⋅fk.Therefore, (64) can be expressed as(66)Ut=12p1z+p1¯z¯⋅fk+Wz,z¯,where Wz,z¯=hz+z¯/2,iz−z¯/2,0. By [30], we can deduce that z satisfies(67)z˙=iω+kτ⌢z+gz,z¯,where(68)gz,z¯=Ψ10−iΨ20fUt,0,fk.Let(69)Wz,z¯=W20z22+W11zz¯+W02z¯22+⋯,(70)gz,z¯=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.From (47) and (64), (66) can be rewritten as(71)Ut=12zeiw+kτ⌢θ1ξ+z¯e−iw+kτ⌢θ1ξ⋅fk+Wz,z¯=12zeiw+kτ⌢θ+z¯e−iw+kτ⌢θcoskxzξeiw+kτ⌢θ+ξ¯z¯e−iw+kτ⌢θcoskx+W1z,z¯W2z,z¯.Then,(72)u1t0=12z+z¯coskx+W2010z22+W1110zz¯+W0210z¯22+⋯;u1t−1=12ze−iω+kτ⌢+z¯eiω+kτ⌢coskx+W201−1z22+W111−1zz¯+W021−1z¯22+⋯;u2t0=12zξ+z¯ξ¯coskx+W2020z22+W1120zz¯+W0220z¯22+⋯;u2t−1=12zξe−iω+kτ⌢+z¯ξ¯eiω+kτ⌢coskx+W202−1z22+W112−1zz¯+W022−1z¯22+⋯.Therefore, we can havefUt,0=τ⌢f11f22, where(73)f11=−a1u1t20−b1y1∗+c1u2t2−1−2b1y2∗u1t0u2t−1−b1u1t0u2t2−1,f22=−a2u2t20+b2y2∗+c2u1t20+2b2y1∗u1t0u2t0+b2u1t20u2t0,fUt,0,fk=τ⌢2M11M121π∫0πcos3kxdxz22+τ⌢4M21M221π∫0πcos3kxdxzz¯+τ⌢2M31M321π∫0πcos3kxdxz¯22+τ⌢M41M42z2z¯2+⋯,where(74)M11=−a1−b1y1∗+c1ξ2e−2iω+kτ⌢−2b1y2∗ξe−iω+kτ⌢,M12=−a2ξ2+b2y2∗+c2+2b2y1∗ξ,M21=−2a1−2b1y1∗+c1ξξ¯−2b1y2∗ξe−iω+kτ⌢+ξ¯eiω+kτ⌢,M22=−2a2ξξ¯+2b2y2∗+c2+2b2y1∗ξ+ξ¯,M31=M11¯,M32=M12¯,M41=−a112W2010+W1110coskx,coskx−b1y1∗+c1ξ¯2W202−1eiω+kτ⌢+ξW112−1e−iω+kτ⌢coskx,coskx−2b1y2∗12W112−1+14W202−1+ξ¯4W2010eiω+kτ⌢+ξ4W1110e−iω+kτ⌢coskx,coskx−b1ξξ¯+ξ22e−2iω+kτ⌢coskx,coskx,(75)M42=−a2ξW1120+ξ¯2W2020coskx,coskx+b2y2∗+c2W1120+12W2010coskx,coskx+2b2y1∗12W1120+14W2020+ξ¯4W2010+ξ2W1110coskx,coskx+b2ξ¯+2ξcoskx,coskx,(76)Wijnθ,coskx=1π∫0πWijnθxcoskxdx,i+j=2,n=1,2.Note that∫0πcos3kxdx=0 for ∀k∈ℕ.Letψ1,ψ2=Ψ10−iΨ20. Comparing the coefficient with that of (70), we can obtain(77)g20=0,k∈ℕ,τ⌢2M11ψ1+M12ψ2,k=0,g11=0,k∈ℕ,τ⌢4M21ψ1+M22ψ2,k=0,g02=g20¯,g21=τ⌢M41ψ1+τ⌢M42ψ2.Next, we compute the termsW20θ and W11θ that appear in the equation for g21. Using a method to the one used in [32, 33] and [27], it follows from (69) that(78)W˙z,z¯=W20zz˙+W11z˙z¯+W11zz¯˙+W02z¯z¯˙+⋯,(79)Aτ⌢W=Aτ⌢W20z22+Aτ⌢W11zz¯+Aτ⌢W02z¯22+⋯.In addition, we can get from [30] that(80)W˙=Aτ⌢W+Hz,z¯,where(81)Hz,z¯=H20z22+H11zz¯+H02z¯22+⋯=X0fUt,0−ΦΨ,X0fUt,0,fk⋅fk.Then, from (67) and (70)–(81), we have(82)2iω+kτ⌢−Aτ⌢W20=H20,−Aτ⌢W11=H11,−2iω+kτ⌢−Aτ⌢W02=H02.Since−Aτ⌢ has only two eigenvalues ±iω+kτ⌢, (82) has a unique solution Wij in Q given by(83)W20=2iω+kτ⌢−Aτ⌢−1H20,W11=−Aτ⌢−1H11,W02=−2iω+kτ⌢−Aτ⌢−1H02.From (81), we can obtain, for θ∈−1,0, that(84)Hz,z¯=−ΦθΨ0fUt,0,fk⋅fk=−p1θ+p2θ2,p1θ−p2θ2iΨ10,Ψ20fUt,0,fk⋅fk=−12p1θΨ10−iΨ20+p2θΨ10+iΨ20fUt,0,fk⋅fk=−12p1θg20+p2θg02¯⋅fkz22−12p1θg11+p2θg11¯⋅fkzz¯+⋯.Hence, forθ∈−1,0, we have(85)H20θ=0,k∈ℕ,−12p1θg20+p2θg02¯⋅f0,k=0,(86)H11θ=0,k∈ℕ,−12p1θg11+p2θg11¯⋅f0,k=0,Hz,z¯0=fUt,0−ΦΨ,X0fUt,0,fk⋅fk,(87)H200=τ⌢2M11M12cos2kx,k∈ℕ,τ⌢2M11M12−p10g20+p20g02¯2⋅f0,k=0,(88)H110=τ⌢4M21M22cos2kx,k∈ℕ,τ⌢4M21M22−p10g11+p20g11¯2⋅f0,k=0.From the definition ofAτ⌢ and (82), we can get(89)W˙20θ=2iω+kτ⌢W20θ+12p1θg20+p2θg02¯⋅fk.Note thatp1θ=p10eiω+kτ⌢θ,θ∈−1,0. Therefore,(90)W20θ=12ig20ω+kτ⌢p1θ+ig02¯3ω+kτ⌢p2θ⋅fk+E1e2iω+kτ⌢θwhere(91)E1=W200,k∈ℕ,W200−12ig20ω+kτ⌢p10+ig02¯3ω+kτ⌢p20⋅f0,k=0.By using the definition ofAτ⌢ and combining (82) and (90), we can obtain the following equation:(92)2iω+kτ⌢12ig20ω+0τ⌢p10⋅f0+ig02¯3ω+0τ⌢p20⋅f0+E1−τ⌢DΔ12ig20ω+0τ⌢p10⋅f0+ig02¯3ω+0τ⌢p20⋅f0+E1−Lτ⌢12ig20ω+0τ⌢p1θ⋅f0+ig02¯3ω+0τ⌢p2θ⋅f0+E1e2iω+0τ⌢θ=τ⌢2M11M12−12p10g20+p20g02¯⋅f0.Note that(93)τ⌢DΔp10⋅f0+Lτ⌢p1θ⋅f0=iω+0τ⌢p10⋅f0,τ⌢DΔp20⋅f0+Lτ⌢p2θ⋅f0=iω+0τ⌢p20⋅f0.Then, fork∈ℕ0, we get(94)2iω+kτ⌢E1−τ⌢DΔE1−Lτ⌢E1e2iω+0τ⌢θ=τ⌢2M11M12cos2kx.By (94), we have(95)E1=12EM11M12cos2kx,where(96)E=2iω+k+a1y1∗+c1+D1k2−2b1y2∗y1∗+c1e−2iω+kτ−2b2y1∗y2∗+c22iω+k+a2y2∗+c2+D2k2−1.Similarly, we can also have(97)−W˙11θ=−12p1θg11+p2θg11¯⋅fk,θ∈−1,0,(98)W11θ=12−ig11ω+kτ⌢p1θ+ig11¯ω+kτ⌢p2θ+E2,where(99)E2=14E′M21M22cos2kx,E′=a1y1∗+c1+D1k2−2b1y2∗y1∗+c1−2b2y1∗y2∗+c2a2y2∗+c2+D2k2−1.Substituting (90) and (98) into the expression of g21, we can compute the following values:(100)c10=i2ω+kτ⌢g20g11−2g112−13g022+12g21,μ2=−Rec10Reλ′τjk,β2=2Rec10,T2=−1ω+kτImc10+μ2Imλ′τjk.Through the above analysis, we can draw the following conclusions.Theorem 3.
Suppose conditions (H1)–(H6) hold. Then, the parametersμ2, β2, and T2 determine the properties of Hopf bifurcation, where μ2 determines the direction of the Hopf bifurcation, β2 determines the stability of the bifurcating periodic solutions on the center manifold, and T2 determines the period of the bifurcating periodic solutions, i.e.,(i)
Ifμ2>0μ2<0, the direction of the Hopf bifurcation is forward (backward)(ii)
Ifβ2<0β2>0, the bifurcating periodic solutions are stable (unstable)(iii)
IfT2>0T2<0, the period increases (decreases)


## 4. Numerical Simulations

Example 1.
We chooseD1=D2=1,c1=c2=1.5,d1=1,d2=0.6,a1=0.4,a2=0.9,b1=0.5, and b2=0.2, which satisfy the above conditions. We take the initial conditions as(101)y1t,x=1.2354×1+0.2sin3.7x+0.3sin1.4x−0.6,y2t,x=1.0058×1+0.2sin2.7x+0.3sin0.74x+0.5.
Then, system (4) with the Neumann boundary has a unique positive constant steady state E=1.2354,1.0058. By performing calculations, we can obtain ω+0=0.8838 and τ00=2.3633. According to Theorem 2, E is locally asymptotically stable if τ∈0,τ00 and is unstable if τ>τ00. When τ crosses through the critical point τ00, E loses its stability, and Hopf bifurcation occurs. The simulation results are shown in Figures 1–4. Based on formulas in Section 3, we can perform computations and obtain that Rec10=−2.7368, μ2=1.8499>0, and β2=−5.4735<0. Hence, when τ crosses through τ00, the bifurcation takes place, and a homogeneous, asymptotically stable periodic solution appears near the trivial steady-state solution E, both of which are shown in Figures 3 and 4.Figure 1
The positive constant steady state of the system with Neumann boundary is locally asymptotically stable forτ1=0.9, τ2=0.9, and τ=τ1+τ2<τ00=2.3633. (a) Spatial distribution of y1t,x. (b) Spatial distribution of y2t,x.
(a)(b)Figure 2
The spatial pattern of the system with Neumann boundary forτ1=0.9,τ2=0.9, and τ=τ1+τ2<τ00=2.3633. (a) Spatial pattern of y1t,x. (b) Spatial pattern of y2t,x.
(a)(b)Figure 3
The positive constant steady state of the system with Neumann boundary loses its stability, and Hopf bifurcation occurs forτ1=0.9, τ2=1.5, and τ=τ1+τ2>τ00=2.3633. (a) Spatial distribution of y1t,x. (b) Spatial distribution of y2t,x.
(a)(b)Figure 4
The spatial pattern of the system with Neumann boundary forτ1=0.9,τ2=1.5, and τ=τ1+τ2>τ00=2.3633. (a) Spatial pattern of y1t,x. (b) Spatial pattern of y2t,x.
(a)(b)Example 2.
In this case, we choose diffusion coefficientsD1=D2=0.001, and we keep other parameter values and initial conditions unchanged. We consider Hopf bifurcation of the system with the Dirichlet boundary condition:(102)y1t,0=y1t,π=1.2,y2t,0=y2t,π=tanh1+t3π.
The results of the simulation are shown in Figures5–8. Apparently, when τ passes through the critical value τ00, Hopf bifurcation occurs, and the spatially homogeneous bifurcating periodic solution appears, which suggests that there are periodic oscillations in the interaction between two enterprises.
Figures9 and 10 show the bifurcation state of the system with and without diffusion for y1 and y2 under the Neumann boundary, respectively. It is found that the diffusion can affect the amplitude of oscillation of the periodic solutions to a certain extent, and the bifurcating periodic solutions are regular and orbitally asymptotically stable.
To illustrate the spatial effect of diffusion on the system, we adopt the method used by references [14, 21] and choose different diffusion coefficients for two different boundary systems under the situations of asymptotically stable periodic solutions. The corresponding spatially dynamic behavior diagrams are shown in Figures 11 and 12.Figure 5
The positive constant steady state of the system with Dirichlet boundary is locally asymptotically stable forτ1=0.9, τ2=0.9, and τ=τ1+τ2<τ00=2.3633. (a) Spatial distribution of y1t,x. (b) Spatial distribution of y2t,x.
(a)(b)Figure 6
The spatial pattern of the system with Dirichlet boundary forτ1=0.9,τ2=0.9, and τ=τ1+τ2<τ00=2.3633. (a) Spatial pattern of y1t,x. (b) Spatial pattern of y2t,x.
(a)(b)Figure 7
The positive constant steady state of the system with Dirichlet boundary loses its stability, and Hopf bifurcation occurs forτ1=0.9, τ2=1.5, and τ=τ1+τ2>τ00=2.3633. (a) Spatial distribution of y1t,x. (b) Spatial distribution of y2t,x.
(a)(b)Figure 8
The spatial pattern of the system with Dirichlet boundary forτ1=0.9,τ2=1.5, and τ=τ1+τ2>τ00=2.3633. (a) Spatial pattern of y1t,x. (b) Spatial pattern of y2t,x.
(a)(b)Figure 9
The bifurcation diagram of the system fory1 with respect to τ. (a) Without diffusion. (b) With diffusion.
(a)(b)Figure 10
The bifurcation diagram of the system fory2 with respect to τ. (a) Without diffusion. (b) With diffusion.
(a)(b)Figure 11
The spatial pattern of the system with Neumann boundary fory1 with different Di at times t=350–380 (from top to bottom, Di=0.01,0.001,0.0005,0.0001).Figure 12
The spatial pattern of the system with Dirichlet boundary fory2 with different Di at times t=350–380 (from top to bottom, Di=0.01,0.001,0.0005,0.0001).


## 5. Conclusions

Since the diffusions of information, technology, and products exist objectively, we deem that the model with reaction-diffusion and delays is more suitable for the interaction between two enterprises, which has no relevant results for system (3) at present. The aim of this study is first to introduce multiple delays and diffusion factors into the traditional two-enterprise interaction model and then to explore the asymptotic stability of the equilibrium and the periodicity of the solution by taking the sum of two time delays as the bifurcation parameter. It is shown that when the time delay passes through a critical value, the system changes from a stable state to an unstable state. The numerical simulations demonstrate that the diffusion effect causes rich, spatially dynamic behaviors.We acknowledge that there are some economic implications. In the enterprise cluster, two enterprises can achieve a stable state and symbiosis of development via dynamic interaction. Although there exists competition, cooperation is very effective, and the enterprises can achieve long-term codevelopment through appropriate positive interactions. It is noteworthy that in addition to the varying economic environments, competition powers, mark-saturation degrees, and delay effects, the diffusion speed of products or technologies is an important factor in the self-development of enterprises. When enterprises are in a state of cyclical fluctuations, appropriate diffusion effects may change the amplitude of these fluctuations, which is beneficial to the enterprises to a certain extent. Furthermore, it is easy to find that differences in the initial resources or positions of an enterprise have a certain impact on its future self-development. Therefore, it is obvious that with the development of the digital economy and the emergence of industrial clusters, the internal relationships between enterprises are increasingly in line with the ecological model. As a consequence, it is of great practical significance to study the dynamic properties of these systems for understanding the evolution of market structure and taking corresponding regulatory strategies. Benefiting from the latest research studies, we will further consider the dynamic behaviors of discrete-time and fractional-order models in the future.



---

*Source: 3466954-2022-02-24.xml*