Second-order metropolitan urban phase transitions

Chaos, Solitons & Fractals 48 (2013) 22–31

Contents lists available at SciVerse ScienceDirect

Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Second-order metropolitan urban phase transitions Roberto Murcio ⇑, Antonio Sosa-Herrera, Suemi Rodriguez-Romo Theoretical Research Center, Mexico’s National Autonomous University (UNAM), Estado de Mexico 54701, Mexico

a r t i c l e

i n f o

Article history: Received 7 August 2012 Accepted 3 January 2013 Available online 31 January 2013

a b s t r a c t The morphology evolution of Metropolitan Urban Areas constituted by different Central Business Districts is studied in this paper. For this matter, we propose a stochastic model which combines an initial percolation setting followed by a diffusion-limited aggregation mechanism. Our model mimics better than either case (percolation or diffusion-limited aggregation) the Metropolitan Urban Areas formation progress. We argue that the Metropolitan Urban Areas case introduced in this paper, grows in such a way that undergoes a non-equilibrium second-order phase transition during this process. This conclusion is supported by a fractal dimension and configurational entropy analysis, as well as by studying an empirical case.  2013 Elsevier Ltd. All rights reserved.

1. Introduction Many structures in natural and artificial systems arise as a result of complex processes. In some cases it is possible to infer the process from the final structure. Besides, the research of human settlements has attracted much attention recently since its results are useful in sciencebased social and political decisions and efficient use of urban resources. The way people establish in a number of settlements is an obvious complex process; it involves several individual decisions such as; government policy, geographical restrictions, historical events (foundations and disasters, for instance), among others. All of these factors shape cities as time goes by; for instance the Metropolitan Urban Area (MUA) studied in this paper (called hereafter MUA-Mexico), built by three different initial Central Business Districts (CBDs); Valle de Mexico Metropolitan Area, Puebla Metropolitan Area and Toluca Metropolitan Area, is one of the largest metropolitan area of the world. We consider some features of this case as parameters in the simulation performed in this paper to test our model, and the empirical data available from the evolution of this Metropolitan Urban Area to be compared with the model. ⇑ Corresponding author. Tel.: +52 15532224044. E-mail address: [email protected] (R. Murcio). 0960-0779/$ - see front matter  2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chaos.2013.01.001

Recent studies on human agglomeration suggest the existence of a universal mechanism governing the population growth at different scales [1]. In this paper we intend to understand a leading process which drives to MUAMexico urban changes all the way through; from foundation dates (by late XV or early XVI century) until nowadays. Our studies have revealed that the topography of the land does not play an important role with respect to historical and economic conditions accompanying the development of MUA-Mexico. Some regularity in our simulations indicates that a common generic mechanism underlies the urbanization processes of our case. To understand this generic mechanism we propose here a model, perform experiments, compare them with empirical data and provide results and conclusions. Our model is based on two well-known fractal growth processes present in different phenomena in nature; diffusion-limited aggregation [2] and percolation [3]. Experiments in our model are performed with the final goal of studying near non-equilibrium phase transitions in complex social phenomena by using a method inspired on MUA-Mexico and based on fractal analysis and the configurational entropy of fractal aggregates, in our case. It is clear that different aggregation states in nature are indeed different phases, but not the other way around, since we can find a number of phases with the same aggregation state, for example the ice is solid, and has a number

23

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

of different phases. Even more, rearrangement of structure in solids and liquids produce different phases; for example there are organic materials with high molecular weight and a long stretched form of the molecules where normal liquid phase and crystal liquid phases occur. Usually, due to the number of independent thermodynamic variables (via the Gibbs’s and Gibbs–Duhem relations) a two dimensional space or phase diagram is good enough to contain all the information related to the phases of a substance; areas are one phase states and lines are phase transitions (two phases coexistence of phase transitions). It is important to remark that this is fully applicable in thermodynamics, therefore, stationary regimes, and this is not the case studied in this paper. However, given the conditions of the matter studied in this paper, we consider our approach good enough to provide some insight over important features in the field. For example, in nature some well-known examples of phase transitions are: the transition of no null viscosity of the 4Helium liquid to superfluid 4Helium, order–disorder phase transitions; the transition to ferroelectricity in certain materials and rearrangements of the orientation in crystal lattices, among others. Even more, there are a number of cases where phase diagrams are constructed with interesting results in games [4,5], human mobility [6], traffic flow [7], urban evolution [8] and granular media [9], among others. Phase transitions are classified in physics basically in two groups: first and higher order transitions. Phase transitions which are connected with an entropy discontinuity are called discontinuous or phase transitions of first order. On the other hand, phase transitions where the entropy is continuous are called continuous or of second or higher order. Like evaporation, melting and sublimation are also phase transitions of first order. Note, however, that some processes which are melting are not phase transitions in our sense. For instance, glass becomes viscous as it is heated up, and then becomes liquid. This corresponds to a continuous change of viscosity, in which no entropy discontinuities appear. Furthermore, glass is not crystalline in the solid state, but shows an amorphous structure which does not change suddenly under heating. Thus, even at room, temperature glass has actually to be considered to be a liquid (however, a liquid with an extremely large viscosity). The degree of order of glass does not change, in contrast to the solid–liquid phase transition water.

2. The model At small geographical scales (1:1000,000–1:500,000), the urban morphologies observed for Metropolitan Urban Area all over the world are not emerging and evolving only as a result of the hundreds of human activities, but as a result of the complex relationships between a few macro variables [10–12]. The model introduced in this paper is based in the work already done in references [13,14]. Five macro variables or urban drivers are proposed as the origin of a particular Metropolitan Urban Area [13]. These drivers can briefly be defined as follows:

(a) Random events ei. We assume that individual human decisions and particular activities are of such nature that at the geographic scale they are merely random and cannot be predicted at scales we are searching for. This random condition could be seen as the equivalent of adding some ‘‘noise’’ to a physical system. We take into account this driver in our MUA-Mexico model, thus we state from the beginning that our approach is not at equilibrium. (b) Historical accidents hi. Particular historical events such as cities foundations or its corresponding massive destructions are considered. We also consider this particular driver as can be seen in Table 1. This is done by choosing a particular set of initial sites which form three different seeds (one per CBD’s; Mexico, Puebla and Toluca Cities) and giving them a potential which is kept constant all the way through the simulation and reflects the relative socio economical importance of each urban area; 20 for Mexico, 10 for Toluca and 4.8 for Puebla Cities. (c) Physical constraint ci. The physical terrain in which cities are developing is clearly a physical restriction for its growth. There are areas where it is at present times impossible to build urban structures such as roads or bridges. We already mentioned that this driver does not influence significantly our morphological results, perhaps because our MUA-Mexico case is located on a valley. This driver is not considered by our model. (d) Natural advantage ni. These advantages, related to the natural resources are, historically, probably the main factor over human group’s decision to settle somewhere. We do not consider this driver since natural resources are not an issue in the case we are studying; the evolution of the three metropolitan areas considered in this paper is mainly driven by other factors, related to economy. (e) Comparative advantage eci. Economy is present in our lives at daily bases. Land prices, access to basic services, security, amenities, parks, schools, among others are considered by people and governments as important issues to urbanize a particular area, this variable is an economic equivalent of natural advantages. Our model is almost completely based on this driver since from foundation, MUA-Mexico evolution has been mostly driven by comparative advantages. Summarizing, three urban drivers changes are considered in our model; random events ei, historical accidents hi and comparative advantage eci. Even though it is clear

Table 1 Initial layout and potential for all the historical accidents considered in our model at time t = 0. MUA

Initial potential

Initial area covered (lattice units)

VMMA PMA TMA

20.0 10.0 4.8

4 3 3

24

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

that the physical terrain and natural resources influence Metropolitan Urban Area spatial morphology and both features should be included in a realistic model, the model presented in this paper reproduces morphologies which are very similar to the actual MUA-Mexico studied, with no inclusion of physical constraints (ci) and natural advantages (ni). This fact leads us to believe that a complex interaction of parameters is involved in the real case and it is exposed in our model. Let us set a 320  240 regular grid and a potential for population growth Pi(t) in the lattice site i at iteration (time) t as a reference state of the MUA-Mexico growing process [15,16]. During evolution this potential undergoes a positive feedback based on a return to scale as a typical urban growth process [16]. The change in potential for urban growth is assumed as a proportion of the existing level of potential in the lattice site i at time t of the grid as follows:

DPi ðtÞ ¼ kPi ðtÞ;

k>0

ð1Þ

Here k is the growth rate. In our construction we use Arthur’s model [17] such that Eq. (1) becomes

Pi ðt þ 1Þ ¼ Pi ðtÞ þ DP i ðtÞ ¼ ð1 þ kÞPi ðtÞ;

k>0

ð2Þ

Eq. (2) can be rewritten in terms of the grid cell potential average in a restricted neighborhood around each grid cell. In particular, we use in this paper Batty’s approach [13]; namely the Von Neumann neighborhood Xi [18] composed of the four adjacent cells: north, south, east and west of a particular cell. The growth rate is set to zero, namely k = 0, then Arthur’s model establish that

Pi ðt þ 1Þ ¼

X Pj ðtÞ j¼Xi

ð3Þ

5

Finally, white noise is added at each period t, to represent the urban driver ei. Thus we have

Pi ðt þ 1Þ ¼

X Pj ðtÞ j¼Xi

5

þ ei ðtÞ with t 2 Zþ ;

since noise is continually being reintroduced into the system. It is clear that the final form of Eq. (4) depends basically on the function Pi(t). Classically the potential measure in social physics depends on an index of economic proximity between zones. Such indexes usually compare the economics of two regions and establish how far are one from each other in economic terms, so potential measures based on those indexes indicates precisely the economic potential for a specific zone or group [19]. As in reference [13] we set Pi(0) = ei, for any i in the grid, so every location in our model has an initial value of 1 or 1. Therefore, we consider a given threshold U to regulate the way Pi(t) is converted to an actual city or urban area, or is not. Metropolitan Urban Area distributions are obtained by iterating equation (4) and translate potential into actual population through U. A variable Di(t) is defined to point out whether or not at position i and time t, an urban configuration (or Metropolitan Urban Area) has undergone some growth. In essence, growth will take place when potential reaches a certain threshold U, once this happens growth persists in the future. Formally we mean that

ð4Þ

Where ei, for simplicity, is set either 1 or 1 with equal probability. The inclusion of random events as an urban driver destroys any tendency towards the equilibrium

If P i ðtÞ > U and Di ðt  1Þ ¼ 0 then Di ðtÞ ¼ 1;

otherwise Di ðtÞ ¼ 0

ð5Þ

As we state above, all locations on the grid have an initial value of 1 or 1, except for some fixed positions representing urban historical accidents. Table 1 shows the initial data chosen. The particular values shown in Table 1 are chosen to represent the spatial morphology of MUA-Mexico. Historical accidents are represented in our model through the position of three clusters, an initial potential and an initial number of urban settlements for each metropolitan area (Valle de Mexico Metropolitan Area, Toluca Metropolitan Area and Puebla Metropolitan Area), a so called initial conditions of the model. The initial potential and number of settlements are a function of the relative economic and cultural importance of the three areas by the date of their foundation, approximately 1500 AC [20]. The positions on the grid for each settlement and the Euclidean distance between them are also necessarily in

Fig. 1. Three typical configurations obtained by our model for t = 500 and thresholds: (a) U = 2.1, (b) U = 4.1, (c) U = 6.5. Actual growth is represented by the black structures in each figure. We can observe how sensitive to the threshold parameter our model is; from 4.1 to 6.5 practically all structures are lost and only three clusters are recognizable.

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

scale with the actual MUA-Mexico. Over this initial scenario we apply Eqs. (4) and (5). See Fig. 1 for instances of our MUA-Mexico model simulations. 3. Experiments At first, we use Eqs. (4) and (5) along with the initial conditions explained above, and then we construct a number of configurations under different thresholds U and different iterations t. In this paper, we select thresholds U 2 [0, 0.1, 0.2,. . ., 19.9] and maximum values tm 2 [10, 20, 30, . . ., 500]. Thus, at time t = 0 all grid positions, except the ones occupied by historical accidents, have a potential Pi(t) = ei, with ei taking values of 1 or 1. To represent the importance of the CBD’s three settlements (Valle de Mexico Metropolitan Area, Toluca Metropolitan Area and Puebla Metropolitan Area) across time, the initial potentials of historical accidents are never changed at any time t, so Eq. (4) is never applied to the cells occupied by the initial conditions set in Table 1 for the MUA-Mexico instance studied in this paper. Up to this point, we have merely used a percolation approach to represent Metropolitan Urban Area and the economic potential they have per site in the urban grid. Following reference [21], we introduce into our model urban migrations attributes through a DLA approach. When our simulation reaches the point t = 0.1tm, we initiate a DLA mechanism over any configuration that is emerging at that particular time t. As in DLA, we define a launching radio which in our case is established at 0.75 times the size of the simulation radio, calculated from the grid’s center to the farthest unit from this geometrical center at iteration t. We replaced the regular random walk in the DLA algorithm with a directed random walk. Normally, any location

25

over the grid has the same probability to be visited by a particle at motion. In our model, we modified this probability in a way that grid zones with greater potential Pi(t) have greater probability to be visited. In a way, we reflect the tendency that regions in urban areas with greater economic potential would attract more people than the ones with smaller economic potentials. We create an ensemble of 9950 different configurations over which we perform 100 runs per configuration, so we can derive a robust statistics in our model. To learn from our MUA-Mexico evolution model, simulations are performed to calculate for each configuration the fractal dimension using the Box counting method (Db) and its configurational entropy H(N) [22,23] by an algorithm based on separation over classes through a Hamming distance [24]. The box counting method to obtain Db is a well-known algorithm. In Section 5 we elaborate on the algorithm to obtain H(N) parameterized by t and U. 4. Box counting fractal dimension By measuring the fractal dimensions of urban patterns a quantitative description for actual urban morphologies can be obtained. The Box fractal dimension Db tells us how much area is occupied by the urban structures, so we calculate Db for all the thresholds U used by our model for different iterations t. In Fig. 2 we show some of these dimensions. Below certain U depending on the t parameter, smooth monotonically decreasing function Db = f(U) is obtained and above U = 8 an oscillating decreasing Db = f(U) is obtained. Please note that, around 4.0 6 U 6 8.0 an interesting phenomenon appears; there is an interval of U where Db = f(U) is a fluctuating increasing function. As we can see, the decreasing tendency is reverted between U = 3.8 and U = 3.9; while the smooth increasing

Fig. 2. Fractal dimension Db (calculated by the box-counting method) as a function of the threshold U for t values of 10, 100, 200, 300 and 500. As threshold begins to increase, dimension follows a decreasing pattern, until a particular point in which suddenly increase its value and then experiment an oscillating pattern and decreasing again. This behavior is independent of t and suggests that larger clusters, as threshold goes up, begin to absorb satellite smaller clusters, compacting the structures generated.

26

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

Fig. 3. Configurations obtained for t = 500, and thresholds (a) U = 3.8; (b) U = 3.9; (c) U = 4.3; (d) U = 4.4. It can be noticed how with a slight increase in the threshold (0.1 units) some clusters disappear, while some others changes size. For example, when threshold goes from 3.8 to 3.9, the clusters inside the square boxes (a) are practically gone in (b); when threshold goes from 4.3 to 4.4, the cluster at the center of the image in (c) clearly gains the size showed enclosed in the triangle in (d). These effects are the physical manifestations of the behavior of function Db = f(U).

pattern is stopped when the threshold changes from 4.3 to 4.4, and then begin to oscillate with a decreasing tendency. In Fig. 3 we can see some examples of configurations with these parameters. We have perceived that, at a small scale, these oscillations in Db = f(U) appear every time some clusters (usually one of the three main central clusters or Metropolitan Areas) gain spatial size, while some of the smaller clusters disappear, mainly the ones far from our MUA-Mexico core. See for instance, Table 2 where the values for Db with t = 500 fixed and U between 3.6 and 4.6, and 5.6 are shown. Generally speaking, the dimension Db ? 0 as U ? 1, and in the limit U ? 1 this behavior is expected since as the threshold U goes larger, less urban growth takes place. For all t (at least the values tested here) the behavior described in Fig. 2 seems to be iteration scale invariant and follows some power function.

5. Configuration entropy based on separation of classes using the Hamming distance In this section we describe a computationally efficient algorithm introduced in this paper to calculate the config-

Table 2 Some examples of fluctuations in H(N), for t = 500 fixed and U between 3.6 and 5.6.

U

H(N)

3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 5.6

0.3086 0.2383 0.2172 0.2188 0.2307 0.2323 0.2696 0.2844 0.2801 0.2821 0.2774 0.2576

urational entropy of aggregates, for instance the ones shown in Fig. 1 obtained as outputs of the simulation of our MUA-Mexico. The algorithm is based on the concept of entropy in the context of information theory introduced by Shannon [24]. The basic idea of our algorithm is to transform the aggregates generated by our model into a 2-dimensional bits sequence for an l  l aggregate with liner size l. With

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

this modification, we define a set A(U, t) which elements ai, i = 1, 2, . . ., l are aggregates represented as a l  l matrix of bits, in which the value of each position is set to 1, if the corresponding position at the aggregated is occupied, and 0 otherwise. Generally speaking the configurational entropy H(A) over a set A of n-bit sequences is defined following Shannon’s definition: n X HðAÞ ¼  pi log2 ðpi Þ

ð6Þ

i¼1

Here pi is the probability that the final spatial distribution of a particular aggregate gives as a result the configuration ai. One aspect to be considered in our case is the practical difficulty to calculate pi since we have a large number of possible configurations for the space of linear size l, namely 2ll. Each pi is going to be very low, so, even with the aid of high-performance computing, it is not possible to directly measure H(A) from the number of occurrences of ai, because at a practical scale, it is unlikely that ai occurs more than once. This is a straightforward consequence of the stochastic model we are proposing to study MUA-Mexico. To get around this problem, instead of calculating pi per configuration, a collection of state classes N = {nk:k = 1, . . ., k} is considered. Each class is a set of similar aggregates generated under the same threshold U and time t. There are many ways to define similarity among different geometric object [25]. In this paper, we consider a similarity measure based on the entropy itself, following an Entropy-based Fuzzy Clustering [26] procedure with a modified version of Eq. (6) that measures the similarity of a particular ai generated with a Uw and a tw with respect to all the others aj with the same Uw, tw. The similarity function used in this paper between elements ai and aj is defined as follows. Hmði;jÞ N

Sij ¼ ea

ð7Þ

where Hm(i, j) is the Hamming distance between two elements ai, aj 2 A(U, t), N is the average of the number of particles for both aggregates ai,aj; namely the average number of occupied space positions. The parameter a = 0.5 is for normalization purposes [27]. Notice that now we are restrict to only two possible outputs; either ai is similar or not to aj, up to certain concept of similarity previously defined by the precision we want to impose in our MUAMexico model. Shannon proposes Eq. (8) to calculate the entropy in the case of two possibilities in a stochastic process; with probabilities p and q = p  1, as follows:

E ¼ ðplogp þ qlogqÞ

ð8Þ

From Eqs. (7) and (8) we define the entropy Ei for and element ai 2 A(U, t) by Eq. (9). A formal derivation of Eq. (9) can be found in [26].

Ei ¼ 

N X

ðSij log2 Sij þ ð1  Sij Þlog2 ðSij ÞÞ

ð9Þ

j¼1;j–i

If two aggregates are completely equal or different to each other, then configurational entropy would be zero; if these aggregates are very similar or dissimilar between

27

them, configurational entropy would be low. From Eq. (9) we observe that the maximum configurational entropy Ei = 1 is obtained when Sij = 0.5. In our paper Ei is not the final goal of our analysis since this is giving information related to one element of A(U,t), recall that we need the entropy involving of the full set A(U, t). However it is a necessary step to find what we are searching for: H(N). We begin by constructing an algorithm to generate the classes N in order to calculate the configurational entropy H(N) parameterized by U and t. Classes N generating algorithm 1. Start with k = 1,X = A(U, t), nk = ;"k where A(U, t) = {ai: i = 1, . . ., n}. 2. Calculate the entropy Ei for each xi 2 X using Eq. (9). 3. Define a tolerance similarity value b. This value regulates the accuracy of the similarity measure chosen. Therefore sets the precision of the H(N) calculated below. 4. Select the ximin with the minimum entropy Ei. 5. Remove ximin and the xj 2 X with a similarity Sxi xj < b min from X. S S Then nk ¼ nk fximin g fxj remov edg. 6. If X – ;, then k = k + 1 and repeat from step 4. Otherwise, finish. Thus, we end up with a set of classes N = {nk} which contains the more representative elements for each configuration threshold U and time t. We calculate the configurational entropy H for N inspired on Eq. (6) as follows.

HðNÞ ¼ 

k X

qi log2 qi

ð10Þ

i¼1

Here qi ¼ Pknkknk k

kk

with knkk being the number of ele-

ments of the nk class. For instance, for t = 500 and b = 0.9 we obtain the values shown in Table 2. At this point it is obvious that b is another factor that influences the results (by accuracy) of the simulations of our MUA-Mexico model. (see Table 3) As expected, the determination of a set of configurational states N partitioned by classes {nk} is quite demanding from a computational point of view. At first, as already mentioned, we had to perform 100 experiments per analyzed configuration (in total 9950). The algorithm to obtain the set of similarity matrixes Si,j in Eq. (7) has a computational complexity of order O(n2) for each i,j since every element has to be compared with all the rest in order to be properly classified. So, to handle the amount of obtained statistical data to obtain conclusions with statistical significance in this paper, we have executed our simulations with parallel computation hardware and computing techniques [28]. To simulate our model we used a general purpose graphic processing unit GPU with 336 processing elements working along with a six-core central processing unit CPU. At each stage of the simulation, the potentials Pi for each cell were calculated in parallel by one processing element inside the GPU, as well as the state of the random walks involved in the DLA model. Evaluation of the

28

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

Table 3 In the second column the entropy H(N) values for t = 250, b = 0.9 are shown for different thresholds U (from 2 to 8) shown in the first column. The number of classes obtained for a particular set U, H(N), t = 250, b = 0.9 is introduced in the third column. Finally the probability (actually frequency) per class is provided in the remaining column.

U

H(N)

knkk

qi

2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8

2.892 3.244 3.789 3.481 3.619 2.309 1.684 0.770 3.5 4 0.65 0.353 0.566

8 10 14 12 14 6 5 3 2 3 2 2 2

0.16 0.06 0.06 0.06 0.1 0.36 0.6 0.83 0.9 0.8 0.83 0.93 0.86

0.06 0.1 0.06 0.06 0.1 0.2 0.16 0.13 0.1 0.16 0.16 0.06 0.13

0.13 0.06 0.06 0.06 0.06 0.16 0.13 0.03

0.06 0.1 0.06 0.03 0.16 0.1 0.03

0.16 0.06 0.06 0.06 0.03 0.13 0.06

0.16 0.1 0.06 0.13 0.1 0.03

0.06 0.06 0.06 0.13 0.1

0.16 0.13 0.06 0.1 0.03

0.16 0.06 0.06 0.06

0.13 0.06 0.13 0.06

0.1 0.066 0.033

0.06 0.06 0.06

0.1

0.066

0.03

0.03

0.03

similarity matrixes and box-counting dimensions were performed by using the six-core CPU. This process was carried on in a pipelined fashion in such way that the computational load was almost equally balanced between the CPU and the GPU. To achieve this, the program was implemented using hybrid multithreading OpenMP [29] and parallel CUDA [30] programming. The computational time required by all the simulations and analysis described here was about 200 h measured in wall-clock time.

6. Results and conclusions The simulation of our model clearly reflects the morphologies of very large urban settlements; this simulation resembles better the empirical data than the ones produce by percolation itself or DLA itself. Actually, our model reproduces, for some parameter setting, fairly well the goal to be studied in this paper for MUA-Mexico. We start by analyzing better the box-counting fractal dimension Db = f(U) behavior related to the cluster aggregation morphology for t = 500 and different threshold U values. The Db = f(U) oscillatory pattern above U = 4.0 and the behavior presented by this function in between U = 4.0 and U = 7.0 is certainly unusual, see Fig. 2. At the scale Fig. 2 is constructed, this effect resembles the behavior of a state equation similar to Van der Waals state equation [31] that, by the way, undergoes a gas–liquid phase transition. Besides, the oscillatory pattern that produces little ‘‘bumps’’ in Db for U > 4 appear every time some cluster gains size using some other cluster particles which disappear due to this process. To understand this behavior better we choose some particular experiments and compare the simulation resulting morphologies as can be seen in Fig. 3. We observe some square boxes framing clusters that in the next threshold (with a difference of 0.1) are not present, while triangles specify the clusters which have increased its spatial size. It is not possible to know if the clusters lost by some aggregates are indeed precisely the ones gained by the greater aggregates in the next threshold, but certainly this interchange of clusters in between

aggregates is affecting the potential field in the whole system, even if this is a small change. This effect recalls also the ‘‘evaporation’’ of particles from an aggregate and the simultaneous ‘‘deposition’’ of particles in another aggregate in an open system, recall that in our model we are launching particles as in DLA, at some point of the growth. In terms of the MUA-Mexico system is clear, besides the obvious conclusions, that our model describes a non-equilibrium balance in the middle of two phases with different order near to equilibrium, this is going to be clearer below, when we analyze the entropy evolution. However, in between the two near-to-equilibrium phases, we have a ‘‘kinetic’’ process similar to what happen in usual chemical reactions, where there are intermediate states. In our case this intermediate states correspond to urban settlements with different size treading urban inhabitants among themselves. We conclude that the box counting fractal dimension Db is very sensitive to the disappearing of aggregates. Theoretically, the configurational entropy H(N) as a function of the threshold U should be zero in both, the ordered and vacuum phases, since the former is constituted mostly by a single cluster of the size of the system and the latter is formed by a set of few small clusters of similar sizes that in the limit of U ? 1 is empty. H(N) should have a maximum at the phase transition where, by the way, the diversity of cluster sizes is also maximum. This is exactly our result as can be seen in Fig. 4 for the model here introduced; thus, a plateau develops as the phase transition takes place. Notice that H(N) = g(U) is a fluctuating function since we are simulation a stochastic model. Cluster-size entropy has been used in the study of problems such as percolation [32–35] and complex systems [32,35–37]. In these cases the cluster entropy shows a maximum at the percolation threshold, where a group of neighboring occupied sites forms a cluster that expands from one edge of the two-dimensional lattice to the opposite one. This resembles our case since we transit from a not percolating configuration (phase I) to a percolating one (phase II), however this is done throughout a wide critical region with elevated configurational randomness where percolation is neither guarantee nor excluded. It is

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

29

Fig. 4. Configurational entropy H(N) as a function of the threshold U. We show a number of configurational examples parameterized by U and t = 500 from the configuration ensemble with a particular H(N) value. For configurations that almost completely cover the study space (thumbnails at the left of the curve) as for configurations that hardly do it (thumbnails at the right of the curve) we obtain configurational entropy values very similar, as expected. At the maximum entropy interval (thresholds values between 3.0 and 4.5) we have many different possible configurations with just a minimum change in threshold value.

clear from Fig. 4 that in our model we have a continuous phase transition, in spite of the fluctuations or noise in the entropy for some threshold values. Even though the precise influence of the switching between percolation and DLA aggregation processes was not considered in this paper, we think that in our instance, the process is somehow influenced by the fact that our model starts as percolation and at a particular point of the process a DLA mechanisms is switched on and the clusters growth to resemble better the MUA-Mexico morphology. This fact might explain the width phase transition region (namely the large number of possible configuration within the H(N) versus U graph plateau) allowing us to conclude that our model for MUA-Mexico based in the Valle de Mexico Metropolitan Area, Toluca Metropolitan Area, Puebla Metropolitan Area foundation and economic power undergoes a non-equilibrium phase transition. We use empirical data for this MUA-Mexico system to support our ansatz. Actually the H(N) versus U graph plateau that can be seen in Fig. 4, corresponds to the state of maximum cluster disorder that occurs in our experimental ensemble. The width of the maximum H(N) entropy plateau may give a clue of the order of the transition; the narrower plateaus, transforming into a peak, indicates that the system crosses over from phase I to phase II in a sudden manner, which is typical of first-order-like transitions. Broader plateaus, instead, suggest that the system moves from one regime to the other in a smooth manner, this is the behavior expected in second-order-like phase transitions and this is exactly what we are obtaining. The cluster-size entropy has hardly been used in the analysis of phase transitions by other authors. Here, we show that the entropy H(N) is a valuable tool that can be utilized to understand the driving processes that determine morphologies of very large Metropolitan Urban Area,

as the instance shown in this paper. We can even show how the state equations Db = f(U) and H(N) = g(U) coincide in predicting a non-equilibrium phase transition for t = 500 (we remark that this effect is replied by other t values), see Fig. 5. In this Figure, Db = f(U) has an inflection point precisely located in the H(N) = g(U) plateau, we explain this as a confirmation of a non-equilibrium second-order phase transition of our MUA-Mexico model as an open system. By a coarse analysis, we also observe that the boxcounting dimension Db = f(U) starts by forming a plateau for small U values; it suddenly decreases in a small interval of U values, afterwards very slowly gets asymptotic to 0. On the other hand, H(N) = g(U) is increasing up to a maximum for small U values, afterwards remains nearly constant with few fluctuations and it decreases asymptotically to 0, with large fluctuations for large U values. Please notice that the H(N) = g(U) quasi-plateau matches (both occur by the same U value interval) with the abrupt ending of the Db = f(U) plateau. For us this fact remarks that, from the viewpoint of both approaches, fractal dimension and entropy, we have a phase transitions for intermediate U values; from an order-1 phase (small U values) to an order-2 phase (large U values). We think that even though our model produces Metropolitan Urban Area morphologies whose behavior is overwhelmingly driven by percolation theory, the DLA approach switch on at a particular time in the stochastic process performed is important in the non-equilibrium second-order phase transition being determinant in the plateau width of the H(N) = g(U) plateau and of the ‘‘evaporation’’ vs ‘‘condensation’’ balance happening in the Metropolitan Urban Area formation process. Our model assumes that two processes (percolation and DLA) lead to the observed MUA-Mexico morphology. We propose that percolation gives rise to the universality exhibited by our model simulation (scale invariance) while DLA establishes

30

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

Fig. 5. Overlap of state equations Db = f(U) and H(N) = g(U) for t = 500. The rapid decay of dimension Db coincides with the maximum entropy interval, and the very oscillating left part of H(N) coincides also with the oscillating section of Db = f(U).

the particular values of the corresponding laws and reflects small-scale ‘‘details’’ that are specific for a given Metropolitan Urban Area system. Finally, let us now face an important matter about this paper: is there any empirical evidence showing that any urban pattern in the real world evolves as the one modeled in this paper within the limitations of the available actual data at the scale studied in this paper? Of course our first hint is the MUA-Mexico, the very same case that inspired our model. It can be shown that most rearrangements of the structure in phase transitions can be described by a so-called order parameter. This order parameter should represent the main difference between the various phases; it is a

measure of the degree of order in the system. Looking for some sort of ‘‘order parameter in urban evolution’’ we study the evolution of urban density (population/urban area) of MUA-Mexico from 1895 until 2010, which are the full set of actual data available for this metropolitan urban area at the scale studied in this paper. We smooth the function urban density vs. year by a ‘‘Smoothing spline’’ technique [38] and found that this function clearly undergoes a change of concavity. Then, we calculate the second derivative where a clear discontinuity appears by the middle 90’s, see Fig. 6. It is clear that at a particular year a specific potential threshold would be established. From the MUA-Mexico actual data we conclude that our model might represent real cases and be coherent with

Fig. 6. Second derivative for the smoothing function obtained from investigate the urban density (population/urban area) of MUA-Mexico from 1895 to 2010. A change on the sign of the derivative is observed after 1990, fact that suggest the second-order phase transition (see text).

R. Murcio et al. / Chaos, Solitons & Fractals 48 (2013) 22–31

empirical data, however straightforward connection is not clear so far due, among others, to the available data at the scale used in this paper in the mankind history. Other order phase transitions at the scale used in this paper are not ruled out, since a variety of real settlements have been built to face different situations (given by nature, socioeconomic, political conditions, among others). However, our ansatz is that they are normally second or higher order due to the fact that only a disaster or aggressive urban policies (of any nature) would lead to a discontinuity in the entropy (as in first-order phase transitions) as a function of the threshold at a particular time in the evolution of the urban areas. References [1] Rozenfeld HD, Rybski D, Andrade Jr JS, Batty M, Stanley HE, Makse HA. Proc Nat Acad Sci USA 2008;105:18702. [2] Witten TA, Sander LM. Phys Rev Lett 1981;47:1400. [3] Makse HA, Havlin S, Stanley HE. Nature 1995;377:608–12. [4] Szolnoki A, Perc M, Szabó G. Phys Rev Lett 2012;109:078701. [5] Szolnoki A, Szabó G, Perc M. Phys Rev E 2011;83:036101. [6] Song Ch, Koren T, Wang P, Barabási AL. Nature Phys 2010;6:818–23. [7] Mei CQ, Liu YJ. Commun Theor Phys 2011;56:945–51. [8] Wilson A, Dearden J. J Geo Syst 2011;13:1–16. [9] Song Ch, Wang P, Makse HA. Nature 2008;453:629–32. [10] Allen PM. Cities and regions as self-organizing systems: models of complexity. NY: Gordon and Breach Science Publishers; 1997. [11] Alexander A. Archit Forum 1965;122(1):58–62 (Part I). 1965;122(2):58–62 (Part II). [12] Fragkias M, Seto KC. Environ Plan B: Plan Des 2007;34(5):858–83. [13] Batty M. Cities and complexity: understanding cities with cellular automata, agent-based models, and fractals. MA: The MIT Press; 2007. [14] Murcio R, Rodríguez-Romo S. Physica A 2011;390:2895–903. [15] Korotayev A, Malkov A, Khaltourina D. Introduction to social macro dynamics: secular cycles and millennial trends. Moscow: URSS; 2006.

31

[16] Duncan A. Feed or feedback: agriculture, population dynamics and the state of the planet. Tuross Head: International Books; 2003. [17] Arthur WB. Increasing returns and path dependence in the economy. MI: The University of Michigan Press; 1994. [18] E.W. Weisstein, ‘‘Von Neumann Neighborhood’’ from mathworld—a wolfram web resource. . [19] Fujita M, Thisse JF. Economics of agglomeration. Cities, industrial location, and regional growth. MA: Cambridge University Press; 2002. [20] Murcio R, Rodríguez-Romo S. Physica A 2009;388:2689–98. [21] Pauling L. J Am Chem Soc 1935;57:2680–4. [22] Craig NC. Entropy analysis. New York: VCH; 1992. p. 86–92. [23] Hamming RW. Error detecting and error correcting codes. Bell Syst Tech J 1950:149–51. [24] Shannon CE. A mathematical theory of communication. Bell Syst Tech J 1948;27:379–423. 623–56. [25] Xu Rui, Wuncsh II Donald C. Clustering. IEEE Press. Wiley; 1996. [26] Yao J, Dash M, Tan ST, Liu H. Entropy-based fuzzy clustering and fuzzy modeling. Fuzzy Sets Syst 2000;113:381–8. [27] Klir GJ. Uncertainty and information foundations of generalized information theory, 9. Hoboken, NJ: John Wiley & Sons; 1950. p. 147–60 [2005]. [28] Herlihy M, Shavit N. The art of multiprocessor programming. Morgan Kaufmann; 2008. [29] Chapman B, Jost G, Van der Pas R. Using openmp: portable shared memory parallel programming. The MIT press; 2007. [30] CUDA C Programming Guide, Nvidia Corporation, 2010. [31] Hill TL. Statistical thermodynamics. Reading: Addison-Wesley; 1960. p. 280. [32] Tsang IR, Tsang IJ. Phys Rev E 1999;60:2684. [33] Stauffer D, Aharony A. Introduction to percolation theory. 2nd ed. London: Taylor and Francis; 1994. [34] Hoshen J, Kopelman R. Phys Rev B 1976;14:3438. [35] Tsang IJ, Tsang IR, Van Dyck D. Phys Rev E 2000;62:6004. [36] Radillo-Díaz A, Pérez LA, del Castillo-Mussot M. Phys Rev E 2009;80:066107. [37] Tsang IR, Tsang IJ. J Phys A 1997;30:L239. [38] De Boor C. A practical guide to splines. NY: Springer; 2001 [Revised edition].