ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 12, с. 2110-2121
УДК 519.676
МОДИФИКАЦИИ ВЕСОВЫХ АЛГОРИТМОВ МЕТОДА МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ НЕЛИНЕЙНЫХ
КИНЕТИЧЕСКИХ УРАВНЕНИЙ1)
© 2007 г. М. А. Коротченко, Г. А. Михайлов, С. В. Рогазинский
(630090 Новосибирск, пр-т Акад. Лаврентьева, 6, ИВМиМГ СО РАН) e-mail: kmary@rol.ru; gam@osmf.sscc.ru; svr@osmf.sscc.ru Поступила в редакцию 22.05.2007 г.
На примере решения тестовых задач для нелинейных кинетических уравнений Больцмана и Смолуховского исследована эффективность различных вариантов весовых "ценностного" моделирования эволюции многочастичных ансамблей. Существенный выигрыш в трудоемкости при решении задач коагуляции достигается путем приближенного ценностного моделирования длины "свободного пробега" ансамбля наряду с ценностным моделированием номера пары взаимодействующих частиц. Для модельного решения уравнения Больцмана наиболее эффективной оказалась весовая модификация моделирования начального распределения скоростей. Разработанная методика может быть полезной при решении реальных задач коагуляции и релаксации системы частиц, для которых рассмотренные модельные решения являются приближенными. Библ. 10. Табл. 3.
Ключевые слова: весовые алгоритмы метода Монте-Карло, решение нелинейных кинетических уравнений Больцмана, уравнение Смолуховского.
1. ВВОДНАЯ ИНФОРМАЦИЯ
В настоящей работе рассматриваются ценностные весовые модификации статистического моделирования для численного решения нелинейных кинетических уравнений Смолуховского и Больцмана в пространственно-однородном случае.
Уравнение Смолуховского описывает широкий класс процессов коагуляции. Для определенности будем рассматривать случай парного слипания частиц. Этот процесс является причиной нелинейности уравнения коагуляции
^Пр 2 ki]n1{t)n]{ t) knn(t)n (1.1)
где n¡(t) - частота (числовая плотность) частиц размера ¡ в момент времени t, причем размер частицы ¡ принимает натуральные значения; k¡j - коэффициенты коагуляции.
Уравнение Больцмана, описывающее процесс однородной релаксации простого однокомпо-нентного газа, имеет вид
= J| |{ - ( i, t)}o(|v- v',х)sinxdXdedvj. (1.2)
Здесьf(v, t) - "одночастичная" плотность распределения по скоростям v в момент времени t е (0, T]; скорости (v', v') и (v, Vj) удовлетворяют законам сохранения импульса и энергии при столкновении: v + v1 = v' + v', v2 + v2 = v'2 + v'2; c(|v - vx|, %) - дифференциальное сечение рассеяния. Присоединяя к уравнению начальные данные f (v, t = 0) = f0(v), получаем задачу Коши для нелинейного уравнения Больцмана.
Известно, что решение кинетического уравнения можно оценивать, статистически моделируя марковский процесс эволюции соответствующего ^-частичного ансамбля, фазовые состояния
j) Работа выполнена при финансовой поддержке РФФИ (код проекта 06-01-00046), Президентской программы "Ведущие научный школы" (грант НШ-4774.2006.1) и "Фонда содействия отечественной науке".
2110
которого меняются вследствие парных взаимодействий (см., например, [1], [2]). Связанное с этим процессом базовое интегральное А-частичное уравнение Каца (см. [1]) использовать непосредственно для построения весовых модификаций статистического моделирования затруднительно, так как его ядро представляет собой сумму взаимно сингулярных слагаемых. Это затруднение было преодолено в [2] посредством модификации фазового пространства системы путем введения номера взаимодействующей пары в число координат фазового пространства системы. Такой важный исходный момент приводит к "расслоению" распределения взаимодействий в системе по номеру пары взаимодействующих частиц п = (7, у). Данный прием позволил в [2] сформулировать новое интегральное уравнение, на основе которого удобно строить стандартные весовые модификации статистического моделирования рассматриваемой многочастичной системы, поскольку ядро уравнения содержит сингулярности лишь в виде сомножителей. Это интегральное уравнение имеет вид
í
Т, г) = ЛZ^, г')К(Z^, Г —- 2, г)йТЛ' + ^о(2)8(г), (1.3)
о Т
где Т = (X, п), йТ = йХйц0(л), причем интегрирование по мере означает суммирование по всем различным парам п = (7, у).
Для уравнения Смолуховского, соответственно (см. [3]), X = (А, 11, ..., 1А), N < А0, ^ е 1, А0, -размер 7-й частицы, а для уравнения БольцманаX = (уь ..., уА), у7 - скорость 7-й частицы. Получаемая из фазовой плотности г) одночастичная плотность является приближенным решением соответствующего кинетического уравнения в предположении молекулярного хаоса (см. [1]).
Для целей настоящей работы целесообразно фиксировать X - состояние А-частичного ансамбля в начале свободного временного пробега. Ядро соответствующего уравнения (1.3) имеет вид
К(Тг^ Т, г) = к(г^ г|Х')К1 (Х| п),
где
к(г' —- ^) = А (X')Е(X, г - г'),
причем Е^, г) = ехр^А^}. Условная переходная плотность К^Я —► X|п) определяет преобразование ансамбля в результате взаимодействия пары с номером п = (7, у). Для уравнения Смолуховского это - замена взаимодействующих частиц на одну частицу размера 17 + 1у, причем N = N - 1, а
для уравнения Больцмана - замена (у), уУ) —- (у7, уу-) соответственно заданной индикатрисе.
Обычно при решении уравнения Каца вычисляют функционалы JH(f) от потока частиц Ф^, г) в заданные моменты времени. На основе теории локальных оценок (см. [4, разд. 6.9]) получается равенство
г
JH(г) = ^^Ф^, г)dX = Ця(X, г - Г)Т, Г)йТйе,
о Т
где Н(X) е ЬН (X, г) = Н^Е^, г).
Рассмотрим цепь Маркова (Тп, гп), п е 0, 1, ..., с нормированной плотностью перехода
Р(Т, г' —► Т, г) = р(г' — )Рх(п|Xр(X — X]п)
и нормированной плотностью начального состояния Р0(Т, г) = Р0(Т)5(0. Случайные веса вводятся по формулам
бо = Т)/Ро(Т), бп = б„-гб(Тп-1, гпгп), б(Тп-1, гп-1; Тп, гп) = К(Т\ г^ Т, г)/Р(Т, г^ Т, г). Величина б очевидным образом факторизуется. Для оценки функционала JH(t) можно исполь-
зовать весовой вариант оценки "по столкновениям" (в начале пробегов):
V
% = X Т - хп), V = тах{ п: хп < Т, п = 0, 1.....N>-1}, (1.4)
г = о
где Т _ момент времени, в который оценивается рассматриваемый функционал. Моделирование цепи Маркова и вычисление значения % происходят по следующей схеме:
1) для t0 = 0 моделируется начальное состояние 20 соответственно плотности Р0(2 и вычисляется б,,;
2) согласно плотности Р(2п _ ь хп _ 1 —» 2п, хп) выбирается хп и моделируется 2п, причем вес преобразуется после каждого элементарного перехода;
3) если хп < Т, то вычисляется слагаемое бпН (Хп, Т _ tn), иначе цепь обрывается.
Интегральное уравнение (1.3) может быть получено непосредственно из этой марковской модели как соотношение интегрального баланса для суммарной плотности состояний.
Заметим, что оценку по столкновениям вида (1.4) можно построить и в случае, когда фиксируется фазовое состояние Х'п ансамбля в конце пробега. Это осуществляется путем очевидного
осреднения значений Н (Хп, Т _ хп) соответственно распределению новых состояний Хп после взаимодействия выбранной пары частиц. Получаемое в результате выражение имеет вид
V
%1 = X хп, т - хп)+бо Н( Хо, Т),
г = 1
где
Н (Хп, Т - Г„) = | К (Х^^^Хп) я) Н (Хп, Т - Г„) Фо (п) йХп, (1.5)
а бШ) _ значение веса непосредственно после выбора хп. Поскольку Н является условным средним величины Н, то можно ожидать, что D%1 < D%. Однако вычисление интеграла (или суммы _ для уравнения Смолуховского) может быть слишком трудоемким. Кроме того, в случае ценностного моделирования значение Н фактически получается методом "выборки по важности" для (1.5) и D% ~ D%1. Это было подтверждено расчетами для рассмотренной в разд. 3 модельной задачи коагуляции.
2. ЦЕННОСТНОЕ МОДЕЛИРОВАНИЕ ДЛЯ КИНЕТИЧЕСКИХ УРАВНЕНИЙ
Согласно терминологии, принятой в теории весовых методов Монте-Карло, функцией ценности ф*(Х, х) называется значение функционала JH для источника с плотностью 5(Х' _Х)5(х' _ 0, т.е.
V
ф* (Х, х) = Е%(Х, 0, где %(Х, о = Н(Х, Т - х) + X бпН(Хп, Т - 1п).
п = 1
Функция ф*(Х, х) является решением интегрального уравнения, сопряженного к уравнению (1.3), и тем самым в предположении молекулярного хаоса может быть приближена решением сопряженного нелинейного кинетического уравнения (см. [5]). Однако, как показано ниже, иногда практически удовлетворительное приближение ф* (Х, х) может быть получено непосредственно на основе оценки функционала JH для модельного кинетического уравнения.
С учетом вида переходной плотности Р для моделирования перехода (2', х') —► (2, х) целесообразно реализовать три элементарных перехода. Согласно [4, разд. 4.5], моделирование называется ценностным, если каждый элементарный переход реализуется соответственно номиро-ванной плотности, полученной путем умножения физической плотности на соответствующую вспомогательную функцию ценности, которая равна ф*(Х, х) лишь для финального элементарного перехода. Для первого и второго переходов в рассматриваемой задаче необходимо использо-
вать следующие ценностные множители:
ф*(г; X) = Цо(п)А_1(X)Kг(X —- XIп)ф*(X, t)dXdц(п),
ф* (п; X', г) = |К1 (X —»- XIп)ф*(X, г)dX.
Следующее утверждение является следствием общей теоремы о ценностном моделировании (см. [6]) всех элементарных переходов.
Теорема 1. Если реализуются нормированные плотности элементарных переходов вида
Cok(tt\X')ф*(t; X), Q-^-ф*(п; X, t), C2K(X ^ X\п)ф*(X, t),
a(n) M(X)'
mo справедливо равенство D^ = 0.
Теорема показывает, что если используются достаточно точные приближенные функции ф* , ф* и ф* , то величина D^ может быть малой.
3. ЦЕННОСТНОЕ МОДЕЛИРОВАНИЕ ДЛЯ РЕШЕНИЯ ТЕСТОВОЙ ЗАДАЧИ КОАГУЛЯЦИИ
3.1. Рассмотрим ценностное моделирование для оценки решения тестовой задачи для уравнения (1.1) при
kj = 1 и и,(0) = 1, n(0) = 0, I > 2. Это решение известно (см. [7]):
(а 1 ( 0.51 V-1
n'( t) =-21 i + 051 , 1 > 1.
(1 + 0.51) V1 + 0 5ty
Для такой задачи а(п) = N- , A(X) = 0.5N(N - 1)/N0, где N0 - начальное количество частиц в ансамбле.
Решается задача оценки величины jH (T) ~ n1(T), т.е. среднего числа (частоты) мономеров в момент времени T, при этом
Н (X) = Н (X) = Nj( X) / N0.
Теорема 2. Для тестовой задачи коагуляции имеем ф*^ t) = C(N, t)N1 при
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.