1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 13
Текст из файла (страница 13)
. . , ξ˜(m) , . . . −(6.1)это последовательность случайных состояний (случайных величин, случайных векторов), принимающих значения в одной и той же областиX ⊆ Rd , при этом число этих состояний бесконечно (см., например,[14]). С точки зрения компьютерных вычислений это является определенным недостатком, так как, с одной стороны, требуется исследоватьсходимость связанных с цепью (6.1) функционалов (усредненных характеристик) при m → ∞, а с другой стороны, численно (на компьютере)удается смоделировать конечное число выборочных значений состоянийцепи Маркова (6.1).В связи с этим появляется содержательная проблема правильной организации обрыва траекторий цепи Маркова, реализуемых на компьютере. Возможен «искусственный» обрыв цепи в фиксированный моментm = s; при этом приближаются только допредельные (неточные) значения требуемых функционалов.Более продуктивной является идея использования прикладнойцепи Маркова или цепи Маркова, обрывающейся с вероятностью80единицаξ (0) , ξ (1) , ..., ξ (N )(6.2)(см., например, [9, 13], а также подраздел 2.8 данного пособия, в частности, соотношение (2.28)) с начальной плотностью π(x) (это плотностьслучайной величины ξ (0) ) и переходной функциейp(x0 , x) = 1 − p(a) (x0 ) × r(x0 , x)(6.3)(см.
также соотношениездесь p(a) (x0 ) – вероятность обрыва, а (j−1) (2.29));00r(x , x) = rξ(j) xξ= x – условная вероятностная плотность распределения случайной величины ξ (j) при фиксированном ξ (j−1) (одинаковая для всех j = 1, ..., N ).Для моделирования траекторий прикладной цепи Маркова (6.2) следует использовать алгоритм 2.8.В дальнейшем будем полагать, чтоp(s) (x0 ) = 1 − p(a) (x0 ) ≤ q < 1 или p(a) (x0 ) ≥ 1 − q > 0;при этомZp(x0 , x) dx = p(s) (x0 ) ≤ q < 1;(6.4)Xэто объясняет использование здесь термина переходная функция дляp(x0 , x) вместо термина переходная плотность для «классической» цепи Маркова (6.1).Для дальнейших рассуждений важным будет рассмотрение последовательности функцийZψ (m) (x) =p(x0 , x)ψ (m−1) (x0 ) dx0 = K̂ψ (m−1) (x); ψ (0) (x) = π(x);X(6.5)здесь m = 1, 2, ..., а K̂ обозначает интегральный оператор с ядромp(x0 , x) (см., например, [19]).Важным является то обстоятельство, что при выполнении условия(6.4) оператор K̂ из (6.5) является сжимающим операторомв проRстранстве функций y из L1 (X) с нормой kykL1 (X) = X |y(x)| dx:kK̂y1 − K̂y2 kL1 (X) = kK̂ (y1 − y2 ) kL1 (X) ≤ qky1 − y2 kL1 (X) ,так как его норма в этом пространстве меньше единицы (см., например,[19]).81В свою очередь, рассмотрение функций (6.5) из пространства L1 (X)вполне естественно для многих практически важных приложений, таккак функции π(x), p(x0 , x) часто содержат интегрируемые особенности(вплоть до дельта-функций).Учитывая соотношения (2.4)–(2.7), можно трактовать функциюp(x0 , x)×ψ (m−1) (x0 ) как аналог (с точностью до соотношений (6.3), (6.4))совместной плотности распределения вектора ξ (m−1) , ξ (m) , а функциюψ (m) (x) из (6.5) – как аналог «безусловной» плотности распределенияслучайной величины ξ (m) .В дальнейшем будем называть функции (6.5) плотностями распределения состояний прикладной цепи Маркова (т.
е. слово «аналог» мы опускаем), а их бесконечную суммуψ(x) = ψ (0) (x) + ψ (1) (x) + ... + ψ (m) (x) + ...(6.6)назовем суммарной плотностью состояний прикладной цепиМаркова.В разделе 7 будет показано, что для функции (6.6) удается выписатьмарковское интегральное уравнение Фредгольма второго рода (см., например, [19]), что, в свою очередь, позволяет строить весовые оцениватели (монте-карловские оценки), нужные для приближенного вычисления линейных функционалов от решения этого уравнения.Идеологию, терминологию и прикладные аспекты представленной вразделах 6–8 теории весовых оценивателей линейных функционалов отрешений интегральных уравнений Фредгольма второго рода во многомопределили реальные расчеты для практически важных «физических»конструкций, среди которых особое место занимает классическая численная модель переноса малых частиц (см., например, [9, 13, 20]).Эта модель достаточно подробно представлена в данном разделе.Определенным дополнительным «оправданием» рассмотрения этоймодели является то обстоятельство, что соответствующие вычислительные алгоритмы были в свое время использованы для содержательныхрасчетов, связанных с разработкой атомных устройств (бомб, реакторови т.
п.), сыгравших важную роль в развитии экономики и укрепленииобороноспособности нашей страны.6.2. Простейшая численная модель переноса частиц. Пустьимеется выпуклая область G (мы рассмотрим случай G ⊆ R3 ), вообщеговоря, неоднородного вещества, в подобласти которой G0 ⊆ G (илина границе ∂G) расположен источник излучения – направленный илистохастический (рис. 6.1).82МОДЕЛЬ 6.1.
Малые частицы излучения («фотоны» – в дальнейшем будем использовать этот термин без кавычек), двигаясь попрямой, сталкиваются с хаотично расположенными более крупнымичастицами вещества в слое.В точке столкновения они либо поглощаются этими частицами,либо рассеиваются согласно некоторому стохастическому (вероятностному) закону распределения.После рассеяния малая частица вновь движется по прямой вплотьдо следующего столкновения.Процесс продолжается вплоть до поглощения малой частицы; приэтом внешность области G считается абсолютно поглощающим веществом (то есть вылет фотона за пределы области означает обрывтраектории).Нас будут интересовать различные усредненные характеристики описанного процесса.
В частности, весьма актуальной нередко являетсяпроблема вычисления вероятности P того, что частица излучения поглотится в области G (другими словами, это доля излучаемого потокамалых частиц, остающаяся внутри области G). Именно эту характеристику мы возьмем в качестве основного примера в данном разделе.Рис. 6.1. Схема простейшего процесса переноса частиц (модель 6.1)Важно отметить, что специально описанные (см. далее) точечные состояния столкновений фотонов с частицами вещества образуют прикладную цепь Маркова.83Соответствующий алгоритм прямого моделирования подразумевает реализацию n траекторий фотонов от «рождения» до «гибели»(поглощения) согласно алгоритму 2.8 и подсчет нужных характеристикпроцесса согласно стандартной схеме метода Монте-Карло (1.1) (см.
алгоритм 1.1).Для алгоритма 2.8 требуется определить начальную плотность и переходную функцию моделируемой прикладной цепи Маркова.В дальнейших рассуждениях выборочные значения точек столкновений (состояний прикладной цепи Маркова) будем снабжать нижниминдексом i, который будет обозначать номер траектории.Каждая точка столкновения описывается трехкомпонентной случайной величиной ρ (с соответствующими индексами) в системе координатr и единичным случайным вектором ω (в системе координат w), определяющим направление прилета малой частицы в точку столкновения ρ.Для описания вектора ω в принципе можно обойтись двумя компонентами (в сферических координатах), но мы для простоты будем использовать три декартовых компоненты w. Таким образом, точка столкновения ξ = (ρ, ω) (снабженная соответствующими индексами) имеетшесть компонент.Особо подчеркнем, что мы будем описывать простейшую односкоростную «столкновительную» модель.Отметим также, что, с одной стороны, описываемый здесь алгоритмпрямого (имитационного) моделирования позволяет учитывать дополнительные эффекты – изменение скоростей и типов частиц, делениефотонов в точках столкновений, слипание фотонов и др.
– с помощьюсоответствующего расширения фазового пространства (а точнее – рандомизации; см. подраздел 3.2).С другой стороны, и в рассматриваемой далее относительно простойситуации нам потребуется определенная изобретательность при реализации и описании численного алгоритма. В частности, будут использованы:– соображения рандомизации (см. подраздел 3.2),– введение дельта-плотностей для описания особенностей распределений (см. подраздел 4.3),– «геометрические» соображения,– стандартные (такие, как метод обратной функции распределения– см.
подраздел 2.3) и специальные методы численного моделированияслучайных величин и векторов (см. разделы 2, 10–13) и др.846.3. Моделирование точки «рождения» и начальногонаправления движения фотона. Случайный изотропныйвектор. Для формирования плотности распределения π(x) (здесьx = (r, w)) начального состояния прикладной цепи Марковаξ (0) = ρ(0) , ω (0) , соответствующего первому столкновению фотона счастицей вещества, следует описать ситуацию «рождения» малой частицы и ее прямолинейный пробег до первого столкновения.Определенные трудности могут возникнуть уже на стадии описания(−1)(0)и численной реализации выборочных значений ρiи ω i точки ρ(−1)«рождения» фотона и начального случайного единичного вектора направления ω (0) .Рассмотрим (для примера) достаточно простой случай, когда источник частиц представляет собой трехмерныйшар B (3,c,R) ⊂ G радиуса(1) (2) (3)R с центром в точке c = c , c , c:nB (3,c,R) = r = r(1) , r(2) , r(3) :qr(1) − c(1)2+ r(2) − c(2)2+ r(3) − c(3)2o≤R ;при этом точка ρ(−1) распределена в этом шаре равномерно.
Предположим также, что вектор ω (0) является изотропным.ОПРЕДЕЛЕНИЕ6.1 (см., например, [9]). Случайный векторξ = ξ (1) , ..., ξ (d) называется изотропным, если точка ω̃ = ξ/kξkравномерно распределена на поверхности d-мерной сферы единичногорадиусав начале координат и не зависит от длиныq с центром22(1)kξk =ξ+ ... + ξ (d) .УТВЕРЖДЕНИЕ 6.1 (см., например, [1]). Если точка ξ (1) , ..., ξ (d)равномерно распределена в d-мерном шаре радиуса Rq22(d,0,R)(1)(d)(1)(d)B= x = x , ..., x: kxk =x+ ... + x≤R(6.7)согласно плотности распределения1при(1)(d)B̄ (d,0,R)=fξ x , ..., x0 иначеx(1) , ..., x(d) ∈ B (d,0,R) ,(6.8)(здесьB̄ (d,0,R) =π d/2 Rd−Γ(d/2 + 1)85(6.9)объем шара B (d,0,R) ), то вектор ξ = ξ (1) , ..., ξ (d) является изотропным и ω = ξ/kξk.В силу утверждения 6.1 и сделанных нами предположений, для мо(−1)(0)делирования выборочных значений ρiи ω i требуется по сути одини тот же алгоритм моделирования случайной точки ξ из утверждения 6.1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
