1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 16
Текст из файла (страница 16)
также формулу(13.6) из раздела 13 данного пособия); при этом g(w|x0 ) ≡ 1/(4π).В случае неизотропного рассеяния пересчет координат направленияпробега осуществляется по формулам, получающимся из специальныхгеометрических соображений (см., например, [5, 9, 20]):(m,1)(m−1,1) (m)(m−1,2)(m)ωi= ωiµi − ωisin θi +(m,3)ωiv(m) 2uu(m−1,1) (m−1,3)(m) t 1 − µi,+ωiωicos θi(m−1,3) 21 − ωi(m,2)(m−1,2) (m)(m−1,1)(m)ωi= ωiµi + ωisin θi −v(m) 2uu(m−1,2) (m−1,3)(m) t 1 − µi,−ωiωicos θi(m−1,3) 21 − ωivu(m) 2hiu(m−1,3) (m)(m−1,3) 2(m) t 1 − µi= ωiµi + 1 − ωi.cos θi(m−1,3) 21 − ωiЗдесь µ(m) – косинус угла между ω (m−1) и ω (m) ; а θ(m) – азимутальный угол рассеяния, который в данном случае отсчитывается отплоскости, проходящейчерез ω (m−1) и Oz; т.
е. θ – угол между плоско(m−1)стями ω; Oz и ω (m−1) ; ω (m−1) .Угол θ(m) распределен равномерно в интервале (0, 2π) и его можно(m)(m)моделировать по=2παi ; при этом пара формуле θi(m)(m)cos θ , sin θявляется двумерным изотропным вектором (см. соотношения (6.29), (6.30) при ρ0 = 1, а также формулу (13.3) из раздела 13 данного пособия).Для моделирования случайных величин µ(m) чаще всего используется индикатриса Хеньи – Гринстейна (см., например, [9]), представляющая собой элементарную плотность видаfµ(m) (u) =1 − c203/22 (1 + c20 − 2 c0 u)99при u, c0 ∈ (−1, +1).(6.44)Эта плотность подобрана таким образом, чтоZ 1Eµ(m) =u fµ(m) (u) du = c0 .−1Для моделирования «рассеяния вперед» принимают c0 .
1, а для«рассеяния назад» берут c0 & −1.(m)Моделирующую формулу для выборочных значений µiполучаем, последовательно решая уравнение (2.16) метода обратной функциираспределения для плотности (6.44) [24]:(m)µiZ−1Zили(m)µi1 − c20 du3/22 (1 + c20 − 2 c0 u)= αi ,1 − c20 d 1 + c20 − 2 c0 u= αi ,3/2−4c0 (1 + c20 − 2 c0 u)(m)2Z1 − c20 1+c0 −2c0 µi−v −3/2 dv = αi ,4c0(1+c0 )2−1или1 − c20или2c0Z(m)1+c20 −2c0 µidv −1/2 = αi ,(1+c0 )2112c0 αi,или q−=1+c1− c20(m)01 + c20 − 2c0 µi12c0 αi + 1 − c0,или q=1 − c20(m)1 + c20 − 2c0 µi"2 #11 − c20(m)2и, окончательно, µi =1 + c0 −.2 c02 c0 αi + 1 − c0(m)Проверка 2.1 для αi = 0 дает µi = 1/(2c0 )× 1+c20 − (1−c20 )/(2c0 ×2 ×0 + 1 − c0 )= 1/(2c0 ) × 1 + c20 − 1 − c20 − 2c0 = −1, а для αi = 12 (m)имеем µi = 1/(2c0 ) × 1 + c20 − (1 − c20 )/(2c0 × 1 + 1 − c0 )= 1/(2c0 )×× 1 + c20 − 1 − c20 + 2c0 = 1.(m)После того, как направление ω iвыбрано, происходит моделирование длины свободного пробега по формуле (6.42).100Используя соображения о рандомизации (увеличении числа случайных параметров) из подраздела 3.2, с помощью формул вида (2.4)–(2.7),(2.10) для совместной плотности компонент случайного вектора получаем следующую переходную функцию p(x0 , x) из соотношения (6.3):0r − r0 Σe−Σ|r−r |×.p(x0 , x) = p(s) × g(w|w0 ) × δ w −|r − r0 ||r − r0 |2(6.45)Особо отметим, что при p(s) < 1 для функции (6.45) выполнено соотношение (6.4), а также то, что функция (6.45) содержит интегрируемуюособенность в виде дельта-функции по направлению, что подтверждает целесообразность рассмотрения соответствующих плотностей состояний (6.5) как функций из L1 (X) (см.
рассуждения из подраздела 6.1).6.6. Простейшая схема прямого численного моделирования.На основании проведенного анализа прикладной цепи Маркова, описывающей перенос частиц, можно приближенно вычислять на ЭВМ значения нужных характеристик (функционалов) этого процесса.По предложению, сформулированному в подразделе 6.2, в качествеосновного примера возьмем вычисление вероятности P того, что малаячастица поглотится внутри области G.Введем случайную величину η̂, которая равна единице в случае, если частица излучения поглотилась в области G, и равна нулю во всехостальных случаях. Заметим, что Eη̂ = P .Рассмотрим алгоритм прямого численного моделирования для приближения этой величины.АЛГОРИТМ 6.2.
(см., например, [9, 13]). Численно реализуем n траекторий частиц и для i-й траектории (i = 1, . . . , n) делаем следующее.1. Полагаем m := 0 и моделируем начальную координату точки(m−1)(−1)«рождения» частицы ρi= ρiи начальное направление движе(m)(0)ния частицы ω i = ω i по формулам (6.10), (6.34), (6.35).(m)2.
Моделируем длину свободного пробега λi частицы по формуле(6.42).3. Проверяем, не вылетела ли частица из области G. Если вылетпроизошел, то заканчиваем траекторию и полагаем η̂i = 0.(m)4. Вычисляем координату r = ρiочередной точки столкнове(m−1)0ния, зная предыдущую координату r = ρi, направление движения(m)(m)w0 = ω i и длину свободного пробега λi :(m)ρi(m−1)= ρi(m)+ λi101(m)ωi.5.
Определяем тип столкновения (поглощение или рассеяние) согласно вероятности p(a) и алгоритму 2.7. Если реализовалось поглощение, то обрываем траекторию и полагаем η̂i = 1.6. Если происходит рассеяние, то полагаем m := m + 1, а затемсогласно индикатрисе(m−1) g w w 0 = ω i(m)выбираем новое направление движения частицы w = ω i и переходимна п.
2 для дальнейшего моделирования траектории.Приближением искомой вероятности P поглощения частицы в области G будет среднее арифметическое вида (1.1):P = Eη̂ ≈η̂1 + . . . + η̂n.nДля дальнейших рассуждений важно заметить, чтоh ξ (N )η̂ =для h(x) = p(a) × χ(G) (r);(6.46)p(a)здесь ξ (N ) = ρ(N ) , ω (N ) – случайная точка столкновения, в которойпроисходит поглощение моделируемой частицы (при этом N – случайный номер обрыва траектории частицы), а χ(G) (r) – индикатор множества G:χ(G) (r) = 1 при r ∈ G и χ(G) (r) = 0 иначе.В теории методов Монте-Карло (во многом благодаря приводимомуздесь примеру блуждания малых частиц) случайная величина (6.46)называется оценивателем (монте-карловской оценкой) по поглощениям для величины P (см., например, [5, 9], а также подраздел 7.5 данногопособия).1027.
Марковское интегральное уравнение.Линейные функционалы от решенияинтегрального уравнения Фредгольмавторого рода. Основной весовой оцениватель.Оцениватель по поглощениям7.1. Интегральное уравнение для суммарной плотностистолкновений. Для расширения спектра оценивателей (монте-карловских оценок) типа (6.46) необходимо учесть следующие важные соображения.Вернемся к материалам подраздела 6.1. Рассмотрим прикладнуюцепь Маркова (6.2) с начальной плотностью π(x) и переходной функцией (6.3) и связанные с ней соотношения (6.5), (6.6), в которых функцииψ, ψ (0) , ψ (1) , ... принадлежат пространству L1 (X); X ⊆ Rd .Индукцией по m из этих соотношений несложно получить, чтоψ(x) = π(x) + K̂π(x) + .
. . + K̂ m π(x) + ... =∞XK̂ m π(x),(7.1)m=0гдеmZK̂ π(x) =Z...Xπ y (0) p y (0) , y (1) × .. × p y (m−1) , x dy (0) ..dy (m−1) .X(7.2)Формулу (7.2) можно пояснить следующим образом. Аналог совместной(2.27) состояний прикладной цепи Маркова (m−1) плотности(m−1)(m)ξ̂,ξ, где ξ̂= ξ (0) , ξ (1) , . . . , ξ (m−1) , имеет видψ̂ˆ (m−1)ξ,ξ (m) y (0) , y (1) , .., y (m) = π y (0) p y (0) , y (1) ×..×p y (m−1) , y (m) .Тогда, согласно формуле (2.5), «безусловная» плотностьψ (m) (x) = K̂ψ (m−1) (x) = K̂ m π(x) компоненты ξ (m) имеет вид (7.2).При выполнении условия (6.4) интегральный оператор K̂ являетсясжимающим в пространстве L1 (X), а функция (7.1) представляет собой функциональный ряд Неймана, сходящийся к единственномурешению интегрального уравнения Фредгольма второго родаZψ(x) =p(x0 , x)ψ(x0 ) dx0 + π(x) или ψ = K̂ψ + π;(7.3)X103здесь p(x0 , x) – ядро интегрального оператора K̂, а π(x) – свободныйчлен уравнения (7.3) (см., например, [19]).То, что функция (7.1) является решением уравнения (7.3), подтверждается непосредственной подстановкой:!∞∞∞XXXK̂ψ + π = K̂K̂ m π + π =K̂ m π + π =K̂ m π = ψ.
(7.4)m=0m=1m=0Уравнение (7.3), соответствующее прикладной цепи Маркова (6.2)будем называть марковским интегральным уравнением.Практически важные величины, связанные с прикладной цепью Маркова (6.2), удается записать в виде линейного функционала от решения уравнения (7.3) вида:ZdefIˆh = (ψ, h) =ψ(x)h(x) dx,(7.5)Xгде h(x) – заданная, определяющая функционал (7.5) функция из пространства L∞ (X), сопряженного пространству L1 (X) с нормойkhkL∞ (X) = vrai supx∈X |h(x)|.В том числе можно представить исследованную в разделе 6 величинуP (вероятность того, что малая частица поглотится внутри области G)в виде функционала вида (7.5): P = Iˆh = (ψ, h), где ψ(x) – решениемарковского интегрального уравнения (7.3) с ядром (6.45) и свободнымчленом (6.43), а функция h(x) имеет вид (6.46).Заметим, что обозначение (ψ, h) из соотношения (7.5) «похоже» наопределение скалярного произведения в гильбертовом функциональномпространстве L2 (X) (см., например, [19]).
Однако использование термина «скалярное произведение» для дальнейших рассуждений некорректно, так как пространство L1 (X) не является гильбертовым.7.2. Произвольное интегральное уравнение Фредгольмавторого рода. Линейный функционал от его решения каксумма интегралов бесконечно возрастающей кратности. Приведенная здесь теория марковских интегральных уравнений (7.3) позволяет строить весовые оцениватели для линейных функционалов вида(7.5)ZdefIh = (ϕ, h) =ϕ(x)h(x) dx(7.6)X(для заданной функции h(x) ∈ L∞ (X)) от неизвестного решения104ϕ(x) ∈ L1 (X) общего (не обязательно марковского) интегрального уравнения Фредгольма второго рода вида (7.3)Zϕ(x) =k(x0 , x)ϕ(x0 ) dx0 + f (x) или ϕ = Kϕ + f ;(7.7)Xздесь k(x0 , x) – ядро интегрального оператора K, а f (x) ∈ L1 (X) –свободный член уравнения (7.7).Наличие практически важных примеров уравнений (7.7), не являющихся марковскими, подтверждает, в частности, следующий пример.ПРИМЕР 7.1 (см., например, [9]).
Рассмотрим трехмерную задачуДирихле для уравнения Гельмгольца (см., например, [25]):∆ϕ(x) − cϕ(x) = −g(x), ϕ(x)|Γ = ψ(x)(7.8)в области x ∈ D ⊂ R3 с границей Γ; здесь c = const > 0, а ∆ϕ(x) =∂ 2 ϕ(x))∂ 2 ϕ(x))∂ 2 ϕ(x))(1), x(2) , x(3) .2 +2 +2 – оператор Лапласа; x = x(∂x(1) )(∂x(2) )(∂x(3) )Предполагаются выполненными условия регулярности заданныхфункций g(x), ψ(x) и границы Γ, обеспечивающие существование и единственность решения ϕ(x) краевой задачи (7.8), а также возможность егопредставления с помощью функции Грина для шара (см., например,[25]).Для функции ϕ(x) можно записать интегральное уравнение Фредгольма второго рода (7.7) с ядром√(cd(x)√× δx (x0 ), x ∈/ Γε ;0sinh( cd(x))k(x , x) =(7.9)0, x ∈ Γεи свободным членом(R sinh(√c(d(x)−|x0 −x|))100√/ Γε ;4π D |x0 −x| sinh( cd(x)) g(x ) dx , x ∈f (x) =ϕ(x), x ∈ Γε .Здесь d(x) – расстояние от точки x до границы Γ, а Γε – ε-окрестностьграницы Γ, т.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
