ДОКЛАДЫ АКАДЕМИИ НАУК, 2007, том 415, № 1, с. 26-30
МАТЕМАТИКА
УДК 519
ЦЕННОСТНЫЕ АЛГОРИТМЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ ДЛЯ РЕШЕНИЯ НЕЛИНЕЙНЫХ КИНЕТИЧЕСКИХ УРАВНЕНИЙ
© 2007 г. Член-корреспондент РАН Г. А. Михайлов, М. А. Коротченко, С. В. Рогазинский
Поступило 31.01.2007 г.
ВВОДНАЯ ИНФОРМАЦИЯ
В настоящей работе рассматриваются "ценностные" весовые модификации статистического моделирования для численного решения нелинейных кинетических уравнений Смолуховского и Больцмана в пространственно-однородном случае.
Уравнение Смолуховского описывает широкий класс процессов коагуляции. Для определенности мы будем рассматривать случай парного слипания частиц. Этот процесс является причиной нелинейности уравнения коагуляции
М( 0 =
2
' +1 = I ' г ]
д t
2 х k'Jn' (t) ni (t) - х кпп> (t) ni (t),
(i)
П (0) = n(0),
где n(t) - частота (числовая плотность) частиц размера I в момент времени t, причем размер частицы I принимает натуральные значения; ку - коэффициенты коагуляции.
Уравнение Больцмана, описывающее процесс однородной релаксации простого однокомпо-нентного газа, имеет вид
v, t) = J|v - vi { f(v', t)f(vi, t) -
-f(v, t)f(vi, t)}o(|v- vi, X)sinX dx de dvi.
Здесь f(v, t) - "одночастичная" плотность распределения по скоростям v в момент времени t е (0, T];
скорости (v', vi) и (v, v1) удовлетворяют законам сохранения импульса и энергии при столкновении; c(|v - v1|, x) - дифференциальное сечение рассеяния. Присоединяя к уравнению начальные данные f(v, t = 0) = fo(v), получим задачу Коши для нелинейного уравнения Больцмана.
Институт вычислительной математики и математической геофизики Сибирского отделения Российской Академии наук, Новосибирск
Известно, что решение кинетического уравнения можно оценивать, статистически моделируя марковский процесс эволюции соответствующего N-частичного ансамбля, фазовые состояния которого меняются вследствие парных взаимодействий частиц (см., например, [1, 2]). Связанное с этим процессом базовое интегральное N-ча-стичное уравнение Каца [1] использовать непосредственно для построения весовых модификаций статистического моделирования затруднительно, так как его ядро представляет собой сумму взаимосингулярных слагаемых. Это затруднение преодолено в [2] посредством модификации фазового пространства системы путем введения номера взаимодействующей пары в число координат фазового пространства системы. Такой важный исходный момент приводит к расслоению распределения взаимодействий в системе по номеру пары взаимодействующих частиц п = (i, j). Данный прием позволил [2] сформулировать новое интегральное уравнение, на основе которого удобно строить стандартные весовые модификации статистического моделирования рассматриваемой многочастичной системы, поскольку ядро уравнения содержит сингулярности лишь в виде сомножителей. Это интегральное уравнение имеет вид
F(Z, t) =
t
= JJF(Z\ t')K(Z, t ^ Z, t)dZdt + F0(Z)5(t), (2)
0 Z
где Z = (X, п), dZ = dXd^0(n), причем интегрирование по мере ц0 означает суммирование по всем различным парам п = (i, j).
Для уравнения Смолуховского [3] X = (N, l1, l2, ..., lN), N < N0, I, е 1, 2, ..., N0 - размер i-й частицы, а для уравнения БольцманаX = (vj, v2, ..., vN), vi - скорость i-й частицы. Получаемая из фазовой плотности F(% t) одночастичная плотность является приближенным решением соответствующего кинетического уравнения в предположении молекулярного хаоса [1].
Для целей настоящей работы целесообразно фиксировать состояние X Л-частичного ансамбля в начале свободного временного пробега. Ядро соответствующего уравнения (2) имеет вид
К( У, ? — 2, г) - к (г'- ЦХ) ^^ К (X1 — Х| п),
где
A (X)
k(t ^ t\X) = A(X1)E(X, t - t),
причем Е(Х, г) = ехр{-г^(Х)}. Условная переходная плотность К1(Х' — Х|п) определяет преобразование ансамбля в результате взаимодействия пары с номером п = (/, у). Для уравнения Смолухов-ского это замена взаимодействующих частиц на одну частицу размера ^ + I, причем N=N _ 1, а для
уравнения Больцмана - замена (V , у,) — (уг-, у,) соответственно заданной индикатрисе.
Обычно при решении уравнения Каца вычисляют функционалы JH(t) от потока "частиц" Ф(Х, г) в заданные моменты времени. На основе теории локальных оценок [4, раздел 6.9] получается равенство
цепи Маркова и вычисление значения % происходят по следующей схеме.
1. Для г0 = 0 моделируется начальное состояние 20 соответственно плотности Р0(2Т) и вычисляется Q0.
2. Согласно плотности Р(2п _ 1, гп _ 1 — 2п, гп) выбирается гп и моделируется 2п, причем вес преобразуется после каждого элементарного перехода.
3. Если гп < Т, то вычисляется слагаемое
QnH (Хп, Т _ гп), иначе цепь обрывается.
Заметим, что интегральное уравнение (2) может быть получено непосредственно из этой марковской модели как соотношение интегрального баланса для суммарной плотности состояний.
ЦЕННОСТНОЕ МОДЕЛИРОВАНИЕ ДЛЯ КИНЕТИЧЕСКИХ УРАВНЕНИЙ
Согласно терминологии, принятой в теории весовых методов Монте-Карло, функцией ценности ф*(Х, г) называется значение функционала JH для источника с плотностью 5(Х' _Х)5(г' _ г), т.е.
г(г) = |Н(Х)Ф(X, г)ёХ -
г
- Л Н (X, г - Г) О (2, ?) М,
0 2
где Н(X) е Н (X, г) = Н^Е^, г).
Рассмотрим цепь Маркова (2п, гп), п е 0, 1, ..., с нормированной плотностью перехода
Р(г — 2, г) - р(г — )Р1(п|X1)Р2(X —XIп)
и нормированной плотностью начального состояния Р0(2, г) = Р0(2)5(0. Случайные веса вводятся по формулам
Р (2)
Q0 - Р (у)' Qn - Qn-1Q(2п-1, гп-1. 2п, гп)? (Л ( 7 г . у г ) - К ( 7 - — г)
Q(2п-1'1п-1; ^1п) Р(2' г — 2, 0'
Величина Q очевидным образом факторизуется. Для оценки функционала JH(t) можно использовать весовой вариант оценки по столкновениям (в начале пробегов):
V
% - X QnH(Xn, Т - гп),
V = max{n: tn < T,n = 0, 1, ..., N0 - 1},
где T - момент времени, в который оценивается рассматриваемый функционал. Моделирование
ф* (X, t) = E£,
( x, t),
Lx, t) = H(X, T - t)■
+
£ QnH(Xn, T - tn).
п - 1
Функция ф*(X, г) является решением интегрального уравнения, сопряженного к уравнению (2), и тем самым в предположении молекулярного хаоса может быть приближена решением сопряженного нелинейного кинетического уравнения [5]. Однако, как показано ниже, иногда практически удовлетворительное приближение фг) может быть получено непосредственно на основе оценки функционала JH для модельного кинетического уравнения.
С учетом вида переходной плотности Р для моделирования перехода (У, О — (2, г) целесообразно реализовать три элементарных перехода. Согласно [4, раздел 4.5], моделирование называется ценностным, если каждый элементарный переход реализуется соответствующей нормированной плотностью, полученной путем умножения физической плотности на соответствующую вспомогательную функцию ценности, которая равна ф*(А", г) лишь для финального элементарного перехода. Для первого и второго переходов в рассматриваемой задаче необходимо использовать следующие ценностные множители:
ф*(г; X1) - Ла(п)X)К1 (X — XIп)х
Хф* (X, г)dX ф(п), Ф*(п; X1, г) - |К1 (X — XIп)ф*(X, г)dX.
V
i = 0
28
МИХАИЛОВ и др.
Следующее утверждение является следствием общей теоремы о ценностном моделировании [6] всех элементарных переходов.
Теорема 1. Если реализуются нормированные плотности элементарных переходов вида
С^(Г ^ ¿IX)ф* (г; X), C1AXФ*(п; X, г),
C2К1 (X ^ XIП)Ф*(X, г),
то = 0.
Теорема 1 показывает, что если используются достаточно точные приближенные функции ф* , ф* и ф* , то величина может быть малой.
ЦЕННОСТНОЕ МОДЕЛИРОВАНИЕ ДЛЯ РЕШЕНИЯ ТЕСТОВОЙ ЗАДАЧИ КОАГУЛЯЦИИ
Далее будет рассмотрено ценностное моделирование для оценки решения тестовой задачи (1) при kij = 1 и п1(0) = 1, п(0) = 0, I > 2. Это решение известно [7]:
п (г) =
1
(1 + 0.5 г )2 V 1 + 0.5 х
0.5 г
I-1
I > 1.
0.5 N1N -1)
Для такой задачи а(п) = N0 , Л^ =-----
N0
Решается задача оценки величины (Т) ■■ ■■ п1(Т), т.е. среднего числа (частоты) мономеров
момент времени Т, при этом Н(X) = H1(X) =
N1 ( X) N '
ф1 (г) =
Ф* (г) = 1Е( г), 1
(1 + 0.5(Т- г))2 (2 + Т + е)2
х
х ехр
2 + Т + е
г • Iе(г),
Т2 = (Л(X) - Л0)-11п (1- а^), ^ = 1 - ехр(-(Л(X) - Л0)(Т + е- г')).
После этого вес домножается соответственно на один из множителей
41 = 42 = -
Л (х)
Л(X) - Л0
(л (x) - л0)х
ел (X) л (X) - Л0)(Т + е)
- е
Соответственно этому (с учетом теорем 1, 2 и нормированное™ переходных плотностей) рассматриваются два приближения к временной функции ценности:
Далее рассматривается ценностное моделирование номера пары. Обозначим: N - общее количество частиц, N - количество мономеров в ансамбле в конце свободного пробега (перед выбором значения п).
Все множество возможных пар для столкновения разобьем на три взаимоисключающих непересекающихся типа:
1) пары вида {мономер, мономер}, число пар такого типа равно М1 = N1 (N - 1)/2;
2) пары вида {мономер, полимер}, число пар такого типа равно Л2 = N1N - N1);
3) пары вида {полимер, полимер}, число пар такого типа равно = 0.5(А - N1 )(К - N1 - 1).
Представим физическое равномерное распре-2
деление Р0(', 1) = —р- номера взаимодействующей пары в следующем рандомизированном виде:
X р(г, 1) = р X fl(1)+
1 < i < 1 < N
(', 1) - пара типа 1
+ Р2 X ^2 ( 1) + Рз X ( ' ' 1),
(3)
(г', 1) - пара типа 2
1
(1) - пара типа 3
где !е(0 - индикатор отрезка [0, Т + е], е - величина продолжения временного промежутка, на котором строится цепь Маркова. В [3] показано, что
|/Н) (Т) - П1(Т)| = 0(N01). С учетом этих приближений свободный пробег т = г - г' моделируется по формулам
т1 = -Л_1(X1) 1п (1- а51), 51 = 1 - ехр(-Л(X1)(Т + е - г')),
где /к(г, 1) = —— равномерная вероятность выбо-
^ k
ра пары (', 1) из пар типа k, а рк - вероятность выбора ^го типа, причем
Р1
N1 (N1-1)
Р2
2 N1 (N - N1)
N(N-1)' N(N-1)
(N - N1)(N - N - 1)
Р3
N (N - 1)
В целях "сохранения" мономеров моделирование пары будем производить соответственно представлению (3), в котором заменим вероятно-
сти pk на величины qk, пропорциональные числу оставшихся мономеров:
qi
qs
p 1 ( N1 - 2)
C '
Рз N1 C :
q2
с =
= p 2 ( N1 - 1 ) C '
N1 (N-2)
N •
(4)
Такая модификация учитывается при построении
веса
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.