Reinterpreting the Origin of Bifurcation and Chaos by Urbanization Dynamics

Chaos associated with bifurcation makes a new science, but the origin and essence of chaos are not yet clear. Based on the well-known logistic map, chaos used to be regarded as intrinsic randomicity of determinate dynamics systems. However, urbanization dynamics indicates new explanation about it. Using mathematical derivation, numerical computation, and empirical analysis, we can explore chaotic dynamics of urbanization. The key is the formula of urbanization level. The urbanization curve can be described with the logistic function, which can be transformed into 1-dimensional map and thus produce bifurcation and chaos. On the other hand, the logistic model of urbanization curve can be derived from the rural-urban population interaction model, and the rural-urban interaction model can be discretized to a 2-dimensional map. An interesting finding is that the 2-dimensional rural-urban coupling map can create the same bifurcation and chaos patterns as those from the 1-dimensional logistic map. This suggests that the urban bifurcation and chaos come from spatial interaction between rural and urban population rather than pure intrinsic randomicity of determinate models. This discovery provides a new way of looking at origin and essence of bifurcation and chaos. By analogy with urbanization models, the classical predator-prey interaction model can be developed to interpret the complex dynamics of the logistic map in physical and social sciences.


Introduction
Chaos is one of the important subjects of science in the twentieth century. However, the problems of origin and essence of chaos were not really solved in last century, and they are passed on to the new century. The simplest model for understanding chaos is the well-known logistic map. The complicated behavior of the logistic growth brought to light by May [1] led to a profound insight into complex dynamics. Thus, chaos is always regarded as intrinsic randomicity of determinate dynamical systems. A pending question is how and why determinate systems have complicated behavior. Many studies are devoted to this problem, and many interesting conjectures are proposed. But the essence of bifurcation and chaos is still puzzling. In fact, a revealing research can be made from the viewpoint of urban dynamics. Urbanization provides a new way of understanding the origin and essence of chaos. Urban systems are complex systems, and the process of urbanization and urban evolution are nonlinear process associated with chaos and fractals [2][3][4][5][6][7]. Using mathematical derivation, numerical computation, and empirical analysis, we can reveal new knowledge about bifurcation and chaos based on the nonlinear dynamics of urban evolution.
New progress may be made by a simple formula of urbanization ratio. A basic and important measurement of urbanization is the proportion of urban population to the total population, which is termed "level of urbanization" in urban geography. The curve of urbanization level is termed "urbanization curve" and can be described with sigmoid functions such as logistic function, which can be discretized to a one-dimensional map. Using the formula of urbanization level, we can derive the logistic equation from the rural-urban population interaction model, which can be discretized to a two-dimensional map. Thus the one-dimensional logistic map can be associated with the two-dimensional rural-urban interaction map. As will be shown below, the two-dimensional rural-urban map can create the bifurcation and chaos that are identical in patterns to those produced by the one-dimensional logistic map. This suggests that the origin of bifurcation and chaos is two-population coupling and interaction rather than intrinsic randomicity of determinate models [8].
The study of chaos associated with bifurcation can help us understand natural and social systems deeply. This paper is a development based on a series of previous studies [8][9][10][11][12][13][14]. The rest of this work is organized as follows. In Section 2, the bifurcation and chaos from ruralurban population interaction dynamics are illustrated by using a two-dimensional map, and a phase portrait analysis of rural-urban interaction is performed. In Section 3, an empirical analysis is made by means of American census data to verify the rural-urban interaction model. The case study lays the foundation of experiments for the urbanization model. In Section 4, several related questions are discussed. First, the two-population interaction model is generalized to explain the ecological phenomena including logistic growth and oscillations of population. Second, the scaling laws of period-doubling cascade are compared with those of hierarchy of cities. Third, the nonlinear dynamics of urbanization curve is further generalized to the fractal dimension curve of urban growth. Fourth, the nonlinear replacement dynamics is outlined. Finally, the discussion is concluded with a brief summary.
world [4]. Based on several assumptions, the spatial interaction model for rural-urban migration can be expressed as below [8]: in which r(t) and u(t) denote rural and urban populations at given time t, respectively, and a, b, and c represent three parameters of population transition. Please note that r(t) > 0, u(t) > 0. This model indicates that the rural-urban population interaction results in urbanization. According to Eq. (1), the growth rate of rural population depends on rural population size and the twopopulation interaction, while that of urban population growth rate only depends on the ruralurban interaction. If the study region is a close system, then the parameters b and c are equal to one another, i.e., b = c, or else they are not. Eq. (1) has a firm basis of statistical analysis. The model can be verified with the population data set of American census since 1790.
It can be proved that the system of differential equations on rural-urban interaction is equivalent to the logistic equation of urbanization curve. For simplicity, Eq. (1) can be rewritten as follows: in which.
In urban geography, the level of urbanization is formulated as where L(t) denotes urbanization level at time t (obviously 0 ≤ L(t) ≤ 1). The level of urbanization is an important measurement in urban study. Just because of the definition of urbanization level, the one-dimension map of logistic growth can be associated with the two-dimension map of rural-urban interaction. In fact, taking the derivative of Eq. (4) yields Substituting Eq. (2) into Eq. (5) gives For simplicity, we can postulate that the region is a close system, which has no population exchanged with outside. In this case, we have b = c and b * = c * , and thus we have As indicated above, c * = c/[r(t) + u(t)], Eq. (7) can be reduced to Based on the level of urbanization defined by Eq. (4), the logistic equation is readily derived as below: In literature, Eq. (9) is always expressed as follows: in which k is just the original rate of growth in the logistic model (k = b-a = c-a). Discretizing Eq. (10) yields a one-dimension logistic mapping. By means of the one-dimension mapping, May [1] created the period-doubling bifurcation and chaotic patterns, which are familiar to many scientists of chaos and complexity.

Bifurcation and chaos based on two-dimensional map
Discretizing the rural-urban population interaction model yields a two-dimensional map, which can be employed to make numerical analysis. Since Eq. (10) can be derived from Eq. (1) through mathematical transformations, we expect that the complicated dynamical behaviors such as period-doubling oscillation and chaos can also be created by the two-dimension maps based on Eq. (1). Discretizing Eq. (1) yields a pair of iterative functions such as in which the parameters α, β, and γ for discrete form correspond to a, b, and c for continuous form in Eq. (1), respectively. The parameters in Eq. (11) will vary slightly after continuousdiscrete transformation. If the region is a close system, we will have β = γ. In light of the US urbanization model as well as the American census in 1790, the parameter α can be set as α = 0.025, and the initial values can be set as r(0) = 3.727559 and u(0) = 0.201655. Thus the behavior of Eq. (11) will depend on the values of the parameters β and γ. The numerical computation can be fulfilled through a common computer with mathematical software. It is expected that the numerical behavior of the two-dimension maps on the base of Eq. (11) is really identical in form to the complicated performance of the one-dimension map based on the logistic function in ecology ( Figure 1).
The numerical iterations can be fulfilled by mathematical software such as MATLAB or even by the well-known spreadsheet, Microsoft Excel. In order to correspond the two-dimension rural-urban maps to the one-dimension logistic map, a limiting condition is set as β = γ.
The iterative values represent the rural and urban population in different times. Using Eq. (4),  we can convert the rural and urban population into the level of urbanization and obtain the urbanization curve. The main results can be summarized as follows. (1) Logistic decay and growth. When β = γ < 0.025, the urbanization curve takes on a monotonous decreasing graph, and the final value is L min = 0; when 0.025 < β = γ < 1.032, it displays a monotonous increasing graph, and the final value is L max = 1. The latter represents the common logistic curve that is familiar to geographers. (2) Steady-state behavior. When 1.033 < β = γ < 2.025, the urbanization curve exhibits an alternating change and finally changes to unit ( Figure 1a). (3) Perioddoubling bifurcation. When 2.025 < β = γ < 2.475, the urbanization curve shows an oscillation of period 2 ( Figure 1b); when 2.475 < β = γ < 2.571, the curve displays an oscillation of period 4 ( Figure 1c); and with β and γ increasing, an oscillation of period 8 (b = c > 2.571) and period 16 (b = c > 2.591) and period 2 n (the natural number n > 4) emerges step by step. Finally, when β = γ > 2.61, the urbanization curve evolves into chaotic state, in which no 2 n cycle can be detected. The upper limit of the parameters β and γ is about 3.033. That is, if β = γ ≥ 3.033, the numerical iteration will break down [8].
A comparison can be drawn between the results from the one-dimension logistic mapping and those from the two-dimension rural-urban interaction mappings. An interesting finding comes from the comparative analysis. In fact, Eq. (10) can be discretized as a one-dimension mapping such as L(t + 1) = (1 + K)L(t)-KL(t) 2 , where the parameter K corresponds to the parameter k in Eq. (10). Then we have K~β-α. Using this one-dimension logistic mapping, we can obtain various urbanization curves. The common characters of the bifurcation and chaos from the onedimension mapping and the two-dimension mappings are the same with one another. What is more, the critical values of the model parameters for the period-doubling bifurcation and chaos are approximate to those of the logistic model. In particular, according to the process of numerical experiments, if the parameter α value becomes small enough, the critical value of the periodic oscillation to chaos based on the rural-urban interaction is close to the value based on the logistic growth. This discovery suggests a new way of looking at the origin and essence of bifurcation and chaos. It is the interaction rather than the intrinsic randomicity that causes the complicated behaviors of a simple dynamic system.
Another finding is the inherent relation between order and chaos. There are narrow ranges of periodic solutions in the chaotic "band." If β = γ > 2.857, the urbanization curve takes on an oscillation of period 3. Further, when β = γ > 2.871, an oscillation of period 5 or period 6 or period 7 appears. All in all, a non-2 n cycle can be found in the chaotic state. Finally, when β = γ > 2.88, the urbanization curve will get into a random state once again. However, the periodic oscillations in the chaotic belt are different from the cycles in the process of perioddoubling bifurcation. The period in bifurcation is of 2 n cycle, while the oscillation in chaos is of non-2 n cycle. The varied non-2 n cycles such as the period 5 and the period 6 indicate chaos rather than bifurcation. The typical non-2 n cycle is period 3 [15]. In short, the limited chaos can be regarded as the sum of randomicity and the cycles of non-doubling period [8]. One the other hand, some slight disorder can be found during the 2 n cycles, which can be revealed by spectral analysis. This suggests that chaos and order cannot be absolutely separated, and they contain one another or are interwoven with each other.
One of properties of chaos is the sensitive dependence on the initial conditions. The property can be testified by the urbanization curve based on the rural-urban interaction mapping. Suppose that the parameter values are given as α = 0.025 and β = γ = 2.785. Let us change the initial rural population from r(0) = 3.727559 to r(0) = 3.727558 (million) but keep the initial urban population u(0) = 0.201655 unchanged. In this case, the numerical iterative curve shows a new urbanization trace. At the beginning, the new urbanization curve and the old urban curve almost coincide with one another; but gradually, the new curve deviates from the old one ( Figure 2). The difference between the two urbanization curves becomes bigger and bigger over time. It should be noted that only one person is reduced from the initial rural population. A minimal error results in wide divergence. This just reflects the sensitive dependence on the initial conditions of the rural-urban interaction.
Urban chaos is an interesting issue, but it seems to appear in the mathematical world instead of the physical world. The model parameter values such as α = 0.025 are determined by the US census data. In terms of Eq. (4), the level of urbanization ranges from 0 to 1, i.e., 0 ≤ L(t) ≤ 1. It will make no sense if L(t) < 0 or L(t) > 1. On the other hand, according to the observational data from the real world, the parameter β and γ values should come between 0.025 and 1.032. Otherwise, the urbanization level L(t) will be less than 0 or greater than 1. Unfortunately, if the parameter β and γ values are confined into 0.025 and 1.032, no bifurcation and chaos will emerge in the rural-urban mapping process. If and only if β = γ > 1.032, we can create the period-doubling bifurcation and chaos, and thus the level of urbanization will exceed 1. This is absurd in both mathematics and urban geography. It is one of the necessary conditions of urban complicated behaviors. This suggests that the urban bifurcation and chaos emerge in the possible world, and we cannot find them in the real world for the time being.

Phase portraits of two-dimension map
Using the two-dimensional map, we can draw the phase portraits of the logistic process based on the one-dimensional map. The spatiotemporal feature of urbanization dynamics can be revealed with the phase portraits. Taking rural population r(t) as x-axis and urban population u(t) as y-axis, we can create a set of scatterplots for the period-doubling bifurcation and chaotic behavior of urbanization. The plots show the rural-urban relationships defined in the phase space based on Eq. (11). Consequently, the period-doubling process can be characterized by 2 n radials (n = 1, 2, 3, …), and the cross point of the rural and urban radials is just the origin of coordinate (Figure 3a-c). If the level of urbanization evolves from bifurcation into chaos, all the scattered points are randomly confined in the triangular region defined by the intersectant rural and urban radials ( Figure 3d). For the chaotic state, the radials indicative of non-2 n cycle Figure 3. The phase portrait of the period-doubling bifurcation and chaos of rural-urban population interaction (Note: the times of iterations are 2500. The four subplots in Figure 3 correspond to the four subplots in Figure 1, respectively. See [8]). may appear in the phase portraits. The feature of phase portrait is independent of the times of iteration. The spatial distribution of scattered points never converge, and this suggests that there is no strange attractor in the phase space of urban chaos. This inference differs from the traditional understanding on the chaotic dynamics based on logistic growth.
Despite the fact that no chaotic attractor can be found, these scatter points follow certain mathematical rule. The distance from a data point (r(t), u(t)) to the origin (0, 0), i.e., the cross point of radicals which act as boundaries of these points, can be formulated as which quantifies the spatial relationships of the scattered points. Thus the distribution of the scattered points in the phase space meets a logarithmic relation as below: in which N(d) refers to the cumulative number of the scattered points within the distance d, and A and B are two parameters representing the slope and intercept, respectively. Now, let us examine mathematical structure of the phase space from the perspective of statistics. The distance is taken as d = 4 n , where n is a natural number, and the number of iterations is set as 5000. As an example, if the parameter value β = γ = 2.785 as given, the estimated values of the parameters in Eq. (13)   The derivative of the logarithmic function is a hyperbolic function. This implies that the density of the points in the phase portrait of chaotic state decays gradually from the origin, and the density change can be characterized by a hyperbolic curve. Despite the fact that both a city and a system of cities bear fractal structure [3-5, 16, 17], the phase portrait of the urban chaos does not display self-similar pattern. The reciprocal function of the logarithmic function is just an exponential function. This suggests that the basic property of the logarithmic distribution can be understood through the exponential distribution. Compared with the Gaussian distribution, the exponential distribution implies complexity [18], while compared with the exponential distribution, the power-law distribution implies complexity [19,20]. This suggests that complexity seems to be a relative concept. Exponential distribution falls between the simplicity based on normal distribution and the complexity based on power-law distribution.
According to the dual relation between the exponential function and the logarithmic function, the logarithmic distribution of the scattered points in the phase space of urban chaos indicates a process appearing between simplicity and complexity.

Data and method
The above-shown numerical iterations are based on the two-dimensional map from the ruralurban population interaction model. It is necessary to make empirical analysis using the dynamic equations of urbanization and observational data. There are two central variables in the study of spatial dynamics of city development: population and wealth [21]. According to the aim of this study, only the first variable, population, is chosen to test the models on urban chaos. In fact, population represents the first dynamics of urban evolution [22]. Generally speaking, the population measure falls roughly into four categories: rural population r(t), urban population u(t), total population P(t), and the level of urbanization indicative of the ratio of urban population to the total population, L(t). The measure relations are as follows-P (t) = r(t) + u(t) and L(t) = u(t)/P(t)-which can be found in Eq. (4).
The American rural and urban data comes from the US ten-yearly population censuses. There are 23 times of census data from 1790 to 2010 available on the website of American population census. However, only the data from 1790 to 1960 are adopted in this work ( Table 1). In fact, the definition of cities in America was changed in 1950, and the new definition came into use since 1970. The US urban population caliber after 1970 may be inconsistent with that before 1960. The observational data can be fitted to the discretization expressions of the United Nations model [23] and the generalized Lotka-Volterra model [24][25][26], respectively. The parameters of models are estimated by the multiple linear regression based on the ordinary least squares (OLS) method. The advantage of the OLS method is to keep the key parameters, slopes, come into a proper range. Two sets of tests can be made after parameter estimation: one is statistical tests and, the other, logical tests. The latter is usually neglected in literature. First, failing to pass the statistical tests indicates that it has some problems like incomplete or redundant variables, inaccurate parameter values, and so on. If so, the modeling process should be reconsidered. Second, failing to pass the logical tests indicates some structural problem. In this instance, the model cannot explain the situation at present and cannot predict the tread of development in the future. Statistical tests bear conventional procedure. However, the logical tests must be made by means of mathematical reasoning and numerical analyses.

Parameter estimation and model selection
The above-stated model on rural-urban interaction is an equation system coming from empirical analysis. One of the general forms of urbanization dynamics models can be expressed as This is in fact the urbanization model of United Nations [23], in which a, b, c, ϕ, ψ, and ω are parameters. In order to make statistical analysis based on the observational data, we must Reinterpreting the Origin of Bifurcation and Chaos by Urbanization Dynamics http://dx.doi.org/10.5772/intechopen.71035 discretize differential equations, Eq. (14), into difference expressions. As a result, the analytical process based on continuous dynamics is converted into the process based on discrete dynamics. Given that Δx/Δt~dx/dt, in which the time difference is Δt = 10. The independent variables include r(t), u(t), and r(t)*u(t)/[r(t) + u(t)], and the dependent variables are Δr(t)/Δt and Δu(t)/Δt, respectively. The model can be fitted to the American census data of rural and urban population. A multivariate stepwise regression analysis based on the least squares calculation gives the following model: which corresponds to Eq. (1). Clearly, the model parameters φ = ψ = ω = 0. Eq. (15) is a pair of difference equations. Given a significance level of α = 0.01, the important statistics such as F value, T value, variance inflation factor (VIF) value, and Durbin-Watson (DW) value can pass the common tests. In theory, as indicated above, we have b = c. However, in the empirical modeling, the two parameters are not equal to one another. The main reasons are as below. First, America is not a truly closed system. It has mass foreign migration. Second, the natural growth of urban population relies heavily on the rural-urban interaction. The latter reason seems to be more important than the former one. All things considered, as a special case of the United Nations model, Eq. (15) can describe the rural and urban population migration and transition of America in the recent 200 years in a better way.
To examine the relationship between the one-dimensional map and the two-dimensional mapping of urbanization, we can investigate the US urbanization curve. According to Eq. (9), the level of urbanization should follow the logistic curve. It is easy to calculate the urbanization ratio using the data in Table 1. For convenience, we set time dummy t = year-1790. A least squares computation involving the percentage urban data gives the following results: The goodness of fit is about R 2 = 0.9839. According to Eq. (16), the intrinsic growth rate is about k = 0.02238. On the other hand, according to Eq. (15), the intrinsic growth rate has two estimated values: the first is k 1 = b-a ≈ 0.03615-0.02584 = 0.01031, and the second is k 2 = da ≈ 0.05044-0.02584 = 0.02460. The number comes between 0.01031 and 0.02460. This suggests that the parameter value based on Eq. (15) is consistent with that based on Eq. (16). The subtle difference between different estimated results can be attributed to three reasons, that is, nonclosed geographical region, imprecise observational data, and the computation error stemming from the conversion from continuous function to discrete equation.
As a reference, the American rural and urban data can be fitted to the predator-prey interaction model. The independent variables include r(t), u(t), and r(t)*u(t), while the dependent variables are Δu(t)/Δt and Δr(t)/Δt, respectively. The multivariate stepwise regression based on the OLS method yields an unacceptable result [27]. If the statistical standard for modeling is lowered, then the US urbanization can be described by the Keyfitz-Rogers model [28,29]. Unfortunately, this model bears two vital deficiencies and is not acceptable for urbanization analysis [27]. All in all, both the linear model proposed by Keyfitz and Rogers and the nonlinear model presented by Lotka and Volterra are inferior to the special case of the United Nations model, Eq. (1).

Numerical experiment
As a complement analysis, the US census data of urban, rural, and total population as well as the level of urbanization can be generated using the rural-urban interaction model. A comparison between the simulation value and observed data shows the effect of urban modeling. The numerical simulation results are based on Eq. (15) and are displayed in Figures 5 and 6, respectively. Clearly, the change of the urban and total population takes on of the sigmoid curves, while the rural population takes on a unimodal curve ( Figure 5). What is more, the urbanization level is also an S-shaped curve, which can be described with the logistic function ( Figure 6). The changing trends of four types of curves based on the numerical simulation are supported by the observation data from the real world [4,27]. In the model, the capacity parameter of the urbanization level is evaluated as 100%, and this does not accord with reality of urban evolution. Nevertheless, the basic characters of the rural and urban development can be brought to light by Eq. (15). Anyway, there is no logical contradiction in the results from the numerical computation based on the rural-urban mapping.
So far, we have finished the building work of the model of urbanization based on the population observation in the real world. To sum up, the calculation results lend empirical support to the theoretical models and relations. First, the rural-urban population interaction model is testified, at least for a number of developed countries. The American model of rural-urban population interaction can be expressed by Eq. (1). This is the experimental foundation of theoretical analysis Figure 5. The predicted curves of the US rural, urban, and total population based on the two-dimension mapping of rural-urban interaction (Notes: the numerical iteration is fulfilled by Eq. (15), and the population unit is 10,000 persons. See [27]).
of discrete urbanization dynamics. Second, the relationship between the one-dimensional map of logistic growth and the two-dimensional map of rural-urban interaction is verified. By using the system of rural-urban interaction models, we can produce the logistic curve of urbanization. What is more, the curves of urban population, rural population, and total population are empirically acceptable. In the following section, I will discuss the related questions about bifurcation, chaos, complexity, and scaling law from the theoretical angle of view.

Generalization and supposition
According to the theoretical derivation, numerical experiments, and empirical analysis, an inference can be reached that chaos originates from nonlinear interaction between two coupling elements. The reasons are as below. First, a one-dimensional logistic map is actually based on a two-dimensional interaction map between two populations. Second, both the onedimensional map and the two-dimensional map processes can create the same patterns of bifurcation and chaos. Further, the theoretical findings can be generalized to the other scientific fields. Where dynamical behaviors are concerned, urban systems bear analogy with ecosystems [21]. Both the logistic equation and the predator-prey interaction model coming from ecology and can be applied to urban studies [24]. The predator-prey system can be modeled by different mathematical expressions, which can produce period-doubling bifurcation and chaos [9,[30][31][32]. On the one hand, the bifurcation and chaos proceeding from the twodimension mapping of urbanization dynamics remind us of the complicated behaviors shown by the one-dimension logistic mapping of insect population. On the other, the model of ruralurban interaction reminds us of the Lotka-Volterra model for the predator-prey interaction [25,26]. Therefore, the conclusions drawn from urban studies may be generalized to ecological  15), and the capacity parameter is 1. The curve is identical in shape to that of logistic growth indicated by Eq. (16). See [27]).
field and vice versa ( Table 2). A speculation is that the logistic growth in ecology can be interpreted by the two-population interaction, and the Lotka-Volterra model can be revised as below [8]: in which x(t) and y(t) denote the numbers of prey and predators at time t, respectively. The (1) can be treated as a special case of Eq. (17). Suppose that the percentage of predator population is defined by Thus, we can derive a logistic equation from Eqs. (17) and (18) as follows: Discretizing Eq. (19) yields a one-dimension mapping of logistic growth as below: where the parameter k~c-a-d. Both the one-dimension mapping based on Eq. (19) and the twodimension mapping based on Eq. (17) can create the same complicated dynamics as those displayed in Figure 1. This implies that the two-population interaction leads to the logistic growth, periodic oscillations, and chaotic behavior in ecosystems.

Model Dynamical equation Urban system Ecosystem
Allometric growth dx t ð Þ=dt ¼ ax t ð Þ dy t ð Þ=dt ¼ by t ð Þ Allometric scaling relations Two-population competition The rural-urban interaction The predator-prey interaction The rural-urban interaction and logistic growth Two-population competition and predator-prey interaction Notes: In these equations, a, b, c, and d are parameters. The equations of allometric growth suggest simplicity, while the two-population interaction models indicate complexity. A conjecture is that the logistic growth of population in ecology is just an approximate expression. It is the ratio of the predator population to the total population rather than the predator population itself that follows the law of logistic growth. Using the two-dimension mapping based on Eq. (17), we can carry out a numerical simulation experiment. The results show that if the percentage of predator population z(t) takes on a logistic growth, the predator population y(t) will grow according to an J-shaped curve in form (Figure 7). However, the latter is a quadratic or fractional logistic growth rather than the conventional logistic growth. What is more, the oscillations of population and the total population can mirror the perioddoubling bifurcation and chaos of percentage population. All in all, the generalized predatorprey interaction can account for more ecological phenomena than the classical Lotka-Volterra model does [8]. Anyway, the studies on urban chaos can help us understand John Holland's question. After discussing the Lotka-Volterra model, Holland [33] said: "In the long run, extensions of such models should help us understand why predator-prey interaction exhibit strong oscillations, whereas the interactions that form a city are typically more stable."

Scaling in bifurcation diagrams
Chaos and fractals are often placed in the same category in literature, although there is no essential correlation between them. A fractal is a hierarchy with cascade structure, which can be testified by urban systems. In fact, a period-doubling bifurcation diagram contains selfsimilar hierarchy. So, the period-doubling bifurcation route to chaos of urbanization dynamics can be compared with the hierarchical structure of cities. The general character of varied bifurcation diagrams can be reflected by Feigenbaum's number, which is a universal constant found by Feigenbaum [34]. This constant can also be figured out through the rural-urban interaction mapping. Based on a bifurcation diagram, we can draw the tent map [35] (Figure 8). If we give up the hypothesis of regional closure, the parameter equation β = γ in Eq. (11) will break. Then more complex and plentiful dynamics can be revealed by the rural-urban interaction mapping.
The period-doubling bifurcation process of urbanization and the cascade structure of systems of cities share the same hierarchical scaling. The bifurcation can be described with three exponential functions such as Figure 8. Tent map: From steady state to chaos (the initial value is L 0 = 0.01) (Note: the graph of tent map is also termed "spider diagram," which can be seen in literature such as ref. [35]. The diagrams are created by using the two-dimensional rural-urban interaction map based on Eq. (11). These subplots correspond to the subplots in Figures 1 and 3).
where m denotes the order of hierarchy of bifurcation, N m refers to period number (or bifurcation number), L m is the range for the stable periodicity, and W m is the span between two bifurcation points of order m. As for the parameters, N 1 = 1, L 1 and W 1 are constants, r = 2, δ ≈ 4.6692, and a ≈ 2.5029 [34,36,37]. In fact, Eq. (22) represents "the bifurcation-rate scaling law" and Eq. (23) "the fork-width scaling law" [38]. Accordingly, Eq. (21) represents "periodnumber scaling law." These what is called scaling laws are linear scaling laws, but they can be transformed into nonlinear scaling laws, i.e., power laws [4]. From Eqs. (22) and (23), it follows an allometric scaling relation between the periodical range (L m ) and the fork width (W m ), and the expression is where the proportionality constant is μ = L 1 W 1 Àb and b denotes a scaling exponent as below: b ¼ ln δ ln a ≈ 1:6796: The physical meaning of this number is not yet clear for the time being and remains to be brought to light in future studies.
The three exponential equations reflect the universal cascade structure of nature and society. An analogy can be drawn between the cascade structure of the bifurcation diagram and the hierarchical structure of urban systems ( Table 3). The scaling laws behind the period-doubling bifurcation can be employed to describe the nonlinear process of urbanization, and the variants of the scaling laws can be adopted to characterize the cascade structure of a hierarchy of cities [10][11][12]39]. Moreover, the allometric scaling relation, Eq. (24), bears an analogy with the fractal relation between urban area and population. The allometric growth law asserts that the rate of relative growth of an organ is a constant fraction of the rate of relative growth of the total organism [40][41][42]. In urban studies, the allometric scaling law can be utilized to describe the measure relation between the urban area (A m ) of a city and its population (P m ) in the urbanized area [4,16,42]. The similarity between urban scaling and bifurcation scaling lends further support to the inference that urban evolution falls between order and chaos [43].

Linear scaling law Period-doubling bifurcation Hierarchy of cities
The first law-number law The second law-length/size law The third law-width/area law Notes: (1) The scaling laws of hierarchy of cities are illuminated by [4]. (2) The period-doubling bifurcation in this work comes from the two-dimension mapping based on the rural-urban interaction model, which differs from the onedimension logistic mapping in form. Table 3. A comparison between the linear scaling laws of period-doubling bifurcation and the exponential laws of hierarchy of cities.

Dynamics of fractal dimension evolution of urban growth
The nonlinear dynamics of urbanization corresponds to the complex dynamics of urban growth and morphology. Urban growth can be measured with the time series of fractal dimension of urban form. The common fractal dimension can be obtained by box-counting method.
In theory, the box dimension of urban form ranges from 0 to 2. However, in practice, the box dimension always comes between 1 and 2. Boltzmann's equation can be employed to describe the fractal dimension growth of cities [13]. In fact, Boltzmann's equation was used to model urban population evolution by Benguigui et al. [44]. Urban population is associated with urban form and urbanization. The Boltzmann model of fractal dimension evolution is as follows: where D(t) refers to the fractal dimension of urban form in time of t; D (0) to the fractal dimension in the initial time/year; D max ≤ 2 to upper limit of fractal dimension, i.e., the capacity of fractal dimension; D min ≥ 0 to the lower limit of fractal dimension; p is a scaling parameter associated with the initial growth rate k; and t 0 a temporal translational parameter indicative of a critical time, when the rate of fractal dimension growth indicating city growth reaches its peak. The scale and scaling parameters can be, respectively, defined by p = 1/k, t 0 = ln[(D max -D (0) )/(D (0) -D min )] p . For the normalized variable of fractal dimension, Eq. (26) can be reexpressed as a logistic function: where D (0) * = (D (0) -D min )/(D max -D min ) denotes the normalized result of D (0) , the original value of fractal dimension. Empirically, Eqs. (26) and (27) can be supported and thus validated by the dataset of London from Batty and Longley [16], the datasets of Tel Aviv from Benguigui et al. [45], and the dataset of Baltimore from Shen [46]. The derivative of Eq. (27) is just the logistic equation: which is actually based on the normalized fractal dimension. Without loss of generality, let the time interval Δt = 1. Thus, discretizing Eq. (28) yields a one-dimensional map such as Defining D t * = (1 + k)x t /k, we can transform Eq. (29) into the following form: where x t is the substitute of D t * and μ = k + 1 is a growth rate parameter. Eq. (30) is just the wellknown logistic map [1]. If the fractal dimension of urban form can be fitted to Boltzmann's equation, it implies that urban evolution can be associated with spatial chaotic dynamics.
The process of urban growth is a dynamic process of urban space filling. An urban region falls into two parts: filled space and unfilled space. We can define a spatial filled-unfilled ratio (FUR) for urban growth [13], that is: Thus we have where U refers to the filled space area with various buildings (space-filling area), measured by the pixel number of built-up land on digital maps, and V to the unfilled space area without any construction or artificial structures (space-saving area). Thus the total space of urbanized region is S = U + V. Obviously, the higher the O value is, the higher the degree of urban spatial filling will be. The normalized fractal dimension can be termed level of space filling (SFL) of cities, implying the degree of spatial replacement.
Based on a digital map with given resolution, the filled space can be measured with the pixels indicating urban and rural built-up area such as structures, outbuildings, and service areas. In contrast, the unfilled space is the complement of the filled space of built-up area. On the digital map, the unfilled space is just the blank space of an urban figure. If a region is extensively developed and is already occupied by various urban infrastructures and superstructures, it is transformed, and the unfilled space is replaced by filled space. This spatial replacement dynamics can be described by a pair of differential equations: where α, β, and λ are parameters. This implies that the growth rate of filled space, dU(t)/dt, is proportional to the size of filled space, U(t), and the coupling between filled and unfilled space, but not directly related to unfilled space size; the growth rate of unfilled space, dV(t)/dt, is proportional to the size of unfilled space, V(t), and the coupling between unfilled and filled space, but not directly related to filled space size. From Eq. (33), we can derive Eq. (28). Discretizing Eq. (33) yields a two-dimensional map of urban growth, which can be used to created periodic oscillation chaos similar to the patterns shown in Figure 1 [13].

Replacement dynamics
The logistic growth model and the rural-urban interaction model can be employed to develop the theory of replacement dynamics. Dynamical replacement is one of the ubiquitous general empirical observations across the individual sciences, which cannot be understood in the set of references developed within the certain scientific domain. We can find the replacement processes associated with competition everywhere in nature and society. The theory of replacement dynamics should be developed in the interdisciplinary perspective. It deals with the replacement of one activity by another. One typical substitution is the replacement of old technology by new; another typical substitution is the replacement of rural population by urban population. Urbanization is a process of population replacement, that is, the urban population substitutes for the rural population [47,48]. The components in a self-organized system, generally speaking, can be distributed into two classes, and the process of a system's evolution is a process of discarding one kind of component in favor of another kind of component. This process is termed "replacement" [13,14]. For example, the population in a geographical region can be divided into urban population and rural population, and urbanization is a process of rural-urban replacement of population [48]; the technologies can be divided into new ones and old ones, and technical innovation is a process of new-old technology replacement [49,50]. In fact, people can be divided into the rich and the poor, the geographical space can be divided into natural space and human space, and so on. Where there are self-organized systems, there is evolution, and where there is evolution, there is replacement. Replacement results from competition and results in evolution. Replacement analysis is a good approach to understanding complex systems and complexity.
The basic and simplest mathematical model of replacement is the logistic function, which can be employed to describe the processes of growth and conversion. Besides, other sigmoid functions such as the quadratic logistic function and Boltzmann's equation may be adopted to model the replacement dynamics. A number of mathematical methods such as allometric scaling can be applied to analyzing various types of replacement. In fact, the allometric scaling can be used to analyze the relationships between the one thing/group (e.g., urban population) and another thing/group (e.g., rural population). A replacement process is always associated with the nonlinear dynamics described by two-group interaction model. The discrete expression of the nonlinear differential equation of replacement is a one-dimensional map, which is equivalent to a two-dimensional map. The maps can generate various simple and complex behaviors including S-shaped growth, periodic oscillations, and chaos. If the rate of replacement is lower, the growth curve is a sigmoid curve. However, if the replacement rate is too high, periodic oscillations or even chaos will arise. This suggests, no matter what kind of replacement it is-virtuous substitution or vicious substitution-the rate of replacement should be befittingly controlled. Otherwise, catastrophic events may take place, and the system will likely fall apart. The studies on the replacement dynamics are revealing for us to understand the evolution in nature and society, and the relationship between the one-dimension map and the two-dimension map is revealing for our understanding of the replacement dynamics.

Conclusions
Researching the origin and essence of bifurcation and chaos in urbanization process offers a new way of looking at complicated dynamics of simple systems. The pattern of phase space cannot be revealed by the one-dimension mapping diagram based on ecological systems, but it can be displayed by the two-dimension mapping diagram based on the rural-urban population migration and transition. This suggests that urban evolution is a good window for examining bifurcation and chaos. Moreover, the similarity between urban dynamics and ecological dynamics will inspire us to explore the implicit substance of natural laws. By the study of urbanization dynamics, we can obtain three aspects of new knowledge about bifurcation and chaos. First, it is interaction rather than the intrinsic randomicity of dynamic systems that leads to bifurcation and chaos. Period-doubling bifurcation and chaos used to be regarded as inherent randomness of determinate systems due to the complicated behaviors of the one-dimension logistic mapping. The people with this viewpoint ignore the following fact: the logistic growth is always based on two-population interaction. However, because of the absence of effective measurement linking the logistic function and the two-population interaction model, the relationships between chaos and interaction cannot be revealed in ecological fields. Second, the chaotic behaviors of the logistic model do not indicate a chaotic attractor, and the relationship between chaos and fractals is scaling. A strange attractor with fractal structure in the phase space used to be treated as a typical sign of chaos. The phase portrait of the logistic growth cannot be demonstrated by the one-dimension mapping. The two-dimension mapping based on rural-urban interaction can be employed to illustrate the phase space of the logistic process. The result shows that the whole trajectory fails to converge into a limited area. No strange attractor or even no fractal structure can be found in the phase portrait of the two-population mapping. However, both fractal structure and the route from bifurcation to chaos can be characterized by hierarchical scaling law. Third, the predator-prey interaction model can be developed to interpret the logistic growth and sigmoid curves. By analogy, we can infer that the predator-prey interaction causes the complicated behaviors of the logistic process in ecological field. In fact, the classical Lotka-Volterra model can be restructured by referring to the expression of the rural-urban interaction model. Consequently, we can get a normalized predator-prey interaction model. Using the revised predator-prey model, we can derive the logistic function for population growth. Finally, where urban geography is concerned, the models of urbanization dynamics can be generalized to describe the spatial dynamics of urban morphology by means of fractal dimension growth. Moreover, both the models of urbanization and urban form evolution can be applied to developing the theory of spatial dynamics of replacement.