Об использовании "частичного ценностного" моделирования в полупространстве фазовых координат Текст научной статьи по специальности «Математика»

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

Том 13, Специальный выпуск 4, 2008

Об использовании "частичного ценностного" моделирования в полупространстве фазовых

координат*

И.Н. Медведев Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: medvedev79@ngs.ru

Effectiveness of modifications of the partial "value" modelling is considered in this paper. These types of modelling are related to a construction of the simulated distribution for some auxiliary random variable by multiplying the initial density by the "value" function, which usually corresponds to a solution of an adjoint integral equation of the second kind. Using our criterion based on the majorant adjoint equation, we focus on an implementation of the partial "value" modeling in the half-space of phase coordinates.

Введение

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

Рассмотрим интегральное уравнения второго рода

где X — т-мерпое евклидово пространство; <*, к Е СЬ(Х); СЬ(Х) — множество неотрицательных непрерывных функций. Здесь и далее предполагается, что К * Е [СЬ(Х)

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00046-а) и программы "Ведущие научные школы" (грант № НШ-4774.2006.1)

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

или ф* = K*ф* + h,

(0.1)

СЬ(Х)], а для субстохастической обобщенной плотности перехода к(х,х') имеет место неравенство

J к(х, х') ¿х' = д(х) < 1 — 8, 8 > 0.

Для построения оценки метода Монте-Карло величины <^*(х) введем цепь Маркова хо,х1,..., xN с плотностью перехода р(х, х'), взаимно регулярпой с к(х',х), причем р(х,х') = 0 па носителе функции к(х,х'). Последнее замечание позволяет нам ввести вспомогательные веса:

= 1, <3га = Яп-1-т-'—г

Р(хп—1, хп)

и оценку по столкновениям'

N

Н(х) + £ ЯпЦхп), (0.2)

п=1

где N — случайный номер последнего состояния цепи.

Функцию <£*(•) в теории весовых методов Монте-Карло принято называть "функцией ценности" в связи с ее вероятностным представлением:

N

^*(х) = Е£х = к(х) + Е^2,ЯпЬ(хп), хо = х, (0.3)

п=1

причем для оценки (0,2) величина Е£% определяется рядом Неймана для уравнения [2]

д = к(2ц)* — к) + К*рд, (0.4)

к2(х,х')

где К * — оператор с ядром —--—,

Известно „|, ™ если с,^ра1ый радиУс ,(*,) < 1. го ^ < +оо. Также ИЗВес,

но [1, 2] что если Н(х) > 0 и

, к(х,х')р*(х')

р(х,х) = —г—-п. , , (0,5)

1 7 [К *р*](х) К

ТО Б^х = 0.

Как правило, переход х ^ х' осуществляется в результате выбора совокупности значений вспомогательных случайных величин, например, номера типа столкновений, углов рассеяния и длины свободного пробега частицы при моделировании процесса переноса, Несомненным плюсом данного подхода является то, что при соответствующей весовой модификации моделирования только части вспомогательных переменных можно определить критерий ограниченности дисперсии оценки по столкновениям без исследования значения спектрального радиуса [3], Пусть 1 = (¿[,¿2) € Т = Т1 х Т2 — набор из двух вспомогательных величин (возможно, векторных), выбор которых осуществляется для реализации перехода в цепи Маркова, В модифицированном [4] фазовом пространстве Т х X = {(1', х)} субстохастическое ядро имеет вид

к((1, х), (1', х')) = 8(х' — х'(х, Ь'))к1(х, ¿[)к2((х, V1), ¿2),

или

где

где х'(х, 1) — функция, определяющая новое значение стандартных евклидовых координат через х и значения вспомогательных переменных 1 . Рассмотрим цепь Маркова с субстохастической плотностью перехода

р((1, х), (1', х')) = ¿(х' - х'(х, 1'))Р1(х, г'М(х, г'1),г'2).

Пусть случайное значение величины t/2 моделируется согласно заданному распределению, а случайное значение величины ¿1 — с использованием соответствующей вспомогательной функции ценности, т, е, переходные плотности имеют вид

Р2 ( (х; ¿1) ) ¿2) = к2((х; ¿1) ) ¿2) ' ЫХ[К*и](х) и*(х)-к(х) [ }

М*Л) = , (0-7)

и(х)

и1(х,Ь'1) = J У 8(х' — х'(х, t'))k2((x,t'1),t'2)u(x')dx'dt'2 =

Т2 X

= ! к2((х, ¿1), ¿2)и(х'(х, t'))dt'2,

Т2

а функция и € СЬ(Х) удовлетворяет мажорантному уравнению

и = К *и + к, (0.8)

причем

- к -вирр к С вирр к, — < С2 < сю Ух € вирр к. (0,9)

к

В этом случае было доказано [3], что дисперсия оценки по столкновениям при "частичном ценностном" моделировании первой вспомогательной величины ^ (т. е. соответственно (0,6)) будет конечна. При моделировании согласно (0,7) дисперсия оценки будет также ограничена, если для элемента ядра к1(х,^1) имеют место неравенства

J k1(x,t'1)dt'1 = 1 — а(х) < 1 — е < 1 Т1

или

к(х)

—-— < С < сю Ух £ вирр к. а(х)

Отметим, что имеют место аналогичные утверждения относительно дисперсии весовой оценки и для частичного ценностного моделирования второй вспомогательной переменной t/2 [3],

Известно, что если для каждого элементарного перехода использовать плотность, равную произведению исходного ядра и соответствующей вспомогательной функции ценности [4], то = 0, Такая глобальная оптимизация моделирования весьма затруднительна, поэтому на практике осуществляется весовая модификация лишь части

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

1. Частичное ценностное моделирование в полупространстве фазовых координат

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

где функция и удовлетворяет (0,8) и (0,9), а Х0 — полупространство фазовых координат, где будет применяться частичное ценностное моделирование.

Теорема 1. Дисперсия оценки по столкновениям £х при частичном ценностном моделировании первой вспомогательной случайной величины ъ\ в полупространстве Х0 (т. е. соответственно (2,1) ) конечна.

из которого в силу неотрицательности всех функций вытекает сходимость ряда

для х Е Х0, для х Е X\Х0,

(1.1)

и = К* х и + к + а (и — к)1 (х Е Хо) Ух Е X,

к

Поскольку к(2<^* — к) < Н2^>* < 2кС, то будет сходиться ряд

п=о

значение которого в силу неотрицательности всех функций [2] совпадает с величиной Е^2, что и требовалось доказать. □

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

Рассмотрим задачу об оценке вероятности вылета, частиц из полупространства —то < г < И. Будем считать, что вне этого полупространства находится абсолютный поглотитель, причем средний свободный пробег а-1 = 1 во всем пространстве. Рассеяние частицы в точке столкновения описывается симметричной нормированной плотностью ), где ц — косинус угла между направлением пробега и осью г. Вероятность выживания при столкновении в точке г < И равна ^ < 1,

Заметим, что уравнение (0,1) для описанной модели с учетом вспомогательной переменной I можно записать в виде

А 1

<^*(г, ¡) = ^У е-1! ') ¿1 + к(г,ц), г < И, (2,1)

0 -1

где А = +то при л < 0 и А = (И — г)/л при л > 0,

Известно, что в схеме по рассеяниям вероятность вылета частицы, стартовавшей в точке х = (г, ¡), определяется вели чиной Е^х (см, [2], разд. 2,4) при

) Г ехр{ —(И — г)/л} при г < И и л > ,

( , | 0 иначе. ^ '

Рассмотрим дополнительно уравнение (2.1) со свободным элементом ка(г, ¡) = а(л)к(г, ¡), где а(л) совместно с параметром с = 1/Ь удовлетворяет характеристическому уравнению Милна [6]

1

(1 — сц)а(ц) = qJ )а(ц')d¡'.

-1

Для реальной индикатрисы рассеяния имеем а(л) > е > 0, и подстановкой нетрудно проверить (см. например [2], разд. 2.4), что

<р*а(г, ¡¡) = ЕСх(а) = а(л) ехр{ —(И — г)/Ь}. (2.3)

Известно, что моделирование длины пробега согласно плотности е-1ес^1 приводит к экепоненциональному преобразованию, для которого а' = 1 — сл [2]. При этом, если обрыв траектории моделируется физически и с = 1/Ь, то экспоненциальное преобразование эквивалентно частичному ценностному моделированию первой (см. (0.6) и (0.7)) вспомогательной переменной I — длине пробега с использованием (2.3). Отметим, что предложенный в работе [3] метод исследования конечности дисперсии позволил установить важный факт ограниченности дисперсии в методе экспоненциального преобразования при 0 < с < 1/Ь.

И=0

из точки ( г, 1) с использован нем ка при q = 0.7 для изотропного рассеяния. В раИ

чаетичное ценностное моделирование для к = ка существенно уменьшает дисперсию по сравнению с прямым моделированием. Однако при малых расстояниях (| г0| = 1, 2, 3) частичное ценностное моделирование увеличивает дисперсию. Для коррелирования результатов использовался одинаковый набор псевдослучайных чисел таким образом, что

106

II. H. Медведев

Ценностное моделирование для оценки

(*, 1) (-20,1) (-10,1) (-5,1) (-2,1) (-1,1)

1) 3.34 • 10"Y 1.48 • 10"3 9.28 • 10"2 1.114 2.553

<70 1.47 • 10~Y 5.57- 10"6 4.04 • 10"5 1.08 • 10"4 1.29 • 10"4

00.829 1.41 • 10"8 4.59 • 10"6 7.09 • 10"5 2.61 • 10"4 3.45 • 10"4

00.829,Х0 9.87- 10"9 5.72 • 10"6 4.05 • 10"5 1.08 • 10"4 1.29 • 10"4

моделирование T = 5 ■ 107 траекторий одного номера в разных вариантах задачи начиналось с одинаковых псевдослучайных чисел из мультипликативного конгруэнтного генератора с параметрами M = 517, m = 240 [1], Здесь и далее а — соответствующая оценка среднеквадратичной вероятностной погрешности,

В таблице представлены результаты расчетов оценки функционала (2,3) с использованием оценки по рассеяниям при прямом (<т0), частичном ценностном (â0.829) и частичном ценностном (â0.829,Xo ) только в полупространетве X0 фазовых переменных моделировании длины свободного пробега, В последнем случае длина свободного пробега моделировалась согласно следующей плотности:

PiXo(x,l)=[ еслихеХо ={(z,ß): z < H - 5, ß > О},

I e-1, если иначе.

Как видно из последней строки таблицы частичное ценностное моделирование только в полупространстве X0, так же как и ценностное (â0.82g) во всем пространстве, существенно уменьшает дисперсию по сравнению с прямым моделированием при больших расстояниях источника от границы H. Однако при небольших расстояниях от границы

X0

точным ценностным моделированием.

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

[1] Ермаков С.М., Михайлов Г.А. Статистическое моделирование. М.: Наука, 1982.

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

[3] Михайлов Г.А., Медведев И.H. Метод исследования дисперсии весовой оценки численного статистического моделирования // Докл. РАН. 2006. Т. 406, № 2.

[4] Михайлов Г.А. Построение весовых методов Монте-Карло на основе увеличения размерности фазового пространства // Докл. РАН. 2003. Т. 389, № 4. С. 461-464.

[5] Михайлов Г.А., Медведев И.Н. Эффективность "ценностного" моделирования в методе Монте-Карло // Докл. РАН. 2005. Т. 401, № 1. С. 16-20.

[6] Дэвисон Б. Теория переноса нейтронов. М.: Атомиздат, 1960.

Поступила в редакцию 20 февраля 2008 г.