# Local Stability in 3D Discrete Dynamical Systems: Application to a Ricker Competition Model

Authors: Rafael Luís; Elias Rodrigues
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A survey on the conditions of local stability of fixed points of three-dimensional discrete dynamical systems or difference equations is provided. In particular, the techniques for studying the stability of nonhyperbolic fixed points via the centre manifold theorem are presented. A nonlinear model in population dynamics is studied, namely, the Ricker competition model of three species. In addition, a conjecture about the global stability of the nontrivial fixed points of the Ricker competition model is presented.


---

## Body

## 1. Introduction

In the study of a dynamical system, specially of nonlinear models, one of the objectives is to determine the stability conditions of equilibrium points or fixed points.According to Osipenko [1], in his excellent paper titled “Center Manifolds”, Lyapunov [2], in his dissertation “General Problem on Stability of Motion,” was perhaps one of the pioneers in the study of the stability of dynamical systems, creating the foundations of the modern theory of the stability in dynamical systems. He considered a continuous dynamical system of the form x˙=F(x), x∈Rn, where F:Rn→Rn is a Ck-smooth vector field, for some k≥1. Near the zero equilibrium O, the system of differential equations can be written as x˙=Ax+f(x), O={x=0}, A=JF(O), and f(O)=Jf(O)=0, where J is the Jacobian matrix of the function. He proved that O is stable if all eigenvalues of the matrix A have negative real parts and O is unstable if there exists an eigenvalue with positive real part. He also studied the case when some eigenvalues have negative real part while the others have zero real part. Moreover, he showed that there exists a two-dimensional invariant surface M passing in the origin when the matrix A has a pair of pure imaginary eigenvalues and the others have negative real parts, in which the equilibrium O of the system is stable whenever O restricted to M is stable. This is precisely what we today know as “centre manifold.” Therefore, Lyapunov built the foundations of stability analysis, proving the existence of a centre manifold, formulating the reduction principle and giving stability conditions for equilibrium points.In discrete dynamical systems of the formxn+1=G(xn), for n=0,1,2,…, where G:Rn→Rn is a Ck smooth mapping, for some k≥1, which can be written near the origin as x→Ax+f(x), O={x=0}, A=JG(O), and f(O)=Jf(O)=0, the equivalent conclusion to Lyapunov’s one is as follows: the fixed point O is locally stable if all the eigenvalues of the matrix A are less than one in absolute value, and O is unstable if there exists at least one of the eigenvalues of A with absolute value greater than one. The centre manifold appears when some of the eigenvalues are one in absolute value, while the others are less than one in absolute value.We should mention that the reduction principle in the study of the qualitative properties of equilibrium points was completely established in 1964 by Pliss [3]. Three years later, Kelley [4] introduced the modern term of centre manifold and gave a complete proof of the existence of manifolds through equilibrium points in general case.There exists a parallel between the study of the stability of fixed points in continuous systems and that in discrete systems. A natural extension of the results for continuous systems has been deeply investigated for discrete systems. See, for instance, the book from 1979 by Iooss [5], where the author presents a complete list of such results, the book by Carr [6], based on a series of lectures that gives an introduction to the applications of centre manifold theory, and the books by Guckenheimer and Holmes [7] and by Kuznetsov [8], where an exploration of the methods of local bifurcation theory for flows and maps, including centre manifolds and normal forms, is provided.Some pedagogical examples explaining how to check the stability of fixed points via the centre manifold may be found in several books, as is the example of the book by Elaydi [9], where the author studies a two-dimensional quadratic model and a two-dimensional cubic model, the book by Wiggis [10], where a three-dimensional quadratic mapping is studied, and the book by Zhang [11], where a two-dimensional quadratic map is deeply studied.In the end of the past century and in the beginnings of the current century, the theory of discrete dynamical systems has developed quicker than ever. The reason may be in the evolution of the computers, as seen in specific software such as Mathematica or Maple, which facilitate the algebraic computations. Nowadays, some nonlinear models can be studied in a much simpler way. Citing some examples, Dektešević et al. [12] studied a two-dimensional quadratic map with two parameters; Guzowska et al. [13] studied a two-dimensional rational model with four parameters denominated by logistic competition model; Karydas and Schinas [14] gave a new proof of the center manifold theorem for a system of difference equations; Kulenović and Merino [15] studied a competitive planar map and established conditions that guarantee the existence of an invariant curve when both eigenvalues of the Jacobian of the map are nonzero and at least one of them has absolute value less than one. These results can be applied to many hyperbolic and nonhyperbolic cases. Several applications to planar systems of difference equations are provided in this paper. Luis et al. [16] studied a nonlinear planar model with four parameters denominated planar Ricker competition model and Psarros et al. [17] studied a three-dimensional model with six parameters.As a consequence of the complexity of the techniques involved in the computations of the centre manifold, the study of the stability of nonhyperbolic fixed points in higher dimensional nonlinear models is usually avoided. The main objective of this survey is precisely to provide the necessary techniques for determining such set in the case of three-dimensional models and illustrate them with a concrete nonlinear example, as is the case of the Ricker competition model.This survey is organized as follows: in Section2, the principal definitions and results concerning the stability of fixed points in one-dimensional and in two-dimensional discrete models are given. In Section 3, the stability conditions for three-dimensional discrete models using the trace, the determinant, and the sum of the principal minors of the Jacobian matrix are studied. In Section 4, the main tools in order to determine the centre manifold of the system are presented. In Section 5, the stability of the fixed points of the Ricker competition model of three species is studied. Finally, in Section 6, a conjecture is presented that the local stability conditions of the nontrivial fixed points of the Ricker competition model imply its global stability.


## 2. Local Stability in 1D and 2D Systems

In this section we provide a brief review of the basic results from the literature, related to the study of stability of fixed points in one- and in two-dimensional discrete models.


### 2.1. One-Dimensional Maps

Consider an intervalI⊆R and a map f:I→I. A point x∗∈R is said to be a fixed point (or equilibrium point) of f if f(x∗)=x∗, and given x0∈R, we define its orbit O(x0) as the set of points O(x0)={x0,f(x0),f2(x0),f3(x0),…}, where fn=f∘fn-1, for n=1,2,….One of the main objectives of the theory of discrete dynamical systems and, in particular, of the stability theory is the study of the behavior of orbits near fixed points, that is, the behavior of solutions of difference equations when the starting points are near equilibrium points. Hence, a basic definition in this field is needed. LetN denote the set of nonnegative integers.Definition 1 (local stability).
Letf:I→I be a map and x∗ be a fixed point of f, where I is an interval of real numbers. Then(1)
the fixed pointx∗ is said to be locally stable if, for any ϵ>0, there exits δ>0 such that, for all x0∈I with |x0-x∗|<δ, we have |fn(x0)-x∗|<ϵ, for all n∈N. Otherwise, the fixed point x∗ will be called unstable;(2)
the fixed pointx∗ is said to be attracting if there exists η>0 such that |x0-x∗|<η implies limn→∞⁡fn(x0)=x∗;(3)
the fixed pointx∗ is said to be locally asymptotically stable if it is both stable and attracting. If in the previous item η=∞, then x∗ is said to be globally asymptotically stable.One of the most effective graphical iteration methods to depict stability of fixed points is the cobweb diagram (also known as stair-step diagram) since it is possible to observe the behavior of the orbit of an initial point near the fixed point. It consists in a sequence of line segments connecting the points in the orbit of an initial pointx0, that is, a sequence of line segments with end points (x0,x1) and (x1,x1), (x1,x1) and (x1,x2), (x1,x2) and (x2,x2), (x2,x2) and (x2,x3), and so forth. On the cobweb plot, a locally stable fixed point corresponds to an inward spiral, while an unstable fixed point is an outward one. It follows from the definition of a fixed point that these spirals will center at a point where the diagonal y=x line crosses the function graph. This technique may be found in many books on discrete dynamical systems or difference equations as are the books by Alligood et al. [18], by Devaney [19], by Elaydi [9, 20], by Zhang [11], and by others, as well. In Figure 1 we present an illustration of the cobweb diagram for the Ricker equation xn+1=xner-xn, for n=0,1,2,… and r>0.Figure 1
Cobweb animation for the Ricker model given byxn+1=xner-xn. Here we are using r=1.8 and one can observe two orbits, the orbit of an initial point x0=0.2 and the orbit of an initial point x0=2.3. Both orbits converge to the fixed point x∗=1.8 meaning that the fixed point x∗=1.8 is locally stable.However, cobweb diagrams are not the most efficient tool to study local stability. There exists a simple but powerful criterion for knowing the local stability of fixed points. We may divide the fixed points into two categories: hyperbolic and nonhyperbolic. A fixed pointx∗ of a map f is said to be hyperbolic if |f′(x∗)|≠1. Otherwise, it is nonhyperbolic.Theorem 2 (see Elaydi [9], page 25).
Letx∗ be a hyperbolic fixed point of a map f, where f is continuous and differentiable at x∗. The following statements hold true:(1)
If|f′(x∗)|<1, then x∗ is locally asymptotically stable.(2)
Iff′x∗>1, then x∗ is unstable.The stability criteria for nonhyperbolic fixed points are more complex and are summarized in the following theorem (for a complete classification of nonhyperbolic fixed points cf. [9], page 33). Before presenting the criteria we introduce the notion of Schwarzian derivative.Definition 3 (Schwarzian derivative).
The Schwarzian derivative,Sf, of a function f, is defined by (1)Sfx=f′′′xf′x-32f′′xf′x2.In particular, when f′(x∗)=-1, we have Sf(x∗)=-f′′′x∗-(3/2)f′′x∗2.Theorem 4 (see Elaydi [9], pages 28–30).
Letx∗ be a fixed point of a map f and f′, f′′, and f′′′ be continuous at x∗.(1)
Letf′(x∗)=1.
(a)
Iff′′(x∗)>0, then x∗ is unstable but semistable from the left.(b)
Iff′′(x∗)<0, then x∗ is unstable but semistable from the right.(c)
Iff′′(x∗)=0 and f′′′(x∗)>0, then x∗ is unstable.(d)
Iff′′(x∗)=0 and f′′′(x∗)<0, then x∗ is locally asymptotically stable.(2)
Letf′(x∗)=-1.
(a)
IfSf(x∗)<0, then x∗ is locally asymptotically stable.(b)
IfSf(x∗)>0, then x∗ is unstable.


### 2.2. Two-Dimensional Maps

Definition1 (local stability) can now be extended to two-dimensional models (or higher dimensional models as well), using an appropriate norm. Since it will follow the same guidelines, we will omit it here.The local stability of a fixed point(x∗,y∗) of the equation(2)xn+1,yn+1=Fxn,yn,F=f,gis guaranteed when all the eigenvalues of the Jacobian, JF(x∗,y∗), evaluated at the fixed point (x∗,y∗) have magnitude less than one. This is precisely the conclusion of the following result; compare Theorem4.11 in Elaydi [9], page 221.Theorem 5 (see Elaydi [9], page 221).
LetF:D⊂R2→R2 be a C1 map, where D is an open subset of R2, (x∗,y∗) a fixed point of F, and JF=JF(x∗,y∗) the Jacobian matrix of F evaluated at the fixed point (x∗,y∗) and let ρ(JF) be the spectral radius of the matrix JF. Then the following statements hold true: (1)
Ifρ(JF)<1, then (x∗,y∗) is locally asymptotically stable.(2)
Ifρ(JF)>1, then (x∗,y∗) is unstable.(3)
Ifρ(JF)=1, then (x∗,y∗) may or may not be stable.The necessary and sufficient conditions can be determined by using the trace,tr, and the determinant, det, of JF as is summarized in the following result; compare Theorem4.4 in Elaydi [9], page 200:Theorem 6 (see Elaydi [9], page 200).
Letx∗,y∗ be a fixed point of the map F=f,g and JF=JF(x∗,y∗) be the Jacobian matrix evaluated at the fixed point (x∗,y∗). Then, ρJF<1 if and only if(3)tr⁡JFx∗,y∗-1<det⁡JFx∗,y∗<1,where ρ(JF) is the spectral radius of the matrix JF.Remark 7.
Notice that Theorem6 states that the fixed point (x∗,y∗) is locally asymptotically stable if conditions (3) are verified.The case of nonhyperbolic fixed point is more involved. There are several scenarios according to the eigenvalues which are1 or -1. When one of the eigenvalues lies in the unit circle and the other eigenvalue is inside the unit circle, it is common to use centre manifold theory [6, 8, 20] in order to determine the stability of the fixed point, since the dynamics on the centre manifold determines the dynamics of the original system.In Section4 we will give the principal definitions and techniques behind the study of nonhyperbolic fixed points for a general dimension. In this specific case of dimension two, a centre manifold is a graph of type y=h(x) (see Theorem 11) on which the dynamics of (2) is given locally by the map on the centre manifold (see Theorem 12).Since the map on the centre manifold, in this case, is a one-dimensional map, one can then use the results presented in the previous subsection and conclude the stability of the nonhyperbolic fixed point. At this point we refer the pedagogical examples present in Elaydi’s book [9, pages 244–248] and in Zhang’s book [11, pages 286–288].


## 2.1. One-Dimensional Maps

Consider an intervalI⊆R and a map f:I→I. A point x∗∈R is said to be a fixed point (or equilibrium point) of f if f(x∗)=x∗, and given x0∈R, we define its orbit O(x0) as the set of points O(x0)={x0,f(x0),f2(x0),f3(x0),…}, where fn=f∘fn-1, for n=1,2,….One of the main objectives of the theory of discrete dynamical systems and, in particular, of the stability theory is the study of the behavior of orbits near fixed points, that is, the behavior of solutions of difference equations when the starting points are near equilibrium points. Hence, a basic definition in this field is needed. LetN denote the set of nonnegative integers.Definition 1 (local stability).
Letf:I→I be a map and x∗ be a fixed point of f, where I is an interval of real numbers. Then(1)
the fixed pointx∗ is said to be locally stable if, for any ϵ>0, there exits δ>0 such that, for all x0∈I with |x0-x∗|<δ, we have |fn(x0)-x∗|<ϵ, for all n∈N. Otherwise, the fixed point x∗ will be called unstable;(2)
the fixed pointx∗ is said to be attracting if there exists η>0 such that |x0-x∗|<η implies limn→∞⁡fn(x0)=x∗;(3)
the fixed pointx∗ is said to be locally asymptotically stable if it is both stable and attracting. If in the previous item η=∞, then x∗ is said to be globally asymptotically stable.One of the most effective graphical iteration methods to depict stability of fixed points is the cobweb diagram (also known as stair-step diagram) since it is possible to observe the behavior of the orbit of an initial point near the fixed point. It consists in a sequence of line segments connecting the points in the orbit of an initial pointx0, that is, a sequence of line segments with end points (x0,x1) and (x1,x1), (x1,x1) and (x1,x2), (x1,x2) and (x2,x2), (x2,x2) and (x2,x3), and so forth. On the cobweb plot, a locally stable fixed point corresponds to an inward spiral, while an unstable fixed point is an outward one. It follows from the definition of a fixed point that these spirals will center at a point where the diagonal y=x line crosses the function graph. This technique may be found in many books on discrete dynamical systems or difference equations as are the books by Alligood et al. [18], by Devaney [19], by Elaydi [9, 20], by Zhang [11], and by others, as well. In Figure 1 we present an illustration of the cobweb diagram for the Ricker equation xn+1=xner-xn, for n=0,1,2,… and r>0.Figure 1
Cobweb animation for the Ricker model given byxn+1=xner-xn. Here we are using r=1.8 and one can observe two orbits, the orbit of an initial point x0=0.2 and the orbit of an initial point x0=2.3. Both orbits converge to the fixed point x∗=1.8 meaning that the fixed point x∗=1.8 is locally stable.However, cobweb diagrams are not the most efficient tool to study local stability. There exists a simple but powerful criterion for knowing the local stability of fixed points. We may divide the fixed points into two categories: hyperbolic and nonhyperbolic. A fixed pointx∗ of a map f is said to be hyperbolic if |f′(x∗)|≠1. Otherwise, it is nonhyperbolic.Theorem 2 (see Elaydi [9], page 25).
Letx∗ be a hyperbolic fixed point of a map f, where f is continuous and differentiable at x∗. The following statements hold true:(1)
If|f′(x∗)|<1, then x∗ is locally asymptotically stable.(2)
Iff′x∗>1, then x∗ is unstable.The stability criteria for nonhyperbolic fixed points are more complex and are summarized in the following theorem (for a complete classification of nonhyperbolic fixed points cf. [9], page 33). Before presenting the criteria we introduce the notion of Schwarzian derivative.Definition 3 (Schwarzian derivative).
The Schwarzian derivative,Sf, of a function f, is defined by (1)Sfx=f′′′xf′x-32f′′xf′x2.In particular, when f′(x∗)=-1, we have Sf(x∗)=-f′′′x∗-(3/2)f′′x∗2.Theorem 4 (see Elaydi [9], pages 28–30).
Letx∗ be a fixed point of a map f and f′, f′′, and f′′′ be continuous at x∗.(1)
Letf′(x∗)=1.
(a)
Iff′′(x∗)>0, then x∗ is unstable but semistable from the left.(b)
Iff′′(x∗)<0, then x∗ is unstable but semistable from the right.(c)
Iff′′(x∗)=0 and f′′′(x∗)>0, then x∗ is unstable.(d)
Iff′′(x∗)=0 and f′′′(x∗)<0, then x∗ is locally asymptotically stable.(2)
Letf′(x∗)=-1.
(a)
IfSf(x∗)<0, then x∗ is locally asymptotically stable.(b)
IfSf(x∗)>0, then x∗ is unstable.


## 2.2. Two-Dimensional Maps

Definition1 (local stability) can now be extended to two-dimensional models (or higher dimensional models as well), using an appropriate norm. Since it will follow the same guidelines, we will omit it here.The local stability of a fixed point(x∗,y∗) of the equation(2)xn+1,yn+1=Fxn,yn,F=f,gis guaranteed when all the eigenvalues of the Jacobian, JF(x∗,y∗), evaluated at the fixed point (x∗,y∗) have magnitude less than one. This is precisely the conclusion of the following result; compare Theorem4.11 in Elaydi [9], page 221.Theorem 5 (see Elaydi [9], page 221).
LetF:D⊂R2→R2 be a C1 map, where D is an open subset of R2, (x∗,y∗) a fixed point of F, and JF=JF(x∗,y∗) the Jacobian matrix of F evaluated at the fixed point (x∗,y∗) and let ρ(JF) be the spectral radius of the matrix JF. Then the following statements hold true: (1)
Ifρ(JF)<1, then (x∗,y∗) is locally asymptotically stable.(2)
Ifρ(JF)>1, then (x∗,y∗) is unstable.(3)
Ifρ(JF)=1, then (x∗,y∗) may or may not be stable.The necessary and sufficient conditions can be determined by using the trace,tr, and the determinant, det, of JF as is summarized in the following result; compare Theorem4.4 in Elaydi [9], page 200:Theorem 6 (see Elaydi [9], page 200).
Letx∗,y∗ be a fixed point of the map F=f,g and JF=JF(x∗,y∗) be the Jacobian matrix evaluated at the fixed point (x∗,y∗). Then, ρJF<1 if and only if(3)tr⁡JFx∗,y∗-1<det⁡JFx∗,y∗<1,where ρ(JF) is the spectral radius of the matrix JF.Remark 7.
Notice that Theorem6 states that the fixed point (x∗,y∗) is locally asymptotically stable if conditions (3) are verified.The case of nonhyperbolic fixed point is more involved. There are several scenarios according to the eigenvalues which are1 or -1. When one of the eigenvalues lies in the unit circle and the other eigenvalue is inside the unit circle, it is common to use centre manifold theory [6, 8, 20] in order to determine the stability of the fixed point, since the dynamics on the centre manifold determines the dynamics of the original system.In Section4 we will give the principal definitions and techniques behind the study of nonhyperbolic fixed points for a general dimension. In this specific case of dimension two, a centre manifold is a graph of type y=h(x) (see Theorem 11) on which the dynamics of (2) is given locally by the map on the centre manifold (see Theorem 12).Since the map on the centre manifold, in this case, is a one-dimensional map, one can then use the results presented in the previous subsection and conclude the stability of the nonhyperbolic fixed point. At this point we refer the pedagogical examples present in Elaydi’s book [9, pages 244–248] and in Zhang’s book [11, pages 286–288].


## 3. Local Stability Conditions in 3D Systems

An autonomous three-dimensional difference equation is usually given by(4)xn+1=fxn,yn,znyn+1=gyn,yn,znzn+1=hxn,yn,zn,where n=0,1,2,… and the functions f, g, and h are of class Ck, with k being as big as we want. For simplicity, let us write system (4) as xn+1=F(xn), where x=(x,y,z) and F=(f,g,h).As mentioned before, we will consider an extension of the stability Definition1 taking a convenient norm in R3. Since one of the main objectives in the field of discrete dynamical systems is the study of the dynamics near the fixed points of F, that is, the local stability of the fixed points, a suitable generalization of Theorem 5 states that the fixed point x∗ of F is locally asymptotically stable if all the eigenvalues of the Jacobian, JF(x∗), evaluated at the fixed point x∗ have magnitude less than one. These conditions can be expressed in terms of the determinant of the Jacobian, det⁡(JF(x∗)), the trace of the Jacobian, tr⁡(JF(x∗)), and the sum of principal minors of the Jacobian, ∑i=13Mii(JF(x∗)). Moreover, these three values can be expressed as functions of the eigenvalues of JF(x∗), λi, for i=1,2,3 (cf. Section2 of Brooks [21]). Namely,(5)det⁡JFx∗=λ1λ2λ3,tr⁡JFx∗=λ1+λ2+λ3,∑i=13MiiJFx∗=λ1λ2+λ1λ3+λ2λ3.The characteristic equation ofJF(x∗) in this case can be given by(6)λ3-tr⁡JFx∗λ2+∑i=13MiiJFx∗λ-det⁡JFx∗=0,which has the three eigenvalues as roots. In order to determine conditions that guarantee that the characteristic roots are inside the unit disk, one can use the Schur-Cohn Criterion; compare Theorem5.1 in Elaydi [20, page 248].Theorem 8 (see Elaydi, [20] page 248).
The zeros of the characteristic polynomial(7)pλ=λk+p1λk-1+⋯+pklie inside the unit disk if and only if the following hold: (i)
p(1)>0.(ii)
(-1)kp(-1)>0.(iii)
The(k-1)×(k-1) matrices(8)Bk-1±=10⋯00p11⋯00⋮⋱⋮pk-3pk-4⋯10pk-2pk-3⋯p11±00⋯0pk00⋯pkpk-1⋮⋱⋮0pk⋯p4p3pkpk-1⋯p3p2are positive innerwise (the determinants of all of its inners are positives), fork=2,3,….Hence, using the Schur-Cohn Criterion, it follows that the necessary and sufficient conditions for which all the roots of (6) lie inside the unit disk (cf. relation (5.1.17) in Elaydi [9], page 248) are(9)tr⁡JFx∗+det⁡JFx∗<1+∑i=13MiiJFx∗∑i=13MiiJFx∗-tr⁡JFx∗det⁡JFx∗<1-det2JFx∗.Thus, conditions (9) imply that the fixed point x∗ of (4) is locally asymptotically stable. A remark about conditions (9) is in order.Remark 9.
Since in most of the concrete cases we deal with linear stability of three-dimensional discrete models depending on parameters, the elements of the Jacobian are most of the times expressed in terms of different parameters. Hence, direct calculations of the eigenvalues through the classical methods are not practical. The main feature of conditions (9) relies in the fact that one can obtain stability conditions, depending on the parameters, without computing the eigenvalues.Before ending this section, let us look into the region of local stability given by conditions (9). This region, where all the eigenvalues are less than one in magnitude, represents the region of local stability of a fixed point in a three-dimensional discrete dynamical system given by (4). Plotting the region defined by these four conditions in a coordinate system representing the trace (Tr), the determinant (Det), and the sum of the minors (SMin), the result is the picture depicted in Figure 2.Figure 2
Region of local stability of a fixed point in the triple(tr,det,∑Mii). Each box represents the same region with different points of view.


## 4. Bifurcation: Stability Using the Centre Manifold Theory

If some of the eigenvalues ofJF(x∗) have absolute value equal to one and the other eigenvalues have absolute value less than one, at least one of the conditions of (9) will be replaced by an equality. When this happens we say that a bifurcation occurs. In this case a small smooth change in the parameter values (the bifurcation parameters) of a system causes a sudden qualitative or topological change in its behavior. For a general framework in this area we refer to Kuznetzov’s book [8], where the author studies in depth the bifurcation phenomenon in both continuous and discrete dynamical systems.In certain cases it is possible to study the stability of the nonhyperbolic fixed point via the centre manifold theory since the dynamics on the centre manifold determines the dynamics of the original system.In the sequel, we will summarize the principal ideas behind the study of stability via the centre manifold theory in discrete dynamical systems. For general framework in the field of difference equations or discrete dynamical systems we refer to the books by Carr [6], Elaydi [9], and Kuznetzov [8] and for specific results and terminology the papers by Guzowska et al. [13], Luis et al. [16], and Psarros et al. [17].We will split our analysis into two cases according to the Jacobian matrix at the fixed point which is a block diagonal matrix or an upper (lower) triangular matrix.


### 4.1. Invariant Manifolds

LetF:Rk→Rk be a map such that F∈C2 and F(0)=0. Then one may write F as a perturbation of a linear map L,(10)Fx=Lx+Rx,where L is a k×k matrix defined by L=JF(0), R(0)=0, and JR(0)=0, and J denotes the Jacobian matrix. Now we will introduce special subspaces of Rk, called invariant manifolds (cf. Wiggis [10], page 28), that will play a central role in our study of stability and bifurcation.An invariant manifold is a manifold embedded in its phase space that it is invariant under the dynamical system generated byF. A subspace M of Rk is an invariant manifold if, whenever x∈M, Fn(x)∈M, for all n∈Z+. For the linear map L, one may split its spectrum σ(L) into three sets σs, σu, and σc, such that λ∈σs if |λ|<1, λ∈σu if |λ|>1, and λ∈σc if |λ|=1. Associated with these sets, there exist three invariant manifolds (linear subspaces) Es, Eu, and Ec which are the generalized eigenspaces corresponding to σs, σu, and σc, respectively. It should be noted that some of these subspaces may be trivial subspaces.The main question here is how to extend this linear theory to nonlinear maps. Corresponding to each of the linear subspacesEs, Eu, and Ec is an invariant manifold: the stable manifold Ws, the unstable manifold Wu, and the centre manifold Wc.The centre manifold theory is interesting only ifWu={0}, for, in this case, the dynamics on the centre manifold Wc determines the dynamics of the system.LetEs⊂Rs, Eu⊂Ru, and Ec⊂Rt, with s+u+t=k. Then one may formally define the above-mentioned invariant manifolds as follows: (11)Ws=x∈Rk∣Fnx⟶0asn⟶∞,Wu=x∈Rk∣Fnx⟶0asn⟶-∞.Since the stability on the centre manifoldWc is not a priori known, we will define it as a manifold of dimension t whose graph is tangent to Ec at the origin. It is noteworthy to mention that the centre manifold isnot unique, while the stable and unstable manifolds areunique.The next result summarizes the basic invariant manifolds theory.Theorem 10 (see Marsden and McCracken [22] page 28).
Suppose thatF∈C2. Then there exist C2 stable Ws and unstable Wu manifolds tangent to Es and Eu, respectively, at X=0 and C1 centre manifold Wc tangent to Ec at X=0. Moreover, Wc, Ws, and Wu are all invariant manifolds.


### 4.2. Centre Manifolds: Block Diagonal Matrix

By suitable change of variables, one may represent the mapF from the k-dimensional equation xn+1=F(xn) as the following system of difference equations:(12)xn+1=Axn+fxn,ynyn+1=Byn+gxn,yn.First we assume that the absolute value of all eigenvalues of thet×t matrix A is one, and all the eigenvalues of the s×s matrix B are less than one in magnitude, where t+s=k. Moreover, (13)f0,0=0,g0,0=0,Jf0,0=0,Jg0,0=0.The following theorem guarantees the existence of a centre manifold, which is a curve y=h(x), on which the dynamics of system (12) is given by the map on the centre manifold. The result may be found in Carr’s book [6], in Elaydi’s book [9, page 243] (cf. Theorem5.1), or in Zhang’s book [11, page 284] (cf. Theorem6.7.1).Theorem 11.
There exists a centre manifold for system (12) tangent to x,y∈Rt×Rs∣y=0 that can be represented locally as the graph of a function h:Rt→Rt such that (14)Wc=x,y∈Rt×Rs∣y=hx,h0=0,Jh0=0,x<δfor  a  sufficiently  smallδ.Furthermore, the dynamics restricted to Wc is given locally by the equation(15)xn+1=Axn+fxn,hxn,x∈Rt.The main feature of (15) is that its dynamics determine the dynamics of (12). So if x∗=0 is a stable, asymptotically stable, or unstable fixed point of (15), then the fixed point (x∗,y∗)=(0,0) of (12) possesses the corresponding property. This is precisely the conclusion of the following result (cf. Theorem5.2 in Elaydi [9, page 243] or cf. Theorem6.7.2 in Zhang [11, page 285]).Theorem 12.
Letx∗=0→t be a t-dimensional zero fixed point of (15). If x∗=0→t is a stable, asymptotically stable, or unstable fixed point of (15), then the fixed point (x∗,y∗) of (12) is stable, asymptotically stable, or unstable fixed point, respectively, where y∗=0→s is a s-dimensional zero vector.To find the mapy=h(x), we substitute for yn for h(xn) in (12) and obtain(16)xn+1=Axn+fxn,hxnyn+1=hxn+1=hAxn+fxn,hxn.But(17)yn+1=Byn+gxn,yn=Bhxn+gxn,hxn.Equating (16) and (17) yields the centre manifold equation(18)hAxn+fxn,hxn=Bhxn+gxn,hxn.Analogously, if all eigenvalues ofA are less than one in magnitude and the absolute value of all the eigenvalues of B is one, one may define the centre manifold Wc and obtain the equation (19)yn+1=Byn+ghyn,yn,where x=h(y).We remark that, in most cases, it is not possible to find an explicit expression forH. Since the main interest is to study locally the properties of the system, we may use some approximation. The most common and useful approximation is the power series as it is established in the following result that may be found in Carr’s book [6], in Zhang’s book [11, page 286], or in a most recent paper of Psarros et al. [17] (cf. Theorem1.3).Theorem 13.
Letψ:Rt→Rs be a C1 map with ψ(0→t)=0→s and Jψ(0→t)=0→s×t, where Jψ is the Jacobian matrix of ψ, and 0→t, 0→s, and 0→s×t are the 1×t, 1×s, and s×t zero matrices, respectively. Suppose that F(ψ(x))=xr as x→0, for some r>1. Then(20)hx=ψx+xr.


### 4.3. Centre Manifolds: Upper (Lower) Triangular Matrix

In working with concrete maps, it is beneficial in certain cases to deal with system where the derivative evaluated at the fixed point is not a block diagonal matrix. Let us now consider, without loss of generality, that the system of difference equationsXn+1=F(Xn) is written of the form(21)xn+1yn+1=AC0Bxnyn+fxn,yngxn,yn.There are two cases to consider: (i) the absolute value of all the eigenvalues ofA is one and all the eigenvalues of B are less than one in absolute value and (ii) the absolute value of all the eigenvalues of B is one and all the eigenvalues of A are less than one in absolute value. We will present details only in the first case, since the second follows in a similar way and can be found in Luis et al. [16].From the general theory in linear algebra, there exists a nonsingular matrixP of the form(22)P=P1P30P2such that(23)AB0C=PA00BP-1.By suitable changes of variables,x=P1u+P3v and y=P2v, one may write(24)un+1vn+1=A00Bunvn+P-1fP1u+P3v,P2vgP1u+P3v,P2v.One can apply the centre manifold theorem and obtain a mapv=h(u) with h(0)=0=h′(0). The dynamics of (24) is completely determined by the dynamics of the equation(25)un+1=Aun+P~1fP1un+P3hu,P2hu+P~3gP1un+P3hu,P2hu,where P~1, P~2, and P~3 are entries of the matrix(26)P-1=P~1P~30P~2.Hence the curvev=h(u) is the implicit solution of (25) and is the equation of the centre manifold. To find the map h we use the centre manifold equation(27)hAu+Chu+fu,hu=Bhu+gu,hu,and using Taylor approximation, we are able to determine locally the map h, as we will illustrate in the next section with the Ricker competition model.


## 4.1. Invariant Manifolds

LetF:Rk→Rk be a map such that F∈C2 and F(0)=0. Then one may write F as a perturbation of a linear map L,(10)Fx=Lx+Rx,where L is a k×k matrix defined by L=JF(0), R(0)=0, and JR(0)=0, and J denotes the Jacobian matrix. Now we will introduce special subspaces of Rk, called invariant manifolds (cf. Wiggis [10], page 28), that will play a central role in our study of stability and bifurcation.An invariant manifold is a manifold embedded in its phase space that it is invariant under the dynamical system generated byF. A subspace M of Rk is an invariant manifold if, whenever x∈M, Fn(x)∈M, for all n∈Z+. For the linear map L, one may split its spectrum σ(L) into three sets σs, σu, and σc, such that λ∈σs if |λ|<1, λ∈σu if |λ|>1, and λ∈σc if |λ|=1. Associated with these sets, there exist three invariant manifolds (linear subspaces) Es, Eu, and Ec which are the generalized eigenspaces corresponding to σs, σu, and σc, respectively. It should be noted that some of these subspaces may be trivial subspaces.The main question here is how to extend this linear theory to nonlinear maps. Corresponding to each of the linear subspacesEs, Eu, and Ec is an invariant manifold: the stable manifold Ws, the unstable manifold Wu, and the centre manifold Wc.The centre manifold theory is interesting only ifWu={0}, for, in this case, the dynamics on the centre manifold Wc determines the dynamics of the system.LetEs⊂Rs, Eu⊂Ru, and Ec⊂Rt, with s+u+t=k. Then one may formally define the above-mentioned invariant manifolds as follows: (11)Ws=x∈Rk∣Fnx⟶0asn⟶∞,Wu=x∈Rk∣Fnx⟶0asn⟶-∞.Since the stability on the centre manifoldWc is not a priori known, we will define it as a manifold of dimension t whose graph is tangent to Ec at the origin. It is noteworthy to mention that the centre manifold isnot unique, while the stable and unstable manifolds areunique.The next result summarizes the basic invariant manifolds theory.Theorem 10 (see Marsden and McCracken [22] page 28).
Suppose thatF∈C2. Then there exist C2 stable Ws and unstable Wu manifolds tangent to Es and Eu, respectively, at X=0 and C1 centre manifold Wc tangent to Ec at X=0. Moreover, Wc, Ws, and Wu are all invariant manifolds.


## 4.2. Centre Manifolds: Block Diagonal Matrix

By suitable change of variables, one may represent the mapF from the k-dimensional equation xn+1=F(xn) as the following system of difference equations:(12)xn+1=Axn+fxn,ynyn+1=Byn+gxn,yn.First we assume that the absolute value of all eigenvalues of thet×t matrix A is one, and all the eigenvalues of the s×s matrix B are less than one in magnitude, where t+s=k. Moreover, (13)f0,0=0,g0,0=0,Jf0,0=0,Jg0,0=0.The following theorem guarantees the existence of a centre manifold, which is a curve y=h(x), on which the dynamics of system (12) is given by the map on the centre manifold. The result may be found in Carr’s book [6], in Elaydi’s book [9, page 243] (cf. Theorem5.1), or in Zhang’s book [11, page 284] (cf. Theorem6.7.1).Theorem 11.
There exists a centre manifold for system (12) tangent to x,y∈Rt×Rs∣y=0 that can be represented locally as the graph of a function h:Rt→Rt such that (14)Wc=x,y∈Rt×Rs∣y=hx,h0=0,Jh0=0,x<δfor  a  sufficiently  smallδ.Furthermore, the dynamics restricted to Wc is given locally by the equation(15)xn+1=Axn+fxn,hxn,x∈Rt.The main feature of (15) is that its dynamics determine the dynamics of (12). So if x∗=0 is a stable, asymptotically stable, or unstable fixed point of (15), then the fixed point (x∗,y∗)=(0,0) of (12) possesses the corresponding property. This is precisely the conclusion of the following result (cf. Theorem5.2 in Elaydi [9, page 243] or cf. Theorem6.7.2 in Zhang [11, page 285]).Theorem 12.
Letx∗=0→t be a t-dimensional zero fixed point of (15). If x∗=0→t is a stable, asymptotically stable, or unstable fixed point of (15), then the fixed point (x∗,y∗) of (12) is stable, asymptotically stable, or unstable fixed point, respectively, where y∗=0→s is a s-dimensional zero vector.To find the mapy=h(x), we substitute for yn for h(xn) in (12) and obtain(16)xn+1=Axn+fxn,hxnyn+1=hxn+1=hAxn+fxn,hxn.But(17)yn+1=Byn+gxn,yn=Bhxn+gxn,hxn.Equating (16) and (17) yields the centre manifold equation(18)hAxn+fxn,hxn=Bhxn+gxn,hxn.Analogously, if all eigenvalues ofA are less than one in magnitude and the absolute value of all the eigenvalues of B is one, one may define the centre manifold Wc and obtain the equation (19)yn+1=Byn+ghyn,yn,where x=h(y).We remark that, in most cases, it is not possible to find an explicit expression forH. Since the main interest is to study locally the properties of the system, we may use some approximation. The most common and useful approximation is the power series as it is established in the following result that may be found in Carr’s book [6], in Zhang’s book [11, page 286], or in a most recent paper of Psarros et al. [17] (cf. Theorem1.3).Theorem 13.
Letψ:Rt→Rs be a C1 map with ψ(0→t)=0→s and Jψ(0→t)=0→s×t, where Jψ is the Jacobian matrix of ψ, and 0→t, 0→s, and 0→s×t are the 1×t, 1×s, and s×t zero matrices, respectively. Suppose that F(ψ(x))=xr as x→0, for some r>1. Then(20)hx=ψx+xr.


## 4.3. Centre Manifolds: Upper (Lower) Triangular Matrix

In working with concrete maps, it is beneficial in certain cases to deal with system where the derivative evaluated at the fixed point is not a block diagonal matrix. Let us now consider, without loss of generality, that the system of difference equationsXn+1=F(Xn) is written of the form(21)xn+1yn+1=AC0Bxnyn+fxn,yngxn,yn.There are two cases to consider: (i) the absolute value of all the eigenvalues ofA is one and all the eigenvalues of B are less than one in absolute value and (ii) the absolute value of all the eigenvalues of B is one and all the eigenvalues of A are less than one in absolute value. We will present details only in the first case, since the second follows in a similar way and can be found in Luis et al. [16].From the general theory in linear algebra, there exists a nonsingular matrixP of the form(22)P=P1P30P2such that(23)AB0C=PA00BP-1.By suitable changes of variables,x=P1u+P3v and y=P2v, one may write(24)un+1vn+1=A00Bunvn+P-1fP1u+P3v,P2vgP1u+P3v,P2v.One can apply the centre manifold theorem and obtain a mapv=h(u) with h(0)=0=h′(0). The dynamics of (24) is completely determined by the dynamics of the equation(25)un+1=Aun+P~1fP1un+P3hu,P2hu+P~3gP1un+P3hu,P2hu,where P~1, P~2, and P~3 are entries of the matrix(26)P-1=P~1P~30P~2.Hence the curvev=h(u) is the implicit solution of (25) and is the equation of the centre manifold. To find the map h we use the centre manifold equation(27)hAu+Chu+fu,hu=Bhu+gu,hu,and using Taylor approximation, we are able to determine locally the map h, as we will illustrate in the next section with the Ricker competition model.


## 5. Illustrative Example: Ricker Competition Model

The three-dimensional Ricker competition model [23, page 31] is usually written as(28)x1n+1=x1nexp⁡r1-∑i=13c1ixinx2n+1=x2nexp⁡r2-∑i=13c2ixinx3n+1=x3nexp⁡r3-∑i=13c3ixin,where xi(n), for i=1,2,3, represents the population sizes of species xi at time unit n, the parameters ri, for i=1,2,3, are the inherent exponential growth rates at low densities, and ci,j, for i,j=1,2,3, are the competition intensity coefficients measuring the effects of intraspecific competition and interspecific competition. More precisely, cii are the intraspecific competition parameters while cij, for i≠j, are the interspecific competition parameters. Notice that these twelve parameters are assumed to be positive.For the sake of simplicity, throughout this section we assume that the intraspecific competition parameters are equal to one; that is,cii=1, and the interspecific competition parameters lie on the unit interval; that is, 0<cij<1, for i≠j.System (28) can be represented by the difference equation xn+1=F(xn), where the map F=(f1,f2,f3) is given by(29)Fx1,x2,x3=x1er1-x1-c12x2-c13x3,x2er2-x2-c21x1-c23x3,x3er3-x3-c31x1-c32x2.The mapF has eight fixed points, namely, the origin O=(0,0,0), three fixed points on the axes E1=(r1,0,0), E2=(0,r2,0), and E3=(0,0,r3), three fixes points in each of the planes E12=(r1-c12r2)/(1-c12c21),(r2-c21r1)/(1-c12c21),0, E13=(r1-c13r3)/(1-c13c31),0,(r3-c31r1)/(1-c13c31), and E23=0,(r2-c23r3)/(1-c23c32),(r3-c32r2)/(1-c23c32), and a possible interior fixed point given by (30)E∗=c23c32-1r1+c13r3-c32r2+c12r2-c23r3c12c21-c23c31+c23c32+c13c31-c21c32-1,c13c31-1r2+c23r3-c31r1+c21r1-c13r3c12c21-c23c31+c23c32+c13c31-c21c32-1,c32r2+c31r1-c12r2+c21c12r3-c32r1-r3c12c21-c23c31+c23c32+c13c31-c21c32-1.Since we are mainly interested in population dynamics, we always assume that these fixed points are nonnegative.First, let us have a look in the dynamics in lower dimensions. In the absence of two species, sayx2 and x3, the dynamics of system (28) is the same as the one given by the one-dimensional equation xn+1=xner1-xn=f(xn). This equation has two fixed points, the origin, which is always unstable, and a positive fixed point x∗=r1, which is globally asymptotically stable when 0<r1≤2. At r1=2, a period-doubling bifurcation occurs since f′(2)=-1. This means that x∗=r1 becomes an unstable fixed point and a new locally asymptotically stable 2-periodic cycle is born. At r1≈2.512, a new period-doubling bifurcation occurs. The 2-periodic cycle becomes unstable and a new locally asymptotically stable 4-periodic cycle is born. Moreover, this scenario of period-doubling bifurcation route to chaos will continue as r1 increases.Now, in the absence of one species, say speciesx3, the dynamics of system (28) is the same as the one of the two-dimensional system(31)x1n+1=x1nexpr1-x1n-c12x2nx2n+1=x2nexpr2-c21x1n-x2n,forn=0,1,2,….The local dynamics and principal properties of the four fixed points of (31) have been studied in depth by Luis et al. [16] and is shortly as follows: (1)
The origin is always an unstable fixed point.(2)
The point(r1,0) is an exclusion fixed point in the x-axis which is locally asymptotically stable if and only if 0<r1≤2 and r2<c21r1.(3)
The point(0,r2) is an exclusion fixed point in the y-axis which is locally asymptotically stable if and only if 0<r2≤2 and r1<c12r2.(4)
The point(x1∗,x2∗)=(r1-c12r2)/(1-c12c21),(r2-c21r1)/(1-c12c21) is a positive fixed point in the interior of the first quadrant, which is locally asymptotically stable, if and only if(32)4c12c21-1+21-c12r2+21-c21r1≤c12r2-r1c21r1-r2<1-c12r2+1-c21r1.The stability regions, in the parameter spaceOr1r2, of the nontrivial fixed points, are denoted by R1, Q1, and S1 and are depicted in Figure 3. The region R1 is the stability region of (r1,0), the region Q1 is the stability region of (0,r2), and the region S1 is the stability region of (x1∗,x2∗).More precisely, the regionS1 is the region in the r1Or2 plane bounded by the lines r2=r1/c12 and r2=c21r1 and the curve γ1, where γ1 is part of a branch of the hyperbola defined by(33)c21r12+21-c21r1-1+c12c21r1r2+21-c12r2+c12r22+4c12c21-1=0.Notice that the positive equilibrium point exists if and only if(34)c12r2<r1,c21r1<r2or(35)c12r2>r1,c21r1>r2.It follows that (34) implies that c12c21<1 while (35) implies c12c21>1. When c12c21>1, the coexistence equilibrium point is a saddle and the asymptotic attractor of an orbit of (31) depends on its initial conditions (see [24] for specific details of this phenomenon in population dynamics). If c12c21=1, system (31) has no coexistence equilibrium point. Henceforth, we shall assume that c12c21<1 and it means that the interspecific competition is less than the intraspecific competition since c12c21<c11c22=1.Equivalently, in the region defined by (32) and under assumption (34), the coexistence fixed point is locally asymptotically stable if and only if (r,s)∈int(S1)∪γ1, where int(S1) denotes the interior of the region S1.Figure 3
The stability regions and the bifurcation scenario of the two-dimensional Ricker competition model (31), in the parameter space Or1r2, when the competition parameters c12 and c21 are fixed such that c12c21<1. The region S1 is the stability region of the coexistence equilibrium point (x1∗,x2∗), the region R1 is the stability region of the exclusion fixed point (r1,0), and the region Q1 is the stability region of the exclusion fixed point (0,r2). A period-doubling bifurcation scenario occurs (in the coexistence case) as we cross from region S1 to region S2, from region S2 to region S3, and so on. Similarly, one has a period-doubling scenario (in the exclusion case) in the x-axis as we cross from region R1 to region R2, from region R2 to region R3, and so on. Similar phenomenon occurs in regions Qi, for i=1,2,….As mentioned before, in order to study the stability of the fixed points of system (28), we have to determine the Jacobian of F which is given by(36)JFx1,x2,x3=-eΔ1x1-1-eΔ1c12x1-eΔ1c13x1-eΔ2c21x2-eΔ2x2-1-eΔ2c23x2-eΔ3c31x3-eΔ3c32x3-eΔ3x3-1,where Δ1=r1-x1-c12x2-c13x3, Δ2=r2-c21x1-x2-c23x3, and Δ3=r3-c31x1-c32x2-x3.It is a straightforward computation to see that the origin is an unstable fixed point since the eigenvalues ofJF(O) are eri>1, for i=1,2,3.


### 5.1. Fixed Points on the Coordinate Axis

We will reduce our analysis to the fixed pointE1=(r1,0,0). The stability of the fixed points Ei, for i=2,3, will follow in a similar way.The Jacobian evaluated at the fixed pointE1 is given by(37)JFE1=1-r1-c12r1-c13r10er2-c21r1000er3-c31r1.Clearly, E1 is locally asymptotically stable if 0<r1<2 and rj<cj1r1, for j=2,3. This region, in the parameter space of the growth rates (triple (r1,r2,r3)), is depicted in Figure 4.Figure 4
(a) Stability region of the fixed pointE1, in the parameter space Or1r2r3 when cj1=0.5, for j=2,3. The projection of this region in the r1Or2 plane is the region R1 depicted in Figure 3. (b) The stability regions of the fixed points Ei, for i=1,2,3, together in the parameter space Or1r2r3 when cii=1 and cij=0.5, for i≠j and i,j=1,2,3.
(a)
(b)Now, ifr1=2 and rj<2cj1, j=2,3 undergoes a period-doubling bifurcation on the x-axis. Under this scenario λ1=-1 is an eigenvalue of JF(E1) and for the other eigenvalues we have |erj-2cj1|<1, for j=2,3. Hence, we are in position to study the stability of the fixed point E1 via the centre manifold theory.Shifting the fixed pointE1 to the origin, it is a straightforward computation to show that system (28) is equivalent to(38)y1n+1=-y1n-2c12y2n-2c13y3n+f~1y1n,y2n,y3ny2n+1=er2-2c21y2n+f~2y1n,y2n,y3ny3n+1=er3-2c31y3n+f~3y1n,y2n,y3n,where (39)f~1y1,y2,y3=y1+2e-c12y2-c13y3-y1+2c12y2+2c13y3+y1-2,f~2y1,y2,y3=er2y2e-c21y1+2-c23y3-y2-e-2c21,f~3y1,y2,y3=er3y3e-c31y1+2-c32y2-y3-e-2c31.The centre manifold in this case can be represented locally as the graph of a function Φ:R2→R2 such that (40)Wc=y1,y2,y3∈R×R2∣y2,y3=Φ1y1,Φ2y1,Φ0=0,JΦ0=0,y<δfor  a  sufficiently  smallδ.The centre manifold equation (27) is therefore given by(41)Φ1-y1-2c12Φ1y1-2c13Φ2y1+f~1y1,Φ1y1,Φ2y1=er2-2c21Φ1y1+f~2y1,Φ1y1,Φ2y1Φ2-y1-2c12Φ1y1-2c13Φ2y1+f~1y1,Φ1y1,Φ2y1=er3-2c31Φ2y1+f~3y1,Φ1y1,Φ2y1,where Φi, i=1,2, can be represented locally by Φi(u)=αiu2+βiu3+O[u]4.Using Mathematica or Maple, we are able to simplify system (41), and using power series one can conclude that (locally) system (41) is equivalent to(42)1-er2-2c21α1y12+c21er2-2c21α1+4c13α1α2+4c12α12-1+er2-2c21β1y13+Oy14=0.α1-α2er3-2c31y12+-β1+4c12α12+4c13α1α2+c31er3-2c31α2-er3-2c31β2y13+Oy14=0In order to have equality, from this system we have to conclude that (43)1-er2-2c21α1=0c21er2-2c21α1+4c13α1α2+4c12α12-1+er2-2c21β1=0α1-α2er3-2c31=0-β1+4c12α12+4c13α1α2+c31er3-2c31α2-er3-2c31β2=0,yielding the trivial solution α1=α2=β1=β2=0.A straightforward computation shows that locally, the map in (25), in this case, can be given by(44)fy1=e-y1y1+2-2.Since f′(0)=-1, we have to compute the Schwarzian derivative at the fixed point 0 (cf. Theorem 4); that is,(45)f′′′0f′0-32f′′0f′02=-1<0. Hence, the fixed point 0 is locally asymptotically stable.Consequently, the fixed pointE1 of system (28) is locally asymptotically stable if r1=2 and rj<2cj1, for j=2,3. These specific values of the parameters undergo a period-doubling bifurcation. As r1 increases, the exclusion fixed point E1 becomes unstable and an asymptotically stable 2-periodic cycle of the form(46)r¯10,0,0,r¯11,0,0,r¯10<r1<r¯11is born in the x-axis.From now on, computations are only possible via numerical simulations. One can observe that atr1≈2.512 a new period-doubling bifurcation occurs. Hence, the exclusion 2-periodic cycle on the x-axis becomes unstable and a new 4-periodic cycle in the x-axis is born. Further simulations suggest that this scenario of period-doubling bifurcation, in the x-axis, will continue route to chaos.We remark that this period-doubling scenario, in three-dimensional mappings, is just a direct observation from the simulations since there is no extension of Sharkovsky’s theorem [25, 26] for dimension k≥2.It remains to study the stability of the fixed pointE1 when r1=2, r2=2c21, and r3<2c31. Under this restriction on the parameters we have λ1=-1, λ2=1, and |λ3|<1. The reduced equation on the centre manifold will be a two-dimensional system. Unfortunately, this scenario in two-dimensional discrete systems (λ~1=-1,λ~2=1) remains as an open problem, and thus we are not able to discuss the stability of the fixed point here.


### 5.2. Fixed Points on the Planes

Let us now consider the fixed pointE12 on the plane x1Ox2. The stability analysis of E13 and E23 follows in a similar way.The Jacobian evaluated at the fixed point is given by(47)JFE12=1-c12c21+c12r2-r11-c12c21c12c12r2-r11-c12c21c13c12r2-r11-c12c21c21c21r1-r21-c12c211-c12c21+r1c21-r21-c12c21c23c21r1-r21-c12c2100ec21c32-c31r1+c31c12-c32r2+1-c21c12r3/1-c12c21.It is a straightforward computation to see that the conditions of local stability are given by(48)4c12c21-1+21-c12r2+21-c21r1<c12r2-r1c21r1-r2<1-c12r2+1-c21r1,(49)c21c32-c31r1+c31c12-c32r2+1-c21c12r3<0.In Figure 5 a prototype of this region is depicted in the triple (r1,r2,r3) when cji=0.5, for i≠j.Figure 5
(a) Stability region of the exclusion fixed pointE12 in the parameter space Or1r2r3 when cji=0.5 for i≠j. The projection of this region in the plane r1Or2 is the region S1 depicted in Figure 3. (b) Stability regions of E12, E13, and E23, together in the parameter space Or1r2r3 when cii=1 and cij=0.5, for i≠j and i,j=1,2,3.
(a)
(b)Considering the equality(50)4c12c21-1+21-c12r2+21-c21r1=c12r2-r1c21r1-r2and keeping the other relations, the eigenvalues of JF(E12) are such that λ1=-1 and |λi|<1, for i=2,3. Hence, we are in position to use the centre manifold theorem. Following the techniques employed in the previous subsection, one can conclude that E12 is a locally asymptotically stable fixed point on the plane x1Ox2 under these restrictions on the parameters. Notice that, in the parameter space, the projection of this region in the r1Or2 plane is part of a branch of an hyperbola as described in the beginning of this section and can be observed in Figure 3 by γ1.This means that this specific region undergoes a period-doubling bifurcation on thex1Ox2 plane. Hence, as the parameters increase, the fixed point E12 becomes unstable and a locally asymptotically stable 2-periodic cycle of the form(51)x¯0,y¯0,0,x¯1,y¯1,0is born on the x1Ox2 plane. As in the precedent case, from now on, computations are only possible via numerical simulations. These simulations suggest a scenario of period-doubling bifurcation on the x1Ox2 plane route to chaos.We should mention that if we replace the inequality signs by equality on (48) and condition (49) is satisfied, then λ1=-1, λ2=1, and |λ3|<1 are the eigenvalues of JF(E12). Unfortunately, we are not able to use here the centre manifold theorem, as mentioned at the end of the previous subsection.


### 5.3. Positive Fixed Point

We now consider the positive fixed pointE∗ known in population dynamics as the coexistence fixed point. Due to the great number of parameters involved and for the sake of simplicity, we will study the symmetric case, that is, when cij=c, for all i≠j. The fixed point E∗=(x1∗,x2∗,x3∗) becomes in this case(52)E∗=cr2+r3-c+1r1c-12c+1,cr1-c+1r2+cr3c-12c+1,cr1+cr2-c+1r3c-12c+1,and the derivative at E∗ is given by(53)JFE∗=-2c2+r2c+r3c+c-c+1r1+1-2c2+c+1cc+1r1-cr2+r3c-12c+1cc+1r1-cr2+r3c-12c+1c-cr1+c+1r2-cr3c-12c+1-2c2+r1c+r3c+c-c+1r2+1-2c2+c+1c-cr1+c+1r2-cr3c-12c+1c-cr1-cr2+c+1r3c-12c+1c-cr1-cr2+c+1r3c-12c+1-2c2+r1c+r2c-r3c+c-r3+1-2c2+c+1.For this matrix we have(54)tr⁡JFE∗=-6c+r1+r2+r3-32c+1,∑i=13MiiJFE∗=3P0+2P1r1+r2+r3+P2r12+r22+r32P0+P4r1r2+r1r3+r2r3P0,det⁡JFE∗=-P0+∑i=13Pir1i+r2i+r3i+P4r1r2+r1r3+r2r3P0-P5r1r22+r32+r2r12+r32+r3r12+r22+P6r1r2r3P0,where (55)P0=4c3-3c-1,P1=-2c2+c+1,P2=c3+3c2+2cP3=c3+c2,P4=-2c3+2c2+c+1,P5=-c3+c2+cP6=1+3c+6c2+2c3.The conditions of local stability ofE∗ now follow from relation (9). A prototype of this region, in the parameter space Or1r2r3, is represented in Figure 6. Notice that this region is a three-dimensional solid with planar faces, with the exception of the top, where the surface is given by(56)tr⁡JFE∗+det⁡JFE∗-∑i=13MiiJFE∗=1.Figure 6
Stability region of the positive fixed pointE∗ in the parameter space Or1r2r3 when cii=1 and cji=0.5, for i≠j and i,j=1,2,3.Now, in conditions (9), when tr⁡(JF(E∗))+det⁡(JF(E∗))-∑i=13Mii(JF(E∗))=1 and assuming the others inequalities, we have λ1=-1 and |λi|<1, for i=2,3. Hence, we have to use the centre manifold theory in order to determine the stability of E∗.Making the change of variables(57)y1=x1-x1∗,y2=x2-x2∗,y3=x3-x3∗, we shift the fixed point to the origin. System (28) is now equivalent to(58)Yn+1=JYn+GYn,where(59)Yn=y1ny2ny3n,J=J11J12J13J21J22J23J31J32J33,GYn=g1Yng2Yng3Yn,J=-2c2+cr2+r3+c+11-r1-2c2+c+1cc+1r1-cr2+r3c-12c+1cc+1r1-cr2+r3c-12c+1-ccr1-c+1r2+cr3c-12c+1-2c2+r1c+r3c+c-c+1r2+1-2c2+c+1-ccr1-c+1r2+cr3c-12c+1-ccr1+cr2-c+1r3c-12c+1-ccr1+cr2-c+1r3c-12c+1-2c2+r1c+r2c-r3c+c-r3+1-2c2+c+1,giy1,y2,y3=fiy1+x1∗,y2+x2∗,y3+x3∗-xi∗-∑j=13Jijyj,i=1,2,3.SinceJ is nonsingular, there exist a diagonal matrix Λ and the matrices P and P-1 such that J=PΛP-1. Making a new change of variables Y=PZ, where Z=(z1,z2,z3), one can write (58) as(60)Zn+1=ΛZn+HZn,where(61)Zn=z1nz2nz3n,Λ=-1000λ2000λ3,HZn=P-1GP×Zn.The centre manifold in this case can be represented locally as the graph of a functionΦ:R2→R2 such that (62)Wc=z1,z2,z3∈R×R2∣z2,z3=Φ1z1,Φ2z1,Φ0=0,DΦ0=0,z1<δfor  a  sufficiently  smallδ.Therefore, the centre manifold equation (27) is given by (63)Φ1-z1+h1z1,Φ1z1,Φ2z1=λ2Φ1z1+h2z1,Φ1z1,Φ2z1Φ2-z1+h1z1,Φ1z1,Φ2z1=λ3Φ2z1+h3z1,Φ1z1,Φ2z1,where Φi, i=1,2, can be represented locally by Φi(u)=αiu2+βiu3+O[u]4 and (64)hiz1,z2,z3=∑j=13P~ijgjz1,z2,z3,i=1,2,3,where P~ij are the entries of the matrix P-1.Unfortunately, at this stage the computations are long and we are not able to find explicitly the matricesP, Λ, and P-1 such that J=PΛP-1. However, we are able to find it numerically. Hence, we will illustrate the study with a concrete example. Let c=0.5 and r2=r3=1.5 (the choice of these values is made considering the region obtained in Figure 6). Substituting in(65)tr⁡JFE∗+det⁡JFE∗-∑i=13MiiJFE∗=1, we obtain r1≈2.1547. Under these values of the parameters, the matrices are given by(66)P=-0.97690-0.6676-0.1511-0.70710.5265-0.15110.70710.5265,Λ=-1.0000.78870000.634,P-1=-0.8558-0.5425-0.54250-0.70710.7071-0.24570.7940.794.And the functions hi(z1,z2,z3), for i=1,2,3, are (67)h1z1,z2,z3=1.9409+z1+0.836z1+0.5713z3-1.4823e1.128z1+0.1411z3+0.082z1-0.3836z2-0.2856z3-0.2293e0.7151z1-0.3536z2-0.4559z3+0.082z1+0.3836z2-0.2856z3-0.2293e0.7151z1+0.3536z2-0.4559z3,h2z1,z2,z3=-0.7887z2+0.1069z1+0.5z2-0.3723z3-0.2989e0.7151z1+0.3536z2-0.456z3+-0.1069z1+0.5z2+0.3723z3+0.2989e0.7151z1-0.3536z2-0.4559z3,h3z1,z2,z3=-0.634z3-0.2456+0.24z1+0.164z3-0.4255e1.128z1+0.1411z3+0.3356-0.12z1+0.5614z2+0.418z3e0.7151z1-0.3536z2-0.456z3+-0.12z1-0.5614z2+0.418z3+0.336e0.7151z1+0.354z2-0.456z3.Using Taylor series and solving system (63), one can conclude that the dynamics of the system is given locally by the dynamics of z1(n+1)=h~(z1(n)), where(68)h~z1=e0.7151z10.164z1-0.4586+e1.128z10.836z1-1.4823+1.9409.Notice that this corresponds to the equation on the centre manifold that is generically given by (15).Sinceh~′(0)=-1, we have to use Theorem 4. The Schwarzian derivative of h~ at 0 is ≈-1.1476<0. It follows from Theorem 4 that, in this particular case, the fixed point is locally asymptotically stable.We remark that, due to the difficulties encountered in finding the matrixP, one can write a script in Mathematica or Maple to automatically obtain, via numerical approximation, the map h~ in each case.


## 5.1. Fixed Points on the Coordinate Axis

We will reduce our analysis to the fixed pointE1=(r1,0,0). The stability of the fixed points Ei, for i=2,3, will follow in a similar way.The Jacobian evaluated at the fixed pointE1 is given by(37)JFE1=1-r1-c12r1-c13r10er2-c21r1000er3-c31r1.Clearly, E1 is locally asymptotically stable if 0<r1<2 and rj<cj1r1, for j=2,3. This region, in the parameter space of the growth rates (triple (r1,r2,r3)), is depicted in Figure 4.Figure 4
(a) Stability region of the fixed pointE1, in the parameter space Or1r2r3 when cj1=0.5, for j=2,3. The projection of this region in the r1Or2 plane is the region R1 depicted in Figure 3. (b) The stability regions of the fixed points Ei, for i=1,2,3, together in the parameter space Or1r2r3 when cii=1 and cij=0.5, for i≠j and i,j=1,2,3.
(a)
(b)Now, ifr1=2 and rj<2cj1, j=2,3 undergoes a period-doubling bifurcation on the x-axis. Under this scenario λ1=-1 is an eigenvalue of JF(E1) and for the other eigenvalues we have |erj-2cj1|<1, for j=2,3. Hence, we are in position to study the stability of the fixed point E1 via the centre manifold theory.Shifting the fixed pointE1 to the origin, it is a straightforward computation to show that system (28) is equivalent to(38)y1n+1=-y1n-2c12y2n-2c13y3n+f~1y1n,y2n,y3ny2n+1=er2-2c21y2n+f~2y1n,y2n,y3ny3n+1=er3-2c31y3n+f~3y1n,y2n,y3n,where (39)f~1y1,y2,y3=y1+2e-c12y2-c13y3-y1+2c12y2+2c13y3+y1-2,f~2y1,y2,y3=er2y2e-c21y1+2-c23y3-y2-e-2c21,f~3y1,y2,y3=er3y3e-c31y1+2-c32y2-y3-e-2c31.The centre manifold in this case can be represented locally as the graph of a function Φ:R2→R2 such that (40)Wc=y1,y2,y3∈R×R2∣y2,y3=Φ1y1,Φ2y1,Φ0=0,JΦ0=0,y<δfor  a  sufficiently  smallδ.The centre manifold equation (27) is therefore given by(41)Φ1-y1-2c12Φ1y1-2c13Φ2y1+f~1y1,Φ1y1,Φ2y1=er2-2c21Φ1y1+f~2y1,Φ1y1,Φ2y1Φ2-y1-2c12Φ1y1-2c13Φ2y1+f~1y1,Φ1y1,Φ2y1=er3-2c31Φ2y1+f~3y1,Φ1y1,Φ2y1,where Φi, i=1,2, can be represented locally by Φi(u)=αiu2+βiu3+O[u]4.Using Mathematica or Maple, we are able to simplify system (41), and using power series one can conclude that (locally) system (41) is equivalent to(42)1-er2-2c21α1y12+c21er2-2c21α1+4c13α1α2+4c12α12-1+er2-2c21β1y13+Oy14=0.α1-α2er3-2c31y12+-β1+4c12α12+4c13α1α2+c31er3-2c31α2-er3-2c31β2y13+Oy14=0In order to have equality, from this system we have to conclude that (43)1-er2-2c21α1=0c21er2-2c21α1+4c13α1α2+4c12α12-1+er2-2c21β1=0α1-α2er3-2c31=0-β1+4c12α12+4c13α1α2+c31er3-2c31α2-er3-2c31β2=0,yielding the trivial solution α1=α2=β1=β2=0.A straightforward computation shows that locally, the map in (25), in this case, can be given by(44)fy1=e-y1y1+2-2.Since f′(0)=-1, we have to compute the Schwarzian derivative at the fixed point 0 (cf. Theorem 4); that is,(45)f′′′0f′0-32f′′0f′02=-1<0. Hence, the fixed point 0 is locally asymptotically stable.Consequently, the fixed pointE1 of system (28) is locally asymptotically stable if r1=2 and rj<2cj1, for j=2,3. These specific values of the parameters undergo a period-doubling bifurcation. As r1 increases, the exclusion fixed point E1 becomes unstable and an asymptotically stable 2-periodic cycle of the form(46)r¯10,0,0,r¯11,0,0,r¯10<r1<r¯11is born in the x-axis.From now on, computations are only possible via numerical simulations. One can observe that atr1≈2.512 a new period-doubling bifurcation occurs. Hence, the exclusion 2-periodic cycle on the x-axis becomes unstable and a new 4-periodic cycle in the x-axis is born. Further simulations suggest that this scenario of period-doubling bifurcation, in the x-axis, will continue route to chaos.We remark that this period-doubling scenario, in three-dimensional mappings, is just a direct observation from the simulations since there is no extension of Sharkovsky’s theorem [25, 26] for dimension k≥2.It remains to study the stability of the fixed pointE1 when r1=2, r2=2c21, and r3<2c31. Under this restriction on the parameters we have λ1=-1, λ2=1, and |λ3|<1. The reduced equation on the centre manifold will be a two-dimensional system. Unfortunately, this scenario in two-dimensional discrete systems (λ~1=-1,λ~2=1) remains as an open problem, and thus we are not able to discuss the stability of the fixed point here.


## 5.2. Fixed Points on the Planes

Let us now consider the fixed pointE12 on the plane x1Ox2. The stability analysis of E13 and E23 follows in a similar way.The Jacobian evaluated at the fixed point is given by(47)JFE12=1-c12c21+c12r2-r11-c12c21c12c12r2-r11-c12c21c13c12r2-r11-c12c21c21c21r1-r21-c12c211-c12c21+r1c21-r21-c12c21c23c21r1-r21-c12c2100ec21c32-c31r1+c31c12-c32r2+1-c21c12r3/1-c12c21.It is a straightforward computation to see that the conditions of local stability are given by(48)4c12c21-1+21-c12r2+21-c21r1<c12r2-r1c21r1-r2<1-c12r2+1-c21r1,(49)c21c32-c31r1+c31c12-c32r2+1-c21c12r3<0.In Figure 5 a prototype of this region is depicted in the triple (r1,r2,r3) when cji=0.5, for i≠j.Figure 5
(a) Stability region of the exclusion fixed pointE12 in the parameter space Or1r2r3 when cji=0.5 for i≠j. The projection of this region in the plane r1Or2 is the region S1 depicted in Figure 3. (b) Stability regions of E12, E13, and E23, together in the parameter space Or1r2r3 when cii=1 and cij=0.5, for i≠j and i,j=1,2,3.
(a)
(b)Considering the equality(50)4c12c21-1+21-c12r2+21-c21r1=c12r2-r1c21r1-r2and keeping the other relations, the eigenvalues of JF(E12) are such that λ1=-1 and |λi|<1, for i=2,3. Hence, we are in position to use the centre manifold theorem. Following the techniques employed in the previous subsection, one can conclude that E12 is a locally asymptotically stable fixed point on the plane x1Ox2 under these restrictions on the parameters. Notice that, in the parameter space, the projection of this region in the r1Or2 plane is part of a branch of an hyperbola as described in the beginning of this section and can be observed in Figure 3 by γ1.This means that this specific region undergoes a period-doubling bifurcation on thex1Ox2 plane. Hence, as the parameters increase, the fixed point E12 becomes unstable and a locally asymptotically stable 2-periodic cycle of the form(51)x¯0,y¯0,0,x¯1,y¯1,0is born on the x1Ox2 plane. As in the precedent case, from now on, computations are only possible via numerical simulations. These simulations suggest a scenario of period-doubling bifurcation on the x1Ox2 plane route to chaos.We should mention that if we replace the inequality signs by equality on (48) and condition (49) is satisfied, then λ1=-1, λ2=1, and |λ3|<1 are the eigenvalues of JF(E12). Unfortunately, we are not able to use here the centre manifold theorem, as mentioned at the end of the previous subsection.


## 5.3. Positive Fixed Point

We now consider the positive fixed pointE∗ known in population dynamics as the coexistence fixed point. Due to the great number of parameters involved and for the sake of simplicity, we will study the symmetric case, that is, when cij=c, for all i≠j. The fixed point E∗=(x1∗,x2∗,x3∗) becomes in this case(52)E∗=cr2+r3-c+1r1c-12c+1,cr1-c+1r2+cr3c-12c+1,cr1+cr2-c+1r3c-12c+1,and the derivative at E∗ is given by(53)JFE∗=-2c2+r2c+r3c+c-c+1r1+1-2c2+c+1cc+1r1-cr2+r3c-12c+1cc+1r1-cr2+r3c-12c+1c-cr1+c+1r2-cr3c-12c+1-2c2+r1c+r3c+c-c+1r2+1-2c2+c+1c-cr1+c+1r2-cr3c-12c+1c-cr1-cr2+c+1r3c-12c+1c-cr1-cr2+c+1r3c-12c+1-2c2+r1c+r2c-r3c+c-r3+1-2c2+c+1.For this matrix we have(54)tr⁡JFE∗=-6c+r1+r2+r3-32c+1,∑i=13MiiJFE∗=3P0+2P1r1+r2+r3+P2r12+r22+r32P0+P4r1r2+r1r3+r2r3P0,det⁡JFE∗=-P0+∑i=13Pir1i+r2i+r3i+P4r1r2+r1r3+r2r3P0-P5r1r22+r32+r2r12+r32+r3r12+r22+P6r1r2r3P0,where (55)P0=4c3-3c-1,P1=-2c2+c+1,P2=c3+3c2+2cP3=c3+c2,P4=-2c3+2c2+c+1,P5=-c3+c2+cP6=1+3c+6c2+2c3.The conditions of local stability ofE∗ now follow from relation (9). A prototype of this region, in the parameter space Or1r2r3, is represented in Figure 6. Notice that this region is a three-dimensional solid with planar faces, with the exception of the top, where the surface is given by(56)tr⁡JFE∗+det⁡JFE∗-∑i=13MiiJFE∗=1.Figure 6
Stability region of the positive fixed pointE∗ in the parameter space Or1r2r3 when cii=1 and cji=0.5, for i≠j and i,j=1,2,3.Now, in conditions (9), when tr⁡(JF(E∗))+det⁡(JF(E∗))-∑i=13Mii(JF(E∗))=1 and assuming the others inequalities, we have λ1=-1 and |λi|<1, for i=2,3. Hence, we have to use the centre manifold theory in order to determine the stability of E∗.Making the change of variables(57)y1=x1-x1∗,y2=x2-x2∗,y3=x3-x3∗, we shift the fixed point to the origin. System (28) is now equivalent to(58)Yn+1=JYn+GYn,where(59)Yn=y1ny2ny3n,J=J11J12J13J21J22J23J31J32J33,GYn=g1Yng2Yng3Yn,J=-2c2+cr2+r3+c+11-r1-2c2+c+1cc+1r1-cr2+r3c-12c+1cc+1r1-cr2+r3c-12c+1-ccr1-c+1r2+cr3c-12c+1-2c2+r1c+r3c+c-c+1r2+1-2c2+c+1-ccr1-c+1r2+cr3c-12c+1-ccr1+cr2-c+1r3c-12c+1-ccr1+cr2-c+1r3c-12c+1-2c2+r1c+r2c-r3c+c-r3+1-2c2+c+1,giy1,y2,y3=fiy1+x1∗,y2+x2∗,y3+x3∗-xi∗-∑j=13Jijyj,i=1,2,3.SinceJ is nonsingular, there exist a diagonal matrix Λ and the matrices P and P-1 such that J=PΛP-1. Making a new change of variables Y=PZ, where Z=(z1,z2,z3), one can write (58) as(60)Zn+1=ΛZn+HZn,where(61)Zn=z1nz2nz3n,Λ=-1000λ2000λ3,HZn=P-1GP×Zn.The centre manifold in this case can be represented locally as the graph of a functionΦ:R2→R2 such that (62)Wc=z1,z2,z3∈R×R2∣z2,z3=Φ1z1,Φ2z1,Φ0=0,DΦ0=0,z1<δfor  a  sufficiently  smallδ.Therefore, the centre manifold equation (27) is given by (63)Φ1-z1+h1z1,Φ1z1,Φ2z1=λ2Φ1z1+h2z1,Φ1z1,Φ2z1Φ2-z1+h1z1,Φ1z1,Φ2z1=λ3Φ2z1+h3z1,Φ1z1,Φ2z1,where Φi, i=1,2, can be represented locally by Φi(u)=αiu2+βiu3+O[u]4 and (64)hiz1,z2,z3=∑j=13P~ijgjz1,z2,z3,i=1,2,3,where P~ij are the entries of the matrix P-1.Unfortunately, at this stage the computations are long and we are not able to find explicitly the matricesP, Λ, and P-1 such that J=PΛP-1. However, we are able to find it numerically. Hence, we will illustrate the study with a concrete example. Let c=0.5 and r2=r3=1.5 (the choice of these values is made considering the region obtained in Figure 6). Substituting in(65)tr⁡JFE∗+det⁡JFE∗-∑i=13MiiJFE∗=1, we obtain r1≈2.1547. Under these values of the parameters, the matrices are given by(66)P=-0.97690-0.6676-0.1511-0.70710.5265-0.15110.70710.5265,Λ=-1.0000.78870000.634,P-1=-0.8558-0.5425-0.54250-0.70710.7071-0.24570.7940.794.And the functions hi(z1,z2,z3), for i=1,2,3, are (67)h1z1,z2,z3=1.9409+z1+0.836z1+0.5713z3-1.4823e1.128z1+0.1411z3+0.082z1-0.3836z2-0.2856z3-0.2293e0.7151z1-0.3536z2-0.4559z3+0.082z1+0.3836z2-0.2856z3-0.2293e0.7151z1+0.3536z2-0.4559z3,h2z1,z2,z3=-0.7887z2+0.1069z1+0.5z2-0.3723z3-0.2989e0.7151z1+0.3536z2-0.456z3+-0.1069z1+0.5z2+0.3723z3+0.2989e0.7151z1-0.3536z2-0.4559z3,h3z1,z2,z3=-0.634z3-0.2456+0.24z1+0.164z3-0.4255e1.128z1+0.1411z3+0.3356-0.12z1+0.5614z2+0.418z3e0.7151z1-0.3536z2-0.456z3+-0.12z1-0.5614z2+0.418z3+0.336e0.7151z1+0.354z2-0.456z3.Using Taylor series and solving system (63), one can conclude that the dynamics of the system is given locally by the dynamics of z1(n+1)=h~(z1(n)), where(68)h~z1=e0.7151z10.164z1-0.4586+e1.128z10.836z1-1.4823+1.9409.Notice that this corresponds to the equation on the centre manifold that is generically given by (15).Sinceh~′(0)=-1, we have to use Theorem 4. The Schwarzian derivative of h~ at 0 is ≈-1.1476<0. It follows from Theorem 4 that, in this particular case, the fixed point is locally asymptotically stable.We remark that, due to the difficulties encountered in finding the matrixP, one can write a script in Mathematica or Maple to automatically obtain, via numerical approximation, the map h~ in each case.


## 6. Conjecture

One of the main questions in population dynamics isunder what conditions does local stability of a fixed point imply its global asymptotic stability?One of the most influential models in population dynamics is the one-dimensional Ricker equation given byxn+1=xner-xn, r>0. The local stability condition of the positive fixed point x∗=r, 0<r≤2, implies global stability of x∗; that is, all solutions of the difference equation xn+1=xner-xn, starting at an initial condition x0>0, converge to x∗=r whenever 0<r≤2. In order to establish this conclusion one can use Coppell’s Theorem since x∗=r is the unique positive fixed point of the map.Theorem 14 (see Coppell [27]).
LetI=[a,b]⊆R and f:I→I be a continuous map. If the equation ffx=x has no roots, with the possible exception of the roots of the equation fx=x, then every orbit under the map f converges to a fixed point.In planar models this question of global stability is more involved since we do not have an analogue of Coppell’s Theorem, but alternative methods do exist. For instance, Balreira et al. [28], used the Singularity Theory of Whitney [29] and proved that thelocal stability conditions (32) of the positive fixed point (69)x1∗,x2∗=r1-c12r21-c12c21,r2-c21r11-c12c21in the two-dimensional Ricker competition model (31) imply its global stability provided that 1<r1, r2<2, and the image of a critical curves [30] do not intersect. Later, Ryals and Sacker [31] improved the result by reducing the infinite number of topological conditions on the critical curves to a finite number.We should mention that Smith [32] uses monotonicity to prove the global stability of the fixed points of the system(70)un+1=unexpr1-un-Bvnvn+1=vnexps1-Cun-vn,when r,s≤1, in which case the invariant set is [0,r-1]×[0,s-1]. Notice that, by the changes of variables ru=x and sv=y, system (70) is equivalent to(71)xn+1=xnexpr-xn-Brsynyn+1=ynexps-yn-Csrxn.Clearly, r1=r, r2=s, c12=Br/s, and c21=Cs/r. Therefore, Smith’s global results cover the local analysis when the carrying capacities are taken in the unit interval and, in this case, global stability is obtained when ri≤1, for i=1,2, such that ri∈S1 (Figure 3).It remains as an open problem to show the global stability in the mixing cases, that is, whenr1≤1 and 1<r2≤2, or r2≤1 and 1<r1≤2, and also in the region between r1=2 or r2=2 and the branch of hyperbola γ1 as is depicted in Figure 3, as part of the region S1. This is precisely part of the conjecture raised by Elaydi and Luís [33].We should mention that the study of global stability can also be extended to the other nontrivial fixed points (by trivial fixed point we mean the origin) of system (31), as well.Now, motivated by the precedent cases in lower dimension, we have the following conjecture concerning the nontrivial fixed points of the three-dimensional Ricker competition model.Conjecture 15.
Consider the Ricker competition model (28) where the growth rates ri are positive and the competition parameters are taken as cii=1 and 0<cji<1, for i≠j and i,j=1,2,3. Let x∗ be a nontrivial fixed point of the map F defined in (29) and assume that we have(72)-1-∑i=13MiiJFx∗<tr⁡JFx∗+det⁡JFx∗≤1+∑i=13MiiJFx∗∑i=13MiiJFx∗-tr⁡JFx∗det⁡JFx∗<1-det2⁡JFx∗,where tr⁡(JF(x∗)), det⁡(JF(x∗)), and ∑i=13Mii(JF(x∗)) are the trace, the determinant, and the sum of the principal minors of the Jacobian matrix evaluated at the fixed point x∗, respectively.
Then, the fixed pointx∗ is globally asymptotically stable with respect to the interior of the first octant.Notice that, for each one of the nontrivial fixed pointsx∗ mentioned in the conjecture (three fixed points on the axes E1, E2, and E3, three fixed points on the planes E12, E13, and E23, and the positive fixed point E∗), the local stability conditions (72) correspond to those developed in the previous section for the respective fixed point. Moreover, geometrically it corresponds to each one of the seven regions depicted in Figure 7 as mentioned before.Figure 7
Stability regions, in the parameter spaceOr1r2r3, of the fixed points of the three-dimensional Ricker competition model when cii=1 and cji=0.5, for i≠j and i,j=1,2,3.Before ending this section, we should mention that in a forthcoming paper, Balreira et al. [34] use a new notion of monotonicity in higher dimensional maps and show that Conjecture 15 is true in the case of the positive fixed point E∗, when ri<1 for i=1,2,3. The result was established since the three-dimensional Ricker competition model is monotone (more specifically, normally monotone) if conditions (72) are valid such that ri<1 for i=1,2,3. This is precisely a subset of the set depicted in Figure 6.Finally, if one considers the remaining region in the parameter space, the map is not monotone. Since it is possible to show global stability, in the planar model, for certain cases (Balreira et al. [28] and Ryals and Sacker [31]), using the Singularity Theory of Whitney [29], we believe that an extension of the ideas may be useful in proving Conjecture 15. Unfortunately, such theory of singularity does not exist for three-dimensional systems. It either is necessary to develop such theory of singularity in three-dimensional models to show Conjecture 15 or may use another different tool such as Lyapunov method (cf. Theorem4.2 in Elaydi [20, page 205]).


## 7. Conclusion

In this paper we illustrate how to apply the centre manifold theorem to three-dimensional discrete dynamical systems. We recall the necessary and sufficient conditions of local stability of three-dimensional systems and give, in depth, the tools to study the stability of nonhyperbolic fixed points in the terminology of difference equations.As illustrative example, we study the Ricker competition model of three species. This nonlinear model has seven nontrivial fixed points. Each fixed point has a region, in the parameter space, where it is locally asymptotically stable. When one determines these regions, the result is the seven regions depicted in Figure7.In addition, we state a conjecture in the global stability of each nontrivial fixed point of the Ricker competition model, in the respective region of local stability.



---

*Source: 6186354-2017-07-06.xml*