Доклады Российской академии наук. Математика, информатика, процессы управления, 2021, T. 498, № 1, стр. 55-58

1. Рассматривается односкоростной (для краткости изложения) процесс переноса частиц (квантов излучения), математическая модель которого определяется кинетическим уравнением [1, 2]

Здесь ${\mathbf{r}} \in {{R}^{3}}$, ω – единичный вектор направления скорости, $\Phi ({\mathbf{r}},\omega )$ – плотность потока частиц-квантов; $\rho ({\mathbf{r}})$ – плотность среды, $\rho ({\mathbf{r}})\sigma ({\mathbf{r}})$ – коэффициент ослабления (полное сечение), $\rho ({\mathbf{r}}){{\sigma }_{s}}({\mathbf{r}})$ – сечение рассеяния, ${{\Phi }_{0}}({\mathbf{r}},\omega )$ – плотность распределения частиц в источнике, $w(\omega ,\omega {\kern 1pt} ';{\mathbf{r}})$ – индикатриса рассеяния кванта в точке ${\mathbf{r}}{\kern 1pt} $. Для решения задач переноса численно-статистически моделируется цепь Маркова столкновений частицы с элементами вещества, свободный пробег l между которыми распределен с плотностью f(l) = = $\rho ({\mathbf{r}}(l))\sigma ({\mathbf{r}}(l)){\text{exp}}( - \tau (l))$, где ${\mathbf{r}}(l) = {\mathbf{r}}{\kern 1pt} '\; + l\omega $, а τ(l) = = $\int\limits_0^l \,\rho ({\mathbf{r}}(t))\sigma ({\mathbf{r}}(t))dt$ [3, 2 ]. Если $\rho \sigma \equiv {\text{const}}$, то пробег можно моделировать по формуле $l = \frac{{ - \ln \alpha }}{{\rho \sigma }}$, где α – случайное число, равномерно распределенное в интервале (0, 1). Если же плотность среды существенно меняется, то может быть полезным метод максимального сечения (выравнивания, дельта-рассеяния) [4–6]. Этот метод основан на том, что в уравнении (1) $\rho ({\mathbf{r}})\sigma ({\mathbf{r}})$ заменяется на ${{\sigma }_{{{\text{max}}}}} \geqslant \rho ({\mathbf{r}})\sigma ({\mathbf{r}})$ и в точке столкновения с вероятностью $\frac{{{{\sigma }_{{{\text{max}}}}} - \rho ({\mathbf{r}})\sigma ({\mathbf{r}})}}{{{{\sigma }_{{{\text{max}}}}}}}$ частица не меняет скорость, т.е. реализуется индикатриса $w(\omega ,\omega {\kern 1pt} ';{\mathbf{r}}) \equiv \delta (\omega - \omega {\kern 1pt} ')$. Адекватность метода следует из уравнения (1), а также из известного свойства сохранения пуассоновости при случайном прореживании пуассоновского потока столкновений [6].

2. Рассматриваемый метод применяется для решения стохастических задач теории переноса со случайным полем $\rho ({\mathbf{r}})$, которое для простоты изложения в основном будем предполагать изотропным с нормированной корреляционной функцией k(r) и корреляционной длиной (корреляционным радиусом) $L = \int\limits_0^\infty \,k(r)dr$. Предполагается, что при этом $\sigma ({\mathbf{r}}) = \rho ({\mathbf{r}})\sigma $. Для решения стохастических задач теории переноса эффективен “метод двойной рандомизации” [7], в котором для каждой реализации среды моделируется одна или несколько (для уменьшения трудоемкости) траекторий частиц и строятся несмещенные оценки линейных функционалов, например, осредненной по реализациям среды вероятности прохождения. Трудоемкость (сложность) такого алгоритма может быть очень большой и даже неограниченной при L → 0, например, для “мозаичных” моделей $\rho ({\mathbf{r}})$ (пуассоновских, Вороного и т.п.) (см., например, [8–10]). В связи с этим далее формулируется несложный корреляционно рандомизированный алгоритм (КР-алгоритм) для оценки влияния стохастичности поля $\rho $ на перенос частиц, эвристически основанный на том факте, что вероятность прохождения в значительной степени определяется корреляционной длиной L и одномерным распределением поля $\rho $ [8].

3. КР-алгоритм определяется пунктами S1, S2:

S1) В процессе переноса пробег из ${\mathbf{r}}{\kern 1pt} '$ в ${\mathbf{r}}$ моделируется по формуле $l = - \frac{{ln\alpha }}{{{{\sigma }_{{{\text{max}}}}}}}$;

S2) Если $l = \left| {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '} \right| < L$, то фиксируется ρ(r) = = $\rho ({\mathbf{r}}{\kern 1pt} ')$, иначе значение $\rho ({\mathbf{r}})$ выбирается случайно из одномерного распределения поля $\rho $.

Как показывают приведенные выше эвристические соображения, а также модельные исследования и численные эксперименты, такой алгоритм дает удовлетворительные результаты при малых значениях L, когда стандартное моделирование слишком трудоемко. При этом значение L можно корректировать на основе упрощенных тестовых задач (см. далее пункт 4). Ясно, что трудоемкость КР-алгоритма ограничена при $L \to 0$ и его можно использовать для исследования соответствующей скорости сходимости в стохастических задачах переноса. Отметим, что в случае неизотропного однородного поля $\rho $ можно также использовать КР-алгоритм, заменив L на $L(\omega ),\omega = \frac{{{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '}}{{\left| {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '} \right|}}$.

4. Далее эффективность КР-алгоритма изучается для модельной задачи о переносе частиц через плоский слой: ${\mathbf{H}} = \{ {\mathbf{r}} = (x,y,z){\kern 1pt} :\;0 \leqslant z \leqslant H\} $, заполненный веществом с единичной плотностью ${{\rho }_{{\mathbf{V}}}} = 1$, “разбавленным” случайным ансамблем пустых (ρ = 0) шаров радиуса ${{r}_{b}}$, центры которых $\{ {{{\mathbf{r}}}_{i}}\} $ образуют пуассоновский точечный поток интенсивности $\lambda $. Процесс переноса является предельно анизотропным, т.е. $w(\omega ,\omega {\kern 1pt} ') = \delta (\omega - \omega {\kern 1pt} ')$. Пусть ${\mathbf{V}}$ – случайная непустая подобласть слоя. Предположим, что ${{\sigma }_{{\mathbf{V}}}} = \sigma \equiv 1$, вероятность выживания при столкновении ${{q}_{{\mathbf{V}}}}({\mathbf{r}}) \equiv q = \frac{{{{\sigma }_{s}}}}{\sigma }$ и Φ₀(z, ω) = $\delta (z)\delta (\omega - {{\omega }_{z}})$. Тогда средняя (при ${{\rho }_{{\mathbf{V}}}} = 1$) плотность ρ среды в слое равна вероятности $P({\mathbf{r}} \in {\mathbf{V}})$, $\forall {\mathbf{r}} \in {\mathbf{H}}$. Согласно распределению Пуассона отсюда имеем

Введем обозначение: $\xi (z)$ для однородной случайной функции со значениями: $\xi (z) = 1$ при ${{{\mathbf{r}}}_{z}} = (0,0,z) \in {\mathbf{V}}$ и $\xi (z) = 0$ при ${{{\mathbf{r}}}_{z}} \notin {\mathbf{V}}$ и, следовательно, с одномерным распределением: $P(\xi = 1)$ = ρ, $P(\xi = 0) = 1 - \rho $. При этом длина части отрезка $[{{{\mathbf{r}}}_{0}},{{{\mathbf{r}}}_{H}}]$, находящейся в веществе, равна ζ = $\int\limits_0^H \,\xi (z)dz$. Хорошо известно (см., например, [2]), что соответствующая вероятность прохождения (при ${{\rho }_{{\mathbf{V}}}} = 1$) равна ${{P}_{\zeta }} = exp( - (1 - q)\zeta )$, и, следовательно,

Элементарные геометрические соображения с учетом свойств пуассоновского потока показывают, что корреляционная функция для $\xi (z)$ равна $k(t) = \frac{{{{e}^{{ - \lambda {{V}_{b}}}}}({{e}^{{\lambda {{V}_{t}}}}} - 1)}}{{1 - {{e}^{{ - \lambda {{V}_{b}}}}}}}$, где V_b = $\frac{4}{3}\pi r_{b}^{3}$ и V_t = = $2\tfrac{{\pi {{h}^{2}}}}{3}\left( {3{{r}_{b}} - \frac{{2{{r}_{b}} - t}}{2}} \right)$ – объем общей части пересекающихся шаров радиуса r_b с центрами, отстоящими на расстоянии $t \leqslant 2{{r}_{b}}$.

Перейдем теперь к построению оценок осредненной вероятности P прохождения для рассматриваемой задачи. Согласно сказанному выше, вероятность прохождения через осредненную среду равна ${{P}_{0}} = exp( - (1 - q)\rho H)$. Далее, на основе (3), используя предельную теорему для однородного случайного процесса $\xi (z)$ [11], получаем следующую асимптотическую (при $H{\text{/}}L \to \infty $) оценку: $\tfrac{P}{{{{P}_{0}}}} \asymp \exp \left( {\tfrac{{{{d}^{2}}}}{2}{{{(1 - q)}}^{2}}} \right)$, где ${{d}^{2}} = 2L\rho (1 - \rho )H$. КР-алгоритм с $L = \int\limits_0^\infty \,k(t)dt$ при ${{\sigma }_{{{\text{max}}}}} = 1$ реализуется согласно сказанному выше последовательным моделированием пробегов по формуле $l = - ln\alpha $. Используя весовое моделирование, следует при столкновении в области V вес домножать на q. Результативный вес дает несмещенную статистическую оценку величины P_s.

Верификация оценки P_s осуществлялась сравнением с несмещенной оценкой величины P, получаемой стандартным методом двойной рандомизации при ${{\sigma }_{{{\text{max}}}}} = 1$ с реализацией пуассоновского ансамбля шаров в цилиндре $\{ 0 \leqslant z \leqslant H,{{x}^{2}} + {{y}^{2}} \leqslant r_{b}^{2}\} $. Как и в КР-алгоритме, при столкновении в области ${\mathbf{V}}$ вес частицы домножался на $q$. Трудоемкость при этом соответственно (2) определяется средним числом $\tfrac{{3H}}{{4{{r}_{b}}}}ln\rho $ шаров в цилиндре.

В табл. 1 приведены для H = 200, $q = \rho = 0.9$ значения ${{P}_{{as}}}{\text{/}}{{P}_{0}}$ и статистические оценки $P{\text{/}}{{P}_{0}}$, ${{P}_{s}}{\text{/}}{{P}_{0}}$, среднеквадратические погрешности которых составляют единицу последнего приведенного разряда.

Эти оценки показывают удовлетворительность здесь КР-алгоритма при $L \lesssim 0.4{\text{/}}{{\sigma }_{{{\text{max}}}}}$. Видно также, что при $L \gtrsim 0.4{\text{/}}{{\sigma }_{{{\text{max}}}}}$ в КР-алгоритме параметр L целесообразно уменьшать, например, в случае ${{r}_{b}} = 1$, $L = 0.735$ использовать $L = 0.579$. Отметим, что трудоемкость оценки ${{P}_{s}}$ ограничена при ${{r}_{b}} \to 0$, в то время как для несмещенной оценки P она, согласно сказанному выше, растет до бесконечности. Например, при ${{r}_{b}} = 0.1$ соответствующее отношение составляет примерно100.

Решалась также задача [12] о переносе (в направлении Oz) гамма-квантов с комптоновским рас-сеянием и поглощением через призму (–50 см ≤ x, y  ≤ 50 см, $0 \leqslant z \leqslant 200$ см), заполненную водой, “разбавленной” пуассоновским ансамблем пустых шаров радиуса ${{r}_{b}}$ (см. пункт 4) со средней плотностью среды $\rho = 0.9$ г/см³; начальная энергия кванта ${{E}_{0}} = 1$ Мэв. Полученные КР-алгоритмом оценки P_s для ${{r}_{b}} = 5$ см и ${{r}_{b}} = 2$ см совпали с соответствующими несмещенными оценками P в пределах достаточно малой погрешности. Это объясняется тем, что, как оказалось, в этой задаче практически выполняется соотношение $L \leqslant 0.4{\text{/}}{{\sigma }_{{{\text{max}}}}}(E)$. Расчеты для меньших значений ${{r}_{b}}$ продемонстрировали сходимость P_s (и тем самым P) к ${{P}_{0}}$ при ${{r}_{b}} \to 0$, которую доказать теоретически здесь не представляется возможным. Отношение трудоемкостей несмещенной оценки P и оценки P_s при этом оказалось близким к ${{(10r_{b}^{{ - 1}})}^{3}}$. Ясно, что КР-алгоритм можно использовать для аналогичных расчетов в неограниченных областях.