Double Hopf bifurcation in delayed reaction-diffusion systems
aa r X i v : . [ m a t h . D S ] N ov Double Hopf bifurcation in delayed reaction-diffusion systems
Yanfei Du
Shaanxi University of Science and Technology, Xi’an 710021, China.
Ben Niu*, Yuxiao Guo, Junjie Wei
Department of Mathematics, Harbin Institute of Technology, Weihai 264209, China.*Corresponding author, [email protected] (Dated: November 27, 2018)
Abstract
Double Hopf bifurcation analysis can be used to reveal some complicated dynamical behavior ina dynamical system, such as the existence or coexistence of periodic orbits, quasi-periodic orbits,or even chaos. In this paper, an algorithm for deriving the normal form near a codimension-twodouble Hopf bifurcation of a reaction-diffusion system with time delay and Neumann boundarycondition is rigorously established, by employing the center manifold reduction technique and thenormal form method. We find that the dynamical behavior near bifurcation points are provedto be governed by twelve distinct unfolding systems. Two examples are performed to illustrateour results: for a stage-structured epidemic model, we find that double Hopf bifurcation appearswhen varying the diffusion rate and time delay, and two stable spatially inhomogeneous periodicoscillations are proved to coexist near the bifurcation point; in a diffusive predator-prey system,we theoretically proved that quasi-periodic orbits exist on two- or three-torus near a double Hopfbifurcation point, which will break down after slight perturbation, leaving the system a strangeattractor.
Keywords: Reaction-diffusion model, delay, double Hopf bifurcation, coexistence, strange attractor . INTRODUCTION A Hopf bifurcation refers to the phenomenon that steady states lose their stability andgive rise to periodic solutions, as a parameter crosses critical values [2, 28, 32, 42, 52]. Sincethe condition of existence of Hopf bifurcation can be easily verified, and both the directionof Hopf bifurcation and the stability of the bifurcating periodic solution can be determinedby formulas derived from center manifold reduction technique, Hopf bifurcation analysis hasbeen an effective method of the study on the existence of periodic solutions of differentialequations. In recent years, Hopf bifurcation has been intensely studied to investigate the dy-namics of differential equations in various fields [6, 31, 43, 48]. Normal forms theory providespowerful tools to bifurcation analysis, whose basic idea is to employ near-identity nonlin-ear transformations that lead the original system to a qualitatively equivalent differentialequation with the simplest form.In a reaction-diffusion system, especially reaction-diffusion system with delays, there mayexist several Hopf bifurcation critical points with parameters’ varying. Thus, like the casein ordinary differential equations [36], double Hopf bifurcations will be characterized byintersections of two Hopf bifurcation curves in a two-parameter plane. In fact, a doubleHopf bifurcation occurs from a critical point at which the Jacobian evaluated there involvestwo conjugate pairs of pure imaginary eigenvalues. Guckenheimer and Holmes [25] gavethe dynamics in the neighborhood of a codimension-two double Hopf point in an ordinarydifferential equation. When bifurcation parameters are closed to the double Hopf bifurcationpoint, the system may exhibit rich dynamics, such as periodic oscillations, quasi-periodicoscillations, coexisting of several oscillations, two- or three-dimensional invariant torus, andeven chaos [25, 36, 60]. For delayed systems governed by delay differential equations, dueto the fact that delay differential equations are of infinite dimensional, the center manifoldmethod was adopted to reduce delay differential equations into finite-dimensional systems,and then normal forms can be calculated on the center manifold [8, 9, 17, 21, 22, 29].Correspondingly, the multiple timescale method has also been successfully applied to analyzecodimension-two non-resonant double Hopf bifurcations [39, 59]. Research on the doubleHopf bifurcation and complex dynamics of delay differential equations is also carried out bythese methods we mentioned above [7, 10, 40, 61, 62].In recent years, the analysis of Hopf bifurcation in a reaction-diffusion systems has at-2racted much attention. Based on the method provided in [30, 53], Yi et al. [58] carried outHopf bifurcation analysis in a reaction-diffusion system without delays. By decomposing theequation with the orthonormal Fourier basis corresponding to the eigenvalue problem, anexplicit expression is given for several key parameters in determining the properties of bifur-cation according to the method given in [30]. Faria [19] investigated the effect of both timedelay and diffusion on Hopf bifurcation, and improved the method given in [20, 21] to studynormal forms with perturbation parameters of Hopf bifurcation in delayed reaction-diffusionequations with Neumann boundary condition. The results in [19] have been further appliedto equations with various practical backgrounds, such as the predator-prey model [12] and soon. Recently, Hopf bifurcation analysis of delayed reaction-diffusion equations with Dirichletboundary condition has also made considerable progress, such as in [27, 51, 57], the exis-tence of Hopf bifurcation and the stability and direction of bifurcating periodic solutionsfor some population models with Dirichlet boundary condition were considered; Further-more, Guo [26], Chen and Yu [13], combining the non-local phenomenon, considered Hopfbifurcation of the delayed reaction-diffusion equations with Dirichlet boundary condition.Comparing the research on Dirichlet boundary conditions with the research on Neumannboundary conditions, the main difficulty when investigating such models is the existence ofpositive steady-state solutions and the complexity of the corresponding eigenvalue problem,so practical analysis is usually carried out by combining the method of phase space analysis,topological degree theory or Liapunov-Schmidt reduction theory [35]. In this paper, onlyhomogeneous Neumann boundary condition is considered, as we can assume that zero is anequilibrium of the system.For analysis on high codimensional bifurcations of reaction-diffusion equations, so far aswe know, there are very few published results. Only a few results have been conductedmainly on the research of Turing-Hopf bifurcation. Research on interaction between theTuring bifurcation and Hopf bifurcation in reaction-diffusion equations can be traced backto the work of De Wit et al. [14]. Later, Meixner et al. [41] conducted a preliminaryanalysis of the problem through codimension-two bifurcation analysis. Parallel research wasalso seen in literature [5] on the predator-prey system. Song et al. [49] investigated theTuring-Hopf bifurcation from the point of view of analysis of bifurcation and normal forms,in which they extended the method given in [19, 58], and deduced the normal form withparameters of Turing-Hopf bifurcation. Xu and Wei [56] applied this method to study the3uring-Hopf bifurcation of a predator-prey model. Turing-Hopf bifurcation can be seenas a special case of Hopf-zero bifurcation of dynamic system from the point of view ofclassification of bifurcation. Results about Hopf-zero bifurcations and normal form analysisin delayed reaction-diffusion equations can be found in An and Jiang’s preprint [1]. Fordouble Hopf bifurcation, Lewis and Nagata [37] studied the transitions from axisymmetricsteady solutions to nonaxisymmetric waves in a Navier-Stokes model of the differentiallyheated rotating annulus experiment, and an analytical-numerical center manifold reductionis used to analyze the double Hopf bifurcation points that occur at this transition. Besides[37], to our best knowledge, there have not been any results on double Hopf bifurcationanalysis and the corresponding normal form calculation in the reaction-diffusion equationswith or without delays.Based on the general method of normal form simplification in ordinary differential equa-tion or delay differential equation, the ratio of two pairs of imaginary roots appearing atthe double Hopf bifurcation point is important to determine the final simplest normal form.When the ratio is not a rational number, the normal form has a universal form as thoseprovided in [7, 8, 10, 39, 40, 59, 61, 62]. However, if the ratio is rational, there must existsome additional terms which couldn’t be eliminated because of the resonance. Precisely,if we are about to obtain the topological classification of the dynamical behavior near thedouble Hopf points by analyzing the third order normal form, this only require that the ratiois not m : n for m, n = 1 , , ,
4, namely, weak resonance case. For all weakly resonant dou-ble Hopf bifurcation, we can treat them in the same method as that used for non-resonantcases. Otherwise, if the ratio is m : n for m, n = 1 , , ,
4, we say that this is a stronglyresonance case, which has also been investigated in ordinary differential equations and delaydifferential equations [3, 11, 23, 24, 34, 44, 45, 50].In this paper, we aim at the development of a method of computing normal forms oncenter manifold near an equilibrium of a system of a general delayed reaction-diffusionequations with non-resonant or weakly resonant double Hopf singularity. Consider a generalreaction-diffusion model ∂u ( x,t ) ∂t = d ( µ )∆ u ( x, t ) + f ( µ, u t , u t , · · · , u tn ) , ∂u ( x,t ) ∂t = d ( µ )∆ u ( x, t ) + f ( µ, u t , u t , · · · , u tn ) , ... ∂u n ( x,t ) ∂t = d n ( µ )∆ u n ( x, t ) + f n ( µ, u t , u t , · · · , u tn ) , x ∈ (0 , lπ ) , (1)4ith homogeneous Neumann boundary conditions and u ti ( θ )( x ) = u i ( t + θ, x ). Define f i ( µ, , · · · ,
0) = 0, i = 1 , , · · · , n , such that (1) has zero equilibrium, and assume thatevery f i is at least C .We will extend the normal form method given by Faria and Magalh˜aes [18, 21] to investi-gate non-resonant or weakly resonant double Hopf bifurcations in reaction-diffusion systemswith delays. The first step of this method is to decompose the phase space, and then todecompose the equations we investigate. Then we can calculate the second and third ordernormal form on the center manifold, and get the normal form for double Hopf bifurcation,which is a four-dimensional system up to the third order with unfolding parameters re-stricted to the center manifold. By polar coordinates transformation, the four-dimensionalsystem can be reduced to a two-dimensional amplitude system, in which the unfolding pa-rameters can be expressed by the parameters in the original system. Due to the unfoldinganalysis given in [25], we can determine the type of twelve distinct kinds of unfoldings bythe coefficients we calculated in the normal form, and get corresponding bifurcation sets.By detailed analysis of bifurcation sets, all the dynamical behaviors near the double Hopfbifurcation point can be figured out. In this paper, a series of explicit calculation formulas ofthe normal form for non-resonant or weakly resonant double Hopf bifurcation is given, withthree different cases of wave number: n = 0 , n = 0, n = 0 , n = 0, and n = 0 , n = 0.The wave number is related to the spatial profile of bifurcating periodic solutions.To illustrate the calculation process of the normal forms, we perform two examples. Thefirst one is a diffusive, stage-structured epidemic model with two delays. Xiao and Chen [55]proposed an SIS epidemic model with stage structure and time delays in the following form dSdt = αy ( t ) − dS ( t ) − αe − dτ y ( t − τ ) − µS ( t − ω ) I ( t ) + γI ( t ) ,dIdt = µS ( t − ω ) I ( t ) − dI ( t ) − γI ( t ) ,dydt = αe − dτ y ( t − τ ) − βy ( t ) , (2)where S ( t ) and I ( t ) are the population of susceptible and infected immature individuals, y ( t ) represents the population of mature individuals, τ is the maturity delay, and ω is thefreely moving delay. Hopf bifurcation results for such a model are given in Du et al. [15].Adding random diffusion of individuals into (2), we have a diffusive epidemic model, whichis proposed in section 4. We take the time delay ω and a diffusion coefficient as bifurcationparameters. We show that double Hopf bifurcation can appear with two different cases of5ave numbers: n = 0 , n = 0, or n = 0 , n = 0. Thus, the spatio-temporal dynamics turnout to be very complicated near the double Hopf bifurcation point. In some regions, thereare two stable and spatially nonhomogeneous periodic solutions or a homogeneous and anonhomogeneous periodic solution coexisting.The second example we consider is a predator-prey system with time delay. Predator-prey systems have been widely investigated [33, 46, 54]. Here, we choose a simple systemwhich was studied in [48] to illustrate the procedure of normal form calculationd u ( t )d t = u ( t )[ r − a u ( t − τ ) − a v ( t )] , d v ( t )d t = v ( t )[ − r + a u ( t ) − a v ( t )] . (3)Incorporating diffusion terms into (3), the system becomes a reaction-diffusion system, whichis proposed in section 5. Taking r and τ as bifurcation parameters, we will show that such asimple system will exhibit complicated dynamical behavior near the double Hopf bifurcationpoint such as the existence of quasi-periodic solution on a 2-torus and quasi-periodic solutionon a 3-torus. Generally, a vanishing 3-torus might accompany strange attractors, and lead asystem into chaos, which is called the “Ruelle-Takens-Newhouse” scenario [4, 16, 47]. Thusa strange attractor is found near the double Hopf bifurcation point in this predator-preysystem.The paper is organized as follows. In section 2, we discuss eigenfunctions and decomposi-tion of the phase space. Section 3 is devoted to calculating the normal forms of non-resonantor weakly resonant double Hopf bifurcation of (1), and explicit formula of normal form trun-cated to the third order is obtained in three cases of different wave numbers. As applicationsof the method, a diffusive epidemic model with two delays and a diffusive predator-prey sys-tem with a delay are considered in sections 4 and 5, respectively, where the conditions forthe existence of double Hopf bifurcation are obtained, the normal forms are calculated usingthe method and formulas in section 3. Two-parameter unfoldings and bifurcation diagramsnear the critical point are analyzed, and some numerical simulations about the rich dynam-ics near these bifurcation points are demonstrated. Finally, we close the paper with someconclusions. 6 . EIGENFUNCTIONS AND DECOMPOSITION OF THE PHASE SPACE Normal form analysis usually depends on a decomposition of the corresponding phasespace [36, 52]. In case of reaction-diffusion model with time delay, it also requires thedecomposition with respect to spatial variables which can be accomplished by finding theeigenfunctions of Laplacian operator. Thus, we first rewrite system (1) into an abstractordinary differential equation in an appropriate phase space.
Using the general approach to put a partial differential equation into an abstract formintroduced in [53], we define the real-valued Hilbert space X := (cid:26) ( u , u , · · · , u n ) T ∈ ( H (0 , lπ )) n : ∂u i ∂x (0 , t ) = ∂u i ∂x ( lπ, t ) = 0 , i = 1 , , · · · , n (cid:27) . Since we are about to deal with the double Hopf bifurcations with two pairs of purely imag-inary eigenvalues, we usually extend the space X into the corresponding complexificationspace of X in a natural way by X C := X ⊕ iX = { U + iU : U , U ∈ X } with the general complex-value L inner product h u, v i = Z lπ ( u v + u v + · · · + u n v n ) dx, for u = ( u , u , · · · , u n ) T , v = ( v , v , · · · , v n ) T ∈ X C . Let C := C ([ − r, , X C ) denote thephase space with the sup norm. We write u t ∈ C for u t ( θ ) = u ( t + θ ), − r ≤ θ ≤ U ( t )d t = D ( µ )∆ U ( t ) + L ( µ )( U t ) + F ( µ, U t ) , t > , (4)where U ( t ) = u ( t ) u ( t )... u n ( t ) ∈ X C , u t = u t u t ... u tn ∈ C , D ( µ ) = d ( µ ) 0 · · · d ( µ ) · · · · · · d n ( µ ) ,7 ( µ, U t ) = F (1) ( µ, U t ) F (2) ( µ, U t )... F ( n ) ( µ, U t ) , and d i ( µ ) > i = 1 , , · · · , n , µ ∈ R . L : R × C → X C is abounded linear operator. F ( µ, φ ) ∈ X C for ( µ, φ ) ∈ R × C , and F is a C k ( k ≥
3) functionsuch that F ( µ,
0) = 0, D φ F ( µ,
0) = 0, where D φ F ( µ,
0) stands for the Fr´echet derivative of F ( µ, φ ) with respect to φ at φ = 0.Linearizing system (4) at (0 , , · · · , U ( t )d t = D ( µ )∆ U ( t ) + L ( µ )( U t ) . (5)It is well known that the eigenvalue problem∆ ϕ = λϕ, x ∈ (0 , lπ ) , ϕ x | x =0 ,lπ = 0 , has eigenvalues λ m = − m l , m ∈ N = N ∪ { } , with corresponding normalized eigenfunctions γ m ( x ) = cos ml x k cos ml x k L = q lπ , m = 0 , q lπ cos ml x, m ≥ . Let β ( j ) m ( x ) = γ m ( x ) e j , where e j is the j -th unit coordinate vector of R n . Then { β ( j ) m } m ≥ areeigenfunctions of D ( µ )∆ with corresponding eigenvalues − d i ( µ ) m l ( i = 1 , , · · · n ). Applyingthe general theory about elliptic operators, we know { β ( j ) m } m ≥ form an orthnormal basis of X . Define B m the subspace of C by B m := span (cid:8) h v ( · ) , β ( j ) m i β ( j ) m | v ∈ C , j = 1 , , · · · , n (cid:9) . For simplification of notations, in the following we write h v ( · ) , β m i = h v ( · ) , β (1) m ih v ( · ) , β (2) m i ... h v ( · ) , β ( n ) m i . λy − D ( µ )∆ y − L ( µ )( e λ · y ) = 0 , (6)where e λ · ( θ ) y = e λθ y , for θ ∈ [ − r, y = ∞ X m =0 a m a m ... a nm γ m , we find that the characteristic equation (6) is equivalent to a sequence of characteristicequations det (cid:20) λI + m l D ( µ ) − L ( µ )( e λ · I ) (cid:21) = 0 , m ∈ N . (7)In order to consider double Hopf bifurcation, we assume that the following conditions holdfor some µ ∈ R [36]:( H ) There exists a neighborhood O of µ such that for µ ∈ O , (5) has two pairs ofcomplex simple conjugate eigenvalues α ( µ ) ± iω ( µ ) and α ( µ ) ± iω ( µ ), all continuouslydifferential in µ with α ( µ ) = 0 , ω ( µ ) = ω > α ( µ ) = 0 , ω ( µ ) = ω >
0, and theremaining eigenvalues of (5) have non-zero real part for µ ∈ O .( H ) Assume that ω < ω and ω : ω = i : j for i, j ∈ N and 1 ≤ i ≤ j ≤
4, i.e., we donot consider the strongly resonant cases.( H ) The conjugate eigenvalues α k ( µ ) ± iω k ( µ ) are obtained by (7 n k ) and the correspond-ing eigenvalues belong to B n k , k = 1 ,
2, and n , n ∈ N . Assume that n ≤ n .To figure out the spatiotemporal dynamical behavior near the critical point µ = µ , let µ = α + µ , where α ∈ R . Now the system (5) is equivalent tod U ( t )dt = D ∆ U ( t ) + L ( U t ) + e F ( α, U t ) , (8)where D = D ( µ ), L ( · ) = L ( µ )( · ) is a linear operator from C to X C , and e F ( α, ϕ ) =[ D ( α + µ ) − D ]∆ ϕ (0) + [ L ( α + µ ) − L ]( ϕ ) + F ( α + µ , ϕ ).9 .2. Decomposition of the Phase Space In order to adapt the center manifold reduction technique, we have to operate on anequation about U t , which requires an enlarged phase space BC defined by BC := { ψ : [ − r, → X C : ψ is continuous on [ − r, , ∃ lim θ → − ψ ( θ ) ∈ X C } . Eq. (8) can be rewritten as an abstract ordinary differential equation on BC [18]:d U t d t = AU t + X e F ( α, U t ) , (9)where A is the infinitesimal generator of the C -semigroup of solution maps of the linearequation (5), defined by A : C ∩ BC → BC , Aϕ = ˙ ϕ + X [ D ∆ ϕ (0) + L ( ϕ ) − ˙ ϕ (0)] , (10)with dom( A ) = { ϕ ∈ C : ˙ ϕ ∈ C , ϕ (0) ∈ dom(∆) } and X is given by X ( θ ) = , − r ≤ θ < ,I, θ = 0 . Then on B m , the linear equationd U ( t )d t = D ∆ U ( t ) + L ( U t )is equivalent to the retarded functional differential equation on C n :˙ y ( t ) = − m l D y ( t ) + L y t . (11)Define functions of bounded variation η k ∈ BC − r, , R ) for n k ( k = 1 ,
2) mentioned in( H ), such that − n k l D ϕ (0) + L ( ϕ ) = Z − r dη k ( θ ) ϕ ( θ ) , ϕ ∈ C . Let A k ( k = 1 ,
2) denote the infinitesimal generator of the semigroup defined by (11) with µ = µ , m = n , n , and A ∗ k denotes the formal adjoint of A k under the bilinear form( ψ, φ ) k = ψ (0) φ (0) − Z − r Z θ ψ ( ξ − θ ) dη k ( θ ) φ ( ξ ) dξ. Let Φ ( θ ) = ( φ ( θ ) , φ ( θ )) , Φ ( θ ) = ( φ ( θ ) , φ ( θ )) , Ψ ( s ) = ψ ( s ) ψ ( s ) , Ψ ( s ) = ψ ( s ) ψ ( s ) ,
10e the basis of the generalized eigenspace of A k , A ∗ k corresponding to the eigenvalues iω and iω , respectively, and satisfying A k Φ k = Φ k B k , A ∗ k Ψ k = B k Ψ k , (Ψ k , Φ k ) k = I, k = 1 , , with B = diag( iω , − iω ), B = diag( iω , − iω ). Denote Φ( θ ) = (Φ ( θ ) , Φ ( θ )), and Ψ( s ) =(Ψ ( s ) , Ψ ( s )) T . Now, we can decompose BC into a center subspace and its orthocomplement,i.e., BC = P ⊕
Ker π, (12)where π : BC → P is the projection defined by π ( ϕ ) = X k =1 Φ k (Ψ k , h ϕ ( · ) , β n k i ) k · β n k , with β n k = (cid:16) β (1) n k , β (2) n k , · · · β ( n ) n k (cid:17) , h ϕ ( · ) , β n k i = h ϕ ( · ) , β (1) n k ih ϕ ( · ) , β (2) n k i ... h ϕ ( · ) , β ( n ) n k i , c · β n k = c β (1) n k + c β (2) n k + · · · + c n β ( n ) n k for c = ( c , c , · · · , c n ) T ∈ C .Decompose U t ∈ C according to (12) as U t ( θ ) = X k =1 Φ k ( θ ) e z k ( t ) · β n k + y t ( θ ) , (13)where e z k ( t ) = (Ψ k , h U t , β n k i ) k , and y t ∈ C T Ker π := Q for any t . Then in BC the system(9) is equivalent to the system˙ z = Bz + Ψ(0) h e F ( α, P k =1 (Φ k e z k ) · β n k + y ) , β n ih e F ( α, P k =1 (Φ k e z k ) · β n k + y ) , β n i , d y d t = A y + ( I − π ) X e F ( α, P k =1 (Φ k e z k ) · β n k + y ) , (14)where z = ( e z , e z ) T := ( z , z , z , z ) T , B = diag(i ω , − i ω , i ω , − i ω ), and A is the restric-tion of A on Q ⊂ Ker π → Ker π , A ϕ = Aϕ for ϕ ∈ Q .11 . CENTER MANIFOLD REDUCTION AND NORMAL FORM In the previous section, we have calculated the basis of the center subspace and defineda projection onto it. Applying the center manifold theory [38, 53], we know there exists aninvariant local center manifold of the equilibrium, which will be calculated approximatelyand be used to obtain the normal form.
Consider the formal Taylor expansion e F ( α, ϕ ) = X j ≥ j ! e F j ( α, ϕ ) , where e F j is the j th Fr´echet derivation of e F . Then (14) can be written as˙ z = Bz + P j ≥ j ! f j ( z, y, α ) , d y d t = A y + P j ≥ j ! f j ( z, y, α ) , (15)where z = ( z , z , z , z ) T ∈ C , y ∈ Q , and f j = ( f j , f j ) , j ≥
2, are defined by f j ( z, y, α ) = Ψ(0) h e F j ( α, ( P k =1 Φ k e z k ) · β n k + y ) , β n ih e F j ( α, ( P k =1 Φ k e z k ) · β n k + y ) , β n i ,f j = ( I − π ) X e F j ( α, ( P k =1 Φ k e z k ) · β n k + y ) . (16)Let us introduce the following notations as those in [18]: for a normed space Y , V j ( Y )denotes the space of homogeneous polynomials of degree j in 6 variables z = ( z , z , z , z ) T , α = ( α , α ) T with coefficients in Y , V j ( Y ) = X | ( q,l ) | = j c ( q,l ) z q α l : ( q, l ) ∈ N , c ( q,l ) ∈ Y and the norm | P | ( q,l ) | = j c ( q,l ) z q α l | = P | ( q,l ) | = j | c ( q,l ) | Y . Define the operator M j = ( M j , M j ), j ≥ M j : V j ( C ) → V j ( C ) , ( M j p )( z, α ) = D z p ( z, α ) Bz − Bp ( z, α ) ,M j : V j ( Q ) ⊂ V j (Ker π ) → V j (Ker π ) , ( M j h )( z, α ) = D z h ( z, α ) Bz − A h ( z, α ) . (17)12ccording to [18], the normal forms are obtained by a recursive transformations of vari-ables of the form ( z, y, α ) = ( b z, b y, α ) + 1 j ! ( U j ( b z, α ) , U j ( b z, α ) , , (18)with U j = ( U j , U j ) ∈ V j ( C ) × V j ( Q ). This recursive process transforms (15) into thenormal form ˙ z = Bz + P j ≥ j ! g j ( z, y, α ) , dydt = A y + P j ≥ j ! g j ( z, y, α ) , where g j = ( g j , g j ) , j ≥
2, are the new terms of order j , given by g j ( z, y, α ) = f j ( z, y, α ) − M j U j ( z, α ) , and U j ∈ V j ( C ) × V j ( Q ) are given by U j ( z, α ) = ( M j ) − P Im( M j ) × Im( M j ) ◦ f j ( z, , α ) , (19)where P is the projection operator, and f j = ( f j , f j ) denote the terms of order j in ( z, y )obtained after the computation of normal forms up to order j − For simplification of notations, we write z = ( z , z , z , z ) T , Φ( θ ) = (Φ ( θ ) , Φ ( θ )), z x =( z γ n , z γ n , z γ n , z γ n ) T . (13) can be written as U t ( θ )= Φ ( θ ) z ( t ) z ( t ) · β n + Φ ( θ ) z ( t ) z ( t ) · β n + y t ( θ )= Φ ( θ ) z ( t ) z ( t ) γ n + Φ ( θ ) z ( t ) z ( t ) γ n + y t ( θ )= Φ( θ ) z ( t ) γ n z ( t ) γ n z ( t ) γ n z ( t ) γ n + y t ( θ )... y nt ( θ ) = φ ( θ ) z ( t ) γ n + φ ( θ ) z ( t ) γ n + φ ( θ ) z ( t ) γ n + φ ( θ ) z ( t ) γ n + y t ( θ ) △ = Φ( θ ) z x ( t ) + y t ( θ ) . (20)13hen, by (14), (4) is equivalent to the system˙ z = Bz + Ψ(0) h e F ( α, Φ z x + y ) , β n ih e F ( α, Φ z x + y ) , β n i , d y d t = A y + ( I − π ) X e F ( α, Φ z x + y ) . When ω : ω = m : n for m, n ∈ N and 1 ≤ m, n ≤
3, by (17), it is easy to verify that M j ( z p α ι e ξ ) = D z ( z p α ι e ξ ) Bz − Bz p α ι e ξ = ( iω p − iω p + iω p − iω p + ( − ξ iω ) z p α ι e ξ , ξ = 1 , , ( iω p − iω p + iω p − iω p + ( − ξ iω ) z p α ι e ξ , ξ = 3 , , (21)where ξ = 1 , , , e = (1 , , , T , e = (0 , , , T , e = (0 , , , T , e = (0 , , , T , z p = z p z p z p z p , α ι = a ι a ι , p , p , p , p , ι , ι ∈ N , p + p + p + p + ι + ι = j .Therefore,Im( M ) c = span α i z , α i z , α i z , α i z , i = 1 , . (22)Thus, the second order normal form of (4) on the center manifold of the origin near µ = µ has the following form ˙ z = Bz + 12! g ( z, , α ) + · · · , with g ( z, , α ) = Proj Im( M ) c f ( z, , α ).Write the Taylor expansions of L ( µ + α ) and D ( µ + α ) as follows L ( α + µ ) = L + α L (1 , + α L (0 , + 12 (cid:16) α L (2 , + 2 α α L (1 , + α L (0 , (cid:17) + · · · ,D ( α + µ ) = D + α D (1 , + α D (0 , + 12 (cid:16) α D (2 , + 2 α α D (1 , + α D (0 , (cid:17) + · · · . Thus, the second order term of e F ( α, U t ) is e F ( α, U t ) = 2 (cid:16) α D (1 , + α D (0 , (cid:17) ∆ U t (0) + 2 (cid:16) α L (1 , + α L (0 , (cid:17) U t + F ( α, U t ) . By (20), we can write e F ( z, y, α ) = e F ( α, Φ z x + y )= 2 (cid:16) α D (1 , + α D (0 , (cid:17) ∆(Φ(0) z x + y (0))+ 2 (cid:16) α L (1 , + α L (0 , (cid:17) (Φ z x + y ) + F ( α, Φ z x + y ) . (23)14ecalling assumptions F ( µ,
0) = 0, D φ F ( µ,
0) = 0, we have F ( α, Φ z x + y ) = F (0 , Φ z x + y ),which can be expressed by F (0 , Φ z x + y ) = X q + q + q + q =2 F q q q q γ q + q n ( x ) γ q + q n ( x ) z q z q z q z q + S (Φ z x , y ) + o ( | y | ) . (24)Here S (Φ z x , y ) represents the linear terms of y , which can be calculated by D φ F (0 , Φ z x + y ) | y =0 ( y ).From (16), (20), and (23), we have12! f ( z, , α ) = 12! Ψ(0) h e F ( z, , α ) , β n ih e F ( z, , α ) , β n i . According to (22) and the fact that Z lπ γ n dx = Z lπ γ n dx = 1 , we obtain 12! g ( z, , α ) = 12! Proj Im( M ) c f ( z, , α ) = B α z + B α z B α z + B α z B α z + B α z B α z + B α z , (25)with B = ψ (0) (cid:16) − n l D (1 , φ (0) + L (1 , φ (cid:17) ,B = ψ (0) (cid:16) − n l D (0 , φ (0) + L (0 , φ (cid:17) ,B = ψ (0) (cid:16) − n l D (1 , φ (0) + L (1 , φ (cid:17) ,B = ψ (0) (cid:16) − n l D (0 , φ (0) + L (0 , φ (cid:17) . (26) Remark 1.
For many systems with discrete time delays, if there are m discrete delaysin system (4), we can write L ( α, ϕ ) = A ( α ) ϕ (0) + G ( α ) ϕ ( − r ) + G ( α ) ϕ ( − r ) + · · · + G m ( α ) ϕ ( − r m ) , and thus (26) can be calculated explicitly as follows B = ψ (0) (cid:16) − n l D (1 , φ (0) + A (1 , φ (0) + G (1 , φ ( − r ) + · · · + G (1 , m φ ( − r m ) (cid:17) ,B = ψ (0) (cid:16) − n l D (0 , φ (0) + A (0 , φ (0) + G (0 , φ ( − r ) + · · · + G (0 , m φ ( − r m ) (cid:17) ,B = ψ (0) (cid:16) − n l D (1 , φ (0) + A (1 , φ (0) + G (1 , φ ( − r ) + · · · + G (1 , m φ ( − r m ) (cid:17) ,B = ψ (0) (cid:16) − n l D (0 , φ (0) + A (0 , φ (0) + G (0 , φ ( − r ) + · · · + G (0 , m φ ( − r m ) (cid:17) . .3. Third Order Normal Form To find the third order normal form of the double Hopf singularity, neglecting the highorder terms of the perturbation parameters and noticing the assumption ( H ), we haveIm( M ) c = span z z , z z z , z z , z z z , z z , z z z , z z , z z z . According to [18], the normal form up to the third order is given by˙ z = Bz + 12! g ( z, , α ) + 13! g ( z, ,
0) + · · · , where g ( z, ,
0) = Proj
Ker( M ) f ( z, , f ( z, ,
0) = f ( z, ,
0) + [ D z f ( z, , U ( z, D y f ( z, , U ( z, − D z U ( z, g ( z, , , (27)with ( U ( z, α ) , U ( z, α )) ∈ V j ( C ) × V j ( Q ) given by (19). After the calculation of thesecond order normal forms, we obtain the third order f ( z, ,
0) given by (27). Noticing that g ( z, ,
0) = 0, we divide the computation of the third term g ( z, ,
0) = Proj
Im( M ) c f ( z, , Im( M ) c ( f ( z, , , Proj
Im( M ) c ( D z f ( z, , U ( z, , and Proj Im( M ) c ( D y f ( z, , U ( z, , which will be calculated in the following part. Step 1
The calculation of Proj
Im( M ) c ( f ( z, , e F (0 , Φ z x ) as follows e F (0 , Φ z x ) = X q + q + q + q =3 F q q q q γ q + q n ( x ) γ q + q n ( x ) z q z q z q z q ,
16e have f ( z, ,
0) = Ψ(0) h e F (0 , Φ z x ) , β n ih e F (0 , Φ z x ) , β n i = Ψ(0) P q + q + q + q =3 F q q q q R lπ γ q + q +1 n ( x ) γ q + q n ( x ) dxz q z q z q z q P q + q + q + q =3 F q q q q R lπ γ q + q n ( x ) γ q + q +1 n ( x ) dxz q z q z q z q . Thus, 13! Proj
Im( M ) c ( f ( z, , C z z + C z z z C z z + C z z z C z z + C z z z C z z + C z z z , where C = 16 ψ (0) F γ , C = 16 ψ (0) F γ ,C = 16 ψ (0) F γ , C = 16 ψ (0) F γ . (28)Here γ ij = R lπ γ in ( x ) γ jn ( x ) dx , and Z lπ γ n k ( x ) dx = lπ , n k = 0 , lπ , n k = 0 , ( k = 1 , , Z lπ γ n ( x ) γ n ( x ) dx = lπ , n = 0 , n = 0 , lπ , n = 0 , n = 0 , lπ , n = 0 , n = 0 , n = n , lπ , n = n = 0 . Step 2
Proj
Im( M ) c ( D z f ( z, , U ( z, . From (16), (20), (23) and (24), we can write f ( z, ,
0) = ( f , f , f , f ) T as f ( z, ,
0) = Ψ(0) h F (0 , Φ z x ) , β n ih F (0 , Φ z x ) , β n i = Ψ(0) P q + q + q + q =2 F q q q q R lπ γ q + q +1 n ( x ) γ q + q n ( x ) dxz q z q z q z q P q + q + q + q =2 F q q q q R lπ γ q + q n ( x ) γ q + q +1 n ( x ) dxz q z q z q z q , (29)where f = ψ (0)( F z β + F z z β + F z z β + F z z β + F z β + F z z β + F z z β + F z β + F z z β + F z β ) , = ψ (0)( F z β + F z z β + F z z β + F z z β + F z β + F z z β + F z z β + F z β + F z z β + F z β ) ,f = ψ (0)( F z β + F z z β + F z z β + F z z β + F z β + F z z β + F z z β + F z β + F z z β + F z β ) ,f = ψ (0)( F z β + F z z β + F z z β + F z z β + F z β + F z z β + F z z β + F z β + F z z β + F z β ) , and f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (30) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (31) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (32) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (33) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (34) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (35) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (36) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (37) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (38) f = ψ (0) F β , f = ψ (0) F β ,f = ψ (0) F β , f = ψ (0) F β , (39)18ith β ij = R lπ γ in ( x ) γ jn ( x ) dx , and Z lπ γ n k ( x ) dx = q lπ , n k = 0 , , n k = 0 , ( k = 1 , , Z lπ γ n ( x ) γ n ( x ) dx = q lπ , n = 0 , n = 0 , , n = 0 , n = 0 , q lπ , n = 0 , n = 0 , n = 2 n , , n = 0 , n = 0 , n = 2 n , Z lπ γ n ( x ) γ n ( x ) dx = q lπ , n = 0 , n = 0 , q lπ , n = 0 , n = 0 , , n = 0 , n = 0 . Combining with (21) and (29), we can calculate U ( z,
0) = ( U , U , U , U ) T fromthe following formula U ( z,
0) = ( M ) − Proj
Im( M ) f ( z, , . Then, we have U = 1 iω f z + 1 − iω f z z + 1 iω f z z + 1 − iω f z z + 1 − iω f z + 1 − iω + iω f z z + 1 − iω − iω f z z + 12 iω − iω f z + 1 − iω f z z + 1 − iω − iω f z ,U = 13 iω f z + 1 iω f z z + 12 iω + iω f z z + 12 iω − iω f z z + 1 − iω f z + 1 iω f z z + 1 − iω f z z + 12 iω + iω f z + 1 iω f z z + 1 − iω + iω f z ,U = 12 iω − iω f z + 1 − iω f z z + 1 iω f z z + 1 iω − iω f z z + 1 − iω − iω f z + 1 − iω f z z + 1 − iω − iω f z z + 1 iω f z + 1 − iω f z z + 1 − iω f z , = 12 iω + iω f z + 1 iω f z z + 1 iω + 2 iω f z z + 1 iω f z z + 1 − iω + iω f z + 1 − iω + 2 iω f z z + 1 − iω f z z + 13 iω f z + 1 iω f z z + 1 − iω f z . Thus, 13! Proj
Im( M ) c ( D z f ( z, , U ( z, D z z + D z z z D z z + D z z z D z z + D z z z D z z + D z z z , (40)where D = 16 (2 f − iω f + f iω f + f iω f + 2 f iω f + f − iω f + f iω − iω f + f iω f + f iω + iω f ) , (41) D = 16 (2 f − iω f + f − iω f + f iω f + f iω f + f iω − iω f + f iω + iω f + f − iω f + 2 f iω − iω f + f iω f + f iω f + f iω f + 2 f iω + 2 iω f ) , (42) D = 16 ( f − iω f + f iω − iω f + f iω f + f iω + iω f + 2 f − iω f + f iω f + f iω f + 2 f iω f ) , (43) D = 16 (2 f − iω + iω f + f iω f + f − iω f + f iω f + 2 f iω + iω f + f iω f + f − iω f + f iω f + 2 f − iω f + f − iω + 2 iω f + f iω + 2 iω f + f iω f ) . (44) Step 3.
The calculation of Proj
Im( M ) c ( D y f ( z, , U ( z, D y f ( z, ,
0) : Q → X C . From (23) and(24), e F ( z, y,
0) can be written as e F ( z, y,
0) = S (Φ z x , y ) + o ( z , y )= S yz ( y ) z γ n + S yz ( y ) z γ n + S yz ( y ) z γ n + S yz ( y ) z γ n + o ( z , y ) , (45)where S yz i ( i = 1 , , ,
4) are linear operators from Q to X C . Remark 2.
Again if there are m discrete delays in the system (4), we can get the explicitformula S yz i ( ϕ ) = ( F y (0) z i , F y (0) z i , · · · , F y n (0) z i ) ϕ (0)+ ( F y ( − r ) z i , F y ( − r ) z i , · · · , F y n ( − r ) z i ) ϕ ( − r ) + · · · + ( F y ( − r m ) z i , F y ( − r m ) z i , · · · , F y n ( − r m ) z i ) ϕ ( − r m ) . Now, we have D y e F ( z, , ϕ ) = S yz ( ϕ ) z γ n + S yz ( ϕ ) z γ n + S yz ( ϕ ) z γ n + S yz ( ϕ ) z γ n . Let U ( z,
0) = h ( z ) = P j ≥ h j ( z ) · β j ( x ) = P j ≥ h j ( z ) γ j ( x ) with h j ( z ) = h (1) j ( z ) h (2) j ( z )... h ( n ) j ( z ) = X q + q + q + q =2 h (1) jq q q q ( z ) h (2) jq q q q ( z )... h ( n ) jq q q q ( z ) z q z q z q z q . Thus, D y f ( z, , U ( z, h D y e F ( z, , U ( z, , β n ih D y e F ( z, , U ( z, , β n i = Ψ(0) h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z + h S yz ( P j ≥ h j γ j ) γ n , β n i z = Ψ(0) P j ≥ [ b jn n S yz ( h j ) z + b jn n S yz ( h j ) z + b jn n S yz ( h j ) z + b jn n S yz ( h j ) z ] P j ≥ [ b jn n S yz ( h j ) z + b jn n S yz ( h j ) z + b jn n S yz ( h j ) z + b jn n S yz ( h j ) z ] .
21o give a clear expression of our derivation about the normal form, the rest part ofderivation are divided into three cases: n = n = 0, n = 0 , n = 0, and n = 0 , n = 0. Case 1
When n = n = 0, in fact b jn n = b jn n = b jn n = b jn n = Z lπ γ j ( x ) γ ( x ) γ ( x ) dx = √ lπ , j = 0 , , j = 0 . Then, obviously D y f ( z, , U ( z, √ lπ Ψ(0) S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z = √ lπ ψ (0)( S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z ) ψ (0)( S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z ) ψ (0)( S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z ) ψ (0)( S yz ( h ) z + S yz ( h ) z + S yz ( h ) z + S yz ( h ) z ) . Thus, 13! Proj
Im( M ) c ( D y f ( z, , U ( z, E z z + E z z z E z z + E z z z E z z + E z z z E z z + E z z z , where E = √ lπ ψ (0) [ S yz ( h ) + S yz ( h )] ,E = √ lπ ψ (0) [ S yz ( h ) + S yz ( h ) + S yz ( h )] ,E = √ lπ ψ (0) [ S yz ( h ) + S yz ( h )] ,E = √ lπ ψ (0) [ S yz ( h ) + S yz ( h ) + S yz ( h )] . (46)Clearly, we still need to calculate h , h , h , h , h , h , and h .From (17), (10) and (14), we have M U ( z, θ ) = M h ( z )( θ ) = D z h ( z ) Bz − A ( h ( z ))= D z h ( z ) Bz − D ∆ h (0) − L ( h ( z )) , θ = 0 ,D z h ( z ) Bz − D θ h ( z ) , θ = 0 , = P j ≥ [ D z h j ( z ) γ j ( x ) Bz − D ∆ h j (0) γ j ( x ) − L ( h j ( z ) γ j ( x ))] , θ = 0 , P j ≥ [ D z h j ( z ) γ j ( x ) Bz − D θ h j ( z ) γ j ( x )] , θ = 0 . f ( z, , e F ( z, , − φ (0) f ( z, , γ n − φ (0) f ( z, , γ n − φ (0) f ( z, , γ n − φ (0) f ( z, , γ n , θ = 0 , − φ ( θ ) f ( z, , γ n − φ ( θ ) f ( z, , γ n − φ ( θ ) f ( z, , γ n − φ ( θ ) f ( z, , γ n , θ = 0 . From (16), (18), and (19), we have (cid:10) M ( U ( z, , β j (cid:11) = (cid:10) f ( z, , , β j (cid:11) . (47)Matching the coefficients of z q z q z q z q in (47) when j = 0, we can get the results of h , h , h , h , h , h , and h . We take h ( θ ) as an example, and theothers can be calculated in the same method.When θ = 0, solve the following equation2 iω h ( θ ) − ˙ h ( θ ) = − h φ ( θ ) f γ n , β i − h φ ( θ ) f γ n , β i− h φ ( θ ) f γ n , β i − h φ ( θ ) f γ n , β i , and we get h ( θ ) = e iω θ h (0) + 1 − iω ( e iω θ − e iω θ ) φ (0) f + 1 − iω ( e − iω θ − e iω θ ) φ (0) f + 1 − iω + iω ( e iω θ − e iω θ ) φ (0) f + 1 − iω − iω ( e − iω θ − e iω θ ) φ (0) f . When θ = 0,2 iω h (0) − D ∆ h (0) − L ( h ) = h F γ n , β i − φ (0) f h γ n , β i− φ (0) f h γ n , β i − φ (0) f h γ n , β i − φ (0) f h γ n , β i . Combining with[ iω − n l D − L ( e iω · I d )] φ (0) = 0 , [ − iω − n l D − L ( e − iω · I d )] φ (0) = 0 , [ iω − n l D − L ( e iω · I d )] φ (0) = 0 , [ − iω − n l D − L ( e − iω · I d )] φ (0) = 0 ,
23e can get h (0) = 1 − iω φ (0) f + 1 − iω φ (0) f + 1 − iω + iω φ (0) f + 1 − iω − iω φ (0) f − [ − iω + L ( e iω · I d )] − h F γ n , β i , and h ( θ ) = 1 − iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 − iω + iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F . Using the same method mentioned above, we can work out h ( θ ) = 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ [ L ( I d )] − F ,h ( θ ) = 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ [ L ( I d )] − F ,h ( θ ) = 1 iω φ ( θ ) f + 1 − iω + iω φ ( θ ) f + 1 − iω + 2 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) + L ( e ( iω − iω ) · I d )] − F ,h ( θ ) = 1 − iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) + L ( e ( iω + iω ) · I d )] − F ,h ( θ ) = 12 iω − iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 iω φ ( θ ) f + 1 iω − iω φ ( θ ) f − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω ) + L ( e ( − iω + iω ) · I d )] − F , and h ( θ ) = 1 iω − iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f + 1 − iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F . ase 2 When n = 0 , n = 0, we have b jn n = R lπ γ j ( x ) γ ( x ) γ ( x ) dx = √ lπ , j = 0 , , j = 0 ,b jn n = b jn n = R lπ γ j ( x ) γ ( x ) γ n ( x ) dx = √ lπ , j = n , , j = n ,b jn n = R lπ γ j ( x ) γ n ( x ) γ n ( x ) dx = √ lπ , j = 0 , √ lπ , j = 2 n , , otherwise, and D y f ( z, , U ( z, P j ≥ [ b jn n ( S yz ( h j ) z + S yz ( h j ) z ) + b jn n ( S yz ( h j ) z + S yz ( h j ) z )] P j ≥ [ b jn n ( S yz ( h j ) z + S yz ( h j ) z ) + b jn n ( S yz ( h j ) z + S yz ( h j ) z )] , = Ψ(0) √ lπ ( S yz ( h ) z + S yz ( h ) z ) + √ lπ ( S yz ( h n ) z + S yz ( h n ) z ) √ lπ ( S yz ( h n ) z + S yz ( h n ) z ) + √ lπ ( S yz ( h ) z + S yz ( h ) z )+ √ lπ ( S yz ( h n ) z + S yz ( h n ) z ) . Thus, we obtain13! Proj
Ker( M ) c ( D y f ( z, , U ( z, E z z + E z z z E z z + E z z z E z z + E z z z E z z + E z z z , where E = √ lπ ψ (0) [ S yz ( h ) + S yz ( h )] ,E = √ lπ ψ (0) [ S yz h ) + S yz ( h n ) + S yz ( h n )] ,E = ψ (0) h √ lπ ( S yz ( h ) + S yz ( h )) + √ lπ ( S yz ( h n ) + S yz ( h n )) i ,E = ψ (0) h √ lπ ( S yz ( h n ) + S yz ( h n ) + S yz ( h )) + √ lπ S yz ( h n ) i . (48)Clearly, we still need to calculate h , h , h , h n , h n , h , h n , h n , h n , and h n . Using the same method used in case 1, we can get the25ollowing results h ( θ ) = 1 − iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F ,h ( θ ) = 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ [ L ( I d )] − F ,h ( θ ) = 1 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ [ L ( I d )] − F ,h n ( θ ) = 1 − iω + 2 iω φ ( θ ) f + 1 − iω φ ( θ ) f − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) − n l D + L ( e ( iω − iω ) · I d )] − F ,h n ( θ ) = 1 − iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) − n l D + L ( e ( iω + iω ) · I d )] − F ,h ( θ ) = 1 iω − iω φ ( θ ) f + 1 − iω − iω φ ( θ ) f − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F ,h n ( θ ) = 1 iω φ ( θ ) f + 1 iω − iω φ ( θ ) f − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω ) − n l D + L ( e ( − iω + iω ) · I d )] − F ,h n ( θ ) = − √ lπ [ − (2 n ) l D + L ( I d )] − F ,h n ( θ ) = − √ lπ e iω θ [ − iω − (2 n ) l D + L ( e iω · I d )] − F , and h n ( θ ) = 0 . ase 3 When n = 0 , n = 0, from b jn k n k = R lπ γ j ( x ) γ n k ( x ) γ n k ( x ) dx = √ lπ , j = 0 , √ lπ , j = 2 n k , , otherwise, ( k = 1 , ,b jn n = b jn n = R lπ γ j ( x ) γ n ( x ) γ n ( x ) dx = √ lπ , j = n + n , √ lπ , j = n − n = 0 , √ lπ , j = n − n = 0 , , otherwise, we have D y f ( z, , U ( z, P j ≥ [ b jn n ( S yz ( h j ) z + S yz ( h j ) z ) + b jn n ( S yz ( h j ) z + S yz ( h j ) z )] P j ≥ [ b jn n ( S yz ( h j ) z + S yz ( h j ) z ) + b jn n ( S yz ( h j ) z + S yz ( h j ) z )] = Ψ(0) b n n ( S yz ( h ) z + S yz ( h ) z ) + b n n n ( S yz ( h n ) z + S yz ( h n ) z )+ b ( n + n ) n n ( S yz ( h n + n ) z + S yz ( h n + n ) z )+ b n − n n n ( S yz ( h n − n ) z + S yz ( h n − n ) z ) b ( n + n ) n n ( S yz ( h n + n ) z + S yz ( h n + n ) z )+ b ( n − n ) n n ( S yz ( h n − n ) z + S yz ( h n − n ) z )+ b n n ( S yz ( h ) z + S yz ( h ) z ) + b n n n ( S yz ( h n ) z + S yz ( h n ) z ) = Ψ(0) √ lπ ( S yz ( h ) z + S yz ( h ) z ) + √ lπ ( S yz ( h n ) z + S yz ( h n ) z )+ √ lπ ( S yz ( h n + n ) z + S yz ( h n + n ) z )+ b n − n n n ( S yz ( h n − n ) z + S yz ( h n − n ) z ) √ lπ ( S yz ( h n + n ) z + S yz ( h n + n ) z )+ b ( n − n ) n n ( S yz ( h n − n ) z + S yz ( h n − n ) z )+ √ lπ ( S yz ( h ) z + S yz ( h ) z ) + √ lπ ( S yz ( h n ) z + S yz ( h n ) z ) , thus, Proj Im( M ) c D y f ( z, , U ( z, E z z + E z z z E z z + E z z z E z z + E z z z E z z + E z z z , (49)27here E = 16 ψ (0)[ 1 √ lπ ( S yz ( h ) + S yz ( h )) + 1 √ lπ ( S yz ( h n ) + S yz ( h n ))] ,E = 16 ψ (0)[ 1 √ lπ ( S yz ( h ) + 1 √ lπ ( S yz ( h n ) + S yz ( h ( n + n )1001 )+ S yz ( h ( n + n )1010 )) + b ( n − n ) n n ( S yz ( h ( n − n )1001 ) + S yz ( h ( n − n )1010 ))] ,E = 16 ψ (0)[ 1 √ lπ ( S yz ( h ) + S yz ( h )) + 1 √ lπ ( S yz ( h n ) + S yz ( h n ))] ,E = 16 ψ (0)[ 1 √ lπ ( S yz ( h ) + 1 √ lπ ( S yz ( h n ) + S yz ( h ( n + n )0110 )+ S yz ( h ( n + n )1010 )) + b ( n − n ) n n ( S yz ( h ( n − n )0110 ) + S yz ( h ( n − n )1010 ))] . (50)Clearly, we still need to calculate h , h , h n , h n , h , h n , h ( n + n )1001 , h ( n + n )1010 , h ( n − n )1001 , h ( n − n )1010 , h , h n , h n , h n , h ( n + n )0110 , h ( n + n )1010 , h ( n − n )0110 , and h ( n − n )1010 . In fact, we have h ( θ ) = − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F ,h ( θ ) = − √ lπ [ L ( I d )] − F ,h n ( θ ) = − √ lπ e iω θ [ − iω − (2 n ) l D + L ( e iω · I d )] − F , n = 2 n , − iω + iω φ ( θ ) f + − iω − iω φ ( θ ) f − √ lπ [ − iω − (2 n ) l D + L ( e iω · I d )] − F , n = 2 n ,h n ( θ ) = − √ lπ [ − D n ) l + L ( I d )] − F , n = 2 n , − √ lπ [ − D n ) l + L ( I d )] − F + iω φ ( θ ) f + − iω φ ( θ ) f , n = 2 n ,h ( θ ) = − √ lπ [ L ( I d )] − F ,h n ( θ ) = − √ lπ [ − D n ) l + L ( I d )] − F , n = n , iω φ ( θ ) f , + − iω φ ( θ ) f , , n = 2 n , , otherwise,h n + n ( θ ) = − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) − ( n + n ) l D + L ( e ( iω − iω ) · I d )] − F , n + n ( θ ) = − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) − ( n + n ) l D + L ( e ( iω + iω ) · I d )] − F ,h n − n ( θ ) = − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) + L ( e ( iω − iω ) · I d )] − F , n = n , iω φ ( θ ) f + − iω + iω φ ( θ ) f − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) − ( n − n ) l D + L ( e ( iω − iω ) · I d )] − F , n = 2 n , − √ lπ e ( iω − iω ) θ [ − ( iω − iω ) − ( n − n ) l D + L ( e ( iω − iω ) · I d )] − F , otherwise,h n − n ( θ ) = − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) + L ( e ( iω + iω ) · I d )] − F , n = n , − iω φ ( θ ) f + − iω − iω φ ( θ ) f − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) − ( n − n ) l D + L ( e ( iω + iω ) · I d )] − F , n = 2 n , − √ lπ e ( iω + iω ) θ [ − ( iω + iω ) − ( n − n ) l D + L ( e ( iω + iω ) · I d )] − F , otherwise,h ( θ ) = − √ lπ e iω θ [ − iω + L ( e iω · I d )] − F ,h n ( θ ) = − √ lπ [ − (2 n ) l D + L ( I d )] − F ,h n ( θ ) = − √ lπ e iω θ [ − iω − (2 n ) l D + L ( e iω · I d )] − F ,h n ( θ ) = − √ lπ [ − (2 n ) l D + L ( I d )] − F , n = n , , n = n ,h n + n ( θ ) = − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω ) − ( n + n ) l D + L ( e ( − iω + iω ) · I d )] − F , n − n ( θ ) = − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω )+ L ( e ( − iω + iω ) · I d )] − F , n = n , iω − iω φ ( θ ) f + − iω φ ( θ ) f − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω ) − ( n − n ) l D + L ( e ( − iω + iω ) · I d )] − F , n = 2 n , − √ lπ e ( − iω + iω ) θ [ − ( − iω + iω ) − ( n − n ) l D + L ( e ( − iω + iω ) · I d )] − F , otherwise. Hence, by (27), (40), (46), and (49), we have13! g ( z, ,
0) = 13! Proj
Im( M ) c f ( z, ,
0) = B z z + B z z z B z z + B z z z B z z + B z z z B z z + B z z z , (51)with B = C + 32 ( D + E ) , B = C + 32 ( D + E ) ,B = C + 32 ( D + E ) , B = C + 32 ( D + E ) . (52)Therefore, by (25) and (51), the normal form truncated to the third order for double Hopfbifurcation reads as˙ z = iω z + B α z + B α z + B z z + B z z z , ˙ z = − iω z + B α z + B α z + B z z + B z z z , ˙ z = iω z + B α z + B α z + B z z + B z z z , ˙ z = − iω z + B α z + B α z + B z z + B z z z . (53) Remark 3.
To sum up, the whole calculation process above can be accomplished by fol-lowing three steps.(1) Obtain the double Hopf bifurcation point by analyzing the associate characteristicequation, and find n , n . Rewrite the original system into the form as (8), and calculate D , L , D (1 , , D (0 , , L (1 , , and L (0 , . Calculate the eigenfunctions φ i and ψ i ( i = 1 , , , .(2) Calculate B , B , B , and B from (26).(3) Calculate F mnij ( m + n + i + j = 3) , and we can get C , C , C , and C from (28); Calculate F mnij ( m + n + i + j = 2) , and we obtain D , D , D , and D from (41)-(44); Calculate f k ) mnij ( m + n + i + j = 2 , k = 1 , , , by (30)-(39), establish the ABLE I. The twelve unfoldings of system (54).Case Ia Ib II III IVa IVb V VIa VIb VIIa VIIb VIII d +1 +1 +1 +1 +1 +1 –1 –1 –1 –1 –1 –1 b + + + – – – + + + – – – c + + – + – – + – – + + – d − b c + – + + + – – + – + – – linear operators S yz i ( i = 1 , , , , we can get E , E , E , and E by the formulasin step3 in three different cases. Finally, by (52), we can get B , B , B , and B . Make double polar coordinates transformation by z = r cos θ + ir sin θ ,z = r cos θ − ir sin θ ,z = r cos θ + ir sin θ ,z = r cos θ − ir sin θ , where r , r >
0. Define ǫ = Sign(Re B ), ǫ = Sign(Re B ), carry out the rescaling b r = r p | B | , b r = r p | B | , b t = tǫ , and drop the hats, then system (53) becomes˙ r = r ( c + r + b r ) , ˙ r = r ( c + c r + d r ) , (54)where c = ǫ (Re B α + Re B α ) = ǫ (Re B ( µ − µ , ) + Re B ( µ − µ , )) ,c = ǫ (Re B α + Re B α ) = ǫ (Re B ( µ − µ , ) + Re B ( µ − µ , ))) ,b = ǫ ǫ Re B Re B , c = Re B Re B , d = ǫ ǫ . By chapter 7.5 in [25], Eq. (54) has twelve distinct kinds of unfoldings (see Table 1).
Remark 4.
In section 4, case Ib appears, thus we draw bifurcation set and phase portraitsfor the unfolding of case Ib in Figure 1 [25]. For case Ib, near the double Hopf bifurcationpoint, the ( α , α ) plane is divided into six regions D1-D6. In region D1, the equilibrium isa sink; when the parameters vary and enter the region D2 (or D6), the stable equilibriumbifurcates into a stable periodic solution via supercritical Hopf bifurcations. For parameters IG. 1. Phase portraits for the unfoldings of case Ib with ǫ = − in D3 (or D5), periodic solutions appear via secondary supercritical Hopf bifurcations, butthey are saddle type and only stable on the center manifold; when the parameters cross theHopf bifurcation curve c = c c (or c = c /b ) into the region D4, there are two stableperiodic solutions coexisting in D4. Remark 5.
In section 5 case VIa arises, thus we draw bifurcation set and phase portraitsfor the unfolding of case VIa in Figure 2 [25]. For case VIa, near the double Hopf bifurcationpoint, the ( α , α ) plane is divided into eight regions D1-D8. In region D2, the positiveequilibrium is a sink; In region D3, there is a stable periodic solution. In D4, there is aquasi-periodic solution on the two-dimensional torus; In region D5, there is a quasi-periodicsolution on the three-dimensional torus. When the parameters vary and enter D6, three-dimensional torus vanish. Generally, a vanishing 3-torus might accompany the phenomenonof chaos [4, 16, 47], so near the double Hopf bifurcation point, strange attractors may exist. IG. 2. Phase portraits for the unfoldings of case VIa with ǫ = −
1, where L : c = c c , L : c = c − b +1 c + o ( c ), L : c = c − b +1 c , L : c = − b c .
4. APPLICATION TO A DIFFUSIVE EPIDEMIC MODEL WITH DELAY ANDSTAGE STRUCTURE
In this section, a diffusive epidemic model with delay and stage structure is considered.Taking the time delay ω and a diffusive coefficient as bifurcation parameters, we show thatdouble Hopf bifurcation can undergo with two wave numbers in different cases, such as n = 0 , n = 0, and n = 0 , n = 0. Following the steps of Remark 3, the normal form canbe calculated, and the unfolding system can be got. Simulations illustrate that the spatio-temporal dynamics turn out to be very complicated near the double Hopf bifurcation point.In some regions, there are two stable nonhomogeneous periodic solutions or a homogeneousand a nonhomogeneous periodic solution coexisting.33 .1. Model and the Existence of Double Hopf Bifurcation The stage-structured epidemic model with the maturation delay and freely-moving delayis given by ∂S ( x, t ) ∂t = d ∆ S ( x, t ) + αy ( x, t ) − dS ( x, t ) − αe − dτ y ( x, t − τ ) − µS ( x, t − ω ) I ( x, t ) + γI ( x, t ) ,∂I ( x, t ) ∂t = d ∆ I ( x, t ) + µS ( x, t − ω ) I ( x, t ) − dI ( x, t ) − γI ( x, t ) ,∂y ( x, t ) ∂t = d ∆ y ( x, t ) + αe − dτ y ( x, t − τ ) − βy ( x, t ) , x ∈ (0 , lπ ) ,∂S ( x, t ) ∂x = 0 , ∂I ( x, t ) ∂x = 0 , ∂y ( x, t ) ∂x = 0 , at x = 0 and lπ, (55)where S ( x, t ) and I ( x, t ) represent the population of susceptible and infected immatureindividuals, and y ( x, t ) represents the population of mature individuals. τ is the time takenfor the immature individuals to maturity, ω is the time taken for the immature individualsfrom birth to moving freely. α is the natural birth rate, d is the natural death rate of theimmature stage, β is the death rate of the mature individuals of logistic nature, µ is thedisease transmission rate, and γ is the recovery rate.Denote the basic reproduction ratio by R = µα e − dτ (1 − e − dτ ) dβ ( d + γ ) . We can easily verify that if R >
1, system (55) has a positive constant equilibrium E ∗ ( S ∗ , I ∗ , y ∗ ) = ( d + γµ , ( d + γ ) µ ( R − , αe − dτ β ) [15].Linearizing system (55) at the positive equilibrium E ∗ ( S ∗ , I ∗ , y ∗ ), we have ∂∂t S ( x, t ) I ( x, t ) y ( x, t ) = ( D ∆ + A ) S ( x, t ) I ( x, t ) y ( x, t ) + G S ( x, t − ω ) I ( x, t − ω ) y ( x, t − ω ) + G S ( x, t − τ ) I ( x, t − τ ) y ( x, t − τ ) , where D = diag { d , d , d } , A = − d − µS ∗ + γ α µS ∗ − d − γ
00 0 − βy ∗ , G = − µI ∗ µI ∗ , G = − αe − dτ αe − dτ . The corresponding characteristic equation isdet( λI − M n − A − G e − λω − G e − λτ ) = 0 , I is the 3 × M n = − n l D , n ∈ N . That is, each characteristicvalue λ is a root of an equation( λ − αe − dτ e − λτ + 2 βy ∗ + d n l ) · ∆ n ( λ, τ ) = 0 , (56)where∆ n ( λ, τ ) = λ + ( d + d n l + d n l ) λ + d n l ( d + d n l ) + e − λω ( µI ∗ λ + µI ∗ d n l + µI ∗ d ) , (57)with n ∈ N .We can easily prove that the roots of λ − αe − dτ e − λτ + 2 βy ∗ + d n l = 0 have negative realpart. To investigate the location of the roots, it remains to consider the roots of ∆ n ( λ, τ ) = 0.When ω = 0, Eq. (57) becomes λ + T n λ + J n = 0 , n ∈ N . (58)Since T n = d + d n l + d n l + µI ∗ > , J n = d n l ( d + d n l ) + µI ∗ d n l + µI ∗ d >
0, we knowthat all roots of Eq. (58) have negative real part, and thus so do the roots of Eq. (56) for ω = 0 when R > iz ( z >
0) into Eq. (57), and separating the real and imaginary parts,we have − z + d n l ( d + d n l ) = − µI ∗ z sin ωz − ( µI ∗ d n l + µI ∗ d ) cos ωz, ( d + d n l + d n l ) z = − µI ∗ z cos ωz + ( µI ∗ d n l + µI ∗ d ) sin ωz, (59)which is solved bysin z n ω = ( d + d n l + d n l )( µI ∗ d n l + µI ∗ d ) z n − [ d n l ( d + d n l ) − z n ] µI ∗ z n ( µI ∗ z n ) + ( µI ∗ d n l + µI ∗ d ) △ = S n ( z n ) , cos z n ω = − ( d + d n l + d n l ) µI ∗ z n + [ d n l ( d + d n l ) − z n ]( µI ∗ d n l + µI ∗ d )( µI ∗ d n l + µI ∗ d ) + ( µI ∗ z n ) △ = C n ( z n ) . Then, we obtain G ( z ) = z + P n z + Q n = 0 , (60)where P n = ( d n l ) + ( d + d n l + µI ∗ )( d + d n l − µI ∗ ) ,Q n = J n K n = J n ( d n l ( d + d n l ) − µI ∗ d n l − µI ∗ d ) . (61)35oticing that J n >
0, the sign of Q n coincides with that of K n . Since K = − µI ∗ d <
0, and K n is a quadratic polynomial with respect to n , we can conclude that there exists k ∈ N ,such that K n < ≤ n ≤ k ,K n > n ≥ k + 1 , n ∈ N . (62)Denote the positive real root of the equation K n = 0 by n ( k < n < k + 1), then wehave K n = d d l n + ( d − µI ∗ ) d l n − µI ∗ d = 0 . Since − µI ∗ d <
0, we have d d l n +( d − µI ∗ ) d l n = ( d + d n l − µI ∗ ) d l n >
0. It means that d + d n l − µI ∗ > . By (61), we have P n = ( d n l ) + ( d + d n l + µI ∗ )( d + d n l − µI ∗ ) > . (63)Noticing k + 1 > n and by (62) and (63), we get P n > , for n ≥ k + 1 , k ∈ N . (64)From (62) and (64), we can conclude that for n ∈ N and n ≤ k , (60) has only onepositive real root z n = s − P n + p P n − Q n . (65)For n ∈ N and n ≥ k + 1, (60) has no positive real roots.In fact, S n ( z n ) = z n { ( d + d n l + d n l )( µI ∗ d n l + µI ∗ d ) − [ d n l ( d + d n l )] µI ∗ + z n µI ∗ } ( µI ∗ z n ) + ( µI ∗ d n l + µI ∗ d ) , where ( d + d n l + d n l )( µI ∗ d n l + µI ∗ d ) − [ d n l ( d + d n l )] µI ∗ = µI ∗ [ d ( d + d n l + d n l ) +( d n l ) ] > . Thus, when n ∈ { , , · · · , k } , S n ≥
0, define ω jn = arccos C n ( z n ) + 2 jπz n . (66)Differentiating the two sides of Eq. (57) with respective to ω , Using (57) and (59), weobtain dRe λ ( ω )d ω | ω = ω jn = p P n − Q n ( µI ∗ ) z n + ( µI ∗ d n l + µI ∗ d ) > . (67)36rom (66), we have the very first critical value as ω ∗ = ω n = min n ∈{ , , ··· ,k } { ω n } , z ∗ = z n . Due to the general Hopf bifurcation theorem [18, 53], we have the following theorem.
Theorem 1.
Suppose R > .(1) The equilibrium E ∗ of system (55) is locally asymptotically stable for ≤ ω < ω n andis unstable for ω > ω n .(2) System (55) undergoes a Hopf bifurcation at the equilibrium E ∗ when ω = ω jn , for j ∈ N and n ∈ { , , · · · , k } . Now we need to give a condition under which a double Hopf bifurcation occurs. From(66) and (67), we have that when n ∈ { , , · · · , k } , ω n = arccos C n ( z n ) z n , and dRe λ ( ω )d ω | ω = ω n > d and ω , and fix other coefficients in Eq. (55),there are k + 1 Hopf bifurcation curves on the d − ω plane. On every Hopf bifurcation curve ω n , the characteristic equation (56) always has one pair of eigenvalues ± iz n , which crossestransversely the imaginary axis when the parameters cross the Hopf bifurcation curve, andthe rest eigenvalues have non-zero real part. If we can find the intersection of two certainHopf bifurcation curves ω n and ω n , there must exist two pairs of eigenvalues ± iz n and ± iz n at the intersection point, and all the other eigenvalues have non-zero real part. Thus, adouble Hopf singularity can be found by searching for the intersection of the Hopf bifurcationcurves, which can be done by the following process. Firstly, for n , n ∈ { , , · · · , k } , weregard z n and z n as functions of d from Eq. (65). Secondly, the expression of ω n and ω n are obtained from (66). Finally, from ω n = ω n , we can solve the value of d , denoted by d ∗ ,such that ω n = ω n . Thus, we have that when d = d ∗ , ω = ω n = ω n , the Hopf bifurcationcurves ω n and ω n intersect. Thus, system (55) undergoes double Hopf bifurcation at theintersection. Theorem 2.
Suppose that R > and there exists d ∗ , n , n , such that when d = d ∗ , ω n = ω n . Then system (55) undergoes a double Hopf bifurcation at E ∗ when d = d ∗ , ω = ω n = ω n △ = ω ∗ . .2. Normal Form of Double Hopf Bifurcation Let u ( x, t ) = S ( x, ωt ) − S ∗ , u ( x, t ) = I ( x, ωt ) − I ∗ , u ( x, t ) = y ( x, ωt ) − y ∗ . Thus,these transformations not only transform the equilibrium ( S ∗ , I ∗ , y ∗ ) into (0 , , ω to 1, and transform the other delay τ into τ /ω . Denote U ( t ) =( u ( x, t ) , u ( x, t ) , u ( x, t )) T , and ω = ω ∗ + α , d = d ∗ + α , then system (55) can betransformed intod U ( t )d t = D ( ω ∗ + α , d ∗ + α )∆ U ( t ) + L ( ω ∗ + α , d ∗ + α )( U t )+ F ( ω ∗ + α , d ∗ + α , U t ) , (68)where D ( ω ∗ + α , d ∗ + α ) = ( ω ∗ + α ) d d ∗ + α
00 0 d = D + α D (1 , + α D (0 , ,L ( ω ∗ + α , d ∗ + α ) U t = ( L + α L (1 , + α L (0 , ) U t ,F ( ω ∗ + α , d ∗ + α , U t ) = ( ω ∗ + α ) − µU t ( − U t (0) µU t ( − U t (0) − βU t (0) , with D = ω ∗ d d ∗
00 0 d , D (1 , = d d ∗
00 0 d , D (0 , = ω ∗ ,L U t = ω ∗ (cid:2) AU t (0) + G U t ( −
1) + G U t ( − τ ∗ ) (cid:3) ,L (1 , U t = (cid:2) AU t (0) + G U t ( −
1) + G U t ( − τ ∗ ) (cid:3) ,L (0 , U t = 0 , and τ ∗ = τ /ω. Eq. (68) can be rewritten as dUdt = D ∆ U ( t ) + L U t + e F ( α, U t ) , (69)38here e F ( α, U t ) = α D (1 , ∆ U + α D (0 , ∆ U + α L (1 , U t + F ( ω ∗ + α , d ∗ + α , U t ) . (70)Consider the linearized system of (69) dUdt = D ∆ U ( t ) + L U t . (71)From the previous discussion, we know that system (71) has pure imaginary eigenvalues {± iz n ω ∗ , ± iz n ω ∗ } at the double Hopf bifurcation point and the other eigenvalues withnon-zero real part.Assume that the non-resonant condition holds true and use the algorithm we give inSection 3. After a few calculations, we have that the bases of P Λ and P ∗ , respectively, are Φ ( θ ) = ( φ ( θ ) , ¯ φ ( θ ) , φ ( θ ) , ¯ φ ( θ )), and Ψ ( s ) = ( ψ ( s ) , ¯ ψ ( s ) , ψ ( s ) , ¯ ψ ( s )) T , with φ ( θ ) = (1 , p , p ) T e iz n ω ∗ θ , φ ( θ ) = (1 , p , p ) T e iz n ω ∗ θ ,ψ ∗ ( s ) = D (1 , q ∗ , q ∗ ) e − iz n ω ∗ s , ψ ∗ ( s ) = D (1 , q ∗ , q ∗ ) e − iz n ω ∗ s , (72)where p = µI ∗ e − iz n ω ∗ d n l + iz n , p = 0 , p = µI ∗ e − iz n ω ∗ d n l + iz n , p = 0 ,q ∗ = − − µS ∗ + γ − d n l − iz n , q ∗ = − − µS ∗ + γ − d n l − iz n ,q ∗ = − α − αe − dτ e − iz n ω ∗ τ ∗ − βy ∗ − d n l + αe − dτ e − iz n ω ∗ τ ∗ − iz n ,q ∗ = − α − αe − dτ e − iz n ω ∗ τ ∗ − βy ∗ − d n l + αe − dτ e − iz n ω ∗ τ ∗ − iz n ,D = q ∗ p + ω ∗ µI ∗ e − izn ω ∗ ( q ∗ − , D = q ∗ p + ω ∗ µI ∗ e − izn ω ∗ ( q ∗ − . Consider the Taylor expansion e F ( α, U t ) = 12! e F ( α, U t ) + 13! e F ( α, U t ) , where e F ( α, U t ) = 2 α D (1 , ∆ U + 2 α D (0 , ∆ U + 2 α L (1 , U t + F ( ω ∗ + α , d ∗ + α , U t ), and39 e F ( α, U t ) = 0. By a few calculations, we have F = 2 ω ∗ − µe − iz n ω ∗ p µe − iz n ω ∗ p − βp ,F = 2 ω ∗ − µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µ ( e − iz n ω ∗ p + e − iz n ω ∗ p ) µ ( e − iz n ω ∗ p + e − iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µe iz n ω ∗ p µe iz n ω ∗ p − βp ,F = 2 ω ∗ − µ ( e iz n ω ∗ p + e − iz n ω ∗ p ) µ ( e iz n ω ∗ p + e − iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µ ( e iz n ω ∗ p + e iz n ω ∗ p ) µ ( e iz n ω ∗ p + e iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µe − iz n ω ∗ p µe − iz n ω ∗ p − βp , = 2 ω ∗ − µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) µ ( e − iz n ω ∗ p + e iz n ω ∗ p ) − βp p ,F = 2 ω ∗ − µe iz n ω ∗ p µe iz n ω ∗ p − βp . By Remark 2 in Section 3, define the linear operators S yz ( y ) = F y (0) z y (0) + F y ( − z y ( −
1) + F y ( − τω ∗ ) z y ( − τω ∗ ) ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( −
1) + F y ( − τω ∗ ) z y ( − τω ∗ ) ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( −
1) + F y ( − τω ∗ ) z y ( − τω ∗ ) ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( −
1) + F y ( − τω ∗ ) z y ( − τω ∗ ) , where F y (0) z = 2 ω ∗ − µe − iz n ω ∗ µe − iz n ω ∗
00 0 − βp , F y ( − z = 2 ω ∗ − µp µp , F y ( − τω ∗ ) z = 0 ,F y (0) z = 2 ω ∗ − µe iz n ω ∗ µe iz n ω ∗
00 0 − βp , F y ( − z = 2 ω ∗ − µp µp , F y ( − τω ∗ ) z = 0 ,F y (0) z = 2 ω ∗ − µe − iz n ω ∗ µe − iz n ω ∗
00 0 − βp , F y ( − z = 2 ω ∗ − µp µp , F y ( − τω ∗ ) z = 0 ,F y (0) z = 2 ω ∗ − µe iz n ω ∗ µe iz n ω ∗
00 0 − βp , F y ( − z = 2 ω ∗ − µp µp , F y ( − τω ∗ ) z = 0 . Following the steps in Remark 3, we can get all the coefficients in (53), and then obtain thenormal forms up to third order. 41 .3. Simulations
In the following, we take α = 2 . d = 0 . µ = 0 . γ = 0 . β = 0 . τ = 1, d = 0 . d = 0 .
06, and let ω and d be the bifurcation parameters.By (66), we can draw the curves of Hopf bifurcation values when d varies. As shownclearly in Fig. 3 a), every two Hopf bifurcation curves intersect at a double Hopf bifurcationpoint such as HH1, HH2 and HH3. When d = 5 . ω intersects ω , and we denotethe double Hopf bifurcation point HH2. It means that the wave numbers are the case of n = 1 , n = 2. For HH2, we have z n = z = 2 . z n = z = 3 . ω = ω = 0 . B = 0 . . i , B =0 . . i , B = 0 . . i , B = − . . i , B = − . − . i , B = − . − . i , B = − . . i , and B = − . − . i .Thus, ǫ = − ǫ = − b = 1 . c = 0 . d = 1, d − b c = − . c = − . α − . α , c = − . α + 0 . α . Thus, the bifurcation set near HH2 is shown in Figure 3 b), in whichthe two black lines are two pitchfork bifurcation curves ω = ( d − . / (28 . . ω = ( d − . / ( − . . ǫ . Thus, when ν = ν ( ǫ ) and ǫ are near 0, (i.e. when ω is near ω ∗ ), two periodic solutions have the following representations respectively U t ( ν, θ )( x ) = ǫReφ ( θ ) e iz n ω ∗ t cos n l x + O ( ǫ ) ,U t ( ν, θ )( x ) = ǫReφ ( θ ) e iz n ω ∗ t cos n l x + O ( ǫ ) , where φ ( θ ) and φ ( θ ) is defined in (72). Figure 4 illustrates that the two solutions havetotally different spatial shape.For double Hopf bifurcation points HH1 and HH3, we can calculate the parameters inTable 2. Obviously, these two points are also of type Ib. We notice that during the calcu-lation of step 3 in third order normal form of the three double Hopf bifurcation points, thecases n = 0 , n = 0 (e.g. HH3), and n = 0 , n = 0 (e.g. HH1, HH2) both occur.42 ) d ω n ω ω ω ω ω b) d ω FIG. 3. a) When α = 2 . d = 0 . µ = 0 . γ = 0 . β = 0 . τ = 1, d = 0 . d = 0 .
06, thebifurcation set on the d − ω plane is drawn, and double Hopf bifurcation points are marked. b)Complete bifurcation set near the double Hopf point HH2.TABLE II. Parameter values at double Hopf bifurcation points. P oint d ∗ ω ∗ n z n n z n b c d d − b c HH1 1 .
62 0 .
5. APPLICATION TO A DIFFUSIVE PREDATOR-PREY MODEL WITH DELAY
In this section, a diffusive predator-prey model with delay is considered. Taking r and τ as bifurcation parameters, we can find the double Hopf bifurcation point. By the steps inRemark 3, we can obtain the normal form and get the unfolding system. We will show thatthe system will exhibit complex dynamical behavior near the double Hopf bifurcation point:the existence of quasi-periodic solution on a 2-torus, quasi-periodic solution on a 3-torus,and even chaos. 43 .1. Model and the Existence of Double Hopf Bifurcation Let us consider the predator-prey system ∂X ( x, t ) ∂t = d ∆ X ( x, t ) + X ( x, t )[ r − a X ( x, t − τ ) − a Y ( x, t )] ,∂Y ( x, t ) ∂t = d ∆ Y ( x, t ) + Y ( x, t )[ − r + a X ( x, t ) − a Y ( x, t )] , x ∈ (0 , lπ ) , (73)equipping with homogeneous Neumann boundary condition, where X ( x, t ) and Y ( x, t ) rep-resent the densities of prey and predator populations at time t and location x , respectively, d , d > τ is the generation time ofthe prey species, and r i , a ij ( i, j = 1 ,
2) are positive constants.It is obvious that if r a − r a >
0, (73) has a unique positive equilibrium E ∗ = ( X ∗ , Y ∗ )[48] where X ∗ = r a + r a a a + a a , Y ∗ = r a − r a a a + a a . By the tanslation u = X − X ∗ , v = Y − Y ∗ , (73) can be written as ∂u ( x, t ) ∂t = d ∆ u ( x, t ) + ( u ( x, t ) + X ∗ )[ − a u ( x, t − τ ) − a v ( x, t )] ,∂v ( x, t ) ∂t = d ∆ v ( x, t ) + ( v ( x, t ) + Y ∗ )[ a u ( x, t ) − a v ( x, t )] . (74)The linearization of (74) at the origin is ∂u ( x, t ) ∂t = d ∆ u ( x, t ) − a X ∗ u ( x, t − τ ) − a X ∗ v ( x, t )] ,∂v ( x, t ) ∂t = d ∆ v ( x, t ) + a Y ∗ u ( x, t ) − a Y ∗ v ( x, t ) , whose characteristic equation is λ + A n λ + B n + e − λτ ( Cλ + D n ) = 0 , (75)with n ∈ N , A n = d n l + d n l + a Y ∗ ,B n = d n l ( a Y ∗ + d n l ) + a a X ∗ Y ∗ ,C = a X ∗ ,D n = ( a Y ∗ + d n l ) a X ∗ . Clearly, λ = 0 is not a root of (75), which excludes the existence of Turing bifurcation.When τ = 0, Eq. (75) becomes the following sequence of quadratic polynomial equations λ + ( A n + C ) λ + ( B n + D n ) = 0 , n ∈ N , (76)44here A n + C = d n l + d n l + a Y ∗ + a X ∗ > , B n + D n = d n l ( a Y ∗ + d n l ) + a a X ∗ Y ∗ + ( a Y ∗ + d n l ) a X ∗ > . We know that all roots of Eq. (76) have negativereal part.We would like to seek critical values of τ such that there exists a pair of imaginaryeigenvalues. Let ± iω ( ω >
0) be the solution of Eq. (75), then we have − ω + B n = − Cω sin ωτ − D n cos ωτ,A n ω = − Cω cos ωτ + D n sin ωτ, which is solved by sin ωτ = A n ωD n − ( B n − ω ) Cω ( Cω ) + D n = S n ( ω ) , cos ωτ = − A n Cω + ( B n − ω ) D n D n + ( Cω ) = C n ( ω ) . Then, we obtian G ( ω ) = ω + ( A n − B n − C ) ω + B n − D n = 0 . (77)Suppose that( H ) A n − B n − C < , B n − D n > , and ( A n − B n − C ) − B n − D n ) > ω ± n = s − ( A n − B n − C ) ± p ( A n − B n − C ) − B n − D n )2 .τ j ± n = arccos C n ( ω )+2 jπω ± n , if S n ( ω ± n ) > , π − arccos C n ( ω )+2 jπω ± n , if S n ( ω ± n ) < . (78)By calculation and the results in [48], we can verify the transversality conditionSignRe dλdτ | λ = iω + n > , and SignRe dλdτ | λ = iω − n < . Theorem 3.
Suppose that ( H ) holds, system (73) undergoes a Hopf bifurcation at theorigin when τ = τ j − n or τ j + n . We fix the other parameters, and choose r and τ as double Hopf bifurcation parameters. Theorem 4.
Suppose that ( H ) holds, and there exists r ∗ such that τ j − n = τ j + n . Thensystem (73) undergoes a double Hopf bifurcation at the origin when r = r ∗ , τ = τ j − n = τ j + n △ = τ ∗ . .2. Normal Form of Double Hopf Bifurcation Let u ( x, t ) = X ( x, τ t ) − X ∗ , v ( x, t ) = Y ( x, τ t ) − Y ∗ , τ = τ ∗ + α , r = r ∗ + α , then (73)can be written as ∂u ( x, t ) ∂t = τ ∗ d ∆ u ( x, t ) + τ ∗ ( u ( x, t ) + X ∗ )[ − a u ( x, t − − a v ( x, t )] ,∂v ( x, t ) ∂t = τ ∗ d ∆ v ( x, t ) + τ ∗ ( v ( x, t ) + Y ∗ )[ a u ( x, t ) − a v ( x, t )] . Denote U ( t ) = ( u ( x, t ) , v ( x, t ) T , then we haved U d t = D ( τ ∗ + α , r ∗ + α )∆ U ( t ) + L ( τ ∗ + α , r ∗ + α ) U t + F ( τ ∗ + α , r ∗ + α , U t ) . (79)Here D ( τ ∗ + α , r ∗ + α ) = ( τ ∗ + α ) d d = D + α D (1 , + α D (0 , ,L ( τ ∗ + α , r ∗ + α ) U t = ( L + α L (1 , + α L (0 , ) U t ,F ( ω ∗ + α , r ∗ + α , U t ) = ( τ ∗ + α ) − a U t ( − U t (0) − a U t (0) U t (0) a U t (0) U t (0) − a U t , with D = τ ∗ d d , D (1 , = d d , D (0 , = 0 ,L U t = τ ∗ ( AU t (0) + G U t ( − , L (1 , U t = AU t (0) + G U t ( − ,L (0 , U t = τ ∗ a a + a a − a a a a − a a U t (0) + − a a
00 0 U t ( − , where A = − a X ∗ r ∗ a Y ∗ r ∗ − a Y ∗ r ∗ , G = − a X ∗ r ∗
00 0 , and X ∗ r ∗ = r ∗ a + r a a a + a a , Y ∗ r ∗ = r ∗ a − r a a a + a a .Eq. (79) can be rewritten asd U d t = D ∆ U ( t ) + L U t + e F ( α, U t ) , (80)where e F ( α, U t ) = ( α D (1 , + α D (0 , )∆ U + ( α L (1 , + α L (0 , ) U t + F ( τ ∗ + α , r ∗ + α , U t ) . U d t = D ∆ U ( t ) + L U t . (81)From the previous discussion, we know that system (81) has pure imaginary eigenvalues ± iω + n τ ∗ , ± iω − n τ ∗ at the double Hopf bifurcation point and the other eigenvalues with non-zero real part. After a few calculations, we have that the bases of P and P ∗ , respectively,are Φ ( θ ) = ( φ ( θ ) , ¯ φ ( θ ) , φ ( θ ) , ¯ φ ( θ )), and Ψ ( s ) = ( ψ ( s ) , ¯ ψ ( s ) , ψ ( s ) , ¯ ψ ( s )) T , with φ ( θ ) = (1 , p ) T e iω − n τ ∗ θ , φ ( θ ) = (1 , p ) T e iω + n τ ∗ θ ,ψ ∗ ( s ) = D (1 , q ) e − iω − n τ ∗ s , ψ ∗ ( s ) = D (1 , q ) e − iω + n τ ∗ s , where p = − a X ∗ r e − iω − n τ ∗ − iω − n − d n l a X ∗ r ∗ , p = − a X ∗ r e − iω + n τ ∗ − iω + n − d n l a X ∗ r ∗ ,q = − − a X ∗ r e − iω − n τ ∗ − iω − n − d n l a Y ∗ r ∗ , q = − − a X ∗ r e − iω + n τ ∗ − iω + n − d n l a Y ∗ r ∗ ,D = 11 + p q − a X ∗ r ∗ τ ∗ e − iω − n τ ∗ , D = 11 + p q − a X ∗ r ∗ τ ∗ e − iω + n τ ∗ . Consider the Taylor expansion e F ( α, U t ) = 12! e F ( α, U t ) + 13! e F ( α, U t ) , where e F ( α, U t ) = 2 α D (1 , ∆ U + 2 α D (0 , ∆ U + (2 α L (1 , + 2 α L (0 , ) U t + 2 F ( ω ∗ + α , d ∗ + α , U t ), and e F ( α, U t ) = 0.By a few calculations, we have F = 2 τ ∗ − a e − iω − n τ ∗ − a p a p − a p ,F = 2 τ ∗ − a ( e iω − n τ ∗ + e − iω − n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p ,F = 2 τ ∗ − a ( e − iω + n τ ∗ + e − iω − n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p ,F = 2 τ ∗ − a ( e iω + n τ ∗ + e − iω − n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p , = 2 τ ∗ − a e iω − n τ ∗ − a p a p − a p ,F = 2 τ ∗ − a ( e − iω + n τ ∗ + e iω − n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p ,F = 2 τ ∗ − a ( e iω + n τ ∗ + e iω − n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p ,F = 2 τ ∗ − a e − iω + n τ ∗ − a p a p − a p ,F = 2 τ ∗ − a ( e iω + n τ ∗ + e − iω + n τ ∗ ) − a ( p + p ) a ( p + p ) − a p p ,F = 2 τ ∗ − a e iω + n τ ∗ − a p a p − a p . By Remark 2 in Section 3, establish the linear operators S yz ( y ) = F y (0) z y (0) + F y ( − z y ( − ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( − ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( − ,S yz ( y ) = F y (0) z y (0) + F y ( − z y ( − , where F y (0) z = 2 τ ∗ − a e − iω − n τ ∗ − a p − a a p a − a p , F y ( − z = 2 τ ∗ − a
00 0 ,F y (0) z = 2 τ ∗ − a e iω − n τ ∗ − a p − a a p a − a p , F y ( − z = 2 τ ∗ − a
00 0 ,F y (0) z = 2 τ ∗ − a e − iω + n τ ∗ − a p − a a p a − a p , F y ( − z = 2 τ ∗ − a
00 0 ,F y (0) z = 2 τ ∗ − a e iω + n τ ∗ − a p − a a p a − a p , F y ( − z = 2 τ ∗ − a
00 0 . Following the steps in Remark 3, we can get all the coefficients in (53), and then obtain thenormal forms up to the third order. 48 .3. Simulations
In this section, we fix r = 1, a = 1, a = 1 . a = 2 . a = 1, d = 0 . d = 0 . l = 3, and take r and τ as bifurcation parameters. We can find that when 0 . < r < . H ) holds for all n = 0, and when 0 . < r < .
19, ( H ) holds for all n = 1. By (78), wecan draw the curves of Hopf bifurcation values when r varies, which is shown in Figure 5(a).When r = 0 . τ − intersects τ , and we denote the double Hopf bifurcationpoint HH. For HH, we have ω +0 = 0 . ω − = 0 . τ − = τ = 10 . B = 0 . . i , B = 1 . . i , B = − . . i , B = − . . i , B = − . − . i , B = − . . i , B = 3 . − . i , B = 0 . . i , ǫ = − ǫ = 1, b = 2 . c = − . d = −
1, and d − b c = 0 . X (0 , t ) , Y (0 , t )), i.e., x = 0. Moreover, because the periodic solutions oscillate in an infinitedimensional phase space [53], we give simulations near the double Hopf bifurcation point onthe Poincar´e section X (0 , t − τ ) = X ∗ . We can see that the system exhibits rich dynamicalbehavior near the bifurcation point HH in Figure 6. System (73) has a quasi-periodic so-lution on a 2-torus which becomes a quasi-periodic solution on a 3-torus and then vanishesthrough a heteroclinic orbit, which is shown in Figure 6 (a) and (b) respectively. After thevanishing of 3-torus, system (73) has a strange attractor, and exhibits chaotic behavior,which is shown in 6 (c). 49 . CONCLUDING REMARKS In this paper, we have extended the center manifold reduction and normal form methodto analyze the dynamical behavior near the double Hopf bifurcation point in a delayedreaction-diffusion system. The method appears to be quite complicated as proceeded in thispaper, but it is still an explicit algorithm and is not difficult to be implanted into a computerprogram. What we should care more about in this method is to calculate the double Hopfpoint and all the rest part can be exactly obtained by following the procedure we have given.After normal form derivation with respect to a delayed reaction-diffusion system, thedynamics near the bifurcation point, like the case in ordinary differential equations, is alsogoverned by twelve distinct kinds of unfolding systems, and the bifurcation set about or-dinary differential equations for each of the twelve types of unfoldings still applies to ourmodel. To show that our method is a powerful method for analyzing local behavior arounda double Hopf bifurcation point, we give two examples: in the stage-structured epidemicmodel, we show that two spatially inhomogeneous periodic oscillations coexist near the sin-gular point; in the predator-prey system, we show the existence of quasi-periodic solutionon a 2-torus, quasi-periodic solution on a 3-torus, and strange attractor appears near thebifurcation point.
ACKNOWLEDGMENTS
The authors are grateful to the handling editor and anonymous referees for their carefulreading of the manuscript and valuable comments, which improve the exposition of the papervery much. This research is supported by National Natural Science Foundation of China(11701120, 11771109) and Shaanxi Provincial Education Department grant (18JK0123). [1] An, Q., Jiang, W: Spatiotemporal attractors generated by the Turing-Hopf bifurcation ina time-delayed reaction-diffusion system. Disctete Cont. Dyn-B doi: 10.3934/dcdsb.2018183(2018)[2] Andronov, A.A.: Application of Poincar´e theorem on bifurcation points and change in stabilityto simple auto-oscillatory systems. C. R. Acad. Sci. Paris 189, 559-561 (1929)
3] Bajaj, A.K., Sethna, P.R.: Bifurcations in three-dimensional motions of articulated tubes. I -Linear systems and symmetry. II - Nonlinear analysis. J. Appl. Mech. 49, 606-618 (1982)[4] Battelino, P.M., Grebogi, C., Ott, E., Yorke, J.A.: Chaotic attractors on a 3-torus, and torusbreak-up. Physica D 39, 299-314 (1989)[5] Baurmann, M., Gross, T., Feudel, U.: Instabilities in spatially extended predator-prey sys-tems: spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations. J. Theor.Bio. 245, 220-229 (2007)[6] Belair, J., Campbell, S.A., Driessche, P.V.D.: Frustration, stability, and delay-induced oscil-lations in a neural network model. SIAM. J. Appl. Math. 56, 245-255 (1996)[7] Bi, P., Ruan, S.: Bifurcations in delay differential equations and applications to tumor andimmune system interaction models. SIAM J. Appl. Dyn. Syst. 12, 1847-1888 (2013)[8] Buono, P.L., B´elair, J.: Restrictions and unfolding of double Hopf bifurcation in functionaldifferential equations. J. Differ. Equ. 189, 234-266 (2003)[9] Campell, S.A., B´elair, J.: Analytical and symbolically-assisted investigation of Hopf bifurca-tions in delay-differential equations. Can. Appl. Math. Q. 3, 137-154 (1995)[10] Campell, S.A., B´elair, J., Ohira, T., Milton, J.: Limit cycles, tori, and complex dynamics ina second-order differential equations with delayed negative feedback. J. Dyn. Differ. Equ. 7,213-236 (1995)[11] Campell, S.A., LeBlanc, V.G.: Resonant Hopf-Hopf interaction in delay differential equations.J. Dyn. Differ. Equ. 10, 327-346 (1998)[12] Chen, S., Shi, J., Wei, J.: Global stability and Hopf bifurcation in a delayed diffusive Leslie-Gower predator-prey system, Int. J. Bifurcat. Chaos 22, 331-517 (2012)[13] Chen, S., Yu, J.: Stability and bifurcations in a nonlocal delayed reaction-diffusion populationmodel. J. Differ. Equations 260, 218-240 (2016)[14] De Wit, A., Dewel, G., Borckmans, P.: Chaotic Turing-Hopf mixed mode. Phys. Rev. E 48,R4191-R4194 (1993)[15] Du, Y., Guo, Y., Xiao, P.: Freely-moving delay induces periodic oscillations in a structuredSEIR model. Int. J. Bifurcat. Chaos 27, 1750122 (2017)[16] Eckmann, J.P.: Roads to turbulence in dissipative dynamical systems. Rev. Modern Phys. 53,643-654 (1981)[17] Elphick, C., Tiraopegui, E., Brachet, M.E., Coullet, P., Iooss, G.: A simple global character- zation for normal forms of singular vector fields. Physica D 29, 95-127 (1987)[18] Faria, T.: Normal forms and Hopf bifurcation for partial differential equations with delays.Trans. Amer. Math. Soc. 352, 2217-2238 (2000)[19] Faria, T., Stability and bifurcation for a delayed predator-prey model and the effect of diffu-sion. J. Math. Anal. Appl. 254, 433-463 (2001)[20] Faria, T., Huang, W.: Stability of periodic solutions arising from Hopf bifurcation for areaction-diffusion equation with time delay. Fields Inst. Comm. 31, 125-141 (2002)[21] Faria, T., Magalh˜aes, L.T.: Normal forms for retarded functional differential equations withparameters and applications to Hopf bifurcation. J. Differ. Equations 122, 181-200 (1995)[22] Faria, T., Magalh˜aes, L.T.: Normal form for retarded functional differential equations andapplications to Bogdanov-Takens singularity. J. Differ. Equations 122, 201-224 (1995)[23] Gils, S.A.V., Krupa, M., Langford, W.F.: Hopf bifurcation with non-semisimple 1:1 resonance.Nonlinearity 3, 825-850 (1990)[24] Govaerts, W., Guckenheimer, J., Khibnik, A.: Defining functions for multiple Hopf bifurca-tions. SIAM J. Numer. Anal. 34, 1269-1288 (1997)[25] Guckenheimer, J., Holmes, P.: Nonlinear Oscillations, Dynamical Systems, and Bifurcationsof Vector Fields. Springer, New York (1983)[26] Guo, S.: Stability and bifurcation in a reaction-diffusion model with nonlocal delay effect. J.Differ. Equations 259, 1409-1448 (2015)[27] Guo, S., Ma, L.: Stability and bifurcation in a delayed reaction-diffusion equation with Dirich-let boundary condition. J. Nonl. Sci. 26, 545-580 (2016)[28] Hale, J.K., Kocak, H.: Dynamics and Bifurcations. Springer, New York (1991)[29] Hale, J.K., Lunel, S.M.V.: Introduction to Functional Differential Equations. Springer, NewYork (1993)[30] Hassard, B.D., Kazarinoff, N.D., Wan, Y.H.: Theory and Applications of Hopf Bifurcation.Cambridge Univ. Press, New York (1981)[31] Hethcote, H.W., Lewis, M.A., Driessche, P.V.D.: An epidemiological model with a delay anda nonlinear incidence rate. J. Math. Biol. 27, 49-64 (1989)[32] Hopf, E.: Abzweigung einer periodischen l¨osung eines Differential Systems. Berichen Math.Phys. Kl. S¨ach. Akad. Wiss. Leipzig 94, 1-22 (1942)[33] Hsu, S.B., Huang, T.W.: Global stability for a class of predator-prey systems. SIAM J. Appl. ath. 55, 763-783 (1995)[34] Ji, J., Li, X., Luo, Z.: Two-to-one resonant Hopf bifurcations in a quadratically nonlinearoscillator involving time delay. Int. J. Bifurcat. Chaos 22, 1250060 (2012)[35] Kielh¨ofer, H.: Bifurcation Theory: An Introduction with Applications to Partial DifferentialEquations. Springer, New York (2011)[36] Kuznetsov, Y.A.: Elements of Applied Bifurcation Theory. Springer, New York (2011)[37] Lewis, G.M., Nagata, W.: Double Hopf bifurcations in the differentially heated rotating an-nulus. SIAM J. Appl. Math. 63, 1029-1055 (2003)[38] Lin, X., So, J.W.H., Wu, J.: Centre manifolds for partial differential equations with delays.P. Roy. Soc. Edinb. A 122, 237-254 (1992)[39] Luongo, A., Paolone, A.: Perturbation methods for bifurcation analysis from multiple non-resonant complex eigenvalues. Nonlinear Dynam. 14, 193-210 (1997)[40] Ma, S., Lu, Q., Feng, Z.: Double Hopf bifurcation for van der Pol-Duffing oscillator withparametric delay feedback control. J. Math. Anal. Appl. 338, 993-1007 (2008)[41] Meixner, M., De Wit, A., Bose, S., Sch¨oll, E.: Generic spatiotemporal dynamics nearcodimension-two Turing-Hopf bifurcations. Phys. Rev. E 55, 6690-6697 (1997)[42] Poincar´e, H.: Les M´ethodes Nouvelles de la M´ecanique C´eleste. Cauthier-Villars, Paris (1892)[43] Reddy, D.V.R., Sen, A., Johnston, G.L.: Time delay effects on coupled limit cycle oscillatorsat Hopf bifurcation. Physica D 129, 15-34 (1999)[44] Revel, G., Alonso, D.M., Moiola, J.L.: Interactions between oscillatory modes near a 2:3resonant Hopf-Hopf bifurcation. Chaos 20, 113-129 (2010)[45] Revel, G., Alonso, D.M., Moiola, J.L.: Numerical semi-global analysis of a 1:2 resonant Hopf-Hopf bifurcation. Physica D 247, 40-53 (2013)[46] Ruan, S., Xiao, D.: Global analysis in a predator-prey system with nonmonotonic functionalresponse. SIAM J. Appl. Math. 61, 1445-1472 (2000)[47] Ruelle, D., Takens, F.: On the nature of turbulence. Comm. Math. Phys. 20, 167-192 (1971)[48] Song, Y., Wei, J.: Local Hopf bifurcation and global periodic solutions in a delayed predator-prey system. J. Math. Anal. Appl. 301, 1-21 (2005)[49] Song, Y., Zhang, T., Peng, Y.: Turing-Hopf bifurcation in the reaction-diffusion equationsand its applications. Commun. Nonlinear Sci. Numer. Simulat. 33, 229-258 (2016)[50] Steen, P.H., Davis, S.H.: Quasiperiodic bifurcation in nonlinearly-coupled oscillators near a oint of strong resonance. SIAM J. Appl. Math. 42, 1345-1368 (1982)[51] Su, Y., Wei, J., Shi, J.: Hopf bifurcations in a reaction-diffusion population model with delayeffect. J. Differ. Equations 247, 1156-1184 (2009)[52] Wiggins, S.: Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, NewYork (2003).[53] Wu, J.: Theory and Applications of Partial Functional-Differential Equations. Springer, NewYork (1996)[54] Xiao, D.: Bifurcations of a ratio-dependent predator-prey system with constant rate harvest-ing. SIAM J. Appl. Math. 65, 737-753 (2005)[55] Xiao, Y., Chen, L.: An SIS epidemic model with stage structure and a delay. Acta Math.Appl. Sin. E. 18, 607-618 (2002)[56] Xu, X., Wei, J.: Turing-Hopf bifurcation of a class of modified Leslie-Gower model withdiffusion. Disc. Continu. Dyn. Sys. B 23, 765-783 (2018)[57] Yan, X., Li, W.: Stability of bifurcating periodic solutions in a delayed reaction-diffusionpopulation model. Nonlinearity 23, 1413-1431 (2010)[58] Yi, F., Wei, J., Shi, J.: Bifurcation and spatiotemporal patterns in a homogenous diffusivepredator-prey system. J. Differ. Equations 246, 1944-1977 (2009)[59] Yu, P.: Analysis on double Hopf bifurcation using computer algebra with the aid of multiplescales. Nonlinear Dynamics 27, 19-53 (2002)[60] Yu, P., Bi, Q.: Analysis of non-linear dynamics and bifurcations of a double pendulum. J.Sound Vib. 217, 691-736 (1998)[61] Yu, P., Yuan, Y., Xu, J.: Study of double Hopf bifurcation and chaos for oscillator with timedelay feedback. Commun. Nonlinear Sci. Numer. Simul. 7, 69-91 (2002)[62] Zhang, Y., Xu, J.: Classification and computation of non-resonant double Hopf bifurcationsand solutions in delayed van der Pol-Duffing system. Int. J. Nonlinear Sci. Numer. Simul. 6,67-74 (2005) IG. 4. When d = 5 .
23, and ω = 0 .
53 in D
4, two stable nonhomogeneous periodic solutionscoexist when we choose two different initial conditions. For (a-c) initial conditions are S ( x, t ) =1 . .
01 cos x , I ( x, t ) = 5 . − .
06 cos x , y ( x, t ) = 4 . − .
05 cos x , t ∈ [ − τ, S ( x, t ) = 1 . .
01 cos 2 x , I ( x, t ) = 5 . − .
06 cos 2 x , y ( x, t ) = 4 . − .
05 cos 2 x , t ∈ [ − τ, ) r τ HH τ τ τ τ τ τ τ τ b) r τ FIG. 5. When r = 1 , a = 1 , a = 1 . , a = 2 . , a = 1 , d = 0 . , d = 0 . , l = 3, (a) thepartial bifurcation set on the r − τ plane and (b) the complete bifurcation set around HH areshown. X(0,t) Y ( ,t ) b)X(0,t) Y ( ,t ) a) X(0,t) Y ( ,t ) c) FIG. 6. (a) When τ = 10 . r = 0 .
69, there exists a stable quasi-periodic solution on a 2-torusof system (73); (b) when τ = 10 . r = 0 .
71, there exists a quasi-periodic solution on a 3-torus ofsystem (73); (c) τ = 10 . r = 0 .726, there is a strange attractor of system (73), and the systemexhibits chaotic behavior.