Xue Xia^1*, Liangli Chen^2
1School of Mathematics, Aba Teachers College, Wenchuan, China.
²School of Mathematics, Southwestern University of Finance and Economics, Chengdu, China.
DOI: 10.4236/jamp.2026.144063   PDF    HTML   XML   52 Downloads   202 Views   Citations

Abstract

This study adopts a kinetic model incorporating a value function to explore the evolution of the city size distribution. Specifically, the value function consistent with prospect theory is utilized to characterize the internal population changes of cities. Based on the interaction rules of population dynamics, we establish the Boltzmann equation and its corresponding Fokker-Planck equation; thereafter, we theoretically analyze population movements by solving these equations. Assuming that the urban population consists of working and non-working people, we employ a utility function to describe cities’ preferences for these two population groups. Numerical analysis is performed, and some key parameters are estimated using real-world data to verify the rationality of the proposed model.

Keywords

City-Size Distribution, Value Function, Boltzmann-Type Equation, Fokker-Planck Equation, Utility Function

Share and Cite:

1. Introduction

There is a large body of research on urban population distribution, employing a variety of research methods. Studies in References [1]-[3] indicate that the urban population distribution follows a power law, while those in References [4] [5] show that it conforms to other regular distributions. Kinetic models have been widely applied in the study of social activities and have achieved certain results [6] [7].

City-size distribution follows power laws, which are shown in a lot of works. Using a stochastic model, Zanette and Manrubia [1] verify that the city-size distribution follows a power law. Benguigui and Blumenfeld-Liebertha [2] set up several models to illustrate that the city-size evolution follows the power law distribution or concave distribution. According to the city data of the southeast and southwest in the United States, Garmestani et al. [3] test the power law behavior of the city-size distribution. Calderín-Ojeda [8] finds that the upper quartile of the commune size data in the early period (1962-1975) obeys the power law distribution. Taking Chinese city-size data, the empirical research on the power law distribution of urban rank-size is presented in [9]. The power law with the index one is Zipf’s law [10], which is a classical explanation of the city-size distribution. Gabaix [11] utilizes a model of random growth and concludes that the city-size distribution obeys Zipf’s law. Marsili and Zhang [12] construct master equations to discuss city distribution, which follows Zipf’s law. Ghosh et al. [13] explain Zipf’s law about city size by using a resource utilization model. Based on the data of the Indian population census and the Chinese population census, Gangopadhyay and Basu [14] prove that the urban size distributions in India and China satisfy Zipf’s law.

Many studies find that urban sizes exhibit different patterns. There are also other works on the distribution of city size. Eeckout [4] [5] concludes that the distribution of urban population obeys the lognormal distribution. Focusing on the role of distance, González-Val [15] analyzes the size distribution of Spanish cities. Li and Zhang [16] find that compared with the power law distribution, the triangle law distribution better fits the Chinese city-size distribution through statistical methods. Levy [17], Zhang et al. [18], Xu and Zhu [19] study the urban distribution from different angles.

Scholars never reach a consensus on the definition of the scope of a city. Brenner and Theodore [20] explore the instability of the geographical patterns of urban space. Jacobs [21] reviews urban geographical studies influenced by relational theory and finds that current research on relationality in urban geography operates with irreconcilable grammars.

Gualandi and Toscani [6] create a kinetic model, which is applied to wealth distribution [22]-[25], opinion formation [26]-[29] and other fields [30]-[32], to study the formation of urban agglomerations when considering the internal mechanism, external mechanism and fluctuation factor. The number of the urban population is represented by the variable $x\left(x\ge 0\right)$ . The interaction equation for the change of urban population in [6] reads

$x^{\text{*}}=x-E\left(x/x_L\right)x+b\left(x\right)z+\eta x,$ (1.1)

where $\eta$ is a random variable. The function $E\left(x/x_L\right)$ describes the population change caused by internal mechanism. When $x<x_L$ , $E\left(x/x_L\right)<0$ . Namely the urban population increases. When $x>x_L$ , $E\left(x/x_L\right)>0$ . Namely the urban population decreases. $b\left(x\right)$ represents the rate of people who move in the city. The variable $z\left(z\ge 0\right)$ denotes the quantity of people who are able to migrate towards the city from certain environment. When the urban population alters, it goes from $x$ to $x^{\text{*}}$ . Considering a certain functions $E\left(x/x_L\right)$ and $b\left(x\right)$ , and using the Boltzmann-type equation and corresponding Fokker-Planck equation, Gualandi and Toscani [6] find that the equilibrium density follows a power law for large cities and a lognormal law for middle and low cities. Yu and Liao [7] employ the non-Maxwellian collision kernel kinetic model. By utilizing the grazing asymptotic method, the corresponding kinetic Fokker-Planck equation is derived. The equilibrium state of this Fokker-Planck equation belongs to the family of generalized Gamma distributions. Numerical tests indicate that the generalized Gamma distribution fits well with the urban size distribution in China.

Inspired by the work in [6] [7], we extend various value functions to explain population changes caused by internal mechanisms. The value function, which satisfies prospect theory, describes the asymmetric changes in population within cities. Using an asymptotic limit approach, we derive the Fokker-Planck equation corresponding to the resulting Boltzmann-type equation, from which we obtain the explicit form of the steady-state solution. Furthermore, we introduce a utility function to describe the changes in the working and non-working populations in cities.

In our work, we establish the ratio of working population to non-working population based on the utility function. Then, according to the interaction rules of the population, the population evolution process is expressed by the Boltzmann equation, and the Fokker-Planck equation is derived through scale transformation. The steady-state solution is obtained under certain conditions, from which steady-state density functions of the working and non-working populations can be derived. Finally, these functions are used to fit real population data.

The structure of this work is as follows. Section 2 introduces the kinetic model. Section 3 discusses the properties of Boltzmann-type equation. Section 4 uses an asymptotic approach to derive the Fokker-Planck equation, which includes a value function. The steady-state distribution of urban population is derived in Section 5. Section 6 gives the numerical simulation. Section 7 summarizes the work.

2. The Kinetic Description

We study the city-size distribution with a kinetic model. Following the idea in [6], we consider that the city population change is related to internal mechanism, external mechanism and fluctuation factor. Suppose that the population distribution in the external environment is known. $g\left(z\right)$ is the probability distribution of variable $z$ . Assume that $g\left(z\right)$ is bounded and satisfies

${\displaystyle {\int }_{R_{+}}g\left(z\right)\text{d}z=1},{\displaystyle {\int }_{R_{+}}{z}^{n}g\left(z\right)\text{d}z}=v_n.$

The phenomenon investigated is the growth process of city size, in which we use a variable $x$ to express the number of urban people. $f\left(x,t\right)$ is the probability density at time $t\ge 0$ . It is generally assumed that the integral of the distribution function $f\left(x,t\right)$ is one, that is

${\displaystyle {\int }_{R_{+}}f\left(x,t\right)\text{d}x}=1,{\displaystyle {\int }_{R_{+}}{x}^{n}f\left(x,t\right)\text{d}x}=m_n\left(t\right).$

For any given value $x$ of urban people, the interaction equation for the change of urban population enloves as

$x^{*}=x+\Phi \left(x/x_L\right)x+b\left(x\right)z+\eta x,$ (2.1)

where $x_L$ is the characteristic city size, e.g. the optimal city size. The left-hand side denotes the population after the change. First item on the right-hand side indicates the initial population change of the city. The second item denotes population change caused by urban internal mechanism. The third item indicates people moving in the city from other place. The last item shows the random factors. Assume that $\eta$ is a random variable with mean zero and variance $\sigma$ .

The core of a city’s optimal population size lies in “balance”. It refers to finding a population scale that aligns with the city’s development stage, on the premise that resources are not overloaded, the environment is not damaged, the economy remains vibrant, and society stays stable. There is no universal formula for it. Instead, it needs to be dynamically calculated based on natural resources, economic structure, technological level, and policy goals, and continuously adjusted as the city develops.

Following the idea in [32], we consider the value function $\Phi \left(x/x_L\right)$ as

$\Phi \left(x/x_L\right)=\mu \frac{1-{\left(x/x_L\right)}^{\delta }}{\left(1+\lambda \right){\left(x/x_L\right)}^{\delta }+1-\lambda },$

in which constants $0<\delta <1$ , $0<\mu <1$ and $0\le \lambda <1$ , which describe the intensity of population change rates near the ideal city size. The value function shows the change rate of urban population caused by psychological factors, so the value function may take positive or negative values. It is positive when the city size $x$ is below the ideal size $x_L$ , and negitive when $x>x_L$ . That is to say, when $x$ is less than $x_L$ , the city’s population is increasing. The other case occurs in the opposite situation. The value fuction $\Phi \left(x/x_L\right)$ equals to zero where $x=x_L$ , namely the total population of the city remains unchanged. Both the value function $\Phi \left(x/x_L\right)$ and function $-E\left(x/x_L\right)$ in Equation (1.1) describe the population change caused by internal mechanism. In agreement with the previous obervations, for any given value of these constants, the function $\Phi \left(x/x_L\right)$ is decreasing and convex on $R_{+}$ . In fact, when $x/x_L\in \left(0,+\infty \right)$ , it takes values in the interval

$-\frac{\mu }{1+\lambda }\le \Phi \left(x/x_L\right)\le \frac{\mu }{1-\lambda },$

which does not depend on the parameter $\delta$ . Based on previous assumptions, we obtain

$0<\frac{\mu }{1+\lambda }<1,$

meaning that the population after migration is greater than 0. In addition, for $0<\Delta x<x_L$ , it holds that

$\Phi \left(1-\Delta x/x_L\right)>-\Phi \left(1+\Delta x/x_L\right).$

This inequality says that given the ideal size of a city, the increase in population in a city with a population below the ideal size is greater than the decrease in the number of residents in a city with a population above the ideal size. When $x<x_L$ , as $x$ increases, the value function $\Phi \left(x/x_L\right)$ grows more and more slowly, namely, the curve is concave. When $x>x_L$ , as $x$ increases, the value function $\Phi \left(x/x_L\right)$ decreases faster and faster, namely, the curve is convex. We compare the value function of different parameters in Figure 1.

Figure 1. The value function $\Phi \left(x/x_L\right)$ .

According to the interaction rule (2.1) and the classical dynamic theory [33], the change of probability density $f\left(x,t\right)$ satisfies a weak form of Boltzmann-type model, namely, for all smooth function $\phi \left(x\right)$ , the $f\left(x,t\right)$ obeys

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{R_{+}}f\left(x,t\right)\phi \left(x\right)\text{d}x}=\langle {\displaystyle {\int }_{R_{+}^2}\left(\phi \left(x^{\text{*}}\right)-\phi \left(x\right)\right)f\left(x,t\right)g\left(z\right)\text{d}x\text{d}z}\rangle ,$ (2.2)

where the bracket $\langle \rangle$ represents the mathematical expectation. Equation (2.2) reflects the population distribution evolution of the city caused by the interaction at $t>0$ .

3. Properties of the Kinetic Model

Since the interaction rule (2.1) is nonlinear, it is difficult to study the time-behavior of the probability density $f\left(x,t\right)$ . Because of this, we analyze the main features of the kinetic model.

In the Boltzmann-type model (2.2), letting $\phi \left(x\right)=1$ , we acquire

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{R_{+}}f\left(x,t\right)\text{d}x}=0,$

from which we obtain that the Equation (2.2) satisfies the conservation property. We study population changes caused by internal mechanisms of cities. In the latter part of the paper, we assume that $b\left(x\right)=b$ in interaction rule (2.1). Since $\langle \eta \rangle =0$ , from the interaction rule (2.1), we have

$x^{*}-x=\Phi \left(x/x_L\right)x+bz+\eta x.$ (3.1)

Then, considering $\phi \left(x\right)=x$ in Equation (2.2), mean value meets

$\frac{\text{d}}{\text{d}t}{m}_1\left(t\right)=\langle {\displaystyle {\int }_{R_{+}^2}\left(x^{\text{*}}-x\right)f\left(x,t\right)g\left(z\right)\text{d}x\text{d}z}\rangle ,$ (3.2)

where

$m_1\left(t\right)={\displaystyle {\int }_{R_{+}}xf\left(x,t\right)\text{d}x}.$

Substituting (3.1) into (3.2) yields

$\frac{\text{d}}{\text{d}t}{m}_1\left(t\right)=\langle {\displaystyle {\int }_{R_{+}^2}\left(\Phi \left(x/x_L\right)x+bz+\eta x\right)f\left(x,t\right)g\left(z\right)\text{d}x\text{d}z}\rangle .$ (3.3)

Choosing $\Phi \left(x/x_L\right)=s$ in (3.3), which is the same as in [34], we have

$\frac{\text{d}{m}_1\left(t\right)}{\text{d}t}=s{m}_1\left(t\right)+b{v}_1.$

Using the constant variation method of differential equations, yields

$m_1\left(t\right)=-\frac{b{v}_1}{s}+\left(m_1\left(0\right)+\frac{b{v}_1}{s}\right)\exp \left(st\right).$

If a city population $x>x_L$ , namely urban population $x$ is more than ideal urban population $x_L$ . Since $s<0$ , the $m_1\left(t\right)$ satisfies

$m_1\left(t\right)\to -\frac{b{v}_1}{s},\left(t\to \infty \right),$

implying that the first moment $-\frac{b{v}_1}{s}$ remains bounded if $t\to \infty$ .

Considering $\phi \left(x\right)=x^2$ in Equation (2.2) gives rise to

$\frac{\text{d}}{\text{d}t}{m}_2\left(t\right)=\langle {\displaystyle {\int }_{R_{+}^2}\left({\left(x^{\text{*}}\right)}^2-x^2\right)f\left(x,t\right)g\left(z\right)\text{d}x\text{d}z}\rangle ,$

where

$m_2\left(t\right)={\displaystyle {\int }_{R_{+}}{x}^{2}f\left(x,t\right)\text{d}x}.$

Since $\langle \eta \rangle =0$ and ${\eta }^2=\sigma$ , from the interaction rule (2.1), we obtain

$\langle {\left(x^{*}\right)}^2-x^2\rangle =x^2\left[2\Phi \left(x/x_L\right)+{\Phi }^2\left(x/x_L\right)+\sigma \right]+2xbz\Phi \left(x/x_L\right)+b^2{z}^2,$

which leads to

$\frac{\text{d}}{\text{d}t}{m}_2\left(t\right)={\displaystyle {\int }_{R_{+}^2}\left[x^2\left(2s+s^2+\sigma \right)+2xbz\left(1+s\right)+b^2{z}^2\right]f\left(x,t\right)g\left(z\right)\text{d}x\text{d}z}$

and

$\frac{\text{d}}{\text{d}t}{m}_2\left(t\right)=\left(2s+s^2+\sigma \right)m_2\left(t\right)+2b\left(1+s\right)v_1{m}_1\left(t\right)+b^2{v}_2.$

If $t\to \infty$ , $m_1\left(t\right)$ and $b^2{v}_2$ remain bounded. Therefore the second moment $m_2\left(t\right)$ keeps bounded as $t$ tends to infinite.

4. Fokker-Planck Equation

As discussed above, it is difficult to find the analytical solution of Equation (2.2). A feasible approach is to investigate its large-time behavior, which can be found in [34]. Considering the kinetic model, we use the asymptotic limit method to derive the Fokker-Planck model relating to the Boltzmann-type equation.

In terms of the value function, letting $\epsilon =\delta$ associated with $0<\epsilon \ll 1$ , $\Phi \left(x/x_L\right)$ is expressed as

$\Phi \left(x/x_L\right)\to {\Phi }^{\epsilon }\left(x/x_L\right),$

where

${\Phi }^{\epsilon }\left(x/x_L\right)=\mu \frac{1-{\left(x/x_L\right)}^{\epsilon }}{\left(1+\lambda \right){\left(x/x_L\right)}^{\epsilon }+1-\lambda },$ (4.1)

from which we obtain

$\underset{\epsilon \to 0^{+}}{\lim }{\Phi }^{\epsilon }\left(x/x_L\right)=\frac{-\mu }{1+\lambda }.$

At this point, it indicates that the rate of population change within a city is only related to two parameters $\mu$ and $\lambda$ . When $\epsilon \to 0$ , it means that the internal mechanism of the city leads to a decrease for the population of the city.

We state that $f\left(x,\tau \right)$ with $t\to \frac{\tau }{\epsilon }$ is independent of $\epsilon$ . Therefore, for all

smooth function $\phi \left(x\right)$ , the distribution function $f\left(x,t\right)$ obeys the Boltzmann-type model, which reads as

$\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{R_{+}}f\left(x,\tau \right)\phi \left(x\right)\text{d}x}=\frac{1}{\epsilon }\langle {\displaystyle {\int }_{R_{+}^2}\left(\phi \left(x^{*}\right)-\phi \left(x\right)\right)f\left(x,\tau \right)g\left(z\right)\text{d}x\text{d}z}\rangle .$ (4.2)

Reconsidering the interaction rule (2.1), $\eta \to \sqrt{\epsilon }\eta$ and $b\to \epsilon b$ , we obtain

$\langle x^{*}-x\rangle =\epsilon {\Phi }^{\epsilon }\left(x/x_L\right)x+\epsilon bz$

and

$\langle {\left(x^{*}-x\right)}^2\rangle ={\epsilon }^2{\left[{\Phi }^{\epsilon }\left(x/x_L\right)x\right]}^2+\epsilon \sigma x^2+{\epsilon }^2{b}^2{z}^2+{\epsilon }^2{\Phi }^{\epsilon }\left(x/x_L\right)bz.$

Using Taylor expansion of $\phi \left(x\right)$ around $x$ , we have

$\phi \left(x^{\text{*}}\right)-\phi \left(x\right)=\left(x^{\text{*}}-x\right){\phi }^{\prime }\left(x\right)+\frac{1}{2}{\left(x^{\text{*}}-x\right)}^2{\phi }^{″}\left(x\right)+\frac{1}{3!}{\left(x^{\text{*}}-x\right)}^3{\phi }^{‴}\left(\tilde{x}\right),$(4.3)

where $\tilde{x}=x+\theta \left(x^{\text{*}}-x\right)$, $0\le \theta \le 1$ .

Therefore, using (4.2), (4.3) and noticing the power of $\epsilon$ yield

$\begin{array}{l}\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{R_{+}}f\left(x,\tau \right)\phi \left(x\right)\text{d}x}\\ ={\displaystyle {\int }_{R_{+}^2}\left[\left({\Phi }^{\epsilon }\left(x/x_L\right)+bz\right]x{\phi }^{\prime }\left(x\right)+\frac{1}{2}\sigma x^2{\phi }^{″}\left(x\right)\right)f\left(x,\tau \right)g\left(z\right)\text{d}x\text{d}z}+R_{ϵ}\left(\tau \right),\end{array}$

where

$\begin{array}{c}{R}_{ϵ}\left(\tau \right)=\frac{1}{\epsilon }\langle {\displaystyle {\int }_{R_{+}}({\epsilon }^2\left({\left({\Phi }^{\epsilon }\left(x/x_L\right)x\right)}^2+{\epsilon }^2{b}^2{z}^2+{\epsilon }^2{\Phi }^{\epsilon }\left(x/x_L\right)bz\right){\phi }^{″}\left(x\right)\stackrel{}{}}\\ +\frac{1}{6}{\left(x^{\text{*}}-x\right)}^3{\phi }^{‴}\left(x\right))f\left(x,\tau \right)g\left(z\right)\text{d}x\text{d}z\rangle .\end{array}$

Taking the limit of the term $R_{ϵ}\left(\tau \right)$ gives rise to

$\underset{\epsilon \to 0}{\text{lim}}{R}_{ϵ}\left(\tau \right)=0.$

Thus, we obtain

$\begin{array}{l}\frac{\text{d}}{\text{d}\tau }{\displaystyle {\int }_{R_{+}}f\left(x,\tau \right)\phi \left(x\right)\text{d}x}\\ ={\displaystyle {\int }_{R_{+}^2}\left(\left({\Phi }^{\epsilon }\left(x/x_L\right)x+bz\right){\phi }^{\prime }\left(x\right)+\frac{1}{2}\sigma x^2{\phi }^{″}\left(x\right)\right)f\left(x,\tau \right)g\left(z\right)\text{d}x\text{d}z}.\end{array}$ (4.4)

When $\epsilon$ is small enough, the value function (4.1) is a particular case of the general class of value functions defined by

${\Phi }^{\epsilon }\left(x/x_L\right)=\mu \frac{1-{\text{e}}^{\epsilon \left({\left(x/x_L\right)}^{\delta }-1\right)/\delta }}{\left(1+\lambda \right){\text{e}}^{\epsilon \left({\left(x/x_L\right)}^{\delta }-1\right)/\delta }+1-\lambda }.$ (4.5)

Letting $\delta \to 0$ , value function (4.5) is the function (4.1). Moreover, we consider that the quasi-invariant scaling [31]. Indeed, $\epsilon \to 0$ , the value function satisfies

${\Phi }^{\epsilon }\left(x/x_L\right)\approx -\mu \epsilon \frac{\left({\left(x/x_L\right)}^{\delta }-1\right)/\delta }{\left(1+\lambda \right)\epsilon \left(1-{\left(x/x_L\right)}^{\delta }\right)/\delta +2},$

which implies

$\underset{\epsilon \to 0}{\text{lim}}\frac{{\text{Φ}}^{\epsilon }\left(x/x_L\right)}{\epsilon }=\frac{\mu }{2\delta }\left(1-{\left(x/x_L\right)}^{\delta }\right).$

Furthermore, we consider the quasi-invariant scaling for the introduced value function

$\underset{\epsilon \to 0}{\lim }\frac{{\Phi }^{\epsilon }\left(x/x_L\right)}{\epsilon }=\Phi \left(x/x_L\right).$ (4.6)

Hence, if we consider the value function $\frac{\mu }{2\delta }\left(1-{\left(x/x_L\right)}^{\delta }\right)$ , for any given $0<\delta <1$ , the kinetic Equation (2.2) in the quasi-invariant limit is given by the Fokker-Planck model.

$\frac{\partial f\left(x,\tau \right)}{\partial \tau }=\frac{\mu }{2\delta }\frac{\partial \left[\left(\frac{\mu }{2\delta }\left({\left(x/x_L\right)}^{\delta }-1\right)x+b{v}_1\right)f\left(x,\tau \right)\right]}{\partial x}+\frac{\sigma }{2}\frac{\partial \left[x^{2}f\left(x,\tau \right)\right]}{\partial x^2}.$ (4.7)

5. The Stationary Distribution

Generally speaking, the analytical solutions of the Fokker-Planck Equation (4.7) is difficult to be obtained. Thus, we use the steady-state of the model (4.7) to explain the asymptotic features of the city-size evolution. The stationary distribution of model (4.7) is found by solving the differential equation.

5.1. Urban Population Size

Moreover, the equilibrium state of Equation (4.7) is given by

$\left[-\frac{\mu }{2\delta }\left({\left(x/x_L\right)}^{\delta }-1\right)x+b{v}_1\right]f\left(x,\tau \right)=\frac{\sigma }{2}\frac{\partial \left[x^{2}f\left(x,\tau \right)\right]}{\partial x},$

from which, we obtain

$f_{\infty }\left(x\right)=c{x}^{\frac{\mu }{\delta \sigma }-2}\exp \left[-\frac{\mu }{{\delta }^2\sigma }{\left(x/x_L\right)}^{\delta }-\frac{2b{v}_1}{\sigma x}\right],$ (5.1)

where the constant $c$ guarantees that the mass of the steady-state density is $1$ .

When $\frac{\mu }{\delta \sigma }-2<0$ , we find that the equilibrium distribution is a polynomial rate decay with respect to $x$ at infinity. The steady distribution (5.1) has a power law tail.

When $\frac{\mu }{\delta \sigma }-2\ge 0$ , the steady-state distribution is exponential decay with respect to $x$ at infinity.

Note that, if $x$ denotes the number of cancer cells, different values about $-1\le \delta \le 1$ are considered in [32].

5.2. Working People and Non-Working People

We assume that a city has a certain proportion of the working population, and the rest people of city are not working. The working people are the persons who are paid for their labor. The city has a population of $x$ . Although the number of urban population $x$ is obviously a natural number, in order to avert unnecessary difficulties, we hypothesise $x\in R^{+}$ in the rest of the paper (see [6]). The working population $q$ who provide labor in exchange for remuneration is given by

$q=\gamma x,$

where $\gamma$ represents the proportion of the working population in the city, $0<\gamma <1$ .

The non-working population in cities includes infants and young children who need care, students who are studying, adults who are not employed, and elderly people without jobs. The non-working population of the city $h$ evolves as

$h=\kappa x,$

where $\kappa$ represents the proportion of the non-working population in the city, $0<\kappa <1$ .

A city’s population is divided into working people and non-working people. Therefore

$\gamma +\kappa =1.$

The non-working people do not work, but they increase the consumption of the entire city. The development of a city is related to both working and non-working population. We use utility function to represent the preference of urban development. The working population can become non-working people, and some non-working people can become working people. We think about the constant elasticity of substitution (CES) utility function to describe the preference of urban development. Following the work in [35] [36], we use the following CES utility function

$U={\left(\alpha q^{\xi }+\beta h^{\xi }\right)}^{\frac{1}{\xi }},\alpha +\beta =1,$ (5.2)

where $0<\xi <1,0<\alpha ,\beta <1$ . When $\xi \to 1$ , the utility function (5.2) denotes that the working population and non-working population are completely interchangeable. Namely, the utility function (5.2) converges to the linear utlity function. When $\xi \to 0$ , the utility function (5.2) indicates that the working population and non working population can be converted into each other with a certain proportion. Namely, the utility function (5.2) converges to the Cobb-Douglas utility function. Different cities have different preferences for the development of cities. The optimal employment structure when a city achieves a balance between production attractiveness and residential livability. For simplicity, we assume that $\xi =1/2$ for the rest of the study. The values $\alpha$ and $\beta$ are linked to the preferences that the city assigns to the two kinds of people. If $\alpha >\beta$ , the city prefers to possess working people. Otherwise, we have an opposite result. The choice $\alpha =\beta =1/2$ means that the two kinds of people are equally important for a city. The total population of the city is $x$ ,

$q+h=x.$

Suppose that $\xi =\frac{1}{2}$ , the utility function (5.2) becomes

$U={\left(\alpha q^{\frac{1}{2}}+\beta {\left(x-q\right)}^{\frac{1}{2}}\right)}^2,\alpha +\beta =1.$

When the utility function is maximized, its derivative is zero.

$\frac{\text{d}U}{\text{d}q}=2\left(\alpha q^{\frac{1}{2}}+\beta {\left(x-q\right)}_{}^{\frac{1}{2}}\right)\left(\frac{1}{2}\alpha q^{-\frac{1}{2}}-\frac{1}{2}\beta {\left(x-q\right)}_{}^{-\frac{1}{2}}\right)=0,$

from which, we get

$\frac{1}{2}\alpha q^{-\frac{1}{2}}-\frac{1}{2}\beta {\left(x-q\right)}_{}^{-\frac{1}{2}}=0.$

Obtaining the maximum utility function, the ratio of working population to non-working population is as follows

$h=\frac{{\beta }^2}{{\alpha }^2}q.$

When the utility function is maximized, the working population $q$ and no-working population $h$ are

$\xi =\frac{1}{2}q=\frac{{\alpha }^2}{{\alpha }^2+{\beta }^2}xh=\frac{{\beta }^2}{{\beta }^2+{\alpha }^2}x.$ (5.3)

From Equation (5.3), when urban utility reaches its maximum, the proportion of urban working population in the whole urban population is $\frac{{\alpha }^2}{{\alpha }^2+{\beta }^2}$ . And the ratio of no-working population in the whole population is $\frac{{\beta }^2}{{\alpha }^2+{\beta }^2}$ .

Following the previous method, we get the distribution of the two types of population when the city’s utility is at its maximum. From Equation (5.3), we obtain

$x=\frac{{\beta }^2+{\alpha }^2}{{\alpha }^2}q=\frac{{\beta }^2+{\alpha }^2}{{\beta }^2}h.$

The optimal distribution of working population is $q=\frac{{\alpha }^2}{{\alpha }^2+{\beta }^2}x$ and non-working population is $h=\frac{{\beta }^2}{{\beta }^2+{\alpha }^2}x$ . The probability densities $f_Q\left(q\right)$ and $f_H\left(h\right)$ of random variables $Q=\frac{{\alpha }^2}{{\alpha }^2+{\beta }^2}X$ and $H=\frac{{\beta }^2}{{\beta }^2+{\alpha }^2}X$ take the form¹

$\begin{array}{l}{f}_Q\left(q\right)\\ =c\frac{{\beta }^2+{\alpha }^2}{{\alpha }^2}{\left(\frac{{\beta }^2+{\alpha }^2}{{\alpha }^2}q\right)}^{\frac{\mu }{\delta \sigma }-2}\exp \left[-\frac{\mu }{{\delta }^2\sigma }{\left(\left(\frac{{\beta }^2+{\alpha }^2}{{\alpha }^2}q\right)/x_L\right)}^{\delta }-\frac{2{\alpha }^{2}b{v}_1}{\sigma \left({\beta }^2+{\alpha }^2\right)q}\right]\end{array}$ (5.4)

and

$\begin{array}{l}{f}_H\left(h\right)\\ =c\frac{{\beta }^2+{\alpha }^2}{{\beta }^2}{\left(\frac{{\beta }^2+{\alpha }^2}{{\beta }^2}h\right)}^{\frac{\mu }{\delta \sigma }-2}\exp \left[-\frac{\mu }{{\delta }^2\sigma }{\left(\left(\frac{{\beta }^2+{\alpha }^2}{{\beta }^2}h\right)/x_L\right)}^{\delta }-\frac{2{\beta }^{2}b{v}_1}{\sigma \left({\beta }^2+{\alpha }^2\right)h}\right].\end{array}$ (5.5)

From (5.1), (5.4) and (5.5), we obtain that the city-size distribution follows a power law. The optimal distribution of working population $f_Q\left(q\right)$ and non-working population $f_H\left(h\right)$ also follow a power law.

6. Numerical Experiments

We conduct the numerical experiments in two parts. In the first part, we observe the image of the steady-state solution of urban population distribution; in the second part, we carry out an analysis in combination with China’s urban population data.

6.1. Test 1

The empirical research on the power law distribution of urban city-size is presented in [9]. When utility function is at its maximum, both the distribution of working population and non-working population follow a power law. Thus, in this part, we mainly consider the impact of urban internal mechanisms on the urban population, such as the psychological state of the urban population. Therefore, we ignore the influence of external factors in the city and assume that the immigration rate of the external environment of the city $b$ is zero. We choose Equation (5.1) to describe the stationary city-size distribution, namely,

$f_{\infty }\left(x\right)=c{x}^{\frac{\mu }{\delta \sigma }-2}\exp \left[-\frac{\mu }{{\delta }^2\sigma }{\left(x/x_L\right)}^{\delta }\right].$ (6.1)

Since ${\displaystyle {\int }_{R_{+}}{f}_{\infty }\left(x\right)\text{d}x}=1$ . Considering $y=\frac{\mu }{{\delta }^2\sigma }{\left(x/x_L\right)}^{\delta }$ , we have $x={\left[\frac{{\delta }^2\sigma }{\mu }{x}_L^{\delta }y\right]}^{\frac{1}{\sigma }}$ and $\text{d}x=\frac{\delta \sigma }{\mu }\left({\left(\frac{{\delta }^2\sigma }{\mu }{x}_L^{\delta }y\right)}^{\frac{1}{\delta }-1}{x}_L^{\delta }\right)\text{d}y$ . We acquire

$c\frac{1}{\delta }{\left[\frac{{\delta }^2\sigma }{\mu }{x}_L^{\delta }\right]}^{\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }}{\displaystyle {\int }_{R_{+}}{y}^{\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }-1}{\text{e}}^{-y}\text{d}y}=1.$

Using $\Gamma \left(\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }\right)={\displaystyle {\int }_0^{+\infty }{y}^{\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }-1}{\text{e}}^{-y}\text{d}y}$ yields

$c=\frac{\delta }{{\left[\frac{{\delta }^2\sigma }{\mu }{x}_L^{\delta }\right]}^{\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }}\Gamma \left(\frac{\mu }{{\delta }^2\sigma }-\frac{1}{\delta }\right)}.$

We choose different values of the parameters of the stationary distribution of urban people, which include working people, and compare several numerical results.

Firstly, we consider different stationary city-size distribution with different values of $\sigma$ , which describes the contribution of diffusion. We suppose that the ideal urban population $x_L=1000000$ and $x_L=2000000$ . The probability density distribution under the steady distribution of the urban population size is shown in Figure 2 when $\mu =0.5$ , $\delta =0.5$ . The steady-state profile $f_{\infty }\left(x\right)$ corresponds to $\sigma =0.4$ , $\sigma =0.3$ , $\sigma =0.2$ . For each value of $\sigma$ , we find that each steady-state profile $f_{\infty }\left(x\right)$ has a unique vertex. The steady distribution of different $\sigma$ with $x_L=1000000$ in Figure 2(a), and the steady distribution of different $\sigma$ with $x_L=2000000$ in Figure 2(b). We conclude that any decrease of $\sigma$ leads to the stationary distribution with lower vertex and more fatter tails.

Secondly, we consider the probability density of urban population size distribution corresponding to different $\delta$ , which represents the contribution of drift. In this case, we also assume that the ideal urban population $x_L=1000000$ and $x_L=2000000$ . The steady distribution of the urban population size is shown in Figure 3 when $\sigma =0.4$ , $\mu =0.5$ . The steady-state profile $f_{\infty }\left(x\right)$ corresponds to $\delta =0.5$ , $\delta =0.4$ , $\delta =0.3$ . For every $\sigma$ , we find that the corresponding stationary profile $f_{\infty }\left(x\right)$ possesses a unique vertex. The steady distribution of different $\sigma$ with $x_L=1000000$ in Figure 3(a), and the steady distribution of different $\sigma$ with $x_L=2000000$ in Figure 3(b). We conclude that any decrease of

Figure 2. Behaviors of the steady-state probability density distribution of urban population size for $\sigma =0.4$ (blue line), $\sigma =0.3$ (black line), $\sigma =0.2$ (red line). (a) $x_L=1000000$ , (b) $x_L=2000000$ .

Figure 3. Behaviors of the steady-state probability density distribution of urban population size for $\delta =0.5$ (blue line), $\delta =0.4$ (black line), $\delta =0.3$ (red line). (a) $x_L=1000000$ , (b) $x_L=2000000$ .

$\sigma$ leads to the stationary distribution with lower vertex and more fatter tails. Thirdly, we consider the probability density of urban population size distribution corresponding to different $\mu$ , which represents the contribution of drift. In this case, we also let ideal urban population $x_L=1000000$ and $x_L=2000000$ . The steady distribution of the urban population size is shown in Figure 4 when $\delta =0.5$ , $\sigma =0.4$ . The stationary profile $f_{\infty }\left(x\right)$ corresponds to $\mu =0.5$ , $\mu =0.48$ , $\mu =0.46$ . For every $\mu$ , we find that the corressponding profile $f_{\infty }\left(x\right)$ possesses a unique vertex. The steady distribution of different $\sigma$ with $x_L=1000000$ in Figure 4(a), and the steady distribution of different $\sigma$ with $x_L=2000000$ in Figure 4(b). We conclude that any decrease of $\mu$ leads to the stationary distribution with lower vertex and more fatter tails.

Figure 4. Behaviors of the steady-state probability density distribution of urban population size for $\mu =0.5$ (blue line), $\mu =0.48$ (black line), $\mu =0.46$ (red line). (a) $x_L=1000000$ , (b) $x_L=2000000$ .

Fourthly, we consider the stationary profile $f_{\infty }\left(x\right)$ when $\mu =0.5$ , $\delta =0.5$ and $\sigma =0.5$ . Seeing from Figure 5, the steady distribution has no vertex. This extreme case implies that stationary distribution possesses thin tail.

Figure 5. Behaviors of the steady-state probability density distribution of urban population size for $\mu =0.5$ , $\sigma =0.5$ and $\delta =0.5$ .

Finally, when utility function is at its maximum, we consider the probability density of working population corresponding to different $\alpha$ and $\beta$ . The steady distribution of working population with $\delta =0.5$ , $\sigma =0.4$ , $\mu =0.5$ and $x_L=2000000$ is shown in Figure 6. We conclude that any increase of $\alpha$ leads to the stationary distribution with lower vertex and more fatter tails.

Figure 6. Behaviors of the steady-state probability density distribution of urban working population for $\alpha =1/2$ , $\alpha =1/3$ , $\alpha =1/4$ , $\alpha =2/3$ , $\alpha =3/4$ .

6.2. Test 2

A city’s development is inseparable from its education, medical care, and industries. Therefore, the working population and non-working population are complementary to each other. If a city is rich in educational resources, it will attract many families with children. If a city has excellent medical resources, it will attract many families with members suffering from illnesses. If a city has a well-developed financial industry, it will be highly attractive to financial talents.

We use the 2024 population data from China Statistical Yearbook, which contains the population data of 2023. We list the total population (unit: 10,000 people), working population (unit: 10,000 people), and gross regional product (GRP unit: 100 million yuan) of several cities in Table 1.

Table 1. Urban population and gross regional product.

Region	Total population	Working people	Gross regional product
Beijing	2186	1129	43760.7
Shanghai	2467	1345	47218.7
Chongqing	3043	1662	30145.8
Hebei	7323	3623	43944.1
Shanxi	3528	1704	25698.2
Liaoning	4358	2091	30209.4

The three cities are municipalities directly under the Central Government of China, each with its own characteristics. Beijing is the political and cultural center, Shanghai is the economic and financial center, and Chongqing is a hub in western China. After calculation, we derive the ratio of parameter $\alpha$ to parameter $\beta$ as shown in Table 2.

Table 2. The ratio of the working population to the non-working population.

Region	$\alpha /\beta$	Average GRP
Beijing	1.03	20.02
Shanghai	1.09	19.14
Chongqing	1.1	9.91
Hebei	0.99	6.00
Shanxi	0.97	7.28
Liaoning	0.96	6.93

Through observation, it is found that the ratio of the working population to the non-working population is approximately 1. Cities with a relatively high average gross regional product (Average GRP) have a ratio of working population to non-working population greater than 1; furthermore, the higher the average GRP, the closer this ratio is to 1.

7. Conclusion

In this work, we employ a kinetic model to investigate the evolution of the urban population. We assume that the urban population is made up of working and non-working people. The model contains a utility function to represent the preference for city development. When the working population and the non-working population reach an appropriate proportion, the city reaches its optimal development. We consider that a value function describes the change rate of the urban population caused by psychological factors. We describe the interaction rules of urban population change, and give the Boltzmann-type equation and the corresponding Fokker-Planck equation. We obtain the stationary city-size distribution, and find that when urban utility reaches its maximum, both the working population and the non-working population distribution obey the power law. Nevertheless, this study has several limitations: the simplification of the population into two groups, the ($b=0$ ) assumption adopted in Test 1, and the lack of empirical identification of the external distribution g(z). These issues warrant further investigation and improvement in future research.

Funding

This work is supported by the Key Research Base Project of Philosophy and Social Sciences in Aba Prefecture (Grant No.ABKT-CXBLY2025-34).

NOTES

¹If the probability density of a continuous random variable $X$ is $f_X\left(x\right)$ , ($-\infty <x<+\infty$ ), the function $y=s\left(x\right)$ is differentiable everywhere and strictly monotonic, then $Y=s\left(X\right)$ is a continuous random variable and the probability density $f_Y\left(y\right)$ is $f_Y\left(y\right)=\{\begin{array}{ll}{f}_X\left[s^{-1}\left(y\right)\right]\left|\frac{\text{d}{s}^{-1}\left(y\right)}{\text{d}y}\right|,\hfill & a<y<r,\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$ where $a=\min \left\{s\left(-\infty \right),s\left(+\infty \right)\right\}$ , $r=\max \left\{s\left(-\infty \right),s\left(+\infty \right)\right\}$ .

Conflicts of Interest

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

References

[1]	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.[CrossRef]
[2]	Benguigui, L. and Blumenfeld-Lieberthal, E. (2007) Beyond the Power Law—A New Approach to Analyze City Size Distributions. Computers, Environment and Urban Systems, 31, 648-666.[CrossRef]
[3]	Garmestani, A.S., Allen, C.R. and Gallagher, C.M. (2008) Power Laws, Discontinuities and Regional City Size Distributions. Journal of Economic Behavior & Organization, 68, 209-216.[CrossRef]
[4]	Eeckhout, J. (2004) Gibrat’s Law for (All) Cities. American Economic Review, 94, 1429-1451.[CrossRef]
[5]	Eeckhout, J. (2009) Gibrat’s Law for (All) Cities: Reply. American Economic Review, 99, 1676-1683.[CrossRef]
[6]	Gualandi, S. and Toscani, G. (2019) Size Distribution of Cities: A Kinetic Explanation. Physica A: Statistical Mechanics and Its Applications, 524, 221-234.[CrossRef]
[7]	Yu, L. and Liao, J. (2021) Non-Maxwellian Kinetic Modelling of City Size Distribution. Journal of Applied Mathematics and Physics, 9, 1329-1339.[CrossRef]
[8]	Calderín-Ojeda, E. (2016) The Distribution of All French Communes: A Composite Parametric Approach. Physica A: Statistical Mechanics and its Applications, 450, 385-394.[CrossRef]
[9]	Wu, G.Z. and Yang, J.W. (2017) The Empirical Study of Chinese City Rank-Size Power Law Distribution. Economic Geography, 37, 59-65. (In Chinese)
[10]	Zipf, G.K. (1949) Human Behavior and the Principle of Least Effort, Addison-Wesley.
[11]	Gabaix, X. (1999) Zipf’s Law for Cities: An Explanation. The Quarterly Journal of Economics, 114, 739-767.[CrossRef]
[12]	Marsili, M. and Zhang, Y. (1998) Interacting Individuals Leading to Zipf’s Law. Physical Review Letters, 80, 2741-2744.[CrossRef]
[13]	Ghosh, A., Chatterjee, A., Chakrabarti, A.S. and Chakrabarti, B.K. (2014) Zipf’s Law in City Size from a Resource Utilization Model. Physical Review E, 90, Article ID: 042815.[CrossRef] [PubMed]
[14]	Gangopadhyay, K. and Basu, B. (2009) City Size Distributions for India and China. Physica A: Statistical Mechanics and Its Applications, 388, 2682-2688.[CrossRef]
[15]	González-Val, R. (2020) The Spanish Spatial City Size Distribution. Environment and Planning B: Urban Analytics and City Science, 48, 1609-1631.[CrossRef]
[16]	Li, W. and Zhang, Y. (2018) Triangle Law or Power Law. City Size Distribution in Sub-National Levelled Administrative Areas in China. International Review for Spatial Planning and Sustainable Development, 6, 203-215.[CrossRef]
[17]	Levy, M. (2009) Gibrat’s Law for (All) Cities: Comment. American Economic Review, 99, 1672-1675.[CrossRef]
[18]	Zhang, W., Yang, D. and Huo, J. (2016) Studies of the Relationship between City Size and Urban Benefits in China Based on a Panel Data Model. Sustainability, 8, Article No. 554.[CrossRef]
[19]	Xu, Z. and Zhu, N. (2009) City Size Distribution in China: Are Large Cities Dominant? Urban Studies, 46, 2159-2185.[CrossRef]
[20]	Brenner, N. and Theodore, N. (2002) Cities and the Geographies of “Actually Existing Neoliberalism”. Antipode, 34, 349-379.[CrossRef]
[21]	Jacobs, J.M. (2012) Urban Geographies I: Still Thinking Cities Relationally. Progress in Human Geography, 36, 412-422.[CrossRef]
[22]	Cordier, S., Pareschi, L. and Piatecki, C. (2009) Mesoscopic Modelling of Financial Markets. Journal of Statistical Physics, 134, 161-184.[CrossRef]
[23]	During, B. and Toscani, G. (2008) International and Domestic Trading and Wealth Distribution. Communications in Mathematical Sciences, 6, 1043-1058.[CrossRef]
[24]	Ghosh, A., Chatterjee, A., Inoue, J. and Chakrabarti, B.K. (2016) Inequality Measures in Kinetic Exchange Models of Wealth Distributions. Physica A: Statistical Mechanics and its Applications, 451, 465-474.[CrossRef]
[25]	Torregrossa, M. and Toscani, G. (2018) Wealth Distribution in Presence of Debts: A Fokker-Planck Description. Communications in Mathematical Sciences, 16, 537-560.[CrossRef]
[26]	Brugna, C. and Toscani, G. (2015) Kinetic Models of Opinion Formation in the Presence of Personal Conviction. Physical Review E, 92, Article ID: 052818.[CrossRef] [PubMed]
[27]	Toscani, G. (2006) Kinetic Models of Opinion Formation. Communications in Mathematical Sciences, 4, 481-496.[CrossRef]
[28]	Lallouache, M., Chakrabarti, A.S., Chakraborti, A. and Chakrabarti, B.K. (2010) Opinion Formation in Kinetic Exchange Models: Spontaneous Symmetry-Breaking Transition. Physical Review E, 82, Article ID: 056112.[CrossRef] [PubMed]
[29]	Boudin, L., Monaco, R. and Salvarani, F. (2010) Kinetic Model for Multidimensional Opinion Formation. Physical Review E, 81, Article ID: 036109.[CrossRef] [PubMed]
[30]	Furioli, G., Pulvirenti, A., Terraneo, E. and Toscani, G. (2017) Fokker-Planck Equations in the Modeling of Socio-Economic Phenomena. Mathematical Models and Methods in Applied Sciences, 27, 115-158.[CrossRef]
[31]	Gualandi, S. and Toscani, G. (2018) Call Center Service Times Are Lognormal: A Fokker-Planck Description. Mathematical Models and Methods in Applied Sciences, 28, 1513-1527.[CrossRef]
[32]	Preziosi, L., Toscani, G. and Zanella, M. (2021) Control of Tumor Growth Distributions through Kinetic Methods. Journal of Theoretical Biology, 514, Article ID: 110579.[CrossRef] [PubMed]
[33]	Pareschi, L. and Toscani, G. (2013) Interacting Multiagent Systems: Kinetic Equations and Monte Carlo Methods. Oxford University Press.
[34]	Zhong, Y., Lai, S. and Hu, C. (2021) Investigations to the Dynamics of Wealth Distribution in a Kinetic Exchange Model. Applied Mathematics and Computation, 404, Article ID: 126231.[CrossRef]
[35]	Brugna, C. and Toscani, G. (2018) Kinetic Models for Goods Exchange in a Multi-Agent Market. Physica A: Statistical Mechanics and its Applications, 499, 362-375.[CrossRef]
[36]	Hu, C., Lai, S. and Lai, C. (2020) Investigations to the Price Evolutions of Goods Exchange with CES Utility Functions. Physica A: Statistical Mechanics and Its Applications, 549, Article ID: 123938.[CrossRef]