Non-Maxwellian Kinetic Modelling of City Size Distribution

Lijuan Yu, Jie Liao^*

Department of Mathematics, East China University of Science and Technology, Shanghai, China.

DOI: 10.4236/jamp.2021.96090                

Abstract

To extend the kinetic formulation of city size distribution introduced in [1], the non-Maxwellian kinetic modeling is introduced in the present study, in which a variable collision kernel is used in the underlying kinetic equation of Boltzmann type. By resorting to the well-known grazing asymptotic, a kinetic Fokker-Planck counterpart is obtained. The equilibrium of the Fokker-Planck equation belongs to the class of generalized Gamma distributions. Numerical test shows good fit of the generalized Gamma distribution with the city size distribution of China.

Keywords

Variable Collision Kernel, Kinetic Fokker-Planck Equation, Generalized Gamma Distribution

1. Introduction

The city size distribution has been long recognised to satisfy a very simple distribution law since Zipf, which is attributed to the generic least effort principle of human behavior [2]. Denote the number of cities having a population size between v and $v+\text{d}v$ by $h\left(v\right)\text{d}v$, and the associated cumulative probability distribution function by

$R\left(v\right)={\displaystyle {\int }_v^{\infty }}h\left(w\right)\text{d}w\mathrm{.}$

Zipf found that empirically, $R\left(v\right)\sim 1/v^{\gamma }$, with $\gamma \approx 1$, especially when focusing on large cities, that is, when v is big. As well-known in city size distribution literatures, this empirical law is quite close to reality for most societies across time [3]. However, existing empirical evidence suggests that Zipf’s law is not always observable even for the upper-tail cities of a territory [3] - [9]. The controversy with empirical findings arises, may due to sample selection biases, methodological weaknesses or data limitations. The hypothesis of Zipf’s law is more likely to be rejected for the entire city size distribution hence alternative distributions may be suggested in such cases.

Moreover, although very accurate in recovering the city size distribution, these contributions do not, in general, give the mechanism for the formation in time of this behavior, thus a theoretical derivation of Zipf’s law for cities has been the object of many studies, see for example [10] [11] [12] [13] and references therein, and more recently, by using kinetic modeling [1] [14], in which Boltzmann type and Fokker-Planck type equations for the size distribution of cities are obtained, by introducing interactions based on some migration rule among cities. This model demonstrated that the city size distribution is kinetically related to some factors such as the rate or the tendency of migration of the inhabitants. As noticed in [1] [14], the reasons behind migration are very complex, and it is quite difficult to select one or another reason as dominant. The different choices in the parameters of the kinetic interactions may explain the origin of different effects or even a mixture of effects, which give in the limit a distribution that can be closer to a Pareto or Zipf law, or a lognormal density or others. In any case, the kinetic modeling considered in this framework is useful to clarify the formation of various distributions in terms of various different microscopic interactions.

We further mention that kinetic models were originally used to describe the dynamics of rarefied gas by constructing a Boltzmann-type equation to analyse the effects of the discrete structure of gas molecules [15]. In recent years, various kinetic models have been developed to study the social and economic interactions in multi-agent systems, for example, in social sciences, the statistical description of wealth distribution [16] [17] [18], opinion formation [19] [20] [21], knowledge formation [22], belief formation [23], criminality [24] and so on. Note that a multi-agent system is often composed of “agents” rather than particles; a kinetic model is used to describe the collective behaviour of individuals in a multi-agent system.

At the kinetic level in [1] [14], the Boltzmann collision operator has been selected to be of Maxwellian type as in the classical kinetic theory, that is, the collision kernel is chosen as a constant that does not depend on the “relative velocity of the molecules”. In the context of city size distribution, the Maxwellian hypothesis corresponds to the strong assumption that the migration rate between agents (cities) does not depend on the amount of inhabitants, thus a constant collision kernel is used. This is a simplification of the sophisticated problem, such that, it could be more easily handled from the mathematical point of view. To extend the kinetic formulation in [1] [14], in the present study, we introduce in the underlying kinetic equation of Boltzmann type a variable collision kernel as in [25].

The arrangement of the rest of the paper is as follows. We will show the detailed kinetic modeling of the problem in Section 2, and derive the quasi-invariant limit in Section 3. Finally, we will carry out some numerical tests to validate the model in Section 4. Note that for the quasi-invariant limit, we will show the asymptotic procedure leading from the kinetic description of Boltzmann type to the Fokker-Planck equation. The equilibrium of the Fokker-Planck equation belongs to the class of generalized Gamma distributions. The model test is based on a collection of 332 cities (prefecture-level administrative regions) in China in the 2019 statistical Yearbook, which fits well with the generalized Gamma distribution.

2. Kinetic Modeling of City Size Distribution

To study the evolution of the city size distribution by kinetic models [1], one first needs to specify the microscopic “collision” rules to describe the change of the population of a city. Consider a multi-agent system in which all agents (cities) are assumed to be indistinguishable [26]. A city’s state at any instant of time $t\ge 0$ is completely characterized by its number of inhabitants v. To avoid inessential difficulties, we can simply assume that $v\in {ℝ}^{+}$ although it is clear that v is a natural number. Consequently, the distribution of the multi-agent system, the city size distribution, can be fully characterized by an unknown probability density function $f=f\left(v\mathrm{,}t\right)$.

Follow [1], we assume that the number of residents of a city will essentially increase with the inflow of immigrants and decrease with the outflow of emigrants. At the same time, due to some uncontrollable factors, the population will change for some other uncertain reasons and show random fluctuations. Hence, the microscopic variation of the city size v is the result of three different contributes

$v^{\mathrm{<mo>\; *\; </mo>}}=v-E\left(v\right)v+I_E\left(v\right)z+\eta v\mathrm{,}$ (2.1)

where

• $v\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}$ : the number of inhabitants of a city before and after a microscopic interaction process, respectively;

• $z\in {ℝ}_{+}$ : the amount of population which can migrate towards a city from the environment (the multi-agent system). This value is usually sampled by a certain given distribution function $\mathcal{E}\left(z\right)$, which characterizes the environment itself;

• $\eta$ : a random variable with zero mean and bounded variance, that is, $\langle \eta \rangle =0$, $\langle {\eta }^2\rangle =\sigma$ with $\sigma >0$ suitably small (such that $v^{\mathrm{<mo>\; *\; </mo>}}$ to be positive);

• $E\left(v\right)\mathrm{,}{I}_E\left(v\right)$ : the rate of variation of the city size v consequent to internal and external mechanism, respectively. More precise description of them will be prescribed in below.

Internal mechanism. For $E\left(v\right)$ related to the internal mechanism, we use the concept of “value function” originally used in the study of the distribution of wealth by Kahneman and Twersky [27]: losses weigh heavier than gains in the change of the value function, that is, the value function is concave in the domain of gains and convex in the domain of losses, thus considerably steeper for losses than for gains. In the collision (2.1), the function $E(\cdot )$ plays the role of the value function, which can be taken with the form [28]

$E_{\epsilon }\left(v\right)=\lambda \frac{{\text{e}}^{\epsilon \left(s^{\delta }-1\right)/\delta }-1}{{\text{e}}^{\epsilon \left(s^{\delta }-1\right)/\delta }+1}\mathrm{,}s=\frac{v}{\bar{v}}\mathrm{,}$ (2.2)

in which the value $\bar{v}$ defines an ideal city size, $\epsilon$ is a small positive parameter introduced to represent the strength of the interaction, $0<\lambda <1$ and $0<\delta \le 1$ are used to quantify the intensity of migration rates near the ideal city size $\bar{v}$. For more explanation on the choice of the value function, we refer [28]. Here, in order to simplify the model, we first consider the ideal size of all cities in the whole system as a given value. However, for different countries, this value may depend on the history, political system or cultural background of different countries, or other factors. It is obvious that E is bounded, $-\lambda \le E\left(v\right)\le \lambda$, and E is negative when the city size v is below the ideal size $\bar{v}$, and positive in the opposite situation. Hence, this quantity describes the tendency of the population to reach the ideal size $\bar{v}$ if $v\ne \bar{v}$, with the reason that people prefer to live in a city of population $\bar{v}$.

External mechanism. A non-negative function $I_E\left(v\right)$ can be used to describe a measure of the immigration rate. For simplicity, following [1], a non-negative constant, i.e., $I_E\left(v\right)=I_E>0$, is chosen in the present study. More general choice of $I_E\left(v\right)$ can be [1]

$I_E\left(v\right)=\mu \frac{v^{\alpha }}{1+v^{\alpha }}$

for some positive parameters $\mu$ and $0<\alpha \le 1$ characterizing the intensity of immigration rate.

With the interaction rule (2.1), the variation in time of $f\left(v\mathrm{,}t\right)$ satisfies a linear Boltzmann-like equation [26], which can be written in weak form, for all smooth functions $\phi \left(v\right)$ (the observable quantities),

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{{ℝ}_{+}}}\phi \left(v\right)f\left(v\mathrm{,}t\right)\text{d}v=\langle {\displaystyle {\int }_{{ℝ}_{+}^2}}\chi \left(v\right)\left(\phi \left(v^{\mathrm{<mo>\; *\; </mo>}}\right)-\phi \left(v\right)\right)f\left(v\mathrm{,}t\right)\mathcal{E}\left(z\right)\text{d}v\text{d}z\rangle \mathrm{,}$ (2.3)

In (2.3), the notation $\langle \cdot \rangle$ denotes mathematical expectation taking into account the presence of the random variable $\eta$ in (2.1). And the function $\chi \left(v\right)$ denotes the collision kernel, which assigns to the interaction a certain probability to occur. Note that in [1], the simplification of the Maxwell molecules, leading to a constant interaction kernel $\chi$, has been assumed. To extend, we notice that the distribution of city size in a country has a close relationship with the national conditions of the country. It is also closely related to many factors such as the speed of urban economic development, urban construction and development conditions.

According to the National Bureau of Statistics (NBS) of China, the size of population of China’s cities has been constantly expanding over the past 70 years, with large, small and medium-sized cities distributed across the country. Among them, small cities attract people because of the low threshold of Hukou. Due to people’s pursuit for better quality of life, many people are willing to live in big cities such as Beijing or Shanghai which can provide people with more employment opportunities, better economic income and higher education, etc. However, big cities are already too crowded, in recent years, with China’s urbanisation, medium-sized cities are also attractive for their new opportunities and living conditions. The migration between different cities is much often than ever. There is strong evidence that the population mobility is greater in cities with a large population, such as Beijing, Shanghai, compared with smaller cities such as Yibin, Xining etc. On the other side, due to geographic and historic reasons, the number of cities (prefecture-level administrative regions) is relatively fixed, so the interactions for v small or near zero should be excluded. Thus, the collision frequency may proportional to the city size v. Hence, to elaborate this behavior, it seems natural to consider a variable collision kernel that

$\chi \left(v\right)=\kappa \cdot v^{\alpha }\mathrm{,}$ (2.4)

where the constants $\kappa >0$ and $0<\alpha <1$. This kernel assigns a high probability of interactions for cities with large population, and low probability of interactions for cities with population v close to zero. By taking into account this new assumption, we consider in the following that $f\left(v\mathrm{,}t\right)$ satisfies the linear kinetic model

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{{ℝ}_{+}}}\phi \left(v\right)f\left(v\mathrm{,}t\right)\text{d}v=\kappa \langle {\displaystyle {\int }_{{ℝ}_{+}^2}}{v}^{\alpha }\left(\phi \left(v^{\mathrm{<mo>\; *\; </mo>}}\right)-\phi \left(v\right)\right)f\left(v\mathrm{,}t\right)\mathcal{E}\left(z\right)\text{d}v\text{d}z\rangle \mathrm{.}$ (2.5)

3. Quasi-Invariant Limit: The Fokker-Planck Equation

In order to describe the development of city size distribution more accurately and intuitively, we carry out the quasi-invariant limit. In this Section, we illustrate the main steps leading from Equation (2.5) to its Fokker-Planck limit. To avoid inessential difficulties, we will assume that the environmental distribution $\mathcal{E}$ has a certain number of bounded moments, more precisely.

$M_{\beta }\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_{{ℝ}_{+}}}{z}^{\beta }\mathcal{E}\left(z\right)\text{d}z\le \infty \mathrm{,}0\le \beta \le \mathrm{4,}$ (3.1)

meanwhile, we introduce the notation

$m_a\left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_{{ℝ}_{+}}}{v}^{a}f\left(v\mathrm{,}t\right)\text{d}v\mathrm{,}a>0.$ (3.2)

It’s obvious that the kinetic equation is mass preserving by taking $\phi \left(v\right)=1$ in (2.5).

For the quasi-invariant limit, one assumes that a single interaction determines only an extremely small change of the value v. Therefore, a small parameter $\epsilon$ is introduced and we consider the scaling

$\begin{array}{l}{I}_E\to \epsilon I_E\mathrm{,}\eta \to \sqrt{\epsilon }\eta \mathrm{.}\hfill \end{array}$ (3.3)

At this point, under the effect of $\epsilon$, the interaction will only produce a very small change to the population size of a city. Obviously, the conservation of “mass” of the system still holds under the scaling. To observe the evolution of the mean value, in (2.5), take $\phi \left(v\right)=v$, there is

$\frac{\text{d}{m}_1\left(t\right)}{\text{d}t}=\kappa \left[-{\displaystyle {\int }_{{ℝ}_{+}}}{v}^{\alpha +1}{E}_{\epsilon }\left(v\right)f\left(v\mathrm{,}t\right)\text{d}v+\epsilon I_E{M}_1{m}_{\alpha }\left(t\right)\right]\mathrm{.}$ (3.4)

Next, denote $A_{\epsilon }\left(v\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{\epsilon }{E}_{\epsilon }\left(v\right)$, then

$A_{\epsilon }\left(v\right)=\frac{\lambda }{\epsilon }\frac{{\text{e}}^{\epsilon \left({\left(v/\bar{v}\right)}^{\delta }-1\right)/\delta }-1}{{\text{e}}^{\epsilon \left({\left(v/\bar{v}\right)}^{\delta }-1\right)/\delta }+1}\to \frac{\lambda }{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{as}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\epsilon \to 0.$ (3.5)

Now, we can resort to a scaling of time to observe an evolution of the average value independent of $\epsilon$. Setting $\tau =\epsilon t$, $f_{\epsilon }\left(v,\tau \right)=f\left(v,t\right)$, then the evolution of the average value for $f_{\epsilon }\left(v\mathrm{,}\tau \right)$ satisfies

$\begin{array}{c}\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{{ℝ}_{+}}}v{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v=\kappa {\displaystyle {\int }_{{ℝ}_{+}}}\left(I_E{M}_1-\frac{\lambda v}{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)\right)v^{\alpha }{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v\\ +\kappa {\displaystyle {\int }_{{ℝ}_{+}}}\left(\frac{\lambda }{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)-\frac{1}{\epsilon }{E}_{\epsilon }\left(v\right)\right)v^{\alpha +1}{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v\mathrm{.}\end{array}$ (3.6)

Since (3.5) means that the second term vanishes as $\epsilon \to 0$, one obtains in the limit a closed form for the evolution of the mean value.

$\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{{ℝ}_{+}}}v{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v=\kappa {\displaystyle {\int }_{{ℝ}_{+}}}\left(I_E{M}_1-\frac{\lambda v}{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)\right)v^{\alpha }{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v\mathrm{.}$ (3.7)

It can be observed that the evolution of the mean ${\displaystyle {\int }_{{ℝ}_{+}}}v{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v$ does not depend on $\epsilon$. Since for $\epsilon \ll 1$ the microscopic interactions produce a very small change of the value v, a finite variation of the mean density can be observed only if agents in the system undergo a huge number of interactions in a fixed period of time to restore the original evolution. Similarly, with this scaling, one obtains in the limit a closed form for the evolution of the second moment

$\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{{ℝ}_{+}}}{v}^2{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v=\kappa {\displaystyle {\int }_{{ℝ}_{+}}}\left(\sigma v^2+2I_E{M}_{1}v-\frac{\lambda v^2}{\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)\right)v^{\alpha }{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v$ (3.8)

The above analysis can be used to justify the passage from the kinetic model (2.5) to its continuous counterpart given by a Fokker-Planck type equation. Given a smooth function $\phi \left(v\right)$, let us expand in Taylor series $\phi \left(v^{\mathrm{<mo>\; *\; </mo>}}\right)$ around $\phi \left(v\right)$. First, by the scaling (3.3), it holds

$\langle v^{\mathrm{<mo>\; *\; </mo>}}-v\rangle =-\epsilon A_{\epsilon }\left(v\right)v+\epsilon I_{E}z\mathrm{,}$ (3.9)

$\langle {\left(v^{\mathrm{<mo>\; *\; </mo>}}-v\right)}^2\rangle =\epsilon \sigma v^2+{\epsilon }^2\left(A_{\epsilon }^2\left(v\right)v^2+I_E^2{z}^2-2I_E{A}_{\epsilon }\left(v\right)vz\right)\mathrm{.}$ (3.10)

Therefore, in terms of powers of $\epsilon$, we easily obtain the expression

$\langle \phi \left(v^{\mathrm{<mo>\; *\; </mo>}}\right)-\phi \left(v\right)\rangle =\epsilon \left[{\phi }^{\prime }\left(v\right)\left(I_{E}z-A_{\epsilon }\left(v\right)v\right)+\frac{1}{2}\sigma v^2{\phi }^{″}\left(v\right)\right]+R_{\epsilon }\left(v\mathrm{,}z\right)\mathrm{,}$ (3.11)

where the remainder term $R_{\epsilon }\left(v\mathrm{,}z\right)$ vanishes at the order ${\epsilon }^{\frac{3}{2}}$ as $\epsilon \to 0$. Therefore, as $\epsilon \to 0$, we can obtain that in consequence of the scaling (3.3) the weak form of the kinetic model (2.5) is well approximated by the weak form of a linear Fokker-Planck equation

$\begin{array}{l}\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{{ℝ}_{+}}}\phi \left(v\right)f_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v\\ =\kappa {\displaystyle {\int }_{{ℝ}_{+}}}\left[{\phi }^{\prime }\left(v\right)\left(I_E{M}_1-\frac{\lambda v}{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)\right)+\frac{1}{2}\sigma v^2{\phi }^{″}\left(v\right)\right]v^{\alpha }{f}_{\epsilon }\left(v\mathrm{,}\tau \right)\text{d}v\mathrm{.}\end{array}$ (3.12)

Providing the boundary terms produced by the integration by parts vanish, Equation (3.12) coincides with the weak form of the Fokker-Planck equation

$\frac{\partial g\left(v\mathrm{,}\tau \right)}{\partial \tau }=\frac{\kappa \sigma }{2}\frac{{\partial }^2}{\partial v^2}\left(v^{2+\alpha }g\left(v\mathrm{,}\tau \right)\right)+\kappa \frac{\partial }{\partial v}\left[\left(\frac{\lambda v}{2\delta }\left({\left(\frac{v}{\bar{v}}\right)}^{\delta }-1\right)-I_E{M}_1\right)v^{\alpha }g\left(v\mathrm{,}\tau \right)\right]\mathrm{.}$ (3.13)

Without loss of generality, we will simplify Equation (3.13) by assuming

$\theta :=\frac{\kappa \sigma }{2},\mu :=\frac{\kappa \lambda }{2\delta {\bar{v}}^{\delta }},\xi :=\frac{\kappa \lambda }{2\delta },\varsigma :=\kappa I_E{M}_1.$ (3.14)

Thus, the resulting Fokker-Planck equation takes the form

$\frac{\partial g\left(v,\tau \right)}{\partial \tau }=\theta \frac{{\partial }^2}{\partial v^2}\left(v^{2+\alpha }g\left(v,\tau \right)\right)+\frac{\partial }{\partial v}\left(\left(\mu v^{1+\delta }-\xi v-\varsigma \right)v^{\alpha }g\left(v,\tau \right)\right).$ (3.15)

As exhaustively discussed in Ref. [29] [30], the right boundary conditions that guarantee mass conservation are the so-called no-flux boundary conditions given by

$\theta \frac{\text{d}}{\text{d}v}\left(v^{2+\alpha }g\left(v,\tau \right)\right)+{\left(\left(\mu v^{1+\delta }-\xi v-\varsigma \right)v^{\alpha }g\left(v,\tau \right)\right)|}_{v=0,+\infty }=0,t>0.$ (3.16)

With these no-flux boundary conditions, we can obtain the explicit stationary solution of the Fokker-Planck Equation (3.13) by solving the ordinary differential equation of first order

$\theta \frac{\text{d}}{\text{d}v}\left(v^{2+\alpha }g\left(v\right)\right)+\left(\mu v^{1+\delta }-\xi v-\varsigma \right)v^{\alpha }g\left(v\right)=0.$ (3.17)

Using $h\left(v\right)=v^{2+\alpha }g\left(v\right)$ in (3.17) as unknown function, separation of variables gives as unique solution to (3.17) the function

$\begin{array}{c}{g}_{\infty }\left(v\right)=C\cdot v^{\frac{\xi }{\theta }-2-\alpha }\cdot \exp \left(-\frac{\varsigma }{\theta v}-\frac{\mu }{\theta \delta }{v}^{\delta }\right)\\ =C\cdot v^{\frac{\lambda }{\delta \sigma }-2-\alpha }\cdot \exp \left(-\frac{2I_E{M}_1}{\sigma v}-\frac{\lambda }{\sigma {\delta }^2}{\left(\frac{v}{\bar{v}}\right)}^{\delta }\right),\end{array}$ (3.18)

where the positive constant C has been chosen to normalize the equilibrium distribution. It is not difficult to discover that $g_{\infty }\left(v\right)$ tends to 0 as $v\to 0$ and $v\to \infty$. In other words, the city population size distribution obtained under the Non-Maxwellian collision does not exist with too little or too much population, which is more consistent with the real situation, and $g_{\infty }\left(v\right)$ is close to the generalized Gamma density as $v\to \infty$.

4. Numerical Tests

In this section, we will use statistical data to verify the validity of the model. Here, we chose data from the Statistical Yearbook 2019 (27 provinces and 4 municipalities directly under the Central Government) released by the National Bureau of Statistics of China with the population of 332 cities (prefecture-level cities) in 2018. The histogram of the city size distribution is shown in Figure 1.

From this probability distribution, we noticed that cities with a population of 1 million to 2 million are the majority. The number of cities with a population of more than 3 million decreases with the increase of the number of people contained, and the rate of decrease also changes from a sharp decline to a slow convergence to zero with the increase of the number of people. To fit the data with our model, we take a set of parameters

$\delta =0.5,\lambda =0.8,I_E=0.6,M_1=20,\sigma =0.5,\bar{v}=400,\alpha =\mathrm{0.2.}$

The equilibrium distribution of both Maxwellian model [1] and our non-Maxwellian model are shown in Figure 2. This result shows that the non-Maxwellian model fits the city size distribution of China better than the Maxwellian model.

Figure 1. Probability distribution histogram of China’s cities size.

Figure 2. Theoretical steady-state distribution and the real data.

5. Conclusion and Perspectives

In this paper, we introduced non-Maxwellian kinetic modeling, in which a variable collision kernel is used in the underlying kinetic equation of Boltzmann type, to explain the evolution of city size in China. By resorting to the well-known quasi-invariant asymptotic, we obtain a kinetic Fokker-Planck counterpart and the steady-state of city size which is defined as the generalized Gamma distribution. Numerical test shows good fit of the generalized Gamma distribution with the city size distribution of China. However, further understanding of the role of each parameter, for example, the ideal city size $\bar{v}$, is still open. It would also be interesting to investigate the trend of the city size distribution under the effect of fast urbanisation of China in recent and next several years.

Acknowledgements

The research is partially supported by the National Science Foundation of China (Grant Nos. 11871335 and 11971008).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Gualandi, S. and Toscani, G. (2019) Size Distribution of Cities: A Kinetic Explanation. Physica A: Statistical Mechanics and Its Applications, 524, 221-234. https://doi.org/10.1016/j.physa.2019.04.260
[2]	Zipf, G.K. (1949) Human Behavior and the Principle of Least Effort. Addison-Wesley Press, Inc., Cambridge, MA.
[3]	Arshad, S., Hu, S. and Ashraf, B.N. (2018) Zipfs Law and City Size Distribution: A Survey of the Literature and Future Research Agenda. Physica A: Statistical Mechanics and Its Applications, 492, 75-92. https://doi.org/10.1016/j.physa.2017.10.005
[4]	Bee, M., Riccaboni, M. and Schiavo, S. (2013) The Size Distribution of US Cities: Not Pareto, Even in the Tail. Economics Letters, 120, 232-237. https://doi.org/10.1016/j.econlet.2013.04.035
[5]	Benguigui, L. and Blumenfeld-Lieberthal, E. (2007) Beyond the Power Law—A New Approach to Analyse City Size Distributions. Computers, Environment and Urban Systems, 31, 648-666. https://doi.org/10.1016/j.compenvurbsys.2006.11.002
[6]	Gangopadhyay, K. and Basu, B. (2009) City Size Distributions for India and China. Physica A: Statistical Mechanics and Its Applications, 388, 2682-2688. https://doi.org/10.1016/j.physa.2009.03.019
[7]	Ioannides, Y. and Skouras, S. (2013) US City Size Distribution: Robustly Pareto, but Only in the Tail. Journal of Urban Economics, 73, 18-29. https://doi.org/10.1016/j.jue.2012.06.005
[8]	Malevergne, Y., Pisarenko, V. and Sornette, D. (2011) Testing the Pareto against the Lognormal Distributions with the Uniformly Most Powerful Unbiased Test Applied to the Distribution of Cities. Physical Review E, 83, Article ID: 036111. https://doi.org/10.1103/PhysRevE.83.036111
[9]	Rozenfeld, H., Rybski, D., Gabaix, X. and Makse, H. (2011) The Area and Population of Cities: New Insights from a Different Perspective on Cities. American Economic Review, 101, 2205-2225. https://doi.org/10.1257/aer.101.5.2205
[10]	Gabaix, X. (1999) Zipf’s Law for Cities: An Explanation. The Quarterly Journal of Economics, 114, 739-767. https://doi.org/10.1162/003355399556133
[11]	Ghosh, A., Chatterjee, A., Chakrabarti, A.S. and Chakrabarti, B.K. (2014) Zipfs Law in City Size from a Resource Utilization Model. Physical Review E, 90, Article ID: 042815. https://doi.org/10.1103/PhysRevE.90.042815
[12]	Marsili, M. and Zhang, Y.-C. (1998) Interacting Individuals Leading to Zipf’s Law. Physical Review Letters, 80, 2741-2744. https://doi.org/10.1103/PhysRevLett.80.2741
[13]	Zanette, D.H. and Manrubia, S.C. (1997) Role of Intermittency in Urban Development: A Model of Large-Scale City Formation. Physical Review Letters, 79, 523-526. https://doi.org/10.1103/PhysRevLett.79.523
[14]	Gualandi, S. and Toscani, G. (2019) Human Behavior and Lognormal Distribution: A Kinetic Description. Mathematical Models and Methods in Applied Sciences, 29, 717-753. https://doi.org/10.1142/S0218202519400049
[15]	Cercignani, C., Illner, R. and Pulvirenti, M. (1994) The Mathematical Theory of Dilute Gases. Vol. 106, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4419-8524-8
[16]	Chatterjee, A., Chakrabarti, B.K. and Manna, S.S. (2004) Pareto Law in a Kinetic Model of Market with Random Saving Propensity. Physica A: Statistical Mechanics and its Applications, 335, 155-163. https://doi.org/10.1016/j.physa.2003.11.014
[17]	Cordier, S., Pareschi, L. and Toscani, G. (2005) On a Kinetic Model for a Simple Market Economy. Journal of Statistical Physics, 120, 253-277. https://doi.org/10.1007/s10955-005-5456-0
[18]	Düring, B., Matthes, D. and Toscani, G. (2008) Kinetic Equations Modelling Wealth Redistribution: A Comparison of Approaches. Physical Review E, 78, 056-103. https://doi.org/10.1103/PhysRevE.78.056103
[19]	Ben-Naim, E. (2005) Opinion Dynamics: Rise and Fall of Political Parties. Europhysics Letters, 69, 671-677. https://doi.org/10.1209/epl/i2004-10421-1
[20]	Boudin, L., Mercier, A. and Salvarani, F. (2012) Conciliatory and Contradictory Dynamics in Opinion Formation. Physica A: Statistical Mechanics and Its Applications, 391, 5672-5684. https://doi.org/10.1016/j.physa.2012.05.070
[21]	Toscani, G. (2006) Kinetic Models of Opinion Formation. Communications in Mathematical Sciences, 4, 481-496. https://doi.org/10.4310/CMS.2006.v4.n3.a1
[22]	Pareschi, L. and Toscani, G. (2014) Wealth Distribution and Collective Knowledge. A Boltzmann Approach. Philosophical Transactions of The Royal Society A, 372, Article No. 0396. https://doi.org/10.1098/rsta.2013.0396
[23]	Brugna, C. and Toscani, G. (2015) Kinetic Models of Opinion Formation in the Presence of Personal Conviction. Physical Review E, 92, Article ID: 052818. https://doi.org/10.1103/PhysRevE.92.052818
[24]	Bellomo, N., Colasuonno, F., Knopoff, D. and Soler, J. (2015) From a Systems Theory of Sociology to Modeling the Onset and Evolution of Criminality. Networks & Heterogeneous Media, 10, 421-441. https://doi.org/10.3934/nhm.2015.10.421
[25]	Furioli, G., Pulvirenti, A., Terraneo, E. and Toscani, G. (2019) Non-Maxwellian Kinetic Equations Modeling the Evolution of Wealth Distribution. Mathematical Models and Methods in Applied Sciences, 30, 685-725. https://doi.org/10.1142/S0218202520400023
[26]	Pareschi, L. and Toscani, G. (2013) Interacting Multiagent Systems. Kinetic Equations and Monte Carlo Methods. Oxford University Press, Oxford.
[27]	Kahneman, D. and Tversky, A. (Eds.) (2000) Choices, Values, and Frames. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511803475
[28]	Dimarco, G. and Toscani, G. (2019) Kinetic Modeling of Alcohol Consumption. Journal of Statistical Physics, 177, 1022-1042. https://doi.org/10.1007/s10955-019-02406-0
[29]	Toscani, G. (2019) Statistical Description of Human Addiction Phenomena. (Preprint)
[30]	Furioli, G., Pulvirenti, A., Terraneo, E. and Toscani, G. (2017) Fokker-Planck Equations in the Modelling of Socio-Economic Phenomena. Mathematical Models and Methods in Applied Sciences, 27, 115-158. https://doi.org/10.1142/S0218202517400048