# Bifurcation in a Discrete Competition System

Authors: Li Xu; Lianjun Zou; Zhongxiang Chang; Shanshan Lou; Xiangwei Peng; Guang Zhang
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

A new difference system is induced from a differential competition system by different discrete methods. We give theoretical analysis for local bifurcation of the fixed points and derive the conditions under which the local bifurcations such as flip occur at the fixed points. Furthermore, one- and two-dimensional diffusion systems are given when diffusion terms are added. We provide the Turing instability conditions by linearization method and inner product technique for the diffusion system with periodic boundary conditions. A series of numerical simulations are performed that not only verify the theoretical analysis, but also display some interesting dynamics.


---

## Body

## 1. Introduction

Interactions of different species may take many forms such as competition, predation, parasitism, and mutualism. One of the most important interactions is the competition relationship. The dynamic relationship between the two competition species is one of the dominant subjects in mathematical ecology due to its universal existence and importance. Lotka-Volterra competition systems are ecological models that describe the interaction among various competing species and have been extensively investigated in recent years (see [1–3] and the references therein). In the earlier literature, the two-competing species competition models are often formulated in the form of ordinary differential systems as follows:
(1)u′(t)=u(t)(r1-a11u(t)-a12v(t)),v′(t)=v(t)(r2-a21v(t)-a22v(t)),
for t∈[0,+∞)aij≥0,i, j=1,2, where u(t) and v(t) are the quantities of the two species at time t, r1>0 and r2>0 are growth rates of the respective species, a11 and a22 represent the strength of the intraspecific competition, and a12 and a21 represent the strength of the interspecific competition.The discrete time models governed by difference equation are more realistic than the continuous ones when the populations have nonoverlapping generations or the population statistics are compiled from given time intervals and not continuously. Moreover, since the discrete time models can also provide efficient computational models of continuous models for numerical simulations, it is reasonable to study discrete time models governed by difference equations.Applying forward Euler scheme to the first equation of system (1) and obtaining a discrete analog of the second equation by considering a variation with piecewise constant arguments for certain terms on the right side (exponential discrete form) [4], we obtain the following equation:
(2)ut+1=ut(r1-a11ut-a12vt),vt+1=vtexp⁡(r2-a21ut-a22vt).For the sake of simplicity, let(3)ut+1=ut(r1+1-a11ut-a12vt),vt+1=vtexp⁡(r2-a21ut-a22vt).By settingUt=a11ut and Vt=a22vt and r1=r2, a11=a22, and a12=a21, we have the following form:
(4)ut+1=ut(r+1-ut-avt)=f(ut,vt),vt+1=vtexp⁡(r-aut-vt)=g(ut,vt).Although numerical variations of system (1) have been extensively studied (see, e.g., the work in [5–8]), some discrete analogs may be found in [9–15], regarding attractivity, persistence, global stabilities of equilibrium, and other dynamics. Up to now, to the best of our knowledge, the discrete system (4) has not been investigated.Since the pioneering theoretical works of Skellam [15] and Turing [16], many works have focused on the effect of spatial factors which play a crucial role in the stability of populations [17–19]. Many important epidemiological and ecological phenomena are strongly influenced by spatial heterogeneities because of the localized nature of transmission or other forms of interaction. Thus, spatial models are more suitable for describing the process of population development. Impact of spatial component on system has been widely investigated (e.g., see [20–22]). It may be a case in reality that the motion of individuals is random and isotropic; that is, without any preferred direction, the individuals are also absolute ones in microscopic sense, and each isolated individual exchanges materials by diffusion with its neighbors [19, 23]. Thus, it is reasonable to consider a 1D or 2D spatially discrete reaction diffusion system to explain the population system. Corresponding to the above analysis, we can obtain the following one-dimensional diffusion systems:
(5)uit+1=uit(r+1-uit-avit)+d1∇2uit,vit+1=vitexp⁡(r-auit-vit)+d2∇2vit,
for i∈{1,2,…,m}=[1,m], t∈Z+ and ∇2uit=ui+1t-2uit+ui-1t, ∇2vit=vi+1t-2vit+vi-1t, and two-dimensional diffusion systems:
(6)uijt+1=uijt(r+1-uijt-avijt)+d1∇2uijt,vijt+1=vijtexp⁡(r-auijt-vijt)+d2∇2vijt,
for i,j∈{1,2,…,m}=[1,m], t∈R+=[0,∞) and ∇2uijt=ui+1,jt+ui,j+1t+ui-1,jt+ui,j-1t-4uijt, ∇2vijt=vi+1,jt+vi,j+1t+vi-1,jt+vi,j-1t-4vijt.In this paper, we will study the dynamical behaviors of models (4), (5), and (6). By using the theory of difference equation, the theory of bifurcation, and the center manifold theorem we will establish the series of criteria on the existence and local stability of equilibria, flip bifurcation for the system (4). For the one- or two-dimensional diffusion systems, with periodic boundary conditions, the Turing instability (or Turing bifurcation) theory analysis will be given. Turing instability conditions can then be deduced combining linearization method and inner product technique. Furthermore, by means of the numerical simulations method, we will indicate the correctness and rationality of our results.The paper is organized as follows. In Section2, we study the existence and stability of equilibria points and the conditions of existence for flip bifurcation are verified for system (4). Turing instability conditions will be illustrated by linearization method and inner product technique for the system (5) and (6) with periodic boundary conditions in Section 3. A series of numerical simulations are performed that not only verify the theoretical analysis, but also display some interesting dynamics. For the system (4), the bifurcation diagrams are given. The impact of the system parameters and diffusion coefficients on patterns can also be observed visually for the given diffusion systems. Finally, some conclusions are given.


## 2. Analysis of Equilibria and Flip Bifurcation

Clearly, the system (4) has four possible steady states; that is, E0=(0,0), exclusion points E1=(r,0), E2=(0,r), and nontrivial coexistence point E3=(u*,v*), where
(7)u*=v*=ra+1.The linearized form of (4) is then
(8)ut+1=fuut+fvvt,vt+1=guut+gvvt,
which has the Jacobian matrix
(9)JEi=[fufvgugv]Pi=[1-u-au-av1-v]Pi,i=0,1,2,3.The characteristic equation of the Jacobian matrixJ can be written as
(10)λ2+pλ+q=0,
where p=-(fu+gv) and q=fugv-fvgu.In order to discuss the stability of the fixed points of (4), we also need the following definitions [20]:(1)
if|λ1|<1 and |λ2|<1, then steady state E is called a sink and E is locally asymptotically stable;(2)
if|λ1|>1 and |λ2|>1, then E is called a source and E is unstable;(3)
if|λ1|>1 and |λ2|<1 (or |λ1|<1 and |λ2|>1), then E is called a saddle;(4)
if either|λ1|=1 and |λ2|≠1 or |λ2|=1 and |λ1|≠1, then E is called nonhyperbolic.Case 1 (the fixed pointE0=(0,0)).
The linearization of (4) about E0 has the Jacobian matrix
(11)JE0=[1001],
which has two eigenvalues
(12)λ1=λ2=1.
The fact means that the system is resonance at the fixed point P0.Case 2 (the fixed pointE1=(0,r)).
At the fixed point, the Jacobian matrix has the form(13)JE1=[10-ar1-r],
and the corresponding eigenvalues of (13) are
(14)λ1=1,λ2=1-r.r is a bifurcation parameter. And r≠2 implies λ2≠-1, and the fixed point E1 is nonhyperbolic.Case 3 (the third fixed pointE2=(r,0)).
The linearization of (4) about E2 has the Jacobian matrix
(15)JE1=[1-r-ar01],
and the eigenvalues of (15) are
(16)λ1=1-r,λ2=1.r is a bifurcation parameter. And r≠2 implies λ1≠-1, and the fixed point E2 is nonhyperbolic.Case 4 (the fixed pointE3=(r/(a+1),r/(a+1))).
The linearization of (4) about E3 has the Jacobian matrix
(17)JE3=[1-ra+1-ara+1-ara+11-ra+1],
and the eigenvalues of (17) are
(18)λ1=1-r,λ2=1+(a-1)ra+1;
then, we have the following results:(1)
|λ1|<1, |λ2|<1 if and only if 0<r<2, 0<a<1;(2)
λ1=-1, |λ2|≠1 if and only if r=2, a≠1;(3)
λ2=1, |λ1|≠1 if and only if a=2, r≠2;(4)
λ2=-1, |λ1|≠1 if and only if a=(r-2)/(r+2), r≠2;(5)
|λ1|<1, |λ2|>1 if and only if 0<r<2, a>1;(6)
|λ1|>1, |λ2|<1 if and only if 2<r<(a+1)/(a-1);(7)
|λ1|>1, |λ2|>1 if and only if r>2, a>1.
The following theorem is the case that the fixed pointE3 is a flip bifurcation point.Theorem 1.
The positive fixed pointE3 undergoes a flip bifurcation at the threshold rF=2.Proof.
Letζn=un-u*, ηn=vn-v*, and μn=r-2, and parameter μn is a new and dependent variable; the system (4) becomes
(19)(ζn+1ηn+1μn+1)=((ζn+μn+2a+1)(1-ζn-aηn)-μn+2a+1(ηn+μn+2a+1)exp⁡(-aζn-ηn)-μn+2a+1μn).
Let(20)T=[1-10110001];
then
(21)T-1=[12120-12120001].
By the following transformation:(22)(ζnηnμn)=T(xnynδn),
the system (19) can be changed into
(23)(xn+1yn+1δn+1)=(-10003a-1a+10001)(xnynδn)+(f(xn,yn,δn)g(xn,yn,δn)0),
where
(24)f(xn,yn,δn)=-(a+1)xn2-2(a+1)xnyn-xnδn+(a-1+(a-1)2a+1)yn2+o((|xn|+|yn|+|δn|)3),g(xn,yn,δn)=(a+1)xn2-2(a+1)xnyn-2(a-1)a+1ynδn+(a-1)2a+1yn2+o((|xn|+|yn|+|δn|)3).
Then, we can consider
(25)yn=h(xn,δn)=a1xn2+a2xnδn+a3δn2+o((|xn|+|δn|)3),
which must satisfy
(26)h(-xn+f(xn,yn,δn),δn+1)=3a-1a+1h(xn,δn)+(a+1)xn2-2(a+1)xnh(xn,δn)-2(a-1)a+1h(xn,δn)δn+(a-1)2a+1h2(xn,δn)+o((|xn|+|yn|+|δn|)3).
By calculating, we can get that
(27)a1=(a+1)22(1-a),a2=0,a3=0,
and the system (19) is restricted to the center manifold, which is given by
(28)f:xn+1=-xn+a3+4a2+3aa-1xn2+(a+1)xnδn+14(a+1)3xn4+o(|xn|4).
Since
(29)(∂f∂δ∂2f∂x2+2∂2f∂x∂δ)|(0,0)=2(a+1)≠0,(12(∂2f∂x2)2+13∂3f∂x3)|(0,0)=2(a3+4a2+3aa-1)2>0,
system (4) undergoes a flip bifurcation at E3. The proof is completed.


## 3. Turing Bifurcation

In this section, we discuss the Turing bifurcation. Turing's theory shows that diffusion could destabilize an otherwise stable equilibrium of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of instability is usually called Turing instability or diffusion-driven instability [16].


### 3.1. One-Dimensional Case

We consider the following diffusion system:(30)uit+1=uit(r+1-uit-avit)+d1∇2uit,vit+1=vitexp⁡(r-auit-vit)+d2∇2vit,
with the periodic boundary conditions
(31)u0t=umt,u1t=um+1t,v0t=vmt,v1t=vm+1t,
for i∈{1,2,…,m}=[1,m] and t∈Z+, where m is a positive integer,
(32)∇2uit=ui+1t-2uit+ui-1t,∇2vit=vi+1t-2vit+vi-1t.In order to study Turing instability of (30) and (31), we firstly consider eigenvalues of the following equation:
(33)∇2Xi+λXi=0
with the periodic boundary conditions
(34)X0=Xm,X1=Xm+1.By calculating, the eigenvalue problem (33)-(34) has the eigenvalues
(35)λs=4sin2(s-1)πmfors∈[1,m].We linearise at the steady state, to get(36)wit+1=Jwit+D∇2wit,D=(d100d2)
with the periodic boundary conditions
(37)w0t=wmt,w1t=wm+1t,
where
(38)wit=(uit-u*vit-v*)=(xityit).Then, respectively, taking the inner product of (36) with the corresponding eigenfunction Xsi of the eigenvalue λs, we see that
(39)∑i=1m‍Xsixit+1=fu∑i=1m‍Xsixit+fv∑i=1m‍Xsiyit+d1∑i=1m‍Xsi∇2xit,∑i=1m‍Xsiyit+1=gu∑i=1m‍Xsixit+gv∑i=1m‍Xsiyit+d2∑i=1m‍Xsi∇2yit.
Let Ut=∑i=1m‍Xsixit and Vt=∑i=1m‍Xsiyit and use the periodic boundary conditions (34) and (37); then we have
(40)Ut+1=fuUt+fvVt-d1λsUt,Vt+1=guUt+gvVt-d2λsVt
or
(41)Ut+1=(fu-d1λs)Ut+fvVt,Vt+1=guUt+(gv-d2λs)Vt.Thus, the following fact can be obtained.Proposition 2.
If(uit,vit) is a solution of the problem of (30) and (31), then
(42)(Ut=∑i=1m‍Xsixit,Vt=∑i=1m‍Xsiyit)
is a solution of (41), where λs is some eigenvalue of (33)-(34) and Xsi is the corresponding eigenfunction. For some eigenvalue λs of (33)-(34), if (Ut,Vt) is a solution of the system (41), then
(43)(uit=UtXsi,vit=VtXsi)
is a solution of (30) with the periodic boundary conditions (31).Proposition 3.
If there exist positive numbersd1, d2 and the eigenvalue λs of the problem (33)-(34) such that one of the conditions
(44)h(λs)<(λs(d1+d2)-(fu*+gv*))-1,h(λs)<-(λs(d1+d2)-(fu*+gv*))-1
or
(45)h(λs)>1
holds, then the problem (30) and (31) at the fixed point (u*,v*) is unstable, where
(46)h(λs)=d1d2λs2-(d1gv*+d2fu*)λs+(fu*gv*-fv*gu*).For the system (5), we have the following results about instability of the positive equilibrium of system.Theorem 4.
0
<
r
<
2, 0<a<1 and Proposition 3 mean or show that the problem (30) and (31) is diffusion-driven unstable or Turing unstable.


### 3.2. Two-Dimensional Case

In this subsection, we will pay our attention to the Turing instability analysis for the following two-dimensional system:(47)uijt+1=uijt(r+1-uijt-avijt)+d1∇2uijt,vijt+1=vijtexp⁡(r-auijt-vijt)+d2∇2vijt,
with the periodic boundary conditions
(48)ui,0t=ui,mt,ui,1t=ui,m+1t,u0,jt=um,jt,u1,jt=um+1,jt,vi,0t=vi,mt,vi,1t=vi,m+1t,v0,jt=vm,jt,v1,jt=vm+1,jt,
for i,j∈{1,2,…,m}=[1,m] and t∈Z+, where m is a positive integer,
(49)∇2uijt=ui+1,jt+ui,j+1t+ui-1,jt+ui,j-1t-4uijt,∇2vijt=vi+1,jt+vi,j+1t+vi-1,jt+vi,j-1t-4vijt.The following theorem will show that the system (47) also undergoes Turing instability. Since the analysis is very similar to the one-dimensional case, the proof is omitted.Theorem 5.
If there exist positive numbersd1, d2 and the eigenvalue kls2 of the corresponding characteristic equation such that one of the conditions
(50)h(kls2)<(kls2(d1+d2)-(fu*+gv*))-1,h(kls2)<-(kls2(d1+d2)-(fu*+gv*))-1
or
(51)h(kls2)>1
holds and 0<r<2, 0<a<1 then the problem (47)-(48) at the fixed point (u*,v*) is diffusion-driven unstable or Turing unstable, where
(52)h(kls2)=d1d2kls4-(d1gv*+d2fu*)kls2+(fu*gv*-fv*gu*),kls2=λl,s=4(sin2((l-1)πm)+sin2((s-1)πm))hhhhhhhhhhhhhhhhhhhhhhhhihhforl,s∈[1,m].


## 3.1. One-Dimensional Case

We consider the following diffusion system:(30)uit+1=uit(r+1-uit-avit)+d1∇2uit,vit+1=vitexp⁡(r-auit-vit)+d2∇2vit,
with the periodic boundary conditions
(31)u0t=umt,u1t=um+1t,v0t=vmt,v1t=vm+1t,
for i∈{1,2,…,m}=[1,m] and t∈Z+, where m is a positive integer,
(32)∇2uit=ui+1t-2uit+ui-1t,∇2vit=vi+1t-2vit+vi-1t.In order to study Turing instability of (30) and (31), we firstly consider eigenvalues of the following equation:
(33)∇2Xi+λXi=0
with the periodic boundary conditions
(34)X0=Xm,X1=Xm+1.By calculating, the eigenvalue problem (33)-(34) has the eigenvalues
(35)λs=4sin2(s-1)πmfors∈[1,m].We linearise at the steady state, to get(36)wit+1=Jwit+D∇2wit,D=(d100d2)
with the periodic boundary conditions
(37)w0t=wmt,w1t=wm+1t,
where
(38)wit=(uit-u*vit-v*)=(xityit).Then, respectively, taking the inner product of (36) with the corresponding eigenfunction Xsi of the eigenvalue λs, we see that
(39)∑i=1m‍Xsixit+1=fu∑i=1m‍Xsixit+fv∑i=1m‍Xsiyit+d1∑i=1m‍Xsi∇2xit,∑i=1m‍Xsiyit+1=gu∑i=1m‍Xsixit+gv∑i=1m‍Xsiyit+d2∑i=1m‍Xsi∇2yit.
Let Ut=∑i=1m‍Xsixit and Vt=∑i=1m‍Xsiyit and use the periodic boundary conditions (34) and (37); then we have
(40)Ut+1=fuUt+fvVt-d1λsUt,Vt+1=guUt+gvVt-d2λsVt
or
(41)Ut+1=(fu-d1λs)Ut+fvVt,Vt+1=guUt+(gv-d2λs)Vt.Thus, the following fact can be obtained.Proposition 2.
If(uit,vit) is a solution of the problem of (30) and (31), then
(42)(Ut=∑i=1m‍Xsixit,Vt=∑i=1m‍Xsiyit)
is a solution of (41), where λs is some eigenvalue of (33)-(34) and Xsi is the corresponding eigenfunction. For some eigenvalue λs of (33)-(34), if (Ut,Vt) is a solution of the system (41), then
(43)(uit=UtXsi,vit=VtXsi)
is a solution of (30) with the periodic boundary conditions (31).Proposition 3.
If there exist positive numbersd1, d2 and the eigenvalue λs of the problem (33)-(34) such that one of the conditions
(44)h(λs)<(λs(d1+d2)-(fu*+gv*))-1,h(λs)<-(λs(d1+d2)-(fu*+gv*))-1
or
(45)h(λs)>1
holds, then the problem (30) and (31) at the fixed point (u*,v*) is unstable, where
(46)h(λs)=d1d2λs2-(d1gv*+d2fu*)λs+(fu*gv*-fv*gu*).For the system (5), we have the following results about instability of the positive equilibrium of system.Theorem 4.
0
<
r
<
2, 0<a<1 and Proposition 3 mean or show that the problem (30) and (31) is diffusion-driven unstable or Turing unstable.


## 3.2. Two-Dimensional Case

In this subsection, we will pay our attention to the Turing instability analysis for the following two-dimensional system:(47)uijt+1=uijt(r+1-uijt-avijt)+d1∇2uijt,vijt+1=vijtexp⁡(r-auijt-vijt)+d2∇2vijt,
with the periodic boundary conditions
(48)ui,0t=ui,mt,ui,1t=ui,m+1t,u0,jt=um,jt,u1,jt=um+1,jt,vi,0t=vi,mt,vi,1t=vi,m+1t,v0,jt=vm,jt,v1,jt=vm+1,jt,
for i,j∈{1,2,…,m}=[1,m] and t∈Z+, where m is a positive integer,
(49)∇2uijt=ui+1,jt+ui,j+1t+ui-1,jt+ui,j-1t-4uijt,∇2vijt=vi+1,jt+vi,j+1t+vi-1,jt+vi,j-1t-4vijt.The following theorem will show that the system (47) also undergoes Turing instability. Since the analysis is very similar to the one-dimensional case, the proof is omitted.Theorem 5.
If there exist positive numbersd1, d2 and the eigenvalue kls2 of the corresponding characteristic equation such that one of the conditions
(50)h(kls2)<(kls2(d1+d2)-(fu*+gv*))-1,h(kls2)<-(kls2(d1+d2)-(fu*+gv*))-1
or
(51)h(kls2)>1
holds and 0<r<2, 0<a<1 then the problem (47)-(48) at the fixed point (u*,v*) is diffusion-driven unstable or Turing unstable, where
(52)h(kls2)=d1d2kls4-(d1gv*+d2fu*)kls2+(fu*gv*-fv*gu*),kls2=λl,s=4(sin2((l-1)πm)+sin2((s-1)πm))hhhhhhhhhhhhhhhhhhhhhhhhihhforl,s∈[1,m].


## 4. Numerical Simulation

As is known to all, the bifurcation diagram provides a general view of the evolution process of the dynamical behaviors by plotting a state variable with the abscissa being one parameter. As a parameter varies, the dynamics of the system we concerned change through a local or global bifurcation which leads to the change of stability at the same time.Now,r is considered as a parameter with the range 0.5–3.5 for the system (4). Since the bifurcation diagrams of r-ut are similar to the bifurcation diagrams of r-vt, we will only show the former which can be seen from Figure 1.Figure 1
Bifurcation diagram forr-ut with a=0.511.Next, we performed a series of simulations for the reaction-diffusion systems, and, in each, the initial condition was always a small amplitude random perturbation 1% around the steady state. As a numerical example, we consider the bifurcation of the two-dimensional system (47)-(48). It is well known that Turing instability (bifurcation) is diffusion-driven instability; thus the diffusion rate is vital to the pattern formation. To investigate the effect of diffusion coefficients on patterns, by keeping all the other parameters of the system fixed (a=0.55, r=0.71, and d1=0.22), we change a diffusion coefficient in the Turing instability region (parameter space which satisfies Turing instability). Figures 2(a)–2(f) exhibit in detail the different distribution of patterns with varying values of d. If we let d2=0.19, a stable pattern of square shapes, namely, stationary wave, is observed in Figure 2(a). With the increase of d2, some tips appear, which attempt to form self-centered spiral waves (Figure 2(b)). Then some tips vanish and some tend to evolve to regular spiral waves, shown in Figure 2(c). Then, with the parameter evolution proceeding, the size of spiral waves rises, but the density of them decreases (see Figure 2(d)). If d2 is further increased, we observe that, in parts of patterns, spiral waves begin to break up. Hardly can spiral waves be seen, disorder and chaotic structure are depicted in Figures 2(e) and 2(f).Pattern selection with the increase ofd2 in the Turing instability region when a=0.55, r=0.71, and d1=0.22. (a) d2=0.195. (b) d2=0.20. (c) d2=0.205. (d) d2=0.21. (e) d2=0.215. (f) d2=0.217.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Discussion and Conclusion

In this paper, we have applied different discrete schemes to convert the continuous Lotka-Volterra competition model to a new discrete model and studied the dynamical characteristic of the discrete model. Our theoretical analysis and numerical simulations have demonstrated that the discrete competition model undergoes flip bifurcation. Furthermore, when the effects of spatial factors are considered, we discuss the Turing instability conditions combining linearization method and inner product technique. The impact of the diffusion coefficients on patterns can also be observed visually, and some interesting situations can be observed. Indeed, the new discrete model can result in a rich set of patterns and we expect that it is more effective in practice.



---

*Source: 193143-2014-04-28.xml*