Periodic solutions in an epidemic model with diffusion and delay

Applied Mathematics and Computation 265 (2015) 275–291

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Periodic solutions in an epidemic model with diffusion and delay Pan-Ping Liu∗ Department of Mathematics, North University of China, Taiyuan, Shanxi 030051, People’s Republic of China

a r t i c l e

i n f o

PACS: 87.23.Cc 89.75.Kd 89.75.Fb 47.54.-r Keywords: Hopf bifurcation Epidemic models Time delay Spatial diffusion

a b s t r a c t A spatial diffusion SI model with delay and Neumann boundary conditions are investigated. We derive the conditions of the existence of Hopf bifurcation in one dimension space. Moreover, we analyze the properties of bifurcating period solutions by using the normal form theory and the center manifold theorem of partial functional differential (PFDs) equations. By numerical simulations, we found that spatiotemporal periodic solutions can occur in the epidemic model with spatial diffusion, which verifies our theoretical results. The obtained results show that interaction of delay and diffusion may induce outbreak of infectious diseases. © 2015 Elsevier Inc. All rights reserved.

1. Introduction From the 2nd century A.D., the Roman Empire occurred Antonine plague [1], then the last few years there were outbreaks of SARS, H1N1 avian influenza, H7N9 avian influenza in the global world. These new emerging infectious diseases continuously increase every year, and the traditional infectious diseases recurrently outbreak [2–4]. Based on pathology of the disease, the route of transmission, and characteristics of the outbreak, we can establish the appropriate mathematical models and in turn, some reasonable suggestions can be provided for the prevention and effective control of infectious diseases. With the study on mechanism of infectious diseases, many researchers found that the corresponding symptoms did not appear immediately after an individual was infected, and the individual took a period of time (namely exposed period) to show symptoms [5–8], such as rabies, AIDS, cholera and so on. In this case, time delay is needed to describe such phenomenon, which reflects that the changes of t moment not only depend on the state at t moment, but also are influenced by some factors before t moment [9–11]. Since an individual does not always stay in one place, conversely, it diffuses around. Therefore, the original reaction model becomes reaction–diffusion model. Moreover, the reaction–diffusion epidemic model with delay is able to better reflect the reality of disease [12–19]. For this kind of model, Thieme and Zhao analyzed a reaction–diffusion system with delay which described the spatial spread of a series bacterium and virus epidemic, then they obtained the traveling wave solutions and the existence of asymptotic propagation speed [20–22]. A spatial diffusion SEI epidemic model with delay was investigated by Kim and Lin, they obtained the sufficient conditions of the local and global asymptotical stability by using upper and lower solutions and its associated monotone iterations [23]. Li et al. studied pattern formation on an epidemic model with time delay, and concluded that the delay and diffusion can induce Turing pattern [10]. However, the studies of bifurcation behaviors on epidemic models with time delay and diffusion are still limited. Hence, in this paper, we studied the existence and the properties of Hopf bifurcation. ∗

Tel.: +8613513656782. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.amc.2015.05.028 0096-3003/© 2015 Elsevier Inc. All rights reserved.

276

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

The structure of this paper is as follows. In Section 2, we derive the characteristic equation of spatial diffusion SI model with delay, then the existence of the Hopf bifurcation is analyzed. In Section 3, by utilizing the normal form theory and the center manifold theorem, we get some properties of Hopf bifurcation. In Section 4, the results of the numerical simulation indicate that time delay can induce periodic outbreak of disease in spatial diffusion SI model with delay. Finally, we give some conclusions and discussions. 2. Mathematical modeling and the existence of Hopf bifurcation 2.1. Model formulation A simple SI epidemic system produced surprising dynamics [24]. However, the authors did not consider the SI system includes both spatial factor and latent case. Thus, we will study the Hopf bifurcation on the basis of this system. The simple SI system is





N dS SI = rN 1 − − β − (μ + m )S, dt K N dI SI = β − (μ + d )I. dt N

(1)

The total population (N) is divided into susceptible groups (S) and infectious groups (I), where r is the intrinsic growth rate, K is the carrying capacity of logistic equation, β denotes the contact transmission rate, μ is the natural mortality, d denotes the disease induced morality, m is the per-capital emigration rate of uninfectives. The basic reproduction number is

R0 =

β , μ+d

and the basic demographic reproductive number is

Rd =

r

μ+m

.

By re-scaling the above system, we have an epidemic system with spatial effect [25]:

∂S SI = vRd (S + I )(1 − (S + I )) − R0 − vS + d1 ∇ 2 S, ∂t S+I ∂I SI = R0 − I + d2 ∇ 2 I, ∂t S+I

(2)

+m where v = μ μ+d is the ratio of the average life-span of susceptibles to that of infections, d1 , d2 are diffusion coefficients. In this

paper, we consider one-dimension space domain, then the Laplacian operator ∇ 2 = ∂∂x2 . For most diseases, they have latent periods [26–28]. For example, Rabies may take several days or months to reach the infectious stage [29–32]. For such diseases, we need to introduce the time delay into the infected population. Moreover, we assume that the population can-not pass the boundary of the domain, and the outside population cannot enter this domain. As a result, we have the following system with Neumann boundary conditions: 2

⎧ SI (t − τ ) ∂S ⎪ 2 ⎪ ⎪ ∂ t = vRd [S + I (t − τ )][1 − (S + I (t − τ ))] − R0 S + I (t − τ ) − vS + d1 ∇ S, ⎨  ∂I SI (t − τ ) ∂S  − I (t − τ ) + d2 ∇ 2 I, t ≥ 0, x ∈ (0, l π ), = 0, = R0 ⎪ ⎪ S + I (t − τ ) ∂ x x=0,lπ ⎪ ⎩ ∂t S(x, t ) = φ1 (x, t ) ≥ 0, I (x, t ) = φ2 (x, t ) ≥ 0, (x, t ) ∈ [0, l π ] × [−τ , 0].

 ∂I  = 0, t ≥ 0, ∂ x x=0,lπ

(3)

Assuming φ = (φ1 , φ2 )T ∈ ℘ = C ([−τ , 0], X ), τ > 0 and X is defined as



X=

(S, I ) : S, I ∈ W T

2,2

  ∂S  ∂I  (0, l π );  = =0 ∂ x x=0,lπ ∂ x x=0,lπ

with the inner product  ·, ·. 2.2. Existence of Hopf bifurcation At first, in the absence of diffusion and delay, system (3) is corresponding to the following system:

∂S SI = vRd (S + I )(1 − (S + I )) − R0 − vS, ∂t S+I ∂I SI = R0 − I. ∂t S+I

(4)

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

According to [24], for Rd > 1 and R0 < 1, the system (4) has a disease-free equilibrium E0 = (S0 , I0 ) = (1 − this system also has an endemic equilibrium E ∗ = (S∗ , I∗ ), which is satisfied

(A1 ) Rd >

R0 + v − 1 , vR0

277 1 Rd

, 0 ). Moreover,

R0 > 1,

then

S∗ =

vR0 Rd − R0 + 1 − v , I∗ = (R0 − 1 )S∗ . vR20 Rd

Let s = S − S∗ , h = I − I∗ , then system (3) can be translated into the following system:

∂s (s + S∗ )(h(t − τ ) + I∗ ) = vRd (s + S∗ + h(t − τ ) + I∗ )[1 − (s + S∗ + h(t − τ ) + I∗ )] − R0 − v(s + S∗ ) + d1 ∇ 2 s, ∂t s + S∗ + h(t − τ ) + I∗ ∂h (s + S∗ )(h(t − τ ) + I∗ ) = R0 − (h(t − τ ) + I∗ ) + d2 ∇ 2 h. ∂t s + S∗ + h(t − τ ) + I∗

(5)

Defining

f (1) (s, h ) = vRd (s + S∗ + h(t − τ ) + I∗ )[1 − (s + S∗ + h(t − τ ) + I∗ )] − R0 f (2) (s, h ) = R0

(s + S∗ )(h(t − τ ) + I∗ ) − (h(t − τ ) + I∗ ), s + S∗ + h(t − τ ) + I∗

(s + S∗ )(h(t − τ ) + I∗ ) − v ( s + S ∗ ), s + S∗ + h(t − τ ) + I∗

and for i, j = 0, 1, 2, . . . , let

fi(j1) =

∂ i+ j f (1) (0, 0 ), ∂ si ∂ h j

fi(j2) =

∂ i+ j f (2) (0, 0 ), i + j ≥ 1. ∂ si ∂ h j

Further we get

dU (t ) = DU (t ) + L(Ut ) + F (Ut ), dt

(6)

where U = (u1 , u2 )T , s(t, x ) = u1 (t, x ), h(t, x ) = u2 (t, x ), D = (d01 d0 ). We set Ut (θ ) = U (t + θ ), φ = (φ1 , φ2 )T ∈ ℘, φ (θ ) = Ut (θ ), and θ ∈ [−τ , 0]. Let L : ℘ → X and F : ℘ → X are given by



L (φ ) =

2



a11 φ1 (0 ) + a12 φ2 (−τ ) , a21 φ1 (0 ) + a22 φ2 (−τ )

(1 ) (1 ) (2 ) (2 ) where a11 = ∂ ∂f s |(0,0 ) , a12 = ∂ ∂f h |(0,0 ) , a21 = ∂ ∂f s |(0,0 ) , a22 = ∂ ∂f h |(0,0 ) , and

⎛ 1 ⎞ fi(j1) (0, 0 )φ1i (0 )φ2j (−τ ) ⎜i+ j≥2 i! j! ⎟ ⎟. F (φ ) = ⎜  1 ⎝ ⎠ fi(j2) (0, 0 )φ1i (0 )φ2j (−τ ) i+ j≥2

i! j!

We obtain the linearing part of (6)

dU (t ) = DU (t ) + L(Ut ), dt

(7)

then we set U (t ) = yeλt , and y = (y1 , y2 )T , the characteristic equation is

λy − D

∂ 2y − L(eλ. y ) = 0, ∂ x2

(8)

where y ∈ dom( ∂∂x2 ) and y = 0, dom( ∂∂x2 ) ∈ X. 2

2

2 2 According to the properties of the Laplacian operator, it is clear that ∂∂x2 on X have eigenvalues − nl 2 with the corresponding

eigenfunctions βn1 = (γ0n ), βn2 = (γ0 ), where γn = cos nl x, n ∈ N0 = {0, 1, 2 . . .}. Clearly, {βn1 , βn2 }∞ construct a basis of the phase n=0 n space X, therefore any element y in X can be expanded as Fourier series in the following form:

y=

∞  n=0



YnT

 1  y, βn βn1  . , and Yn =  βn2 y, βn2

(9)

278

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

Besides, by some simple computations, we get

  1  1 βn βn T T L φ = L (φ ) , n ∈ N0 . βn2 βn2

(10)

Based on the above (9) and (10), (8) becomes ∞ 

YnT

n=0

   1 βn n2 a11 a12 e−λτ = 0, λI2 + D 2 − a21 a22 e−λτ l βn2

(11)

we deduce the eigenpolynomial associated with λ :

 n  n 2 4 λ2 − −(d1 + d2 ) + a11 + a22 e−λτ λ + d1 d2 l l n 2 − (d1 a22 e−λτ + d2 a11 ) + (a11 a22 − a12 a21 )e−λτ = 0,

(12)

l

where

R20 + vR0 Rd + vR0 − 4R0 − 2v + 3 , R0 2R0 − vR0 Rd + 2v − 3 = , R0

a11 = − a12

(R0 − 1 )2

a21 =

R0

> 0,

1 − 1 < 0, R0

a22 =

for R0 > 1.

Moreover, a11 + a22 < 0, and a11 a22 − a12 a21 > 0. When n = 0, the eigenpolynomial is

λ2 − (a11 + a22 e−λτ )λ + (a11 a22 − a12 a21 )e−λτ = 0.

(13)

Let λ = iw(w > 0 ), then

(iw )2 − (a11 + a22 e−iwτ )iw + (a11 a22 − a12 a21 )e−iwτ = 0, − w2 − [a11 + a22 (coswτ − isinwτ )]iw + (a11 a22 − a12 a21 )(coswτ − isinwτ ) = 0. Separating the real and imaginary parts, one can have



−w2 = a22 wsinwτ − (a11 a22 − a12 a21 )coswτ ,

(14)

−a11 w = a22 wcoswτ + (a11 a22 − a12 a21 )sinwτ . Further, by squaring and adding the two parts of Eq. (14), then we have





w4 + a211 − a222 w2 − (a11 a22 − a12 a21 )2 = 0. So we have



2

w =



− a211 − a222 +



a211 − a222

2

(15)

+ 4(a11 a22 − a12 a21 )2

2

.

Obviously, w2 > 0 holds, so Eq. (13) has a pair of simple imaginary roots which is ±iw0 . From (14), we can obtain

cos(w0 τ ) = sin(w0 τ ) =

−a12 a21 w20

a222 w20

+ (a11 a22 − a12 a21 )2

= P ( w0 ),

−a11 w0 (a11 a22 − a12 a21 ) − a22 w30 a222 w20 + (a11 a22 − a12 a21 )2

= Q ( w0 ),

furthermore, some simple derivations show that

τ0j =

⎧ 1 ⎪ ⎨ arccos(P (w0 ) + 2 jπ ), ⎪ ⎩

w0 1 (2π − arccos(P (w0 )) + 2 jπ ), w0

when Q (w0 ) ≥ 0, when Q (w0 ) ≤ 0,

j = 0, 1, 2, . . ..

Let λ(τ ) = α (τ ) + iw(τ ) be the root of (13) near τ0 which satisfies α (τ0 ) = 0 and w(τ0 ) = w0 , where j = 0, 1, 2, . . .. j

j

j

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

279

Lemma 2.1. For n = 0, assuming

(A2 ) a11 − a22 < 0, then the transversality condition satisfies

 Re

dλ dτ



τ =τ0j > 0.

Proof. We take the derivation with respect to τ on Eq. (13), then we derive

dλ {2λ − (a11 + a22 e−λτ ) + a22 λτ e−λτ − (a11 a22 − a12 a21 )τ e−λτ } + [a22 λ − (a11 a22 − a12 a21 )]λe−λτ = 0, dτ  dλ −1 ⇒ Re dτ τ =τ0j

−(a11 a22 − a12 a21 )τ e−λτ 2λ − (a11 + a22 e−λτ ) + a22 λτ e−λτ = Re + [−a22 λ + (a11 a22 − a12 a21 )]λe−λτ [−a22 λ + (a11 a22 − a12 a21 )]λe−λτ τ =τ0j

−(a11 a22 − a12 a21 )τ 2λeλτ − (a11 eλτ + a22 ) + a22 λτ = Re + [−a22 λ + (a11 a22 − a12 a21 )]λ [−a22 λ + (a11 a22 − a12 a21 )]λ τ =τ0j

2λeλτ − (a11 eλτ + a22 ) τ  − = Re  λ τ =τ j − a22 λ + (a11 a22 − a12 a21 ) λ 0 

2 2 2 2a22 w20 − a11 (a11 a22 − a12 a21 ) w0 sinw0 τ − (a11 a22 − 2a12 a21 )w0 cosw0 τ − a22 w0  2 2  2 +  2 2  2 = a22 w0 + (a11 a22 − a12 a21 )2 w0 a22 w0 + (a11 a22 − a12 a21 )2 w0 τ =τ0j      a211 − a222 (a11 a22 − a12 a21 )2 2a222 w40 + 2(a11 a22 − a12 a21 )2 + a222 a211 − a222 w20 = +  2 2 2 2 . a22 w0 + (a11 a22 − a12 a21 )2 a222 w20 + (a11 a22 − a12 a21 )2 From the assumption (A2), and combining a11 + a22 < 0, we can get a211 − a222 > 0. Consequently,

 Re

dλ dτ



−1



τ =τ > 0 ⇒ j 0

Re

dλ dτ



τ =τ0j > 0.



Theorem 2.1. If (A1), (A2) are all established , in the absence of diffusion, system (3) undergoes a spatially homogeneous Hopf bifurj cation at equilibrium E ∗ = (S∗ , I∗ ) when τ = τ0 , and period solution will emerge. Lemma 2.3. Suppose that

(S1 ) R0 + v − 2 ≥ 0, and (S2) if there exists a certain n0 ∈ N = {1, 2, . . .}, which makes

d1 d2

 n 4 0

l

− (d2 a11 − d1 a22 )

 n 2 0

l

− (a11 a22 − a12 a21 ) < 0.

(16)

Then equation (12) has a pair purely imaginary roots ±iwn0 , and

√

w n0 =

2 2



−Bn0 +



B2n0 − 4Cn0 ,

(17)

this moment, τ is satisfied τn0 . j

Proof. Supposing n = n0 ∈ N, let λ = iw(w > 0 ) be a root of equation (12), and setting k = by using the same way as before:

n , then getting the following equation l

(iw )2 − [a11 + a22 e−iwτ − (d1 + d2 )k2 ]iw + d1 d2 k4 − (d1 a22 e−iwτ + d2 a11 )k2 + (a11 a22 − a12 a21 )e−iwτ = 0.

(18)

Separating the real and imaginary parts, one can have



−w2 + d1 d2 k4 − d2 a11 k2 = a22 wsinwτ + d1 a22 k2 coswτ − (a11 a22 − a12 a21 )coswτ ,

(d1 + d2 )k2 w − a11 w = a22 wcoswτ − d1 a22 k2 sinwτ + (a11 a22 − a12 a21 )sinwτ .

(19)

280

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

By squaring and adding the two parts of Eq. (19), then we obtain

w4 +







d12 + d22 k4 − 2d1 a11 k2 + a211 − a222



w2 + [d1 d2 k4 − d2 a11 k2 ]2 − [d1 a22 k2 − (a11 a22 − a12 a21 )]2 = 0,

(20)

namely w4 + Bw2 + C = 0, ∀n ∈ N. Further, we can deduce 2

w =

−B ±

√ B2 − 4C . 2

(21)

If w > 0 exists, then w2 needs to be satisfied w2 > 0.

C = C1 × C2 = [d1 d2 k4 − d2 a11 k2 − d1 a22 k2 + (a11 a22 − a12 a21 )] × [d1 d2 k4 − d2 a11 k2 + d1 a22 k2 − (a11 a22 − a12 a21 )]. For assumptions (A1) and (S1), we have a11 < 0, then −d2 a11 k2 > 0, so C1 > 0, ∀n ∈ N. n For assumption (S2), we obtain C2 (k0 ) < 0 for k0 = l0 . Consequently, w2n0 > 0, we deduce

w n0

√   2 = −Bn0 + B2n0 − 4Cn0 . 2

On the basis of two equations of (19),

cos(wn0 τ ) =

−d12 d2 a22 ( nl0 )6 + d1 d2 (2a11 a22 − a12 a21 )( nl0 )4 a222 w2n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2

+ sin(wn0 τ ) =

−d2 [a11 (a11 a22 − a12 a21 ) + a22 w2n0 ]( nl0 )2 + a12 a21 w2n0 a222 w2n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2

= X ( w n0 ) ,

−d12 a22 wn0 ( nl0 )4 + d1 (2a11 a22 − a12 a21 )wn0 ( nl0 )2 a222 w2n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2

+

−d2 a12 a21 wn0 ( nl0 )2 − (a11 a22 − a12 a21 )a11 wn0 − a22 w3n0 a222 w2n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2

= Y ( w n0 ) ,

then we get

τnj0 =

1 wn0

arccos(X (wn0 ) + 2 jπ ),

1 wn0

(2π − arccos(X (wn0 )) + 2 jπ ), when Y (wn0 ) ≤ 0,

when Y (wn0 ) ≥ 0, j = 0, 1, 2, . . ..

Let λ(τ ) = α (τ ) + iw(τ ) be the root of (12) near τn0 which satisfies α (τn0 ) = 0 and w(τn0 ) = wn0 , where j = 0, 1, 2, . . .. j

j

j

Lemma 2.4. For a certain n = n0 ∈ N, the transversality condition

 Re

dλ dτ



τ =τnj0 > 0

holds. Proof. Taking the derivation on Eq. (12) with respect to τ , then we get





2λ − a11 + a22 e

−λτ

 n 2  dλ



 n 2

dλ − (a11 a22 − a12 a21 ) τ e−λτ dτ

+ a22 λ + d1 a22 l   n 2 + a22 λ + d1 a22 − (a11 a22 − a12 a21 ) λe−λτ = 0, l 

−2λeλτ + a11 eλτ + a22 − (d1 + d2 )( nl )2 eλτ dλ −1 τ ⇒ Re = Re − dτ τ =τnj0 λ τ =τ j [a22 λ + d1 a22 ( nl )2 − (a11 a22 − a12 a21 )]λ n0

n 2 − ( d + d )( ) (coswτ + isinwτ ) a −2iw(coswτ + isinwτ ) + a11 (coswτ + isinwτ ) 22 1 2 l = Re + iw[a22 iw + d1 a22 ( nl )2 − (a11 a22 − a12 a21 )] iw[a22 iw + d1 a22 ( nl )2 − (a11 a22 − a12 a21 )] τ =τ j − ( d 1 + d2 )



l

dτ

n0

= =

−2w2n0 [−w2n0 + d1 d2 ( nl0 )4 − d2 a11 ( nl0 )2 ]

a222 w4n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2 w2n0

(

)( )

a222 w4n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2 w2n0

) − a22 w2n0 ( ) − (a11 a22 − a12 a21 )]2 w2n0

d1 + d2 nl0 2 w2n0 [−a11 a222 w4n0 + [d1 a22 nl0 2

+ (d1 + d2 )(

+

n0 2 ] l

−a11 w2n0 [−a11 + (d1 + d2 )( nl0 )2 ]



P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

= = = =

281

[2w2n0 + (d12 + d22 )( nl0 )4 − 2d1 a11 ( nl0 )2 + (a211 − a222 )]w2n0 a222 w4n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2 w2n0 [2w2n0 + Bn0 ]w2n0

a222 w4n0 + [d1 a22 ( nl0 )2 − (a11 a22 − a12 a21 )]2 w2n0 w2n0 [−Bn0 + a222 w4n0

+



( ) − (a11 a22 − a12 a21 )]2 w2n0  B2n0 − 4Cn0

w2n0 a222 w4n0 +

B2n0 − 4Cn0 + Bn0 ]

[d1 a22 nl0 2

( ) − (a11 a22 − a12 a21 )]2 w2n0

[d1 a22 nl0 2

> 0.

Lemma 2.4 is proved.  Theorem 2.2. If the conditions (A1), (S1) and (S2) are satisfied, in the presence of space, system (3) experiences a spatially inhomogej neous Hopf bifurcation at E ∗ = (S∗ , I∗ ) when τ = τn0 , and period solution will appear. 3. Direction and stability of Hopf bifurcation From the above part, we obtain conditions of the existence of Hopf bifurcation, next we will investigate properties of these bifurcating periodic solutions from the positive constant steady state E∗ (I∗ , H∗ ) of system (3) which includes the direction, stability and period. Here we will use the normal form theory and the center manifold theorem of partial functional differential j equations (PFDEs) [33,34]. For the sake of simplicity, we set τc = τn (n = 0, n0 ; j = 0, 1, 2, . . . ) and k = nl . Let s¯(t, x ) = s(τ t, x ), h¯ (t, x ) = h(τ t, x ), and drop the bar which do not affect the results, then system (3) is transformed into

 1 (1 ) ∂s ∂ 2 s(t ) = τ d1 f (0, 0 )si (x, t )h j (x, t − 1 ), + τ [a11 s + a12 h(t − 1 )] + τ ∂t i! j! i j ∂ x2 i+ j≥2  1 (2 ) ∂h ∂ 2 h(t ) = τ d2 f (0, 0 )si (x, t )h j (x, t − 1 ). + τ [a21 s + a22 h(t − 1 )] + τ ∂t i! j! i j ∂ x2 i+ j≥2

(22)

We make a little disturbance on τ , namely τ = τc + a, a ∈ R. Setting U = (u1 , u2 )T , s(t, x ) = u1 (t, x ), h(t, x ) = u2 (t, x ), then (22) can be rewritten in the space ℘ = C ([−1, 0], X ),

dU (t ) = τc DU (t ) + L(τc )Ut + F (Ut , a ). dt

(23)

We set Ut (θ ) = U (t + θ ), φ = (φ1 , φ2 )T ∈ ℘, and φ (θ ) = Ut (θ ), and θ ∈ [−1, 0]. L(b)(· ) : R × ℘ −→ X (b is τ c or a) and F : ℘× R −→ X are respectively defined as





a φ (0 ) + a12 φ2 (−1 ) L(b)(φ ) = b 11 1 , a21 φ1 (0 ) + a22 φ2 (−1 ) and

F (φ , a ) = aDφ (0 ) + L(a )(φ ) + f (φ , a ), where

⎛! f ( φ , a ) = ( τc + a ) ⎝ !

i+ j≥2 i+ j≥2

1 f (1 ) i! j! i j

(0, 0 )φ1i (0 )φ2j (−1 )

1 f (2 ) i! j! i j

(0, 0 )φ1i (0 )φ2j (−1 )

⎞ ⎠.

Now consider the linearing part of (23),

U˙ (t ) = τc DU (t ) + L(τc )Ut .

(24)

Based on the discussion of Section 2, it is obvious that the equilibrium of Eq. (22) is the origin (0,0), and then we derive that the characteristic equation of Eq. (24) has a pair of purely imaginary eigenvalues 0 = {iwn τc , −iwn τc }, (n = 0, n0 ). Considering the ordinary functional differential equation:

X˙ (t ) = −τc Dk2 X (t ) + L(τc )(Xt ).

(25)

By the Riesz representation theorem, for φ ∈ C ([−1, 0], X ), there exists a 2 × 2 matrix function η (θ , τc )(−1 ≤ θ ≤ 0 ), then we have:

−τc Dk2 φ (0 ) + L(τc )(φ ) =

"

0 −1

d[η (θ , τc )]φ (θ ).

(26)

282

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

This moment,

⎧  −d1 k2 + a11 0 ⎪ ⎪ , θ = 0, τc ⎪ ⎪ −d2 k2 a21 ⎪ ⎨ θ ∈ (−1, 0 ), 0, η ( θ , τc ) = ⎪  ⎪ ⎪ 0 −a12 ⎪ ⎪ , θ = −1. ⎩ τc

(27)

−a22

0

For φ ∈ C ([−1, 0], X ), defining semigroup induced by the solution of the linear equation (25), and the infinitesimal generator A(τ c ) of the semigroup is

⎧ d φ (θ ) ⎪ , θ ∈ [−1, 0 ), ⎨ dθ A ( τc ) φ ( θ ) = " 0 ⎪ ⎩ dθ η (θ , τc )φ (θ ), θ = 0.

(28)

−1

For ψ ∈ C([0, 1], X), the formal adjoint operators of A(τ c ) is A∗ (τ c ) which denotes [35]

A∗ ( τc ) ψ ( s ) =

⎧ d ψ (s ) ⎪ , ⎨− ⎪ ⎩

s ∈ (0, 1],

ds

"

0 −1

(29)

ds η (s, τc )ψ (−s ),

s = 0.

Here, the bilinear pairing form associated A(τ c ) with A∗ (τ c ) is

( ψ , φ ) = ψ ( 0 )φ ( 0 ) −

"

" θ

0

ψ ( ξ − θ )d η ( θ , τc )φ ( ξ )d ξ  " 0 0 −a12 = ψ ( 0 ) φ ( 0 ) − τc ψ (ξ + 1 ) φ ( ξ )d ξ . ξ =0

−1

−a22

0

−1

(30)

By the discussion of Section 2, it is easy to know that A(τ c ) has a pair purely imaginary eigenvalues ±iwn τ c , in fact they are also eigenvalues of A∗ (τ c ). We set P and P∗ as the center subspace, namely the generalized eigenspace of A(τ c ) and A∗ (τ c ) associated with 0 , respectively. Moreover, P∗ is the adjoint space of P and dimP = dimP∗ = 2 [34]. Next we give the following results by some computations. Lemma 3.1. Let

ξ=

(iwn + d1 k2 − a11 )eiwn τc a12

,

η=

iwn + d1 k2 − a11 , a21

(31)

then a basis of P with 0 is

q1 (θ ) = eiwn τc θ (1, ξ )T ,

q2 (θ ) = q1 (θ ),

−1 ≤ θ ≤ 0,

and a basis of P∗ with 0 is

q∗1 (s ) = (1, η )e−iwn τc s ,

q∗2 (s ) = q∗1 (s ),

0 ≤ s ≤ 1.

 = (1 , 2 ) is obtained by separating the real and imaginary parts of q1 (θ ), and  is also the basis of P. Similarly, ∗ = (∗1 , ∗2 )T is also the basis of P∗ . Then we give the following calculation: q (θ ) + q2 (θ ) = 1 (θ ) = 1 2



=



Re{eiwn τc θ }



Re{ξ eiwn τc θ }



coswn τc θ

,

(d1 k2 −a11 )coswn τc (1+θ )−wn sinwn τc (1+θ ) a12

q (θ ) − q2 (θ ) = 2 (θ ) = 1 2i



=



Im{eiwn τc θ }



Im{ξ eiwn τc θ }

sinwn τc θ

(d1 k −a11 )sinwn τc (1+θ )+wn coswn τc (1+θ ) 2

a12

.

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

 ∗1 (s ) =

283



coswn τc s

(d1 k2 −a11 )coswn τc s+wn sinwn τc s

T

,

a21

  (s ) = ∗ 2



−sinwn τc s

T

−(d1 k2 −a11 )sinwn τc s+wn coswn τc s a21

.

According to the bilinear pairing form (30), we can compute:

(∗1 , 1 ) = 1 +

(d1 k2 − a11 )2 coswn τc − wn (d1 k2 − a11 )sinwn τc

1 + τc 2 +

(∗1 , 2 ) =



a12 a21



sin2wn τc 1 1+ (d1 k2 − a11 ) − sin2 wn τc 2wn τc τc





1 sin2wn τc a τc 22 −w2n 1 − 2 a12 a21 2wn τc



+ (d1 k2 − a11 )2 1 +

sin2wn τc 2wn τc



,

(d1 k2 − a11 )wn coswn τc + (d1 k2 − a11 )2 sinwn τc



a12 a21





1 sin2wn τc τc w n 1 + 2 2wn τc

+

  sin wn τc a 1 τc 22 2wn (d1 k2 − a11 ) + w2n + (d1 k2 − a11 )2 , 2 a12 a21 w n τc

+

(d1 k2 − a11 )sin2 wn τc w n τc



+



(∗2 , 1 ) =





2

wn (d1 k2 − a11 )coswn τc − (wn )2 sinwn τc a12 a21



+



1 sin2wn τc (d k2 − a11 )sin wn τc τc − 1 + wn 1 − 2 w n τc 2wn τc 2









  sin wn τc 1 a22 + τc 2wn (d1 k2 − a11 ) − w2n + (d1 k2 − a11 )2 , 2 a12 a21 w n τc (∗2 , 2 ) =

2

(wn )2 coswn τc + wn (d1 k2 − a11 )sinwn τc 1 − τc 2



a12 a21



sin2wn τc sin wn τc + (d1 k − a11 ) 1 − τc 2wn τc 1

2



2



1 sin2wn τc a22 + τc w2n 1 + 2 a12 a21 2wn τc



− (d1 k − a11 ) 2



2



sin2wn τc 1− 2wn τc



.

In the following, we can define (∗ , ) = (∗ı ,  ), (ı,  = 1, 2 ). For P∗ , we construct a new basis  , and  = (1 , 2 )T = cos nl x coskx (∗ , )−1 ∗ , this moment,  and  are satisfied ( , ) = I2 . Besides, fn = (βn1 , βn2 ), where βn1 = ( )=( ), βn2 = 0 0 0 0 ( n )=( ). For c = (c1 , c2 ) ∈ C ([−1, 0], X ), we define c · fn = c1 βn1 + c2 βn2 . cos l x coskx # lπ Define  ·, · as the complex-valued L2 inner product on X, for U = (u1 , u2 )T , V = (v1 , v2 )T ∈ X, then U, V  = l1π 0 (u1 v¯ 1 + u2 v¯ 2 )dx, furthermore, for ∀φ ∈ ℘,

φ , fn  = (φ , βn1 , φ , βn2  )T . On the basis of the decomposition of the phase space, we have ℘ = PCN℘+ Ps℘, where PCN℘ is the center subspace of linear equation (24),

PCN℘(φ ) = ( , φ , fn  ) · fn ,

φ ∈ ℘,

and Ps℘ is the complement subspace of PCN℘.

(32)

284

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

Due to the infinitesimal generator A(τ c ) is deduced by the solution of (24), so (22) can be rewritten as

U˙ t = A(τc )Ut + X0 F (Ut , a ),

where X0 (θ ) =

(33)

θ ∈ [−1, 0 ),

0,

I, θ = 0. By using the phase space decomposition ℘ = PCN℘+ Ps℘ and (32), the solution of (23) is expressed as

Ut = (x1 (t )x2 (t )) · fn + h(x1 , x2 , a ),

(34)

x1 (t ) ) = ( , Ut , fn  ), and h(x1 , x2 , a) ∈ Ps℘, h(0, 0, 0 ) = 0, Dh(0, 0, 0 ) = 0. x2 (t ) Especially, the solution of (23) on center manifold is

where (

Ut = (x1 (t )x2 (t )) · fn + h(x1 , x2 , 0 ).

(35)

Let z = x1 − ix2 , z¯ = x1 + ix2 and (0 ) = (1 (0 ), 2 (0 ))T , and let q1 = 1 + i2 , so

(x1 (t )x2 (t )) · fn = (1 , 2 )

 z+z¯ i(z−z¯ )  2

2

· fn =

1 (q1 z + q¯1 z¯ ) · fn . 2

(36)

Further (34) becomes

Ut =

1 (q1 z + q¯1 z¯ ) · fn + W (z, z¯ ), 2

(37)

z¯ i(z−z¯ ) where W (z, z¯ ) = h( z+ 2 , 2 , 0 ), and setting

W (z, z¯ ) = W20

z2 z¯2 + W11 zz¯ + W02 + . . .. 2 2

(38)

Furthermore, according to [34] z satisfies

z˙ = iwn τc z + g(z, z¯ ),

(39)

g(z, z¯ ) = (1 (0 ) − i2 (0 ))F (Ut , 0 ), fn  = (1 (0 ) − i2 (0 )) f (Ut , 0 ), fn ,

(40)

where

and setting

g(z, z¯ ) = g20

z2 z2 z¯ z¯2 + g11 zz¯ + g02 + g21 . . .. 2 2 2

(41)

From f(φ , a), (37) and (38), it is easy to compute that

 f (Ut , 0 ), fn  =

+

+

τc 4

τc 4

τc 4

$ ·

$ ·

$ ·

(1 ) (1 ) −iwn τc (1 ) 2 −2iwn τc f20 + 2 f11 ξe + f02 ξ e

%

(2 ) (2 ) −iwn τc (2 ) 2 −2iwn τc f20 + 2 f11 ξe + f02 ξ e

·

1

" lπ

lπ

0

(1 ) (1 ) ¯ iwn τc (1 ) ¯ f20 + f11 (ξ e + ξ e−iwn τc ) + 2 f02 ξξ

% ·

(2 ) (2 ) ¯ iwn τc (2 ) ¯ f20 + f11 (ξ e + ξ e−iwn τc ) + 2 f02 ξξ (1 ) (1 ) ¯ iwn τc (1 ) ¯ 2 2iwn τc f20 + 2 f11 ξe + f02 ξ e (2 ) (2 ) ¯ iwn τc (2 ) ¯ 2 2iwn τc f20 + 2 f11 ξe + f02 ξ e

%

·

(coskx )3 dx ·

1

" lπ

lπ

0

1

" lπ

lπ

0

z2 2

(coskx )3 dx · zz¯

(coskx )3 dx ·

z¯2 2

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

  (1 )

⎛



(1 ) (1 ) (1 ) f20 W11 (0 ) + 12 W20 (0 ) · coskx, coskx







⎞

(2 ) (2 ) ⎟ ⎜ + f11 W11 (−1 ) + 12 W20 (−1 ) · coskx, coskx ⎟ ⎜    ⎟ ⎜ (1 ) (1 ) (0 ) · coskx, coskx ⎟ ξ e−iwn τc W11(1) (0 ) + 12 ξ¯ eiwn τc W20 ⎜ + f11 ⎜ (1)   ⎟ (2 ) (2 ) −iwn τc ⎜+ f W11 (−1 ) + 12 ξ¯ eiwn τc W20 (−1 ) · coskx, coskx ⎟ ⎟ ⎜ 02 ξ e ⎟ ⎜ 1 (1 ) (1 ) ¯ iwn τc (1 ) (ξ e + 2ξ e−iwn τc ) + f12 (2ξ ξ¯ + ξ 2 e−2iwn τc ) ⎟ ⎜ + 8 f30 + f21 ⎟ ⎜   (1 ) 2 ¯ ⎟ ⎜ + f03 ξ ξ · coskx, coskx ⎟ z2 z¯ ⎜ ⎟· + τc · ⎜ ⎟ 2 + ..., ⎜  (1)   (2 ) (1 ) 1 ⎟ ⎜ f W ( 0 ) + W ( 0 ) · coskx, coskx 2 20 20 ⎟ ⎜  (11   ⎟ ⎜ (2 ) 2) (2 ) 1 + f11 W11 (−1 ) + 2 W20 (−1 ) · coskx, coskx ⎟ ⎜ ⎟ ⎜    ⎜ + f (2) ξ e−iwn τc W (1) (0 ) + 1 ξ¯ eiwn τc W (1) (0 ) · coskx, coskx ⎟ 2 11 11 20 ⎜  ⎟ ⎟ ⎜ (2)  −iwn τc (2) (2 ) W11 (−1 ) + 12 ξ¯ eiwn τc W20 (−1 ) · coskx, coskx ⎟ ⎜+ f02 ξ e ⎟ ⎜ 1  (2) (2 ) ¯ iwn τc (2 ) ⎝ + 8 f30 + f21 (ξ e + 2ξ e−iwn τc ) + f12 (2ξ ξ¯ + ξ 2 e−2iwn τc ) ⎠   (2 ) 2 ¯ + f03 ξ ξ · coskx, coskx

where





Wi(jm) (θ ), coskx =

1

" lπ

lπ

0

Wi(jm) (θ )(x )coskxdx,

and i, j = 0, 1, 2, . . . , m = 1, 2. # lπ # lπ And 0 (coskx )3 dx = 0, for ∀n ∈ N, 0 (coskx )3 dx = l π , for n = 0. Let (ψ1 , ψ2 ) = 1 − i2 , then we can get

g20 =

g11 =

⎧ ⎪ ⎪ ⎨ τ  c 4 ⎪ ⎪ ⎩ 

⎧ ⎪ ⎪ ⎨

0, (1 )

(1 )

f20 + 2 f11

ξe

−iwn τc

(1 ) 2 −2iwn τc

+ f02

ξ e

(2 ) (2 ) −iwn τc (2 ) 2 + f20 + 2 f11 ξe + f02 ξ e

τc 

n ∈ N,



ψ1  ψ2 ,

 −2iw τ n c

0, (1 ) (1 ) ¯ iwn τc (1 ) ¯ f20 + f11 (ξ e + ξ e−iwn τc ) + 2 f02 ξξ



n = 0, n ∈ N,

ψ1   (2 ) (2 ) ¯ iwn τc (2 ) ¯ + f20 + f11 (ξ e + ξ e−iwn τc ) + 2 f02 ξ ξ ψ2 , n = 0,

4 ⎪ ⎪ ⎩ 

g02 = g20 ,

⎛







(1 ) (1 ) (1 ) f20 W11 (0 ) + 12 W20 (0 ) · coskx, coskx

⎞

 (2)   ⎟ ⎜ (1 ) (2 ) + f11 W11 (−1 ) + 12 W20 (−1 ) · coskx, coskx ⎟ ⎜ ⎟ ⎜ ⎜ + f (1) ξ e−iwn τc W (1) (0 ) + 1 ξ¯ eiwn τc W (1) (0 ) · coskx, coskx ⎟ ⎟ ⎜ 2 11 11 20 ⎟ g21 = τc ⎜ ⎜+ f (1) ξ e−iwn τc W (2) (−1 ) + 1 ξ¯ eiwn τc W (2) (−1 ) · coskx, coskx⎟ψ1 2 11 20 ⎟ ⎜ 02 ⎟ ⎜ 1  (1) (1 ) ¯ iwn τc (1 ) ⎜ + 8 f30 + f21 (ξ e + 2ξ e−iwn τc ) + f12 (2ξ ξ¯ + ξ 2 e−2iwn τc ) ⎟ ⎠ ⎝   (1 ) 2 ¯ + f03 ξ ξ · coskx, coskx ⎛





(2 ) (1 ) (1 ) f20 W11 (0 ) + 12 W20 (0 ) · coskx, coskx



⎞

 (2)   ⎟ ⎜ (2 ) (2 ) ⎟ ⎜ + f11 W11 (−1 ) + 12 W20 (−1 ) · coskx, coskx ⎟ ⎜  −iw τ (1)   ⎟ ⎜ (2 ) ( 1 ) 1 ⎜ + f11 ξ e n c W11 (0 ) + 2 ξ¯ eiwn τc W20 (0 ) · coskx, coskx ⎟ ⎟ + τc ⎜ ⎜+ f (2) ξ e−iwn τc W (2) (−1 ) + 1 ξ¯ eiwn τc W (2) (−1 ) · coskx, coskx⎟ψ2 , n ∈ {0, N}. ⎟ ⎜ 02 2 11 20 ⎟ ⎜  ⎜ + 1 f (2) + f (2) (ξ¯ eiwn τc + 2ξ e−iwn τc ) + f (2) (2ξ ξ¯ + ξ 2 e−2iwn τc ) ⎟ 21 12 ⎠ ⎝ 8 30   (2 ) 2 ¯ + f03 ξ ξ · coskx, coskx

285

286

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

For g21 contains W20 (θ ) and W11 (θ ), so it is necessary to compute them. From (38), we can deduce

W˙ (z, z¯ ) = W20 zz˙ + W11 z˙ z¯ + W11 zz¯˙ + W02 z¯z¯˙ + . . . , A(τc )W = A(τc )W20

z2 z¯2 + A(τc )W11 zz¯ + A(τc )W02 . 2 2

(42)

(43)

Besides, from [34],

W˙ = A(τc )W + H (z, z¯ ),

(44)

and

H (z, z¯ ) = H20

z2 z¯2 + H11 zz¯ + H02 + . . . = X0 f (Ut , 0 ) − ( , X0 f (Ut , 0 ), fn  ) · fn , 2 2

(45)

with Hij ∈ P∗ , i, j = 0, 1, 2 . . .. Therefore, from (39) and (41)–(45), we can obtain the following form

(2iwn τc − A(τc ))W20 = H20 ,

(46)

− A(τc )W11 = H11 . Since A(τ c ) has only two eigenvalues ±iwn τ c , then (46) has unique solution Wij given by

W20 = (2iwn τc − A(τc ))−1 H20 ,

(47)

W11 = −A(τc )−1 H11 . From (45), for −1 ≤ θ < 0,

H (z, z¯ ) = −(θ )(0 ) f (Ut , 0 ), fn  · fn

 q (θ ) + q (θ ) q (θ ) − q (θ )  1 2 1 2

=−

,

2

2i

(1 (0 ), 2 (0 )) ·  f (Ut , 0 ), fn  · fn

z2 1 1 = − (q1 (θ )g20 + q2 (θ )g02 ) · fn − (q1 (θ )g11 + q2 (θ )g11 ) · fn zz¯ + . . .. 2 2 2 Thus, for −1 ≤ θ < 0,

H20 (θ ) =

H11 (θ ) =

⎧ ⎨

0,

n ∈ N,

⎩ − 1 (q (θ )g + q (θ )g ) · f , n = 0, 1 20 2 02 0 2 ⎧ 0, n ∈ N, ⎨

(48)

⎩ − 1 (q (θ )g + q (θ )g ) · f , n = 0. 1 11 2 11 0

(49)

2

For θ = 0, H (z, z¯ )(0 ) = f (Ut , 0 ) − ( ,  f (Ut , 0 ), fn  ) · fn , then

% ⎧ $ (1 ) (1 ) (1 ) ⎪ τc f20 + 2 f11 ξ e−iwn τc + f02 ξ 2 e−2iwn τc ⎪ ⎪ · (coskx )2 , n ∈ N, ⎪ (2 ) (2 ) −iwn τc (2 ) 2 −2iwn τc ⎪ ⎨ 4 f20 + 2 f11 ξe + f02 ξ e H20 (0 ) = $ (1 ) % ⎪ (1 ) (1 ) ⎪ ⎪ 1 τc f20 + 2 f11 ξ e−iwn τc + f02 ξ 2 e−2iwn τc ⎪ ⎪ − (q1 (0 )g20 + q2 (0 )g02 ) · f0 , n = 0, ⎩ 4 (2 ) (2 ) −iwn τc (2 ) 2 −2iwn τc 2 f20 + 2 f11 ξ e + f02 ξ e % $ ⎧ (1 ) (1 ) ¯ iwn τc (1 ) ¯ + ξ e−iwn τc ) + 2 f02 ξξ ⎪ ⎪ τc f20 + f11 (ξ e ⎪ · (coskx )2 , n ∈ N, ⎪ ( 2) ( 2 ) ¯ iwn τc ( 2) ¯ ⎪ 4 −iw τ n c ⎨ f20 + f11 (ξ e + ξe ) + 2 f02 ξ ξ H11 (0 ) = % $ (1 ) ⎪ (1 ) ¯ iwn τc (1 ) ¯ −iwn τc ⎪ f + f ( ξ e + ξ e ) + 2 f ξ ξ ⎪ 1 τ 20 11 02 c ⎪ ⎪ − (q1 (0 )g11 + q2 (0 )g11 ) · f0 , n = 0. ⎩ 4 (2 ) (2 ) ¯ iwn τc (2 ) ¯ 2 −iwn τc f20 + f11 (ξ e + ξe ) + 2 f02 ξ ξ

(50)

(51)

By the definition of infinitesimal generator A(τ c ), then (46) is transformed

W˙ 20 (θ ) = 2iwn τcW20 (θ ) +

1 (q1 (θ )g20 + q2 (θ )g02 ) · f0 , 2

(52)

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

287

and −1 ≤ θ < 0. From q1 (θ ) = q1 (0 )eiwn τc θ , −1 ≤ θ ≤ 0, we have





1 ig20 ig02 q1 (θ ) + q2 (θ ) · fn + C1 e2iwn τc θ , 2 w n τc 3wn τc

W20 (θ ) =

(53)

further we also get

⎧ ⎨

W20 (0 ), n ∈ N,   C1 = ⎩ W20 (0 ) − 1 ig20 q1 (0 ) + ig02 q2 (0 ) · f0 , n = 0. 2 w n τc 3wn τc

(54)

For n = 0, θ = 0, by using the definition of A(τ c ) and (53), the first equation of (46) is transformed

 1  ig   ig02 20 q1 (0 ) + q2 (0 ) · f0 + C1 2 w n τc 3wn τc     ig02 ∂ 2 1 ig20 − τc D 2 q1 (0 ) + q2 (0 ) · f0 + C1 3wn τc ∂ x 2 w n τc  1  ig   ig02 20 − L ( τc ) q1 (θ ) + q2 (θ ) · f0 + C1 e2iwn τc θ 2 w n τc 3wn τc $ (1 ) % (1 ) (1 ) 1 τc f20 + 2 f11 ξ e−iwn τc + f02 ξ 2 e−2iwn τc = − (q1 (0 )g20 + q2 (0 )g02 ) · f0 . 4 f (2) + 2 f (2) ξ e−iwn τc + f (2) ξ 2 e−2iwn τc 2 20 11 02

2iwn τc

Therefore, we can deduce

$ (1 ) % (1 ) (1 ) ∂2 τc f20 + 2 f11 ξ e−iwn τc + f02 ξ 2 e−2iwn τc 2iwn τc θ 2iwn τcC1 − τc D 2 C1 − L(τc )(C1 e )= (coskx )2 . 4 f (2) + 2 f (2) ξ e−iwn τc + f (2) ξ 2 e−2iwn τc ∂x 20

We can see from the above expression,



C1 = where

1 2iwn + d1 k2 − a11 −a21 4

$ J1 =

−a12 e−2iwn τc 2iwn + d2 k2 − a22 e−2iwn τc

(1 ) (1 ) −iwn τc (1 ) 2 −2iwn τc f20 + 2 f11 ξe + f02 ξ e

11

(55)

02

−1

J1 ,

(56)

%

(2 ) (2 ) −iwn τc (2 ) 2 −2iwn τc f20 + 2 f11 ξe + f02 ξ e

(coskx )2 .

Similar to the method of solving W20 (θ ), for −1 ≤ θ < 0,

W˙ 11 (θ ) =

1 (q1 (θ )g11 + q2 (θ )g11 ) · f0 , 2

1 W11 (θ ) = 2



ig11 −ig11 q1 (θ ) + q2 (θ ) w n τc w n τc

further we obtain



C2 = where

2 1 d1 k − a11 4 −a21

$ J2 =

−a12

(57)

· fn + C2 ,

−1

d2 k2 − a22

J2 ,

(59)

(1 ) (1 ) ¯ iwn τc (1 ) ¯ f20 + f11 (ξ e + ξ e−iwn τc ) + 2 f02 ξξ (2 )

(2 )

f20 + f11

(58)

%

(2 ) ¯ (ξ¯ eiwn τc + ξ e−iwn τc ) + 2 f02 ξξ

(coskx )2 .

By the above computation, we can get the expression of g21 . Further we can determine the properties of Hopf bifurcating period solutions at the threshold value τ c by computing the following values:

c1 (0 ) =

i 2wn τc



g20 g11 − 2|g11 |2 −

1 |g02 |2 3



+

1 g21 , 2

288 P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291 Fig. 1. The initial condition is that S (i, j ) = 0.30 + 0.001 sin(ix + jt ) and I (i, j ) = 0.27 + 0.001 sin(ix + jt ). (a) and (b) when τ = 1.7, the constant steady state E∗ is asymptotically stable, namely the epidemic becomes the endemic; (c) and (d) when τ = 2.2, the spatially inhomogeneous bifurcating periodic solutions are asymptotically stable, namely the epidemic outbreak recurrently.

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291 Fig. 2. S (i, j ) = 0.30 + 0.001 sin(ix + jt ) and I (i, j ) = 0.27 + 0.001 sin(ix + jt ) are initial conditions. (a) and (b) when τ = 1.4, the constant steady state E∗ is asymptotically stable, that is, the disease develops the endemic in this region; (c) and (d): when τ = 2.1, the bifurcating periodic solutions are spatially heterogeneous and stable, which describes that the disease occurs periodically.

289

290

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

μ2 = −

Re(c1 (0 )) Re(λ (τnj ))

,

(60)

β2 = 2Re(c1 (0 )), T2 = −

1  (Im(c1 (0 )) + μ2 Im(λ (τnj ))). w n τc

If μ2 > 0(μ2 < 0), then the direction of the Hopf bifurcation is forward(backward) which corresponds to the existence of the bifurcating period solutions for τ > τ c (τ < τ c ); the stability of the bifurcating period solutions on center manifold are determined by β 2 , that is, the bifurcating period solutions are stable (unstable) if β 2 < 0(β 2 > 0); moreover, T2 determines the period of the bifurcating period solutions, if T2 < 0(T2 > 0), then the period decreases (increases). 4. Simulation In this section, we perform a series of simulations of the spatial epidemic model with time delay (3) in one-dimensional space. The reaction–diffusion system is solved in a discrete domain. The Laplacian describing diffusion is approximated by using finite differences, and we also discrete the time evolution. Thus, we can solve the numerical solutions of system (3) by using Matlab. Moreover, we find that the Hopf bifurcating periodic solutions with space is spatially inhomogeneous, which is in accord with the previous conclusions [34]. For n = 80, we fix l = 50, d1 = 0.01, d2 = 0.25, R0 = 1.9, Rd = 5, v = 0.3, then the equilibrium is E ∗ = (S∗ , I∗ ) = 0 = 1.7590, c (0 ) = −5.1344 − 3.9002i, and λ (τ 0 ) = 0.2484 + 0.0972i, by the for(0.3047, 0.2742 ). Further, we obtain τ80 1 80 mulae derived in Section 3, μ2 > 0, β 2 < 0 and T2 > 0 are obtain. E∗ is asymptotically stable for 0 ≤ τ < τ c . Next, E∗ loses its stability and Hopf bifurcation occur at τ c . μ2 > 0 shows that the direction of bifurcation is τ > τ c , these bifurcating periodic solutions are stable which are presented in Fig. 1. For n = 100, setting l = 50, d1 = 0.008, d2 = 0.16, R0 = 1.9, Rd = 5, v = 0.3, then the equilibrium is E ∗ = (S∗ , I∗ ) = 0 = 1.7574, c (0 ) = −5.1865 − 4.0337i and (0.3047, 0.2742 ). Moreover, by utilizing the formulae derived in Section 3, we get τ100 1  0 λ (τ100 ) = 0.2481 + 0.0984i. Calculating the formulae (60), μ2 > 0, β 2 < 0 and T2 > 0 are obtain. Therefore, the direction of Hopf bifurcation is τ > τ c , these bifurcating periodic solutions are stable which are showed in Fig. 2. 5. Discussion and conclusion In a one-dimension spatial diffusion SI model with delay, we derive the characteristic equation at the positive constant steady state E∗ (S∗ , I∗ ), and choose the time delay τ as a bifurcation parameter. Furthermore, we obtain the two classes conditions of the existence of Hopf bifurcation: one is the absence of diffusion n = 0, the other is the presence of diffusion n = n0 ∈ N. When τ j j crosses through the critical value τ c (τ0 or τn0 ), the system (3) occurs the Hopf bifurcation. When n = 0, the system (3) occurs spatially homogeneous Hopf bifurcation. For n = n0 ∈ N, the Hopf bifurcating period solutions of the system (3) is spatially heterogeneous, which is consistent with our theoretical results. Next we get the properties of bifurcating period solutions including direction, stability and period by utilizing the normal formal theory and the center manifold theorem of PFD equations. In the real world, many infectious disease exhibit seasonal fluctuations, such as measles, whooping cough, infuenza, chickenpox, rabies [30,36,37]. Based on the studies of the Hopf bifurcating period solutions, these results may provide some new insights on the seasonal outbreaks of diseases. In the future, we will study whether time delay or diffusion also induce the seasonal outbreaks of diseases combining real data. References [1] F. Brauer, C. Castillo-Chavez, Mathematical models in population biology and epidemioloy, Test Appl. Math. 40 (2001) 171–229. [2] L. Wang, Z. Wang, Y. Zhang, X. Li, How human location-specific contact patterns impact spatial transmission between population? Sci. Rep. 3 (2013) 146801– 146810. [3] L. Wang, X. Li, Spatial epidemiology of networked metapopulation: An overview, Chin. Sci. Bull. 59 (2014) 3511–3522. [4] S. Boccaletti, The structure and dynamics of multilayer networks, Phys. Rep. 544 (2014) 1–122. [5] E. Beretta, Y. Takeuchi, Global stability of an sir epidemic model with time delays, J. Math. Biol. 33 (1995) 250–260. [6] A. Kadder, On the dyanmics of a delayed sir epiemic model with a modified saturated incidence rate, Electronic J.D.E. 1 (2009) 1–7. [7] J. Zhang, Z. Jin, J. Yan, Stability and hopf bifurcation in a delayed competition system, Nonlinear Anal.: T.M.A. 70 (2009) 658–670. [8] A. Abta, A. Kadder, T.H. Alaoui, Global stability for delay sir and seir epidemic models with saturated incidece rates, Electronic J.D.E. 2012 (2012) 1–13. [9] Z. Ma, Y. Zhou, W. Wang, Mathematical Models and Dynamics of Infectious Diseases, China Sciences Press, Beijing, 2004. [10] J. Li, G.-Q. Sun, Z. Jin, Pattern formation of an epidemic model with time delay, Phys. A 403 (2014) 100–109. [11] G.-Q. Sun, A. Chakrabort, Q.-X. Liu, Z. Jin, K.E. Anderson, B.-L. Li, Influence of time delay and nonlinear diffusion on herbivore outbreak, Commun. Nonlinear Sci. Numer. Simulat. 19 (2014) 1507–1518. [12] G. Sun, Z. Jin, Q.-X. Liu, L. Li, Pattern formation in a s-i model with nonlinear incidence rates, J. Stat. Mech. 11 (2007) P11011. [13] G.-Q. Sun, Z. Jin, Q.-X. Liu, L. Li, Chaos induced by breakup of waves in a spatial epidemic model with nonlinear incidence rate, J. Stat. Mech. 8 (2008) P08011. [14] K. Wang, W. Wang, S. Song, Dynamics of an hbv model with diffusion and delay, J. Theor. Biol. 253 (2008) 36–44. [15] R. Xu, Z. Ma, An hbv model with diffusion and time delay, J. Theor. Biol. 257 (2009) 499–509. [16] G.-Q. Sun, Q.-X. Liu, Z. Jin, A. Chakraborty, B.-L. Li, Influence of infection rate and migration on extinction of disease in spatial epidemics, J. Theor. Biol. 264 (2010) 95–103. [17] G.-Q. Sun, Pattern formation of an epidemic model with diffusion, Nonlinear Dynam. 69 (2012) 1097–1104. [18] C. Xia, L. Wang, S. Sun, J. Wang, An sir model with infection delay and propagation vector in complex networks, Nonlinear Dynam. 69 (2012) 927–934.

P.-P. Liu / Applied Mathematics and Computation 265 (2015) 275–291

291

[19] J. Sanz, C.-Y. Xia, S. Meloni, Y. Moreno, Dynamics of interacting diseases, Phys. Rev. X 4 (2014) 041005. [20] V. Capasso, L. Maddalena, A nonlinear diffusion system modeling the spread of oro-faecal diseases, in: V. Lakshmikanthan (Ed.), Nonlinear phenomena in mathematical science, Academic Press, New York, 1981. [21] V. Capasso, L. Maddalena, Convergence to equilibrium states for a reaction-diffusion system modeling the spatial spread of a class of bacterial and viral diseases, J. Math. Biol. 13 (1981b) 173–184. [22] H.R. Thieme, X.-Q. Zhao, Asymptotic speeds of spread and traveling waves for integral equations and delayed rection-diffusion models, J. Differ. Equat. 195 (2003) 430–470. [23] K.I. Kim, Z. Lin, Asymptotic behavior of an sei epidemic model with diffusion, Math. Comput. Model. 47 (2008) 1314–1322. [24] F. Berezovsky, G. Karev, B. Song, A simple epidemic model with surprising dynamics, Math. Biosci. Eng. 2 (2005) 133–152. [25] W. Wang, Y. Cai, Complex dynamics of a reaction-diffusion epidemic model, Nonlinear Anal. RWA 13 (2012) 2240–2258. [26] S. Busenberg, K.L. Cooke, The effect of integral conditions in certain equations modelling epidmic and population growth, J. Math. Biol. 10 (1980) 13–32. [27] D. Klinkenberg, H. Nishiura, The correlation between infectivity and incubation period of measles, estimated from households with two cases, J. Theor. Biol. 284 (2011) 52–60. [28] M.C. White, R.W. Nelson, L.M. Kawamura, Changes in characteristics of inmates with latent tuberculosis infection, Publ. Health 126 (2012) 752–759. [29] J. Zhang, Z. Jin, G.-Q. Sun, Analysis of rabies in china: Transmission dynamics and control, PLoS ONE 6 (2011) e20891. [30] J. Zhang, Z. Jin, G.-Q. Sun, X. Sun, S. Ruan, Modeling seasonal rabies epidemics in china, Bull. Math. Biol. 74 (2012) 1226–1251. [31] R. Xu, Global stability of a delayed epidemic model with latent period and vaccination strategy, Appl. Math. Model. 36 (2012) 5293–5300. [32] R. Xu, Global stability of an hiv-1 infection model with saturation infection and intracellular delay, J. Math. Anal. Appl. 375 (2011) 75–81. [33] T. Faria, L.T. Magalhaes, Normal form for retared functional differential equations with parameters and applications to hopf bifurcation, J. Differ. Equat. 22 (1995) 181–200. [34] J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. [35] J. Hale, Theory of Functional Differential Equations, Springer-Verlag, Berlin, 1977. [36] W.P. London, J.A. Yorke, Recurrent outbreaks of measles, chickenpox and mumps. I. Seasonal variation in contact rates, Amer. J. Epidem. 98 (1973) 453–468. [37] S.F. Dowell, Seasonal variation in host susceptibility and cycles of certain infectious diseases, Emerg. Infect. Dis. 7 (2001) 369–374.