Interesting Features of Three-Dimensional Discrete Lotka-Volterra Dynamics

Yogesh Joshi¹, Micelle Savescu², Musa Syed³, Denis Blackmore⁴

¹Department of Mathematics and Computer Science, Kingsborough Community College, Brooklyn, NY, USA.
²Department of Mathematics, Kutztown University, Kutztown, PA, USA.
³Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ, USA.
⁴Department of Mathematical Sciences and Center for Applied Mathematics and Statistics, New Jersey Institute of Technology, Newark, NJ, USA.

DOI: 10.4236/am.2021.128049                

Abstract

Discrete Lotka-Volterra systems in one dimension (the logistic equation) and two dimensions have been studied extensively, revealing a wealth of complex dynamical regimes. We show that three-dimensional discrete Lotka-Volterra dynamical systems exhibit all of the dynamics of the lower dimensional systems and a great deal more. In fact and in particular, there are dynamical features including analogs of flip bifurcations, Neimark-Sacker bifurcations and chaotic strange attracting sets that are essentially three-dimensional. Among these are new generalizations of Neimark-Sacker bifurcations and novel chaotic strange attractors with distinctive candy cane type shapes. Several of these dynamical are investigated in detail using both analytical and simulation techniques.

Keywords

Discrete Lotka-Volterra Systems, Flip Bifurcations, Higher Dimensional Neimark-Sacker Type Bifurcations, Chaotic Strange Attracting Sets, Horseshoe Type Dynamics

1. Introduction

We shall consider a natural discrete analog of the famous and ubiquitous Lotka-Volterra model, developed independently by Alfred Lotka [1] and Vito Volterra [2] for the dynamics of a population comprising several interacting species, say $x\mathrm{<mo>\; :\; </mo>}=\left(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_m\right)$, as described in such references as [3] [4]. This model may be expressed in the form

${\stackrel{\dot{}}{x}}_i\mathrm{<mo>\; :\; </mo>}=\frac{\text{d}{x}_i}{\text{d}t}={\alpha }_i{x}_i\left(1-{\displaystyle {\sum }_{j=1}^m}{\beta }_{ij}{x}_j\right)\mathrm{,}1\le i\le m\mathrm{,}$ (1)

where the coefficients ${\alpha }_i\mathrm{,}{\beta }_{ij}$ are real constants and typically nonnegative, with ${\beta }_{ii}=1$ for all $1\le i\le m$, $1\le j\le n$. Note that for a single species, (1) is essentially the (continuum) logistic equation. If one employs the Euler method to obtain discrete approximate solutions of (1), the result can be written as the following system of difference equations:

$x_{i,n+1}={\lambda }_i{x}_{i,n}\left(1-{\displaystyle {\sum }_{j=1}^m}{\gamma }_{ij}{x}_{j,n}\right),1\le i\le m,$ (2)

where the birth-rate parameters ${\lambda }_i$ and interaction coefficients ${\gamma }_{ij}$ are usually nonnegative constants as we shall assume in the sequel, with ${\gamma }_{ii}=1$ for all $1\le i\le m$, $1\le j\le n$. Analogous to our remark above, we note that the 1-dimensional version of (2) is the discrete logistic equation that, in contrast to (1) for $1\le m\le 2$, can exhibit chaotic dynamics.

Although our focus here will be on the 3-dimensional version of (2), we first describe a few useful details for the m-dimensional system, which has dynamics defined by iterates of the map $F\mathrm{<mo>\; :\; </mo>}{ℝ}^m\to {ℝ}^m$ comprising a discrete (semi-) dynamical system and is represented as

$\begin{array}{c}F\left(x\right)=F\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_m\right)\mathrm{<mo>\; :\; </mo>}=\left(f_1\left(x\right)\mathrm{,}\cdots \mathrm{,}{f}_m\left(x\right)\right)\\ =\left(f_1\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_m\right)\mathrm{,}\cdots ,f_m\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_m\right)\right)\\ =\left({\lambda }_1{x}_1\left(1-{\displaystyle {\sum }_{j=1}^m}{\gamma }_{1j}{x}_j\right),\cdots ,{\lambda }_m{x}_m\left(1-{\displaystyle {\sum }_{j=1}^m}{\gamma }_{mj}{x}_j\right)\right),\end{array}$ (3)

where the coefficients are as described in (2). Here, the $\lambda$ ’s and $\gamma$ ’s represent the birth-rates and interaction strengths of and among the populations, respectively.

Inasmuch as we are primarily thinking of $x$ as a population vector, it is natural to consider the map F restricted to

${\mathfrak{N}}^m\mathrm{<mo>\; :\; </mo>}=\left\{x\in {ℝ}^m\mathrm{<mo>\; :\; </mo>}{x}_1\mathrm{,}\cdots \mathrm{,}{x}_m\ge 0\right\}\mathrm{.}$

Moreover, since we want the forward iterates of F of points in this set to remain in ${\mathfrak{N}}^m$, we would like to find the largest subset $\mathfrak{X}={\mathfrak{X}}_F$ such that $F\left(\mathfrak{X}\right)\subset \mathfrak{X}$ for sizeable ranges of the birth-rate parameters. For example, if all the off-diagonal interaction coefficients are zero, the obvious choice of $\mathfrak{X}$ is the m-cube $K^m={\left[\mathrm{0,1}\right]}^m\mathrm{<mo>\; :\; </mo>}=\left\{x\in {ℝ}^m\mathrm{<mo>\; :\; </mo>0}\le x_i\le \mathrm{1,1}\le i\le m\right\}$ in which case $F\left(\mathfrak{X}\right)\subset \mathfrak{X}$ for all $0\le {\lambda }_i\le 4$, $1\le i\le m$. For convenience, in this case we say the birth-rate parameter ranges are admissible. When some of the off-diagonal interaction coefficients are positive, the manifest generalization of the invariant set is

$\mathfrak{X}={\mathfrak{X}}_F:=\left\{x\in {ℝ}^m\mathrm{<mo>\; :\; </mo>0}\le x_i\le 1-{\displaystyle \underset{j\left(\ne i\right)=1}{\overset{m}{\sum }}}{\gamma }_{ij}{x}_j\right\}$ (4)

and the admissible ${\lambda }_i$ ranges have to be adjusted accordingly.

It is clear from (3) that a necessary and sufficient condition for $F\left(\mathfrak{X}\right)\subset \mathfrak{X}$ is that the following system of inequalities be satisfied:

$0\le {\phi }_i\left(x\mathrm{,}\lambda \right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\sum }_{j=1}^m}{\gamma }_{ij}{f}_j\left(x\right)\le \mathrm{1,}1\le i\le m\mathrm{,}$ (5)

where $\lambda \mathrm{<mo>\; :\; </mo>}=\left({\lambda }_1\mathrm{,}\cdots \mathrm{,}{\lambda }_m\right)$. One can by lengthy but straightforward calculations find the maximizers and associated maximum values of each of the ${\phi }_i$, which yield a rather complicated system of inequalities determining the optimal admissible values of $\lambda$. However, the system is complex to the point of being very unwieldy, we shall opt for a cruder but much easier to use system of inequalities described in the following result.

Lemma 1. The birth-rate parameters satisfying the system of inequalities

${\displaystyle {\sum }_{j=1}^m}{\gamma }_{ij}{\lambda }_j\le \mathrm{4,}1\le i\le m\mathrm{,}$ (6)

are admissible in $\mathfrak{X}$ for the discrete dynamical system generated by the iterates of the map (3).

Proof. It follows from (4) and (5) that $0\le {\phi }_i\left(x\mathrm{,}\lambda \right)$ for all $x\in \mathfrak{X}$ and $\lambda \in {\mathfrak{N}}^m$, $1\le i\le m$, so only the upper inequalities in (6) require verification. For this, we observe that it is easy to see that each $f_j\left(x\right)\le {\lambda }_j/4$ on $\mathfrak{X}$. Consequently,

${\phi }_i\left(x\mathrm{,}\lambda \right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\sum }_{j=1}^m}{\gamma }_{ij}{f}_j\left(x\right)\le \frac{1}{4}{\displaystyle {\sum }_{j=1}^m}{\gamma }_{ij}{\lambda }_j\mathrm{,}$

which implies that (5) follows from (6), thereby completing the proof.£

The generalities discussed above aside, our almost exclusive focus shall be on the discrete dynamical system generated by the iterates of the map $F\mathrm{<mo>\; :\; </mo>}{ℝ}^3\to {ℝ}^3$ defined as

$\begin{array}{l}F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(f_1\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_2\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_3\left(x\mathrm{,}y\mathrm{,}z\right)\right)\\ =\left({\lambda }_{1}x\left(1-x-{\gamma }_{12}y-{\gamma }_{13}z\right)\mathrm{,}{\lambda }_{2}y\left(1-{\gamma }_{21}x-y-{\gamma }_{23}z\right)\mathrm{,}{\lambda }_{3}z\left(1-{\gamma }_{31}x-{\gamma }_{32}y-z\right)\right).\end{array}$ (7)

Clearly, the x, y-, x, z- and y, z -planes as well as the x-, y- and z-axes are F-invariant, so all of the interesting dynamics for the 1-dimensional and 2-dimensional cases of (2)-(3) including flip bifurcations, period-doubling cascades (leading to chaos), Neimark-Sacker bifurcations, (2-dimensional) transverse heteroclinic orbit induced chaos and chaotic strange attractors of Hausdorff dimension between one and two such as found in [5] [6] [7] [8] are subsumed under the dynamics defined by (7). Naturally then, one can expect even more complex and interesting higher dimensional analogs of such dynamics in the 3-dimensional case, which shall be elucidated in the sequel. In particular, we shall find and analyze flip bifurcations, certain higher dimensional Neimark-Sacker type bifurcations to be described in the sequel, several 3-dimensional chaotic regimes and a few unusual chaotic strange attractors corresponding to long-term population dynamic states. This includes an apparently new type of chaotic strange attractor shaped like a candy cane. As one might expect, there have been several studies of the dynamics of higher dimensional systems of the Lotka-Volterra (L-V) and related types such as [9] - [15], but our investigation is novel in a number of respects, especially with regard to the Neimark-Sacker bifurcation generalizations and the candy cane chaotic strange attractors. Nevertheless, our work shares some common elements with the dynamics literature. For example, Bischi and Tramontana [10] show that their 3-dimensional Lotka-Volterra type model for industrial clusters exhibits standard flip and Neimark-Sacker bifurcations along with interesting chaotic attractors, which are not candy canes. Another example in Yousef et al. [15], although fundamentally different from our model owing to the delay in one of the growth rates, also is shown to exhibit flip and Neimark-Sacker bifurcations along with Marotto chaos.

The investigation is organized as follows: In Section 2, we describe some basic features of the discrete dynamical system generated by the map (7) such as fixed points and invariant manifolds. One of our foci will be various types of bifurcations that are fully 3-dimensional versions of those found in lower dimensional systems, and one of them—the flip bifurcation shall be investigated in some depth in Section 3. This is followed in Section 4 with a definition of certain higher dimensional analogs of Neimark-Sacker bifurcations and the identification of various coefficient/parameter ranges for which they are exhibited for (7). In addition, several simulations are included to illustrate these bifurcations. Next, in Section 5 the focus is on chaotic strange attractors. We show simulation examples indicating the wide variety of such steady-state sets. In particular, we provide examples where the nature of the steady-state attractor is strongly dependent on the parameters and auxiliary data associated to the dynamical system. Most importantly, however, we introduce new types of chaotic strange attractors, which we analyze in detail. Finally, in Section 6, we present conclusions and plans for future related research including an in depth analysis of the higher dimensional Neimark-Sacker type bifurcations and candy cane attractors, which appear to be quite common in population dynamics.

2. Basic Dynamical Properties

It is clear from its definition that the map (7) enjoys the following invariance properties concerning the coordinate axes denoted as $L_x\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,0,0}\right)\mathrm{<mo>\; :\; </mo>}x\in ℝ\right\}$, $L_y\mathrm{<mo>\; :\; </mo>}=\left\{\left(\mathrm{0,}y\mathrm{,0}\right)\mathrm{<mo>\; :\; </mo>}y\in ℝ\right\}$ and $L_z\mathrm{<mo>\; :\; </mo>}=\left\{\left(\mathrm{0,0,}z\right)\mathrm{<mo>\; :\; </mo>}z\in ℝ\right\}$ and the coordinate planes by $P_{xy}\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,0}\right)\mathrm{<mo>\; :\; </mo>}x\mathrm{,}y\in ℝ\right\}$, $P_{xz}\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,0,}z\right)\mathrm{<mo>\; :\; </mo>}x\mathrm{,}z\in ℝ\right\}$ and $P_{yz}\mathrm{<mo>\; :\; </mo>}=\left\{\left(\mathrm{0,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}y\mathrm{,}z\in ℝ\right\}$ :

Theorem 2. The origin $\left(\mathrm{0,0,0}\right)$ is a fixed point of F and all of the coordinate sets defined above are F-invariant, i.e. $F\left(L_x\right)\subset L_x$, $F\left(L_y\right)\subset L_y$, $F\left(L_z\right)\subset L_z$, $F\left(P_{xy}\right)\subset P_{xy}$, $F\left(P_{xz}\right)\subset P_{xz}$ and $F\left(P_{yz}\right)\subset P_{yz}$.

If any of the coordinates of the initial point for the discrete dynamics for (4) is zero, it follows from Theorem 1 that the whole process reduces to one of two dimensions or less, where it has been shown that the system can exhibit flip bifurcations, period-doubling cascades to chaos, horseshoe type chaos and a variety of other dynamical phenomena, so we shall focus on initial points with no zero coordinates, which are sometimes referred to as coexistence points, mainly in the context of population dynamics. Thus, it may be said that we are restricting our attention to fully three-dimensional analogs of such dynamical behavior. In particular, we shall be most interested in fixed points of the map (7) that are also coexistence points. These fixed points can be found by solving the system of linear equations

$\begin{array}{l}x+{\gamma }_{12}y+{\gamma }_{13}z=\frac{{\lambda }_1-1}{{\lambda }_1}\mathrm{,}\\ {\gamma }_{21}x+y+{\gamma }_{23}z=\frac{{\lambda }_2-1}{{\lambda }_2}\mathrm{,}\\ {\gamma }_{31}x+{\gamma }_{32}y+z=\frac{{\lambda }_3-1}{{\lambda }_3}\mathrm{,}\end{array}$ $\iff \Gamma {\left(x\mathrm{,}y\mathrm{,}z\right)}^{\text{T}}={\left(\frac{{\lambda }_1-1}{{\lambda }_1}\mathrm{,}\frac{{\lambda }_2-1}{{\lambda }_2}\mathrm{,}\frac{{\lambda }_3-1}{{\lambda }_3}\right)}^{\text{T}}\mathrm{,}$ (8)

which has a unique solution if the matrix $\Gamma$ is nonsingular.

We note that it is sometimes convenient to recast the discrete dynamical system in terms of (translated coordinates) with respect to a fixed point $\left(x_{\ast }\mathrm{,}{y}_{\ast }\mathrm{,}{z}_{\ast }\right)$ with $x_{\ast }\mathrm{,}{y}_{\ast }$ and $z_{\ast }$ positive, namely with $\xi :=x-x_{\ast }$, $\eta :=y-y_{\ast }$ and $\zeta \mathrm{<mo>\; :\; </mo>}=z-z_{\ast }$, so that the map takes the form

$\hat{F}\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}=F\left(\xi +x_{\ast }\mathrm{,}\eta +y_{\ast }\mathrm{,}\zeta +z_{\ast }\right)-\left(x_{\ast }\mathrm{,}{y}_{\ast }\mathrm{,}{z}_{\ast }\right)\mathrm{,}$ (9)

with corresponding fixed point $\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)=\left(\mathrm{0,0,0}\right)$.

3. Flip Bifurcations and Period-Doubling Cascades

Flip bifurcations and associated period-doubling cascades to chaos are well known types of bifurcations and behaviors thoroughly explicated in such texts as [16] [17], which have been shown to exist in discrete Lotka-Volterra (L-V) dynamical systems associated with the map (7) and its higher and lower dimensional versions. There are numerous parameter sets for which (7) exhibits flip bifurcations and subsequent period-doubling cascades leading to one-dimensional logistic map type chaos, several of which are described in this section. We begin with some rather simple examples and then analyze a less obvious case.

3.1. Some Simple Flip and Period-Doubling Examples

Consider the special case chosen for its ease of analysis of the family of maps (7) given as $F_{\lambda }\mathrm{<mo>\; :\; </mo>}\mathfrak{X}\to \mathfrak{X}$ defined by

$F_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\left(\lambda -1\right)x\left(1-x-\alpha y\right)\mathrm{,}\left(\lambda -1\right)y\left(1-\beta x-y\right)\mathrm{,}\lambda z\left(1-z\right)\right)\mathrm{,}$

$\mathfrak{X}={\mathfrak{X}}_{F_{\lambda }}\mathrm{<mo>\; :\; </mo>}=\left\{x=\left(x\mathrm{,}y\mathrm{,}z\right)\in {ℝ}^3\mathrm{<mo>\; :\; </mo>0}\le x\le 1-\alpha y\mathrm{,0}\le y\le 1-\beta x\mathrm{,0}\le z\le 1\right\}\mathrm{,}$ (10)

with $0\le \alpha ,\beta <0.3$ and the single (bifurcation) parameter $\lambda \in \left(\mathrm{2,4}\right]$. It follows from Lemma 1 that $0\le \lambda \le 4$ is an admissible range for the above map. The coexistence fixed points of this map as functions of the parameter $\lambda$ are readily found to be

$p_{\lambda }\mathrm{<mo>\; :\; </mo>}=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)=\left(\left(\frac{1-\alpha }{1-\alpha \beta }\right)\frac{\lambda -2}{\lambda -1}\mathrm{,}\left(\frac{1-\beta }{1-\alpha \beta }\right)\frac{\lambda -2}{\lambda -1}\mathrm{,}\frac{\lambda -1}{\lambda }\right)\mathrm{.}$ (11)

This class of L-V maps has the bifurcation behavior described in the following result, which has a straightforward proof that is left to the reader. The dynamics are also illustrated in Figure 1.

Theorem 3. The discrete dynamical system generated by the map (10) satisfies the following properties when $0\le \alpha ,\beta <0.3$ and $\lambda \in \left(\mathrm{2,4}\right]$ :

i) The system $F_{\lambda }$ has a flip bifurcation at $p_{\lambda }$ for $\lambda =3$ at which the restriction $f\mathrm{<mo>\; :\; </mo>}=F_{\lambda \mathrm{<mo>\; |\; </mo>}{L}_{\lambda }}$ to the $F_{\lambda }$ -invariant line (segment) $L_{\lambda }\mathrm{<mo>\; :\; </mo>}=\left\{\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\in ℝ\right\}\cap \mathfrak{X}$ has derivative $f^{\prime }\left(z_{\lambda }\right)=-1$.

ii) The restriction f has a period-doubling cascade leading to chaos along $L_{\lambda }$ as $\lambda$ increases from 3 to a limit of approximately 3.57 following that of the standard 1-dimensional logistic map.

iii) The fixed point $p_{\lambda }$ is a sink for $2<\lambda <3$ and is hyperbolic for $3<\lambda <4$ with a 2-dimensional (planar) stable manifold $W^s\left(p_{\lambda }\right)=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$ and 1-dimensional (linear) unstable manifold $W_{\text{lin}}^u\left(p_{\lambda }\right)=\left\{\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\in ℝ\right\}$, where $\stackrel{\circ }{\mathfrak{X}}$ denotes the interior of $\mathfrak{X}$.

iv) When $\lambda =3$, the fixed point $p_{\lambda }=p_3$ has the stable manifold $W^s\left(p_3\right)=\left\{\left(x\mathrm{,}y\mathrm{,}2/3\right)\mathrm{<mo>\; :\; </mo>\; <mo\; stretchy="false">\; (\; </mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">\; )\; </mo>}\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$ and the center manifold $W^c\left(p_3\right)=\left\{\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\in \left(\mathrm{0,1}\right)\right\}$.

v) There is a Milnor attractor ${\mathcal{A}}_{\infty }$ along $W_{\text{lin}}^u\left(p_{\lambda }\right)\cap \mathfrak{X}$, as in [18], at the limit of the period-doubling cascade as well as a dense chaotic (standard tent map-type Milnor) attractor ${\mathcal{A}}_T$ for $\lambda =4$, each with a basin of attraction $\stackrel{\circ }{\mathfrak{X}}\backslash \left\{\text{denumerable}\text{subset}\text{of}\text{periodic}\text{points}\text{on}{W}_{\text{lin}}^u\left(p_{\lambda }\right)\cap \mathfrak{X}\right\}$.

It should be noted that by permuting the coordinates, two analogs of (10) with flips along lines parallel to the x- and y-axes are obtained. Also observe that these examples can be readily extended to higher-dimensional systems comprising an arbitrary number of populations.

Figure 1. Flip bifurcation for (10): (a) $\lambda =2.8$ ; (b) $\lambda =3.4$ ; (c) $\lambda =3.56$ ; (d) $\lambda =4$ with $\alpha =0.2$ and $\beta =0.3$.

3.2. More Complicated Flip and Period-Doubling Examples

There are numerous other flip bifurcation/period-doubling cascades to chaos that can occur in the L-V discrete dynamical systems generated by maps of the form (4), but which are not nearly as recognizable as those described above. We shall confine our attention to just the one example $G_{\lambda }\mathrm{<mo>\; :\; </mo>}\mathfrak{X}\to \mathfrak{X}$, with equal birth rates serving as the bifurcation parameter (with $2<\lambda \le 3.6$ ), shown below. The dynamics shall be analyzed in detail for prescribed ranges of variation of the interaction coefficients for

$G_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{<mo\; stretchy="false">\; (\; </mo>}\lambda x\left(1-x-a_{1}y-b_{1}z\right)\mathrm{,}\lambda y\left(1-a_{2}x-y-b_{2}z\right)\mathrm{,}\lambda z\left(1-cx-dy-z\right)\right)\mathrm{,}$

$\mathfrak{X}={\mathfrak{X}}_{G_{\lambda }}\mathrm{<mo>\; :\; </mo>}=\left\{x\mathrm{<mo>\; :\; </mo>0}\le x\le 1-a_{1}y-b_{1}z\mathrm{,0}\le y\le 1-a_{2}x-b_{2}z\mathrm{,0}\le z\le 1-cx-dy\right\}\mathrm{.}$ (12)

For the coexistence fixed point, we simply solve the following system of equations in which we impose the initial constraints $0\le a_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\le 0.1$ :

$x+a_{1}y+b_{1}z=\sigma \left(\lambda \right):=\frac{\lambda -1}{\lambda },$

$a_{2}x+y+b_{2}z=\sigma \left(\lambda \right):=\frac{\lambda -1}{\lambda },$ (13)

$cx+dy+z=\sigma \left(\lambda \right):=\frac{\lambda -1}{\lambda },$

to obtain the unique solution

$p_{\lambda }\mathrm{<mo>\; :\; </mo>}=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)\mathrm{.}$ (14)

In order to better understand the flip bifurcations that can occur for the map (12), we first consider the following special case:

$g_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x-b_{1}z\right)\mathrm{,}\lambda y\left(1-y-b_{2}z\right)\mathrm{,}\lambda z\left(1-z\right)\right)\mathrm{,}$ (15)

with $0.05\le b_1,b_2<0.08$. It follows from Lemma 1 that $0\le \lambda \le 3.7$ is an admissible range for this map. Then, from (13) we find that the unique coexistence fixed point is

$p_{\lambda }\mathrm{<mo>\; :\; </mo>}=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)=\sigma \left(\lambda \right)\left(1-b_1\mathrm{,1}-b_2\mathrm{,1}\right)\mathrm{.}$ (16)

The type of this fixed point is determined by the eigenvalues of the derivative matrix

$\begin{array}{c}{g^{\prime }}_{\lambda }\left(p_{\lambda }\right)=\left[\begin{array}{ccc}\lambda \left(1-2x_{\lambda }-b_1{z}_{\lambda }\right)& 0& -\lambda b_1{x}_{\lambda }\\ 0& \lambda \left(1-2y_{\lambda }-b_2{z}_{\lambda }\right)& -\lambda b_2{y}_{\lambda }\\ 0& 0& \lambda \left(1-2z_{\lambda }\right)\end{array}\right]\\ =\left[\begin{array}{ccc}\left(2-\lambda \right)+b_1\left(\lambda -1\right)& 0& -b_1\left(1-b_1\right)\left(\lambda -1\right)\\ 0& \left(2-\lambda \right)+b_2\left(\lambda -1\right)& -b_2\left(1-b_2\right)\left(\lambda -1\right)\\ 0& 0& 2-\lambda \end{array}\right]\mathrm{,}\end{array}$ (17)

which are ${\mu }_1=\left(2-\lambda \right)+b_1\left(\lambda -1\right)$, ${\mu }_2=\left(2-\lambda \right)+b_2\left(\lambda -1\right)$ and ${\mu }_3=2-\lambda$, with corresponding eigenfunctions $e_1={\left(\mathrm{1,0,0}\right)}^{\text{T}}$, $e_2={\left(\mathrm{0,1,0}\right)}^{\text{T}}$ and $e_3={\left(1-b_1\mathrm{,1}-b_2\mathrm{,1}\right)}^{\text{T}}$ for $\lambda \in \left(\mathrm{1,3.7}\right]$, respectively. We note that $\left|{\mu }_1\right|,\left|{\mu }_2\right|<1$ for all $\lambda \in \left(\mathrm{2,3.1}\right]$, so that $P\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_{\lambda }\right)\mathrm{<mo>\; :\; </mo>\; <mo\; stretchy="false">\; (\; </mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">\; )\; </mo>}\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$, which is $g_{\lambda }$ -invariant for all $\lambda$, comprises the stable manifold $W^s\left(p_{\lambda }\right)$ over the same parameter interval. Moreover, $-1<{\mu }_3<0$ for $2<\lambda <3$, ${\mu }_3=-1$ when $\lambda =3$ and $-1.7\le {\mu }_3<-1$ for $3<\lambda \le 3.7$. It follows of course, that there is a flip bifurcation at $\lambda =3$ at which value the linearized center manifold is $W_{\text{lin}}^c\left(p_3\right)\mathrm{<mo>\; :\; </mo>}=\left\{z\left(1-b_1\mathrm{,1}-b_2\mathrm{,1}\right)\mathrm{<mo>\; :\; </mo>}z\in ℝ\right\}$ and the ( $g_{\lambda }$ -invariant) stable manifold is $W^s\left(p_3\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_3\right)\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$. When $3<\lambda \le 3.1$, the stable manifold is $W^s\left(p_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$ and the ( $g_{\lambda }$ -invariant) unstable manifold $W^u\left(p_{\lambda }\right)$ has as its linearization the line $W_{\text{lin}}^u\left(p_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}=W_{\text{lin}}^c\left(p_3\right)$, which is not $g_{\lambda }$ -invariant. In this parameter range, the attractor is a 2-cycle on $W^u\left(p_{\lambda }\right)$. The flip bifurcation and period-doubling cascade converging to a Milnor attractor is determined by the restriction of $g_{\lambda }$ to the continuation of $W^u\left(p_{\lambda }\right)$ as a 1-dimensional $g_{\lambda }$ -invariant manifold, which we now analyze.

The determination of the center and unstable manifolds is more conveniently handled by translating the coordinates so that the fixed point is always at the origin as described in (9). In particular, we define the $\lambda$ -dependent translation $h_{\lambda }\mathrm{<mo>\; :\; </mo>}{ℝ}^3\to {ℝ}^3$ by $h_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)=\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}=\left(x-x_{\lambda }\mathrm{,}y-y_{\lambda }\mathrm{,}z-z_{\lambda }\right)$ and the corresponding translated map as

$\begin{array}{c}{\hat{g}}_{\lambda }\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}=h_{\lambda }\circ g_{\lambda }\circ h_{\lambda }^{-1}\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\\ =(\lambda \xi \left\{\left[1-\left(b_1-2\right)z_{\lambda }\right]-\xi -b_1\zeta \right\}-\lambda b_1{x}_{\lambda }\zeta \mathrm{,}\\ \lambda \xi \left\{\left[1-\left(b_2-2\right)z_{\lambda }\right]-\eta -b_2\zeta \right\}-\lambda b_2{y}_{\lambda }\zeta \mathrm{,}\\ \lambda \zeta \left[\left(1-2z_{\lambda }\right)-\zeta \right])\mathrm{,}\end{array}$ (18)

where the coexistence point for the original map in the translated coordinates now corresponds to $\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)=\left(\mathrm{0,0,0}\right)\mathrm{<mo>\; :\; </mo>}=0$ for all $\lambda \in \left(\mathrm{2,3.7}\right]$. We assume that the invariant 1-dimensional invariant manifold that begins as (super) attracting for $2<\lambda <3$, becomes the center manifold $W_{\lambda }^c\left(0\right)$ at $\lambda =3$ and is an unstable manifold $W_{\lambda }^u\left(0\right)$ for $3<\lambda \le 3.1$ can be continued as $W_{\lambda }^{\ast }\left(0\right)$ corresponding to the eigenvalue of largest magnitude for $3.1<\lambda \le 3.7$, is analytic in $\zeta$ and expressible as

$W_{\lambda }^{\ast }\left(0\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}\xi ={\phi }_{\lambda }\left(\zeta \right)\mathrm{,}\eta ={\psi }_{\lambda }\left(\zeta \right)\right\}\mathrm{,}$ (19)

where

${\phi }_{\lambda }\left(\zeta \right)={\displaystyle {\sum }_{m=1}^{\infty }}{\alpha }_m\left(\lambda \right){\zeta }^m\mathrm{,}{\psi }_{\lambda }\left(\zeta \right)={\displaystyle {\sum }_{m=1}^{\infty }}{\beta }_m\left(\lambda \right){\zeta }^m\mathrm{,}$ (20)

with coefficients that are analytic functions of the parameter, converge on a $\zeta$ -interval containing the origin and large enough so that the restriction ${\hat{f}}_{\lambda }\mathrm{<mo>\; :\; </mo>}={{\hat{g}}_{\lambda }|}_{W^{\ast }\left(0\right)}$ exhibits a flip bifurcation at $\lambda =3$ and a period-doubling cascade to converging to a Milnor attractor ${\mathcal{A}}_{\infty }$.

For invariance, we must have

$\begin{array}{l}\left[\lambda -\left(b_1-2\right)\left(\lambda -1\right)\right]{\phi }_{\lambda }\left(\zeta \right)-\lambda {\phi }_{\lambda }\left(\zeta \right)\left({\phi }_{\lambda }\left(\zeta \right)-b_1\zeta \right)-b_1\left(1-b_1\right)\left(\lambda -1\right)\zeta \\ ={\phi }_{\lambda }\left(\zeta \left[\left(2-\lambda \right)-\lambda \zeta \right]\right)\end{array}$ (21)

and

$\begin{array}{l}\left[\lambda -\left(b_2-2\right)\left(\lambda -1\right)\right]{\psi }_{\lambda }\left(\zeta \right)-\lambda {\psi }_{\lambda }\left(\zeta \right)\left({\psi }_{\lambda }\left(\zeta \right)-b_2\zeta \right)-b_2\left(1-b_2\right)\left(\lambda -1\right)\zeta \\ ={\psi }_{\lambda }\left(\zeta \left[\left(2-\lambda \right)-\lambda \zeta \right]\right)\mathrm{,}\end{array}$ (22)

with $0.05\le b_1,b_2<0.08$, which we want to solve for $3\le \lambda \le {\lambda }_{\mathrm{<mo>\; \#\; </mo>}}$. After finding $W_{\lambda }^{\ast }\left(0\right)$, we wish to determine a parameter interval $\left[\mathrm{3,}{\lambda }_{\mathrm{<mo>\; \#\; </mo>}}\right]$ containing in its interior a ${\lambda }_{\infty }$ such that ${\hat{f}}_{{\lambda }_{\infty }}$ has a Milnor attractor that is the period-doubling cascade limit for ${\hat{f}}_{\lambda }$ as $\lambda \to {\lambda }_{\infty }$. Substituting (20) in (21) and (22) and equating coefficients of like powers of $\zeta$, we obtain the values of the leading order coefficients and the following recurrence relations for the higher order terms ( $m\ge 1$ ):

${\alpha }_1=\frac{b_1\left(1-b_1\right)}{4-b_1}\mathrm{,}$ (23)

$\begin{array}{l}{\alpha }_{m+1}=A_m{\left(\lambda \right)}^{-1}\lambda \{\left({\displaystyle \underset{i+j=m-1}{\sum }}{\alpha }_{i+1}{\alpha }_{j+1}\right)-b_1{\alpha }_m\\ +{\displaystyle \underset{i=1}{\overset{\left[\left(m+1\right)/2\right]}{\sum }}}{\left(-1\right)}^i\left(\begin{array}{c}m+1-i\\ i\end{array}\right){\lambda }^{i-1}{\left(2-\lambda \right)}^{m+1-2i}{\alpha }_{m+1-i}\}\mathrm{,}\end{array}$ (24)

and

${\beta }_1=\frac{b_1\left(1-b_1\right)}{4-b_1}\mathrm{,}$ (25)

$\begin{array}{l}{\beta }_{m+1}=B_m{\left(\lambda \right)}^{-1}\lambda \{\left({\displaystyle \underset{i+j=m-1}{\sum }}{\beta }_{i+1}{\beta }_{j+1}\right)-b_2{\beta }_m\\ +{\displaystyle \underset{i=1}{\overset{\left[\left(m+1\right)/2\right]}{\sum }}}{\left(-1\right)}^i\left(\begin{array}{c}m+1-i\\ i\end{array}\right){\lambda }^{i-1}{\left(2-\lambda \right)}^{m+1-2i}{\beta }_{m+1-i}\}\mathrm{,}\end{array}$ (26)

where $A_m\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left[\lambda -\left(b_1-2\right)\left(\lambda -1\right)-{\left(2-\lambda \right)}^{m+1}\right]$, $B_m\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left[\lambda -\left(b_2-2\right)\left(\lambda -1\right)-{\left(2-\lambda \right)}^{m+1}\right]$ and $\left[\right]$ is the usual greatest integer function. Several terms in the series (20) computed using the above recursion formulas are as follows:

$\begin{array}{l}{\phi }_{\lambda }\left(\zeta \right)=\zeta \left[{\alpha }_1+{\alpha }_2\zeta +{\alpha }_3{\zeta }^2+{\alpha }_4{\zeta }^3+\cdots \right]\mathrm{,}\\ {\psi }_{\lambda }\left(\zeta \right)=\zeta \left[{\beta }_1+{\beta }_2\zeta +{\beta }_3{\zeta }^2+{\beta }_4{\zeta }^3+\cdots \right]\mathrm{,}\end{array}$ (27)

where

$\begin{array}{l}{\alpha }_2=\frac{2\lambda b_1\left(1-b_1\right)\left(2-b_1\right){\left(4-b_1\right)}^{-2}}{\left[\left(4-b_1\right)+\left(b_1-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^2\right]},\\ {\beta }_2=\frac{2\lambda b_2\left(1-b_2\right)\left(2-b_2\right){\left(4-b_2\right)}^{-2}}{\left[\left(4-b_2\right)+\left(b_2-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^2\right]},\end{array}$

$\begin{array}{l}{\alpha }_3=\frac{\lambda {\alpha }_2\left[2{\alpha }_1-b_1-2\left(2-\lambda \right)\right]}{\left[\left(4-b_1\right)+\left(b_1-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^3\right]},\\ {\beta }_3=\frac{\lambda {\beta }_2\left[2{\beta }_1-b_2-2\left(2-\lambda \right)\right]}{\left[\left(4-b_2\right)+\left(b_2-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^3\right]},\end{array}$ (28)

$\begin{array}{l}{\alpha }_4=\frac{\lambda \left\{\left[2{\alpha }_1-b_1-3{\left(2-\lambda \right)}^2\right]{\alpha }_3+\left({\alpha }_2+\lambda \right){\alpha }_2\right\}}{\left[\left(4-b_1\right)+\left(b_1-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^4\right]},\\ {\beta }_4=\frac{\lambda \left\{\left[2{\beta }_1-b_2-3{\left(2-\lambda \right)}^2\right]{\beta }_3+\left({\beta }_2+\lambda \right){\beta }_2\right\}}{\left[\left(4-b_2\right)+\left(b_2-3\right)\left(2-\lambda \right)-{\left(2-\lambda \right)}^4\right]}.\end{array}$

Now, it a straightforward matter to prove, for example using the method of majorants, that the power series (27) converge for $\left|\zeta \right|<1$, so the following result is an immediate corollary of Theorem 2 and the fact that the bifurcation behavior is essentially completely determined by the z-coordinate of the map (15).

Lemma 4. The discrete dynamical system generated by the iterates of $g_{\lambda }$ defined by (15) has the following properties for $0.05\le b_1,b_2<0.08$ and $2<\lambda \le 3.7$.

(a) The coexistence fixed point $p_{\lambda }$ given by (16) is a sink for $2<\lambda <3$, has a center manifold $W^c\left(p_3\right)$ at $\lambda =3$ (associated with an eigenvalue $\mu =-1$ of ${g^{\prime }}_{\lambda }\left(p_{\lambda }\right)$ ) along which (15) has a flip bifurcation and there is a stable manifold $W^s\left(p_3\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_3\right)\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$.

(b) For $3<\lambda \le 3.1$, $p_{\lambda }$ is hyperbolic, with stable manifold $W^s\left(p_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}{z}_{\lambda }\right)\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\right\}\cap \stackrel{\circ }{\mathfrak{X}}$ and unstable manifold $W^u\left(p_{\lambda }\right)$ along which there is an 2-cycle attractor.

(c) $W^u\left(p_{\lambda }\right)$ possesses a smooth invariant continuation $W^{\ast }\left(p_{\lambda }\right)$ for $3.1<\lambda \le 3.7$ and $f_{\lambda }\mathrm{<mo>\; :\; </mo>}={g_{\lambda }|}_{W^{\ast }\left(p_{\lambda }\right)}\mathrm{<mo>\; :\; </mo>}{W}^{\ast }\left(p_{\lambda }\right)\to W^{\ast }\left(p_{\lambda }\right)$ has a period-doubling cascade sequence ${\lambda }_n\uparrow {\lambda }_{\infty }\simeq 3.57$ such that $f_{{\lambda }_{\infty }}$ has a Milnor attractor.

It follows from continuity and smoothness considerations that the qualitative behavior characterized in Lemma 4 should persist for sufficiently small additional (nonzero) interaction coefficients in (12). In order to lend some precision to this observation, a decidedly non-optimal example of this is formulated and a proof is sketched in what follows. In addition, the bifurcation behavior of the this system is illustrated by the simulations in Figure 2.

Theorem 5. Suppose the discrete dynamical system generated by the map (12), which is represented in translated form taking the varying coexistence fixed point into the origin by the map (18) when the interaction coefficients $a_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c$

Figure 2. Flip bifurcation for (12): (a) $\lambda =2.8$ ; (b) $\lambda =3.3$ ; (c) $\lambda =3.56$ ; (d) $\lambda =3.7$ with $a_1=0.03$, $a_2=0.01$, $b_1=0.05$, $b_2=0.06$, $c=0.01$ and $d=0.015$.

and d differ from zero, satisfies the following properties:

(i) $0\le a_1,a_2,b_1,b_2,c,d<0.1$.

(ii) $0.05\le b_1\mathrm{,}{b}_2\le 0.08$.

(iii) $0.05\le a_1+b_1\mathrm{,}{a}_2+b_2\le 0.08$.

(iv) $c\mathrm{,}d\mathrm{<mo>\; <\; </mo>}\frac{1}{2}max\left\{a_1\mathrm{,}{a}_2\right\}$.

Then the dynamical system generated by the map (12) has a flip bifurcation at $p_{\lambda }$ for a value of the bifurcation parameter $\lambda ={\lambda }_{\ast }$ near 3 along what is initially a center and then an unstable manifold followed by a smooth invariant extension $W^{\ast }\left(p_{\lambda }\right)$ containing a period-doubling cascade to a (Milnor) attractor ${\mathcal{A}}_{\infty }$ as $\lambda$ increases to a value of about 3.6.

Proof Sketch. We describe the foundational elements of the proof and leave the verification of some of the details which tends to be rather lengthy, albeit routine to the reader. First, it follows from (13) that the unique coexistence fixed point (depending on $a_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c$ and d as well as the bifurcation parameter $\lambda$ ) is

$p_{\lambda }=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)=\sigma \left(\lambda \right)\left(K_1\mathrm{,}{K}_2\mathrm{,}{K}_3\right)\mathrm{,}$ (29)

where

$\begin{array}{l}{K}_3\mathrm{<mo>\; :\; </mo>}=\frac{1-a_1{a}_2-c\left(1-a_1\right)-d\left(1-a_2\right)}{1-a_1{a}_2-c\left(b_1-a_1{b}_2\right)-d\left(b_2-a_2{b}_1\right)}\mathrm{,}\\ K_2\mathrm{<mo>\; :\; </mo>}=\frac{\left[1-a_2\left(1-b_1\right)-b_2\right]}{1-a_1{a}_2}{K}_3\mathrm{,}\\ K_1\mathrm{<mo>\; :\; </mo>}=1-a_1{K}_2-b_1{K}_3\mathrm{,}\end{array}$ (30)

and it is easy to verify that all of the K’s are positive owing to assumptions (i)-(iv). Next, we compute the eigenvalues of the derivative of the map at the fixed points given by the solutions of the characteristic equation

$\Phi \left(\mu \right)\mathrm{<mo>\; :\; </mo>}={\mu }^3+P{\mu }^2+Q\mu +R=\mathrm{0,}$ (31)

which we shall compare with the characteristic equation for (15); namely,

${\Phi }_0\left(\mu \right)\mathrm{<mo>\; :\; </mo>}={\mu }^3+P_0{\mu }^2+Q_0\mu +R_0=\mathrm{0,}$ (32)

where

$P_0=P_0\left(\lambda \mathrm{<mo>\; ;\; </mo>}{a}_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\right)\mathrm{<mo>\; :\; </mo>}=-b_1-b_2-\left(3-b_1-b_2\right)\left(2-\lambda \right)$

$\begin{array}{l}{Q}_0=Q_0\left(\lambda \mathrm{<mo>\; ;\; </mo>}{a}_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\right)\\ \mathrm{<mo>\; :\; </mo>}=b_1{b}_2+\left(b_1+b_2+b_1{b}_2^{\ast }+b_1^{\ast }{b}_2\right)\left(2-\lambda \right)+\left(b_1^{\ast }+b_2^{\ast }+b_1^{\ast }{b}_2^{\ast }\right){\left(2-\lambda \right)}^2\end{array}$

$\begin{array}{l}{R}_0=R_0\left(\lambda \mathrm{<mo>\; ;\; </mo>}{a}_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\right)\\ \mathrm{<mo>\; :\; </mo>}=-\left(2-\lambda \right)\left[b_1+b_2+\left(b_1{b}_2^{\ast }+b_1^{\ast }{b}_2\right)\left(2-\lambda \right)+b_1^{\ast }{b}_2^{\ast }{\left(2-\lambda \right)}^2\right]\end{array}$ (33)

and $b_i^{\ast }\mathrm{<mo>\; :\; </mo>}=1-b_i$, $i=\mathrm{1,2}$. We shall compare these with the coefficients of (31), which are readily found (in terms of the fixed point (29)) to be

$P=P\left(\lambda ;a_1,a_2,b_1,b_2,c,d\right):=P_0+{\Delta }_P=\lambda \left(x_{\lambda }+y_{\lambda }+z_{\lambda }\right)-3,$

$\begin{array}{l}Q=Q\left(\lambda \mathrm{<mo>\; ;\; </mo>}{a}_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\right)\mathrm{<mo>\; :\; </mo>}=Q_0+{\Delta }_Q\\ =3-2\lambda \left(x_{\lambda }+y_{\lambda }+z_{\lambda }\right)+{\lambda }^2\left[\left(1-a_1{a}_2\right)x_{\lambda }{y}_{\lambda }+\left(1-b_{1}c\right)x_{\lambda }{z}_{\lambda }+\left(1-b_{2}d\right)y_{\lambda }{z}_{\lambda }\right]\end{array}$ (34)

$\begin{array}{l}R=R\left(\lambda \mathrm{<mo>\; ;\; </mo>}{a}_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c\mathrm{,}d\right)\mathrm{<mo>\; :\; </mo>}=R_0+{\Delta }_R\\ =-1+\lambda \left(x_{\lambda }+y_{\lambda }+z_{\lambda }\right)-{\lambda }^2[\left(1-a_1{a}_2\right)x_{\lambda }{y}_{\lambda }+\left(1-b_{1}c\right)x_{\lambda }{z}_{\lambda }\\ +\left(1-b_{2}d\right)y_{\lambda }{z}_{\lambda }]+{\lambda }^3\left[1+a_1\left(b_{2}c-a_2\right)+b_1\left(a_{2}d-c\right)-b_{2}d\right]x_{\lambda }{y}_{\lambda }{z}_{\lambda }\mathrm{.}\end{array}$

It is a simple, but rather tedious, matter to use the formulas above to verify that the assumptions (i)-(iv) imply that the following properties obtain:

(p1) $\left|P-P_0\right|,\left|Q-Q_0\right|,\left|R-R_0\right|<0.06$ for all $\lambda \in \left(\mathrm{2.6,3.6}\right]$.

(p2) Two of the eigenvalues (solutions of (31)) of ${G^{\prime }}_{\lambda }\left(p_{\lambda }\right)$, with corresponding eigenvectors nearly parallel to the xy-plane, satisfy $\left|{\mu }_1\right|\mathrm{,}\left|{\mu }_2\right|<1$ for all $\lambda \in \left(\mathrm{2.9,3.1}\right)$.

(p3) The remaining eigenvalue ${\mu }_3$ of ${G^{\prime }}_{\lambda }\left(p_{\lambda }\right)$ is negative for $\lambda \in \left(\mathrm{2.1,3.7}\right)$ and is decreasing as a function of $\lambda$ such that ${\mu }_3=-1$ at $\lambda ={\lambda }_1$, where $\left|{\lambda }_1-3\right|<0.06$, which corresponds to a flip bifurcation.

Using (p1)-(p3) in conjunction with the translation of variables to convert the fixed point to the origin for all $\lambda \in \left(\mathrm{2.6,4}\right)$ as was done for the simpler map (15), it can again be shown that there is an invariant 1-dimensional (real) analytic manifold that is a center manifold for $\lambda ={\lambda }_1\underset{\_}{\approx }3$ and is the unstable manifold $W^u\left(p_{\lambda }\right)$ for $3<\lambda <3.09$. Consequently, the restriction of the map to $W^u\left(p_{\lambda }\right)=\left\{\left(x_{\lambda }+{\tilde{\phi }}_{\lambda }\left(z-z_{\lambda }\right)\mathrm{,}{y}_{\lambda }+{\tilde{\psi }}_{\lambda }\left(z-z_{\lambda }\right)\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\in ℝ\right\}$ denoted as $\mathfrak{g}\mathrm{<mo>\; :\; </mo>}{W}^u\left(p_{\lambda }\right)\to W^u\left(p_{\lambda }\right)$, where we have used a form of parametrization analogous to that in our analysis of (15), can be shown to have the same basic unimodular form as that of (15). Hence, it follows that this is an analog of $W^{\ast }\left(p_{\lambda }\right)$ on which the map has a period-doubling (parameter) cascade ${\lambda }_n\uparrow {\lambda }_{\infty }$, with $\left|{\lambda }_{\infty }-3.57\right|<0.08$, to a Milnor attractor. Thus, the sketch of our proof is complete. $\square$

4. Higher Dimensional Neimark-Sacker Type Bifurcations

The type of higher dimensional analog of the Neimark-Sacker bifurcation that we have in mind is one of codimension 1 for a discrete dynamical system in ${ℝ}^m$, $m\ge 3$, with a fixed point that bifurcates from a sink to a source, giving birth to an invariant homeomorph of an $\left(m-1\right)$ -sphere (which is typically as smooth as the system except on a lower dimensional subset) that grows in size as a particular parameter increases. Although we concentrate mainly on $m=3$, the following is a nice smooth example of the type of bifurcation we interested in for $m\ge 2$ : Consider the smooth (actually real analytic) map $T_{\lambda }\mathrm{<mo>\; :\; </mo>}{ℝ}^m\to {ℝ}^m$ defined as

$T_{\lambda }\left(x\right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{ta{n}^{-1}\left(\lambda \left|x\right|\right)}{\left|x\right|}\right)x\mathrm{,}$ (35)

where $|\cdot |$ denotes the Euclidean norm on ${ℝ}^m$ and the parameter $\lambda \in \left(\mathrm{0,}\infty \right)$. The origin is a fixed point for all $\lambda$, which is a sink for $0<\lambda <1$ and a source when $\lambda >1$. We note that the origin is the only fixed point for $0<\lambda \le 1$, but for $\lambda >1$ it is easy to see that the fixed point set of the map comprises all points on the $\left(m-1\right)$ -sphere $S_{\lambda }\mathrm{<mo>\; :\; </mo>}=\left\{x\in {ℝ}^m\mathrm{<mo>\; :\; </mo>}tan\left|x\right|=\lambda \left|x\right|\mathrm{<mo>\; :\; </mo>}=\lambda r\left(\lambda \right)\right\}$ in addition to the origin. Moreover, it is rather easy to show that $r=r\left(\lambda \right)$ is a smooth increasing function for $1\le \lambda <\infty$ such that $r\left(1\right)=0$ and $r\left(\lambda \right)\to \pi /2$ as $\lambda \to \infty$. Furthermore, $S_{\lambda }$ is an attractor with basin of attraction $B\left(S_{\lambda }\right)={ℝ}^m\backslash \left\{0\right\}$.

In the sequel, we shall almost exclusively consider such bifurcations for $m=3$, where a particularly simple example of the piecewise smooth variety is

$F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x\right)\mathrm{,}\lambda y\left(1-y\right)\mathrm{,}\lambda z\left(1-z\right)\right)\mathrm{,}$ (36)

for the fixed point $P\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)=\left(\left(\lambda -1\right)/\lambda \right)\left(\mathrm{1,1,1}\right)$ in a small neighborhood of $\lambda =3$. Given that each of the coordinate functions is just a logistic map, it is easy to see that $P\left(\lambda \right)$ is a sink for $2<\lambda <3$ and a source for $3<\lambda <4$, so that $\lambda =3$ is a bifurcation value for the parameter. Moreover, the system (36) has the four attracting 2-cycles for $3<\lambda <3.4$

$\begin{array}{l}\left\{\left(x_1\mathrm{,}{y}_1\mathrm{,}{z}_1\right),\left(x_2\mathrm{,}{y}_2\mathrm{,}{z}_2\right)\right\},\left\{\left(x_2\mathrm{,}{y}_1\mathrm{,}{z}_1\right),\left(x_1\mathrm{,}{y}_2\mathrm{,}{z}_2\right)\right\},\\ \left\{\left(x_1\mathrm{,}{y}_2\mathrm{,}{z}_1\right),\left(x_2\mathrm{,}{y}_1\mathrm{,}{z}_2\right)\right\},\left\{\left(x_1\mathrm{,}{y}_1\mathrm{,}{z}_2\right),\left(x_2\mathrm{,}{y}_2\mathrm{,}{z}_1\right)\right\},\end{array}$ (37)

in addition to nine 2-cycles of saddle type obtained by replacing one or two coordinate pairs in these cycles with fixed point coordinates; such as, $\left\{\left(x_{\lambda }\mathrm{,}{y}_1\mathrm{,}{z}_1\right)\mathrm{,}\left(x_{\lambda }\mathrm{,}{y}_2\mathrm{,}{z}_2\right)\right\}$ and $\left\{\left(x_1\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)\mathrm{,}\left(x_2\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)\right\}$. Here

$\begin{array}{l}{x}_1=y_1=z_1=\frac{1}{2\lambda }\left[\left(\lambda +1\right)-\sqrt{{\lambda }^2-2\lambda -3}\right],\\ x_2=y_2=z_2=\frac{1}{2\lambda }\left[\left(\lambda +1\right)+\sqrt{{\lambda }^2-2\lambda -3}\right],\end{array}$ (38)

which define the vertices (37) of the boundary $\partial R\left(\lambda \right)$ of the rectangular solid $R\left(\lambda \right)$ containing $P\left(\lambda \right)$ when $\lambda -3$ is small and positive. For these parameter values, it is not difficult to verify that $\partial R\left(\lambda \right)$, which is homeomorphic and piecewise-linearly diffeomorphic to the 2-sphere ${\mathbb{S}}^2$ and is an F-invariant, locally attracting set for the dynamical system having a diameter that increases with $\lambda$, as indicated by the simulation in Fig. 4.1. This behavior is somewhat like a piecewise-linear analog of Neimark-Sacker bifurcations in ${ℝ}^2$ (see [19] [20] ); one in fact that has several natural variations. For example, if we consider the following minor modification of (36):

$F_{\epsilon }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x\right)\mathrm{,}\left(\lambda -\epsilon \right)y\left(1-y\right)\mathrm{,}\left(\lambda -\epsilon \right)z\left(1-z\right)\right)\mathrm{,}$ (39)

where $0<\epsilon <1$, it is not difficult to verify that this system first has a flip bifurcation at $\lambda =3$, followed by the fixed point of interest changing from a saddle to a source (at $\lambda =3+\epsilon$ ) and generating a locally attracting, F-invariant piecewise-linear isomorph of ${\mathbb{S}}^2$ like that of (36), only elongated in the x-direction. It is also rather easy to envisage even higher dimensional generalization of the bifurcations described above for dynamical system far more varied than discrete Lotka-Volterra types. However, here we shall confine our attention to such bifurcations in 3-dimensional Lotka-Volterra systems, reserving a more extensive, general treatment for a later investigation.

Definition of Neimark-Sacker Analogs and Another L-V Example

We first define the analogs of Neimark-Sacker bifurcations for discrete 3-dimensional dynamical systems briefly described above and then state an existence theorem for certain Lotka-Volterra systems. For this, we consider a parameter-dependent $C^1$ map of the form

$\Theta \mathrm{<mo>\; :\; </mo>}U\times \left(a\mathrm{,}b\right)\to U\subset {ℝ}^3\mathrm{,}$

$\Theta \left(x\mathrm{,}y\mathrm{,}z\mathrm{<mo>\; ;\; </mo>}\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left({\theta }_1\left(x\mathrm{,}y\mathrm{,}z\mathrm{<mo>\; ;\; </mo>}\lambda \right)\mathrm{,}{\theta }_2\left(x\mathrm{,}y\mathrm{,}z\mathrm{<mo>\; ;\; </mo>}\lambda \right)\mathrm{,}{\theta }_3\left(x\mathrm{,}y\mathrm{,}z\mathrm{<mo>\; ;\; </mo>}\lambda \right)\right)\mathrm{,}$ (40)

where U is an open subset of ${ℝ}^3$ and the parameter $\lambda$ ranges over the interval $\left(a\mathrm{,}b\right)\subset ℝ$.

Definition NST. Suppose that (40) has an isolated hyperbolic fixed point $P\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(x_{\ast }\left(\lambda \right)\mathrm{,}{y}_{\ast }\left(\lambda \right)\mathrm{,}{z}_{\ast }\left(\lambda \right)\right)$ for each $\lambda \in \left(a\mathrm{,}b\right)$. Then (40) has an NST bifurcation at $\lambda ={\lambda }_0\in \left(a\mathrm{,}b\right)$ if the following properties obtain:

(i) There is a neighborhood $\left({\lambda }_0-\epsilon \mathrm{,}{\lambda }_0+\epsilon \right)\subset \left(a\mathrm{,}b\right)$ such that at least one eigenvalue $\mu$ of ${{\Theta }^{\prime }}_{\lambda }$ at $P\left(\lambda \right)$ satisfies $\left|\mu \right|<1$ whenever ${\lambda }_0-\epsilon <\lambda <{\lambda }_0$ and all eigenvalues $\mu$ of ${{\Theta }^{\prime }}_{\lambda }$ at $P\left(\lambda \right)$ satisfy $\left|\mu \right|>1$ for ${\lambda }_0<\lambda <{\lambda }_0+\epsilon$.

(ii) There is a locally attracting, ${\Theta }_{\lambda }$ -invariant homeomorph $\partial R\left(\lambda \right)$ of the boundary of a rectangular box $R\left(\lambda \right)$ or 2-sphere ${\mathbb{S}}^2$ enclosing $P\left(\lambda \right)$ for ${\lambda }_0<\lambda <{\lambda }_0+\epsilon$ such that the Lebesgue measure of $R\left(\lambda \right)$ increases with $\lambda$ from zero at $\lambda ={\lambda }_0$.

Observe that the key elements, such as the 2-cycles and invariant 2-cycle pairs of planes, of the NST bifurcation for the simple map (36) are all hyperbolic or normally hyperbolic and therefore locally structurally stable. Consequently, one expects that small changes in the interaction coefficients (which are zero) should still produce NST bifurcations, which is confirmed by our next result.

Theorem 6. If the Lotka-Volterra map $F\mathrm{<mo>\; :\; </mo>}{ℝ}^3\to {ℝ}^3$ defined as

$F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x-a_{1}y-b_{1}z\right)\mathrm{,}\lambda y\left(1-a_{2}x-y-b_{2}z\right)\mathrm{,}\lambda z\left(1-cx-dy-z\right)\right)\mathrm{,}$ (41)

is such that $a_1\mathrm{,}{a}_2\mathrm{,}{b}_1\mathrm{,}{b}_2\mathrm{,}c$ and d are sufficiently small and nonnegative, it has an NST bifurcation in the interior of the first octant $O_{+}:=\left\{\left(x,y,z\right):x,y,z>0\right\}$ such that the corresponding invariant sphere homeomorph $\partial R\left(\lambda \right)$ is the boundary of a curvilinear rectangular box with edges and diagonal vertices comprising a locally attracting set and four attracting 2-cycles, respectively.

Proof. In the interest of simplicity, we shall give the proof only for the case where $a_1=a_2:=a,b_1=b_2:=b$ and $c,d=0$. There is really no loss of generality in this since the general argument is essentially the same as this simpler one, only the calculations are considerably more complicated, albeit routine. Therefore, we assume that

$F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x-ay-bz\right)\mathrm{,}\lambda y\left(1-ax-y-bz\right)\mathrm{,}\lambda z\left(1-z\right)\right)\mathrm{,}$ (42)

and $0\le a,b<0.1$. It is easy to show that the coexistence fixed point of (42) is $p\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{\lambda -1}{\lambda }\right)\left(\frac{1-b}{1+a}\mathrm{,}\frac{1-b}{1+a}\mathrm{,1}\right)\mathrm{<mo>\; :\; </mo>}=\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)$ for $\lambda >1$ and that the eigenvalues of $F^{\prime }\left(p\left(\lambda \right)\right)$ are

${\mu }_1={\mu }_2=\frac{\left(2+a-b\right)-\left(1-b\right)\lambda }{1+a},{\mu }_3=2-\lambda ,$ (43)

so it follows that $p\left(\lambda \right)$ is a hyperbolic sink for $1<\lambda <3$ and there is a bifurcation at $\lambda =3$ across which we observe, for all permissible choices of the interaction coefficients, a center manifold at $\lambda =3$, $W^c\left(p\left(\lambda \right)\right)$, followed by an unstable manifold, $W^u\left(p\left(\lambda \right)\right)$, for $\lambda >3$, which is 1-dimensional and transverse to the plane $z=z_{\lambda }$ at the coexistence fixed point. One should also note that the plane $z=z_{\lambda }$ is F-invariant. Our focus is actually going to be on F² and we shall find it convenient to use the translated form of the map described in (9) for which the coexistence fixed point is moved to the origin. The translated version of (42) is

$\begin{array}{l}\hat{F}\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}=\left(F_1\mathrm{,}{F}_2\mathrm{,}{F}_3\right)\\ \mathrm{<mo>\; :\; </mo>}=\left(\xi -\lambda \left(\xi +x_{\lambda }\right)G\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{,}\eta -\lambda \left(\eta +y_{\lambda }\right)G\left(\eta \mathrm{,}\xi \mathrm{,}\zeta \right)\mathrm{,}\zeta -\lambda \left(\zeta +z_{\lambda }\right)\zeta \right)\mathrm{,}\end{array}$ (44)

from which we find that

$\begin{array}{l}{\hat{F}}^2\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}=\left(F_1^2\mathrm{,}{F}_2^2\mathrm{,}{F}_3^2\right)\\ \mathrm{<mo>\; :\; </mo>}=(F_1-\lambda \left(F_1+x_{\lambda }\right)G\left(F_1\mathrm{,}{F}_2\mathrm{,}{F}_3\right)\mathrm{,}{F}_2-\lambda \left(F_2+y_{\lambda }\right)G\left(F_2\mathrm{,}{F}_1\mathrm{,}{F}_3\right)\mathrm{,}\\ F_3\left[1-\lambda \mathrm{<mo\; stretchy="false">\; (\; </mo>}{F}_3+z_{\lambda }\mathrm{<mo\; stretchy="false">\; )\; </mo>}\right])\mathrm{,}\end{array}$ (45)

where $G=G\left(\xi ,\eta ,\zeta \right):=\xi +a\eta +b\zeta$ and

$F_1^2:=F_1^2\left(\xi ,\eta ,\zeta \right):=A\left(\lambda \right)\xi -2aB\left(\lambda \right)\eta -bC\left(\lambda \right)\zeta -H_1\left(\xi ,\eta ,\zeta ;\lambda \right),$

$F_2^2:=F_2^2\left(\xi ,\eta ,\zeta \right):=-2aB\left(\lambda \right)\xi +A\left(\lambda \right)\eta -bC\left(\lambda \right)\zeta -H_2\left(\xi ,\eta ,\zeta ;\lambda \right),$ (46)

$F_3^2:=F_3^2\left(\xi ,\eta ,\zeta \right):=\zeta \left[{\left(2-\lambda \right)}^2-\lambda \left(\lambda -1\right)\left(4-\lambda \right)\zeta -2\lambda \left(\lambda -2\right){\zeta }^2-{\lambda }^3{\zeta }^3\right],$

where

$A\left(\lambda \right):=\left\{1-\left(\frac{1-b}{1+a}\right)\left(\lambda -1\right)\left[2-\left(\frac{1-b}{1+a}\right)\left(\lambda -1\right)\right]\right\},$

$B\left(\lambda \right)=\left(\frac{1-b}{1+a}\right)\left(\lambda -1\right)\left[1-\left(\frac{1-b}{1+a}\right)\left(\lambda -1\right)\right]\mathrm{,}$ (47)

$C\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{1-b}{1+a}\right)\left(\lambda -1\right)\left\{-2+\left(\lambda -1\right)\left[1+\left(1-b\right)\left(\lambda -1\right)\right]\right\}$

and the higher order terms in (46) are

$\begin{array}{l}{H}_1\mathrm{<mo>\; :\; </mo>}=\lambda \{\left(1-\lambda x_{\lambda }\right)\xi G\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)+a\lambda x_{\lambda }\eta G\left(\eta \mathrm{,}\xi \mathrm{,}\zeta \right)\\ +b\lambda x_{\lambda }{\zeta }^2+{\left(F_1\right)}^2+a{F}_1{F}_2+b{F}_1{F}_3\}\mathrm{,}\end{array}$

$\begin{array}{l}{H}_2\mathrm{<mo>\; :\; </mo>}=\lambda \{a\lambda y_{\lambda }\xi G\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)+\left(1-\lambda y_{\lambda }\right)\eta G\left(\eta \mathrm{,}\xi \mathrm{,}\zeta \right)\\ +b\lambda y_{\lambda }{\zeta }^2+{\left(F_2\right)}^2+a{F}_1{F}_2+b{F}_2{F}_3\}\mathrm{.}\end{array}$ (48)

As a first step in completing the proof, we note that it follows from our previous analysis that there exists a positive $\epsilon$ such that when $0<a,b<\epsilon$, there are parameter values ${\lambda }_1<3<{\lambda }_2$ for which the following property obtains: The coexistence fixed point 0 of (44) is a sink at $\lambda ={\lambda }_1$ and a source when $\lambda ={\lambda }_2$. Moreover, there is a unique ${\lambda }_{\ast }\in \left({\lambda }_1\mathrm{,}{\lambda }_2\right)$ across which 0 transitions from a sink (as in the case when all the interaction coefficients are zero) or a saddle point to a source. What we have also shown in the zero interaction (uncoupled) case, which we denote as ${\hat{F}}_0$, is that there are unique pairs of invariant, locally attracting plane cycles for $\lambda >{\lambda }_{\ast }$ of the form

${\mathcal{P}}_{\xi }\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left\{P_{\xi }\left(\lambda \right)\mathrm{,}{Q}_{\xi }\left(\lambda \right)\right\}\mathrm{<mo>\; =\; </mo>}\left\{\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}\xi \mathrm{<mo>\; =\; </mo>}{A}_{+}\left(\lambda \right)\right\}\mathrm{,}\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}\xi =A_{-}\left(\lambda \right)\right\}\right\}\mathrm{,}$

${\mathcal{P}}_{\eta }\left(\lambda \right):=\left\{P_{\eta }\left(\lambda \right),Q_{\eta }\left(\lambda \right)\right\}=\left\{\left\{\left(\xi ,\eta ,\zeta \right):\eta =B_{+}\left(\lambda \right)\right\},\left\{\left(\xi ,\eta ,\zeta \right):\eta =B_{-}\left(\lambda \right)\right\}\right\},$ (49)

${\mathcal{P}}_{\zeta }\left(\lambda \right):=\left\{P_{\zeta }\left(\lambda \right),Q_{\zeta }\left(\lambda \right)\right\}=\left\{\left\{\left(\xi ,\eta ,\zeta \right):\zeta =C_{+}\left(\lambda \right)\right\},\left\{\left(\xi ,\eta ,\zeta \right):\zeta =C_{-}\left(\lambda \right)\right\}\right\},$

$B_{-}\left(\lambda \right)$ where $A_{-}\left(\lambda \right)<0<A_{+}\left(\lambda \right)$, $B_{-}\left(\lambda \right)<0<B_{+}\left(\lambda \right)$, $C_{-}\left(\lambda \right)<0<C_{+}\left(\lambda \right)$, and $A_{+}\left(\lambda \right)-A_{-}\left(\lambda \right)$, $B_{+}\left(\lambda \right)-B_{-}\left(\lambda \right)$ and $C_{+}\left(\lambda \right)-C_{-}\left(\lambda \right)$ are increasing functions of $\lambda$ that converge to zero as $\lambda \downarrow {\lambda }_{\ast }$. Furthermore, ${\hat{F}}_0\left(P_{\xi }\left(\lambda \right)\right)=Q_{\xi }\left(\lambda \right)$, ${\hat{F}}_0\left(Q_{\xi }\left(\lambda \right)\right)=P_{\xi }\left(\lambda \right)$, ${\hat{F}}_0\left(P_{\eta }\left(\lambda \right)\right)=Q_{\eta }\left(\lambda \right)$, ${\hat{F}}_0\left(Q_{\eta }\left(\lambda \right)\right)=P_{\eta }\left(\lambda \right)$ and ${\hat{F}}_0\left(P_{\zeta }\left(\lambda \right)\right)=Q_{\zeta }\left(\lambda \right)$, ${\hat{F}}_0\left(Q_{\zeta }\left(\lambda \right)\right)=P_{\zeta }\left(\lambda \right)$. With these planes, we have associated lamina defined as

${\mathcal{L}}_{\xi }\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}{A}_{-}\left(\lambda \right)\le \xi \le A_{+}\left(\lambda \right)\right\}\mathrm{,}$

${\mathcal{L}}_{\eta }\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}{B}_{-}\left(\lambda \right)\le \eta \le B_{+}\left(\lambda \right)\right\}\mathrm{,}$ (50)

${\mathcal{L}}_{\zeta }\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}{C}_{-}\left(\lambda \right)\le \zeta \le C_{+}\left(\lambda \right)\right\}\mathrm{,}$

such that for every $\lambda >{\lambda }_{\ast }$, the fixed point 0 of ${\hat{F}}_0$ is in the interior of the invariant rectangular box

$R\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}={\mathcal{L}}_{\xi }\left(\lambda \right)\cap {\mathcal{L}}_{\eta }\left(\lambda \right)\cap {\mathcal{L}}_{\zeta }\left(\lambda \right)\mathrm{,}$ (51)

which is such that every point in a sufficiently small neighborhood, with the exception of 0, is attracted to the four pairs of diagonal 2-cycles comprising the vertices of $\partial R\left(\lambda \right)$. Naturally, one expects that curvilinear analogs of the above results to hold for $\epsilon$ sufficiently small, which what we shall now confirm to complete the proof. We first seek an attracting ${\hat{F}}^2$ -invariant surface of the form

${\mathcal{S}}_{\zeta }^{+}\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}\zeta ={\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right\}\mathrm{,}$

where ${\Phi }_{\lambda }$ is (real) analytic in a neighborhood of $\left(\xi \mathrm{,}\eta \right)=\left(\mathrm{0,0}\right)$ having the power series expansion

${\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)={\displaystyle {\sum }_{m\mathrm{,}n\ge 0}}{c}_{mn}\left(\lambda \right){\xi }^m{\eta }^n\mathrm{,}$ (52)

with $a_{00}\left(\lambda \right)>0$ an increasing function of $\lambda$ for $\lambda >{\lambda }_{\ast }$ that converges to zero as $\lambda \downarrow {\lambda }_{\ast }$. Invariance requires that ${\Phi }_{\lambda }$ satisfy the equation

$F_3^2\left(\xi \mathrm{,}\eta \mathrm{,}{\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right)={\Psi }_{\lambda }\left(F_1^2\left(\xi \mathrm{,}\eta \mathrm{,}{\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right)\mathrm{,}{F}_2^2\left(\xi \mathrm{,}\eta \mathrm{,}{\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right)\right)\mathrm{.}$ (53)

Upon substituting (52) in the above equation, it is straightforward to inductively show that if $\epsilon$ is sufficiently small and $\lambda -{\lambda }_{\ast }\le 0.1$, then

$\left|c_{mn}\right|\le 2^{\left(m+n+1\right)}$

for all integers $m\mathrm{,}n\ge 0$. Hence, the series (52) converges uniformly for ${\xi }^2+{\eta }^2\le 1/5$ over the specified range of the parameter and the sufficiently small value of $\epsilon$. Of course this leads to the companion surface beneath the $\zeta =0$ plane in the 2-cycle given by ${\mathcal{S}}_{\zeta }^{-}\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\hat{F}\left({\mathcal{S}}_{\zeta }^{+}\left(\lambda \right)\right)$. Clearly, this procedure can be repeated for the first two coordinates to obtain analytic ${\hat{F}}^2$ -invariant surfaces

${\mathcal{S}}_{\xi }^{+}\left(\lambda \right):=\left\{\left(\xi ,\eta ,\zeta \right):\xi ={\Theta }_{\lambda }\left(\eta ,\zeta \right)\right\}\text{and}{\mathcal{S}}_{\eta }^{+}\left(\lambda \right):=\left\{\left(\xi ,\eta ,\zeta \right):\eta ={\Phi }_{\lambda }\left(\xi ,\zeta \right)\right\},$

where

${\Theta }_{\lambda }\left(\eta \mathrm{,}\zeta \right)={\displaystyle {\sum }_{m\mathrm{,}n\ge 0}}{a}_{mn}\left(\lambda \right){\eta }^m{\zeta }^n\text{and}{\Phi }_{\lambda }\left(\xi \mathrm{,}\zeta \right)={\displaystyle {\sum }_{m\mathrm{,}n\ge 0}}{b}_{mn}\left(\lambda \right){\xi }^m{\zeta }^n$

converge uniformly for ${\eta }^2+{\zeta }^2\le 1/5$ and ${\xi }^2+{\zeta }^2\le 1/5$, respectively. In addition, one has their respective analytic companion 2-cycle surfaces

${\mathcal{S}}_{\xi }^{-}\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\hat{F}\left({\mathcal{S}}_{\xi }^{+}\left(\lambda \right)\right)\text{and}{\mathcal{S}}_{\eta }^{-}\left(\lambda \right)\mathrm{<mo>\; :\; </mo>\; <mo>\; =\; </mo>}\hat{F}\left({\mathcal{S}}_{\eta }^{+}\left(\lambda \right)\right)\mathrm{.}$

Putting all of this together, we obtain the $\hat{F}$ -invariant curvilinear rectangular box

$\begin{array}{l}\tilde{R}\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\{\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\mathrm{<mo>\; :\; </mo>}{F}_3\left({\Theta }_{\lambda }\left(\eta \mathrm{,}\zeta \right)\mathrm{,}\eta \mathrm{,}\zeta \right)\le \xi \le {\Theta }_{\lambda }\left(\eta \mathrm{,}\zeta \right)\mathrm{,}\\ F_3\left(\xi \mathrm{,}{\Phi }_{\lambda }\left(\xi \mathrm{,}\zeta \right)\mathrm{,}\zeta \right)\le \eta \le {\Phi }_{\lambda }\left(\xi \mathrm{,}\zeta \right)\mathrm{,}\\ \left\{F_3\left(\xi \mathrm{,}\eta \mathrm{,}{\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right)\le \zeta \le {\Psi }_{\lambda }\left(\xi \mathrm{,}\eta \right)\right\}\}\end{array}$

Then, it follows from the definition of the map that all iterates of points in the interior of $\tilde{R}\left(\lambda \right)$, except for the fixed point 0, converge to $\partial \tilde{R}\left(\lambda \right)$, as do all sufficiently close points in the exterior of $\tilde{R}\left(\lambda \right)$. More precisely, the convergence is to the union of the four 2-cycle sinks comprising the diagonal vertex pairs. Thus, the proof is complete.

It should also to be possible to prove the above theorem using the implicit function theorem, but neither this nor the method used in the proof appears to be feasible for higher dimensions and more general parameter ranges. However, it appears some of the iterative image techniques employed in [21] can adapted to extending the results here and we plan to demonstrate this in a forthcoming paper.

The NST bifurcation for (42) in the case where $a=b=0.03$ is illustrated in Figure 3 via orbits starting near the fixed point. Observe how for $\lambda =2.8$ all the iterates tend to the fixed point sink. Note the contrast in behavior for $\lambda =3.3$ where all of the iterates converge to an attracting 2-cycle comprising diagonal vertices of the boundary of the invariant rectangular box.

5. Steady-State Chaotic Strange Attracting Sets

There are numerous types of chaotic strange attracting sets that can occur for the dynamics of the map (1). Naturally, there are the well-known lower dimensional examples embedded in the coordinate lines and planes such as in [4] [5] [8]. More to the point for the purposes here are fully 3-dimensional examples such as those indicated by simulation and shown in Figures 4-6.

In this section we shall focus our attention on a detailed analysis of types of

Figure 3. NST bifurcation for (42) with $a=b=0.03$ shown for (a) $\lambda =2.8$ ; (b) $\lambda =3.2$ via orbits starting at three positions near the fixed point in each case.

Figure 4. Attractor for (1) with ${\lambda }_1=2.2$, ${\lambda }_2=4$, ${\lambda }_3=3.9$ ; ${\gamma }_{12}=0.02$, ${\gamma }_{13}=0.03$, ${\gamma }_{21}=0.07$, ${\gamma }_{23}=0.05$, ${\gamma }_{31}=0.001$, ${\gamma }_{32}=0.01$.

Figure 5. Attractor for (1) with ${\lambda }_1=2.2$, ${\lambda }_2=4$, ${\lambda }_3=3.9$ ; ${\gamma }_{12}=0.2$, ${\gamma }_{13}=0.3$, ${\gamma }_{21}=0.7$, ${\gamma }_{23}=0.5$, ${\gamma }_{31}=0.01$, ${\gamma }_{32}=0.1$.

3-dimensional Lotka-Volterra maps exhibiting steady-state behavior manifesting itself in chaotic strange attracting sets (cf. [22] [23] ).

We consider maps of the form $F=F_{\lambda }\mathrm{<mo>\; :\; </mo>}{\mathfrak{X}}_F\to {\mathfrak{X}}_F$ defined as

$\begin{array}{l}F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(f_1\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_2\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_3\left(x\mathrm{,}y\mathrm{,}z\right)\right)\\ \mathrm{<mo>\; :\; </mo>}=({\lambda }_{1}x\left(1-x-{\gamma }_{12}y-{\gamma }_{13}z\right)\mathrm{,}{\lambda }_{2}y\left(1-{\gamma }_{21}x-y-{\gamma }_{23}z\right)\mathrm{,}{\lambda }_{3}z\left(1-{\gamma }_{31}x-{\gamma }_{32}y-z\right))\mathrm{,}\end{array}$ (54)

$\Lambda$ where ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3$ are nonnegative and represented collectively as $\lambda \mathrm{<mo>\; :\; </mo>}=\left({\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3\right)$, all interaction coefficients are in $\left(\mathrm{0,0.5}\right)$ and

Figure 6. Attractor for (1) with ${\lambda }_1=2.2$, ${\lambda }_2=4$, ${\lambda }_3=3.9$ ; ${\gamma }_{12}=0.04$, ${\gamma }_{13}=0.01$, ${\gamma }_{21}=0.07$, ${\gamma }_{23}=0.01$, ${\gamma }_{31}=0.09$, ${\gamma }_{32}=0.03$.

$\begin{array}{l}{\mathfrak{X}}_F\mathrm{<mo>\; :\; </mo>}=\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>0}\le x\le 1-{\gamma }_{12}y-{\gamma }_{13}z\mathrm{,0}\le y\le 1-{\gamma }_{21}x-{\gamma }_{23}z\mathrm{,}\\ 0\le z\le 1-{\gamma }_{31}x-{\gamma }_{32}y\}\mathrm{.}\end{array}$ (55)

Furthermore, we assume that $\lambda \in \Lambda \left({\mathfrak{X}}_F\right)$ implies that $F_{\lambda }\left({\mathfrak{X}}_F\right)\subset {\mathfrak{X}}_F$, so that the set $\Lambda \left({\mathfrak{X}}_F\right)$ is admissible inasmuch as it contains all parameter values such that $F_{\lambda }$ maps the compact, convex set ${\mathfrak{X}}_F$ into itself. For convenience, in the sequel we denote the attracting sets and attractors (minimal attracting sets) of the map (54) by $\mathcal{A}$ and $\mathfrak{A}$, respectively, possibly with subscripts.

In what follows we shall demonstrate just how complicated the attracting sets and attractors can be and how wildly and dramatically they can change as certain parameters are varied. However, an inkling of the variety and complexity leading to chaotic attractors can readily be obtained by considering the simple uncoupled map $S=S_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left(\lambda x\left(1-x\right)\mathrm{,}\lambda y\left(1-y\right)\mathrm{,}\lambda z\left(1-z\right)\right)$ on $I^3$ as $\lambda$ increases over the interval $\left(\mathrm{0,4}\right]$. In this case it is easy to see that ${\mathfrak{A}}_0=\left\{0=\left(\mathrm{0,0,0}\right)\right\}$ is a global attractor for $0<\lambda \le 1$ and

${\mathfrak{A}}_1=\left\{\sigma \left(\lambda \right)1=\frac{\lambda -1}{\lambda }1\right\}$, where $1\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{1,1,1}\right)$, for ${\lambda }_1\mathrm{<mo>\; :\; </mo>}=\mathrm{1\; <mo>\; <\; </mo>}\lambda \le 3$ is an attractor

with basin of attraction $\mathcal{B}\left({\mathfrak{A}}_1\right)=I^3\backslash \left\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}x\mathrm{,}y\mathrm{,}z=\mathrm{0,1}\right\}$. Then, there is an increasing (period doubling) sequence $\left\{{\lambda }_n\right\}$ with ${\lambda }_2\mathrm{<mo>\; :\; </mo>}=3$ and ${\lambda }_n\uparrow {\lambda }_{\infty }\simeq 3.57$ with attracting set ${\mathcal{A}}_2={\displaystyle {\cup }_{k=1}^4{\mathfrak{A}}_{\mathrm{2,}k}}$ (comprising four disjoint 2-cycle attractors) with basin of attraction $\mathcal{B}\left({\mathcal{A}}_2\right)=I^3\backslash \left(\left\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}x\mathrm{,}y\mathrm{,}z=\mathrm{0,1}\right\}\cup {\mathfrak{A}}_1\right)$ for ${\lambda }_2\mathrm{<mo>\; :\; </mo>}=3<\lambda \le {\lambda }_3$. In general, for ${\lambda }_n<\lambda \le {\lambda }_{n+1}$ there is an attracting set ${\mathcal{A}}_n={\displaystyle {\cup }_{k=1}^{2^{n+1}}{\mathfrak{A}}_{2^n\mathrm{,}k}}$ (comprising $2^{n+1}$ disjoint $2^n$ -cycle attractors) with basin of attraction $\mathcal{B}\left({\mathcal{A}}_n\right)=I^3\backslash \left(\left\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}x\mathrm{,}y\mathrm{,}z=\mathrm{0,1}\right\}\cup {\mathfrak{A}}_1\cup {\mathcal{A}}_2\cup \cdots \cup {\mathcal{A}}_{n-1}\right)$. Moreover, $S_{{\lambda }_{\infty }}$ has a (Milnor) attractor that is the cartesian product ${\mathcal{A}}_{\infty }^3\mathrm{<mo>\; :\; </mo>}={\mathcal{A}}_{\infty }\times {\mathcal{A}}_{\infty }\times {\mathcal{A}}_{\infty }\subset X_{S_{{\lambda }_{\infty }}}$ of the period-doubling cascade limit attractor for each of the coordinate functions and $S_4$ has a dense chaotic (Milnor) attractor $\mathfrak{A}\subset {\mathfrak{X}}_F=I^3$.

5.1. Bifurcation of an Attractor from a Sink to an Attracting Cycle to a Chaotic Strange Attractor

We shall now analyze a L-V map family with a codimension-1 bifurcation sequence beginning with a sink attractor followed by cycle attractors and ending in a true horseshoe type chaotic strange attractor. In particular, for ease of computation we consider the family of maps ${\Phi }_{\lambda }\mathrm{<mo>\; :\; </mo>}\mathfrak{X}\to \mathfrak{X}$ for $\lambda >0$ defined as

${\Phi }_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\lambda \left(\left(x/2\right)\left(1-x-0.2z\right)\mathrm{,}\left(y/2\right)\left(1-y-0.2z\right)\mathrm{,}z\left(1-z\right)\right)\mathrm{,}$ (56)

where

$\mathfrak{X}:={\mathfrak{X}}_{\Phi }\mathrm{<mo>\; :\; </mo>}=\left\{\mathrm{<mo\; stretchy="false">\; (\; </mo>}x\mathrm{,}y\mathrm{,}z\mathrm{<mo\; stretchy="false">\; )\; </mo>}\in {ℝ}^3\mathrm{<mo>\; :\; </mo>0}\le x\le 1-0.2z\mathrm{,0}\le y\le 1-0.2z\mathrm{,0}\le z\le 1\right\}\mathrm{.}$ (57)

It should be noted that $0\le \lambda \le 4$ is an admissible parameter range for this family of maps.

5.1.1. Types of Fixed Points of the Map(56)

The fixed points of (56) are easily computed to be

$\begin{array}{l}{p}_0\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,0,0}\right)\mathrm{<mo>\; ;\; </mo>}{p}_1\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{\lambda -2}{\lambda }\mathrm{,0,0}\right)\mathrm{<mo>\; ;\; </mo>}\\ p_2\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,}\frac{\lambda -2}{\lambda }\mathrm{,0}\right)\mathrm{<mo>\; ;\; </mo>}{p}_3\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,0,}\frac{\lambda -1}{\lambda }\right)\mathrm{<mo>\; ;\; </mo>}\\ p_4\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{\lambda -2}{\lambda }\mathrm{,}\frac{\lambda -2}{\lambda }\mathrm{,0}\right)\mathrm{<mo>\; ;\; </mo>}{p}_5\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{4\lambda -9}{5\lambda }\mathrm{,0,}\frac{\lambda -1}{\lambda }\right)\mathrm{<mo>\; ;\; </mo>}\\ p_6\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,}\frac{4\lambda -9}{5\lambda }\mathrm{,}\frac{\lambda -1}{\lambda }\right)\mathrm{<mo>\; ;\; </mo>}{p}_7\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{4\lambda -9}{5\lambda }\mathrm{,}\frac{4\lambda -9}{5\lambda }\mathrm{,}\frac{\lambda -1}{\lambda }\right)\mathrm{.}\end{array}$ (58)

We can then determine the types of these fixed points from the spectral properties of the derivative matrix

${{\Phi }^{\prime }}_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=\left[\begin{array}{ccc}\frac{\lambda }{2}\left(1-2x-0.2z\right)& 0& -0.1\lambda x\\ 0& \frac{\lambda }{2}\left(1-2y-0.2z\right)& -0.1\lambda y\\ 0& 0& \lambda \left(1-2z\right)\end{array}\right]\mathrm{.}$ (59)

As the above matrix is upper triangular, the eigenvalues at the fixed points are just the diagonal elements evaluated at the fixed $\left(x_{\lambda }\mathrm{,}{y}_{\lambda }\mathrm{,}{z}_{\lambda }\right)$ ; namely,

$\begin{array}{l}{p}_0\mathrm{<mo>\; :\; </mo>}{\mu }_1={\mu }_2=\frac{\lambda }{2}\mathrm{,}{\mu }_3=\lambda \mathrm{<mo>\; ;\; </mo>}{p}_1\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1=\frac{4-\lambda }{2}\mathrm{,}{\mu }_2=\frac{\lambda }{2}\mathrm{,}{\mu }_3=\lambda \mathrm{<mo>\; ;\; </mo>}\\ p_2\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1=\frac{\lambda }{2}\mathrm{,}{\mu }_2=\frac{4-\lambda }{2}\mathrm{,}{\mu }_3=\lambda \mathrm{<mo>\; ;\; </mo>}{p}_3\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1={\mu }_2=\frac{4\lambda +1}{10}\mathrm{,}{\mu }_3=2-\lambda \mathrm{<mo>\; ;\; </mo>}\\ p_4\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1={\mu }_2=\frac{4-\lambda }{2}\mathrm{,}{\mu }_3=\lambda \mathrm{<mo>\; ;\; </mo>}\\ p_5\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1=\frac{19-4\lambda }{10}\mathrm{,}{\mu }_2=\frac{4\lambda +1}{10}\mathrm{,}{\mu }_3=2-\lambda \mathrm{<mo>\; ;\; </mo>}\end{array}$

$\begin{array}{l}{p}_6\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1=\frac{4\lambda +1}{10}\mathrm{,}{\mu }_2\mathrm{<mo>\; =\; </mo>}\frac{19-4\lambda }{10}\mathrm{,}{\mu }_3=2-\lambda \mathrm{<mo>\; ;\; </mo>}\\ p_7\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}{\mu }_1={\mu }_2=\frac{19-4\lambda }{10}\mathrm{,}{\mu }_3=2-\lambda \mathrm{,}\end{array}$ (60)

which determine the types of the fixed points.

5.1.2. Dynamics in $\mathfrak{X}$ for $0<\lambda \le 1$

In this case, it follows from the above that the only fixed point in $\mathfrak{X}$ is the origin $p_0$, which is a hyperbolic sink when $0<\lambda <1$ and still a sink when $\lambda =1$. Moreover, $\left\{p_0\right\}$ is a global attractor in $\mathfrak{X}$ for all $0<\lambda \le 1$.

5.1.3. Dynamics in $\mathfrak{X}$ for $1<\lambda \le 2$

The dynamics in $\mathfrak{X}$ for this case is still rather easy to describe. The only fixed points are $p_0$ and $p_3\left(\lambda \right)\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,0,}\frac{\lambda -1}{\lambda }\right)$. The eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_0\right)$ are

(from (60)) ${\mu }_1={\mu }_2=\lambda /2$ and ${\mu }_3=\lambda$, with corresponding eigenvectors along the x-, y- and z-axes, respectively, so $p_0$ is a hyperbolic saddle point for $1<\lambda <2$, with a 2-dimensional stable $W^s\left(p_0\right)$ and 1-dimensional unstable manifold $W^u\left(p_0\right)$. At $\lambda =2$, the stable manifold becomes a center manifold $W^c\left(p_0\right)$ that still attracts all points in the portion of $\mathfrak{X}$ in the xy-plane. The eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_3\right)$ are ${\mu }_1={\mu }_2=0.4\lambda +0.1$ and ${\mu }_3=2-\lambda$, from which it readily follows that $p_3$ is a hyperbolic sink for all $1<\lambda \le 2$ and the basin of attraction (in $\mathfrak{X}$ ) is $\mathcal{B}\left(\mathrm{<mo>\; \{\; </mo>}{p}_3\mathrm{<mo>\; \}\; </mo>}\right)=\mathfrak{X}\cap \left\{\left(x,y,z\right):0<z<1\right\}$. We note that $\lambda =1$ is actually a bifurcation value if the map is considered on all of ${�}^3$ inasmuch as $p_0$ and $p_3\left(\lambda \right)$ are distinct for $0<\lambda <1$, merge at $\lambda =1$, then re-emerge for $1<\lambda \le 2$ as distinct points, both of which are in $\mathfrak{X}$.

5.1.4. Dynamics in $\mathfrak{X}$ for $2<\lambda \le 3$

For this case, things become a bit more interesting. To begin with, we have to consider two subcases: (i) $2<\lambda \le 9/4$ ; (ii) $9/4<\lambda \le 3$.

(i) $2<\lambda \le 9/4$ : For this subcase, the fixed points in $\mathfrak{X}$ are $p_0$, $p_1\left(\lambda \right)$, $p_2\left(\lambda \right)$, $p_3\left(\lambda \right)$ and $p_4\left(\lambda \right)$. It follows from (60) that the eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_0\right)$ are ${\mu }_1={\mu }_2=\lambda /2$ and ${\mu }_3=\lambda$, so $p_0$ is a hyperbolic source, actually for all $2<\lambda \le 3$. On the other hand, the eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_1\right)$ are

${\mu }_1=\frac{4-\lambda }{2}$, so $\left|{\mu }_1\right|<1$ for all $2<\lambda \le 3$, ${\mu }_2=\frac{\lambda }{2}$, so $\left|{\mu }_2\right|>1$ for all $2<\lambda \le 3$ and ${\mu }_3=\lambda$, so $\left|{\mu }_3\right|>1$ for all $2<\lambda \le 3$. Hence, $p_1$ is a hyperbolic saddle point for all $2<\lambda \le 3$. Next, the eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_2\right)$ are ${\mu }_1=\frac{\lambda }{2}$, so $\left|{\mu }_1\right|>1$ for all $2<\lambda \le 3$, ${\mu }_2=\frac{4-\lambda }{2}$, so $\left|{\mu }_2\right|<1$ for all $2<\lambda \le 3$ and ${\mu }_3=\lambda$, so $\left|{\mu }_3\right|>1$ for all $2<\lambda \le 3$. Consequently, $p_2$ is also (symmetrically) a hyperbolic saddle point for all $2<\lambda \le 3$.

The eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_3\right)$ are ${\mu }_1={\mu }_2=\frac{4\lambda +1}{10}$ and ${\mu }_3=2-\lambda$. Therefore, $0<{\mu }_1={\mu }_2<1$ for $2<\lambda <9/4$ and $-1<{\mu }_3<0$ for all $2<\lambda \le 3$ and

${\mu }_3=-1$ when $\lambda =3$. Hence, there is a flip bifurcation (along the z-axis) at $\lambda =3$. It also follows from the eigenvalues and an examination of the map (56) that $p_3$ is a hyperbolic sink for $2<\lambda <9/4$ and still a sink for $2<\lambda \le 9/4$, but it changes into a hyperbolic saddle point for $9/4<\lambda <3$. Thus, $\lambda =9/4$ is a bifurcation point on account of this behavior, but also due to the emergence and re-emergence of fixed points across this value since $p_3=p_5=p_6=p_7$ when $\lambda =9/4$ and these points are all different for $\lambda >9/4$.

Observe that $\lambda =2$ is also a bifurcation value since $p_0$, $p_1$, $p_2$ and $p_4$ are all equal when $\lambda =2$, but are distinct when $\lambda \ne 2$. The eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_4\right)$ are ${\mu }_1={\mu }_2=\frac{4-\lambda }{2}$ and ${\mu }_3=\lambda$ from which we conclude that $p_4$

is a hyperbolic saddle point for all $2<\lambda \le 3$, since $0<{\mu }_1={\mu }_2<1<{\mu }_3$ in this parameter range. It is worth noting that the generalized eigenspace corresponding to ${\mu }_1={\mu }_2$ lies in the xy-plane, which happens to contain the stable manifold of $p_4$. We also conclude from these considerations that $\left\{p_3\right\}$ is an attractor for all $1<\lambda \le 9/4$ with basin of attraction given by $\mathcal{B}\left(\left\{p_3\right\}\right)=\mathfrak{X}\cap \left\{\left(x,y,z\right):0<z<1\right\}$.

(ii) $9/4<\lambda \le 3$ : In this range, we must add the fixed points $p_5$, $p_6$ and $p_7$, noting once again that $p_3=p_5=p_6=p_7$ when $\lambda =9/4$. Moreover, $p_3$, $p_5$, $p_6$ and $p_7$ are distinct, with only $p_3\in \mathfrak{X}$ for $2<\lambda <9/4$ and all of them are distinct when $9/4<\lambda \le 3$. The fixed points $p_5$ and $p_6$ are essentially

symmetric and ${{\Phi }^{\prime }}_{\lambda }\left(p_5\right)$ and ${{\Phi }^{\prime }}_{\lambda }\left(p_6\right)$ have the eigenvalues ${\mu }_1=\frac{19-4\lambda }{10}$, ${\mu }_2=\frac{4\lambda +1}{10}$, ${\mu }_3=2-\lambda$ and ${\mu }_1=\frac{4\lambda +1}{10}$, ${\mu }_2=\frac{19-4\lambda }{10}$, ${\mu }_3=2-\lambda$, respectively. Consequently, both $p_5$ and $p_6$ are hyperbolic saddle points for all $9/4<\lambda <3$ and a closer examination of the map shows the dynamics in neighborhoods of these points still exhibit saddle point behavior for $\lambda =3$. The eigenvalues of ${{\Phi }^{\prime }}_{\lambda }\left(p_7\right)$ are ${\mu }_1={\mu }_2=\frac{19-4\lambda }{10}$, ${\mu }_3=2-\lambda$, so it follows that $\left\{p_7\right\}$ is a hyperbolic sink for $9/4\le \lambda <3$. In fact, an analysis of (56) shows that $\left\{p_7\right\}$ is an attractor whenever $9/4\le \lambda \le 3$ with $\mathcal{B}\left(\left\{p_7\right\}\right)=\mathfrak{X}$.

5.1.5. Dynamics in $\mathfrak{X}$, Including Chaotic Strange Attractors, for $3<\lambda \le 4$

This is the parameter range where we find the most interesting and complicated dynamics, but some features are rather simple or are subsumed by attributes discussed above. For example, the origin $p_0$ is a hyperbolic source for this range, while $p_1$ and $p_2$ are both hyperbolic saddle points with 1-dimensional stable and 2-dimensional unstable manifolds for $3<\lambda \le 4$ and at $\lambda =4$ these fixed points have a stable manifold that is actually superstable. Furthermore, $p_3$ is a hyperbolic source for all $3<\lambda \le 4$ that generates a period-doubling cascade along the z-axis leading to chaos in the manner of the flip bifurcations described in Section 3. Associated with this cascade is a sequence of hyperbolic fixed points of ${\Phi }_{\lambda }^{2n}$ that change from hyperbolic saddle points to sources as $\lambda$ increases. The fixed point $p_4$ continues to be a hyperbolic saddle point (with stable manifold $W^s\left(p_4\left(\lambda \right)\right)=\mathfrak{X}\cap \left\{xy\text{-plane}\right\}$ ) for this parameter range. Both $p_5$ and $p_6$ are hyperbolic saddle points with a 1-dimensional stable manifold and 2-dimensional unstable manifold.

It is, in fact, the changes in the fixed point $p_7\left(\lambda \right)$ in this parameter range that are the sources of the especially complicated dynamics, which ultimately includes a chaotic strange attractor. As shown in (60), the eigenvalues of ${{\Phi }^{\prime }}_{\lambda }(\; p\; 7\; )$

are ${\mu }_1={\mu }_2=\frac{19-4\lambda }{10}$ and ${\mu }_3=2-\lambda$, so that $p_7\left(\lambda \right)$ is hyperbolic, with a

2-dimensional stable manifold and 1-dimensional unstable manifolds over the whole parameter range. Observe that $\left|{\mu }_1\right|=\left|{\mu }_2\right|$ decreases and $\left|{\mu }_3\right|$ increases with increasing $\lambda$ for $3\le \lambda \le 4$. Therefore, there is a flip bifurcation at $\lambda =3$ with a period-doubling cascade leading to chaos along the unstable manifold $W^u\left(p_7\left(\lambda \right)\right)$ as $\lambda$ increases beyond 3 in the manner described in Section 3. Moreover, the sequence of attracting 2n-cycles $C_n$ in $W_{\text{lin}}^u\left(p_7\left(\lambda \right)\right)$ actually comprise attractors with basins of attraction $\mathcal{B}\left(C_n\right)=\mathfrak{X}$ for all natural numbers $n$. These cycle attractors (periodic sinks) can also be viewed as results of tangent intersections of $W^u\left(p_7\left(\lambda \right)\right)$ and $W^s\left(p_7\left(\lambda \right)\right)$ in the manner of Newhouse theory [24].

More interesting attractors chaotic strange ones rooted in the nature of the 2-dimensional stable and 1-dimensional unstable manifolds of $p_7\left(\lambda \right)$, respectively, occur when $\lambda$ is near four. These horseshoe generated attractors are described in the following result and illustrated geometrically and with simulations of the sequence of bifurcations preceding them, respectively, in Figure 7 and Figure 8.

Theorem 7. Suppose that the special cases of the map (54) $F\mathrm{<mo>\; :\; </mo>}{\mathfrak{X}}_F\to {\mathfrak{X}}_F$ defined as

$\begin{array}{l}F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}{F}_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)=\left(f_1\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_2\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}{f}_3\left(x\mathrm{,}y\mathrm{,}z\right)\right)\\ \mathrm{<mo>\; :\; </mo>}=(\left(\lambda /2\right)x\left(1-x-{\gamma }_{12}y-{\gamma }_{13}z\right)\mathrm{,}\left(\lambda /2\right)y\left(1-{\gamma }_{21}x-y-{\gamma }_{23}z\right)\mathrm{,}\\ \lambda z\left(1-{\gamma }_{31}x-{\gamma }_{32}y-z\right))\mathrm{,}\end{array}$ (61)

where ${\mathfrak{X}}_F$ is as in (55), satisfies the following properties: (i) $3.7<\lambda \le {\lambda }_a$ ; (ii) $0\le {\gamma }_{12}\mathrm{,}{\gamma }_{21}\mathrm{,}{\gamma }_{31}\mathrm{,}{\gamma }_{32}\le 0.05$ and sufficiently small and (iii) $\left|{\gamma }_{13}-0.2\right|\mathrm{,}\left|{\gamma }_{13}-0.2\right|$ are also sufficiently small. Here ${\lambda }_a$ is the maximum of the admissible parameter range $\left[\mathrm{0,}{\lambda }_a\right]$, which is greater than 3.8 in view of Lemma 1. Then, there is a horseshoe-like chaotic strange attractor (CSA) of the generalized attracting horseshoe (GAH) type in ${\mathfrak{X}}_F$ of the form

${\mathfrak{A}}_{\lambda }\mathrm{<mo>\; :\; </mo>}={\displaystyle {\cap }_{n=0}^{\infty }}{F}_{\lambda }^n\left(K\right)\mathrm{,}$ (62)

with basin of attraction $\mathcal{B}\left({\mathfrak{A}}_{\lambda }\right)\supset {\mathfrak{X}}_F\supset K$, where K is compact.

Proof. Obviously, the map depends smoothly on the parameters. Moreover, as we shall see, the principal features of our proof, both analytical and qualitative,

are smoothly persistent under small perturbations of the interaction coefficients. Therefore, owing to the manner in which the theorem is formulated, it suffices to verify the result for the rather simple system just studied above; namely,

$F\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}=F_{\lambda }\left(x\mathrm{,}y\mathrm{,}z\right)=\lambda \left(\left(x/2\right)\left(1-x-0.2z\right)\mathrm{,}\left(y/2\right)\left(1-y-0.2z\right)\mathrm{,}z\left(1-z\right)\right)\mathrm{,}$ (63)

with ${\lambda }_a=4$, where

$\mathfrak{X}\mathrm{<mo>\; :\; </mo>}={\mathfrak{X}}_F\mathrm{<mo>\; :\; </mo>}=\left\{\left(x\mathrm{,}y\mathrm{,}z\right)\in {ℝ}^3\mathrm{<mo>\; :\; </mo>0}\le x\le 1-0.2z\mathrm{,0}\le y\le 1-0.2z\mathrm{,0}\le z\le 1\right\}\mathrm{.}$ (64)

Let Q be the solid prism having the following vertices in the $z=0$ plane: $\left(\frac{3}{20}\mathrm{,}\frac{5}{20}\mathrm{,0}\right)$, $\left(\frac{5}{20}\mathrm{,}\frac{3}{20}\mathrm{,0}\right)$, $\left(\frac{11}{20}\mathrm{,}\frac{9}{20}\mathrm{,0}\right)$ and $\left(\frac{9}{20}\mathrm{,}\frac{11}{20}\mathrm{,0}\right)$ and corresponding vertices in the plane $z=1$ : $\left(\frac{3}{20}\mathrm{,}\frac{5}{20}\mathrm{,1}\right)$, $\left(\frac{5}{20}\mathrm{,}\frac{3}{20}\mathrm{,1}\right)$, $\left(\frac{11}{20}\mathrm{,}\frac{9}{20}\mathrm{,1}\right)$ and

$\left(\frac{9}{20}\mathrm{,}\frac{11}{20}\mathrm{,1}\right)$ as shown in Figure 7. We note that both the upper and lower faces

of Q are mapped into the $x\mathrm{,}y$ -plane ( $z=0$ ). Moreover, it is easy to verify that Q is a trapping set for the map for $3.7<\lambda \le {\lambda }_a$ ; i.e., $F_{\lambda }\left(Q\right)\subset Q$ for all $\lambda \in \left(\mathrm{3.7,}{\lambda }_a\right]$ and the image has the horseshoe shape depicted in Figure 7. This and other properties are most conveniently apprehended by analyzing the images under F of the vertical fibers

$l\left(x_0\mathrm{,}{y}_0\right)\mathrm{<mo>\; :\; </mo>}=\left\{\left(x_0\mathrm{,}{y}_0\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}\left(x_0\mathrm{,}{y}_0\mathrm{,0}\right)\in R\mathrm{,0}\le z\le 1\right\}$

where R is the rectangular bottom face of Q, which lies in the xy-plane. We compute that

$F_{\lambda }\left(l\left(x_0\mathrm{,}{y}_0\right)\right)=\left\{\left(\lambda /2\right)\left(x_0\left(1-x_0-0.2z\right)\mathrm{,}{y}_0\left(1-y_0-0.2z\right)\mathrm{,2}z\left(1-z\right)\right)\mathrm{<mo>\; :\; </mo>0}\le z\le 1\right\}\mathrm{,}$

which is a parabolic curve in Q that begins at

$\left(\lambda /2\right)\left(x_0\left(1-x_0\right)\mathrm{,}{y}_0\left(1-y_0\right)\mathrm{,0}\right)$

Figure 7. Depiction of generators of the GAH: the trapping set Q and $F\left(Q\right)$ for (63).

and ends at

$\left(\lambda /2\right)\left(x_0\left(0.8-x_0\right)\mathrm{,}{y}_0\left(0.8-y_0\right)\mathrm{,0}\right)\mathrm{.}$

For example, the parabolic curve $F_{\lambda }\left(l\left(\frac{1}{5}\mathrm{,}\frac{1}{5}\right)\right)$ begins at $\frac{2\lambda }{25}\left(\mathrm{1,1,0}\right)$, which lies between $0.296\left(\mathrm{1,1,0}\right)$ and $0.32\left(\mathrm{1,1,0}\right)$, while $F_{\lambda }\left(l\left(\frac{1}{2}\mathrm{,}\frac{1}{2}\right)\right)$ begins at $\frac{\lambda }{8}\left(\mathrm{1,1,0}\right)$, which lies between $0.4625\left(\mathrm{1,1,0}\right)$ and $0.5\left(\mathrm{1,1,0}\right)$ for $3.7\le \lambda \le 4$.

Then it is straightforward to verify that combining the images of all the vertical fibers in Q produces the (twisted) horseshoe shaped image $F\left(Q\right)$ shown in Figure 7. Naturally, the hyperbolic fixed point $p_7\left(\lambda \right)$ belongs to $F_{\lambda }\left(Q\right)$ and also to

${\mathfrak{A}}_{\lambda }\mathrm{<mo>\; :\; </mo>}={\displaystyle {\cap }_{n=0}^{\infty }}{F}_{\lambda }^n\left(Q\right)\mathrm{.}$

For $3.7\le \lambda \le 4$ these sets ${\mathfrak{A}}_{\lambda }$ by virtue of their construction and Birkhoff-Moser-Smale theory (cf. [17] [25] [26] [27] [28] ), attracts all points in the ${\stackrel{\text{o}}{\mathfrak{X}}}_{F_{\lambda }}$ and are fractal $F_{\lambda }$ -invariant sets on which the restriction $F_{\mathrm{<mo>\; |\; </mo>}{\mathfrak{A}}_{\lambda }}$ is chaotic. Consequently, the ${\mathfrak{A}}_{\lambda }$ are chaotic, strange attracting sets. To complete the proof, it remains to show that for $\lambda \in \left[\mathrm{3.7,4}\right]$ each ${\mathfrak{A}}_{\lambda }$ is minimally attracting instead of a periodic sink subset generated by tangent intersections of

$W^u\left(p_7\left(4\right)\right)$ and $W^s\left(p_7\left(4\right)\right)=\left\{\left(x\mathrm{,}y\mathrm{,}z\right)\in {ℝ}^3\mathrm{<mo>\; :\; </mo>}z=z_{\lambda }=\left(\lambda -1\right)/\lambda \right\}\cap {\stackrel{\text{o}}{\mathfrak{X}}}_{F_{\lambda }}$ in the manner of Newhouse theory [24]. As a matter of fact, it follows that such ${\mathfrak{A}}_{\lambda }$ are generalized horseshoe attractors of the type described in [23] if we can verify that there are no tangent intersections of $W^u\left(p_7\left(\lambda \right)\right)$ and $W^s\left(p_7\left(\lambda \right)\right)$ for $\lambda \in \left[\mathrm{3.7,4}\right]$ by verifying that the $F_{\lambda }$ -image of the arch of $F_{\lambda }\left(Q\right)$ lies well below

the plane $z=z_{\lambda }=\frac{\lambda -1}{\lambda }$. We shall describe the arch of the horseshoe by finding

estimates for the maximum height (z-value) of the fiber images $F_{\lambda }\left(l\left(x_0\mathrm{,}{y}_0\right)\right)$ as $\left(x_0\mathrm{,}{y}_0\mathrm{,0}\right)$ ranges over R. For this, we define

$Z_{\lambda }\left(z\right):=Z_{\left(x_0,y_0\right)}\left(z\right):=\lambda z\left[1-z\right]$

and find its maximizer $z_M\mathrm{<mo>\; :\; </mo>}=z_M\left(x_0\mathrm{,}{y}_0\right)$ and maximum height $Z\left(z_M\right)$, which are readily found to be $z_M=0.5$ and $Z_{\lambda }\left(z_M\right)=0.25\lambda$. Consequently, $0.925\le Z_{\lambda }\left(z_M\right)\le 1$ on R whenever $3.7\le \lambda \le 4$. Hence, the arch of $F_{\lambda }(\; Q\; )$

certainly lies in $A=A\left(0.9\right)\mathrm{<mo>\; :\; </mo>}=Q\cap \left\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\ge 0.9\right\}$ and it follows that $F_{\lambda }\left(A\right)\subset Q\cap \left\{\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{<mo>\; :\; </mo>}z\le 0.36\right\}$ for all $\lambda \in \left[\mathrm{3.7,4}\right]$. Thus, since $\frac{27}{37}\le z_{\lambda }\le \frac{3}{4}$

for all $3.7\le \lambda \le 4$, there are no tangent intersections of $W^u\left(p_7\left(\lambda \right)\right)$ and $W^s\left(p_7\left(\lambda \right)\right)$ for the given range of the birth-rate parameters, which means that each ${\mathfrak{A}}_{\lambda }$ is a GAH and the proof is complete.

Observe that the CSA shown in the simulation in Figure 8 has a rather distinctive shape that brings to mind a chaotic candy cane. It also should be noted that

Figure 8. Bifurcations for (61): (a) $\lambda =3.0$ ; (b) $\lambda =3.25$ ; (c) $\lambda =3.5$ ; (d) $\lambda =3.75$ with ${\gamma }_{12}=0.05$, ${\gamma }_{13}=0.2$, ${\gamma }_{21}=0.03$, ${\gamma }_{23}=0.21$, ${\gamma }_{31}=0.04$ and ${\gamma }_{32}=0.05$.

the horseshoe of the above theorem is analogous to the twisted horseshoe for the 2-dimensional Lotka-Volterra map [5].

6. Conclusions and Suggestions

In this study of the discrete Lotka-Volterra dynamical system model, we have shown not only that there is an abundance of fully 3-dimensional flip bifurcations and related period-doubling cascades to chaos, but some rather novel analogs of Neimark-Sacker bifurcations and chaotic strange attractors. These bifurcations have been analyzed in considerable detail exploiting, for example, the identification of the strange chaotic attractors as generalized attracting horseshoes, which are shaped like candy canes. The focus in this investigation was on 3-dimensional systems, but in future research we plan to identify and analyze higher dimensional analogs of the Neimark-Sacker-type bifurcations, candy cane attractors and related dynamical phenomena. Moreover, we shall make an effort to find other examples of interesting and useful discrete dynamical system models, besides those of Lotka-Volterra type, that exhibit these and possibly other novel, unusual or unexpected types of behavior.

Acknowledgements

Y. Joshi is grateful to the Department of Mathematics and Computer Science at Kingsborough Community College, CUNY for its support, M. Savescu wishes to thank the Department of Mathematics at Kutztown University for encouragement and support and M. Syed and D. Blackmore are beholden to the Department of Mathematical Sciences and CAMS at NJIT for generous support. Also thanks are due the editorial staff of JAM for its fast and efficient work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lotka, A. (1910) Contribution to the Theory of Periodic Reaction. Journal of Physical Chemistry, 14, 271-274. https://doi.org/10.1021/j150111a004
[2]	Volterra, V. (1926) Variazioni del Numero d’Individuii in Specie Animali Conviventi. Memorie Accademia dei Lincei di Roma, 2, 31-113.
[3]	Murray, J. (1993) Mathematical Biology. Springer, New York. https://doi.org/10.1007/978-3-662-08542-4
[4]	Smith, J.M. (1978) Models in Ecology. Cambridge University Press, Cambridge.
[5]	Blackmore, D., Chen, J., Perez, J. and Savescu, M. (2001) Dynamical Properties of Discrete Lotka-Volterra Equations. Chaos, Solitons and Fractals, 12, 2553-2568. https://doi.org/10.1016/S0960-0779(00)00214-9
[6]	Din, Q. (2013) Dynamics of a Discrete Lotka-Volterra Model. Advances in Difference Equations, 2013, Article No. 95. https://doi.org/10.1186/1687-1847-2013-95
[7]	Liu, D. and Xiao, D. (2006) Bifurcations in a Discrete Lotka-Volterra Predator-Prey System. Discrete and Continuous Dynamical Systems B, 6, 559-572. https://doi.org/10.3934/dcdsb.2006.6.559
[8]	Zhao, M., Yuan, Z. and Li, C. (2016) Dynamics of a Discrete-Time Predator-Prey System. Advances in Difference Equations, 2016, Article No. 191. https://doi.org/10.1186/s13662-016-0903-6
[9]	Roeger, L.-I. and Lahodny, Jr.G. (2013) Dynamically Consistent Discrete Lotka-Volterra Competition Systems. Journal of Difference Equations and Applications, 19, 191-200. https://doi.org/10.1080/10236198.2011.621894
[10]	Bischi, G. and Tramontana, F. (2010) Three-Dimensional Discrete-Time Lotka-Volterra Models with an Application to Industrial Clusters. Journal of Difference Equations and Their Applications, 15, 3000-3014. https://doi.org/10.1016/j.cnsns.2009.10.021
[11]	Choo, S. (2014) Global Stabiliy in n-Dimensional Discrete Lotka-Volterra Predator-Prey Models. Advances in Difference Equations, 2014, Article No. 11. https://doi.org/10.1186/1687-1847-2014-11
[12]	Lyu, X. and Jablonski, P.G. (2012) Four-Dimensional Discrete-Time Lotka-Volterra Models with an Application to Ecology.
[13]	Mickens, R. (2016) A Note on Exact Finite Difference Schemes for Modified Lotka-Volterra Differential Equations. Journal of Difference Equations and Applications, 24, 1016-1022. https://doi.org/10.1080/10236198.2018.1430792
[14]	Mukhamedov, F. and Saburov, M. (2012) On Discrete Lotka-Volterra Type Models. International Journal of Modern Physics: Conference Series, 9, 341-346. https://doi.org/10.1142/S2010194512005405
[15]	Yousef, A., Salmon, S. and Elsadany, A. (2018) Stability and Bifurcation of a Delayed Predator-Prey Model. International Journal of Bifurcation and Chaos, 28, 1850116-1850129. https://doi.org/10.1142/S021812741850116X
[16]	Kuznetsov, Y. (2004) Elements of Applied Bifurcation Theory. 3rd Edition, Springer, New York. https://doi.org/10.1007/978-1-4757-3978-7
[17]	Wiggins, S. (2014) Global Bifurcations and Chaos. Springer, New York.
[18]	Milnor, J. (1985) On the Concept of Attractor. Communications in Mathematical Physics, 99, 177-195. https://doi.org/10.1007/BF01212280
[19]	Neimark, J. (1959) On Some Cases of Periodic Motion Depending on Parameters. Doklady Akademii Nauk SSSR, 129, 36-39.
[20]	Sacker, R. (1964) On Invariant Surfaces and Bifurcations of Periodic Solutions of Ordinary Differential Equations. Report IMM-NYU 333, 1-62.
[21]	Champanerkar, J. and Blackmore, D. (2007) Pitchfork Bifurcations of Invariant Manifolds. Topology and Its Applications, 154, 1650-1663. https://doi.org/10.1016/j.topol.2006.12.014
[22]	Joshi, Y. and Blackmore, D. (2014) Strange Attractors for Asymptotically Zero Maps. Chaos, Solitons and Fractals, 68, 123-138. https://doi.org/10.1016/j.chaos.2014.08.005
[23]	Joshi, Y., Blackmore, D. and Rahman, A. (2016) Generalized Attracting Horseshoes and Chaotic Strange Attractors.
[24]	Newhouse, S. (1974) Diffeomorphisms with Infinitely Many Sinks. Topology, 13, 9-18. https://doi.org/10.1016/0040-9383(74)90034-2
[25]	Katok, A. and Hassleblatt, B. (1995) Introduction to the Modern Theory of Dynamical Systems. Cambridge University Press, Cambridge.
[26]	Moser, J. (1973) Stable and Random Motion in Dynamical Systems. Princeton University Press, Princeton.
[27]	Smale, S. (1980) The Mathematics of Time: Essays on Dynamical Systems, Economic Processes, and Related Topics. Springer, New York. https://doi.org/10.1007/978-1-4613-8101-3
[28]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, New York.