#### Abstract

A Leslie population model for two generations is investigated by qualitative analysis and numerical simulation. For the different parameters *a* and *b* in the model, the dynamics of the system are studied, respectively. It shows many complex dynamic behavior, including several types of bifurcations leading to chaos, such as period-doubling bifurcations and Neimark–Sacker bifurcations. With the change of parameters, attractor crises and chaotic bands with periodic windows appear. The largest Lyapunov exponents are numerically computed and can verify the rationality of the theoretical analysis.

#### 1. Introduction

In 1964, Hénon took KAM theorem as the background and found that there was deterministic random behavior in two-dimensional nonintegrable Hamiltonian system, that is claimed Hénon attractor. Ruelle and Takens proposed the concept of “strange attractor,” which promoted the research of Smale horseshoe attractor [1]. Lorenz [2] pointed out that there must be a connection between the inexact recurrence of climate and the inability of long-term weather forecast. It is found that the chaos phenomenon is “extremely sensitive to the initial conditions.” The word “chaos” is formally used by Li and Yorke [3], which is considered as the first formal expression of chaos theory. In [1], an example of horseshoe mapping is given, which opened the mathematical method of studying chaos. Lasota [4] studied the initial value problem of the first order nonlinear partial differential equation of the wormhole model. Brunovsky [5] gave the definition of chaotic mapping.

A single species continuous-time population model is studied by [6] and later studied further by Pearl and Reed [7].

May [8] studied several discrete-time models with a single-species displaying chaotic behavior for certain parameter values. Beddington et al. [9] studied that discrete-time host, parasitoid models show complex dynamic behavior.

Leslie [10, 11] introduced an age-structured linear population model. Later, this kind of discrete-time models contain linear or nonlinearity forms, which are called Leslie models. Ugarcovici and Weiss [12] studied the dynamical behavior of Ricker model. This model is described by the two-dimensional mapping :where *x* and *y* stand for the density of the first age group and the second age group. *a* and *γa* are the group’s initial fertility rates (*a*, *γ* > 0), *b* is the survival rate from the ﬁrst age group to the second one, and *λ* is the decay index, *λ* > 0. In equation (1), the fertility rate monotonically decreases as a function of the total population size, and the fertility decay is exponential. The other mold is the Hassell model. It is described by the two-dimensional mapping :where *a* and *γa* are the group’s initial fertility rates, *b* is the survival rate from the ﬁrst age group to the second one, and *β* is the decay index, *β* > 1. In equation (2), the fertility rate monotonically decreases as a function of the total population size and the fertility decay is polynomial. Guo [13] studied the chaos control of the population model.

For some parameter values, these models admit an ergodic attractor which supports a unique physical probability measure. This physical measure satisﬁes in the strongest possible sense the population biologist’s requirement for ergodicity in their population models. Wikan and Mjølhus [14] and Ugarcovici and Weiss [12] showed that Ricker mapping and Hassell mapping produce Hénon-like chaotic attractors.

Here, for the different parameters *a* and *b* in the model, the dynamics of the system is studied, respectively. It shows many complex dynamic behavior, including several types of bifurcations leading to chaos, such as period-doubling bifurcations and Neimark–Sacker bifurcations. With the change of parameters, attractor crises and chaotic bands with periodic windows appear. The largest Lyapunov exponents are numerically computed and can verify the rationality of the theoretical analysis. Zhu and Zhao [15] studied the dynamic complexity of a host, parasitoid ecological model with the Hassell growth function for the host. Liu and Xiao [16] studied the complex dynamic behaviors of a discrete-time predator-prey system.

In biology or ecology, the complex chaotic behavior of this mapping shows the relationship among the number, birth rate, and survival rate in a population, whether it survives in a balanced state or makes the population develop in disorder or chaos. This research can provide theoretical basis and help for the research in biology or ecology. Through the study of these complex dynamic properties of the population model, it can guide people to make corresponding meaningful work to keep the ecological balance. For example, it is applied in the marine fishing, or in the reproduction and population growth of a certain species in nature.

#### 2. The Nonlinear Hassell Population Model

Hassell model is described by the two-dimensional mapping :wherewhere are the density of the first age group and the second one. *a* and *γa* are the group’s initial fertility rates, *b* is the survival rate from the ﬁrst age group to the second one, and *β* is the decay index, *β* > 1. In equation (3), the fertility rate monotonically decreases as a function of the total population size and the fertility decay is polynomial.

#### 3. Stability Analysis of Hassell Population Model

The equation (3) has one equilibrium point: , where are positive and satisfy

The Taylor series expansion of equation (3) at the equilibrium point: is written aswhere

Let the matrix bewhere

The characteristic equation of equation (9) is

Equation (11) can be rewritten aswhere

The roots of equation (12) are

Both eigenvalues are real numbers and < 1, if

Equation (15) can be solved as

The eigenvalues become complex numbers and are inside the unit circle in the complex -plane for

That is,

The positive equilibrium point is stable under conditions (17) and (18).

#### 4. Bifurcation Analysis

Equation (3) cannot be solved analytically, and therefore its long-term behavior must be studied by numerical simulation. Figure 1 shows the bifurcation diagrams of equation (3) for the density of the first age group *x* with , and *b* *=* 0.7, as the parameter *a* increases. Figure 1 is magnified as Figure 2.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

When the initial fertility rate *a* changes between 0 and 60, equation (3) shows complicated features. For 0 < *a* < 0.75, (0*,* 0) is a global attractor. There exists a positive fixed point that is asymptotically stable for 0.75 < *a* < 15.8. At *a* = 15.8, it appears a cascade of period-3 orbit. For *a* between 15.8 and 24, it exists a period-3 window which is embedded in the strange attractor for a long time. Then, each of the period-3 orbits begins to undergo a flip bifurcation, leading to chaos. The phase portraits of various *a* corresponding to Figure 2(a) are plotted in Figures 3(a)‐3(b). Tianyan-Li and Yorke [3] proved that if a system has period-3 point, it has all periodic points. This is also verified in the following analysis as *a* increases.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

When *a* passes through the range (24.4, 27), the bifurcation diagram in Figure 2(b) shows that this window is not a periodic window with a cascade of periodic attractors, but that it includes more complex dynamic patterns. It appears that several attractors coexist in this region: period-3 and period-5 attractors. There is a very wide chaotic band with period-6 orbit meanwhile (see Figure 3(b)).

When *a* = 27.2, the chaotic band suddenly disappears in a crisis and the system enters a periodic window with a cascade of period-doubling bifurcations leading to a chaotic attractor with periodic windows. It appears period-11 and starts to appear period-9 when *a* = 27.24. The phase portraits which are associated with Figure 2(b) are disposed in Figures 3(b)‐3(c). When *a* passes through the range (27.24, 30), the detailed bifurcation diagram in Figure 2(c). As *a* > 28.26, there is a cascade of period-doubling bifurcation leading to a wide chaotic region, characterized by tangent bifurcation and attractor crises. A very narrow periodic window appears again when *a* = 29. When *a* increases beyond 28.5, a chaotic attractor abruptly appears and the periodic attractor disappears. The chaotic band has changed into four chaotic bands. This kind of transilience of the chaotic state caused by the continuous change of parameters is called chaos crisis by Grebogi et al. [17]. Because the unstable periodic orbit meets the secondary chaotic band, the orbits in the chaotic band is filled between all levels of orbits, leading to the emergence of chaos crisis. The phase portraits of various *a* corresponding to Figure 2(c) are plotted in Figures 3(d)–3(f). When *a* passes through the range (30.8, 31.7) in Figure 2(d), a large periodic window appears; then, chaos appears because of period-doubling bifurcation. When *a* = 49, it appears as period-4, and it appears as period-9 at *a* = 51.2, at last, equation (6) appears as a chaotic attractor. The phase portraits which are associated with Figure 2(e) are disposed in Figures 3(g) and 3(h).

From the above analysis, it can be seen that the two-dimensional nonlinear mapping is a process of intermittently breaking up in accordance with periodic behavior and chaos phenomenon so that the system appears chaotic motion state, which is called Pomeau–Manneville path (through intermittence chaos) by Eckmann [18]. In the process of parameter change, when the periodic window appears, the one of the roots of equation (11) is −1. The unstable period-doubling bifurcation happens immediately, and the periodic point becomes the periodic saddle point. Once the mapping point falls near the unstable periodic point, it will leave along the unstable manifold, and chaos will appear. Equation (3) also contains Feigenbaum path (through fork bifurcation) to chaos.

Figures 4(a) and 4(b) show the bifurcation diagrams of equation (3) for the density of the first age group *x* and the second one *y* with , and *a* *=* 38, as the parameter *b* changes between 0 and 0.7. Figure 4(c) is the graph magnified version of Figure 4(a).

**(a)**

**(b)**

**(c)**

In Figure 4(a), as *b* increases from 0 to 0.03, a narrow chaotic band appears. When *b* = 0.03, the chaotic regime suddenly disappears in a crisis and the system appears a period-2 orbit for *b* ∈ (0.03, 0.435), and a stable fixed point for *b* ∈ (0.435, 0.5526). A flip bifurcation occurs at *b* = 0.435. Figure 5(a) presents the phase portrait of period-2. As *b* ∈ (0.5526, 0.61) in Figure 4(c), equation (3) goes through a quasi-periodic region (including tangent bifurcation, narrow and wide periodic windows, and frequency lockings which appear as a collapse of the invariant circle to a periodic orbit). We can see that there is a stable fixed point for *b* < 0.5526, and the fixed point loses its stability as *b* increases. A Neimark–Sacker bifurcation occurs at *b* = 0.5526 and an attracting invariant cycle bifurcates from the fixed point. In this case, equation (11) has two conjugate complex roots and , satisfying = 1.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

As *b* is increasing, the smooth invariant loop becomes larger and no longer smooth. At last, the invariant loop breaks. As *b* = 0.61, there is a cascade of period-doubling bifurcation leading to a wide chaotic region, characterized by tangent bifurcation and attractor crises. Furthermore, when *b* ∈ (0.565, 0.57) and *b* ∈ (0.615, 0.619), we can observe the period-5, 8 windows within the chaotic regions, respectively. The phase portraits of various *b* corresponding to Figure 4(c) are plotted in Figures 5(b)–5(f).

According to the above analysis, when *b* is increasing, chaos appears by the Ruelle Takens Newhouse scheme (through Neimark–Sacker bifurcation) and the Pomeau–Manneville path (through intermittence chaos). In the larger region of the parameter space, the path leading to chaos is related to the Neimark–Sacker bifurcation. In these ways, the locking phase and quasiperiodic motion can be observed.

#### 5. The Largest Lyapunov Exponent

In this section, the largest Lyapunov exponent, which has proven to be the most useful dynamic diagnostic tool for chaotic systems, is considered. This quantity represents the average exponential rate of divergence or convergence of nearby orbits in phase space [19]. For a chaotic attractor, the largest Lyapunov exponent LE_{max} must be positive. If LE_{max} is negative, this implies a stable state or a periodic attractor.

The largest Lyapunov exponents corresponding to the cases shown in Figures 1 and 4 have been calculated and plotted in Figures 6(a) and 6(b), respectively. The existence of chaotic regions in the parametric space is clearly visible in these figures. The largest Lyapunov exponents for the strange attractors were found to be LE_{max} = 20.02,31.81, corresponding to the parameters *a* and *b*, respectively.

**(a)**

**(b)**

In Figure 6(a), we can easily see that the maximum Lyapunov exponents are negative for *a*(0, 25.5), that is to say, equation (3) has no chaotic region. For *a*(25.5, 35), some Lyapunov exponents are bigger than 0, some are smaller than 0, so there exists stable fixed point or stable periodic windows in the chaotic region. In Figure 6(b), the maximum Lyapunov exponents are positive for *b*(0, 0.03), that is to say, equation (3) appears in a chaotic region. For *b*(0.03, 0.526), the maximum Lyapunov exponents are negative, equation (3) appears periodic. For *b*(0.526, 0.7), some Lyapunov exponents are bigger than 0, some are smaller than 0, so there exist periodic windows in the chaotic region.

We find that the values of the maximum Lyapunov exponents are consistent with the dynamic properties of the system under different parameters.

#### 6. Conclusion

A Leslie population model for two generations is investigated by qualitative analysis and numerical simulation. For the different parameters *a* and *b* in the model, the dynamics of the system is studied, respectively. It shows many complex dynamic behavior, including several types of bifurcations leading to chaos, such as period-doubling bifurcations and Neimark–Sacker bifurcations. With the change of parameters, attractor crises and chaotic bands with periodic windows appear. The largest Lyapunov exponents are numerically computed and can verify the rationality of the theoretical analysis. Numerical simulation results not only show the consistence with the theoretical analysis but also display the new and interesting dynamical behaviors, including different periodic orbits in chaotic regions and boundary crisis. When the parameter *γ* changes in different value ranges, what more complex dynamic properties of the model will be a problem worthy of further study.

In biology or ecology, the complex chaotic behavior of this mapping shows the relationship among the number, birth rate, and survival rate in a population, whether it survives in a balanced state or makes the population develop in disorder or chaos. This research can provide theoretical basis and help for the research in biology or ecology. For example, it is applied in the marine fishing, or in the reproduction and population growth of a certain species in nature.

#### Data Availability

No data, models, or code were generated or used during the study.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the High Level Scientific Research Project Cultivation Fund of Shandong Women’s University and Zhejiang Provincial Natural Science Foundation of China (Grant no. Y515387).