On stable and quasi-chaotic regimes in a one-dimensional unimodal mapping obtained by modeling the dynamics of a biological population

The paper considers the properties of the difference equation describing the dynamics of the animal population, obtained earlier in the framework of studies of tundra communities. We have considered a special case in which the model is represented by a one-parameter difference equation that defines a one-dimensional unimodal mapping of a segment into itself, similar to the well-known triangular (tent) mapping, supplemented by a region with a constant value. A change in the mapping parameter generates a bifurcation scenario, in which stability zones arise, characterized by orbits of a constant period, interspersed with zones with more complicated, “quasi-chaotic” regimes. Based on the properties of the n-iterated triangular mapping, a necessary and sufficient condition for the localization of cyclic orbits in the considered type of unimodal mappings is formulated, which makes it possible to identify stability regions for any given period n. On the basis of this condition an algorithm for detecting stability zones is proposed. The main subject of the study is the fractal properties of the set, which is the complement of the obtained set of stability regions to the entire domain of mapping definition. The dynamics of the DH (n) value is obtained, the limit of which at n → ∞ is equal to the fractal dimension dH . It is shown that in the studied range of n (2 ≤ n ≤ 22), the DH (n) < 0.9, which suggest that dH < 1. If so, then according to the definition of a fractal set, its topological dimension is dT = 0, which means that the complement of the set of stability regions consists of isolated points.


The origin of the problem. Triangular mapping with a plateau as a model of population dynamics
Mathematical modeling of processes of various nature in discrete time is closely related to the use of first-order difference equations, or succession functions, when the value of the function at the current moment of time depends on its value at the previous moment: x (n+1) = f (x (n) ).One of the main questions that arise in such models is related to the existence of cyclic trajectories in them and their stability.Since the subsequent value of the variable is uniquely determined by the previous one, the repetition of the same value twice means that the entire dynamics is completely repeated and localized in a limited region of the phase space.The opposite is the question of the existence of non-cyclic or chaotic trajectories.If such trajectories exist, then any exogenous shock (a short-term sharp change in external conditions) or a regular change in exogenous parameters can lead to the system being in a chaotic regime.Such models and related mathematical problems are actively studied within the framework of population modeling, but their scope is much wider (including both problems of the natural sciences, such as the cycles of division of microorganisms in a nutrient medium, and problems of the social sciences, in particular, economic cycles, finance, dynamics of public opinion, etc.).The most well known are the works of Robert May [1], where the term "chaos" was used in describing the behavior of a biological population model, and Tien-Yen Lee and James York [2], who connected the existence of a cycle of period three with chaotic behavior.A more general condition for the coexistence of cycles was given earlier by A.N. Sharkovskii for an arbitrary continuous mapping of a segment into itself [3] regardless of stability, while in [1,2] both the presence and stability of cyclic regimes were considered using the example of a logistic mapping, which later became often used to study the properties of difference equations.
An example of a model in discrete time is the one proposed by Sarancha et al. a simplified model of the tundra animal population rate in the form of a difference equation [4,5], which differs from the generally accepted logistic one.The equation specifies a piecewise linear three-parameter mapping of a segment into itself (Fig. 3, equation on page 7 in [4]), close to the well-known triangular (tent) mapping, supplemented by a region with a constant value -a plateau.
From a biological and ecological point of view, autocorrelation modeling of population dynamics (the dependence of the population size at each moment of time on its own size at previous moments) corresponds to the assumption of the invariance of the external conditions in which the population exists (i.e., volume, rate of renewal of resources and etc.).As long as the population size does not exceed a certain threshold level, the population grows at a rate corresponding to the rate of reproduction.Once the threshold is exceeded, the population causes serious damage to the ecosystem (depletes resources), and then its abundance begin to decline until the ecosystem recovers.
In turn, the presence of a plateau means the existence of a certain "optimal biotope" for a given population, i.e. such a combination of external conditions under which a certain level of abundance remains constant [6].Thus, the population dynamics is described by a piecewise linear function with three modes (growth, decline and equilibrium).At the same time, in the computational aspect, the presence of a plateau in the mapping function means the existence of a stable trajectory -an attractor with an attraction region that coincides with the entire domain of the function.
Later, in [7], as a result of calibrating the obtained difference equation using the identification sets method [8], it was possible to reproduce the temporal dynamics close to the dynamics of the number of real populations.
Note that in this paper we are only interested in the mathematical properties of this mapping, while the given population interpretation is only an example of its practical application.When studying the properties of the resulting model in [9,10], using the example of some approximation, which is a triangular mapping with a plateau (see formula (5.2), on p. 73 in [9]), the properties were described that allowed us to formulate a necessary and sufficient condition for the localization of cyclic trajectories in the considered type of unimodal mappings of a segment into itself, based on the following statements.

Properties of cyclic trajectories in a triangular mapping and in a triangular mapping with a plateau
In the triangular mapping of the segment [0, 1] into itself: the set of cyclic trajectories of period n consists of all those and only those fixed points of the niterated mapping f (n) = f ( f (n-1) ) (i.e.intersection points of the graph of the mapping f(n) with the diagonal at [0,1][0,1]) (2), which are not simultaneously fixed points of the mapping f (m) , where m is a divisor of n, including two fixed points of the mapping f (1) (in this case m = 1).
Here   () is the x-coordinate of the i-th intersection of the n-iterated mapping f (n) with the diagonal; i is the number of the linear segment of the piecewise linear function f (n) , whose intersection with the diagonal gives the point   () .
Thus, we know the set of points that make up the entire set of cyclic trajectories of period n, as well as the number of trajectories (p. 25 in [11]).However, the "partition" of this set into separate trajectories is unknown.Consider the question of localization of trajectories of period n on the set [0,1].By the localization of a cyclic trajectory, we mean the range of values enclosed between the points with the smallest and largest coordinates of this trajectory.. Definition: The "boundary" point of the cycle will be called the point of the cyclic trajectory with the largest x-coordinate.
Let's list the properties of the fixed points of an n-iterated triangular mapping.
Prop. 1.The points of intersection of the sides of some tooth of the graph f (n) with the diagonal always belong to different cyclic trajectories of period n in the mapping f (1) .
Prop. 2. The points of cyclic trajectories of period n of the set described above are boundary if and only if the segments on the diagonal bounded by them, enclosed inside the teeth f (n) , are not embedded in any segment formed on the diagonal by the tooth of the mapping f (m) , m < n.
From properties 1 and 2 follows: Prop.3. If the point of intersection of the side of some tooth f (n) with the diagonal is the boundary point of the cyclic trajectory of period n in the mapping f (1) , then the point of intersection of the other side of this tooth with the diagonal is the boundary point of another cyclic trajectory of period n in f (1) .Except for the case when the left boundary of the segment (formed by the intersection of the sides of the tooth f (n) with the diagonal) coincides with the boundary point of the segment corresponding to the cycle, the period of which is a divisor of n.
Let E be a non-trivial fixed point of the triangular mapping, f(x E ) = x E , x E  0. Then any cyclic trajectory in the triangular mapping has at least one point x > x E . Proof: For x  (0, x E ): f(x) > x, therefore, for any x i from (0, x E ) we have x i+1 > x i , and x i+1 belongs to either (0, x E ), or (x E , 1).On the contrary, x  (x E , 1): f(x) < x, moreover, x  (x E , 1): f(x) < x E , therefore, for any x i from (x E , 1) x i+1 < x i ,, and x i+1 always belongs to (0, x E ).▲ Corollary of Statement 1.The boundary point of any cyclic trajectory in the triangular mapping belongs to the interval (x E , 1).
Thus, to determine the localization of cyclic trajectories of some period n in a triangular mapping, it suffices to find all boundary points of cycles of this period in the interval (x E , 1).) ,  +1 () ](x E , 1)} -the set of segments bounded by the teeth of the graph f (n) on the diagonal in the interval (x E , 1), projected onto the X axis, where the points x i , x i+1 are determined by formula (2).Statement 2. All interior points of an arbitrary segment [x i , x i+1 ]  X (n) fall into the half-interval in n steps of the mapping f (1) .

Stable cycles in a triangular mapping with a plateau
Proof: The function f (n) is continuous on the entire interval (x i , x i+1 ), strictly monotonically increases on the interval (x i , ), and on the interval ( ) = 1.In other words, the mapping f (n) : (x i , bijective and increasing, while the mapping f (n) : ( ] is bijective and decreasing.▲Therefore, if the boundary points of the segment [x i , x i+1 ] are the boundary points of the cycles of period n in the mapping f (1) , and a horizontal segment (the "plateau") of the mapping passes through some point x A (x i , x i+1 ) (Fig. 1, formula (3)), then the trajectory that left the point x(x i , x A ) will reach a plateau in k•n cycles: x (kn)  [x A , 1], followed by a transition to the n-cycle, including the points Due to the fact that the derivative f (n) ′(x A ) = 0, such a cycle is stable.
where point А(x A , 2(1-x A )) is the left boundary of the plateau.So, in contrast to the triangular mapping (1), in which all points are unstable, including the fixed 0, E, in the mapping with a plateau (3) stability regions appear, the boundary points of which are the boundary points of n-cycles in the triangular mapping.
Therefore, if in a triangular mapping each period n corresponds to a finite number of unstable cycles [11], then in a triangular mapping with a plateau there is an uncountable set of stable n-cycles in each segment [ x i , x i+1 ].If x A is the transition point of the descending side of the triangular mapping to the plateau (Fig. 1 ), and x A belongs to such a segment: x A [ x i , x i+1 ], then x A corresponds to some stable n-cycle with attraction domain [0, 1].By varying the value of x A from the fixed point x E to 1, one can obtain a successive change of stable cycles satisfying the Sharkovskii order [3], separated by regions of "chaotic dynamics", whose trajectory period cannot be estimated.In this case, the mapping (3) preserves those unstable cycles, all of whose points lie to the left of x A .Thus, adding a plateau region to a triangular mapping reduces the problem of localizing unstable cycles to determining the boundaries of segments from the interval (x E , 1) in which the corresponding stable cycles are realized.This allows solving the problem not only analytically, but also computationally, by implementing the bifurcation scenario, changing the position of the plateau region in the interval (x E , 1), however, this approach is associated with a large amount of calculations.
Based on the foregoing and formula (2), an algorithm is constructed that sifts through the set of segments formed by the teeth f (n) and excludes nested segments.Removing the sets X (n) , n = 2, ... , n max from the initial segment [x E , 1] one by one, we obtain a process similar to the construction of the wellknown Cantor set (for more details on the construction of the set, see Section 2).Note that by deleting the segments corresponding to the cycles of period n, we also delete the segments of cycles of large periods nested in them.Thus, in fact, the process is reduced to the removal of segments bounded by the largest points of the cycles of periods n = 2, ... , n max , in which, if they have a plateau, n-cycles are realized.
However, if we associate these removed segments, in which there are cycles, with the removed intervals of the classical Cantor set, then we get a set whose topological and fractal properties are different from those of the Cantor set.At least the following differences can be distinguished: i. by construction, this set consists of open subsets (that is, the intervals that make up this set do not contain their boundary points); ii. the exhaustibility of the original segment by the removed subsets has not been proved; iii.although a deterministic algorithm is used in the construction (i.e., there is no random parameter), nevertheless, for the resulting object, an analytical description of regularity is not known, as in the case of classical regular fractals.The purpose of this work is to establish a connection between the value of the fractal and topological dimensions of the set under consideration and such a property as the "exhaustibility" of the original set, which is a segment [x E , 1] by segments containing cycles (see the description of the construction of a pseudo-Cantor set).According to the definition of a fractal set given by B. Mandelbrot, the fractal dimension of a set or the Hausdorff-Besikovich dimension d H is always strictly greater than its topological dimension d T [12].Note that the fractal dimension is defined in the limit state of the set as the limit of a continuous function of variable n, while the topological dimension, which is a discrete value, is defined for all intermediate states of the set and is equal to the dimension of the original set for each finite value of n.The topological dimension of the set in the limit state may differ from that of the original set.So, if 0 < d H < 1, then the topological dimension is d T = 0, which means that in the limit state any point of the set has an arbitrarily small neighborhood with an empty boundary ( [13], p. 162).If in the limiting state with the largest period of the cycle n in the removed segments tending to infinity, the set is topologically zero-dimensional, then we will say that the set is "exhausted".

Numerical estimation of the metric dimension based on the total length of intervals and rounding error
One of the basic definitions of the fractal dimension based on covering by balls (in the Euclid or Chebyshev metric) is known as the Minkowski dimension [14].By definition, it is equal to where M(ε) is the minimum number of elements required for the covering, ε is the diameter of the coating element.It should be noted that since when determining the Hausdorff dimension d H , covering elements of different diameters are used, while when determining the dimension d, all elements of the same diameter are used, then the relation d  d H is fulfilled.As a result, for a numerical estimate of the fractal dimension from above, it is permissible to use expression (4).
Since the fractal dimension is estimated for the limit state of the pseudo-Cantor set, for which there is only an algorithmic, but not an analytical description, we need to consider in more detail the process of its construction.The segment [x E , 1] is taken as the initial set, from which, at each n-th step, the segments corresponding to the cycles of period n described in Section 1.3 are removed: X (n) = {[  () ,  +1 () ](x E , 1)}.We obtain an infinite sequence of states of the set of intervals on the segment [x E , 1] remaining at step n after the removal of segments X (n) : [x E , 1]\{ X (k) , k = 2, ..., n}, with n.Visualization of the set structure is presented in subsection 2.2.
Usually, when calculating the metric dimension by formula (4) in the case of classical regular fractals, the diameter of the covering element ε is expressed as a function of the step number n.However, for the set we are considering, unlike most regular fractals, it is difficult to establish such a value of ε, which depends on n, by which the lengths of all existing segments at a given step are divided.To overcome this difficulty, we suggest using the total length of segments and estimating the rounding error when dividing by the value ε(n) corresponding to the given step.Then, as a value characterizing the covered set, we take the total length of intervals at the n-th step that do not contain cycles of periods from 2 to n: () = ∑ δ(, ) () =1 , where δ(n, j) is the length of the j-th interval that does not contain cycles of periods from 2 to n, j max (n) is the number of such intervals at step n.
The number of elements of diameter ε required to cover the set at step n is expressed by the approximate ratio: 4) can be replaced by the limit: Denote by D(n) the value under the limit sign for finite n: Let us determine the rounding error when calculating the number of segments of size ε(n) required to cover the set at the n-th step of constructing the set.In this case, the condition ε(n)  δ min (n) must be satisfied, where δ min (n) is the smallest length of the segment contained in the set at step n.
The number of segments needed to cover the set of intervals in step n: () + θ  (, ε()), where the rounding error is: Here δ(n, j) is the interval number j at the n-th step, j max (n) is the number of intervals at the n-th step.Because each interval δ(n, j) has a size of at least ε(n)  δ min (n), then by the definition of rounding up the quotient: Summing up both parts of inequality ( 8) over the number of covered intervals from 1 to j max (n), we obtain: Substituting equality (7) into inequality (9), we obtain the upper estimate of the rounding error: and, accordingly, the upper bound on the number of segments of the set cover: () ⌋ +   ().The lower bound on the number of coverage segments is reached when  m = 0 (if the length of each interval (n,j) of the covered set is a multiple of the length of the smallest interval  min (n), see equality (7)): Based on the obtained values of the lower and upper bounds of the number of coverage segments  and , we determine the lower and upper bounds of the value D(n) (see formula (6)), the limit of which for n   is equal to the metric dimension of the set: There are various ways to select a coverage element.

Coverage by δ min segments.
As an element of coverage ε(n), we choose a segment equal to the minimum size of the interval that does not contain cycles at the n-th step: δ min (n).Note that δ min (n) decreases monotonically (not strictly) as n increases, and the limit lim δ min (n), as n   is undefined.The values of () and () depending on the step n of constructing the set when covered with segments δ min (n) are presented in Table 1(Supplement 1) and on the graph in Figure 2.  (10,11)), when covered with segments δ min (n), depending on the maximum cycle period n.

Coverage by g min segments.
Let us choose as a coverage element ε(n) a segment with a value equal to the size of the smallest segment containing a cycle of period n, removed at the n-th step:  (n).For a triangular mapping, the value g min (n) is determined from formula (2) for i = 2 n and is equal to: Unlike δ min (n) -the value of the smallest interval that does not contain cycles of period less than or equal to n, the value of g min (n) decreases strictly monotonically and has the limit lim g min (n) = 0, n   .At the same time, expressions (10) and (11) for the lower and upper boundaries of the value D(n) remain valid for g min (n).
Table 2 (Supplement 2) and the graph in Figure 3 present the values of () and () depending on the step n of constructing the set when covered by segments g min (n).(10), (11)), when covered with segments g min (n), depending on the maximum cycle period n.The results of the numerical experiment show that, regardless of the type of coverage element, the rounding error Θ, equal to the difference between the upper and lower estimated boundaries of the value D(n), decreases monotonically with an increase of the maximal cycle of period n, which is a necessary condition for the existence of a limit of D(n) at n  .

Estimating the metric dimension of a "pseudo-Cantor" Set by the Deep Holes Method
Another method for estimating the metric dimension is simulation modeling of covering that is close to optimal, based on the Deep Holes Method.If optimal (or close-to-optimal) metric ε-net [15] is constructed with radius (or diameter) ε  0, approximating a metrically bounded set, then it is possible to determine the metric dimension, since it forms an ε-covering with a certain number of elements (balls in the chosen metric).(An optimal net is a metric net containing the minimum possible number of balls.)However, there are no universal algorithms that provide approximation of any set by optimal coverage, especially for small ε (because the complexity of the distance matrix NxN elements construction increases quadratically with the number of elements in the net).Therefore, in practice research, (ε,δ)-nets [16] are used, which, with metric accuracy (coverage radius) ε/2 and completeness in measure δ (that is, the coverage covers not 100% of the original set, but its given share 1 -δ with a given reliability).
To construct semi-optimal metric networks, the Deep Holes Method (DHM) [17] can be used: an adaptive iterative method for approximating completely bounded sets in chosen metric space based on near-optimal metric ε-networks and ε-distinguishable subsets construction (coverings and packings of metrics balls).The general principle of the method is as follows: the set of balls that form the covering supplementation with new elements that are the most distant from those already present in the covering.This makes it possible to obtain an estimate of the number of balls in the cover M and the radius of the cover ε, the accuracy of which depends on the duration of the simulation experiment and the size of the covered set.
For our experiments, we used a variant of the DHM based on the selection of points from a net covering the original set.A feature of this approach is additional restrictions on the original net (it must reflect the topology of the set, i.e., show emptiness between segments).
We should also note the internal methodological features of the DHM: it gives an estimate of the number of points in the coverage "from above", since the metric net is close to optimal, but not optimal (the number of points in the covering may be slightly too high).This error is the greater, the more primitive the topology of the object (in our case, the more the object looks like a continuous segment without emptiness).Accordingly, for the first steps of the set construction, the dimension estimate turns out to be significantly overestimated.This, however, is not a disadvantage, since the definition of the metric dimension is, a construction of step tending to infinity.Accordingly, our ability to reproduce the exact geometry of the intervals of the set at the latest steps becomes a real limiter in the simulation study of the metric dimension.The number of segments covered by elements and the distance between them characterize the fractionality of the set (for 21-22 steps, these calculations can take several days).
A graphical representation of stepwise process of construction a set in the form of a set of layers corresponding to step n, based on the results of simulation experiments, is shown in Figure 4. Here, gaps denote excluded areas with k-cycles, where k = 2, ... , n.The calculated values of the metric dimension for different construction steps are shown in Figure 5.The analysis shows some nonmonotonicity of the metric dimension estimation, which is explained by the alternation of even and odd set construction steps.The general trend, however, confirms the metric dimension approximately ~0.9 for sufficiently small ε and big n, which corresponds with the estimation based on the total length of the segments.

Conclusion
The problem of calculating the metric dimension of the Pseudo-Cantor set formulated in this work arose in the process of studying the properties of the difference equation (3), which describes the dynamics of the population size of tundra animals, which were previously considered in [4,5].The equation, in contrast to the well-known logistic one, is a triangular mapping, supplemented by an area with a constant value -a plateau.In these works, along with the type of equation, the choice of the bifurcation scenario was also substantiated: a change in the height of the plateau, which characterizes the capacity of the optimal biotope in the case of a specific model [6].At the same time, the mapping remains unchanged in the rest of the domain of definition.
The computational experiment [18] showed the presence of regions in which cyclic trajectories have a constant period, separated by zones with more complicated, "quasi-chaotic" regimes, which are characterized by a drastic change in the period length with a small change in the bifurcation parameter.formulated, which makes it possible to identify stability regions for any given period n, on the basis of which an algorithm for detecting such zones is proposed.
By stepwise removal of regions with stable n-cycles from the initial segment [x E , 1] of the X axis containing all the boundary points of the cycles, as the period n increased, a new set similar to Cantor's set was constructed, which, as we assumed, has fractal properties.
For the complement of the set of cycle stability regions up to the segment [x E , 1] obtained in this way, the metric dimension was estimated in different ways.ii.Based on the methods of simulation modeling of covering close to optimal, using the Deep Holes Method [17].Despite the differences in approaches (due to its specificity, the DHM can somewhat overestimate the required number of points in the coverage), different methods show a close estimate of the metric dimension d H ≈ 0.9.
At the same time, according to the definition of a fractal set given by B. Mandelbrot [12], if the metric dimension of the set is d H < 1, then its topological dimension is d T = 0.This means that, in the limit, the complement of the set of regions of stability up to the entire domain of mapping definition (in this case, this is the segment [x E , 1] ) consists of isolated points.Consequently, in the limit, the entire region [x E , 1] is exhausted by the set of stability regions.

Figure 1 .
Figure 1.A triangular mapping of the segment [0,1] into itself, supplemented by a region with a constant value -a plateau.The transition point of the descending side of the triangular mapping to the plateau x A belongs to the segment [ x i , x i+1 ] -the projection onto the X axis of the diagonal segment bounded by the tooth of the n-iterative mapping function f (n) (x).

Figure 2 .
Figure 2. The lower and upper bounds of the value D (see formulas(10,11)), when covered with segments δ min (n), depending on the maximum cycle period n.

Figure 3 .
Figure 3.The lower and upper bounds of the value D (see formulas(10),(11)), when covered with segments g min (n), depending on the maximum cycle period n.
Upper and lower bounds of value D(n)vs maximal cycle period n, g min -coverage

Figure 5 .
Figure 5. Estimation of metric dimension D(n) based on imitational experiment using the DHM, for set construction steps n from 2 to 22.
i. Numerically, based on the values of the coordinates of regions with stable cycles, the dynamics of the quantity   () = − log (ε()) log ε() is obtained, the limit of which at n   is equal to the fractal dimension: d H = lim →∞   ().Here ε(n) is the diameter of the covering ball, M(ε(n)) is the number of balls in the covering, where n is the longest cycle period.In the studied range of n values (2  n  22), the value D H (n) < 0.9, which suggests that d H < 1.