Содержание к диссертации
Введение
1. Оптимизация весовых методов Монте-Карло по вспомогательным переменным 13
1.1. Модификация фазового пространства и весовой оценки 13
1.2. Оптимизация моделирования по части переменных 16
2. Ценностное моделирование 24
2.1. Эффективность ценностного моделирования начального распределения , 24
2.2. Моделирование начального распределения процесса для некоторых задач теории переноса 30
2.3. Частичное ценностное моделирование элементов траектории 3G
2.4. Метод исследования дисперсии несовой оценки 42
3. Исследование дисперсии весовых оценок в методе "экспоненциального преобразования" 48
3.1. Асимптотическая оптимизация распределения длины свободного пробега 48
3.2. Исследование дисперсии весовых оценок 50
3.2.1. Одномерный вариант 50
3.2.2. Сферический вариант 57
3.2.3. Численные расчеты для сферического варианта 62
Заключение 69
Литература
- Оптимизация моделирования по части переменных
- Моделирование начального распределения процесса для некоторых задач теории переноса
- Метод исследования дисперсии несовой оценки
- Исследование дисперсии весовых оценок
Введение к работе
Введение
Разработка алгоритмов численного статистического моделирования в настоящее время имеет особое значение в связи с возможностью их идеального распараллеливания путем распределения численных статистических испытаний по отдельным процессорам. Практически статистическое моделирование наиболее часто используется для решения задач физики и техники, в основе которых лежат вероятностные модели, связанные с некоторыми цепями Маркова [2J,[15],[16]. В принципе, такие задачи можно решать непосредственно численно моделируя траектории этих цепей. На основе исходного "феномеїіологического" описания проблемы можно даже строить весовые модификации алгоритма, домножая вспомогательный "вес" после каждого элементарного перехода на отношение соответствующих (может быть обобщенных) плотностей исходного и моделируемого распределений. С другой стороны, моделирование траекторий частиц можно интерпретировать как алгоритм оценки функционалов от соответсвующего интегрального уравнения второго рода, ядро которого совпадает с плотностью перехода базовых цепей Маркова [2],[14]. Однако зачастую прямое моделирование не позволяет оценивать изучаемые величины с требуемой точностью. В этом случае можно использовать весовые методы Монте-Карло, состоящие в том, что на ЭВМ моделируется подходящая цепь Маркова, а требуемые функционалы оцениваются с помощью веса, который после очередного перехода в цепи домножается на отношение ядра интегрального уравнения к переходной плотности. Несомненный плюс рассмотрения соответствующих интегральных уравнений - возможность детального изучения эффективности различных весоиых модификаций. В частности, на основе таких уравнений можно разрабатывать "ценностные" алгоритмы с малыми вероятностными погрешностями, что особенно важно для оценок малых вероятностей [2]. В последнее время было выяснено [1], что для построения весовых модификаций может быть особенно эффективным
Введение
увеличение размерности фазового пространства путем включения в число фазовых координат моделируемых вспомогательных случайных величин. Соответствующая факторизация ядра базового интегрального уравнения и вспомогательной плотности перехода иногда дает требуемую весовую модификацию.
Изложим более подробно введенные понятия. Математическая модель ряда прикладных задач строится на основе рассмотрения некоторого скачкообразного обрывающегося с вероятностью единица однородного марковского процесса (см. например, [2]), При этом траектория процесса вполне определяется ее состояниями в моменты скачков, т.е. фактически можно рассматривать обрывающуюся однородную цепь Маркова с заданной переходной функцией P(x',S), где х Є X — тп-мерное евклидово пространство, S С X - измеримое по Лебегу множество. Для построения весовых алгоритмов моделирования целесообразно рассматривать соответствующую условной мере P(x',S) обобщенную субстохастичсскую плотность перехода к(х',х). Обобщенная плотность распределения к(х\х) определяется для любого х' Є X равенством
f h{x)P{x',dx)=- f k{x\x)h(x)dx V/ієОД,
X X
где Cq{X) - множество непрерывных ограниченных функций. Отметим, что здесь и далее рассматриваются и ненормированные (т.е. не обязательно вероятностные) распределения.
Настоящая работа ориентирована на приложения в теории переноса частиц, в которой кроме измеримых плотностей, соответствующих абсолютно непрерывным распределениям, используются также "дельта-функции", означающие интегрирование по некоторым многообразиям меньшей сравнительно с т размерности (см., например, [2] раздел 2.1.1). Использовать обобщенные плотности (вместо интегрирования по соответствующим мерам) в теории статистического моделирования предложил Н.Н. Ченцов [3] в связи с тем, что такой подход упрощает
Введение.
построение и реализацию модификаций моделирования. Это важно для настоящей работы, так как в пей рассматриваются дополнительные фазовые переменные, причем базовые переменные - координаты и скорости - связаны с дополнительными функционально, т.е. их условные распределения определяются соответствующими "дельта-функциями". Предполагается, что
J к{х\х) dx = q(x') <1-S, 5>0.
Величина q(x') имеет смысл вероятности необрыва траектории в заданной точке х'. Вследствие последнего неравенства цепь обрывается с вероятностью единица и среднее число состояний конечно. Пусть
X0,Xi,.. .,xN
- однородная обрывающаяся цепь Маркова, определяемая плотностью f(x) распределения начального состояния xq и субстохастической обобщенной плотностью перехода к(х\х). Здесь N - номер состояния, в котором реализуется обрыв траектории {иначе, момент остановки). Ясно, что обобщенная плотность распределения состояний, непосредственно следующих за начальным, выражается равенством
>і(я) = J f(x')k(x',x)dx' = [Kf](x).
Следовательно, обобщенная плотность распределения среднего числа фазовых состояний цепи
71=0
где <рп{х) - плотность распределения состояний номера п, представляет собой ряд Неймана для следующего интегрального уравнения 2-го рода:
<р(х)= [ k(x',x)
Введение
или tp = Kip 4- f. Это уравнение будем рассматривать в пространстве Ni(X) обобщенных плотностей мер ограниченной вариации, так как ряд
^2 ( ^""Л ^ ) в силу условия q(x') < 1—5 сходится для любой h Є Cq(X).
71=0 \ /
Для заданной функции h Є Cq(X) рассмотрим также сопряженное уравнение
где [К*(р*](х) = j k(x,x')(p*(x') dxr. Предполагается, что х
1С є [Со(Х) -> С0{Х)]. Справедливо представление
П=0 71=0
где 5(-) - плотность локализованного в соответствующей точке источника.
Методы Монте-Карло обычно используются для оценки линейных функционалов вида
/
оо
фЩх) dx = 5D(/f"/, /і), he C0{X).
x "=o
Если реализуется прямое моделирование исходной цепи Маркова, то для оценки величины Ik используется соотношение
" N
Ih = E ^h(xn)
,п=0
Однако можно моделировать и другую, вспомогательную, цепь Маркова с плотностью перехода р(х', х), взаимно регулярной с к(х', х), и начальной плотностью 7г(ж), взаимно регулярной с f(x). Кроме того, предполагается, что р(х',х) ^ 0 и тг(ж) ^ 0 на носителях функций к(х',х) и /(ж) соответственно. При выполнении этих "общих условий
Введение
несмещенности" отношения к(х',х)/р(х',х) и f(x)/iT{x) имеют смысл и определяются отношениями измеримых сомножителей рассматриваемых функций. Это позволяет ввести вспомогательные веса по формулам:
^Ы~ф^У Qn~Q"-1P(^-i,x„)- (а3)
Полагая
= 53 Qnh(xn),
имеем /д — Е (см., например, [2]). Величина представляет собой стандартную весовую оценку "по столкновениям". Справедливы следующие соотношения:
- ^U & = Ms) + Х>ЛЫ- (0.4)
Заметим, что в теории переноса, если в точке столкновения величина h(xn) рассчитывается после выбора нового направления пробега частицы, то оценку (0.4) принято называть оценкой "по рассеяниям" [4].
Известно [2], что для оценки функционала можно построить большое количество различных оценок, среднее значение которых равно /д (например оценка по поглощениям). Общим недостатком таких оценок является то, что они включают в себя, как правило величины, которые сильно меняются в реальных задачах, либо трудно вычисляются. Практические расчеты показывают, что для задач теории переноса значительно более эффективно разрабатывать модификации основной оценки по столкновениям, чем применять другие виды оценок (за исключением задач о вычислении вероятности поглощения какого-либо типа) [2].
Функцию <>*() в теории весовых методов Монте-Карло принято называть "функцией ценности" в связи с ее вероятностным представлением:
Введение
<р*{х) = Щх = h(x) + E^Qn/i(i„), xQ = хл (0.5)
причем для оценки (0.4) величина Е^ определяется рядом Неймана для следующего уравнения [4]
g = k(2y>* - h) + Kfa (0.6)
где К* - оператор с ядром к2(х,х')/р(х,х').
Введем спектральный радиус р{К) по формуле
р(К) = p(IC) = inf || К И1/",
Известно [2], что если р(Кр) < 1 и /2/тт Є N\(X), то D < +оо. Также известно [17],[2], что, если h{x) > 0 и
fc(o/,z)^) №)y^) fQ ?)
то D = 0. Отмстим, что реализация оценки с нулевой дисперсией невозможна по следующим причинам: 1) изначально функция (р* неизвестна; 2) вследствие того, что вероятность обрыва тождественно равна нулю, необходимо моделировать бесконечные цепочки. Поэтому на практике используют приближенное моделирование по ценности, основанное на замене в характеристиках цепи (0.7) функции <р* па функцию р ~ <р* (замечательным свойством получаемого алгоритма является его независимость от постоянного множителя функции р). Кроме того, начиная с т - го состояния вводят поглощение с вероятностью, обеспечивающей обрыв траектории с вероятностью 1. Например, если вместо <р*{х) здесь использовать функцию:
т{х) = Ор*(х)[1 -Ь е{х)}, \е{х)\ < б,
то при малом є величина D мала [4]. Соответствующие алгоритмы объединяются под названием "моделирование по ценности" [2],[4].
Введение
Как правило, при увеличении размерности фазового пространства, переход х —» х' осуществляется в результате выбора совокупности значений вспомогательных случайных величин, причем, если соответствующие функции ценности определяются формулами типа (0.5), то ценностное моделирование всех элементарных вспомогательных переходов дает оценку с нулевой дисперсией [1]. Однако на практике осуществляется весовая модификация моделирования только части вспомогательных переменных, например в теории переноса наиболее часто используется модификация моделирования длины пробега на основе так называемого "экспонециального преобразования" [2]. При этом оказалось, что для некоторых задач частичная ценностная модификация вспомогательных переменных может увеличивать дисперсию по сравнению с прямым моделированием [12], а в некоторых случаях возникает вопрос о конечности дисперсии и, как следствие, величины трудоемкости s = f()D, где () - среднее время расчетов на ЭВМ для получения одного выборочного значения
Далее следует краткое содержание диссертационной работы по главам и параграфам.
Оптимизация моделирования по части переменных
В дальнейшем вспомогательные переменные будут рассматриваться всегда, поэтому будут использоваться более простые обозначения, принятые в теории весовых оценок (см.[2]): , ср, р , Jft, К, Кр.
Здесь будут использованы две, вообще говоря, векторные вспомогательные случайные величины ti,t2, т.е. ядро имеет вид: k((t, х), (t , х )) - S{x - х {х, t )) i(:c, t\)k2{{x, t\), t 2).
Рассмотрим сопряженное уравнение (0.2) в множестве Сь{Х) -неотрицательных функций из CQ(X), дополнительно предполагая, что К Є [Сь{Х) — Сь{Х)]. Поскольку для оценки по столкновениям справедливо: то целесообразно рассмотреть задачу равномерной относительно х G X минимизации величины Е . Согласно теореме 1.1, использование ценностных плотностей вида (1,3) для всех элементарных переходов даст абсолютный минимум : Е — ( (х))2 \/х Є X. Практически такая глобальная оптимизация моделирования весьма затруднительна, поэтому важно рассмотреть возможность уменьшения величины Е 2. путем оптимального подбора плотности распределения части вспомогательных случайных величин, например, ti.
Рассмотрим следующую задачу оптимизации : для фиксированной плотности Р2{{х, t\), t2) определить плотность pi(x, t[) таким образом, чтобы функция Е V.x Є X была минимальной при выполнении общих условий несмещенности. Такую плотность pi(x, t[) будем называть равномерно наилучшей.
Впрочем, и без рассмотрения выше предельного перехода ясно, что соотношение (1.8) является необходимым условием оптимальности алгоритма с Щ% = д(х), так как это соотношение означает его неулучшаемость.
Поэтому K2g д1//2д 5 [. Нетрудно аналогично получить неравенство K glg q(n l 2q \\ g \\ , т.к. при переходе к следующему значению п в соответствующем выражении дополнительно появляется агрегат, оцененный выше величиной q1 2. Таким образом, функция g равна риду Неймана для уравнения (1.6) при выполнении (1.9). Лемма 1.1 доказана.
Теорема 1.2. Пусть h Є Сь(х) и выполнены соотношегшя (1.5). Тогда решение g уравнения (1.8) в Сь{Х) существует и единственно, причем соответствующая плотность перехода (1.9) дает решение поставленной задачи оптимизации.
Заметим, что ранее в работе [4] была решена задача минимизации величины (1.6) путем оптимального подбора плотности распределения t2. Меняя способ факторизации субстохастического ядра (т.е. порядок выбора величин ti, 2)) для любой вспомогательной переменой можно решать соответствующую задачу оптимизации. Отметим, что представленная здесь оптимизация распределения ti ( первой по порядку выбора) вспомогательной переменной существенно упрощает вид уравнения (1.6) для построения соответствующей оптимальной плотности по сравнению с оптимизацией распределения t2
Покажем, как можно использовать полученные результаты. Рассмотрим систему линейных алгебраических уравнений :
Здесь слово частичная" для (1.12) означает, что оптимизируется распределение только одной вспомогательной случайной величины, а именно - номера следующего состояния. Вторая вспомогательная случайная величина: поглощение или выживание, моделируется физически.
Отметим, что разделе 3.1, на основе теоремы 1.2, осуществлена асимптотическая оптимизация распределения длины свободного пробега для практически важной задачи теории переноса .
В этой главе изучается эффективность ценностного и частичного ценностного моделирования, связанного с построением моделируемого распределения какой-либо из вспомогательных случайных величин путем умножения исходной плотности на функцию ценности, обычно соответствующей решению сопряженного уравнения. Получены условия, при выполнении которых ценностное моделирование начального распределения уменьшает дисперсию по сравнению с прямым. Используя полученные условия, проведен анализ эффективности ценностного моделирования начального распределения для некоторых задач теории переноса. Также, па примере модельной задачи теории переноса с "дельта" рассеянием и задачи оценки решения системы линейных алгебраических уравнений, показано, что частичное ценностное моделирование может увеличивать дисперсию оценки по сравнению с прямым моделированием. Доказано, что дисперсия весовой оценки в случае частично-ценностного статистического моделирования конечна. На этой основе предложен метод решения вопроса о конечности дисперсии весовой оценки с использованием мажорантного сопряженного уравнения. Результаты данной главы опубликованы в статьях [8],[12],[13],[18].
Моделирование начального распределения процесса для некоторых задач теории переноса
Как видно из формул дисперсий оценок для прямого, ценностного и оптимального моделирования начальной точки ((2.15), (2.16) и (2.17) соответственно), при моделировании осредняются величины, зависящие от {aij} И {/І}.
Перебирая всевозможные наборы {а } согласно (2.13), и варьируя распределение начальной точки формуле /2 = 1-/1, /i-m 0.1, (m -1,2,.., N), (2.18) мы получили, что при h\ = 1,/Ї2 = 2, в случае прямого моделировании вдоль цепи, ценностное моделирование начальной точки не увеличивает дисперсию по сравнения с прямым для 2618 наборов {ац} из возможных 2916 для любой начальной плотности (2.18), т.е. в среднем для 90 % наборов {dij}. Для частичного ценностного и частичного оптимального моделирования вдоль цепи эта доля составила 92 %. Заметим, для перечисленных выше наборов, было получено, что max -— = 0.93, — 0.23, т.е. ценностное моделирование начальной точки может весьма существенно уменьшать дисперсию сравнительно с прямым, не значительно проигрывая, оптимальному. В тех случаях, когда для какого либо набора {а } и начального распределения (2.18), D D, оказалось, что Г1 0.03, т.е. ценностное моделирование начальной точки увеличивает дисперсию сравнительно с прямым, но весьма незначительно. Замечено, что при уменьшении шага сетки получаются аналогичные результаты.
Из полученных результатов можно сделать вывод, что в основном, в случаях прямого, частичного ценностного и частичного оптимального моделирования вдоль цепи оценки по столкновениям для задачи (1.10), ценностное моделирование начальной точки не увеличивает дисперсию по сравнению с прямым. Причем, чем эффективнее моделирование вдоль цепи, тем чаще ценностное моделирование начальной точки уменьшает дисперсию по сравнению с прямым.
Теоретически важным в этом пункте является установление факта существования задач, для которых частичное ценностное моделирование вдоль цепи увеличивает дисперсию по сравнению с прямым моделированием.
Как было показано в пункте 2.3, при решении некоторых задач использование частичного ценностного моделирования вдоль цепи может быть неэффективным по сравнению с прямым моделированием. Более тогО) оказалось, что в задаче об оценке вероятности вылета частиц из полупространства —со z Н (см. главу 3), при ценностном моделировании части вспомогательных переменных для оценки по столкновениям не выполняется стандартный критерий р(К ) 1 конечности дисперсии (см. параграф 3.1). Далее будут получены утверждения, позволяющие определить ограниченность дисперсии весовой оценки безотносительно значения спектрального радиуса. Пусть іІ = (іі,«/2)єГ = ГіхГ2 - набор из двух вспомогательных величин (возможно векторных), выбор которых осуществляется для реализации перехода в цепи Маркова. В 2.4- Метод исследования дисперсии весовой оценки модифицированном фазовом пространстве TxX = {(t ,x)}, субстохастическос ядро имеет вид (см. главу 1) k((t, х), (tf, О) = 5(х - х (х, t ))fti(s, t[)k2{{x, Єг), t 2), где х (х, t ) - функция, определяющая новое значение стандартных евклидовых координат через х и значения вспомогательных переменных t . Дополнительно предположим, что Ух Є X і ki{x,t\)dl\ = \-а{х) .
Пусть случайное значение величины 1 2 моделируется согласно заданному распределению, а случайное значение величины t\ - с использованием соответствующей вспомогательной функции ценности, т.е. переходные плотности имеют вид
Дисперсия оценки по столкновениям х при частичном ценностном моделировании первой вспомогательной случайной величины t[ (соответственно (2.19),(2.20)j конечна.
Метод исследования дисперсии весовой оценки Доказательство. Поскольку h, ц Є Сь(Х), то справедливо следующее утверждение [4]: n=0 (2.22) Покажем, что ряд (2.22) сходится для любой h Є Сь{Х). Учитывая (2.19)-(2.21) и условие на ki(x,t[), прямой подстановкой можно проверить равенство: [К 1(а0 = г. РіОМ ї)) J J к2{{х, t[),t 2)5(x -x {x,t ))tp (x )dx dt 2 4 , X dt\ = Гі 2((2:, ), ) ( (1, ) &i = = [K y \{x) J hix dt , = [ V](a;)(l-a(a;)) - ( )-/ ))(1-0 )), Гі которое можно переписать в виде Подставив в это равенство под знак оператора К вместо р равную ей функцию K ip + a(tp — h) 4- h, получаем равенство Р = Kfy + К; (а{ р -h) + h\ +a( / - h) + h. Далее сделаем такую же подстановку в последнем равенстве иод знак оператора К 2, и так далее.
Метод исследования дисперсии несовой оценки
Простые вычисления (см.[4], раздел 1.6) показывают, что и для оценки по столкновениям величина р(Кр) определяется соотношением (3.6). Совпадение значений спектрального радиуса соответствует тому очевидному факту, что дисперсии оценок по рассеяниям и по столкновениям конечны или бесконечны одновременно, т.е. для одних и тех же весовых алгоритмов. Поэтому величина р{Кр) инвариантна относительно любого сдвига фазовой точки х вдоль цепочки элементарных вспомогательных переходов. Таким образом, уравнение для критического значения с имеет вид q2 r lnS±i-l 2(1-р)с c/a-1 и при 1 — р — q совпадает с соответствующим характеристическим уравнением Милна (3.4). Поэтому, при использовании экспоненциального преобразования со значением параметра с = 1/L, в случае изотропного рассеяния имеем р(Кр) — 1, то есть при этом не гарантирована даже конечность дисперсии.
Однако предложенный в конце главы 2 метод исследования конечности дисперсии позволяет установить практически важный факт ограниченности дисперсии в методе экспоненциального преобразования при 0 с 1/L. Известно [6], что для реальной индикатрисы w( ,//), уравнение (3.4) устанавливает соответствие с -» {q(c), ас(ц)}, причем q{c) Ч, если с (0,1/L]. Рассмотрим функцию ценности вида: Vc(z / ) =ас(д)ехр{-(Я-г)с}, (3.7) для которой справедливо равенство: причем
Заметим, что, если вероятность выживания в точке столкновения равна #(с), то экспоненциальное преобразование эквивалентно частичному ценностному моделированию с использованием функции ценности вида (3.7). Для исходной вероятности выживания q экспоненциальное преобразование при с (О,1/L] эквивалентно использованию частичной ценностной плотности, помноженной на q{c)/q. Заметим, что предложенный метод исследования позволяет использовать такую плотность, поскольку в этом случае в выражении для K tp c появляется множитель q/q{c) 1 (см. док-во теоремы 2.5). В случае, когда q/q(c) — 1, т.е. при с = 1/L, нетрудно проверить справедливость условия теоремы 2.5 на функции h, а.
Численно решалась задача оценки вероятности вылета частицы за границу Я = 0 из точки (z, ) — (—20,1) с использованием (3.3) при q — 0,7 для изотропного рассеяния, В этом случае экспоненциальное преобразование при с = 0 есть прямое моделирование длины пробега, при с — \/l q 0,548 - асимптотически (для II — оо) оптимальное (см. раздел 3.1), а при с = 1/L « 0.829 -приближенно ценностное моделирование длины пробега, для которого функция ценности определяется выражением (3.5) и, соответственно,
В таблицах 3.1,3.2 представлены результаты расчетов (здесь и далее полученные путем моделирования 107 траекторий частиц) с использованием оценки по рассеяниям (0.4) и модифицированной "бернуллисвой" оценки, которая "подсчитывает" веса Q N вылетающих частиц, причем Q N = QNe-{H-z»)r:
Отметим, что дисперсия модифицирован)юй "берпуллиевой" оценки при ценностном моделировании длины пробега также будет конечна, т.к. эта оценка получается путем стандартной рандомизации рассмотренной оценки по рассеяниями (см. [2], раздел 4.4.1).
Исследование дисперсии весовых оценок
Рассмотрим некоторые задачи, связанные с прохождением частиц через оптически толстый слой.
Пусть внутри исходной сферы находится абсолютно "черный" шар радиуса Ri R, попав на поверхность которого, частица поглощается с вероятность 1. Полагаем, что внутри слоя #i г R рассеяние частицы изотропно и вероятность выживания после столкновения q — 0.7. Пусть на поверхности черного шара помещен точечный источник, испускающий частицы по закону Ламберта [19), для которого где wo - вектор внешней нормали к сфере радиуса R\. Необходимо оцепить вероятность Р = P(R) вылета частицы за пределы сферы радиуса R.
При прохождении через толстый однородный слой вещества Ri г R поток частиц асимптотически по г убывает, как е а\ поэтому, учитывая раздел 3.2.2, в данной задаче длину свободного пробега целесообразно моделировать с использованием модифицированного сечения (ЗЛО) при с 0.
Численно решалась задача оценки функционала P(R) при R = 10, Я]. = 2. В таблице 3.6 представлены результаты расчетов с использованием модифицированной бернуллиевой оценки, которая подсчитывает веса вылетающих частиц:
Приняты обозначения: P{R) - статистическая оценка величины P(R), а - соответствующая оценка среднеквадратичной вероятностной погрешности. Выбраны следующие, значения экспоненциального преобразования: с — 0 - соответствующее прямому моделированию, с = 1 — q = 0.3 - наибольшее значение, при котором имеет место частичное ценностное моделирование со вспомогательной функцией ценности, удовлетворяющей мажорантному сопряженному уравнению, с = у/1 — q л; 0.548 - критическое или максимальное значение с, для которого доказана конечность дисперсии весовой оценки. Был также просчитан вариант для с — 1/L « 0.829, что соответствует критическому значению параметра экспоненциального преобразования для аналогичной задачи на плоскости. Заметим, что здесь и далее для каждого расчета моделировалось 107 траекторий частиц, причем для коррелирования результатов использовался одинаковый набор псевдослучайных чисел таким образом, что моделирование траекторий одного номера в разных вариантах задачи начиналось с одинаковых псевдослучайных чисел из мультипликативного конгруэнтного генератора с параметрами М = 517,т — 2ю [2].
Как видно из таблицы 3.6, частичное ценностное моделирование (с = 0.3) заметно уменьшает дисперсию по сравнению с прямым моделированием, при этом не значительно проигрывая экснонециалыюму преобразованию с максимально допустимым значением с = 0.548. В случае, когда конечность дисперсии не гарантирована (с — 0.829), модифицированное моделирование, тем не менее, незначительно проигрывает прямому.
Пусть теперь точечный ламбертовский источник частиц расположен на внутренней поверхности сферы радиуса R и требуется оценить вероятность Р = -P(-fti) поглощения частицы на поверхности абсолютно черного шара радиуса Ri.
Здесь необходимо рассчитывать прохождение частиц через оптически толстый слой R г Ri на внутреннюю сферу радиуса R\, и, согласно сказанному в разделе 3.2.2, в данной задаче длину свободного пробега целесообразно моделировать с использованием модифицированного сечения (3.10) при с 0.
Численно решалась задача оценки функционала P{R\) при R = 10, R\ — 2. В таблице 3.7 представлены результаты расчетов с использованием модифицированной бернуллиевой оценки, которая подсчитывает веса проходящих частиц: Модифицированная берпуллиева оценка вероятности P(Ri) при использовании экспоненциального преобразования
Как видно из таблицы 3.7, также, как и в предыдущей задаче, частичное ценностное моделирование на основе мажорантного сопряженного уравнения (с = 0.3) заметно уменьшает дисперсию по сравнению с прямым моделированием, при этом не значительно проигрывая экспонециальному преобразованию с максимально допустимым значением с = 0.548. В случае, когда конечность дисперсии теоретически не гарантирована (с = 0.829), модифицированное моделирование также незначительно проигрывает прямому.
Кроме того, были проведены расчеты с использованием оценки по столкновениям. Чтобы оценить вероятность Р(Яі), достаточно в каждой точке столкновения (r,w) вычислять величину [2]: конус направлений с вершиной в точке г, "опирающийся" па внутреннюю сферу радиуса Я . Заметим, что вычисление интеграла h(v,u ) в каждой точке столкновения представляет достаточно трудоемкую процедуру. Для данной задачи целесообразно использовать рандомизированную оценку по столкновениям [2], а именно, после моделирования нового направления и рассчитывать величину:
В таблице (3.8) приведены расчеты функционала P{R\) с использованием рандомизированной оценки по столкновениям при экспоненциальном преобразовании без вылета. В данной модификации длина пробега I внутри слоя R\ г R моделируется согласно плотности C(l + jufJ)) -"" , поэтому траектория может оборваться только к результате поглощения с вероятностью 1 — q. После моделирования рассеяния, в точке (r7t,u/) вес определяется по формуле: где величина 1 п есть I или l R, в зависимости от направления пробега и. Заметим, что, в данной модификации длины пробега, время моделирования одной траектории определяется только заданной вероятностью поглощения 1 — q, поэтому среднее время моделирования одной траектории одинаково для любого значения параметра с.
