Вычислительные технологии
Том 17, № 2, 2012
К обоснованию модификации
>к
метода максимального сечения*
В. С. Антюфеев Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: ant@osmf.sscc.ru
Рассматривается метод максимального сечения, используемый при решении задач переноса излучения методом Монте-Карло в неоднородной среде. Кратко описаны стандартный и модифицированный способы моделирования длины пробега. Представлено вероятностное доказательство модификации этого метода.
Ключевые слова: метод Монте-Карло, моделирование, вероятностное распределение, длина свободного пробега, фиктивное рассеяние, скачкообразный процесс.
Введение
Расчет характеристик поля рассеянного излучения методом Монте-Карло связан с моделированием траекторий частиц, и это занимает основную часть времени расчета. Метод максимального сечения предназначен для ускорения моделирования траекторий.
Уравнение переноса излучения и его параметры рассмотрены в работах [1-3]. Приведем некоторые используемые ниже термины и обозначения. Траектория частицы — это ломаная линия, случайные точки Гк столкновения частицы с частицами среды — ее узлы. Плотность распределения случайной длины ¿к = \г'к-\ — г к | пробега частицы
где ш — единичный вектор направления пробега, и(г) > 0 — переменное сечение ослабления среды. Величина I = ¿к моделируется за два шага.
Сначала моделируется оптическая длина свободного пробега т = т([гк-1, гк]):
Здесь т имеет экспоненциальное распределение, ее моделируют по формуле т = — 1п(а) [1-4] (а — случайная величина, равномерно распределенная на [0, 1]). Затем (2) рассматривается как уравнение относительно I и вычисляется ¿.
* Работа выполнена при финансовой поддержке Междисциплинарного интеграционного проекта СО РАН № 47.
(1)
(2)
о
Для удобства моделирования траекторий а (г) приближают кусочно-постоянной функцией. Область разбивается на подобласти $1, Б2, ... , в каждой из которых значения а (г) считаются постоянными: а (г) = ак при г Е Бк. Тогда (2) принимает вид
т
т = ^2 аи 4, (3)
к= 1
и искомая длина £ = й1 + ... +
Первый шаг — моделирование т — не вызывает затруднений. Однако последующее вычисление £ из предыдущей формулы и из (3) представляет собой трудоемкую процедуру: необходимо заранее вычислить точки пересечения луча {г = гк + £ шк+1} с границами областей Б к и длины отрезков 4, что значительно увеличивает трудоемкость моделирования.
1. Метод максимального сечения и его модификация
В работе [4] предложен метод максимального сечения моделирования длины пробега £ в средах с переменным сечением а (г). Зададим постоянную мажоранту Е : Е > а (г) для а (г). Следующий алгоритм задает моделирование £:
1) полагаем £ = 0;
2) моделируем случайную величину Д£ с экспоненциальной плотностью распределения
/Е(*) = Еехр(-Е£) , £ > 0 (4)
(Д£ — длина свободного пробега в среде с постоянным сечением Е);
3) увеличиваем длину £: £ ^ £ + Д£;
4) моделируем а Е [0, 1] и проверяем, выполняется ли неравенство
а < а(г + (£ + Д£)и)/Е :
если последнее не выполняется, возвращаемся к п. 2, если выполняется — заканчиваем моделирование.
Таким образом, £ представляется в виде суммы
£ = Д£1 + Д£2 + ... + Д£М. (5)
Здесь количество ^ случайных слагаемых Д£ также случайное:
Обычно точки
^ = шт | п : а < Е ^г + ^^ Д£к ^ Е
+ ^Д£и ^к, п = 1, ... ,(^- 1),
8п г к=1
интерпретируют как точки "фиктивного рассеяния" (последняя ^-я точка — точка "истинного рассеяния"). Моделирование слагаемых Д£к в среде с постоянным сечением не является трудоемким, если число фиктивных точек рассеяния невелико. В данном случае не надо вычислять точки пересечения луча пробега с границами областей Бк
и числа ^. В этом преимущество метода максимального сечения перед стандартным методом моделирования. Доказано [4], что моделирование длины пробега методом максимального сечения не изменяет ее вероятностное распределение.
В ряде задач значения функции а (г) сильно варьируют, и постоянная мажоранта Е в некоторых областях значительно превышает функцию а(г). В результате количество фиктивных столкновений и слагаемых в сумме (5) оказывается очень большим, что снижает эффективность метода максимального сечения. В этом случае для увеличения эффективности лучше использовать не постоянную, а кусочно-постоянную мажоранту сечения ослабления, которая лучше аппроксимирует функцию а (г).
В [5] предложена соответствующая модификация метода максимального сечения. Здесь слагаемые являются зависимыми случайными величинами, и их вероятностные распределения различны. Пусть уже смоделированы слагаемые Д^1, Д^2, ... , Д^к-1 и вычислена очередная точка фиктивного рассеяния Гк-1; тогда плотность следующей длины Д4
/s(i) = E(rfc-1 + tuk) exp
- E(rk-i + uuk) du
t > 0. (6)
Краткое обоснование этой модификации состоит в преобразовании уравнения переноса излучения. Точнее говоря, в [5] доказано, что моделирование длины пробега t методом максимального сечения для E(t) = const не изменяет функцию потока излучения, которая является решением уравнения переноса излучения.
В работе [8] приводится доказательство данного факта, не опирающееся на уравнение физического процесса. В следующем разделе дано еще одно доказательство, основанное на использовании уравнения Колмогорова — Феллера для функции распределения чисто разрывного процесса.
2. Вычисление вероятностной плотности распределения длины пробега
В области, где происходит перенос частиц, зафиксируем точку столкновения r = Гк-1, направление пробега ш = Uk после столкновения в r и выделим луч {r + tu : t > 0}. Последующие выкладки относятся к этому лучу. Нас будут интересовать лишь одномерные вероятностные распределения. Вместо точек (r + tu) Е R3 рассмотрим точки t на луче и соответственно функции одной переменной t.
Итак, пусть на луче {t > 0} заданы функция сечения a(t) и ее мажоранта E(t). Отметим на данном луче последовательность случайных точек
sb s2, ... , , Sk =
i=1
Ё Ati, (7)
соответствующих точкам рассеяния частицы на одном звене траектории: точки вм истинному рассеянию, остальные точки — фиктивному рассеянию (рис. 1).
0 Si S2 s
t
Рис. 1. Звено траектории: точки истинного (•) и фиктивного (о) рассеяния
Теорема. Вероятностная плотность случайной величины имеет вид
¡а (г) = а (г) ехр
— а(м)
г > 0.
Доказательство. Запишем уравнение Колмогорова — Феллера [6, 7] для функции распределения чисто разрывного случайного процесса £(г), следуя обозначениям [6]:
У сс
(' ' 'У) ^ р(т, г) 4^(г,х; т,у) + [ р(т, г)С(т, г, у) 4^(г,х; т, г). (8)
дт
— с —с
Здесь
1 — — моменты скачков процесса, х, у, г — значения выборочных функций;
2 — ^(г, ж; т, у) (функция распределения процесса £(г)) задает вероятность события {£(т) < У | £(г) = ж}, г < т; 4^(г, ж; т,у) — вероятностная мера Стилтьеса;
3 — С(г,х,у) = С(г,х ^ г < у) — условная функция распределения величины £(т) после скачка при условии, что в момент г произошел скачок, а непосредственно перед скачком £ (г — 0) = ж;
вероятность скачка в промежутке (г, г + Дг)
4 — р(г,х) = 11т ---; £(т) = х до скачка.
Ниже (см. пункты 1-3) определим функциональные характеристики процесса £ (г) таким образом, чтобы вероятностные распределения соответствующих членов последовательностей {¿1, , ... } и |з1, ... } совпали. Кроме того, фиксируем значения некоторых переменных в (8). Это позволит свести исходную задачу к решению обыкновенного дифференциального уравнения.
1. Пусть каждая выборочная функция рассматриваемого процесса (рис. 2) принимает лишь значения 0 и 1. В таком случае траекторию процесса можно записать в следующем виде:
'1, 0 < г < и ,
£(г) • 0, ( > 1„.
г
■о
¿1 ¿2
Рис. 2. Траектория чисто разрывного процесса £(£): фиктивных скачков
— точка истинного скачка, о — точки
1
г
0
V
Эта функция испытывает скачок (разрыв) лишь в точке ¿^. Однако данный процесс испытывает (фиктивные) скачки в точках слева и справа от .
Итак, траектории процесса £(¿) удовлетворяют двум условиям:
— начальная точка процесса £(£) фиксирована: £ = 0, £(0) = 1 (отсюда следует, что все рассматриваемые далее значения т положительны);
— процесс имеет лишь один истинный скачок; нас интересует вероятностное распределение момента этого скачка; после него процесс перейдет в состояние 0, что эквивалентно утверждению {£(£) < 0 при £ > }.
Данные условия позволяют упростить вид уравнения (8). Фиксируем значения переменных ¿, х, у в соответствующих частях (8) и перепишем его в виде
0 те
дР Г С
— (0, 1; т, 0) = - р(т, г) 4Р(0, 1; т, г) + р(т, г) С(т, г, 0) 4Р(0, 1; т, г).
Множество значений процесса состоит из точек 0 и 1. Следовательно, при фиксированном т > 0 вероятностная мера Р(0, 1; т, г) сосредоточена на множестве {г = 0, 1}. В силу этого оба интеграла по мере Р(•) можно переписать в виде сумм, причем первая сумма содержит лишь одно слагаемое. Обозначим как т(т, г) вероятностную меру точек {г}, (г = 0, 1). Теперь (8) принимает вид
дР
— (0, 1; т, 0) = - [р(т, 0) • т(т, 0)] +
+ [р(т, 1) С(т, 1, 0) • т(т, 1) + р(т, 0) С(т, 0, 0) • т(т, 0)]. (9)
Таким образом, уравнение (8) относительно функции четырех переменных преобразовано в уравнение относительно функции одной переменной т. Здесь следует также правильно задать функции р(-), С(-) в соответствии с целью, поставленной выше, и вычислить значения меры та(-) через значения функции Р.
2. Зададим функцию р(т, г), г = 1, 0, отвечающую за распределение случайных промежутков времени (¿и — 1) между моментами скачков процесса £(£) (фиктивных или истинных). Определим функцию р(-) так, чтобы распределения разностей (¿и — совпадали с распределениями соответствующих разностей в и — = .
Согласно [6], обозначим через рп(£,х, т) вероятность того, что, отправляясь из состояния х в момент ¿, до момента т процесс п раз изменит свое состояние. В работе [6] приведены формулы для всех п > 0. В частности, вероятность того, что от момента £ = ¿и до момента т процесс не совершит ни одного скачка (п = 0), равна
р0(£, х, т) = ехр
— р(м, х)
а вероятность хотя бы одного скачка в промежутке (¿,т)
Р^ (т) = 1 — ехр
— р(м, х)
т
т
Последнее утверждение, очевидно, эквивалентно следующему: вероятность того, что после момента £ скачок произойдет раньше момента т, равна Р^(т). Другими словами, Р^ (т) представляет собой функцию распределения первого (после момента ¿) скачка.
С другой стороны, известно [1-4], что функция распределения случайной длины пробега Д£ для частицы, стартовавшей в точке г, имеет вид
Ре(т ) = 1 — ехр
— Е(м)
Отличие этой формулы от предыдущей лишь в том, что под знаком интеграла стоит другая функция. Напомним, что мы хотим задать функцию р(г, м) так, чтобы распределения соответствующих случайных величин — 1) и (вк — вк_ 1) совпадали. Поэтому пусть
р(г, м) = Е(м), м = 0, 1. (10)
В данном случае функции распределения Р^ и Р^ совпадают. Следовательно, для всех к распределение величины г к совпадает с распределением величины вк. Итак, функция р(м, х) выбрана так, что условие 2 выполнено.
3. Для каждого скачка процесса функция С(т,г,у), г =1, 0, определяет вероятность события, будет ли этот скачок истинным. Определим С(-) так, чтобы распределение момента истинного скачка гv совпало с распределением длины истинного пробега частицы вм. Напомним, что согласно определению С(т,г,у) — это (условная) вероятность того, что после скачка процесс перейдет в состояние {< у}. Нас интересует лишь переход процесса в состояние 0, поэтому значение последнего аргумента в левой части (8) равно 0. Перепишем соответственно выражение для С(т,г,у) (здесь г = 0 или 1):
С(т, г, у) = С(т, г ^ 0) = С(т, г, 0).
Вероятность перехода (0 ^ 0) равна 1, так как процесс, попав в состояние 0, остается там навсегда. В то же время скачок (1 ^ 0) соответствует истинному рассеянию частицы. Однако акт истинного рассеяния происходит с вероятностью а(г)/Е(г) < 1 (см. п. 4 алгоритма моделирования длины Д£). Поэтому следует определить значение С(-) следующим образом:
( 1, г = 0, С(т,г, 0)= { ^ (11)
I Е(г), г = 1.
Тогда распределение величины гv совпадает с распределением вм. Таким образом, функция С(т, г, 0) выбрана так, что условие 3 выполнено.
Учитывая п. 3 и 4 для р, С, перепишем уравнение (9) в виде
дР
— (0, 1; т, 0) = — [Е(т) • т(т, 0)] +
Е(т) • ^ • т(т, 1) + Е(т) • 1 • т(т, 0)
Е(т)
или
или
дР
— (0, 1; т, 0) = — [Е(т) • т(т, 0)] + [а(т) • т(т, 1) + Е(т) • т(т, 0)] , дт
дР
—(0, 1; т, 0) = а(т) • т(т, 1). (12)
Согласно определению, Р(0, 1; т, 0) = Р^ < т), т.е. Р(0, 1; т, 0) — искомая функция распределения величины гv. Обозначим ее Ф(т): Ф(т) = Р(0, 1; т, 0). Мера т(т, г) сосредоточена на множестве {г = 0, 1} и является вероятностной. В силу этого
т(т, 0) + т(т, 1) = 1, или т(т, 1) = 1 — т(т, 0). (13)
т
По определению меры Стилтьеса:
т(т, 0) = lim ( F(0, 1; т, 0) - F(0, 1; т, 0 - Az)) = lim ( Ф(т) - 0) = Ф(т). (14)
Az^ü Az^ü
Пользуясь формулами (13), (14), перепишем уравнение (12):
d
— Ф(т) = а(т)(1 - Ф(т)). (15)
dT
Согласно определению процесса £ (t):
Ф(0) = 0.
Проинтегрировав обыкновенное дифференциальное уравнение (15) с этим начальным условием, находим:
Ф(т) = 1 - exp
— a(u) du
о
Отсюда искомая плотность распределения случайной величины tv имеет вид
d
fa(t) = -Ф(t) = a(t) exp
— J a(u)du о
т. е. совпадает с естественной плотностью распределения длины пробега, что и требовалось доказать.
Список литературы
[1] Spanier J., Gelbard E.M. Monte Carlo principles and Neutron Transport. Massachusetts: Addison-Wesley, Reading, 1969.
[2] Ермаков С.М., Михайлов Г.А. Курс статистического моделирования. М.: Наука, 1976.
[3] Марчук Г.И., Михайлов Г.А., Назаралиев М.А. и др. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976.
[4] Coleman W.A. Mathematical verification of a certain Monte Carlo sampling technique to radiation transport problems // Nucl. Sci. Eng. 1968. Vol. 32, No. 1. P. 76-81.
[5] Михайлов Г.А. Метод моделирования длины свободного пробега частицы // Атомная энергия. 1970. Т. 28, вып. 2. С. 175.
[6] Гнеденко Б.В. Курс теории вероятностей. М.: Наука, 1965.
[7] Feller W. An introduction to probability theory and its applications. New York; Chichester; Brisbane; Toronto: John Wiley & Sons, 1970.
[8] Аверина Т.А., Михайлов Г.А. Алгоритмы точного и приближенного статистического моделирования пуассоновских ансамблей // Журн. вычисл. математики и матем. физики. 2010. Т. 50. С. 1005-1016.
Т
Т
Поступила в редакцию 8 августа 2011 г., с доработки — 26 декабря 2011 г.