1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 15
Текст из файла (страница 15)
(6.28)fϕ (φ|r) =2fρ (r)2π0 πRСледовательно, плотность (6.26) представима в видеf(ρ,ϕ) (r, φ) = fρ (r) × fϕ (φ).91Это означает, что компоненты случайного вектора (ρ, ϕ), определяемого соотношением (6.25), независимы и распределены согласно элементарным плотностям fρ (r) и fϕ (φ) из соотношений (6.27), (6.28) соответственно. Решая соответствующие уравнения (2.16) метода обратнойфункции распределения, получаем моделирующие формулы√(6.29)ρ0 = R α1 , ϕ0 = 2πα2 ,где α1 , α2 ∈ U (0, 1) – стандартные случайные числа.С учетом соотношений (6.24) и (6.29) получаем моделирующие формулы для выборочных значений компонент двумерного случайного вектора ξ, имеющего плотность распределения (6.12):√√(2)(1)ξ0 = R α1 sin 2πα2 , ξ0 = R α1 cos 2πα2 .(6.30)Таким образом, мы показали, что проблема 6.1 для случая d = 2разрешима, а также то, что, согласно определению 6.1, преобразованиедекартовых координат {u, v} к полярным координатам {r, φ} является вычислимым моделируемым для двумерного случайного вектора ξ сраспределением (6.12).Теперь исследуем возможности решения проблемы 6.1 для случаяd = 3.Рассмотрим случайную точку ξ = ξ (1) , ξ (2) , ξ (3) , имеющую плотность распределения (6.11).Выкладки, аналогичные (6.13), (6.14), для плотности (6.11) подтверждают замечание 6.1 о том, что прямое применение алгоритма 2.3 длямоделирования случайной точки, равномерно распределенной в трехмерном шаре B (3,0,R) , не дает экономичных моделирующих формул метода обратной функции распределения для всех компонент этой точки.По аналогии с двумерным случаем, применим утверждение 6.2 ирассмотрим взаимно однозначное преобразование Φ(u, v, w), описывающее переход от декартовых координат {u, v, w} к сферическим координатам {r, φ, θ}.
Обратное к Φ(u, v, w) преобразование определяетсяследующим частным случаем формул (6.19):u = r sin φ sin θ, v = r cos φ sin θ, w = r cos θ; r ≥ 0, 0 ≤ φ < 2π, 0 ≤ θ < π.(6.31)Якобиан преобразования (6.31) равен ∂ (u(r, φ, θ), v(r, φ, θ), w(r, φ, θ)) = r2 sin θ∂(r, φ, θ)92(см. формулу (6.20)) и, согласно утверждению 6.2, плотность случайного вектораρ, ϕ, θ̃ = Φ ξ (1) , ξ (2) , ξ (3)(6.32)(здесь вектор ξ = ξ (1) , ξ (2) , ξ (3) имеет плотность распределения (6.11))имеет вид3 r2 sin θ; 0 ≤ r < R, 0 ≤ φ < 2π, 0 ≤ θ < π.f(ρ,ϕ,θ̃) (r, φ, θ) =4 π R3Реализуя выкладки типа (6.27), (6.28), получаем аналог соотношения (6.23): 2 3r1sin θf(ρ,ϕ,θ̃) (r, φ, θ) = fρ (r) × fϕ (φ) × fθ̃ (θ) =××.R32π2 (6.33)Это означает, что компоненты случайного вектора ρ, ϕ, θ̃ , определяемого соотношением (6.32), независимы и распределены согласно элементарным плотностям fρ (r), fϕ (φ) и fθ̃ (θ) из соотношения (6.33). Решая соответствующие уравнения (2.16) метода обратной функции распределения, получаем моделирующие формулы√ρ0 = R 3 α1 , ϕ0 = 2πα2 , cos θ̃0 = 1 − 2α3 ,(6.34)где α1 , α2 , α3 ∈ U (0, 1) – стандартные случайные числа.Декартовы координаты выборочного значения ξ 0 случайного вектора ξ, имеющего плотность распределения (6.11), равны(1)(2)(3)ξ0 = ρ0 sin ϕ0 sin θ̃0 , ξ0 = ρ0 cos ϕ0 sin θ̃0 , ξ0 = ρ0 cos θ̃0 .(6.35)Формулы (6.34), (6.35) можно использовать при получении выбо(−1)рочных значений точки «рождения» малой частицы ρiи начального(0)направления ее движения ω i по формулам (6.10).Таким образом, мы показали, что проблема 6.1 для случая d = 3 разрешима, а также то, что, согласно определению 6.1, преобразование декартовых координат {u, v, w} к сферическим координатам {r, φ, θ} является вычислимым моделируемым для трехмерного случайного вектора ξ с распределением (6.11).Однако, начиная с d = 4, уже не все частные плотности распределения из разложения (6.23) являются элементарными: только три первые плотности дают моделирующие формулы метода обратной функции93распределения, полученные из уравнений типа (2.16):√(1)(2)ρ0 = R d α1 , ϕ0 = 2πα2 , ϕ0 = arccos(1 − 2α3 )(см.
аналог этих соотношений – формулы (6.34) – для случая d = 3),а вот для следующих компонент таких формул вывести не удается; вчастности, для компоненты ϕ(3) уравнение (2.16) имеет вид2πZ(3)ϕ0(3)(3)sin2 u du = α4 или 2ϕ0 − sin 2ϕ0 = 2πα4 ,0где α4 ∈ U (0, 1). Последнее уравнение не разрешимо в элементарных(3)функциях относительно ϕ0 .Таким образом,преобразование Φ(x) декартовыхкоординат (1)x , . . .
, x(d) к гиперсферическим координатам r, φ(1) , . . . , φ(d−1) неявляется вычислимым моделируемым для случайного вектора ξ размерности d ≥ 4 с распределением (6.8), и проблема 6.1 для этого случаяосталась неразрешимой.Тем не менее, далее в подразделе 13.2 будет представлено решениепроблемы 6.1 для случая d ≥ 4, основанное на специальных свойствахстандартного гауссовского распределения (1.12).В качестве основного вывода этого подраздела заметим, что уже настадии розыгрыша случайной точки «рождения» малой частицы (фо(−1)(0)тона) ρiи начального направления ее движения ω i – даже дляпростейшей геометрии источника и изотропии начального направления– получились некоторые специальные (отличные от стандартного алгоритма 2.3) алгоритмы численного моделирования случайных векторов.6.5.
Моделирование длины свободного пробега и точекстолкновения фотона. После «рождения» малая частица движется(−1)(0)прямолинейно от точки ρiвдоль направления ω i до точки первого(0)столкновения ρi .Для фиксации прямолинейности этого движения в плотность распределения первого столкновения должен быть включен соответствующий «ограничитель возможностей» (см.
подраздел 4.3) – множительr − r0 δ w−–|r − r0 |дельта-функция по направлению.94(6.36)Наконец, должна быть разыграна случайная величинаλ(m) = ρ(m) − ρ(m−1) , m = 0,отражающая длину свободного пробега фотона (т. е. длину отрезка прямолинейного пути частицы, на котором она не испытывает столкновений).При этом используется следующее предположение: вероятность того, что частица, вылетевшая из точки s, испытает столкновение в интервале (s, s+∆s) (ось s совпадает с направлением движения частицы),равнаΣ ρ(m−1) + ω (m) s ∆s + o(∆s).(6.37)Функция Σ(ρ) (ρ обозначает координату частицы во «внешней» системе координат) предполагается заданной и называется полным сечением взаимодействия частицы со средой.Из соотношения (6.37) получаетсяследующее рассуждение [5].
Обозначим через Fλ(m) (u) = P λ(m) < u функцию распределения длинысвободного пробега.Вероятность того, что частица испытает первое столкновение в интервале (u, u + ∆u) равнаhiFλ(m) (u+∆u)−Fλ(m) (u) = [1−Fλ(m) (u)] Σ ρ(m−1) + ω (m) u ∆u + o(∆u) ;здесь 1 − Fλ(m) (u) – это вероятность того, что частица долетит до точкиu. Разделив последнее соотношение на ∆u и перейдя к пределу при∆u → 0, получим дифференциальное уравнениеdFλ(m) (u)= Σ ρ(m−1) + ω (m) u [1 − Fλ(m) (u)]duс начальным условием Fλ(m) (0) = 0.Решением этого уравнения является функция обобщенного экспоненциального распределения Z uFλ(m) (u) = 1 − exp −Σ ρ(m−1) + ω (m) s ds ;0соответствующая плотность распределения равна Z ufλ(m) (u) = Σ ρ(m−1) + ω (m) u × exp −Σ ρ(m−1) + ω (m) s ds .0(6.38)95Весьма интересным является вопрос о численном моделировании выборочных значений случайной величины, имеющей плотность распределения (6.38).
В частности, имеется оригинальный алгоритм, которыйносит название метод мажорантного (максимального) сечения, идеякоторого состоит в следующем (см., например, [9, 23]).Полагаем, что существует мажоранта функции сеченияΣ ρ(m−1) + ω (m) s ≤ σm (s),(6.39)такая, что численное моделирование выборочных значений ηi случайной величины η согласно плотности Z uσm (s) ds , u > 0,(6.40)fη (u) = σm (u) exp −0реализуется достаточно просто.АЛГОРИТМ 6.1 (см., например, [9, 23]).
1. Полагаем сначала i := 1;L := 0.2. Моделируем выборочное значение ηi случайной величины η, распределенной согласно плотности (6.40), а также стандартное случайное число αi . Вычисляем L := L + ηi .3. ЕслиΣ ρ(m−1) + ω (m) L,(6.41)αi <σm (L)(m)то λ0 = L, иначе полагаем i := i + 1 и повторяем пункт 2 и т. д.Напомним, что знак «:=» означает переприсваивание, т. е.
приA := B числа или результаты вычислений в правой части B помещаются в ячейку памяти ЭВМ с именем A.Обоснование алгоритма 6.1 связано со свойством прореживания мо(m)делируемого в этом алгоритме пуассоновского потока τi = η1 + ... + ηiинтенсивности σm (s).(m)УТВЕРЖДЕНИЕ 6.3 [23].
Если точки τi исключать условно неза(m)(m) висимо с вероятностями 1 − σ τi/σm τi(т. е. принимать с ве(m) (m) роятностями p(τi ) = σ τi/σm τi), то отобранные точки τi реализуют пуассоновский поток интенсивности σ(s).Доказательство этого утверждения основано на том, что, с точностью до рассмотрения величин порядка o(ds), вероятность появленияотобранной точки в интервале (s, s + ds) равна p(s) × σm (s) ds = σ(s) ds[23].96Моделирование обобщенного экспоненциального распределения(6.38) будет эффективным, если в алгоритме 6.1 количество проверокнеравенства (6.41) будет относительно мало. Для этого мажоранта σm (u)из соотношения (6.39) должна быть близка к σ(u).В частности, достаточно эффективным может оказаться выбор в качестве σm (u) кусочно-постоянногоприближения «сверху» для функцииΣ ρ(m−1) + ω (m) u . В этом случае функция (6.40) оказывается составной плотностью с частными плотностями усеченного экспоненциального распределения.
Для моделирования выборочных значений ηi с такимраспределением можно использовать соображения, представленные далее в подразделе 11.3 данного пособия.Здесь мы ограничимся рассмотрением простейшего случая, предполагая, что вещество, заполняющее область G, однородно, что означаетвыполнение равенства Σ(ρ) ≡ Σ = const.В этом случае равенство (6.38) определяет плотность экспоненциального (или показательного) распределения (2.18) для λ = Σ, моделирующая формула для которой имеет вид (2.20):(m)λi=−ln αm,iΣ(6.42)для m = 0.Используя соображения о рандомизации (увеличении числа случайных параметров) из подраздела 3.2, с помощью формул вида (2.4)–(2.7), (2.10) для совместной плотности компонент случайного вектора получаем, что плотность распределения вектора ρ(−1) , ξ (0) (здесьξ (0) = (ρ(0) , ω (0) )) равнаf(ρ(−1) ,ξ(0) ) (r0 , x) = fρ(−1) (r0 ) × fξ(0) xρ(−1) = r0 =011r − r0 Σe−Σ|r−r |××δw−×=(4/3)πR34π|r − r0 ||r − r0 |203Σr − r0 =×δw−× e−Σ|r−r | ; x = (r, w); w ∈ S (3,0,1) ,16π 2 R3 |r − r0 |2|r − r0 |=где S (3,0,1) – трехмерная единичная сфера с центром в начале координат(здесь использовано то, что площадь этой сферы равна 4π).
Появлениев знаменателе множителя |r − r0 |2 связано с соответствующей нормировкой вдоль направления движения фотона.97Отсюда плотность распределения первого («нулевого») столкновения ξ (0) = ρ(0) , ω (0) (это начальное состояние соответствующей прикладной цепи Маркова (6.2)) имеет видZπ(x) = f(ρ(−1) ,ξ(0) ) (r0 , x) dr0 =3Σ=16π 2 R3Z0r − r0 e−Σ|r−r | 0δ w−×dr ; x = (r, w).|r − r0 ||r − r0 |2(6.43)Получив распределение первого («нулевого»)столкновения,попытаемся описать следующие состояния ξ (m) = ρ(m) , ω (m) ; m = 1, ..., Nстолкновений фотонов с частицами вещества.Определим переходную функцию p(x0 , x) соответствующей прикладной цепи Маркова (6.2).Для этого проанализируем то, как происходит переход из (m − 1)-госостояния в m-е.(m−1)Во-первых, для реализованного состояния ξiопределяется, непроизойдет ли в этой точке поглощения (обрыва) траектории.
Это моделируется следующим образом. Полное сечение рассеяния разлагаетсяв суммуΣ(r) = Σ(a) (r) + Σ(s) (r),где Σ(a) (r) – сечение поглощения, Σ(s) (r) – сечение рассеяния; все этинеотрицательные функции предполагаются заданными.Величина p(a) (r) = Σ(a) (r)/Σ(r) описывает вероятность того, чтопри столкновении частицы излучения с частицей вещества в точке r ∈ Gпроисходит поглощение фотона частицей вещества (при этом траектория движения частицы излучения прерывается). Для r ∈ G(−) (здесьG(−) – внешность области G) полагаем p(a) (r) ≡ 1. Соответственно,величина p(s) (r) = Σ(s) (r)/Σ(r) = 1 − p(a) (r) – это вероятность «выживания» фотона.Розыгрыш события поглощения происходит согласно алгоритму 2.7.Заметим, что в рассматриваемом нами случае однородной рассеивающей и поглощающей среды величины Σ(r), Σ(a) (r), Σ(s) (r),p(a) (r), p(s) (r) являются константами.Если реализовалось рассеяние, то происходит численное моделиро(m)вание выборочного значения ω i нового направления движения фотона согласно условной плотности распределения g(w|x0 ) (для однороднойсреды – g(w|w0 )), которую называют индикатрисой рассеяния.98Если рассеяние изотропное, то вне зависимости от x0 (или w0 ) можноиспользовать формулы (6.34), (6.35) для ρ0 = 1 (см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
