We shall consider a natural discrete analog of the famous and ubiquitous Lotka-Volterra model, developed independently by Alfred Lotka  and Vito Volterra  for the dynamics of a population comprising several interacting species, say , as described in such references as  . This model may be expressed in the form
where the coefficients are real constants and typically nonnegative, with for all , . 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:
where the birth-rate parameters and interaction coefficients are usually nonnegative constants as we shall assume in the sequel, with for all , . 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 , 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 comprising a discrete (semi-) dynamical system and is represented as
where the coefficients are as described in (2). Here, the ’s and ’s represent the birth-rates and interaction strengths of and among the populations, respectively.
Inasmuch as we are primarily thinking of as a population vector, it is natural to consider the map F restricted to
Moreover, since we want the forward iterates of F of points in this set to remain in , we would like to find the largest subset such that for sizeable ranges of the birth-rate parameters. For example, if all the off-diagonal interaction coefficients are zero, the obvious choice of is the m-cube in which case for all , . 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
and the admissible ranges have to be adjusted accordingly.
It is clear from (3) that a necessary and sufficient condition for is that the following system of inequalities be satisfied:
where . One can by lengthy but straightforward calculations find the maximizers and associated maximum values of each of the , which yield a rather complicated system of inequalities determining the optimal admissible values of . 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
are admissible in for the discrete dynamical system generated by the iterates of the map (3).
Proof. It follows from (4) and (5) that for all and , , so only the upper inequalities in (6) require verification. For this, we observe that it is easy to see that each on . Consequently,
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 defined as
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     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  - , 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  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. , 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 , and and the coordinate planes by , and :
Theorem 2. The origin is a fixed point of F and all of the coordinate sets defined above are F-invariant, i.e. , , , , and .
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
which has a unique solution if the matrix 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 with and positive, namely with , and , so that the map takes the form
with corresponding fixed point .
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  , 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 defined by
with and the single (bifurcation) parameter . It follows from Lemma 1 that is an admissible range for the above map. The coexistence fixed points of this map as functions of the parameter are readily found to be
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 and :
i) The system has a flip bifurcation at for at which the restriction to the -invariant line (segment) has derivative .
ii) The restriction f has a period-doubling cascade leading to chaos along as increases from 3 to a limit of approximately 3.57 following that of the standard 1-dimensional logistic map.
iii) The fixed point is a sink for and is hyperbolic for with a 2-dimensional (planar) stable manifold and 1-dimensional (linear) unstable manifold , where denotes the interior of .
iv) When , the fixed point has the stable manifold and the center manifold .
v) There is a Milnor attractor along , as in , at the limit of the period-doubling cascade as well as a dense chaotic (standard tent map-type Milnor) attractor for , each with a basin of attraction .
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) ; (b) ; (c) ; (d) with and .
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 , with equal birth rates serving as the bifurcation parameter (with ), shown below. The dynamics shall be analyzed in detail for prescribed ranges of variation of the interaction coefficients for
For the coexistence fixed point, we simply solve the following system of equations in which we impose the initial constraints :
to obtain the unique solution
In order to better understand the flip bifurcations that can occur for the map (12), we first consider the following special case:
with . It follows from Lemma 1 that is an admissible range for this map. Then, from (13) we find that the unique coexistence fixed point is
The type of this fixed point is determined by the eigenvalues of the derivative matrix
which are , and , with corresponding eigenfunctions , and for , respectively. We note that for all , so that , which is -invariant for all , comprises the stable manifold over the same parameter interval. Moreover, for , when and for . It follows of course, that there is a flip bifurcation at at which value the linearized center manifold is and the ( -invariant) stable manifold is . When , the stable manifold is and the ( -invariant) unstable manifold has as its linearization the line , which is not -invariant. In this parameter range, the attractor is a 2-cycle on . The flip bifurcation and period-doubling cascade converging to a Milnor attractor is determined by the restriction of to the continuation of as a 1-dimensional -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 -dependent translation by and the corresponding translated map as
where the coexistence point for the original map in the translated coordinates now corresponds to for all . We assume that the invariant 1-dimensional invariant manifold that begins as (super) attracting for , becomes the center manifold at and is an unstable manifold for can be continued as corresponding to the eigenvalue of largest magnitude for , is analytic in and expressible as
with coefficients that are analytic functions of the parameter, converge on a -interval containing the origin and large enough so that the restriction exhibits a flip bifurcation at and a period-doubling cascade to converging to a Milnor attractor .
For invariance, we must have
with , which we want to solve for . After finding , we wish to determine a parameter interval containing in its interior a such that has a Milnor attractor that is the period-doubling cascade limit for as . Substituting (20) in (21) and (22) and equating coefficients of like powers of , we obtain the values of the leading order coefficients and the following recurrence relations for the higher order terms ( ):
where , and is the usual greatest integer function. Several terms in the series (20) computed using the above recursion formulas are as follows:
Now, it a straightforward matter to prove, for example using the method of majorants, that the power series (27) converge for , 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 defined by (15) has the following properties for and .
(a) The coexistence fixed point given by (16) is a sink for , has a center manifold at (associated with an eigenvalue of ) along which (15) has a flip bifurcation and there is a stable manifold .
(b) For , is hyperbolic, with stable manifold and unstable manifold along which there is an 2-cycle attractor.
(c) possesses a smooth invariant continuation for and has a period-doubling cascade sequence such that 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
Figure 2. Flip bifurcation for (12): (a) ; (b) ; (c) ; (d) with , , , , and .
and d differ from zero, satisfies the following properties:
Then the dynamical system generated by the map (12) has a flip bifurcation at for a value of the bifurcation parameter near 3 along what is initially a center and then an unstable manifold followed by a smooth invariant extension containing a period-doubling cascade to a (Milnor) attractor as 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 and d as well as the bifurcation parameter ) is
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
which we shall compare with the characteristic equation for (15); namely,
and , . We shall compare these with the coefficients of (31), which are readily found (in terms of the fixed point (29)) to be
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) for all .
(p2) Two of the eigenvalues (solutions of (31)) of , with corresponding eigenvectors nearly parallel to the xy-plane, satisfy for all .
(p3) The remaining eigenvalue of is negative for and is decreasing as a function of such that at , where , 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 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 and is the unstable manifold for . Consequently, the restriction of the map to denoted as , 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 on which the map has a period-doubling (parameter) cascade , with , to a Milnor attractor. Thus, the sketch of our proof is complete.
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 , , with a fixed point that bifurcates from a sink to a source, giving birth to an invariant homeomorph of an -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 , the following is a nice smooth example of the type of bifurcation we interested in for : Consider the smooth (actually real analytic) map defined as
where denotes the Euclidean norm on and the parameter . The origin is a fixed point for all , which is a sink for and a source when . We note that the origin is the only fixed point for , but for it is easy to see that the fixed point set of the map comprises all points on the -sphere in addition to the origin. Moreover, it is rather easy to show that is a smooth increasing function for such that and as . Furthermore, is an attractor with basin of attraction .
In the sequel, we shall almost exclusively consider such bifurcations for , where a particularly simple example of the piecewise smooth variety is
for the fixed point in a small neighborhood of . Given that each of the coordinate functions is just a logistic map, it is easy to see that is a sink for and a source for , so that is a bifurcation value for the parameter. Moreover, the system (36) has the four attracting 2-cycles for
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, and . Here
which define the vertices (37) of the boundary of the rectangular solid containing when is small and positive. For these parameter values, it is not difficult to verify that , which is homeomorphic and piecewise-linearly diffeomorphic to the 2-sphere and is an F-invariant, locally attracting set for the dynamical system having a diameter that increases with , as indicated by the simulation in Fig. 4.1. This behavior is somewhat like a piecewise-linear analog of Neimark-Sacker bifurcations in (see   ); one in fact that has several natural variations. For example, if we consider the following minor modification of (36):
where , it is not difficult to verify that this system first has a flip bifurcation at , followed by the fixed point of interest changing from a saddle to a source (at ) and generating a locally attracting, F-invariant piecewise-linear isomorph of 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 map of the form
where U is an open subset of and the parameter ranges over the interval .
Definition NST. Suppose that (40) has an isolated hyperbolic fixed point for each . Then (40) has an NST bifurcation at if the following properties obtain:
(i) There is a neighborhood such that at least one eigenvalue of at satisfies whenever and all eigenvalues of at satisfy for .
(ii) There is a locally attracting, -invariant homeomorph of the boundary of a rectangular box or 2-sphere enclosing for such that the Lebesgue measure of increases with from zero at .
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 defined as
is such that and d are sufficiently small and nonnegative, it has an NST bifurcation in the interior of the first octant such that the corresponding invariant sphere homeomorph 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 and . 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
and . It is easy to show that the coexistence fixed point of (42) is for and that the eigenvalues of are
so it follows that is a hyperbolic sink for and there is a bifurcation at across which we observe, for all permissible choices of the interaction coefficients, a center manifold at , , followed by an unstable manifold, , for , which is 1-dimensional and transverse to the plane at the coexistence fixed point. One should also note that the plane is F-invariant. Our focus is actually going to be on F2 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
from which we find that
and the higher order terms in (46) are
As a first step in completing the proof, we note that it follows from our previous analysis that there exists a positive such that when , there are parameter values for which the following property obtains: The coexistence fixed point 0 of (44) is a sink at and a source when . Moreover, there is a unique 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 , is that there are unique pairs of invariant, locally attracting plane cycles for of the form
where , , , and , and are increasing functions of that converge to zero as . Furthermore, , , , and , . With these planes, we have associated lamina defined as
such that for every , the fixed point 0 of is in the interior of the invariant rectangular box
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 . Naturally, one expects that curvilinear analogs of the above results to hold for sufficiently small, which what we shall now confirm to complete the proof. We first seek an attracting -invariant surface of the form
where is (real) analytic in a neighborhood of having the power series expansion
with an increasing function of for that converges to zero as . Invariance requires that satisfy the equation
Upon substituting (52) in the above equation, it is straightforward to inductively show that if is sufficiently small and , then
for all integers . Hence, the series (52) converges uniformly for over the specified range of the parameter and the sufficiently small value of . Of course this leads to the companion surface beneath the plane in the 2-cycle given by . Clearly, this procedure can be repeated for the first two coordinates to obtain analytic -invariant surfaces
converge uniformly for and , respectively. In addition, one has their respective analytic companion 2-cycle surfaces
Putting all of this together, we obtain the -invariant curvilinear rectangular box
Then, it follows from the definition of the map that all iterates of points in the interior of , except for the fixed point 0, converge to , as do all sufficiently close points in the exterior of . 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  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 is illustrated in Figure 3 via orbits starting near the fixed point. Observe how for all the iterates tend to the fixed point sink. Note the contrast in behavior for 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   . 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 shown for (a) ; (b) via orbits starting at three positions near the fixed point in each case.
Figure 4. Attractor for (1) with , , ; , , , , , .
Figure 5. Attractor for (1) with , , ; , , , , , .
3-dimensional Lotka-Volterra maps exhibiting steady-state behavior manifesting itself in chaotic strange attracting sets (cf.   ).
We consider maps of the form defined as
where are nonnegative and represented collectively as , all interaction coefficients are in and
Figure 6. Attractor for (1) with , , ; , , , , , .
Furthermore, we assume that implies that , so that the set is admissible inasmuch as it contains all parameter values such that maps the compact, convex set into itself. For convenience, in the sequel we denote the attracting sets and attractors (minimal attracting sets) of the map (54) by and , 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 on as increases over the interval . In this case it is easy to see that is a global attractor for and
, where , for is an attractor
with basin of attraction . Then, there is an increasing (period doubling) sequence with and with attracting set (comprising four disjoint 2-cycle attractors) with basin of attraction for . In general, for there is an attracting set (comprising disjoint -cycle attractors) with basin of attraction . Moreover, has a (Milnor) attractor that is the cartesian product of the period-doubling cascade limit attractor for each of the coordinate functions and has a dense chaotic (Milnor) attractor .
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 for defined as
It should be noted that 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
We can then determine the types of these fixed points from the spectral properties of the derivative matrix
As the above matrix is upper triangular, the eigenvalues at the fixed points are just the diagonal elements evaluated at the fixed ; namely,
which determine the types of the fixed points.
5.1.2. Dynamics in for
In this case, it follows from the above that the only fixed point in is the origin , which is a hyperbolic sink when and still a sink when . Moreover, is a global attractor in for all .
5.1.3. Dynamics in for
The dynamics in for this case is still rather easy to describe. The only fixed points are and . The eigenvalues of are
(from (60)) and , with corresponding eigenvectors along the x-, y- and z-axes, respectively, so is a hyperbolic saddle point for , with a 2-dimensional stable and 1-dimensional unstable manifold . At , the stable manifold becomes a center manifold that still attracts all points in the portion of in the xy-plane. The eigenvalues of are and , from which it readily follows that is a hyperbolic sink for all and the basin of attraction (in ) is . We note that is actually a bifurcation value if the map is considered on all of inasmuch as and are distinct for , merge at , then re-emerge for as distinct points, both of which are in .
5.1.4. Dynamics in for
For this case, things become a bit more interesting. To begin with, we have to consider two subcases: (i) ; (ii) .
(i) : For this subcase, the fixed points in are , , , and . It follows from (60) that the eigenvalues of are and , so is a hyperbolic source, actually for all . On the other hand, the eigenvalues of are
, so for all , , so for all and , so for all . Hence, is a hyperbolic saddle point for all . Next, the eigenvalues of are , so for all , , so for all and , so for all . Consequently, is also (symmetrically) a hyperbolic saddle point for all .
The eigenvalues of are and . Therefore, for and for all and
when . Hence, there is a flip bifurcation (along the z-axis) at . It also follows from the eigenvalues and an examination of the map (56) that is a hyperbolic sink for and still a sink for , but it changes into a hyperbolic saddle point for . Thus, 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 when and these points are all different for .
Observe that is also a bifurcation value since , , and are all equal when , but are distinct when . The eigenvalues of are and from which we conclude that
is a hyperbolic saddle point for all , since in this parameter range. It is worth noting that the generalized eigenspace corresponding to lies in the xy-plane, which happens to contain the stable manifold of . We also conclude from these considerations that is an attractor for all with basin of attraction given by .
(ii) : In this range, we must add the fixed points , and , noting once again that when . Moreover, , , and are distinct, with only for and all of them are distinct when . The fixed points and are essentially
symmetric and and have the eigenvalues , , and , , , respectively. Consequently, both and are hyperbolic saddle points for all and a closer examination of the map shows the dynamics in neighborhoods of these points still exhibit saddle point behavior for . The eigenvalues of are , , so it follows that is a hyperbolic sink for . In fact, an analysis of (56) shows that is an attractor whenever with .
5.1.5. Dynamics in , Including Chaotic Strange Attractors, for
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 is a hyperbolic source for this range, while and are both hyperbolic saddle points with 1-dimensional stable and 2-dimensional unstable manifolds for and at these fixed points have a stable manifold that is actually superstable. Furthermore, is a hyperbolic source for all 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 that change from hyperbolic saddle points to sources as increases. The fixed point continues to be a hyperbolic saddle point (with stable manifold ) for this parameter range. Both and 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 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
are and , so that is hyperbolic, with a
2-dimensional stable manifold and 1-dimensional unstable manifolds over the whole parameter range. Observe that decreases and increases with increasing for . Therefore, there is a flip bifurcation at with a period-doubling cascade leading to chaos along the unstable manifold as increases beyond 3 in the manner described in Section 3. Moreover, the sequence of attracting 2n-cycles in actually comprise attractors with basins of attraction for all natural numbers . These cycle attractors (periodic sinks) can also be viewed as results of tangent intersections of and in the manner of Newhouse theory .
More interesting attractors chaotic strange ones rooted in the nature of the 2-dimensional stable and 1-dimensional unstable manifolds of , respectively, occur when 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) defined as
where is as in (55), satisfies the following properties: (i) ; (ii) and sufficiently small and (iii) are also sufficiently small. Here is the maximum of the admissible parameter range , 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 of the form
with basin of attraction , 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,
with , where
Let Q be the solid prism having the following vertices in the plane: , , and and corresponding vertices in the plane : , , and
as shown in Figure 7. We note that both the upper and lower faces
of Q are mapped into the -plane ( ). Moreover, it is easy to verify that Q is a trapping set for the map for ; i.e., for all 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
where R is the rectangular bottom face of Q, which lies in the xy-plane. We compute that
which is a parabolic curve in Q that begins at
Figure 7. Depiction of generators of the GAH: the trapping set Q and for (63).
and ends at
For example, the parabolic curve begins at , which lies between and , while begins at , which lies between and for .
Then it is straightforward to verify that combining the images of all the vertical fibers in Q produces the (twisted) horseshoe shaped image shown in Figure 7. Naturally, the hyperbolic fixed point belongs to and also to
For these sets by virtue of their construction and Birkhoff-Moser-Smale theory (cf.      ), attracts all points in the and are fractal -invariant sets on which the restriction is chaotic. Consequently, the are chaotic, strange attracting sets. To complete the proof, it remains to show that for each is minimally attracting instead of a periodic sink subset generated by tangent intersections of
and in the manner of Newhouse theory . As a matter of fact, it follows that such are generalized horseshoe attractors of the type described in  if we can verify that there are no tangent intersections of and for by verifying that the -image of the arch of lies well below
the plane . We shall describe the arch of the horseshoe by finding
estimates for the maximum height (z-value) of the fiber images as ranges over R. For this, we define
and find its maximizer and maximum height , which are readily found to be and . Consequently, on R whenever . Hence, the arch of
certainly lies in and it follows that for all . Thus, since
for all , there are no tangent intersections of and for the given range of the birth-rate parameters, which means that each 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) ; (b) ; (c) ; (d) with , , , , and .
the horseshoe of the above theorem is analogous to the twisted horseshoe for the 2-dimensional Lotka-Volterra map .
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.
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.
 Blackmore, D., Chen, J., Perez, J. and Savescu, M. (2001) Dynamical Properties of Discrete Lotka-Volterra Equations. Chaos, Solitons and Fractals, 12, 2553-2568.
 Roeger, L.-I. and Lahodny, Jr.G. (2013) Dynamically Consistent Discrete Lotka-Volterra Competition Systems. Journal of Difference Equations and Applications, 19, 191-200.
 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.
 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.
 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.