К обоснованию модификации метода максимального сечения Текст научной статьи по специальности «Математика»

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

Том 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 г.