Application of kinetic models and simulation of multi-particle systems in some fields of science Текст научной статьи по специальности «Математика»

Электронный научный журнал "Математическое моделирование, компьютерный и натурный эксперимент в естественных науках" http://mathmod.esrae.ru/ URL статьи: mathmod.esrae.ru/24-79 Ссылка для цитирования этой статьи:

Burmistrov A.V., Korotchenko M.A. Application of Kinetic Models and Simulation of Multi-particle Systems in Some Fields of Science // Математическое моделирование, компьютерный и натурный эксперимент в естественных науках. 2019. №1

Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект 0315-2016-0002) и при частичной поддержке РФФИ (гранты 17-01-00698, 18-01-00356, 18-01-00599)_

УДК 519.6, 519.2

ПРИМЕНЕНИЕ КИНЕТИЧЕСКИХ МОДЕЛЕЙ И МОДЕЛИРОВАНИЯ МНОГОЧАСТИЧНЫХ СИСТЕМ В НЕКОТОРЫХ ОБЛАСТЯХ НАУКИ

Бурмистров А.В.1, Коротченко М.А.2 1Институт вычислительной математики и математической геофизики СО РАН;

Новосибирский госуниверситет, Россия, Новосибирск, burm@osmf.sscc.ru 2Институт вычислительной математики и математической геофизики СО РАН,

Россия, Новосибирск, kmaria@osmf.sscc.ru

APPLICATION OF KINETIC MODELS AND SIMULATION OF MULTI-PARTICLE SYSTEMS IN SOME FIELDS OF SCIENCE

Burmistrov A.V.1, Korotchenko M.A.2 institute of Computational Mathematics and Mathematical Geophysics SB RAS; Novosibirsk State University, Russia, Novosibirsk, e-mail: burm@osmf.sscc.ru 2Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Russia, Novosibirsk, e-mail: kmaria@osmf.sscc.ru

Аннотация. В статье рассмотрен и систематизирован опыт авторов, полученный в области моделирования многочастичных систем для ряда задач, возникающих в различных областях науки, таких как динамика разреженного газа, коагуляция частиц, транспортные потоки и ценообразование. Все эти задачи описываются уравнениями типа Больцмана. Для их численного решения применяется подход, использующий интегральное уравнение второго рода и соответствующую марковскую цепь, которая однозначно определяется коэффициентами интегрального уравнения. Это позволяет распространить хорошо развитую теорию весовых методов Монте-Карло на рассмотренные задачи.

Ключевые слова: метод Монте-Карло; эволюция многочастичной системы; парные взаимодействия; уравнение коагуляции; автотранспортный поток; ценообразование

Abstract. In this paper we cover our experience applied in the field of multi-particle systems simulation for some problems arising in various fields of science such as rarefied gas dynamics, particle coagulation, vehicular traffic flows, and price formation. All these problems are described by the Boltzmann type equations. To solve them numerically, we apply the approach which uses the integral equation of the second kind and the corresponding Markov chain. The latter is uniquely

determined by the coefficients of the integral equation. It allows spreading the well-developed theory of weight Monte Carlo methods to the considered problems.

Keywords: Monte Carlo method; multi-particle system evolution; pairwise interactions; coagulation equation; vehicular traffic flow; price formation

Introduction. A number of problems in mathematics, physics, economics etc. could be reduced to the estimation of some linear functionals of solutions to integral equations. This is mainly because mathematical models of such problems are constructed based on the corresponding stepwise Markov process, which terminates with probability one. The transition density of the latter process is the substochastic kernel of the integral operator, which describes the evolution of the system during a single step. The first example of such equation is the nonlinear kinetic Boltzmann equation, which remains the basis of the kinetic gas theory. This integro-differential equation describes the dynamics of a rarefied gas and was derived by Ludwig Boltzmann in 1872 (it could be found in the original text [1]). The nonlinearity of this equation lies in the collision integral, which describes the pairwise interactions of the particles. Though the kinetic equation was first obtained for the homogeneous relaxation of a single component ideal gas, the range of its applications turned out to be much wider. Boltzmann type equations are used to study the radiation transfer in matter, neutrons transfer in a nuclear reactor, electrons transfer in solids and plasmas, and also to study the growth of droplets in clouds, defects in materials, gas pores in metals, etc.

In this paper we cover our experience applied in the field of multi-particle systems simulation for some problems arising in various fields of science such as rarefied gas dynamics, particle coagulation, vehicular traffic flows, and price formation. All these problems are described by the Boltzmann type equations. To solve them numerically, we apply the approach which uses the integral equation of the second kind and the corresponding Markov chain. The latter is uniquely determined by the coefficients of the integral equation. It allows spreading the well-developed theory of weight Monte Carlo methods to the considered problems. Moreover, this makes it possible to estimate the parametric derivatives of the solution. This is especially important when one studies the influence of various parameters on the solution of nonlinear kinetic equations. For some problems we construct "value" modifications of the weight statistical modeling, which leads to a considerable reduction of computational costs when solving numerically the kinetic equation.

Rarefied Gas Dynamics. In the spatial homogeneous case, the Boltzmann equation for a single-component gas, as shown for example in [2], could be written using the symmetric probability density w = w(v', v\ | v, vj) of the velocities transition of two particles from values (v ', v ' 1) to new velocities (v, vj):

f (v,t) = J{f (v',t)f (v'j,t) - f (v,t)f (vj,t)}wdv'dv'j dvj.

Here f (v, t) is one-particle velocity distribution function at time t. Density w is the

product of the differential cross section for scattering of particles and two 5-functions, which take into account the physical conservation laws of momentum and energy in the case of pairwise particles interaction.

It is often assumed that the number N of interacting particles does not change, when using the probabilistic model of the multi-particle system evolution (see, for example, [3, 4]) for an approximate solution of the nonlinear Boltzmann equation by the Monte-Carlo method. In this case, according to [3], the model process of the stochastic kinetics of N particles system is a homogeneous Markov process. The phase states of the latter change as a result of elementary, namely pairwise, interactions of particles, since the number of interactions of a larger number of particles is negligible. The time between elementary interactions in a system is determined by the state of the entire system and has a generalized exponential distribution. However, since the kernel of the standard basic integral equation is the sum of mutually singular terms, the latter equation associated with the Markov process cannot be directly used to construct weight modifications of statistical modeling.

This obstacle was overcome in [5] with the help of expansion of the system phase space by introducing into the set of phase variables the number n = (i, j) of two particles that interact in the system. Such a procedure leads to a splitting of the distribution of interactions in the system by the number of a pair of interacting particles. This technique allowed in [5] to formulate an integral equation of the second kind of a special type, in which the density of collisions in the extended system serves as the desired quantity F:

Since the kernel of the latter equation contains singularities only as multipliers, it is convenient to construct standard weight modifications of the statistical modeling of the multiparticle system under consideration. The operator K could be considered as acting from Lx into Lx, despite the presence of generalized functions in its kernel, as it is indicated in [6]. Moreover, the Neumann series for this integral equation converges in the L norm.

It is well-known that under the "molecular chaos" assumption (see [3] for details) the normalized particle density in an N -particle system satisfies the Boltzmann equation asymptotically, when N ^ro. Moreover, the corresponding

error has order O( N_1), as a rule (see, for example, [4]).

For a given weight function h, we should estimate the linear functional (h, f)

of a solution f to the original Boltzmann equation. It is proved to be equal to a linear

functional (H, F) of a solution F to the integral equation of a special type and the

new weight function H~ is directly related to the function h (see [5] for more details). For the latter functional, weight modifications of statistical "collision" estimat ors as

F = KF + F.

0

well as "absorption" estimators are constructed in the standard way (see, for example,

In the spatially inhomogeneous case, the phase coordinates are the physical coordinates R = (rj,...,rN) and velocities V = (vj,...,vN) of all N particles as well as the number 7 = (i, j) of the interacting pair. Moreover, the state of the system is fixed immediately before the new free path modeling. In this case, the transition of the system from the current state (Z', t') to the next state (Z, t) is carried out as follows [4]:

1. choose the free path time r of the system according to the generalized exponential density, in this case we have t = t'+r;

2. calculate new coordinates of all particles R = R'+r • V';

3. simulate the number 7 of a pair of particles that implement the next pairwise interaction;

4. determine the velocities of all particles: for particles with numbers (i, j) the new velocities are chosen according to the given distribution density, while the velocities of the other particles are not changed.

Based on the described transition to an integral equation of the second kind, both correlated and weight algorithms were constructed to study the dependence of the solution of the Boltzmann kinetic equation on the time parameter and the initial one [7]. In addition, the weight algorithm of modeling the initial velocity distribution was developed for solving the Boltzmann equation (for more details, see [4]).

Coagulation Equation. The model proposed in the previous section could be adapted with some changes to estimate the solution to Smoluchowski coagulation equation. This kinetic equation describes a wide class of coagulation (or coalescence) processes in various physical systems consisting of particles with integer sizes (we are considering such a case so far, and that is why sums appear instead of integrals in the equation).

Suppose that for given coagulation coefficients Kkl, the probability of interaction (or collision) of particles with sizes k and l over a time interval At is Kw At. We call a particle of size k as a k -mer or k -meric particle. Then, in the

spatially homogeneous case, the concentration of k -mers nk (t) at the time instant t satisfies the following kinetic equation:

In this equation the rate of change in the concentration of k -mers over the time is represented as the sum of two terms: the first term is the rate at which k -mers originate as a result of the coupling of two smaller particles (the factor 1/2 is necessary here to take into account every such interaction only once); the second term describes a decrease in the concentration of k -mers due to their coupling with other particles.

[4]).

£ Klmnl (t)nm (t) - £ Kklnk (t)nl (t).

In the case of Smolukhovsky equation, the integral equation is related to the Kolmogorov equation, which is a probabilistic description of the evolution of N -particle system and gives an approximate solution of the corresponding nonlinear kinetic equation (see [8] for more details). Note that, in contrast to the Boltzmann equation, the number of particles in the system is variable since k -meric and l -meric particles interact to produce a single (k +1) -meric particle.

In the model used to estimate the solution of the coagulation equation, the phase coordinates are the current number N of particles in the system, the set of integer sizes for all these particles LN = (l1,...,lN), and the number ж = (i, j) of the interacting pair: Z = (N, LN ,ж). In this case, the transition of the system from the current state (Z', t') to the next state (Z, t) is carried out as follows (see, for example,

[9]):

1. choose the free path time t of the system according to the generalized exponential density, in this case we have t = t '+t ;

2. simulate the number ж = (i, j) of a pair of particles that implement the next pairwise interaction;

3. transform the system, i.e. replace two interacting particles i and j with a single particle whose size is l = lt + lj. Thus, as a result of the interaction, the

number of particles in the system decreases by one: N = N-1.

We investigated the Smoluchowski equation with linear coagulation coefficients depending on two parameters. For the numerical estimation of linear functionals of the solution to the equation under consideration, both the weight absorption estimator

[10] and collision estimator [11] have been developed. These algorithms allow estimating functionals for various coefficient parameters as well as its parametric derivatives with respect to these parameters using one set of trajectories. Moreover, value modifications of the algorithms for different transition steps are constructed: both for value modeling the free path time [7, 9], and for choosing the number of an interacting pair [10, 11, 12]. As a result we obtain a considerable reduction of computational costs for estimating the desired functionals. In addition, the problems with random parameters in coagulation coefficient are now of great interest [13].

Vehicular Traffic Flow (VTF). The authors succeeded to modify developed algorithms and apply them for solving a number of vehicular traffic problems. To solve them, a kinetic model of VTF proposed in [14], was used. A distinctive feature of the latter model is the introduction of acceleration into the set of phase coordinates. Now the state of the car is described by velocity and spatial coordinate, traditionally used in kinetic models, as well as acceleration. Such a transformation of the phase space allowed the kinetic model to be extended to sufficiently dense VTFs, i.e. to congested flows and synchronized flows. We study the relationship between the number of cars (which is analog of the gas concentration in gas dynamics) and the frequency of acceleration changes made by cars (which is analogous to the reaction speed), therefore the model is called kinetic. In the framework of the model under

consideration, the single-particle distribution density f = f (a, v, t) of cars with acceleration a and velocity v satisfies the Boltzmann-type integro-differential equation

д д

--ь a —

dt dv

f = J (E(a | d, v, a, v ) f (a], v, t)-Z(a' | a, v, a, v ) f}f (a, v, t)dadvdd,

here a and v are acceleration and velocity of the leader car, correspondingly. The leader here is the car straight in front of the "current" car; just between these two cars the interaction takes place (acceleration, braking or overtaking). The function !(•) here is a given weighted interaction density.

There are several differences between the VTF model and the gas dynamics. First one is the consequence of introduction of acceleration into the phase space in this model. So we consider the usage of accelerator and brake pedals by drivers as acts of kinetic pairwise interactions of cars, which lead to jumps in acceleration variables. It is unlike gas dynamics, where the pairwise interactions of molecules lead to jumps in velocity variable. Moreover, the leader does not change its acceleration after interaction of two cars, this explains the asymmetry of the function !(•). Therefore, interacting cars should be described by ordered pairs 7 = (i, j).

For definiteness, we can assume that the first index is the number of the follower, i.e. the car following the leading one, and the second index is the number of the leader. Note also that, in contrast to the Boltzmann equation, in the case of VTF the conservation laws of momentum and energy are not satisfied with pairwise interactions. It is also should be noted that in this model, as in many others VTF kinetic models, the "automobile chaos" assumption is used, which is similar to the "molecular chaos" hypothesis in gas dynamics. This assumption is further confirmed by the fact that changes in acceleration depend not only on the type of car, but also on the behavior and skills of each particular driver, which introduce additional randomness into the model.

In this model, the phase coordinates are the velocities V = (vj,...,vN) and accelerations A = (aj,...,aN) of all particles (cars), as well as the number 7 = (i, j) of the interacting pair [15]. In this case, the transition of the system from the current state (Z ,t ) to the next state (Z,t) is as follows:

1. choose the free path time r of the system according to the generalized exponential density, in this case t = t'+r;

2. calculate all the new velocities of cars V = V+rA';

3. simulate the number 7 = (i, j) of a pair of cars that interact;

4. determine the accelerations A of all cars: for the follower i , the new acceleration is selected according to the obtained distribution density, while the accelerations of the leader j and other vehicles do not change.

We considered the interaction profiles with threshold functions which depend on the velocities of the cars in the leader-follower pair [16]. Algorithms for estimating the distribution of velocity and acceleration in space homogeneous case were

developed. With the help of numerical experiments, the practicability and efficiency of using the integral equation and modeling the corresponding Markov chains in solving VTF problems were demonstrated (for more details, see [15]). Furthermore, the problems with random parameters in properties of cars and drivers are of our interest [17].

Price Formation. For the price formation model proposed in [18] we study the relationship between the number of orders on the exchange (analog of substance concentration) and the frequency of transactions (analog of reaction speed), so this model is also called kinetic. Like the interaction of molecules in gas dynamics, in this model, transactions are considered as acts of kinetic interactions of traders, namely buyers and vendors [19].

We supply all the functions associated with buyers with an auxiliary indicator n = 0, for vendors we put n = 1, hereafter n = 1 -n. In addition to discrete-time transactions between these two groups of traders, the price of an asset continuously changes due to random fluctuations, which are generally described by some stochastic differential equation [20], in the simplest case the Brownian motion with constant volatility g is used. The pricing process is modeled by splitting [21] into two stages: continuous (or dynamic) [22] and discrete (or kinetic) [23].

We denote s(x, t,n) = s the density of price distribution, here x is a bid price

for buyer and ask price for vendor. The transaction takes place with some probability at the time, when the vendor and the buyer reach agreement on the price. After the transaction, the buyer and the vendor change roles (buyer becomes the vendor and wants to sell the asset at a higher price, and vendor visa versa). Kinetic equation has the following form (see [18]):

Ся _2 r\2 \ О JO

s(x, t,ж) = J[К0ж • s(x', t,ж) ■ s(y', t,ж) - К1ж • s ■ s(y, t,ж)]dydx' dy',

ydt 2 dx2 j

The kernel K = Knn( x, y, x', y') shows the number of transactions per unit of time between buyers who wish to buy at a price x and resell after a transaction at a price x , and vendors who want to sell at a price y and after the transaction buy again at a price y . For the numerical solution of the equation, it is proposed to use the modeling of price changes in the market represented by a large number of buyers and vendors. In this case, the transactions are the result of the interaction of the buyer-vendor pair in the corresponding Markov chain.

In this model, the phase coordinate is the set of bid and ask prices X = (x1,... ,xN ) , as well as the number n = (i, j) of the pair involved in the transaction: Z = (X,n). In this case, the transition of the system from the current state (Z , t ) to the next state (Z, t) is as follows:

1. choose the free path time r of the system according to the generalized exponential density, in this case t = t'+r;

2. simulate the number n = (i, j) of a pair of traders involved in the transaction;

3. transformation takes place in the system, taking into account the transaction price and change of xi and Xj in accordance with a given density.

To estimate linear functionals of the solution of equations describing price dynamics on the stock exchange, in the framework of the kinetic model, weight Monte Carlo algorithms [23] were constructed. In addition, algorithms are being developed to study the parametric dependencies of the desired functionals (for example, on the frequency of transactions or transaction prices).

Conclusion. We consider four problems described by the kinetic Boltzmann type equations: rarefied gas dynamics, particle coagulation, vehicular traffic flow and price formation. To solve these problems, we propose to use the integral equation of the second kind approach and the weight statistical simulation of the corresponding Markov chain. For some problems we construct "value" modifications of the numerical algorithms, which lead to a considerable reduction of computational costs.

This work is conducted within the framework of the budget project 0315-20160002 for ICMMG SB RAS and is partly supported by the Russian Foundation for Basic Research (grants 17-01-00698, 18-01-00356, 18-01-00599).

Литература

1. Boltzmann L. Weitere Studien über das Wärmegleichgewicht unter Gasmolekülen. In: Brush S.G. Kinetische Theorie II: Irreversible Prozesse Einführung und Originaltexte, Wiesbaden: Vieweg+Teubner Verlag. 1970. P.115-225.

2. Ivanov M.S., Rogasinsky S.V. The direct statistical simulation method in dilute gas dynamics. Novosibirsk: CC SB USSR AS, 1988, 117pp. [in Russian]

3. Kac M. Probability and Related Topics in Physical Sciences. London ; New York: Interscience Publishers, 1959, 266pp.

4. Ivanov M.S., Korotchenko M.A., Mikhailov G.A., Rogasinsky S.V. New Monte Carlo Global Weight Method for the Approximate Solution of the Nonlinear Boltzmann Equation // Russian Journal of Numerical Analysis and Mathematical Modelling. 2004. V. 19. No.3. P.223-238.

5. Mikhailov G.A., Rogazinskii S.V. Weighted Monte Carlo Methods For Approximate Solution Of A Nonlinear Boltzmann Equation // Siberian Mathematical Journal. 2002. V.43. No.3. P.496-503.

6. Mikhailov G.A. New Monte Carlo methods with estimating derivatives. Utrecht: VSP, 1995, 186pp.

7. Mikhailov G.A., Korotchenko M.A., Rogasinsky S.V. Importance Modeling Algorithms for Solving Nonlinear Kinetic Equations // Doklady Mathematics. 2007. V.76. No.1. P.502-505.

8. Lushnikov A. Some new aspects of coagulation theory. Izvestiya of the Academy of Sciences of the USSR. Atmospheric and Oceanic Physics. 1978. V.14. No.10. P.1046-1054. [In Russian]

9. Korotchenko M.A. Statistical algorithms of value modeling for solving the Smoluchowski equation // Computational Technologies. 2008. V.13. Special issue 4. P.68-74. [In Russian]

10.Burmistrov A.V., Korotchenko M.A. Weight Monte Carlo Algorithms for Estimation and Parametric Analysis of the Solution to the Kinetic Coagulation Equation // Numerical Analysis and Applications. 2014. V.7. No.2. P. 104-116.

11.Burmistrov A.V., Korotchenko M.A. Usage of collision estimator for solving the Smoluchowski equation by Monte Carlo method // Marchuk Scientific Readings-2017, International Conference Proceedings. Novosibirsk: ICMMG SB RAS. 2017. P.126-132. [In Russian]

12.Korotchenko M. Value Monte Carlo Algorithms for Estimating the Solution to the Coagulation Equation // Springer Proceedings in Mathematics and Statistics. 2012. V.23. P.511-522.

13.Burmistrov A.V., Korotchenko M.A. Algorithms of statistical modeling for solving the kinetic coagulation equation with random parameters // Abstracts of International Conference "Computational Mathematics and Mathematical Geophysics". Novosibirsk: Akademizdat, 2018. P.26-27. [In Russian]

14.Waldeer K.T. A Vehicular Traffic Flow Model Based on a Stochastic Acceleration Process // Transport Theory and Statistical Physics. 2004. V.33. No.1. P.7-30.

15.Burmistrov A., Korotchenko M. Weight Monte Carlo Method Applied To Acceleration Oriented Traffic Flow Model // Springer Proceedings in Mathematics and Statistics. 2012. V.23. P.297-311.

16.Burmistrov A., Korotchenko M. Monte Carlo Algorithm for Simulation of the Vehicular Traffic Flow Within the Kinetic Model with Velocity Dependent Thresholds // Springer Proceedings in Mathematics and Statistics. 2014. Vol.114. P.109-117.

17.Burmistrov A.V., Korotchenko M.A. Simulation of vehicular traffic flow within the kinetic model with random parameters // Abstracts of XIXth all-russian conference of young scientists on mathematical modeling and information technologies. Novosibirsk: ICT SB RAS. 2018. P.12-13. [In Russian]

18.Burger M., Caffarelli L., Markowich P., Wolfram M.-T. On a Boltzmann-type price formation model // Proceedings of The Royal Society A. 2013. V.469. No.2157. 20130126.

19.Burmistrov A.V., Novikov A.V. Stochastic kinetic model of price formation // International Research Journal. 2017. No.7-3(61). P.107-112. [In Russian]

20.Novikov A.V., Burmistrov A.V. Improvement of the stochastic differential model of the price series // Marchuk Scientific Readings-2017, International Conference Proceedings. Novosibirsk: ICMMG SB RAS. 2017. P.654-658. [In Russian]

21.Novikov A.V., Burmistrov A.V. Usage of trading algorithms in adaptive intellectual ecosystem // Economy and business: theory and practice. 2018. No.4. P.184-198. [In Russian]

22.Burmistrov A.V., Novikov A.V. Definition of the coefficients in the stochastic differential model of price formation // Modern Science: actual problems of theory and practice. Series: Natural and Technical Sciences. 2018. No.3. P.39-44. [In Russian]

23.Korotchenko M.A. Construction of statistical estimates by modeling the price dynamics in the framework of the kinetic price model // In: Mathematical and computer modeling of natural scientific and social problems. 2018. P.174-179. [In Russian]