Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло тема диссертации и автореферата по ВАК РФ 01.01.07, кандидат физико-математических наук Роженко, Сергей Александрович

• Роженко, Сергей Александрович
• кандидат физико-математических наук
• 2013

Оглавление диссертации кандидат физико-математических наук Роженко, Сергей Александрович

Оглавление

Введение

Глава 1. Вводная информация

1.1. Процесс переноса частиц. Плотность столкновений и поток частиц

1.2. Общая схема моделирования процесса переноса. Оценка функционалов от потока частиц

1.3. Переходная плотность в цепи столкновений

1.4. Интегральное уравнение переноса (с обобщенным ядром)

1.5. Оценка по столкновениям

1.6. Рандомизированное ветвление траекторий

1.7. Коэффициент к^ размножения частиц

1.8. Метод расщепления

Глава 2. Модификация двухэтапных алгоритмов метода Монте-Карло с учетом свойств

симметрии задачи

2.1. Двухэтаипые алгоритмы метода Монте-Карло

2.2. Модификации на основе свойств симметрии первого этапа

2.3. Модификация расщепления траекторий в случае изотропного источника. Постановка задачи

2.4. Оценка трудоемкости на основе упрощённой модели

2.5. Выборка "вращения" по важности

2.6. Описание алгоритма численного моделирования

2.7. Анализ результатов

Глава 3. Модификация метода расщепления для решения задачи с ветвлением траектории

3.1. Постановка задачи

3.2. Оценки на основе упрощенной модели типа процесса Гальтона-Ватсопа в шаре

3.3. Упрощенная модель для всей задачи, "ветвление веса"

3.4. Численные эксперименты

3.5. Дополнительные замечания

Глава 4. Минимаксные параметрические весовые оценки в методе Монте-Карло

4.1. Постановка задачи

4.2. Вспомогательные утверждения

4.3. Минимаксная оптимизация стандартного МПТ

4.4. Минимаксный алгоритм с ветвлением траекторий

4.5. Оптимизация МПТ при вариации индикатрисы

4.6. Параметрическая оценка погрешности транспортного приближения

Литература

Введение диссертации (часть автореферата) на тему «Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло»

Введение

Свойства симметрии и подобия изучаемых процессов позволяют строить эффективные вычислительные алгоритмы с использованием соответствующих аналитических преобразований осреднения, поворота, сдвига и т.п. В методе Моите-Карло это связано с использованием условных математических ожиданий или в рандомизированном варианте—алгоритмов расщепления. Разработка таких приемов восходит к известной работе [22]. Построение более сложных весовых алгоритмов в значительной степени связано с геометрическими соображениями о подобии траекторий моделируемого процесса для различных вариантов физической модели (см., напр., [3,7,18,22,27]).

Свойства симметрии и подобия можно эффективно использовать в "двух-этапных алгоритмах" метода Монте-Карло. В литературе (см., напр., [9]) рассматриваются два варианта стандартного двухэтапного алгоритма: метод математических ожиданий и метод расщепления. На практике чаще всего используется метод расщепления. В диссертации разработаны модификации двухэтапных алгоритмов метода Монте-Карло с учетом свойства симметрии, то есть инвариантности, первого этапа относительно некоторого начального векторного параметра моделируемой траектории. Предлагаемая модификация состоит в формальном переносе моделирования указанного параметра на второй этап алгоритма. В "методе расщепления" это означает рандомизацию начальных точек вспомогательных траекторий. Показано, что такую рандомизацию можно улучшить, фактически применяя принцип Беллмаиа.

На модельных задачах радиационного контроля в диссертации разработан содержательный вариант учёта симметрии на первом этапе моделирования траектории для оптимизации соответствующего алгоритма расщепления. При этом оптимальные параметры алгоритма определяются на основе специально построенных упрощённых моделей задачи.

Можно полагать, что в главах 2, 3 диссертации рассматриваются варианты модельной задачи радиационного контроля, целью которой является определение типа или параметров "внутренней" среды по показаниям детектора частиц, которые сравниваются с заданными "эталонными значениями". В частности, практически важной является задача определения кон-

центрации активного изотопа, в котором реализуется процесс деления. Соответствующий параметрический анализ эффективно осуществляется весовым методом, причем для ограничения дисперсии в некоторых случаях оценки следует строить на ветвящейся цепи Маркова. Другой важной задачей является оценка величины полного сечения для внутренней среды. При этом для параметрического анализа показаний детектора естественно применять так называемый "метод подобных траекторий", который детально изучается и разрабатывается в главе 4.

Метод подобных траекторий (МПТ) позволяет проводить эффективный параметрический анализ исследуемых функционалов, необходимый, в частности, для решения обратных задач восстановления значений параметров по экспериментальным наблюдениям. Решение задач теории переноса излучения с помощью МПТ реализуется методом Монте-Карло путем численно-статистического моделирования траекторий частиц — "квантов излучения" — соответственно заданной радиационной модели среды. С помощью вспомогательного случайного веса при этом получаются оценки исследуемых функционалов для различных значений параметров модели, например, плотности среды.

Стандартный вариант МПТ связан с вариацией плотности или линейного размера среды и даёт возможность эффективной оценки, например, вероятности прохождения кванта через среду в зависимости от этих параметров [3,7,18,27]. Случайные "физические" траектории для различных значений модельной плотности при этом отличаются лишь длинами пробегов — этим и объясняется термин МПТ. В диссертации рассмотрен и соответствующий невесовой алгоритм — "геометрический МПТ". Изучается также весовой алгоритм для случая, когда варьируется параметр индикатрисы, то есть плотности распределения косинуса угла рассеяния частицы. Такой алгоритм можно рассматривать как вариант МПТ, так как соответствующие физические траектории отличаются лишь углами рассеяния.

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

обеспечить конечность дисперсий оценок во всем интервале параметров и тем самым расширить его при использовании МПТ. В диссертации теоретически и численно исследованы различные варианты этой задачи. Существенно уточнены ранее известные и получены новые утверждения о минимаксных алгоритмах МПТ.

Основные цели диссертационной работы:

• Разработка и обоснование модифицированного двухэтапного алгоритма расщепления с учетом свойства симметрии, т.е. инвариантности, первого этапа относительно некоторого начального векторного параметра моделируемой траектории. Сравнение традиционного и модифицированного методов расщепления на основе упрощенной модели и численных экспериментов.

• Разработка двухэтапного алгоритма с переходом от весового моделирования к невесовому в задаче с делением частиц. Исследование его эффективности в сравнении с весовым и аналоговым алгоритмами на основе упрощенной модели и численных экспериментов.

• Исследование задачи минимаксной оптимизации параметрического максимума среднеквадратической погрешности оценок алгоритмов МПТ для различных семейств плотностей. Сравнение соответствующих решений задачи в случае стандартного МПТ и в случае вариации параметра индикатрисы.

• Получение параметрических оценок погрешности транспортного приближения путем сравнения результатов расчетов с помощью МПТ с результатами расчетов в транспортном приближении.

Диссертация состоит из четырех глав, заключения и списка литературы.

В главе 1 приводится вводная информация но теме исследования, взятая из [4,9,10], за исключением описанного в и. 1.2 простого способа моделирования случайной величины щ с распределением

Р{щ = И) = 1 - {у - И), Р(иг = [и] + 1) = и- И с использованием новой формулы

VI = [У а].

Здесь и далее а—равномерно распределённая на промежутке [0, 1) случайная величина.

Глава 2 посвящена модификации двухэтапных алгоритмов метода Монте-Карло с учетом свойства симметрии, то есть инвариантности, первого этапа относительно некоторого начального векторного параметра моделируемой траектории. Предлагаемая модификация состоит в формальном переносе моделирования указанного параметра на второй этап алгоритма.

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

где Р\(-) — вероятностная мера в II, I и) — вероятностная мера в V при и £ и. В методе Монте-Карло

Случайный вектор (£, г/) численно моделируется в два этапа. На первом из них реализуется £ соответственно мере Р\(-), а на втором — г] соответственно условной мере Р2(- | £). Искомое приближённое значение величины 3 строится путем осреднения выборочных значений (, получаемых в результате численного моделирования [9]. Известно, что трудоёмкость алгоритма метода Монте-Карло оценивается величиной Б = tD(, где — дисперсия, а £— средние арифметические затраты на одну реализацию алгоритма [9].

Метод расщепления состоит в том, что для каждого значения £ на втором этапе реализуется К условно-независимых значений 771, г]к соответственно условной мере Р2(- | £) и вычисляется значение случайной величины

В п. 2.2 предлагается модификация двухэтапных алгоритмов метода расщепления, заключающаяся в следующем. Пусть £ однозначно определяется вектором (£1,^2), £1 € ¿Л, £2 € 6^2, причём £1 и £2 независимы и реализация £2

3 = ЕС, С = <7й, »7), пеУ.

= пРичём ЕС(Х) = ЕС

к

сравнительно малотрудоёмка. Тогда эффективным может оказаться преобразование этапов моделирования, которое определяется заменой £ на £1, ц на вектор (£2>т?)> и) соответственно, и на СД, V на и2 х К. Характерной в этом плане является задача переноса частиц с рассеянием от источника, расположенного в центре сферически-однородного шара. Пусть £—внутренняя часть траектории частицы (до вылета из шара), ц — остальная часть траектории. Через обозначим случайные начальное направление частицы и угол поворота траектории вокруг соответствующего луча, а через ^ — внутреннюю часть траектории, "освобождённую" от этих начальных параметров (то есть, для некоторого фиксированного начального направления, например, вдоль какой-либо из координатных осей, с заданным углом поворота). Траектория £ определяется вектором (£1,^2) путем вращения траектории £]_ от фиксированных начальных параметров к случайным параметрам

Следующая теорема показывает целесообразность соответствующей модификации алгоритма расщепления.

Теорема 2.1. Пусть требуется вычислить величину Е#(£1,??); где £1, £2 и V — некоторые случайные величины, ад — вычисляемая по ним статистика.

Пусть и ~ случайные величины, моделируемые с использованием метода расщепления следующим образом:

Если затраты на моделирование £2 пренебрежимо малы по сравнению с затратами на моделирование £1 и г/, то трудоёмкость модифицированного алгоритма не превосходит трудоёмкость традиционного алгоритма при одинаковом значении параметра расщепления К. Это соотношение, очевидно, усиливается при оптимальных параметрах расщепления.

Дополнительно можно использовать весовую модификацию моделирования £2 — "выборку по важности". Известно, что оптимальная плотность вероятности для такой модификации, согласно стохастическому аналогу принципа Беллмана [9,23], определяется выражением

Тогда

> 0.

/ьЫ = с^иЕ(С2 I У)

где /^(^г) — плотность распределения ^2(^2), а — нормирующая константа. Соответствующая весовая модификация оценки такова:

С =Сс*ч/Е(С21^'

В п. 2.3 описана модельная задача, обладающая упомянутым свойством сферической симметрии на первом этапе моделирования. В центре однородного шара, расположенного в бесконечной цилиндрической полости, находится изотропный источник частиц. На некотором расстоянии от шара находится перегородка (защита), за которой расположен детектор частиц. Фиксируется среднее число частиц, попавших в детектор. Расщсплснис траекторий частиц при моделировании производится при пересечении частицей "сферы расщепления", центр которой совпадает с источником частиц, а радиус не превосходит радиуса однородного шара. В качестве £2 здесь выступает вектор (р, <£>), состоящий из координат точки расщепления и угла поворота направления частицы вокруг радиус-вектора р при вылете частицы со сферы расщепления. Модификация метода расщепления с "вращением" состоит в переносе моделирования £2 на второй этап.

В п. 2.4 построена упрощенная модель задачи, в которой £1 и £2 — бер-нул л невские случайные величины, определяющие вылет частицы из шара и попадание частицы в цилиндрическую часть защиты около детектора соответственно. Полученные формулы позволяют получить априорную оценку оптимального числа частиц К в методе расщепления.

Вращение частиц также позволяет использовать "выборку по важности". В п. 2.5 рассматривается простейший вариант выборки по важности: сфера расщепления разбивается на две области плоскостью, перпендикулярной к оси симметрии. Понятно, что из ближайшей к детектору области \¥1 частицы попадут в детектор с большей вероятностью, чем из области И^- Обозначим: ¿1, — площади областей и \¥2; <71, ~~ модифицирующие "ценностные множители". Моделируется случайная величина

С' = УС(1)- + (1-У)С(2)~,

Ч1 42

где

, Р = Я1в1/а,

а = qlSl + д2.ч2, 7 = <

Р = «Дог/а,

с(г) — случайное значение £ при условии, что частица стартует на втором этапе из области г = 1,2.

Оценка всегда остается несмещенной и совпадает с £, когда дх и равны. Оптимальные весовые множители определяются по формулам

<Й = \/Е(С21 6 € ии Й = \/Е(С2 | & е \¥2).

Таким образом, и в дискретном варианте оптимизации "вращения" реализуется стохастический вариант принципа Беллмана [23].

В п. 2.6 приводится описание алгоритма численного моделирования, а в заключительном п. 2.7 приводятся результаты численных расчетов и выполнен анализ полученных результатов. Выигрыш в трудоёмкости в рассматриваемом примере при использовании модифицированного расщеплеиия с выборкой по важности по сравнению со стандартным методом расщепления при оптимальных значениях параметров составил 48 раз. Сравнение оптимальных значений К* для обычного и модифицированного алгоритмов, полученных для упрощенной модели, с расчетными значениями К, дающими наименьшую трудоёмкость, показало, что построенную в работе упрощенную модель целесообразно использовать для исследования и оптимизации расщепления.

В главе 3 рассмотрена модифицированная модельная задача, в которой объемлющий шар разбит па две области — внутренний шар, и внешний слой. Во внутреннем шаре происходит деление частиц, т.е. в точке столкновения с некоторой вероятностью происходит ветвление траектории.

В п. 3.1 приводится постановка задачи, а также описаны варианты моделирования с ветвлением траекторий и весового моделирования, при котором ветвление заменяется увеличением "веса" частицы. Предлагается модификация, комбинирующая весовой алгоритм внутри шара с расщеплением, соответствующим весу частицы, на границе шара.

В п. 3.2 рассматривается упрощенная модель шара, заключающаяся в замене двухслойного шара на модельный однородный при помощи простейшей "гомогенизации", а затем перехода к модели типа процесса Гальтона-Ватсона.

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

Для упрощенной модели получены оценки дисперсии числа вылетевших из шара частиц при аналоговом и весовом моделировании. Параметры модели выражены через сечения сг8, а/, ас и критическое значение и* коэффициента размножения, которое можно оценить с помощью улучшенного диффузионного приближения [16].

В п. 3.3 получены оценки оптимального параметра расщепления К* и трудоёмкости 5* в упрощённой модели для вариантов аналогового моделирования, весового моделирования, а также для нового варианта весового моделирования с "ветвлением веса". Идея последнего варианта заключается в следующем: при вылете частицы с весом -д из объемлющего шара в момент расщепления испускается частиц с весом 1/К, где

Показано, что весовое моделирование с ветвлением веса в методе расщепления выигрывает по трудоемкости у обычного весового моделирования при прочих равных условиях. При некоторых значениях параметров моделирования выигрыш может быть весьма существенным.

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

М [КЩ, Р = 1-{КЩ, к | [Кд] + 1, Р = {Кд}.

кости аналогового моделирования.

Использование "вращения" и выборки по важности при расщеплении существенно снижают трудоемкость алгоритма по сравнению с традиционным алгоритмом расщепления. При увеличении радиуса шара эффективность "вращения" возрастает.

Сравнение результатов расчётов со значениями, полученными для упрощенной модели, показало, что построенную в работе упрощенную модель целесообразно использовать для исследования и оптимизации модифицированного расщепления.

Глава 4 посвящена теоретическому и численному решению задачи выбора модели для минимизации параметрического максимума среднеквадрати-ческой погрешности весовых оценок функционалов в алгоритмах МПТ.

В п. 4.1 ставится задача о минимаксной оптимизации алгоритмов МПТ. Рассматривается односкоростной процесс переноса частиц, вероятностная модель которого определяется плотностью ае~ах (а > 0, х > 0) распределения случайной длины \ "свободного пробега", плотностью w(y) ( — 1 < у < 1) распределения косинуса /х угла рассеяния и вероятностью "выживания" q в точке "столкновения" (см., напр., [6]). Функция w(y) обычно называется индикатрисой рассеяния. Для решения задач теории переноса методом Монте-Карло численно моделируется цепь Маркова столкновений частицы с элементами вещества в соответствии со вспомогательной радиационной моделью и строится "оценка по столкновениям". При этом используется вспомогательный вес, который после очередного случайного перехода домножается на величину

q ае ах w(ß)

Ф) р(х) гО)'

где р{х) — моделируемая плотность распределения х-> Г('У) — моделируемая индикатриса, a г/о— моделируемая вероятность выживания.

В работе используется значение q0 = 1, причем обрыв траектории реализуется вследствие вылета из среды. В рассмотренных условиях средний квадрат оценки но столкновениям (для функционалов от плотности "физических" столкновений) определяется интегральным уравнением 2-го рода с оператором Кр, спектральный радиус которого в случае бесконечной среды равен

Р(кр) = Г ^J^ dx f1 ^ dy, (1)

а в случае конечной среды мажорируется этой величиной [8].

Из (1) следует, что для глобальной оптимизации параметрической весовой оценки по столкновениям целесообразно определять моделируемую плотность путем решения минимаксной задачи вида

9н(х)

р* = are: min max / —;—^ dx, (2)

P€P ß!<ß<ß2 J0 p(x)

где P—некоторое семейство непрерывных положительных плотностей, ß — параметр, gß е Р, ß\ < ß < /32. Обозначим

Лр\9) = ГЩ**-

Jo Pi?)

Решение задачи (2) при определенных условиях совпадает с решением более простой задачи

р* — arg min max J(p; gi), (3)

per 7=1,2

гAegt=9ßt, ¿ = 1,2.

В п. 4.2 рассмотрены варианты решения задачи оптимизации (3) для различных семейств плотностей.

Следующая теорема устанавливает решение задачи (3) на множество Pq всех непрерывных положительных плотностей.

Теорема 4.1. При д\ ф д2 решение минимаксной задачи (3) определяется выражением

р5Сг) = С(Л*)у Х*д1(х) + (1 - А*)д|(:г),

где величина А* является единственным решением относительно А Е (0,1) уравнения

Лр*91) = Лр*92). 0<А<1. (4)

Рассмотрим теперь семейство Р\ удобно моделируемых методом "суперпозиции" плотностей вида

Рг{х) = Хдг{х) + (1 - Х)д2(х), 0 < А < 1, в предположении, что д\ ф д2.

Теорема 4.2. Для семейства Р\ решение минимаксной задачи (3) определяется выражением

Л{х) = (5)

Теоремы 4.1 и 4.2 были приведены в [8,20] только для экспоненциальной плотности без достаточно детализированных доказательств.

В п. 4.3 исследуется минимаксная оптимизация стандартного МПТ. Коэффициент а рассматривается в интервале 0 < о-! < а < а 2 <оо. Очевидным образом задача сводится к интервалу 1 < а < I = а2/аг. Для семейства экспонент Р2 = {5е_®х, 1 < 5 < Ь = сг2/<Т1}, соответствующего простейшему варианту стандартного МПТ, плотность р\(х) — в*е~8*х будет решением

21

задачи (2) при я* = -—- е (1,шт{2^}).

I 1

Отметим, что плотность р\ неэффективна по сравнению с и р\. Так, например,

10) = 1.544, в{р\\ 10) = 1.557, 10) = 3.025.

Здесь — значение величины минимакса на оптимальном решении за-

дачи (3) для Еи = {М}.

Хорошее совпадение величин Ю) и С(р{; 10) объясняется следую-

щим утверждением.

Лемма 4.5. Имеет место предельное соотношение

у*>0.

р*(ж) г-оо

Таким образом, для оптимизации параметрических оценок в задачах теории переноса целесообразно использовать вспомогательные плотности вида р*, которые легко моделируются модифицированным методом суперпозиции:

если а <1/2 то £ := '01 (2а) иначе £ := ф2(2(а — 1/2)),

где 'фг — моделирующая функция для плотности р«, г = 1,2.

На основе достаточно точных расчетов было получено, что решение задачи (3) перестаёт совпадать с решением задачи (2) в случае р* = р£ при Ь > ¿ц = 12.62, а в случае р* = р\ — при £ > ^ = 17.59.

Следующее утверждение даёт условия, гарантирующие, что решение задачи (3) будет решением задачи (2).

Лемма 4.6. Пусть плотность р* {х) минимаксна для Si,2 = причем

w пусть

d2J(p*\ (je~ax) f°° (4а2х2 - 8ах + 2)е'2ах ,

min ----J- = / ^-—-J--dx > 0.

1 <a<t da2 Jq p*(x)

Тогда p* минимаксна для E = [1, t).

С помощью достаточно точных расчетов было получено, что плотности Pq, pi заведомо минимаксны соответственно для о £ [1, 2] и a £ [1, 8].

В п. 4.4 изучается минимаксный алгоритм, в котором из точки рассеяния "испускается" в среднем v условно-независимых частиц, "вес" которых домно-жается на q/v. Для рассматриваемой модели среды при этом имеем [10]:

р(Кр(и)) = = plp-aer").

Здесь для ограничения трудоёмкости моделирования следует использовать V < v*(p), где v*(p) соответствует критическому варианту моделируемой ветвящейся цепи Маркова.

Для приближенной оценки величины v*(pi) можно использовать среду с фиксированным значением сг^ = 2/(1 + Известно, что вероятность вылета частицы при рандомизации радиационной модели увеличивается (см., напр., [8]), поэтому можно ожидать, что соответствующее > v*{p\).

Для проверки данного неравенства был использован шар радиуса 7.72043 с параметрами среды

(т = 1, (Tf = 0.03, (ts = 0.97, и = 2.5,

и изотропным рассеянием, для которого известно v* = 1.045. В численных экспериментах с использованием а1 = 1 /л/2 и соответствующего <т2 = 1/л/2 + 1 получена оценка v*(p\) = 1.049 > v* со среднеквадратической погрешностью 0.00025.

В п. 4.5 рассматривается МПТ для решения задач теории переноса с индикатрисой ws(y), — 1 < у < 1, достаточно регулярно зависящей от некоторого параметра s, причем si < s < S2-

Для численного эксперимента используется стандартная индикатриса Хеньи-Гринстейна:

+ С3'

где й = Е//,— средний косинус угла рассеяния.

Численные расчеты были выполнены для аналогов рассмотренных в п. 4.3 плотностей и так же показали целесообразность использования плотности

г;(у) = „,,(,).„„(,), < = 1,2,

дающей решение задачи (3) в семействе плотностей Р\.

В п. 4.6 численно исследована параметрическая зависимость погрешности "транспортного приближения" для вероятностей прохождения, поглощения и альбедо частицы.

Транспортное приближение состоит в замене индикатрисы (6) на

и}(у) = ^- + 85(у- 1), -1<у<1. (7)

Подстановка равенства (7) в уравнение переноса приводит к следующей модификации параметров радиационной модели:

а а = (1 - дв)а, д д = ^-^ I (8)

1 — дй 2

В численных экспериментах рассматривался плоский слой 0 < г < Н вещества с индикатрисой вида (6). Для различных значений параметров с помощью МПТ оценивались вероятности прохождения Р1 и альбедо Ра кванта излучения для источника частиц, "падающих" извне перпендикулярно к поверхности .г = 0, а также вероятность Рс = 1 — Рг — Ра поглощения кванта внутри слоя. При этом различные значения Н были реализованы геометрическим вариантом МПТ —траектории моделировались для максимальной толщины, а для меньших слоев учитывались вплоть до вылета частицы за их пределы. Было введено описанное в п. 4.4 ветвление траекторий, обеспечивающее конечность дисперсий оценок для всех значений параметров.

Для сравнения были вычислены оценки функционалов Рг, Ра и Рс для транспортного приближения по формулам

Р1 = 1 [ Г1 (1у Ах + е~дН,

2 Л Jo

= Ъ Г Г ¿у <1г,

^ ¿о Jo

Рс = 1 - А - Ра,

где <р(г) — функция плотности столкновений, являющаяся решением "плоского" интегрального уравнения Пайерлса для параметров (8):

ГН

ф) = ад / ф') К(а\г - г'\) <1х' + а еГа\ (9)

Л

1 Г°° е~4 с ядром К(х) — - / — сИ.

Уравнение (9) решалось путём разбиения отрезка [О, на М равных частей со следующей дискретизацией: ^ = </>(^г), 2г — — 1 /2)/г., /г = Н/М,

а/ г*+у2

Щ = / А'(ф7; - г'|) + (те^, г = 1,..., Л/. (10)

•? = 1 2,-/1/2

Следующая теорема показывает, что для решения системы уравнений (10) можно применять высокоточные прямые методы решения систем линейных уравнений, например метод квадратного корня или ЬИЬТ факторизации (см., напр., [2]).

Теорема 4.3. Матрица системы линейных уравнений (10) симметрична. положительно определена и имеет строгое диагональное преобладание. Её спектр содержится в интервале (1 — д, 1 + д).

Отметим, что число обусловленности матрицы системы (10) оценивается величиной т—4, не зависящей от М и Н.

С помощью замены переменных во внутреннем интеграле получены следующие формулы для приближенного вычисления Д и Ра:

- м *'+А/2

У Е2(а(Н - г)) ¿¡г + е-*н, (И)

-7"1 г^-—/г/2

- м г'+л/2

I Е2{аг)<1г. (12)

Здесь Е-2(х) — интегральная показательная функция [19, п. 5.1.45]:

roo -t

Еп{х)=хп-Ч —dt, n = 1,2,...

С учетом тождества E'n{x) = —i?n_i(.x') и равенства i^(x') = 0.5£i(x-), интегралы в (10) и в (11), (12) эффективно вычисляются соответственно через значения интегральных показательных функций Е2 и в узлах сетки.

Заключение содержит перечень основных результатов диссертационной работы.

Основные результаты диссертации опубликованы в работах [11-13,24-26] в рецензируемых изданиях, рекомендуемых ВАК.

Автор диссертации принял активное участие во всех совместных публикациях, в особенности на стадии разработки алгоритмов, вычислительных программ и проведения расчётов. Поэтому его участие можно заведомо оцепить, как 50%.

На защиту выносятся следующие основные результаты:

1. Построена модификация алгоритма расщепления с учетом свойства симметрии, то есть инвариантности, первого этапа относительно некоторого начального векторного параметра моделируемой траектории. Предложенная модификация состоит в формальном переносе моделирования указанного параметра на второй этап алгоритма и эквивалентна рандомизации начальных точек вспомогательных траекторий. Показано, что такую рандомизацию можно улучшить, фактически применяя принцип Беллмана. Эффективность описанной модификации численно проверена для двух вариантов модельной задачи, первый этап которой обладает свойством сферической симметрии. При этом были использованы оценки оптимального параметра расщепления, построенные с использованием упрощённых моделей.

2. Исследована задача минимаксной оптимизации весовых оценок в методе подобных траекторий для различных семейств плотностей распределений, определяющих соответствующие алгоритмы. Уточнены ранее известные и получены новые утверждения о минимаксных алгоритмах МПТ. Для стандартного МПТ и для варианта с вариацией параметра индикатрисы показана целесообразность использования плотности, оптимальной на семействе плотностей типа "смеси", легко моделируемых методом суперпозиции.

3. Получены параметрические оценки погрешности вероятностей прохождения, поглощения и альбедо частицы в транспортном приближении для

задачи теории переноса в слое вещества с индикатрисой Хсньи-Гринстейна. С использованием МПТ значения вероятностей для различных вариантов радиационной модели были получены на одном ансамбле траекторий, а транспортное приближение было реализовано путём решения "плоского" интегрального уравнения Пайерлса.

Результаты, изложенные в диссертационной работе, неоднократно докладывались и обсуждались на семинаре "Методы Монте-Карло в вычислительной математике и математической физике" ИВМиМГ СО РАН, докладывались на Международных научных студенческих конференциях "Студент и научно-технический прогресс" (2008, 2009 гг.), на Конференциях молодых учёных ИВМиМГ СО РАН (2009, 2011, 2013 гг.), на шестом Санкт-Петербургском международном семинаре по стохастическому моделированию (2009 г.), на Всероссийской конференции по вычислительной математике КВМ-2011.

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

Заключение диссертации по теме «Вычислительная математика», Роженко, Сергей Александрович

Заключение

В диссертационной работе получены следующие результаты:

1. Построена модификация алгоритма расщепления с учетом свойства симметрии, то есть инвариантности, первого этапа относительно некоторого начального векторного параметра моделируемой траектории. Предложенная модификация состоит в формальном переносе моделирования указанного параметра на второй этап алгоритма и эквивалентна рандомизации начальных точек вспомогательных траекторий. Показано, что такую рандомизацию можно улучшить, фактически применяя принцип Беллмана. Эффективность описанной модификации численно проверена для двух вариантов модельной задачи, первый этап которой обладает свойством сферической симметрии. При этом были использованы оценки оптимального параметра расщепления, построенные с использованием упрощённых моделей.

2. Исследована задача минимаксной оптимизации весовых оценок в методе подобных траекторий для различных семейств плотностей распределений, определяющих соответствующие алгоритмы. Уточнены ранее известные и получены новые утверждения о минимаксных алгоритмах МПТ. Для стандартного МПТ и для варианта с вариацией параметра индикатрисы показана целесообразность использования плотности, оптимальной на семействе плотностей типа "смеси", легко моделируемых методом суперпозиции.

3. Получены параметрические оценки погрешности вероятностей прохождения, поглощения и альбедо частицы в транспортном приближении для задачи теории переноса в слое вещества с индикатрисой Хеньи-Гринстейна. С использованием МПТ значения вероятностей для различных вариантов радиационной модели были получены на одном ансамбле траекторий, а транспортное приближение было реализовано путём решения "плоского" интегрального уравнения Пайерлса. Для использованной при этом дискретизации доказано, что матрица получающейся системы линейных уравнений положительно определена и имеет строгое диагональное преобладание, а число обусловленности оценивается константой, не зависящей от шага дискретизации.

Список литературы диссертационного исследования кандидат физико-математических наук Роженко, Сергей Александрович, 2013 год

Литература

[1] Боровков A.A. Теория вероятностей. — М.: Наука, 1986.

[2] Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. — М.: Наука, 1984.

[3] Ермаков С. М. Метод Монте-Карло и смежные вопросы. — М.: Наука, 1974.

[4] Ермаков С.М., Михайлов Г. А. Курс статистического моделирования. — М.: Наука, 197G.

[5] Марчук Г.И. Методы расчета ядерных реакторов. — М.: Госатомиздат, 1961.

[6] Марчук Г.И., Михайлов Г.А., Назаралиев М.А. и др. Методы Монте-Карло в атмосферной оптике. — Новосибирск: Наука, 1976.

[7] Михайлов Г.А. Некоторые вопросы теории методов Монте-Карло. — Новосибирск: Наука, 1974.

[8] Михайлов Г.А. Оптимизация весовых методов Монте-Карло. — М.: Наука, 1987.

[9] Михайлов Г.А., Войтишек A.B. Численное статистическое моделирование. Методы Монте-Карло: Учеб. пособие. — М.: Изд. центр "Академия", 2006.

[10] Михайлов Г.А., Медведев И.Н. Оптимизация весовых алгоритмов статистического моделирования. — Новосибирск: Омега Принт, 2011.

[11] Михайлов Г.А., Роженко С.А. Модификация двухэтапных алгоритмов метода Монте-Карло на основе свойств симметрии первого этапа // ЖВМиМФ. — 2009. — Т. 49, № 11, — С. 2010-2019; G.A. Mikhailov, S.A. Rozhenko. Modification of two-step Monte Carlo algorithms based on the symmetry of the first step // Computational Mathematics and Mathematical Physics. —2009. —Vol. 49, № 11. —P. 1921-1929.

[12] Михайлов Г.А., Роженко С.А. Минимаксная оптимизация численно-статистического "метода подобных траекторий" // ДАН. — 2012. — Т. 446, № 1.— С.15-17.

[13] Михайлов Г.А., Роженко С.А. Минимаксные параметрические весовые оценки в методе Монте-Карло // ЖВМиМФ. —2013.— Т. 53, № 9 (в печати).

[14] Роженко С. А. Модификация двух-этаиных алгоритмов метода Монте-Карло с учетом свойств симметрии задачи: Магистерская диссертация / Министерство образования и науки Российской федерации. Федеральное агентство по образованию. Новосибирский государственный университет. Механико-математический факультет. Кафедра вычислительной математики. — Новосибирск, 2010.

[15] Роженко С.А. "Расщепление" и ветвление траекторий в методе Монте-Карло с учетом симметрии источника частиц // Тезисы докл. Всероссийской конференции по вычислительной математике КВМ-2011, 29 июня - 1 июля 2011, Новосибирск.— http://www.sbras.ru/ws/show_abstract.dhtml?ru+220+16150.

[16] Романов Ю.А. Точные решения односкоростного кинетического уравнения и их использование для рассчета диффузионных задач (усовершенствованный диффузионный метод) // Исследование критических параметров реакторных систем. — М.: Госатомиздат, 1960. — С. 3-26.

[17] Смирнов Н.В., Дунин-Барковский И.В. Краткий курс математической статистики для технических приложений. — М.: Физматгиз, 1959.

[18] Соболь И.М. Численные методы Монте-Карло. — М.: Наука, 1973.

[19] Справочник по специальным функциям с формулами, графиками и математическими таблицами / Пер. с англ. под ред. В.А. Диткина и J1.H. Кармазиной; Под ред. М. Абрамовица и И. Стиган. — М.: Наука, 1979.

[20] Старков А.В., Стронгина Т.А. Минимаксная оптимизация на классе плотностей типа смеси // Труды Вычислительного центра СО РАН. Серия Вычислительная математика.—1989.—Вып. 1. — С. 53-62.

[21] Boost С++ Libraries. — http://www.boost.org.

[22] Kahn Н. Use of differential Monte Carlo sampling techniques // Sympos. On Monte Carlo methods / Ed. H.A. Meyer. —New York: Wiley, 1956. —P. 145-191.

[23] Mikhailov G.A. Recurrent formulae and the Bellman principle in the Monte Carlo method // Rus. J. Num. Anal, and Math. Modell. —1994, —Vol. 9, № 3.-P. 281-300.

[24] Rozhenko S.A. Trajectories splitting algorithm modification taking account of system simmetry / S.M. Ermakov, V.B. Melas, and A.N. Pepelyshev, eds. // Proc. 6th St. Petersburg workshop on simulation, St. Petersburg, June 28 - July 4, 2009. — St. Petcrburg VVM со. Ltd., 2009. —Vol. I. —P. 181-186.

[25] Rozhenko S.A. Splitting of trajectories in the Monte Carlo method taking into account the symmetry of the source /'/ RJNAMM. —2012, —Vol. 27, № 2. —P. 175-189.

[26] Rozhenko S.A., Mikhailov G.A. Minimax parametric optimization of numerical-statistical 'method of similar trajectories' for solution of radiation transfer theory problems // RJNAMM.-2013.-Vol. 28, № 2. —P. 201-212.

[27] Spanier J., Gelbard E. Monte Carlo principles and neutron transport problems.— Reading: Addison-Wesley, 1969; Дж. Спанье, 3. Гелбард. Метод Монте-Карло и задачи переноса нейтронов. — М.: Атомизда