Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло

Роженко, Сергей Александрович

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

Роженко, Сергей Александрович АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Новосибирск МЕСТО ЗАЩИТЫ

2013 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

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

Автореферат диссертации на тему "Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло"

На правах рукописи

Роженко Сергей Александрович

Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло

01.01.07— Вычислительная

математика

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

18 \:ж ¿'¡ш

005531525

НОВОСИБИРСК - 2013

005531525

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте вычислительной математики и математической геофизики Сибирского отделения Российской академии наук.

Научный руководитель:

член-корреспондент РАН, Михайлов Геннадий Алексеевич

Официальные оппоненты:

Симонов Николай Александрович, доктор физико-математических наук,

Baker Hughes, Новосибирский технологический центр, старший научный сотрудник.

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

Федеральное государственное бюджетное учреждение науки Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук, старший научный сотрудник.

Ведущая организация:

Федеральное государственное бюджетное учреждение науки Институт теоретической и прикладной механики им. С.А. Христиановича Сибирского отделения Российской академии наук.

Защита состоится 18 сентября 2013 г. в 16 часов 30 минут на заседании диссертационного совета Д 003.061.01 на базе Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (ИВМиМГ СО РАН) по адресу: 630090, г. Новосибирск, проспект Академика Лаврентьева, 6.

С диссертацией .можно ознакомиться в библиотеке ИВМиМГ СО

РАН.

Автореферат разослан 3 июля 2013 г.

Ученый секретарь диссертационного совета д.ф.-м.н.

С.В. Рогазинский

Общая характеристика работы

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

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

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

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

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

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

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

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

Научная новизна. Автором получены новые результаты в следующих направлениях:

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

• Показано, что построенную модификацию алгоритма расщепления можно улучшить, фактически применяя принцип Белл-мана.

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

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

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

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

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

Методология и методы исследования. В диссертации используется математический аппарат методов Монте-Карло, числен-

ные методы линейной алгебры, аппарат математического анализа, теории оптимизации и теории вероятности.

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

Положения, выносимые на защиту:

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

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

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

ло реализовано путём решения "плоского" интегрального уравнения Пайерлса.

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

Публикации. По теме диссертации опубликовано 5 работ и еще одна принята к печати. Все работы опубликованы в рецензируемых изданиях, рекомендуемых ВАК. Список работ автора по теме диссертации приведен в конце автореферата.

Работа состоит из введения, четырех глав, заключения и списка литературы из 27 наименований. Объём диссертации—84 с.

В главе 1 приводится вводная информация по теме исследования.

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

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

где Pi( ) — вероятностная мера в U, Р2(- | и) — вероятностная мера в V при и Є U. В методе Монте-Карло J = Е£, £ = q{^,rf), £ Є U, г] Є V. Случайный вектор (£, rj) численно моделируется в два этапа.

Содержание работы

На первом из них реализуется £ соответственно мере (-), а на втором— Г] соответственно условной мере /"2 (• | £). Искомое приближённое значение величины 3 строится путем осреднения выборочных значений получаемых в результате численного моделирования.

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

случайной величины — ~ Е^-1 г/д-). причём Е ((к'> =

1л. к— -I

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

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

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

Пусть и — случайные величины, вычисляемые методом

расщепления: с[К) = ± £ д(6,6,Л ' = £ Ё Тогда к=1 к=1

В([К) -ВС(2К) = ^ (Е,(9(6,6,77) I 6,6) | б) > 0.

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

Дополнительно можно использовать весовую модификацию моделирования £2 — "выборку по важности". Оптимальная плотность вероятности для такой модификации, согласно стохастическому аналогу принципа Беллмана, определяется выражением /г2 (и2) = с*/<?2(иг)\/Е(С2 I где Д2(и2)— плотность распределения Р2(е?г12), а с* — нормирующая константа.

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

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

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

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

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

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

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

В п. 3.3 получены оценки оптимального параметра расщепления К* и трудоёмкости в упрощённой модели для вариантов аналогового моделирования, весового моделирования, а также для нового варианта весового моделирования с ветвлением, который, как показали расчёты, может дать существенный выигрыш в трудоёмкости.

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

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

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

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

q ае~ах w(fj.) , .

перехода домножается на величину--—--—-, где р (х)— модели-

Яо р(х) r(ß)

руемая плотность распределения х, г (у) — моделируемая индикатриса, а д0 — моделируемая вероятность выживания.

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

„ /*оо _2„ —2сгх Г1 .,,.2/,Л

= (Ч

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

При этом для оптимизации р(х) и г(у) целесообразно рассматривать соответствующие задачи вида

f°° 91{х) 1

р = arg mm max / , ч ах, (2)

реР ßl<ß<02 Jo р(х)

где Р—некоторое семейство непрерывных положительных плотностей, ß—параметр, др £ Р, ßi < ß < /?2-

roo q2(x)

Обозначим J{p\ a) = / dx. Решение задачи (2) при опре-

J о Р(х)

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

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

pEP ¿=1,2

где Qi = gßi, i = 1,2.

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

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

Теорема 4.1. При д\ ф д2 решение минимаксной задачи (3) определяется выражением р5(х) = С{\*)^/\*g\{x) + (1 - \*)д%(х), где величина А* является единственным решением относительно А <Е (0,1) уравнения J(p\;gi) = J{p\\92), 0 < А < 1.

Рассмотрим теперь семейство Р\ удобно моделируемых методом "суперпозиции" плотностей вида р\(х) = \<j\(х) + (1 — Х)д2(х), 0 < А < 1, в предположении, что gi ф д2.

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

Л */ 9l(x)+92(x)

определяется выражением р1{х) = ---.

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

В п. 4.3 исследуется минимаксная оптимизация стандартного МПТ. Коэффициент и рассматривается в интервале 1 < а < t < оо. Для семейства экспонент Р2 = {se-sx, 1 < s < t = сг2/сгi}, соответствующего простейшему варианту стандартного МПТ, плотность р2{х) = s*e~s*x является решением задачи (2) при s* = ^ G (l,min{2,í}).

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

G{p*q\ 10) = 1.544, G(pí; 10) = 1.557, G{p*2; 10) = 3.025.

Здесь G{p*\t) — значение величины минимакса на оптимальном решении задачи (3) для S1)2 = {1, £}. Хорошее совпадение величин G(pq; 10) и G(pl; 10) объясняется следующим утверждением.

р*(х)

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

pj(l) t-+oo

Vx > 0.

Таким образом, для оптимизации параметрических оценок в задачах теории переноса целесообразно использовать вспомогательные плотности вида р\.

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

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

Лемма 4.6. Пусть плотность р* (х) минимаксна для £^2 = {l,i}? причем J(p*;e-X) = J(p*-,te~tx), и пусть

. d2J{p*-ae-"x) f°° (4а2х2 - 8ах + 2)е~2°х ,

mm -—-- — / --——--dx > 0.

l<a<t da"1 J 0 p {x)

Тогда p* минимаксна для £ = [l,i].

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

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

среды при этом имеем: р{Кр[и)) — — —J(p\cre~ax). Здесь для

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

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

Для численного эксперимента используется стандартная индикатриса Хеньи-Гринстейна, причем s = E/i — средний косинус угла рассеяния.

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

*/ \ т(у) + и>2{у) , ч _ ( \ • л о зования плотности г1{у) = ---, ги{(у) = -Шэ^у), г = 1,2,

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

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

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

й!(у)=1-^ + 8б(у-1), -1 < у < 1. (4)

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

о- —> ОТ = (1 — Я -> Я = ^, ю(у) -» (5)

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

Для сравнения были вычислены оценки функционалов Ра и Рс в транспортном приближении соответственно плотности столкновений 1р(г), являющейся решением "плоского" интегрального уравнения Пайерлса для параметров (5):

ф) = щ [" ф')К{д\г-г'\)йг' + (6)

V О

с ядром К(х) = ^ Г° ей.

2 с

Уравнение (6) решалось путём разбиения отрезка [О, Л"] на М равных частей со следующей дискретизацией: (¿зг = гг = (г —

1/2)/г, /г = Н/М, г = 1,...,М,

M Zj + h/2

y>i = crqi£>j I K{ö\Zl-z'\)dz'+ г = 1,..., M. (7)

J = 1 Zj-h/2

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

Отметим, что число обусловленности матрицы системы (7) оценивается величиной +, не зависящей от М и II.

1 - q

На этой основе для решения системы с дальнейшей высокоточной оценкой требуемых функционалов был использован известный метод LDLT факторизации.

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

Заключение

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

с вариацией параметра индикатрисы показана целесообразность использования плотности, оптимальной на семействе плотностей типа "смеси", легко моделируемых методом суперпозиции.

Основные публикации по теме диссертации

[1] Михайлов Г.А., Роженко С.А. Модификация двухэтапных алгоритмов метода Монте-Карло на основе свойств симметрии первого этапа // ЖВМиМФ, —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.

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

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

[4] 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. Peterburg VVM со. Ltd., 2009. —Vol. I. —P. 181-186.

[5] 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.

[6] 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.

Подписано в печать 28.06.2013 Формат 60x84 1\16 Усл. печ. л. 1 Объем 16 стр. Тираж 100 экз. Заказ № 183 Отпечатано Омега Принт 630090, г. Новосибирск, пр. Ак.Лаврентьева,6 email: omegap@yandex.fu

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

РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ НАУКИ

ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ГЕОФИЗИКИ

04201361705 На правах рукописи

Роженко Сергей Александрович

Использование свойств симметрии и подобия в алгоритмах метода Монте-Карло

01.01.07— Вычислительная математика

ДИССЕРТАЦИЯ

на соискание ученой степени кандидата физико-математических наук

Научный руководитель: чл. корр. РАН Г.А. Михайлов

НОВОСИБИРСК - 2013

Оглавление

Введение ........................................................................................3

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Литература......................................................................................83

Введение

Свойства симметрии и подобия изучаемых процессов позволяют строить эффективные вычислительные алгоритмы с использованием соответствующих аналитических преобразований осреднения, поворота, сдвига и т.п. В методе Моите-Карло это связано с использованием условных математических ожиданий или в рандомизированном варианте—алгоритмов расщепления. Разработка таких приемов восходит к известной работе [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) случайная величина.

где Р\(-) — вероятностная мера в 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 — "выборку по важности". Известно, что оптимальная плотность вероятности для такой модификации, согласно стохастическому аналогу принципа Беллмана [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 рассмотрена модифицированная модельная задача, в которой объемлющий шар разбит па две области — внутренний шар, и внешний слой. Во внутреннем шаре происходит деление частиц, т.е. в точке столкновения с некоторой вероятностью происходит ветвление траектории.

Используемую в данной главе для численных экспериментов среду внешнего сферического слоя можно представить, как среду с такими же сечени