Статистические алгоритмы ценностного моделирования для решения уравнения Смолуховского Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 13, Специальный выпуск 4, 2008

Статистические алгоритмы ценностного моделирования для решения уравнения Смолуховекого*

М. А. Коротченко Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: kmaria@osmf . sscc. ru

Modifications of statistical algorithms applied for solution of the Smoluchowski equation are considered. Efficiency of estimations for the total monomer and dimer quantity in an ensemble is studied. The estimate is obtained using the value modeling of the free path, interaction pair number, and their combinations for various value functions.

Введение

В настоящей работе рассматриваются “ценностные” весовые модификации статистического моделирования для численного решения уравнения Смолуховекого в пространственно-однородном случае.

Уравнение Смолуховекого описывает широкий класс процессов коагуляции. Мы будем рассматривать случай парного слипания частиц. Данное уравнение в пространственно-однородном случае имеет вид

П3 (*) “ X! КиЩ Пг(°) = п<10)’ 1 - °>

г+]=1 г>1

где пг(£) — частота (числовая плотность) частиц размера I в момент времени £, причем размер частицы I принимает натуральные значения; К^ — коэффициенты коагуляции (предполагаются заданными). Под численным решением уравнения понимается нахождение линейных функционалов от функции пг(£).

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00046).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

F (Z,t) = J J F (Z',t') K (Z',t' ^ Z,t) dZ'dt' + Fo (Z)5(t), (1)

o Z

где Z = (Х,7г), X = i \. / |. • • • ,lN), N < N0, /.. G 1,N0 — размер i-ж частицы; dZ = dXd^0 (п), причем интегрирование по мере ^0 означает суммирование по всем различным номерам пар п = (i,j). Ядро уравнения (1) имеет вид [2]

K(Z',t' Z,t)= r(t' ^ - I).

Здесь использованы следующие обозначения: r(t' ^ t\X') = A(X') exp{ — (t—t')A(X')}, —

N-1 N

pdf длины “свободного пробега” ансамбля; A (X) = a (N, h, lj), N = N' — 1,

i=l j=i+l ж

. (li,lj ^ l) ^ли N > ^,

a (п) = a (N,li,lj) = ^ =

0, если иначе,

причем a(n)/A(X) — вероятность выбора значения п; k(li, lj ^ l) = N0-1KijSli+lj определяет преобразование ансамбля при взаимодействии пары с номером п = (i,j), в результате которого происходит замена взаимодействующих частиц на одну частицу размера li + lj ‘

Обычно при решении уравнения вычисляют функционалы JH(t) = (Ф, H) от потока “частиц” Ф(X,t) в заданные моменты времени, В [1] получено равенство

t Jh (t) = J j H(X, t — t')F(Z, t')dZdt' = (F, H),

o Z

где H(X) G Lж, H(X,t) = H(X) exp{—A(X)t}.

Рассмотрим цепь Маркова (Zn, tn), n G 0,1,..., с нормированной плотностью перехода

P(Z',t' ^ Z,t)= p(t' ^ t\X')Pi(AX')P2(X' ^ X\п)

и нормированной плотностью начального состояния P0(Z,t) = P0(Z)S(t). Случайные веса вводятся по формулам

Q0 = F0(Z0)/P0(Z0), Qn = Qn-\Q(Zn-l,tn-l; Z^ tn), Q(Zn-l, tn-i; Zn-, tn) = K(Z' , £ ^ Z,t)/P(Z' , £ ^ Z,t).

Для оценки функционала JH (T) можно использовать весовой вариант оценки “по столкновениям” (в начале пробегов):

V ^ = ^2, QnH(Xn,T — tn), v = max{n : tn <T, n =0,1,..., N0 — 1},

i=0

T

Моделирование цепи (2п, гп) и вычисление значения £ происходит по следующей схеме,

1, Для г0 = 0 моделируется начальное с остоянпе плот пости Р0(2) и вычисляется Q0.

2, Согласно плотности Р(2п-1,гп_1 ^ 2п, гп) выбирается гп и моделируется 2п, причем вес преобразуется после каждого элементарного перехода,

3, Если гп < Т, то вычисляется слагаемое QnH(Xn, Т — гп), иначе цепь обрывается,

1. Ценностное моделирование для уравнения Смолуховского

Согласно определению [3], функцией ценности ф*(Х,г) называется значение функционала Зн для источника с плотностью 5(Х' — X)5(£ — г), т, е,

Ф*(Х,г) = Е£(х,г), где £(х,г) = Н(Х,Т — г) + Е QnH(Xn,T — 1п). п= 1

Функция Ф*(Х, г) является решением интегрального уравнения, сопряженного с уравнением (1), следовательно, может быть приближена решением нелинейного сопряженного уравнения [4].

С учетом вида, переходной плотности Р для моделирования перехода (2', г') ^ (2, г) целесообразно реализовать два элементарных перехода (выбор г и п), Согласно [3] моделирование называют ценностным, если каждый элементарный переход реализуется соответственно нормированной плотности, полученной умножением “физической” плотности на соответствующую вспомогательную функцию ценности, которая равна Ф* (Х, г) лишь для финального элементарного перехода. Если при этом используется точная функция ценности, то дисперсия = 0, а если приближенная — то величина может быть малой [2].

2. Ценностное моделирование для решения тестовой задачи коагуляции

Далее будет рассмотрено ценностное моделирование с целью оценки решения тестовой

задачи для уравнения (1) при К^ = 1 и п1(0) = 1, щ(0) = 0, I > 2. Это решение

известно [5]:

1 ( 0.5г У-1 , _

~ (1 + О.Ы)2 (1 + 0.5*) ’ - • П

_1 N Ш — 1)

Для такой задачи а(тг) = , А(Х) = ——-------.

2No

Оценка величины п1(Т). Рассмотрим задачу оценки величины ^(Т) Щ(Т),

т, е, среднего числа (частоты) мономеров в момент времени Т, при этом Н(Х) = Н1(Х) = Хг(Х)

Х0 '

В [2] показано, что ^^(Т) — п1(Т)| = О^_1). Принимая во внимание последнее равенство, рассмотрим “точную” временную функцию ценности ф\(г) = п1(Т£ — г)1£(г), а также приближения ф*1(г) ъ Ф1 к ней:

ф*(г) = І£(г),

Й(*) = ехр (I,(t) « щ(Т, - t)Ut) = (1+0^_())2, (3)

где T = T + е; I=(t) — индикатор отрезка [0,T=], Величина е продолжения временного промежутка, на котором строится цепь Маркова, вводится с целью обрыва цепи Маркова, Функционал H зануляется при пересечении уровня T, а плотность моделируется па расширенном временном интервале [0,T=], На этом интервале цепь Маркова имеет бесконечное число переходов, но функционал H отличен от 0 только в интервале [0,T], е

Для разных случаев выбора ф* свободный пробег моделируетея в [t',T=) соответственно плотностям, пропорциональным различным функциям, приведенным ниже. Во

Ts

всех случаях нормировочные константы С\(Х', t) выбираются из уеловия J f (t) dt = 1.

• При ф* « ф\ плотность fi(t) = Ci(X',t') exp{-A(X')(t — t)}.

2

• При ф* и ф\ плотность f2(t) = C2(X',t') exp{ —(A(X') - A0)(t - t')}, A0 = -——.

2 + T =

• Пр и ф* = ф\ плотность f0(t) = С0(Х', t )nl(T= — t) exp{—A(X' )(t — t)}, и ее можно

fi ( t)

< exp{-A{X’){t - t')} = ^

Со(Х',г')- ^ у " С1(Х',г')

г

величину д = А(Х') ехр{—А(Х')(г — г')}//г(г).

После моделирования времени свободного пробега т = г — г' выполняется ценностное моделирование ном,ера, пары взаимодействующих частиц. Обозначим: N — общее количество частиц, N1 — количество мономеров в ансамбле в конце свободного пробега (перед выбором значения п).

Каждая из всевозможных пар для столкновения попадает в одну из неперееека-ющихся групп в зависимости от изменения количества мономеров после столкновения

1

2

“0-пары”):

“1-пары” — пары вида {мономер, полимер}, их число N = N1(^ — N1);

“2-пары” — пары вида {мономер, мономер}, их число N2 = Щ(Щ — 1)/2;

“0-пары” — пары вида {полимер, полимер}, их число N = (№ — Щ)(№ — N1 — 1)/2.

2

Представим “физическое” равномерное распределение Ро{%,]) = Р0 = номера взаимодействующей пары в следующем рандомизированном виде:

1 = X! Ро ({,3)= Р1 ^ /1 ({,Э)+ Р2 ^ ({,3)+ ро ^ ?о ({,1), (4)

1<i<j<N (г,])_1-пара (г,])_2-пара (г,])-О-пара

где /к (%,]) = 1 /Мк — равномерная вероятность вы бора пары (%,]) из множес тва к-пар, а рк = Р<^к — вероятность выбора типа пары,

В целях “сохранения” мономеров моделирование пары будем производить соответственно представлению (4), в котором заменим вероятности рк на величины дк, пропорциональные числу оставшихся мономеров:

, Р1И-1) мщ-2) : РоЛ" ЛГ;(ЛГ'-2)

<11 - ---~-----> <72 - -----£----, О - - -----—------.

Такая модификация учитывается при построении веса: 0 = 0'рк/дк. В [2] показано, что

использованная “ценность” номера пары соответствует точной функции ф*.

Оценка величины п1(Т)+ п2(Т). Далее рассмотрим задачу оценки функционала (12)

3% ’{Т) п\(Т)+п2(Т), т. е, суммарного числа мономеров и димеров в заданный момент

времени Т, при этом Н(Х) = Н\2(Х) = ———^

Nо

В качестве временной функции ценности будем использовать приближение

8(1 + Т) (( 3 1

—з ехр

і\ & Іє(і)(пі(Тє - і) + п2(Тє - і)).

(2 + Т£)6 ^ [\2 + Т 1 + Т£

В этом случае свободный пробег моделируется в [г',Те) соответственно плотности:

2Т + 1

Ы*) = С12(Х', И) ехр{—(А(Х') - А2Щ, А £ +

После выбора і вес домпожается па величину д =

(2 + Тє)(1+ Тє)‘

А(Х/)ехр{-А(Х/)(*-0} /і2(*)

Для ценностного моделирования ном,ера пары, рассмотрим плотность, пропорциональную величине N1(X) + N2(X). Обозначим N2 — количество димеров в ансамбле в конце свободного пробега. Далее, считая полимерами частицы с размером I > 3, все множество возможных пар для столкновения разобьем на шесть неперееекающих-ся групп в зависимости от изменения общего количества мономеров и димеров после столкновения (это количество может уменьшиться на 1, на 2 и не измениться):

1) “1-пары” вида {мономер, мономер}, их число Ыи = N'^N1 — 1)/2;

2) “1-пары” вида { мономер, полимер}, их число Ы1к = ^1(^ — (N1 + N2))',

3) “1-пары” вида {димер, полимер}, их число Ы2к = Щ(^ — (N1 + N2))',

4) “2-пары” вида {димер, димер}, их число М22 = Щ(Щ — 1)/2\

5) “2-пары” вида {мономер, димер}, их число М12 = N'1N'2^,

6) “0-пары” вида {полимер, полимер}, их число Ыкк = N(^ — 1)/2 — (Ы" + М12 +

N1к + -^22 + N2к )■

В аналогичном (4) представлении физического распределения Р0 номера взаимодействующей пары вероятности выбора типа пары имеют вил рг] = N1]Р0, г] € {1, 2, к}.

В целях “сохранения” мономеров и димеров моделирование пары будем производить в соответствии с представлением, аналогичным (4), где вероятности рг] заменим на д^, пропорциональные сумме оставшихся в системе мономеров и димеров:

«и = (лг; + ^ - 1)^г, «1* = (лг; + - 1)^г, <й* = (лг; + лг;-1)^,

= (лг; + ЛГ' - 2)^2, д22 = (ЛГ[ + ЛГ' - 2)^,

<7и = (лг; + ЛГ' — 0)М

Здесь С

0'Ри /д%з-

нормировочная константа. Такая модификация учитывается в весе: 0

3. Результаты численных экспериментов

В численных экспериментах для задачи коагуляции, результаты которых приведены ниже, использовались следующие расчетные параметры: е = 10_5, начальное число

частиц = 100, число траекторий 106. Выбранное значение е дает почти минимальную (по е) дисперсию, не сильно увеличивая среднее число столкновений по сравнению с прямым моделированием,

В качестве функционалов выбраны З^^1 (Т), З^2 (Т) и З^21(Т) — среднее число мономеров, димеров или суммарное число мономеров и димеров соответственно в момент времени Т. Смещение оценок относительно П\(Т), п2(Т) и П\(Т) + п2(Т) соответственно объясняется конечностью N0. Согласно [2] оно имеет порядок О (Ж-1),

Обозначения в таблицах: а — средпеквадратичеекое уклонение; £с — время расчета на РС РепМит 4 (3,5 ГГц); Ба и Бь — трудоемкости [3] алгоритмов прямого статистического моделирования (подсчет числа мономеров, димеров или суммы мономеров и димеров при £ = Т) и ценностного моделирования соответственно.

Дополнительные расчеты показали, что ценностное моделирование лишь длины пробега не дает выигрыша по сравнению с прямым моделированием, В то же время результаты численных оценок функционалов З^1 (Т) (табл. 1) и З^22 (Т) (табл. 2) показывают, что приближенно ценностное моделирование обоих элементарных переходов радикально уменьшает трудоемкость расчетов. Оценка числа димеров может быть получена как разность З^2 (Т) = З^^22 (Т) — З^2(Т), причем трудоемкость оптимальной оценки величины З(1 (Т) здесь пренебрежимо мала.

Таблица 1. Оценка функционала 3н\т) с использованием ценностного моделирования номера пары и различных “ценностей” свободного пробега для уравнения Смолуховского

<РХ 3^{Т) а с

Т = 1, щ( 1) = 4.4Т 14 • 10“ 1

ф* 4.46928 • КГ1 4.9 • 10“5 74 0.86

ФХ 4.46922 • КГ1 2.5 • 10“6 87 280

Т = 5, щ(5) = 8.1632 • 10“ 2

ф* 8.2380 • 10“2 1.6 • 10“5 113 2.0

ФХ 8.2382 • 10“2 1.1 • 10“6 142 337.3

Т = 10, щ(10) = 2.7777 • 10 _2

ФХ 2.8047 • 10“2 6.6 • 10_6 124 4.3

ФХ 2.8044 • 10-2 5.9 • 10“7 153 432.5

Т = 20, щ(20) = 8.2644 • 10 -з

ФХ 8.3410 • 10_3 2.4 • 10_6 129 12.7

ФХ 8.3405 • 10“3 2.7 • 10“7 161 808.4

(12)

Таблица 2. Оценка функционала 3КН (Т) с использованием ценностного моделирования номера пары и с приближенной функцией ценности свободного пробега ф1

<Р* 3^2\Т) а с Ба/Б,

Т= 1, щ(1)+п2(1) = 5.9259 • 10“1

ФХ 5.9574 • 10-1 6.9 • 10-5 108 0.31

Г = 5, щ(5) +п2(5) = 1.3994 • 10“1

ФХ 1.4161 • 10-1 7.5 • 10_6 197 9.50

Г =10, щ(10) + п2(10) = 5.0926 • 10“2

ФХ 5.1607 • 10-2 4.2 • 10_6 197 12.82

Т = 20, щ(20) +п2(20) = 1.5777- 10“2

ФХ 1.5991-10“2 1.7-10“6 206 26.85

Расчеты е функцией ценности (2) дали неудовлетворительные результаты даже для числа начальных частиц N порядка нескольких тысяч; это объясняется выполнением равенства З#1^) = п1(Т) только асимптотически при М0 ^ ж. Таким образом, приближение (3) в сочетании с ценностным моделированием номера пары оказалось достаточно эффективным,

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

Автор выражает благодарность чл.-корр, РАН Г,А, Михайлову за ценные замечания и предложения в процессе подготовки статьи.

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

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

[2] Михайлов Г.А., Коротченко М.А., Рогазинский С.В. Ценностные алгоритмы статистического моделирования для нелинейных кинетических уравнений // Докл. РАН. 2007. Т. 415, № 1. С. 26-30.

[3] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование (метод Монте-Карло). М.: Изд. центр “Академия”, 2006.

[4] Марчук Г.И., Агошков В.И., Шутяев В.П. Сопряженные уравнения и методы возмущений в нелинейных задачах математической физики. М.: Наука, 1993.

[5] ВОЛОЩУК В.М. Кинетическая теория коагуляции. Л.: Гидрометеоиздат, 1984.

Поступила в редакцию 20 февраля 2008 г.