Prepared by Avinav Sharan, Himanshu Sharma and Gaurav Kshirsagar

**Rank-Size Population Distributions And Dynamics of Cities
**

We analyse how the shift in population of a given city happens with time given their ranks. We try to rank them by their population size.\\\\

We define some parameters to represent the population and their ranks. \\

$P_{1}(t),…..,P_{N}(t)$: populations of N cities at time t

$r_{1},….,r_{N}$: population ranks of cities in decreasing

$P_{i,r_{i}=1}≥P_{j,r_{j}=2}≥P_{m,r_{m}=N}$: population sorted in decreasing order

$\Gamma(t+1)=\sum_{i}|p_{i,r_{i}(t+1)}(t+1) – p_{i,r_{i}(t)}(t)|$: percentage shift in population, where $pi,r=P_{i,r}/\sum_{i}P_{i,r}$

$\wedge(t+1)=\sum_{i}|r_{i}(t+1)-r_{i}(t)|/N$ percentage shift in rank order

insert image 1

insert image 2

We plotted the $\log(rank)$ vs $\log(population)$ for 458 standardized administrative areas from 1901 to 1991. We can see that it shows a power law distribution between rank and population. The curve is very stable as we can see from the plot. As 1901 and 1901 the curve is similar with same coefficient even with population of the cities changes.\\

How can we explain this above pattern based on microscopic agent behaviour and interaction. If we can explain then we can define the model for population development:\\

insert image 3

We will try to explain one of the multiple interpretation.

let increment in population is some random variable with Gaussian distribution.\\

If we have the initial state (population at time t = 0, $P_{i_{0}}$) and increment of population at time t = 0 is $\epsilon_{0}$ and for city 1 increment is $\epsilon_{1}$ and so on.

$\log P_{i}(t) = \log P_{i}(0) + \sum_{\tau=0}^t \epsilon_{i}(\tau)$

Also assume mean is $\mu$ and variance is $\sigma^2$, than we get the log normal distribution. (meaning logarithm of the population is normal distribution)

Now simulate the same number of cities, 458 cities and we start from a uniform population distribution. In each city we simulate the population growth according to Gaussian distribution. We sample the scaling factor of all of the 458 cities and we simulate from 991 AD to 1991 and what we find out is indeed that we have some kind of nice straight line of the log population shear versus log ranks. It is one feasible interpretation, hypothesis about city development.

**Explaining city growth with Ants Model**

How can we explain the above model (growth in population of cities)?

insert image 4

We use this agents based model to explain the growth of cities.

In the picture we have nest in the center, and three food sources. In starting ants move around randomly until they reach the food source. They carry the food and leave the chemical trails while returning for the other ants. So that they can find the food source easily. We can see the simulation of this model in the NetLogo tool.\\

We can identify many things with this model. First is that nearest food resource will be probabilistically explored first and furthest will be explored last.

Second, whenever ant find it and they will leave trails. And It will have large concentration of the chemicals. So we can relate this ant model with the behaviour of the people.

**Formulation of Ants Model**

Following random variables represent the Ants model:

Position of any ant at time t : $A(x_i(t),y_i(t))$

Home of ant i: $A(X_i,Y_i)$

Position of resource k : $R(X_k, Y_k)$

The chemical concentration of scent at point (x,y) : $\psi(x,y)$

There are few things to note in the above representation. ‘Home of ant’ does not vary with time. The chemical concentration of scent is a function of time.

**Process**

When the ant is moving back to home with food, it leaves the trail of chemical scent. This is represented by the following equation:

\begin{equation}

\psi(x_i(t+1),y_i(t+1)) = \psi(x_i(t),y_i(t)) + \phi

\end{equation}

This equation says that concentration of chemical scent depends on the same at the postion in previous time step, with some added noise $\phi$.

Generally, we can represent the ants motion as the following:

\begin{equation}

x_i(t+1) = x_i(t) + \Delta x_i(t) + \epsilon_x(t)

\end{equation}

\begin{equation}

y_i(t+1) = y_i(t) + \Delta y_i(t) + \epsilon_y(t)

\end{equation}

That is, the position of ant depends on the gradient of previous time step plus some added noise, $\epsilon$.

The gradient, $\Delta x$ or $\Delta y$ depends on the two cases:

When the ant is seeking food: We know that the ant is trying at random (given by $\epsilon(t)$, but following any chemical trail if it finds. This is intuitively represented by the equation:

\begin{equation}

(\Delta x_i(t), \Delta y_i(t)) = c \cdot \nabla \psi(x_i(t),y_i(t)) + \epsilon(t)

\end{equation}

When the ant is moving food back: We know that ant follows the path that is dependent on the ant’s home. This is given by the equation:

\begin{equation}

(\Delta x_i(t), \Delta y_i(t)) = (f(X_i),f(Y_i))

\end{equation}

This formulation let us think the ant model probabilistically.

**Simulation and Observation of Ants Model**

First we simulate a \textbf{simplified ants model, with one home (at origin) and one food center} (represented by a square box). In the image below we can observe three time stages. Intitially (at t=0), we observe that the ants move in random directions. There are no definite permanent tracks yet.\par

As we skip to time, t=200, we observe that ants have found the food and the track has started to emerge. But the ants are still randomly distributed in the space.\par

Finally, at time, t=2000, we observe that the track is well defined and the ants have started to cluster around this track.\linebreak

Insert image 5

Next we simulate an \textbf{ant model with random homes and multiple resources}. In the below image, we can clearly observe one of the main properties of cities, that is roads. We can observe star like structures (roads) between homes and resources after sufficient number of iterations.

Insert image 6

Next let us make this formulation more interesting and complex. We add following features instead of resources: ‘Industry’ ($I_{xy}$) and ‘Service’ ($S_{xy}$). We also dynamically add population ($P_{xy}$) to the model. ‘Service’ is dependent on ‘Industry’. The population is dependent of ‘Service’ and ‘Industry’. We can model any kind of interaction between Industry, Service and Population.

The process of the model is defined by:

The Position at any time, depends on previous time step population, the diffusion of population ($\nabla^2P_{xy}(t)$) (laplacian) at previous time step, the interaction of industry and service ($o_{xy}^{(k)}(t)$) and random noise ($\epsilon$).

\begin{equation}

P_{xy}(t+1) = P_{xy}(t) + \theta^{(P)}\nabla^2P_{xy}(t) + \Sigma_k o_{xy}^{(k)}(t) + \epsilon_{xy}^{(P)}(t)

\end{equation}

where, $\theta$ is diffusion constant

Similarly, we can represent the denisty of Indurtrial resouce and Service:

\begin{equation}

I_{xy}(t+1) = I_{xy}(t) + \theta^{(I)}\nabla^2 I_{xy}(t) + \Sigma_k o_{xy}^{(k)}(t) + \epsilon_{xy}^{(I)}(t)

\end{equation}

\begin{equation}

S_{xy}(t+1) = S_{xy}(t) + \theta^{(S)}\nabla^2 S_{xy}(t) + \Sigma_k o_{xy}^{(k)}(t) + \epsilon_{xy}^{(S)}(t)

\end{equation}

**Observations**

In the image below we can observe the simulation. Intially, we see the formulation of nuclear cities (at t=100). We observe cluster of population which represent cities. The log-size of the ‘cities’ is linear with log-rank, which we had seen in the british cities data from 1901 to 1991. We also observe the formation of roads.

Insert image 7

**SandPile model**

Let us perform the following experiment:

Drop sand uniformally onto a table. Gradually, sandpile height increases, and the sand starts to fall in surrounding areas. The key observation is that sand fall can trigger more sand fall and an ‘avalanche’ will occur.

In Sandpile Model we try to model the aboce experiment. Model Representation:

\begin{equation}

X_{i,j} \geq 4: X_{i,j} = X_{i,j} – 4,

\end{equation}

\begin{equation}

X_{i \pm 1,j} = X_{i \pm 1,j} + 1,

\end{equation}

\begin{equation}

X_{i,j \pm 1} = X_{i,j \pm 1} + 1,

\end{equation}

where $X_{i,j}$ is the height of the sand pile.

NetLogo Demo:

Checkout:

http://netlogoweb.org/launch#http://netlogoweb.org/assets/modelslib/Sample%20Models/Chemistry%20&%20Physics/Sandpile.nlogo

You can play with it, and observe that small avalanch occur more often than large ones. We study Sandpile model as it is related to many real-life scenarios. For instance, the migration of the population from one city to another can be said to follow sand pile model. Similarly, we can observe it in forest fires, frequency and size of earthquakes, etc.

**Size of Avalanche based on SandPile model**

Here we are using the SandPile model to understand the working of avalanches.

Attempts to model avalanche behaviour date from the early 20th century, notably the work of Professor Lagotala in preparation for the 1924 Winter Olympics in Chamonix. His method was developed by A. Voellmy and popularised following the publication in 1955 of his Ueber die Zerstoerungskraft von Lawinen (On the Destructive Force of Avalanches). Voellmy used a simple empirical formula, treating an avalanche as a sliding block of snow moving with a drag force that was proportional to the square of the speed of its flow:

\begin{equation}

Force = 1/2 * \rho \nu^2

\end{equation}

Avalanche increases in the following manner having a relation to the previous state

\linebreak

\begin{equation}

A + \phi \rightarrow B

\end{equation}

\begin{equation}

B + \phi \rightarrow C

\end{equation}

\begin{equation}

C + \phi \rightarrow D

\end{equation}

\begin{equation}

D + \phi \rightarrow p_1 D + \phi, p_2 C + 2\phi, B + 3\phi, A + 4\phi

\end{equation}

The kinetic equations for this growth are:

\begin{equation}

\dot{n_A} = n_\phi (p_4 n_D – n_A)

\end{equation}

\begin{equation}

\dot{n_B} = n_\phi (p_3 n_D + n_A – n_B)

\end{equation}

\begin{equation}

\dot{n_C} = n_\phi (p_2 n_D + n_B – n_C)

\end{equation}

\begin{equation}

\dot{n_D} = n_\phi (p_1 n_D + n_C – n_D)

\end{equation}

In general,

\begin{equation}

\dot{n_\phi} = n_\phi (\bar{p} n_D – 1) + \bar{p}\nu\nabla^2 (n_\phi n_D) + \eta(r,t)

\end{equation}

where

\bar{p} = p_1 + 2p_2 + 3p_3 + 4p_4

\nu = diffusion coefficient

\eta = noise

Insert image 8

Avalanches in 2D sandpile model.The dark regions show the size

of avalanches triggered by a single addition.

Though at first glance, it seems that the number of topplings and the area both measure the number of affected sites, a single site may topple more than once in a single avalanche, and hence the two are truly different variables. The duration of the avalanche is defined as the number of updates one must perform on the system before all the sites become stable after the addition of one grain of sand.

Video An avalanche triggered which follows SandPile model

**Avalanche Equilibrium**

For the avalanche to reach its equilibrium we have the probability of the number of grains affected after t steps

Probability generating function:

\begin{equation}

G_t(\lambda_A,\lambda_B,\lambda_C,\lambda_D) = \sum_{N} (N_A,N_B,N_C,N_D) \lambda^{N_A}_A \lambda^{N_B}_B \lambda^{N_B}_B \lambda^{N_D}_D

\end{equation}

At time 1,

\begin{equation}

G_1 = (n_Dp_4)\lambda_A + (n_Dp_3 + n_A)\lambda_B + (n_Dp_2 + n_B)\lambda_C + (n_Dp_1 + n_C)\lambda_D

\end{equation} \linebreak

At time t+1, the recursive relationship with $G_t is

\begin{equation}

G_{t+1} = (n_A\lambda_B + n_B\lambda_C + n_C\lambda_D) (n_Dp_4)\lambda_A + n_D(p_1\lambda_DG_t + p_2\lambda_DG_t^2 + p_3\lambda_BG_t^3 + p_4\lambda_AG_t^4)

\end{equation}

\begin{equation}

G_\infty(\lambda) = \frac{1-\sqrt{a-4\lambda^2n_D(1-n_D)}}{2\lambda n_D}

\end{equation}

Example: Metroploitan Buffalo population growth (1750-1990)

Insert image 9

It has been discusseed and found that when a system diverges the various parameters of the Sandpile model exhibit a power law distribution for scenarios like avalanches or population. For additional information see this http://guava.physics.uiuc.edu/~nigel/courses/563/Essays_2012/PDF/banerjee.pdf

For such a distribution in case of population, we have the following parameters:

\begin{equation}

N(R) = (2R)^D,

\end{equation}

where $\textit{D}$ = fractal dimension of the city 1\leq\textit{D}\leq2

R = effective radius of city

Cumulative populaiton density:

\begin{equation}

P(R) = N(R)/A(R) = (2R)^{D-2}

\end{equation}

where A(R) = Area of square

Incremental Populaion

\begin{equation}

n(R) = \frac{dN(R)}{dR} = D2R^{D-1}

\end{equation}

Incremental density

\begin{equation}

\rho(R) = \frac{dn(R)}{dR} = D^2R^{D-2}

\end{equation}