Nonlinear Analysis: Real World Applications 10 (2009) 116–129
www.elsevier.com/locate/nonrwa
Chaotic dynamics of a discrete prey–predator model with
Holling type II
H.N. Agiza∗, E.M. ELabbasy, H. EL-Metwally, A.A. Elsadany
Department of Mathematics, Faculty of Science, Mansoura University, Mansoura 35516, Egypt
Received 22 August 2006; accepted 15 August 2007
Abstract
A discrete-time prey–predator model with Holling type II is investigated. For this model, the existence and stability of three
fixed points are analyzed. The bifurcation diagrams, phase portraits and Lyapunov exponents are obtained for different parameters
of the model. The fractal dimension of a strange attractor of the model was also calculated. Numerical simulations show that the
discrete model exhibits rich dynamics compared with the continuous model, which means that the present model is a chaotic, and
complex one.
c 2008 Published by Elsevier Ltd
Keywords: Prey–predator model; Holling type II functional response; Chaotic behavior; Layapunov exponents; Fractal dimension
1. Introduction
It is well-known that the Lotka–Volterra prey–predator model is one of the fundamental population models. A
predator–prey interaction has been described firstly by two pioneers Lotka (1924) [1] and Volterra (1926) [2] in two
independent works. After them, more realistic prey–predator models were introduced by Holling suggesting three
kinds of functional responses for different species to model the phenomena of predation [3]. The research dealing
with interspecific interactions has mainly focused on continuous prey–predator models of two variables, where the
dynamics include only stable equilibrium or limit cycles. Nevertheless, some works by Danca et al. [4], Jing and
Yang [5], Liu and Xiao [6] and Elabbasy et al. [7] showed that, for the discrete-time prey–predator models the
dynamics can produce a much richer set of patterns than those observed in continuous-time models. Also Summers
et al. have examined four typical discrete-time ecosystem models under the effects of periodic forcing [8]. They found
that a system which has simplistic behavior in its unforced state can assume chaotic behavior when subjected to
periodic forcing, dependent on the values chosen for the controlling parameters; such a phenomenon is well-known
in the physical sciences in the theory of nonlinear oscillators see [8]. Danca et al. [4] demonstrated that, the chaotic
dynamics in a simple discrete-time prey–predator model with Holling I take place [4]. In previous work, we modified
Danca et al. model [4] and studied the complex dynamics [7], and found that the modified model is more realistic than
∗ Corresponding author. Tel.: +20 105160714; fax: +20 502246781.
E-mail addresses: agizah@mans.edu.eg (H.N. Agiza), emelabbasy@mans.edu.eg (E.M. ELabbasy), helmetwally@mans.edu.eg
(H. EL-Metwally), aelsadany1@yahoo.com (A.A. Elsadany).
1468-1218/$ - see front matter c 2008 Published by Elsevier Ltd
doi:10.1016/j.nonrwa.2007.08.029
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 117
that of Danca et al. [4]. However, the same conclusion was obtained using another technique, in which Euler method
was used by Jing and Yang [5] and Liu and Xiao [6].
One of the real life applications of the prey–predator model is the Lynx and its prey the snowshoe Hare study
documented by the Hudson Bay company for the time interval 1845–1935 [9]. The extension of discrete prey–predator
model to cover the Holling type II had a little attention in the discrete case till now, due to its complexities. Therefore,
the present work aims to shed more light on this subject through analyzing the dynamic complexities in a discrete-time
prey–predator model with the Hollings type II functional response. That is, we shall focus our attention on analyzing
how the Holling type II response [3] affects the dynamic complexities of prey–predator interactions.
This paper is organized as follows: in Section 2, the discrete prey–predator model with Holling type II is formulated,
then the existence and stability of three fixed points are derived. In Section 3, some values of the parameters, such
that the model undergoes the flip bifurcation and the Hopf bifurcation in the interior R2
+, were derived and also
discussed. The numerical simulation of the analytic results, such as the bifurcation diagrams, strange attractors,
Lyapunov exponents and fractal dimension were presented in Section 4. Finally, Section 5 draws the conclusion.
2. Model
The classical prey–predator system always be in the following form
x (t) = xq(x) − αyp(x)
y (t) = (p(x) − β)y
x(0), y(0) > 0,
(1)
where x, y represent the prey and predator density, respectively. p(x) is the so-called predator functional response
and α, β > 0 are the conversion and predator’s death rates, respectively. If p(x) = mx
1+εx , q(x) = ax(1 − x), then Eq.
(1) becomes the following well-known prey–predator model with the Holling type II functional response [10]:



x (t) = ax(1 − x) − α
mxy
1 + εx
y (t) =
mx
1 + εx
− β y,
(2)
where a, m and ε are the positive parameters that stand for prey intrinsic growth parameter, half saturation parameter,
limitation of the growth velocity of the predator population with increase in the number of prey, respectively. The
above model (Eq. (2)) has been studied by many authors [11–13] and it was shown that, the dynamics include only
stable equilibrium or limit cycles.
Another possible way to understand the complex problem of competition between two interacting species is by
using discrete models [4]. In the present work we study the dynamics of discrete prey–predator model with Holling
type II which has the following two difference equations:
T :



xn+1 = axn(1 − xn) −
bxn yn
1 + εxn
yn+1 =
dxn yn
1 + εxn
,
(3)
where a, b, c and d are the nonnegative parameters. The map given by Eq. (3) is a noninvertible map of the plane.
The study of the dynamical properties of the above map allows us to have information about the long-run behavior
of prey–predator populations. Starting from given initial condition (x0, y0), the iteration of (3) uniquely determines a
trajectory of the states of population output in the following form
(x(n), y(n)) = T n
(x0, y0),
where n = 0, 1, 2, . . . .
3. The fixed points and their stability
In this section, we first determine the existence of the fixed points of map (3), then investigate their stability by
calculating the eigenvalues for the variational matrix of (3) at each fixed point. To determine the fixed points we have
118 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
to solve the nonlinear system given by



x = ax(1 − x) −
bxy
1 + εx
y =
dxy
1 + εx
.
By simple computation of the above algebraic system it was found that there are three nonnegative fixed points,
(i) E0(0, 0) is origin,
(ii) E1(a−1
a , 0) is the axial fixed point in the absence of predator (y = 0) exist for a > 1, and
(iii) E2(x∗, y∗) is the interior fixed point, where
x∗
=
1
d − ε
y∗
=
d
d − ε
a
b
1 −
1
d − ε
−
1
b
(4)
exists if and only if the following condition is satisfied:
d > ε +
a
a − 1
.
3.1. The dynamic behavior of the model
Now, we investigate the local behavior of the model (3) around each of the above fixed points. The local stability
analysis of the model (3) can be studied by computing the variation matrix corresponding to each fixed point. The
Jacobian matrix of Eq. (3) at the state variable is given by
J(x, y) =



a(1 − 2x) −
by
(1 + εx)2
−
bx
1 + εx
dy
(1 + εx)2
dx
1 + εx


 . (5)
The characteristic equation of Jacobian matrix can be written as
λ2
− Tr λ + Det = 0, (6)
where Tr is the trace and Det is the determinant of the Jacobian matrix J(x, y) which is defined as Tr = a(1 − 2x) −
by
(1+εx)2 + dx
1+εx and Det = adx(1−2x)
1+εx . Hence the model (3) is a dissipative dynamical system if
adx(1 − 2x)
1 + εx
< 1,
conservative dynamical one, if and only if
adx(1 − 2x)
1 + εx
= 1,
and is an undissipated dynamical system otherwise.
In order to study the stability of the fixed points of the model, we first give the following lemma, which can be
easily proved by the relations between roots and coefficients of a quadratic equation [6].
Lemma 1. Let F(λ) = λ2 − Bλ + C. Suppose that F(1) > 0, λ1 and λ2 are the two roots of F(λ) = 0. Then
(i) |λ1| < 1 and |λ2| < 1 if and only if F(−1) > 0 and C < 1;
(ii) |λ1| < 1 and |λ2| > 1 (or |λ1| > 1 and |λ2| < 1) if and only if F(−1) < 0;
(iii) |λ1| > 1 and |λ2| > 1 if and only if F(−1) > 0 and C > 1;
(iv) λ1 = −1 and λ2 = 1 if and only if F(−1) = 0 and B = 0, 2;
(v) λ1 and λ2 are complex and |λ1| = |λ2| if and only if B2 − 4C < 0 and C = 1.
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 119
Let λ1 and λ2 be the two roots of Eq. (6), which are called eigenvalues of the fixed point (x, y). We recall some
definitions of topological types for a fixed point (x, y). A fixed point (x, y) is called a sink if |λ1| < 1 and |λ2| < 1,
so the sink is locally asymptotically stable. (x, y) is called a source if |λ1| > 1 and |λ2| > 1, so the source is
locally unstable. (x, y) is called a saddle if |λ1| > 1 and |λ2| < 1 (or |λ1| < 1 and |λ2| > 1). And (x, y) is called
non-hyperbolic if either |λ1| = 1 or |λ2| = 1.
Proposition 2. The fixed point E0 is a sink if a < 1; saddle if a > 1 and non-hyperbolic if a = 1.
Proof. One can see that the Jacobian matrix J at E0 is given by
J(E0) =
a 0
0 0
.
Hence the eigenvalues of matrix are λ1 = a and λ2 = 0. Thus it is clear that E0 is a sink point if a < 1; a saddle point
if a > 1 and non-hyperbolic point if a = 1.
Proposition 3. If a > 1 there are at least four different topological types of E1(a−1
a , 0) for all permissible values of
parameters:
(i) E1 is a sink if 1 < a < 3 and d < a+ε(a−1)
a−1 ;
(ii) E1 is a source if a > 3 and d > a+ε(a−1)
a−1 ;
(iii) E1 is non-hyperbolic if a = 3 or d = a+ε(a−1)
a−1 ;
(iv) E1 is a saddle for the other values of parameters except those values in (i)–(iii).
Proof. In order to prove this result, we calculate the eigenvalues of Jacobian matrix for E1 which is given by
J(E1) =



2 − a −
b(1 − a)
a + ε(a − 1)
0
d(a − 1)
a + ε(a − 1)


 .
The eigenvalues of the matrix, J(E1), are λ1 = 2 − a and λ2 = d(a−1)
a+ε(a−1) . By using Lemma 1 it is easy to see that,
E1 is a sink if 1 < a < 3 and d < a+ε(a−1)
a−1 ; E1 is a source if a > 3 and d > a+ε(a−1)
a−1 ; E1 is non-hyperbolic if a = 3
or d = a+ε(a−1)
a−1 and E1 is a saddle point for the other values of parameters.
From the condition (iii) of Proposition 3, it is easy to see that one of the eigenvalues of the fixed point E1(a−1
a , 0) is
−1 and the other eigenvalue is neither 1 nor −1 and then it also implies that all the parameters locate in the following
set:
SE1 = (a, b, d, ε) : d =
a + ε(a − 1)
(1 − a)
, a = 3, a > 1, b, d, ε > 0 .
The fixed point E1 can pass through flip bifurcation when parameters vary in the small neighborhood of SE1 and a
center manifold of Eq. (3) at E1(a−1
a , 0) is y = 0 which is restricted to the center manifold of logistic model [14]. In
this case, the predator becomes extinct and the prey pass through the period-doubling bifurcation to chaos in the sense
of Li–Yorke by the varying bifurcation parameter a.
3.2. Local stability and dynamic behavior around interior fixed point E2
Now, we shall discuss the stability and bifurcations of interior fixed point E2. The fixed point E2 is stable if it
satisfies the following conditions:



1 + Tr (J(E2)) + Det(J(E2)) > 0
1 − Tr (J(E2)) + Det(J(E2)) > 0
1 − Det(J(E2)) > 0,
(7)
120 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
where J(E2) is given by
J(E2) =




a 1 −
2
d − ε
−
d − ε
d
a 1 −
1
d − ε
− 1
−b
d
d − ε
b
a 1 −
1
d − ε
− 1 1



 . (8)
The above conditions are well-known [15] and are sufficient for the local stability of the fixed point and are necessary
for the two roots of the characteristic equation, λ1,2, to be inside the unit circle of the complex plane, since
Tr (J(E2)) = 1 + a 1 −
2
d − ε
−
d − ε
d
a 1 −
1
d − ε
− 1
Det(J(E2)) = a 1 −
2
d − ε
. (9)
By using Lemma 1, it is easy to deduce the following proposition.
Proposition 4. When d > ε and a > d−ε
d−ε−1 , the model (3) has a unique positive fixed point E2(x∗, y∗) and:
(i) it is a sink if a > (3d−ε)(d−ε)
(d−ε)(2d−ε)−2d(d−ε−2) and a < d−ε
d−ε−2 .
(ii) it is a source if a < (3d−ε)(d−ε)
(d−ε)(2d−ε)−2d(d−ε−2) and a > d−ε
d−ε−2 .
(iii) it is non-hyperbolic if a = (3d−ε)(d−ε)
(d−ε)(2d−ε)−2d(d−ε−2) .
(iv) it is a saddle if a < (3d−ε)(d−ε)
(d−ε)(2d−ε)−2d(d−ε−2) .
From Lemma 1, one can see that one of the eigenvalues of the positive fixed point E2 is −1 and the other is neither
1 nor −1 if (iii) of Proposition 4 holds. Rewriting the condition (iii) of Proposition 4 one get the following set:
SE2 = (a, b, d, ε) : a =
(3d − ε)(d − ε)
(d − ε)(2d − ε) − 2d(d − ε − 2)
, d > ε, a >
d − ε
d − ε − 1
.
By solving the characteristic equation (6) at the interior fixed point E2, the roots (eigenvalues at E2) will be:
λ1,2 =
Tr (J(E2)) ± (Tr (J(E2)))2 − 4Det(J(E2))
2
, (10)
where Tr (J(E2)) and Det(J(E2)) are given by Eq. (9). The eigenvalues in numerical simulation can be used to
classify different types of bifurcations.
4. Numerical simulations
To provide some numerical evidence for the qualitative dynamic behavior of the model (3), the phase portraits,
bifurcation diagrams, Lyapunov exponents, sensitive dependence on initial conditions and fractal dimension were
used to illustrate the above analytical results and for finding new dynamics as the parameters vary.
Since the dynamics of discrete prey–predator model has been examined [4–7], we will now mainly focus our
attention on the effect of Holling type II functional response. First we will study the phase portrait of the model (3)
when we change only the parameter ε and fix the others. To study the behavior of the model (3) when the parameter ε
varied in the interval ]0.22, 1[, one can consider the initial condition (x0, y0) situated in the basin of attraction of fixed
point E2. When the control parameter ε varies, the stability of a periodic solution may be lost through various types
of bifurcations.
Without loss of generality:
We fix the parameters a = 4.1, b = 3, d = 3.5 and assume that ε varies. The phase portraits are considered in the
following cases:
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 121
Fig. 1. An attractor fixed point for system (3) which exists for ε = 0.99.
Fig. 2. A stable fixed point for system (3) which exists for ε = 0.9.
An attractor fixed point takes place for ε = 0.99, which means that the system orbit is a fixed point, as shown in
Fig. 1.
Fig. 2 shows that the fixed point E2 is a stable attractor at ε = 0.9. For this parameter value, the fixed point
E2 occurs at x∗ = 0.3846153846, y∗ = 0.6834319527 and the associated complex conjugate eigenvalues are
λ± = 0.4073626372 ± i0.8832947005, then |λ±| = 0.9727043981, and this means that the fixed point E2 is
asymptotically stable.
The behavior of model (3) before a Hopf bifurcation at ε = 0.856 is shown in Fig. 3. Fig. 4 demonstrates the
behavior of the model after a Hopf bifurcation when ε = 0.854. From Figs. 3 and 4, we deduce that the fixed point
E2 loses its stability through a Hopf bifurcation, when the parameter ε varies from 0.856 to 0.854.
Decreasing the control parameter ε forward, (ε = 0.845) leads to the fixed point E2 which has (0.3766478343,
0.6836288233) and the associated eigenvalues are λ± = 0.4156724510 ± i0.9158079345. Since the modulus of the
complex conjugate eigenvalues is |λ±| = 1.005727478, one can conclude that the fixed point became unstable and an
invariant closed curve was created around the fixed point. Fig. 5 demonstrates and confirms the above argument.
122 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
Fig. 3. Phase portrait for system (3) before a Hopf bifurcation which exists for ε = 0.856.
Fig. 4. Phase portrait for system (3) after a Hopf bifurcation which exists for ε = 0.854.
As ε is further decreased, however, the phase portrait starts to fold – interrupted by periodic and then by windows
– a quasi-periodic behavior which changes to a chaotic one. We can see that the circle, after being stretched, shrunk
and folded creates new phenomena due to the breakdown of the invariant closed curve, Fig. 6.
For decreasing ε we obtain the multiple invariant closed curves brought about by a Hopf bifurcation of iterates of
the model (3). In these cases higher bifurcations of the tours occurs as the system moves out of the quasi-periodic
region by decreasing ε. The dynamics move from one cycle to another by a periodic way, but the dynamics in each
closed curve, may be periodic or quasi-periodic. Fig. 7 represents the set of 16 closed curves brought about by a Hopf
bifurcation of the 16th iterate of the model (3), obtained for ε = 0.54. The radius of quasi-periodic solution grows as
ε is further decreased as shown in Fig. 8.
Moreover, these closed curves may break leading to multiple fractal tori on which the dynamics are chaotic. When
0.52 ≤ ε ≤ 0.22, a chaotic attractor appears in the phase portrait and the chaotic orbits evolves upon ε decreasing.
Figs. 9 and 10 represent strange attractors for the model (3) with ε = 0.52 and ε = 0.4, respectively which exhibit
fractal structure.
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 123
Fig. 5. The invariant closed curve for system (3) around the fixed point created after the bifurcation which exists for ε = 0.845.
Fig. 6. The breakdown of the invariant closed curve for system (3) which exists for ε = 0.6.
The strange attractor is produced by the breaking of the invariant circles which results in the appearance of the
sixteen chaotic regions which change as they are linked to a single chaotic attractor, as shown clearly in Figs. 8 and 9.
In addition, in system (3) full developed chaos occurs when ε = 0.22, Fig. 11.
In Fig. 12 the bifurcation diagram for system (3) is plotted as a function of the control parameter ε for 0.22 ≤ ε ≤ 1.
From this figure it is clear that the fixed point is stable for ε > 0.85483871, and loses its stability at a Hopf bifurcation
parameter value ε = 0.85483871. As ε decreases the behavior of this model becomes very complicated, including
many chaotic bands, the Hopf bifurcation, tangent bifurcations and crises. The bifurcation diagram for system (3)
versus the control parameter a is drawn in Fig. 13, while other parameters are kept fixed as follows b = 3, d = 3.5
and ε = 0.4. The nature of this stable coexistence of the fixed changes as a Hopf bifurcation takes place at a ≈ 3,
initiating the rise of complicated dynamics which are not easy to be characterized in the case is illustrated in Fig. 13.
Figs. 14 and 15 represent the dominant Lyapunov exponent with respect to ε and a, respectively. The maximal
Lyapunov exponent corresponding to Fig. 12 are given in Fig. 14, which shows the complete inverse a Hopf bifurcation
exist at ε = 0.85483871 and the existences of chaotic regions as parameter ε decreases. Also, the maximal Lyapunov
124 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
Fig. 7. The existence of multiple invariant closed curves for system (3) which exists for ε = 0.54.
Fig. 8. The existence of multiple invariant closed curves for system (3) which exists for ε = 0.53.
exponent corresponding to Fig. 13 are computed and plotted in Fig. 15 in which we can easily see that the maximal
Lyapunov exponents have negative values for a (1, 3), that is, to say that the stable region is bigger than the chaotic
one (3, 4.1). In general, when the maximal Lyapunov exponent is positive, this can be considered as an evidence for
chaos.
4.1. Sensitive dependence on initial conditions
To demonstrate the sensitivity to initial conditions of the model (3), two orbits with initial points (x0 , y0) and
(x0 + 0.0001, y0) are computed, respectively, and are represented in Fig. 16. At the beginning, the two time series are
overlapped and are indistinguishable; but after a number of iterations, the difference between them builds up rapidly.
Fig. 16 shows a sensitive dependence on initial conditions for x-coordinate of the two orbits for the model (3), which
is plotted against the time with the parameters values (a, b, d, ε) = (4.1, 3, 3.5, 0.4). The difference between the two
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 125
Fig. 9. Chaotic attractor for system (3) which exists for ε = 0.52.
Fig. 10. Chaotic attractor for system (3) which exists for ε = 0.4.
x-coordinates is 0.001, while the other coordinate is kept to have the same value. For this case the two orbits with
initial points (0.2, 0.1) and (0.201, 0.1) are computed respectively and plotted in Fig. 16 as a function of time.
Also, sensitive dependence on initial conditions, y-coordinates of the two orbits, for model (3), are plotted against
the time with the parameter constellation (a, b, d, ε) = (4.1, 3, 3.5, 0.4); in Fig. 17. The y-coordinate of initial
conditions differ by 0.001, while the other coordinate is kept at the same value. In Figs. 16 and 17, it is shown that the
time series of the system (3) sensitively depends on the initial conditions, i.e. complex dynamic behavior occur in this
model.
4.2. Fractal dimension of the map (3)
Strange attractors are typically characterized by fractal dimensions. We examine the important characteristics
of neighboring chaotic orbits to see how rapidly they separate each other. This separation is quantified by using
the concept of Lyapunov exponents [16]. These exponents represent a dynamic measure of chaos that average the
126 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
Fig. 11. Full developed chaos for system (3) which exists for ε = 0.22.
Fig. 12. Bifurcation diagram for system (3) versus ε.
separation of the orbits of nearby initial conditions as system moves forward in time. The Lyapunov dimension
[17,18] is defined by using Lyapunov exponents as follows:
dL = j +
i= j
i=1
Λi
|Λj |
,
with Λ1, Λ2, . . . ., Λn, where j is the largest integer such that
i= j
i=1 Λi ≥ 0 and
i= j+1
i=1 Λi < 0.
Our model is a two-dimensional map which has the Lyapunov dimension in the form
dL = 1 +
Λ1
|Λ2|
, Λ1 > 0 > Λ2.
By the definition of the Lyapunov exponents and with help of the computer simulation one can calculate the Lyapunov
dimension of the strange attractor of system (3). Each one of the chaotic attractors plotted in Figs. 9 and 10 has positive
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 127
Fig. 13. Bifurcation diagram for system (3) versus a.
Fig. 14. Maximal Lyapunov exponent versus ε corresponding to Fig. 12.
Lyapunov exponent. Two Lyapunov exponents are estimated and are found to be Λ1 ≈ 0.0160 and Λ2 ≈ −0.2061 for
ε = 0.52 and Λ1 ≈ 0.10601 and Λ2 ≈ −0.2071 for ε = 0.4. The above results show that the present model (3) has a
fractal dimension in the following form
dL ≈ 1 +
0.0160
0.2061
= 1.0776, for ε = 0.52 and
dL ≈ 1 +
0.10601
0.2071
= 1.5119, for ε = 0.4.
(11)
From the strange attractors presented in Figs. 9 and 10 and corresponding fractal dimension which was calculated in
Eq. (11), one can deduce that, the discrete prey–predator model with Holling type II is of a very complicated behavior.
128 H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129
Fig. 15. Maximal Lyapunov exponent versus a corresponding to Fig. 13.
Fig. 16. Sensitive dependence on initial conditions of the model, x-coordinates of the two orbits, plotted against time; the x-coordinates of initial
conditions differ by 0.001, and the other coordinate kept equal.
5. Conclusions
In this paper, chaotic dynamics and bifurcation of a nonlinear discrete-time prey–predator model with second
Holling type have been investigated. The stability of fixed points and bifurcations are investigated. Basic properties
of the model have been analyzed by means of phase portrait, bifurcation diagrams, Lyapunov exponents and sensitive
dependence on initial conditions. Under certain parametric conditions, the interior fixed point enters a Hopf bifurcation
phenomenon. The suggested model has more rich features and more complicated dynamics than that exist in the
continuous case. Nevertheless, identifying complicated, possibly chaotic dynamics in population data has remained a
major challenge in ecological studies [19]. This could be very useful for the biologists who work with discrete-time
prey–predator models.
H.N. Agiza et al. / Nonlinear Analysis: Real World Applications 10 (2009) 116–129 129
Fig. 17. Sensitive dependence on initial conditions of the model, y-coordinates of the two orbits, plotted against time; the y-coordinates of initial
conditions differ by 0.001, and the other coordinate kept equal.
Acknowledgments
We thank referees for several useful suggestions that have considerably improved the paper.
References
[1] V. Voltera, Opere matematiche: mmemorie e note, vol. V, Acc. Naz. dei Lincei, Roma, Cremon, 1962.
[2] A.J. Lotka, Elements of Mathematical Biology, Dover, New York, 1962.
[3] C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Ent. Soc. Canada
45 (1965) 1–60.
[4] M. Danca, S. Codreanu, B. Bako, Detailed analysis of a nonlinear prey-predator model, J. Biol. Phys. 23 (1997) 11–20.
[5] Z.J. Jing, J. Yang, Bifurcation and chaos discrete-time predator–prey system, Chaos Solitons Fractals 27 (2006) 259–277.
[6] X. Liu, D. Xiao, Complex dynamic behaviors of a discrete-time predator–prey system, Chaos Solitons Fractals 32 (2007) 80–94.
[7] E.M. ELabbasy, H.N. Agiza, H. EL-Metwally, A.A. Elsadany, Chaos and bifurcation of a nonlinear discrete-time prey-predator system, World
J. Model. Simul. (submitted for publication).
[8] D. Summers, C. Justian, H. Brian, Chaos in periodically forced discrete-time ecosystem models, Chaos Solitons Fractals 11 (2000) 2331–2342.
[9] L. Edelstein-Keshet, Mathematical Models in Biology, Random House, New York, 1988.
[10] S.B. Hsu, T.W. Hwang, Global stability for a class predator–prey system, SIAM J. Appl. Math. 55 (1995) 763–783.
[11] N. Kasarinoff, P. van der Deiesch, A model of predator–prey system with functional response, Math. Biosci. 39 (1978) 124–134.
[12] S. Ruan, D. Xiao, Global analysis in a predator–prey system with nonmonotonic functional response, SIAM J. Appl. Math. 61 (2001)
1445–1472.
[13] J. Hainzl, Stability and Hopf bifurcation a predator–prey system with several parameters, SIAM J. Appl. Math. 48 (1988) 170–180.
[14] J.D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1998.
[15] S.N. Elaydi, An Introduction to Difference Equations, Springer-Verlag Publishers, 1996.
[16] V.I. Oseledec, A multiplicative ergodic theorem: Lyapunov characteristic numbers for dynamical systems, Trans. Moscow Math. Soc. 19
(1968) 197–231.
[17] J.L. Kaplan, Y.A. Yorke, A regime observed in a fluid flow model of Lorenz, Comm. Math. Phys. 67 (1979) 93–108.
[18] J.H.E. Cartwright, Nonlinear stiffness, Lyapunov exponents, and attractor dimension, Phys. Lett. A 264 (1999) 298–304.
[19] Y. Xiao, D. Cheng, S. Tang, Dynamic complexities in predator-prey ecosystem models with age-structure for predator, Chaos Solitons Fractals
14 (2002) 1403–1411.