КИНЕТИЧЕСКИЕ МОДЕЛИ И МОДЕЛИРОВАНИЕ МНОГОЧАСТИЧНЫХ СИСТЕМ

KINETIC MODELS AND MULTIPARTICLE SYSTEMS SIMULATION

Aleksandr Burmistrov

candidate of Science, research associate at Institute of Computational Mathematics and Mathematical Geophysics SB RAS, senior teacher at Novosibirsk State University,

Russia, Novosibirsk

Mariya Korotchenko

candidate of Science, research associate at Institute of Computational Mathematics and Mathematical Geophysics SB RAS,

Russia, Novosibirsk

АННОТАЦИЯ

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

ABSTRACT

We consider four problems, three of which are not related to rarefied gas dynamics, but are described by the Boltzmann type equations. To solve these problems, we propose to use the integral equation of the second kind approach and the weight simulation of Markov chain, the latter is uniquely determined by the coefficients of the integral equations. In the paper, authors summarize the developed approaches in the field of multiparticle systems simulation for different problems: rarefied gas dynamics, particle coagulation, vehicular traffic flow and price formation.

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

Keywords: Monte Carlo method; pair interactions; vehicular traffic flow; coagulation equation; price formation.

Введение.

Решение целого ряда задач математической физики и индустриальной математики можно свести к оценке линейных функционалов от решений интегральных уравнений. Это связано, прежде всего, с тем, что математическая модель таких задач строится на основе рассмотрения соответствующего скачкообразного и обрывающегося с вероятностью единица марковского процесса, чья переходная плотность является субстахостическим ядром интегрального оператора, описывающего динамику системы за один шаг. Одним из примеров таких задач является нелинейное кинетическое уравнение Больцмана в динамике разреженного газа. Это нелинейное интегро-дифференциальное уравнение описывает поведение разреженного газа и было выведено Людвигом Больцманом в 1872 году, однако до сих пор остается основой кинетической теории газов. Нелинейность уравнения Больцмана заключается в интеграле столкновений – квадратичном операторе, который описывает парные столкновения (или взаимодействия) частиц.

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

Применение алгоритма для решения задач динамика разреженного газа.

В пространственно-однородном случае уравнение Больцмана для однокомпонентного газа, как было показано в [2], может быть записано с использованием симметричной плотности  вероятности перехода пары скоростей частиц от () к новым скоростям ():

,

здесь  – одночастичная функция плотности распределения по скоростям в момент времени . Плотность  является произведением дифференциального сечения рассеяния частиц и d-функций, учитывающих физические законы сохранения (импульса и энергии) при парном взаимодействии частиц.

При использовании вероятностной модели эволюции многочастичной системы (см., например, [3; 4]) для приближенного решения нелинейного уравнения Больцмана методом Монте-Карло часто предполагается, что количество взаимодействующих частиц  не изменяется. При этом, согласно [3], модельный процесс стохастической кинетики системы из  частиц является однородным марковским процессом, фазовые состояния которого меняются в результате элементарных. а именно парных, взаимодействий частиц, поскольку число взаимодействий большего количества частиц пренебрежимо мало. Время между элементарными взаимодействиями в системе определяется состоянием всей системы и имеет обобщенное экспоненциальное распределение. Однако, связанные с этим процессом стандартные базовые интегральные уравнения невозможно непосредственно использовать для построения весовых модификаций статистического моделирования, так как ядра этих уравнений представляют собой сумму взаимосингулярных слагаемых. Это затруднение было преодолено в [6] посредством расширением фазового пространства системы путем введения в число фазовых переменных номера пары, реализующей взаимодействие в системе. Такая процедура приводит к «расслоению» распределения взаимодействий в системе по номеру пары взаимодействующих частиц . Данный прием позволил в [6] сформулировать интегральное уравнение второго рода специального вида, в котором в качестве искомой величины выступает плотность столкновений в расширенной системе :

.

На основе этого уравнения удобно строить стандартные весовые модификации статистического моделирования рассматриваемой многочастичной системы, поскольку ядро уравнения содержит сингулярности лишь в виде сомножителей. Оператор , несмотря на наличие в ядре обобщенных функций, можно, как указано в [9], рассматривать действующим из  в , причем ряд Неймана для интегрального уравнения сходится по норме . Известно, что в предположении «молекулярного хаоса» (подробнее см. [3]) нормированная плотность частиц в -частичной системе асимптотически, при , удовлетворяет уравнению Больцмана, причем как правило соответствующая погрешность имеет порядок  (см., например, [4]).

Линейный функционал  от решения  исходного уравнения Больцмана ( – заданная весовая функция) оказывается равен линейному функционалу  от решения  интегрального уравнения специального вида (– напрямую связанная с  функция) (детали см. в [6]). Для последнего функционала стандартным образом строятся весовые модификации статистических оценок «по столкновениям» или «по поглощениям» (см., например, [4]). В пространственно-неоднородном случае фазовыми координатами являются физические координаты и скорости всех частиц  и , а также номер взаимодействующей пары . Причем состояние  системы фиксируется непосредственно перед началом нового свободного пробега. При этом переход системы из состояния  в состояние  осуществляется следующим образом:

1.  выбирается время  свободного пробега системы согласно обобщенной экспоненциальной плотности, при этом ;

2.  вычисляются новые координаты всех частиц ;

3.  моделируется номер  пары частиц, реализующей очередное взаимодействие;

4.  определяются скорости  всех частиц: для частиц с номерами  новые скорости выбираются согласно полученной плотности распределения, а скорости остальных частиц не изменяются.

На основе описанного перехода к интегральному уравнению второго рода были построены алгоритмы коррелированного и весового моделирования для изучения зависимости решения кинетического уравнения Больцмана от начального и временнóго параметров. Кроме того, был разработан алгоритм весового моделирования начального распределения скоростей для решения уравнения Больцмана (подробнее о результатах см. в [4]).

Использование алгоритма для решения уравнения коагуляции.

Предложенную в предыдущем разделе модель можно с некоторыми изменениями адаптировать для оценки решения кинетического уравнения коагуляции Смолуховского. Это уравнение описывает широкий класс процессов коагуляции (слипания) в различных физических системах, состоящих из частиц с целочисленными размерами (мы рассматриваем пока такой случай, и именно поэтому вместо интегралов в уравнении появляются суммы).

Предположим, что при заданных коэффициентах коагуляции  вероятность взаимодействия (или столкновения) частиц с размерами  и  в течение временного интервала  равна . Будем называть -мером частицу размера . Тогда, в пространственно однородном случае, концентрация -меров  в момент времени  удовлетворяет следующему кинетическому уравнению:

.

С помощью этого уравнения скорость изменения концентрации -меров по времени представляется в виде суммы двух слагаемых: первое слагаемое – это скорость, с которой -меры появляются в результате слипания частиц меньшего размера (множитель 1/2 необходим для того, чтобы каждое такое взаимодействие было учтено только один раз); второе слагаемое описывает уменьшение концентрации -меров за счет их слипания с другими частицами.

В случае уравнения Смолуховского интегральное уравнение связано с уравнением Колмогорова, которое представляет собой вероятностное описание эволюции системы из  частиц и даёт приближённое решение соответствующего нелинейного кинетического уравнения (подробнее см. в [5]). Заметим, что, в отличие от уравнения Больцмана, количество частиц в системе является переменным.

В модели, используемой для оценки решения уравнения коагуляции, фазовыми координатами являются текущее количество частиц  в системе, набор целочисленных размеров этих частиц , а также номер взаимодействующей пары : . При этом переход системы из состояния  в состояние  осуществляется следующим образом:

1.  выбирается время  свободного пробега системы согласно обобщенной экспоненциальной плотности, при этом ;

2.  моделируется номер  пары частиц, реализующей очередное взаимодействие;

3.  происходит преобразование в системе, которое состоит в замене взаимодействующих частиц  и  одной частицей с размером . Таким образом, в результате взаимодействия число частиц в системе уменьшается, т. е. .

Авторы исследовали уравнение Смолуховского с линейными коэффициентами коагуляции, зависящими от двух параметров. Для численной оценки линейных функционалов от решения рассматриваемого уравнения разработаны весовые алгоритмы, позволяющие на одном наборе траекторий оценивать функционалы для различных параметров, а также параметрические производные по этим параметрам. Кроме того, построены ценностные модификации алгоритма выбора номера взаимодействующей пары частиц в процессе коагуляции, что позволило существенно уменьшить трудоемкости оценок искомых функционалов (более подробно см. [1]).

Обобщение алгоритма на автотранспортные задачи.

Авторам удалось применить с некоторыми модификациями разработанные алгоритмы на решение ряда автотранспортных задач. Для решения этих задач авторы в своей работе использовали кинетическую модель автотранспортного потока (АТП), предложенную в [10]. Отличительной особенностью выбранной модели является внесение ускорения в число фазовых координат, описывающих состояние автомобиля, наряду со скоростью и пространственной координатой, традиционно используемыми в кинетических моделях. Такое преобразование фазового пространства позволило распространить кинетическую модель на достаточно плотные АТП, т. е. на потоки с затрудненным или синхронизированным движением. В рамках рассматриваемой модели, одночастичная плотность  распределения автомобилей с ускорением  и скоростью  удовлетворяет интегро-дифференциальному уравнению Больцмановского типа

здесь  и  – ускорение и скорость лидера, соответственно. Под лидером здесь понимаем автомобиль, находящийся непосредственно перед «текущим»; между этими двумя автомобилями и происходит взаимодействие (ускорение, торможение или обгон). Функция  – заданная взвешенная плотность взаимодействий.

Следует отметить несколько отличий модели от газовой динамики. Прежде всего, в результате введения ускорения в число фазовых координат в данной модели при парных взаимодействиях скачкообразно меняется не скорость, а ускорение автомобиля. Кроме того, при взаимодействии двух автомобилей лидер не меняет своего ускорения, что выражено в несимметричности функции . Поэтому взаимодействующие автомобили нужно представлять упорядоченными парами . Для определенности можно полагать, что первый индекс  – номер преследователя, то есть автомобиля, следующего за лидирующим, а второй индекс  – номер лидера. Заметим также, что, в отличие от уравнения Больцмана, в случае АТП при парных взаимодействиях законы сохранения момента и энергии не выполняются. Также стоит отметить, что в данной модели, как и во многих кинетических моделях АТП, используется предположение об «автомобильном хаосе», которое является аналогичным гипотезе о «молекулярном хаосе» в газовой динамике. Данное предположение дополнительно подтверждается тем, что изменения ускорения зависят не только от типа автомобиля, но также от поведения и навыков каждого конкретного водителя, что вносит дополнительную случайность в модель.

В данной модели фазовыми координатами являются скорости и ускорения всех частиц (автомобилей)  и , а также номер взаимодействующей пары : . При этом переход системы из состояния  в состояние  осуществляется следующим образом:

1.  выбирается время  свободного пробега системы согласно обобщенной экспоненциальной плотности, при этом ;

2.  вычисляются новые скорости всех автомобилей ;

3.  моделируется номер  пары автомобилей, реализующей очередное взаимодействие;

4.  определяются ускорения  всех автомобилей: у преследователя  новое ускорение выбирается согласно полученной плотности распределения, а ускорения лидера  и остальных автомобилей не изменяются.

Авторами были рассмотрены профили взаимодействия с пороговыми функциями, зависящими от скоростей автомобилей в паре лидер-преследователь. Разработаны алгоритмы оценки распределения скорости и ускорения в однородном по пространству случае. С помощью численных экспериментов продемонстрирована практическая целесообразность и эффективность использования интегрального уравнения и моделирования соответствующих цепей Маркова при решении автотранспортных задач (подробнее см. в [8]).

Моделирование цены на рынке.

В заключении данной статьи мы рассмотрим систему двух кинетически уравнений типа Больцмана, которая возникает в соответствующей модели формирования цены, предложенной в [7]. В данной модели, кроме наличия двух типов частиц (покупатель и продавец), есть еще одно отличие от предыдущих уравнений: присутствуют случайные флуктуации цены, которые описываются второй производной и могут моделироваться с помощью расщепления винеровским процессом. Типы частиц характеризуются плотностями  и . Здесь  – цена спроса или предложения для покупателей и продавцов соответственно. Взаимодействие (сделка) между покупателем в состоянии  и продавцом в состоянии  происходит с некоторой вероятностью, зависящей от их состояний. После сделки покупатель и продавец меняются ролями (купивший становиться продавцом и желает продать актив по большей цене, а продавший – наоборот). Система уравнений имеет следующий вид (см. [7]):

Ядро  показывает количество сделок за единицу времени между покупателями, желающими купить по цене  и после сделки перепродать по цене , и продавцами, желающими продать по цене  и после сделки снова купить по цене . Для численного решения системы уравнений предлагается использовать моделирование изменений цен на рынке, представленном большим количеством покупателей и продавцов. При этом сделки являются результатом взаимодействия пары покупатель-продавец в соответствующей цепи Маркова.

В данной модели фазовой координатой является набор цен спроса и цен предложения , а также номер пары , участвующей в сделке: . При этом переход системы из состояния  в состояние  осуществляется следующим образом:

1.  выбирается время  свободного пробега системы согласно обобщенной экспоненциальной плотности, при этом ;

2.  моделируется номер  пары, участвующей в сделке;

3.  происходит преобразование в системе, учитывающее цену сделки и изменение  и  в соответствии с заданной плотностью.

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

Список литературы:

1.  Бурмистров А.В., Коротченко М.А. Весовые алгоритмы метода Монте-Карло для оценки и параметрического анализа решения кинетического уравнения коагуляции // Сиб. журн. вычисл. математики. 2014. Т. 17. № 2. С. 125–138.

2.  Иванов М.С., Рогазинский С.В. Метод прямого статистического моделирования в динамике разреженного газа. Новосибирск: ВЦ СО АН СССР, 1988.

3.  Кац М. Вероятность и смежные вопросы в физике. – М.: Мир, 1965.

4.  Коротченко М.А., Весовые параметрические алгоритмы статистического моделирования для решения нелинейных кинетических уравнений: дис. … канд. физ-мат. наук. – Новосибирск, 2008.

5.  Лушников А.А. Некоторые новые аспекты теории коагуляции // Изв. АН СССР, Физ. атмосферы и океана. 1978. Т. 14. № 10. С. 1048–1055.

6.  Михайлов Г.А., Рогазинский С.В. Весовые методы Монте–Карло для приближенного решения нелинейного уравнения Больцмана // Сиб. мат. журн. 2002. T. 43. № 3. С. 620–628.

7.  Burger M., Caffarelli L., Markowich P., Wolfram M.-T. On a Boltzmann-type price formation model // Proc. R. Soc. A 2013. V. 469. № 2157. 20130126.

8.  Burmistrov A., Korotchenko M. Monte Carlo Algorithm for Simulation of the Vehicular Traffic Flow Within the Kinetic Model with Velocity Dependent Thresholds // Springer Proc. Math. Stat. 2014. Vol. 114. P. 109–117.

9.  Mikhailov G.A. Parametric estimates by the Monte Carlo method. Utrecht: VSP, 1999.

10.  Waldeer K.T. The direct simulation Monte Carlo method applied to a Boltzmann-like vehicular traffic flow model // Comput. Phys. Commun. 2003. V. 156. № 1. P. 1–12.