close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Вычислительные технологии
Том 13, Специальный выпуск 4, 2008
Статистические алгоритмы
ценностного моделирования
для решения уравнения Смолуховского?
М. А. Коротченко
Институт вычислительной математики
и математической геоизики СО АН, Новосибирск, оссия
e-mail: kmariaosmf.ss.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.
Введение
В настоящей работе рассматриваются ценностные весовые модиикации статистического моделирования для численного решения уравнения Смолуховского в пространственно-однородном случае.
Уравнение Смолуховского описывает широкий класс процессов коагуляции. Мы будем рассматривать случай парного слипания частиц. Данное уравнение в пространственно-однородном случае имеет вид
X
?nl (t)
1 X
(0)
=
Kij ni (t) nj (t) ?
Kil ni (t) nl (t), nl (0) = nl , l ? 0,
?t
2
i?1
i+j=l
где nl (t) частота (числовая плотность) частиц размера l в момент времени t, причем
размер частицы l принимает натуральные значения; Kij коэициенты коагуляции
(предполагаются заданными). Под численным решением уравнения понимается нахождение линейных ункционалов от ункции nl (t).
ешение кинетического уравнения можно оценивать с помощью моделирования марковского процесса эволюции N -частичного ансамбля. Введение номера взаимодействующей пары в число координат азового пространства позволило в [1? сормулировать
новое интегральное уравнение, удобное для построения стандартных весовых модиикаций статистического моделирования рассматриваемой многочастичной системы, так
как его ядро содержит сингулярности лишь в виде сомножителей:
?
абота выполнена при инансовой поддержке оссийского онда ундаментальных исследований
(грант ќ 06-01-00046).
Институт вычислительных технологий Сибирского отделения оссийской академии наук, 2008.
68
Ценностное моделирование для решения уравнения Смолуховского
F (Z, t) =
Zt Z
0
F (Z ? , t? ) K (Z ? , t? ? Z, t) dZ ?dt? + F0 (Z)?(t),
69
(1)
Z
где Z = (X, ?), X = (N, l1 , · · · , lN ), N ? N0 , li ? 1, N0 размер i-й частицы; dZ =
dXdµ0 (?), причем интегрирование по мере µ0 означает суммирование по всем различным номерам пар ? = (i, j). Ядро уравнения (1) имеет вид [2?
K(Z ? , t? ? Z, t) = r(t? ? t|X ? )
a(?)
k(li , lj ? l).
A(X ? )
Здесь использованы следующие обозначения: r(t? ? t|X ? ) = A(X ? ) exp{?(t?t? )A(X ? )}, NP
?1 P
N
pdf длины свободного пробега ансамбля; A (X) =
a (N, li , lj ), N = N ? ? 1,
i=1 j=i+1
? ?
? P k (l , l ? l) , если N ? 2,
i j
a (?) ? a (N, li , lj ) =
l=1
?
0, если иначе,
причем a(?)/A(X) вероятность выбора значения ? ; k(li , lj ? l) = N0?1 Kij ?li +lj определяет преобразование ансамбля при взаимодействии пары с номером ? = (i, j), в результате которого происходит замена взаимодействующих частиц на одну частицу размера
li + lj .
Обычно при решении уравнения вычисляют ункционалы JH (t) = (?, H) от потока
частиц ?(X, t) в заданные моменты времени. В [1? получено равенство
JH (t) =
Zt Z
0
Z
e
H?(X, t ? t? )F (Z, t?)dZdt? = (F, H),
e
где H(X) ? L? , H(X,
t) = H(X) exp{?A(X)t}.
ассмотрим цепь Маркова (Zn , tn ), n ? 0, 1, . . . , с нормированной плотностью перехода
P (Z ?, t? ? Z, t) = p(t? ? t|X ? )P1 (?|X ?)P2 (X ? ? X|?)
и нормированной плотностью начального состояния P0 (Z, t) = P0 (Z)?(t). Случайные
веса вводятся по ормулам
Q0 = F0 (Z0 )/P0 (Z0 ), Qn = Qn?1 Q(Zn?1 , tn?1 ; Zn , tn ),
Q(Zn?1 , tn?1 ; Zn , tn ) = K(Z ? , t? ? Z, t)/P (Z ?, t? ? Z, t).
Для оценки ункционала JH (T ) можно использовать весовой вариант оценки по столкновениям (в начале пробегов):
?=
?
X
Qn H?(Xn , T ? tn ),
? = max{n : tn < T, n = 0, 1, . . . , N0 ? 1},
i=0
где T момент времени, в который оценивается рассматриваемый ункционал.
70
М. А. Коротченко
Моделирование цепи (Zn , tn ) и вычисление значения ? происходит по следующей
схеме.
1. Для t0 = 0 моделируется начальное состояние Z0 соответственно плотности P0 (Z)
и вычисляется Q0 .
2. Согласно плотности P (Zn?1, tn?1 ? Zn , tn ) выбирается tn и моделируется Zn , причем вес преобразуется после каждого элементарного перехода.
3. Если tn < T , то вычисляется слагаемое Qn H?(Xn , T ? tn ), иначе цепь обрывается.
1. Ценностное моделирование для уравнения Смолуховского
Согласно определению [3?, ункцией ценности ?? (X, t) называется значение ункционала JH для источника с плотностью ?(X ? ? X)?(t? ? t), т. е.
?
?? (X, t) = E?(X,t) , где ?(X,t) = H?(X, T ? t) + ? Qn H?(Xn , T ? tn ).
n=1
Функция ?? (X, t) является решением интегрального уравнения, сопряженного с уравнением (1), следовательно, может быть приближена решением нелинейного сопряженного
уравнения [4?.
С учетом вида переходной плотности P для моделирования перехода (Z ? , t? ) ? (Z, t)
целесообразно реализовать два элементарных перехода (выбор t и ? ). Согласно [3? моделирование называют ценностным, если каждый элементарный переход реализуется соответственно нормированной плотности, полученной умножением изической
плотности на соответствующую вспомогательную ункцию ценности, которая равна
?? (X, t) лишь для инального элементарного перехода. Если при этом используется
точная ункция ценности, то дисперсия D? = 0, а если приближенная то величина
D? может быть малой [2?.
2. Ценностное моделирование для решения
тестовой задачи коагуляции
Далее будет рассмотрено ценностное моделирование с целью оценки решения тестовой
задачи для уравнения (1) при Kij ? 1 и n1 (0) = 1, nl (0) = 0, l ? 2. Это решение
известно [5?:
l?1
1
0.5t
nl (t) =
, l ? 1.
(2)
(1 + 0.5t)2 1 + 0.5t
N(N ? 1)
.
2N0
(1)
Оценка величины n1 (T ). ассмотрим задачу оценки величины JH (T ) ? n1 (T ),
т. е. среднего числа (частоты) мономеров в момент времени T , при этом H(X) = H1 (X) ?
N1 (X)
.
N0
(1)
В [2? показано, что |JH (T ) ? n1 (T )| = O(N0?1 ). Принимая во внимание последнее
равенство, рассмотрим точную временн
ую ункцию ценности ??1 (t) = n1 (T? ? t)I? (t),
?
?
а также приближения ??1 (t) и ???1 к ней:
Для такой задачи a(?) = N0?1 , A(X) =
???1 (t) = I? (t),
71
Ценностное моделирование для решения уравнения Смолуховского
2
2
I? (t)
2
(3)
=
exp
t I? (t) ? n1 (T? ? t)I? (t) =
,
2 + T?
2 + T?
(1 + 0.5(T? ? t))2
где T? = T + ?; I? (t) индикатор отрезка [0, T? ]. Величина ? продолжения временного
промежутка, на котором строится цепь Маркова, вводится с целью обрыва цепи Маркова. Функционал H зануляется при пересечении уровня T , а плотность моделируется
на расширенном временном интервале [0, T? ]. На этом интервале цепь Маркова имеет
бесконечное число переходов, но ункционал H отличен от 0 только в интервале [0, T ].
Введение ? оставляет оценку ункционала несмещенной [3?, но влияет на дисперсию.
Для разных случаев выбора ?? свободный пробег моделируется в [t? , T? ) соответственно плотностям, пропорциональным различным ункциям, приведенным ниже. Во
ZT?
? ?
всех случаях нормировочные константы Ci (X , t ) выбираются из условия fi (t) dt = 1.
????1 (t)
t?
? При ?? ? ???1 плотность f1 (t) = C1 (X ? , t? ) exp{?A(X ? )(t ? t? )}.
? При ?? ? ????1 плотность f2 (t) = C2 (X ? , t? ) exp{?(A(X ? ) ? A0 )(t ? t? )}, A0 =
2
.
2 + T?
? При ?? = ??1 плотность f0 (t) = C0 (X ? , t? )n1 (T? ? t) exp{?A(X ? )(t ? t? )}, и ее можно
моделировать методом исключения [3? с мажорантой f1 (t):
f0 (t)
f1 (t)
? exp{?A(X ? )(t ? t? )} =
.
?
?
C0 (X , t )
C1 (X ? , t? )
После выбора момента времени следующего столкновения t вес домножается на
величину q = A(X ? ) exp{?A(X ? )(t ? t? )}/fi (t).
После моделирования времени свободного пробега ? = t ? t? выполняется ценностное моделирование номера пары взаимодействующих частиц. Обозначим: N ? общее
количество частиц, N1? количество мономеров в ансамбле в конце свободного пробега
(перед выбором значения ? ).
Каждая из всевозможных пар для столкновения попадает в одну из непересекающихся групп в зависимости от изменения количества мономеров после столкновения
(количество мономеров может уменьшиться на 1 в результате взаимодействия 1-пары,
на 2 в результате взаимодействия 2-пары и не измениться при взаимодействии
0-пары):
1-пары пары вида {мономер, полимер }, их число N1 = N1? (N ? ? N1? );
2-пары пары вида {мономер, мономер }, их число N2 = N1? (N1? ? 1)/2;
0-пары пары вида {полимер, полимер }, их число N0 = (N ? ? N1? )(N ? ? N1? ? 1)/2.
2
Представим изическое равномерное распределение P0 (i, j) ? P0 =
?
N (N ? ? 1)
номера взаимодействующей пары в следующем рандомизированном виде:
X
X
X
X
1?
P0 (i, j) = p1
f1 (i, j) + p2
f2 (i, j) + p0
f0 (i, j),
(4)
1?i<j?N
(i,j)?1-пара
(i,j)?2-пара
(i,j)?0-пара
где fk (i, j) = 1/Nk равномерная вероятность выбора пары (i, j) из множества k -пар,
а pk = P0 Nk вероятность выбора типа пары.
В целях сохранения мономеров моделирование пары будем производить соответственно представлению (4), в котором заменим вероятности pk на величины qk , пропорциональные числу оставшихся мономеров:
q1 =
p1 (N1? ? 1)
p2 (N1? ? 2)
p0 N1?
N ? (N ? ? 2)
, q2 =
, q0 =
, C = E(N1 ) = 1
.
C
C
C
N?
72
М. А. Коротченко
Такая модиикация учитывается при построении веса: Q = Q? pk /qk . В [2? показано, что
использованная ценность номера пары соответствует точной ункции ?? .
Оценка величины n1 (T )+n2 (T ). Далее рассмотрим задачу оценки ункционала
(12)
JH (T ) ? n1 (T )+n2 (T ), т. е. суммарного числа мономеров и димеров в заданный момент
N1 (X) + N2 (X)
времени T , при этом H(X) = H12 (X) ?
.
N0
В качестве временной ункции ценности будем использовать приближение
8(1 + T? )
3
1
?
exp
??12 (t) = I? (t)
?
t ? I? (t)(n1 (T? ? t) + n2 (T? ? t)).
2 + T? 1 + T?
(2 + T? )3
В этом случае свободный пробег моделируется в [t? , T? ) соответственно плотности:
f12 (t) = C12 (X ? , t? ) exp{?(A(X ? ) ? A2 )t},
A2 =
2T? + 1
.
(2 + T? ) (1 + T? )
A(X ? ) exp{?A(X ? )(t ? t? )}
После выбора t вес домножается на величину q =
.
f12 (t)
Для ценностного моделирования номера пары рассмотрим плотность, пропорциональную величине N1 (X) + N2 (X). Обозначим N2? количество димеров в ансамбле
в конце свободного пробега. Далее, считая полимерами частицы с размером l ? 3,
все множество возможных пар для столкновения разобьем на шесть непересекающихся групп в зависимости от изменения общего количества мономеров и димеров после
столкновения (это количество может уменьшиться на 1, на 2 и не измениться):
1) 1-пары вида {мономер, мономер }, их число N11 = N1? (N1? ? 1)/2;
2) 1-пары вида { мономер, полимер }, их число N1k = N1? (N ? ? (N1? + N2? ));
3) 1-пары вида {димер, полимер }, их число N2k = N2? (N ? ? (N1? + N2? ));
4) 2-пары вида {димер, димер }, их число N22 = N2? (N2? ? 1)/2;
5) 2-пары вида {мономер, димер }, их число N12 = N1? N2? ;
6) 0-пары вида {полимер, полимер }, их число Nkk = N ? (N ? ? 1)/2 ? (N11 + N12 +
N1k + N22 + N2k ).
В аналогичном (4) представлении изического распределения P0 номера взаимодействующей пары вероятности выбора типа пары имеют вид pij = Nij P0 , i, j ? {1, 2, k}.
В целях сохранения мономеров и димеров моделирование пары будем производить
в соответствии с представлением, аналогичным (4), где вероятности pij заменим на qij ,
пропорциональные сумме оставшихся в системе мономеров и димеров:
p11
p1k
p2k
, q1k = (N1? + N2? ? 1) , q2k = (N1? + N2? ? 1) ,
C
C
C
p12
p22
pkk
?
?
?
?
?
?
= (N1 + N2 ? 2) , q22 = (N1 + N2 ? 2) , qkk = (N1 + N2 ? 0)
.
C
C
C
q11 = (N1? + N2? ? 1)
q12
Здесь C нормировочная константа. Такая модиикация учитывается в весе: Q =
Q? pij /qij .
3. езультаты численных экспериментов
В численных экспериментах для задачи коагуляции, результаты которых приведены
ниже, использовались следующие расчетные параметры: ? = 10?5, начальное число
Ценностное моделирование для решения уравнения Смолуховского
73
частиц N0 = 100, число траекторий 106 . Выбранное значение ? дает почти минимальную
(по ?) дисперсию, не сильно увеличивая среднее число столкновений по сравнению с
прямым моделированием.
(1)
(2)
(12)
В качестве ункционалов выбраны JH (T ), JH (T ) и JH (T ) среднее число мономеров, димеров или суммарное число мономеров и димеров соответственно в момент
времени T . Смещение оценок относительно n1 (T ), n2 (T ) и n1 (T ) + n2 (T ) соответственно
объясняется конечностью N0 . Согласно [2? оно имеет порядок O(N0?1 ).
Обозначения в таблицах: ? среднеквадратическое уклонение; tc время расчета
на PC Pentium 4 (3.5 ц); Sd и Sv трудоемкости [3? алгоритмов прямого статистического моделирования (подсчет числа мономеров, димеров или суммы мономеров и
димеров при t = T ) и ценностного моделирования соответственно.
Дополнительные расчеты показали, что ценностное моделирование лишь длины пробега не дает выигрыша по сравнению с прямым моделированием. В то же время резуль(1)
(12)
таты численных оценок ункционалов JH (T ) (табл. 1) и JH (T ) (табл. 2) показывают,
что приближенно ценностное моделирование обоих элементарных переходов радикально уменьшает трудоемкость расчетов. Оценка числа димеров может быть получена как
(2)
(12)
(1)
разность JH (T ) = JH (T ) ? JH (T ), причем трудоемкость оптимальной оценки вели(1)
чины JH (T ) здесь пренебрежимо мала.
Т а б л и ц а
1. Оценка ункционала
(1)
JH (T )
с использованием ценностного моделирования
номера пары и различных ценностей свободного пробега для уравнения Смолуховского
??1
???1
????
1
???1
????1
???1
????1
???1
????1
Т а б л и ц а
(1)
J?H (T )
?
tc , с Sd /Sv
T = 1, n1 (1) = 4.4444 · 10?1
4.46928 · 10?1 4.9 · 10?5
74
0.86
4.46922 · 10?1 2.5 · 10?6
87
280
?2
T = 5, n1 (5) = 8.1632 · 10
8.2380 · 10?2 1.6 · 10?5 113
2.0
?2
?6
8.2382 · 10
1.1 · 10
142 337.3
T = 10, n1 (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
T = 20, n1 (20) = 8.2644 · 10?3
8.3410 · 10?3 2.4 · 10?6 129
12.7
?3
?7
8.3405 · 10
2.7 · 10
161 808.4
2. Оценка ункционала
(12)
JH (T )
с использованием ценностного моделирования
номера пары и с приближенной ункцией ценности свободного пробега
??1
????1
????1
????1
????1
(12)
J?H (T )
?
tc , с Sd /Sv
T = 1, n1 (1) + n2 (1) = 5.9259 · 10?1
5.9574 · 10?1 6.9 · 10?5 108
0.31
T = 5, n1 (5) + n2 (5) = 1.3994 · 10?1
1.4161 · 10?1 7.5 · 10?6 197
9.50
T = 10, n1 (10) + n2 (10) = 5.0926 · 10?2
5.1607 · 10?2 4.2 · 10?6 197
12.82
T = 20, n1 (20) + n2 (20) = 1.5777 · 10?2
1.5991 · 10?2 1.7 · 10?6 206
26.85
????1
74
М. А. Коротченко
асчеты с ункцией ценности (2) дали неудовлетворительные результаты даже для
числа начальных частиц N0 порядка нескольких тысяч; это объясняется выполнени(1)
ем равенства JH (t) = n1 (T ) только асимптотически при N0 ? ?. Таким образом,
приближение (3) в сочетании с ценностным моделированием номера пары оказалось
достаточно эективным.
В заключение отметим, что построенные алгоритмы ценностного моделирования
можно использовать с целью улучшения алгоритмов решения реальных задач с переменными коэициентами Kij , для которых рассмотренные задачи дают приближенные решения.
Автор выражает благодарность чл.-корр. АН .А. Михайлову за ценные замечания
и предложения в процессе подготовки статьи.
Список литературы
[1? Михайлов .А., огазинский С.В. Весовые методы Монте-Карло для приближенного
решения нелинейного уравнения Больцмана // Сиб. мат. журн. 2002. Т. 43, ќ 3. С. 620628.
[2? Михайлов .А., Коротченко М.А., огазинский С.В. Ценностные алгоритмы статистического моделирования для нелинейных кинетических уравнений // Докл. АН. 2007.
Т. 415, ќ 1. С. 2630.
[3? Михайлов
.А.,
Войтишек
А.В. Численное статистическое
моделирование (метод
Монте-Карло). М.: Изд. центр Академия, 2006.
[4? Марчук .И., Агошков В.И., Шутяев В.П. Сопряженные уравнения и методы возмущений в нелинейных задачах математической изики. М.: Наука, 1993.
[5? Волощук В.М. Кинетическая теория коагуляции. Л.: идрометеоиздат, 1984.
Поступила в редакцию 20 евраля 2008 г.
Документ
Категория
Без категории
Просмотров
5
Размер файла
210 Кб
Теги
решение, уравнения, моделирование, алгоритм, статистический, смолуховского, ценностного
1/--страниц
Пожаловаться на содержимое документа