Неимитационные оценки и улучшение способов математических ожиданий при статистическом моделировании переноса частиц Текст научной статьи по специальности «Математика»

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

Том 17, № 2, 2012

Неимитационные оценки и улучшение способов математических ожиданий при статистическом моделировании переноса частиц

А. И. Хислмутдинов1, Б. В. БАНЗАРОВ2 1 Институт нефтегазовой геологии и геофизики СО РАН, Новосибирск, Россия 2Новосибирский государственный университет, Россия e-mail: KhisamutdinovAI@ipgg.nsc.ru, banzarov@gmail.com

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

Ключевые слова: статистическое моделирование, методы Монте-Карло, имитационные и неимитационные оценки, способ математических ожиданий (expected value estimator), приближения к решению сопряженной задачи, уменьшение трудоемкости.

Введение

Методы Монте-Карло активно используют для моделирования процессов переноса частиц при решении различных научно-технических задач [1-4]. Вычисляются и анализируются как показания приборов, так и различные распределения частиц: пространственные, временные, энергетические и др. Один из ранних примеров изучения распределений — работа [5], в которой рассматривались пространственно-энергетические распределения от мононаправленного моноэнергетического точечного источника гамма-квантов. Математической моделью исследуемых здесь процессов переноса частиц является марковский скачкообразный процесс, и фазовая плотность в этом процессе рассматривается как удовлетворяющая уравнению переноса. Пара, составленная из процесса и оценки (функционала от траектории процесса), образует метод Монте-Карло. В настоящей работе в представленных методах соответствующий марковский процесс считается заданным и фиксированным и речь идет о выборе и построении оценок. Цель работы — исследование различных оценок для задач, в которых размер области, отвечающей носителю весовой функции в линейном функционале, достаточно мал. В связи с такими задачами принято говорить о проблеме вычисления малых вероятностей, имея

в виду то, что частицы либо редко пересекают данную область, либо редко в ней соударяются. Пусть р — малая вероятность и последовательность "вкладов" траекторий описывается схемой Бернулли. Тогда для относительной стандартной дисперсии на одну траекторию имеем выражение

что показывает неограниченность относительной дисперсии при p ^ 0. Еще в 19501960 гг. был предложен и усовершенствован способ преодоления трудностей "малых вероятностей", введена и усовершенствована неимитационная оценка — "expected value estimator", или оценка способа математических ожиданий (СМО) [1, 3, 6, 7]. В методе с этой оценкой относительную дисперсию удавалось ограничивать. Иногда считают, что в данную оценку введена итерация весовой функции. Интерес к развитию и улучшению СМО сохраняется (см., например, [8, 9]).

Хорошо известная "локальная оценка" для вычисления "потока в точке" родственна оценке способа математических ожиданий, и много усилий было затрачено на усовершенствование этой оценки и улучшение ее порядка сходимости. Неплохие результаты были, в частности, получены, используя неимитационные оценки "нулевого класса" и априорную информацию о решении сопряженной задачи [10-12]. Аналоги данного подхода, в котором (сопряженная) функция ценности (р приближалась весовой функцией, были применены и для улучшения СМО [10, 13]. Вместе с тем следует напомнить, что уже в первых работах об оценках СМО появилось словосочетание "длинные ряды экспонент", отражающее определенную тяжеловесность их конструкции. Часто вклады от определенных областей малы, однако их необходимо вычислять, чтобы оставаться в рамках несмещенности. В основе предлагаемого здесь авторами подхода — стремление к простоте и легкости имитационных оценок, что и реализуется в построенных неимитационных оценках. Конструкция последних позволяет им "работать" как имитационные в областях с "малым вкладом" и как оценки СМО в остальных областях. Дисперсии неимитационных оценок и оценок способа математических ожиданий близки, в то же время первые обладают меньшей трудоемкостью, поскольку их реализация требует меньшего числа вычислительных операций. В сравнении с [13] в настоящей работе исследуется общий нестационарный случай, рассматриваются вычисления как средних "по поверхностям", так и средних по пространственным областям, и здесь, что принципиально, вводятся другого типа эффективные дополнительные случайные выборки и некоторые новые приближения (р к решениям сопряженной задач. Одним из основных объектов внимания авторов является средний поток по объему конечного кругового цилиндра. Эта форма характерна для детекторов частиц; вычисление показаний детекторов, включая случай детекторов малого объема, представляет собой важную задачу. Актуален также вопрос об оптимизации методов вычислений [7, 14].

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

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

1. Терминология, обозначения и основные понятия

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

X — фазовое пространство координат г, направлений П и энергий Е; X = М3 0 0 [0, гс);

х = (г, П, Е), X е X, 51 = {П : П е М3, |П| = 1};

V = v(E) — абсолютная скорость частиц, £ — временная переменная, £ е [0, гс);

а — поглощающее состояние;

{Х4}£=0 — марковский скачкообразный случайный процесс с состояниями в X и {а}, описывающий распространение частиц; процесс считается непрерывным справа;

Р0 — цепь Маркова с состояниями в [0, гс) 0 (X и {а}), соответствующая процессу

ш — траектория цепи Р0, ш = {(£0, х0), (£1, х1),..., (tj, х),..., (£п, хп) = а}, где п — номер состояния обрыва цепи, xj• = х^., Xj-0 = х^.-0; состояние (£п,хп-0) непосредственно предшествует обрыву цепи;

?(£, х) — плотность начальных состояний цепи Р0, сосредоточенная на множестве [0, гс) 0 X;

Е(х) — полное макроскопическое сечение взаимодействий частиц (со средой), (х) — макроскопическое сечение рассеяний; считаем, что для всех х верно соотношение Ез(х) < Е(х) < Се < гс;

С(П, Е ^ П', Е'|г') — индикатриса рассеяния частиц;

т(£, х) — фазовая плотность частиц, Ф = vm, ^ = ЕФ — плотности потока и взаимодействий; все эти функции считаются обобщенными плотностями, удовлетворяющими соответствующим интегральным формам уравнения переноса;

1у, !я — линейные функционалы от Ф, подлежащие вычислению по формулам

Ту = ¿£ / ¿г / ¿П / ¿ЕФ(£,х)Еу(£,х), (1)

^ М J ¿Б J ¿П J ¿ЕФ(£,х)|(П,п5)|Ез(£,х), (2)

где (V), (Б) — заданные область и поверхность в М3, п^ = п^(г) — нормаль к поверхности (Б), Еу(£,х), Ез(£,х) — заданные ограниченные кусочно-непрерывные функции, равные нулю вне (V) и (Б) соответственно;

Т+ и С + — хорошо известные линейные интегральные операторы, ядра которых есть

| г — г/ |

- Г Я(т'+Ш,Е)Л

Т(г',т' ^ г,т\П,Е) = у(Е)У(т,Е)е 0 5(т - т' - у(Е)П(г - г')),

У (т Е')

СоШ',Е' ^ П,Е\т) = / С(П',Е' ^ П,Е\т)

У(т, Е')

соответственно;

К + — композиция операторов Т+ и С+: Iе + = Т+С+; Iе — оператор, сопряженный к 1С +;

р = Iе+р + f — интегральная форма уравнения переноса для р на X, f = Т+д; <^(ш) — оценка "по соударениям" для вычисления 1у,

п Еу (г х)

Ыш) = ^2к(гз,тз,^j-l,Ej-l), где ь(г,х)= у^ ; (з)

Пя(ш) — оценка "по пересечениям" для вычисления 1я,

п

Пя(ш) = ^ Ня(гз-1,тз-1 ^ гз ,тз\^з-1,Ез-1), 3 = 1

ь

Ня (гз-1,тз-1 ^ гз ,тз\Пз-1,Ез-1) = £я (ti,тi, ^з-1,Ез-1), (4)

г=1

где 1з — количество пересечений траектории ш с поверхностью (Б) в пределах полуинтервала [тз-1, тз), гг, тг — момент времени и точка ¿-го пересечения, оценки ^(ш) и пя(ш) принято называть имитационными для вычисления 1у и 1я соответственно;

С1,у, С1,я — оценки способа математических ожиданий с дополнительной случайной выборкой для вычисления 1у и 1я [7]:

(1,у (ш,шр) = Н(г1,х1-о) + кк(гр,хр\Ьз ,хз-0),

з=1

(1,я(ш, шр) = Ня(го, то ^ 11, п\По,Ео) + ^ кНя(гр, х%,хз-о),

з=1

здесь кк, кНя — оценки итераций Iек(г,х), IеНя(г,х), для которых справедливы выражения

Мр^х(кк(гр,хр\г,х)) = ! кк(гр,хр\г,х)ру (гр,хр\г,х)сИрс1хр = К к(г,х),

Мрх(кНя(гр,хр\г,х)) = У кНя(гр,хр\г,х)'р/я(гр}хр\1}х)(Ир(1Ер(тр(1Б = КНя(г,х), где (гр,хр) — переменные дополнительной случайной выборки, шр = ((гр, хр),..., (гп, хП)).

2. Неимитационные оценки

§ 2.1. В работе рассматривается семейство ^^ из "нулевого" класса оценок [10, 12, 13]. Для вычисления /у оценки этого семейства определяются конструкцией

С2,У (ш,шр) = ДМ^ + ^-о) - <^С3+1,Ж3+1 ^ (7)

3=1

где (/? — заданная измеримая интегрируемая функции, при этом <^(а|£,ж) = 0. Оценка считается заданной, если кроме функции </? также определена оценка итерации К</?, обозначенная как к</?(£р,Жр-0|£3-, ж3-0), где ¿р,Жр_0 — переменные дополнительной случайной выборки. Описанные выше оценки <^(ш) и (ш, шр) содержатся в ZPу и соответствуют определенным функциям (/?.

Аналогичным образом введем семейство Z(рS для вычисления функционала /88:

п

(2,8 (Ш^р) = ^[Я^ (^з-1,Гз-1 ^ ¿з ,Гз |Пз-1,Ез-1) + 3=1

+к^р,жр_0^,хз-0) - )], (8)

при этом оценки п8(ш) и £1)8(ш,шр) содержатся в Zр8.

Оба семейства, z2Py и , являются подмножествами более широких классов, в которых содержатся оценки с нулевой дисперсией [10, 12]. Сформулируем условия о конечности дисперсии для оценки (7).

Предложение. Пусть функции Л,(-); <^(ф) и к<^(ф) равномерно ограничены по всем своим переменным. Пусть также существует такая постоянная с < 1, что вероятность события "не менее последовательных соударений" меньше с п.н. Тогда дисперсия оценки (2,у(ш,шр) конечна.

Доказательство. Согласно первому условию предложения существует такая постоянная С, что для каждого звена произвольной траектории ш справедливо неравенство

Мз",х3-0) + ,хз-0) - ^СЗ+ъХ'+^з)| < Сс,

откуда следует, что |^2,у(ш,шр)| < пС, где п — случайная длина траектории. Хорошо известно, что случайная величина в правой части последнего неравенства имеет ко-

оо

нечный второй момент. Действительно, М(пС^)2 = ^^ Рп(пС)2, где Рп — вероятность

п=1

обрыва траектории на п-м звене. Согласно второму условию, для всех целых к > 0 и всех натуральных п Е [к^0 + 1, (к + 1)^0] имеем Рп < ск. Окончательно получим

^ к3о+3о -(2 ^^ ^^ „к„2 < с2^ ~3„к,

к=0 п=к30+1 к=0

МСС(ш) < £ скп2 < ^(к + 1)2.

Ряд в правой половине неравенства сходится, что доказывает предложение. □

Выбор эффективной оценки заключается, во-первых, в нахождении функции </?(£, ж|£', ж'), во-вторых, в определении плотности для дополнительной случайной выборки. Конструкции (7) и (8) позволяют строить оценки с конечной дисперсией и трудоемкостью меньшей, чем у стандартных оценок. Далее в разделе будет рассмотрен выбор функций </?(£, ж|£', ж'); дополнительные случайные выборки анализируются в § 2.2.

Назначение функции р(г,х\г',х') производим согласно следующим двум доводам. Как уже было отмечено, при р(г, х\г', х') = к(г,х) оценка принимает вид оценки математических ожиданий (•). В данном случае итерации Iер(г',х') имеют ненулевые значения для всех точек взаимодействия т', и их нужно как-то оценивать. При этом, если т' находится на большом расстоянии от объема (V) или поверхности (Б), то значения итераций являются малыми величинами. Последнее привело к тому, чтобы полагать функцию р^Ю равной нулю при условии, если точка т' не находится в заданной области О', которую мы считаем близкой к (V) или (Б). При таком назначении р(ф) нет необходимости рассчитывать итерации Iер(г',х') при каждом соударении.

Второй довод касается случаев, когда область (V) или поверхность (Б) имеют большую протяженность по пространственным переменным. Тогда расстояние \т - т'\ может очень сильно варьироваться и возрастает вероятность того, что вклад в оценку будет малым. В данной работе второй довод рассмотрен только для поверхностей (Б). Он заключается в том, что (Б) будет разделена на ближнюю и дальнюю части и р(г, х\г', х') будет полагаться равной нулю для всех т, принадлежащих дальней части. При таком условии дополнительная случайная выборка разыгрывается только по ближней части поверхности.

Формально вышеуказанные доводы принимают следующую реализацию. Пусть О — некоторая область в К3. Тогда для оценок (2-у (•) функции р(ф) выбираются в виде

р(1,х\г',х') = ( к(г,х) при т' 6 (9)

I и во всех остальных случаях.

Пусть дополнительно (Б') — часть поверхности (Б), не включающая в себя точки, отдаленные от т'. Тогда для оценок (2. я(•) функции р(ф) имеют вид

{ Ня- > пРи ' 6 6 (10)

и во всех остальных случаях.

Далее мы переходим к рассмотрению конкретных типов областей (Б) и (V). Детально изучим следующие важные для приложений типы. В качестве поверхностей (Б) рассматриваются цилиндрические и плоские кольца, задаваемые уравнениями

(Бс) = {х, у, г\х2 + у2 = р0, го - Н < г < го + Н}

и

(Бр) = {х,у,г\(ро - Н)2 <х2 + у2 < (ро + Н)2 г = го}. При этом центральной окружностью будем называть окружность во, заданную в виде

во = {х, у, г\х2 + у2 = ро, г = го}.

В качестве области (V) рассматривается конечный цилиндр:

(V) = {х, у, г\х2 + у2 < ро, го - Н < г < го + Н}.

§ 2.2. В соответствии с вышеизложенным определим явный вид функций р(г, х\г', х') в оценках для вычисления 1яс, 1яР и 1у. Для 1у полагаем в (9)

О = {х, у, г\х2 + у2 < (ро + Ь)2, го - Н - Ь < г < го + к + Ь}, 0 < Ь < ж,

ния точки соударения на примере плоского кольца Бр

где параметр Ь должен быть задан. Для 1яс и 1яР полагаем

Б' = {х, у, г\(х2 + у2 + г2 + р\ - Ь2)2 - 4р20(х2 + у2) > 0},

или, другими словами, область Б' задается в виде тора, окружность вращения которого совпадает с центральной окружностью колец, а радиус образующей окружности берется равным параметру Ь, который также должен быть задан. Отметим, что если параметр Ь больше радиуса центральной окружности ро, то это не приводит к недоразумениям при реализации алгоритма выборки.

Поверхность Б' для колец Бс и Бр определяем следующим образом. Пусть расстояние от точки соударения т' до оси OZ больше радиуса центральной окружности ро. Тогда для центральной окружности и точки соударения строится двугранный угол, ребро которого проходит через точку соударения параллельно оси центральной окружности. Грани угла касаются центральной окружности, указывая на ней пару точек, определяющих еще один двугранный угол, ребром которого является ось центральной окружности, а грани проходят через построенные точки. Этот угол делит кольца Бс и Бр на две части — дальнюю и ближнюю по отношению к точке взаимодействия. Именно ближнюю часть будем рассматривать в качестве поверхности Б'. Пусть теперь расстояние от точки соударения до оси OZ не превышает ро, и нельзя построить касательный двугранный угол. В этом случае проводим через точку соударения плоскость, перпендикулярную радиусу-вектору с координатами (х',у', 0), где х' и у' — декартовы координаты точки соударения. Пересечение данной плоскости с центральной окружностью также показывает пару точек, в свою очередь определяющих двугранный угол, делящий кольца на две части, одна из которых есть (Б'). Итак, функции р(1, х\Ь', х') определены для всех случаев. На рис. 1 представлены схемы, поясняющие алгоритмы построения поверхности (Б') для двух положений точки т' на примере плоского кольца БР.

3. Дополнительные случайные выборки

§ 3.1. Построим дополнительные случайные выборки для оценивания итераций от заданных функций р(1,х\1' ,х'). Как правило, индикатриса рассеяния С является смесью парциальных индикатрис, описывающих разные взаимодействия, и итерация Кр представляет собой смесь вкладов от каждого из них. Для простоты изложения будем счи-

тать, что индикатриса описывает только один тип рассеяния. При этом ограничимся важнейшим случаем, когда энергия частицы после взаимодействия определяется только углом ее рассеяния относительно направления до взаимодействия. Тогда = (П, П'), а $(^'|Е,г') — заданная плотность по которую считаем ограниченной. Также будем считать, что функция Л,(£, ж) внутри объема (V) зависит только от переменной Е. С учетом перечисленных предположений выражения для итераций К(/? в оценках и £2,у имеют вид

КК^з,Гз, П-ьЕ^Н / х

х ^р(^'|Ез--1, Гз)е-т^Е^(¿5, , П, Е^) (11)

2п

и

К^з,Гз, Пз-1,Ез-1)^ ¿П з_) —у (М'|Ез-_1 ,Гз) х

./ _(гз ,Ез_1)

те

х У ^/Е(гз + /П, Еу)е_т(г*^^^ ^ + Гз + /П, П, Еу^ =

^П-

(Гз ,Ез_1) 1

■0Ои'|Ез-_ь гз)й(Еу) (е_т(г*) - е_т(г*(12)

Е(Гз ,Ез_1) 2П

соответственно. В выражениях (11), (12) ¿5 и Г5 — момент времени и точка пересечения луча Гз+/П, / > 0, с поверхностью (Б); Е^, Еу — значения энергии, которые выражаются через где = (Пз-1, П); Гз,1, Гз,2 — точки пересечения луча Г3 + /П, / > 0 с границами объема (V).

Пусть в (11), (12) ^П = ^а^, где а и ^ — переменные сферической системы координат с началом в точке Г3, ^ — косинус широтного угла, отсчитываемого от положительного направления оси О^, а — азимутальный угол. Тогда равномерная выборка по этим переменным является наиболее естественной дополнительной случайной выборкой, но с той особенностью, что при ее реализации достаточно просто производится только выборка по переменной а. Определение границ интегрирования по переменной ^ трудоемко, и в настоящей работе предлагаются такие алгоритмы, которые либо с помощью замен переменных, либо путем выбора плотности распределения выборки позволяют обойти трудоемкие операции. В первую очередь необходимо разбить область Б' на две подобласти Д и Д. В качестве Д выбираются сами области (V) и (Б) вместе с некоторыми их окрестностями, и при Г3 Е Д предлагается использовать выборку, равномерную и по а, и по Оставшаяся часть Д' рассматривается в качестве Д, и при Г3 Е Д в качестве выборки используются далее описываемые алгоритмы. Подобласти Д и Д можно выбирать таким образом, что почти все соударения будут происходить в Д, и в этом случае трудоемкость всего алгоритма будет слабо зависеть от того, какая дополнительная случайная выборка используется при Г3 Е Д.

Для кольца (Бс) рассмотрим следующий вариант. Пусть переменная а уже разыграна. Тогда на кольце однозначно определен отрезок, расположенный параллельно оси О^, и для всех его точек координаты ж^ и постоянны. Таким образом, ^ зависит только от и функция ^(г^) является монотонной и дифференцируемой. Теперь можно провести замену переменных с ^ на ^ и использовать дополнительную случайную

выборку по . Замена переменных имеет вид ^а^ = ^а ^

. Использовать для

выборки переменную zs удобно из тех соображений, что границы ее изменения всегда заданы и постоянны. В работе предлагается равномерная выборка по этой переменной.

Совместная плотность по переменным а и zs имеет вид p(a,zs) = -, где а0 еще

4a0H

требуется определить (см. ниже § 3.2).

Рассмотрим случай вычисления функционала IsP. Интеграл (11) представляется

в виде поверхностного интеграла, т. е. dQ = j ( S,—dS. Поверхностный интеграл

\rs - Tj\2

в свою очередь может быть представлен в виде двойного интеграла по переменным Ps и фs цилиндрической системы координат с осью OZ и началом при z = 0, т. е.

, )\ dS = ii^—^psdpsdфs,

\rs - Tj\2 \rs - r при этом фs зависит только от а. Окончательно имеем

dфs

dQ = dS = P

\rs - rj\2

das

\(ns, -\ dpsda.

\rs - rj\2

Итак, мы перешли от интегрирования по а и ^ к интегрированию по а и ря, что проще, поскольку границы изменения для ря фиксированы. В настоящей работе предлагается равномерная выборка по а и ря, и совместная плотность распределения имеет вид, идентичный выборке для 1яс: р(а,ря) = --—. Способ вычисления ао, а также построение

4аоН

моделирующих формул для алгоритма описаны в § 3.2.

При вычислении "объемного" функционала 1у выражение (12) имеет представление в виде интеграла по dQ и dl,в котором интегрирование по с11 проводится аналитически, и остается выражение в виде интеграла по СП. Равномерная выборка по П имеет конечную дисперсию, но при этом остаются сложности при определении границ интегрирования по переменной ^ (П = (а,ц)). Для розыгрыша ^ предлагается использование простой "индуцированной" выборки. Рассмотрим цилиндрическую систему координат с началом в точке Tj и осью Z', параллельной оси данного цилиндра; переменными этой системы координат являются а, р' и г'. Если переменная а выбрана, то тем самым определяется прямоугольник (Р\а), являющийся пересечением ребра а с (V). Пусть Q(а, р', г') есть случайная точка, имеющая равномерное распределение в (Р\а). Рассмотрим теперь переменные (а,ц,1) в сферической системе координат с центром в точке Tj. Распределение точки Q описанным способом индуцирует некоторое распределение по переменным ^ и I. Полученное распределение мы предлагаем использовать для дополнительной случайной выборки. Выражение для плотности этого распределения и моделирующие формулы для переменных выборки будут представлены в § 3.2. Отметим, что индуцированная выборка другого типа ранее рассматривалась в работе [14], где случайная точка выбиралась равномерно в области (V). Сравнение этих двух вариантов приводится ниже.

§ 3.2. Построим моделирующие формулы для дополнительной случайной выборки, используемой при вычислении 1я и 1у. В § 3.1 были определены переменные и намечены плотности распределения для дополнительных случайных выборок. Сконструируем сначала выборку по углу а, алгоритм которой является одинаковым для обоих типов поверхностей (Б) и цилиндра (V). Во всех случаях необходимо построить двугранный угол, касательный к цилиндрической поверхности, разыграть в полученном растворе угол и построить секущую плоскость, соответствующую этому углу.

и

Обозначим координаты точки соударения Tj как (Xj ,yj ,Zj) и введем величину р1 = 'sjx] + У2 • Тогда требуемый двугранный угол имеет раствор величиной 2а0, где

Í. Ро

arcsm — при р0 < р1,

РП

2 пРи Ро > pi.

Зададим далее переменную а* таким образом, что cos а* = — и sin а* = —. Пусть Yi —

Pi PI

стандартное случайное число. Тогда выборочное значение а определяется выражением а = а * + (1 — 2^1)а0. Итак, выборка по а сформулирована.

Координаты (x*,y*) пересечения "разыгранной" секущей плоскости с цилиндром имеют вид

X * = xj + pS cos а, y * = yj + sgn(yj )ps sin а, (13)

где

pS = — (xj cos а + yj sin а) — , если p0 < p1,

pS = —(xj cos а + yj sin а) + , если p0 > p1.

Таким образом, если вычисляются функционалы Is, то определена точка на центральной окружности, имеющая координаты (x*,y*,z0), зная которые, можно разыгрывать

x

переменные zS и pS. Для SP прежде всего определим угол фS так, что cos фS = — и

Ро

y

sin ф3 = — .В этом случае Ро

xs = (ро + (1 — ъ)Н) cos фs, ys = (ро + (1 — ъ)Н) sin фs, йф^ Р2 + р\ — 2pipo cos фs

ZS = z0^ — = -т-2-.

аа р1р0 cos фц; — р0

Для Se имеем

* * , (л . m Ф (xs— xj)2 + (ys— yj)2 (лл.

xs = x , ys = y , zs = Zo + (1 — 2ji)H, — =-¡---. (14)

dzs \ts — Tj

Итак, дополнительные случайные выборки для обоих типов поверхностей построены и по построению их дисперсии конечны.

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

Найдем явный вид плотности распределения для выборки по переменной Переменная а определяет прямоугольное сечение цилиндра (P|а), высота которого равна 2H,

ширина — 2^J~{xjcoSа+yjsiñ'а)2—(Xf+У2—рP0l), а произведение этих величин дает его

площадь |£|. Далее, если мы разыграем точку равномерно внутри данного прямоугольника, то к плотности по ^ приводят следующие преобразования:

6»2 L

dS f dpdz [ . n 1П [ /d/ — ' sin ftdtf '

M2

J |S| J |S|

(S) (S)

sin 0|S |

L 2 L 2

2 sin 0|S |

PM =

L2 - _ (L - Li)(L2 + Li)

:i5)

2 sin 0|S| 2 sin 0|S|

где ^ = cos

Теперь запишем моделирующие формулы для ^ и точек rj-,1 и rj,2 из (12). Для этого вычислим координаты разыгранной точки в прямоугольнике (Р|а). Пусть Y1 и y2 стандартные случайные числа. Подобно (13) имеем

xy = Xj + cosaR(yi), yy = yj + sgn(yj) sinaR(71), zy = zo + (1 - 2y2)H,

:i6)

где

(71) = —(xj cos a + yj sin a) + (1 — 27^^(xj cos a + yj sin a)2 — (ж2 + y2 — p0).

Отсюда

Zy — Zj

^(xy - xj)2 + (yy - yj)2 + (zy - zj)2

Чтобы найти координаты точек rj,1 и rj,2, необходимо в (16) в качестве аргумента функции R подставить соответственно 0 и 1. Итак, выписаны все требуемые величины и алгоритм определен.

Поскольку в выражение для оценки итерации входит множитель, обратный плотности (15), то выражение (L2 — L1)(L2 + L1) привносит в нее две особенности при L2 ^ L1 и при L2 ^ 0, которые тем не менее не приводят к неограниченности оценки. В случае первой из них оценка остается ограниченной, так как разность экспонент в (12) при L2 ^ L1 стремится к выражению E(rj, Ej-1)(L2 — L1)eT(rj) и (L2 — L1) сокращается. Второй предельный переход не рассматривается, поскольку по построению оценки L1 и L2 ограничены снизу некоторым положительным числом.

В отличие от настоящей работы числитель плотности, построенной в [14], содержит выражение (L2 — L1)(L2 + L1L2 + L1) и разность (L2 — L1) является устранимой особенностью, а сумма (L2 + L1L2 + L1) при L2 ^ 0 представляет собой малую величину большего порядка по сравнению с (L2 + L1).

§ 3.3. Рассмотрим еще один подход к построению оценок с конечной дисперсией. В данной работе для того, чтобы оценка была ограниченной, при условии r' Е D'£ мы переходим к равномерным выборкам по угловым переменным. Это вызвано тем фактом,

что в случае использования выборок, построенных в § 3.1 и 3.2, при |rj — rs| ^ 0 оценка

const

итерации становится неограниченной и имеет асимптотический вид

lrj - rs|

. Функции

в этом случае имеют вид (9) и (10). Определим функции следующим образом. Пусть

0 во всех остальных случаях,

ям^х) = < н((',г'-(,г|п'в,)з--^ пригЕ^ гЕ(5,)

0 во всех остальных случаях.

При таком выборе ф(ф) оценка итерации всегда является ограниченной и дисперсия алгоритма конечна. Данный подход по ограничению оценок при дополнительных случайных выборках рассматривался в [4, 10], но по отношению к другим (родственным) проблемам. В этом случае не нужно определять, в какой подобласти области Б' произошло соударение, так как всегда реализуется единый алгоритм дополнительной случайной выборки. Привлекательным моментом в данном подходе является отказ от равномерных выборок по угловым переменным при г' Е .

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

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

§ 4.1. Изучим задачу о пространственно-временном распределении быстрых нейтронов в конфигурации скважина — пласт. Скважина радиусом 9.85 см заполнена водой; состав пласта — песчаник; импульсный моноэнергетический источник нейтронов равномерно распределен на окружности радиусом 7 см, ось которой совпадает с осью скважины и является осью О^ декартовой системы координат, начало которой находится

в центре окружности. Частицы испускаются в момент £ = 0 с энергией 14.1 МэВ и вып

летают под углом — относительно оси О^. Под быстрыми подразумеваются нейтроны

с энергиями выше 1.4 МэВ. В этой области энергий нейтроны испытывают неупругие рассеяния, являющиеся источником гамма-излучения неупругого рассеяния, которое регистрируется детекторами приборов. Для интерпретации данных измерений наряду с показаниями детекторов, размещенных в скважине, необходимы понимание характера и особенностей распространения нейтронного газа и гамма-квантов в системе скважина — пласт, а также знание, какие области горной породы и как влияют на показания детекторов. Кроме того важно понимать, каково пространственно-временное распределение быстрых нейтронов и соответствующих источников гамма-квантов. Результаты изучения этой задачи изложены в работе [16], здесь ограничимся некоторыми из них.

Для функционала (2) рассмотрим весовую функцию (£, х), принимающую значение -г—--, если )| > еп, и 0, если )| < еп, где еп — некоторое

малое положительное число, |£| — площадь заданного плоского кольца (£) вида

5 = {(г,- - к)2 <г2 < (г,- + к)2, г = Z},

где г^ — радиусы центральных окружностей колец. Тогда функционал можно представить в виде функции от параметров весовой функции Е^(£,х). На рис. 2 приведены графики зависимости функционала от г^; к = 2 мм, г3 = 8, 9,... , 40 см, Z = 20 и 40 см.

Рис. 2. Пространственно-временные распределения быстрых нейтронов от импульсного источника, распределенного на окружности: а — 2 = 20, б — 2 = 40 см; при вычислениях применялась оценка С,2,Б

В качестве (Е1, Е2) и (¿1, рассмотрены энергетический интервал (1.4, 14.1) МэВ и временные интервалы (5, 10), (10, 20), (20, 30) и (30, 40) нс. Зависимости на графиках представляют эволюцию первоначального пучка частиц.

§ 4.2. Оценка использовалась также для решения задачи о пространственном распределении нейтронов вблизи границы скважины и пласта. Ранее эта задача рассматривалась в [15] для надтепловых нейтронов. Мы представим результаты для тепловых и частично надтепловых нейтронов.

Точечный изотропный стационарный источник расположен на оси скважины, радиус скважины 9.85 см, распределение начальной энергии определяется спектром 252СГ, скважина заполнена водой, состав пласта — песчаник с пористостью 0 и 30 %. Изучалось радиальное распределение нейтронов на разных уровнях относительно источника, представляемое в виде зависимости от параметра г^. Результаты вычислений для надтепловых нейтронов показали согласие с данными, опубликованными в [15]. Здесь будут представлены результаты, полученные для тепловых нейтронов, т. е. для энергетического интервала (0, 0.215) эВ. Кольца располагались на высоте 20 см от плоскости источника, радиусы центральных окружностей составляли 4, 6, ... , 16 см, пористость среды равнялась нулю. График на рис. 3, а демонстрирует полученное распределение. Отметим, что убывание этой функции в скважине и в пласте происходит по-разному.

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

Рис. 3. Результаты вычислений пространственного распределения тепловых нейтронов: а — потоки нейтронов, б — отношения трудоемкостей методов с имитационной (^1) и неимитационной (^2) оценкой для колец разной ширины

метим, что при вычислении пространственного распределения с оценкой (2 по одним траекториям вклады в функционалы от соседних точек являются коррелированными, что позволяет достигать гладкости кривых при меньшем числе траекторий.

Поставим вопрос — как выигрыш в трудоемкости зависит от количества функционалов по одним и тем же траекториям? При малом числе поверхностей время вычисления вкладов в оценку мало по сравнению со временем моделирования траектории независимо от типа оценки. Однако при большом числе поверхностей это время увеличивается для неимитационной оценки, но остается малым для имитационной. Таким образом, программа с неимитационным методом работает медленнее программы с имитационным методом. При очень большом числе поверхностей имитационная оценка может выигрывать по сравнению с неимитационной, однако в задачах переноса число локальных экстремумов мало и для корректного расчета особенностей пространственного распределения не требуется большого числа поверхностей. Отметим, что наборы колец, рассмотренные в § 4.1 и в задаче о надтепловых нейтронах, описанной в § 4.2, являются избыточными по сравнению с достаточным набором, т. е. в данном случае можно было бы использовать меньшее число колец.

§ 4.3. Теперь сравним методы для вычисления функционалов вида /у. Рассмотрим вычисление средних потоков гамма-квантов по заданному объему (V). На оси скважины расположен цилиндр диаметром 1 см, центр которого находится на расстоянии 40 см от источника. Вычисляется средний поток гамма-квантов в цилиндре. Источник гамма-квантов — неупругие рассеяния нейтронов на атомах кислорода среды. Сравнение неимитационного метода с методом математических ожиданий показало: при одинаковых погрешностях в первом случае вычисление производится в полтора раза быстрее, чем во втором, что дает выигрыш по трудоемкости. Также было определено, что скорости вычислений имитационным и неимитационным методами практически равны. На рис. 4, а представлены абсолютные погрешности (корни из дисперсий расчета) имитационной оценки по длине пробега и неимитационной оценки. На рис. 4, б приведена кривая отношений трудоемкостей методов, т. е., фактически — отношений их дисперсий. Как и ожидалось, выигрыш в трудоемкости растет с уменьшением к.

Рис. 4. Абсолютные погрешности методов с имитационной оценкой по длине пробега (кривая 1) и неимитационной оценкой (кривая 2) при вычислении среднего числа соударений в цилиндре с полувысотой Н (а); отношение трудоемкостей методов с имитационной (^1) и неимитационной оценкой в зависимости от Н (б)

Аналогично кривым на рис. 3, б кривая отношений трудоемкостей не является монотонной, что объясняется случайным характером погрешности.

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

Список литературы

[1] Фано У., Спенсер Л., Бергер М. Перенос гамма-излучения. М.: Госатомиздат, 1963. 284 с.

[2] Метод Монте-Карло в проблеме переноса излучений / Под ред. Г.И. Марчука. М.: Атом-издат, 1967. 256 с.

[3] Спанье Дж. Метод Монте-Карло и задачи переноса нейтронов. М.: Атомиздат, 1972. 271 с.

[4] ХисАмутдинов А.И., Стариков В.Н., Морозов А.А. Алгоритмы Монте-Карло в ядерной геофизике. Новосибирск: ВЦ СО АН СССР, 1986. 157 с.

[5] Berger M.J., Spencer L.V. Some radiological applications of gamma-ray transport theory // Radiation Res. 1959. Vol. 10, No. 5. P. 552-570.

[6] Kahn H. Use of different Monte Carlo sampling techniques // Symp. on Monte Carlo Methods. N.Y.: Wiley, 1956. P. 146-190.

[7] ХисАмутдинов А.И. Об эффективности метода математических ожиданий для задач одного класса // Журн. вычисл. математики и матем. физики. 1967. Т. 7, № 4. С. 946-953.

[8] Banerjee K., Martin W.R. Using kernel density estimation for Monte-Carlo tallies with unbounded variance // 21st Intern. Conf. on Transport Theory. Torino, 2009. P. 8.

[9] Mosher S., Mauceo M., Spanier J. et al. Expected-value techniques for Monte Carlo modeling of well logging problems // Nuclear Instruments and Methods in Phys. Res. 2010. Vol. 613. P. 334-341.

[10] Хислмутдинов А.И. Несмещенные оценки в методах Монте-Карло для интегральных уравнений. Новосибирск: ВЦ СО АН СССР, 1986. 147 с.

[11] Воронин А.Ф., Хислмутдинов А.И. Метод Монте-Карло с дополнительной случайной выборкой для вычисления потока частиц в точке // Журн. вычисл. математики и матем. физики. 1985. Т. 25, № 8. С. 1155-1163.

[12] Khisamutdinov A.I. Unbiased non-simulation estimators in Monte Carlo methods and their applications in particle transport // Transport Theory and Statistical Phys. 1998. Vol. 27, No. 3-4. P. 303-316.

[13] Хислмутдинов А.И., Шишмареба С.О. Неимитационные оценки для методов Монте-Карло // Журн. вычисл. математики и матем. физики. 1992. Т. 32, № 1. С. 115-122.

[14] Хислмутдинов А.И., Шишмареба С.О. Метод с индуцированной плотностью дополнительной случайной выборки в способе математических ожиданий // VII Всесоюзное совещание "Методы Монте-Карло в вычислительной математике и математической физике" (тезисы докл.). Новосибирск, 1985.

[15] Khisamutdinov A.I., Voronin A.F. Spatial distributions of the epithermal neutron flux from a 252Cf source in a "Borehole-formation" system // Nuclear Geophysics. 1990. Vol. 4, No. 4. P. 437-441.

[16] Хислмутдинов А.И., Банзаров Б.В., Федорин М.А. Математическое моделирование нестационарного переноса частиц в задачах импульсного нейтронного-гамма каротажа. Новосибирск, 2008. 54 с. (Препр. ИНГГ СО РАН).

Поступила в 'редакцию 27 июня 2011 г., с доработки — 21 декабря 2011 г.