Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 490, № 1, стр. 47-50

1. Рассматривается для простоты изложения односкоростной нестационарный процесс переноса частиц с изотропным рассеянием. В качестве математической модели процесса используется цепь Маркова {x_n}, n = 0,1, …, N, соответствующих “столкновений” в фазовом пространстве X = R × V × T координат, скоростей (с фиксированным значением ${v}$ модуля) и времени, т.е. x_n = = (r_n, v_n, t_n), причем v_n = $\frac{{{v}({{{\mathbf{r}}}_{n}} - {{{\mathbf{r}}}_{{n - 1}}})}}{{\left| {{{{\mathbf{r}}}_{n}}\, - \,{{{\mathbf{r}}}_{{n - 1}}}} \right|}}$, ${{t}_{n}}\, = \,{{t}_{{n - 1}}}\, + \,\frac{{\left| {{{{\mathbf{r}}}_{n}} - {{{\mathbf{r}}}_{{n - 1}}}} \right|}}{{v}}$. Эта модель определяется плотностью f(x) точки x₀ и плотностью k(x′, x) перехода из состояния x′ в x, причем

и тем самым среднее число переходов конечно, например, в ограниченной системе [1]. Заданы σ_s(x) и σ_f(x) – коэффициенты рассеяния и деления, среднее число ν частиц в точке деления, σ_c(x) – коэффициент поглощения. Полный коэффициент ослабления потока частиц σ(r) = σ_s(r) + σ_f(r) + + σ_c(r) (см. [1]).

Методы Монте-Карло используют для оценки линейных функционалов вида J_h = (φ, h), $h \in {{L}_{\infty }}$, где φ(x) – плотность столкновений, удовлетворяющая уравнению φ = Kφ + f, в котором K – интегральный оператор с ядром k(x′, x). При выполнении условия (1) имеем ${{\left\| K \right\|}_{{{{L}_{1}}}}} < 1$ и, следовательно, спектральный радиус оператора $\rho (K) < 1$ (см. [1]). Для построения весовой модификации алгоритма используют цепь Маркова с начальной плотностью f₀(x) и плотностью перехода p(x', x), вспомогательные веса

и “оценку по столкновениям” $\xi = \sum\limits_{n = 0}^N {} {{Q}_{n}}h({{x}_{n}})$. Если выполняются “условия несмещенности” [1], то Eξ = (φ, h). Если при этом ρ(K_p) < 1, где K_p – оператор с ядром k²(x', x)/p(x', x), и ${{f}^{2}}{\text{/}}{{f}_{0}} \in {{L}_{1}}(X)$, то ${\text{D}}\xi < + \infty $ [1].

2. Пусть ${{\varphi }_{0}}(x;{{{\mathbf{r}}}_{0}},{{{\mathbf{v}}}_{0}})$ – плотность столкновений (по аргументу x) от одного столкновения в точке $({{{\mathbf{r}}}_{0}},{{{\mathbf{v}}}_{0}},0)$, т.е. для $f(x) = \delta ({\mathbf{r}} - {{{\mathbf{r}}}_{0}})\delta ({\mathbf{v}} - {{{\mathbf{v}}}_{0}})\delta (t)$.

Функционал

можно представить [2] в виде

J(t) = $\int\limits_R \int\limits_V \int\limits_0^\infty f({{{\mathbf{r}}}_{0}},{{{\mathbf{v}}}_{0}}$, $t - \tau )F({{{\mathbf{r}}}_{0}},{{{\mathbf{v}}}_{0}},\tau ){\text{d}}{{{\mathbf{r}}}_{0}}{\text{d}}{{{\mathbf{v}}}_{0}}{\text{d}}\tau $,

где F(r₀, v₀, t) = $\int\limits_R \int\limits_V {{\varphi }_{0}}({\mathbf{r}},{\mathbf{v}},t$; r₀, v₀)h(r, v)drdv.

Предполагается, что f(r, v, –t) ≡ 0, и вследствие этого F(r, v, –t) ≡ 0 при t > 0.

Символом $f_{t}^{{(m)}}$ далее обозначается m-кратная производная от функции f по t. В работе [2] фактически доказана следующая

Лемма 1. Пусть точка (r₀, v₀) распределена для t₀ ≡ 0 с плотностью f₀(r, v), $\left| {\frac{{{{f}^{{(m)}}}({\mathbf{r}},{\mathbf{v}},t)}}{{{{f}_{0}}({\mathbf{r}},{\mathbf{v}})}}} \right|$ < C < +∞, m = 0, 1, …, n, $\rho ({{K}_{p}}) < 1$ и $F(x) < C < + \infty $. Тогда выполняется соотношение ${{J}^{{(m)}}}(t) = {\text{E}}\xi _{t}^{{(m)}}$, где

причем ${\text{D}}\xi _{t}^{{(m)}} < + \infty $, m = 0, 1, …, n.

3. Известно, что при выполнении довольно общих условий в случае источника, локализованного в точке (r₀, v₀, 0), имеет место асимптотическое при t → ∞ соотношение [3, 4]:

где λ – ведущее характеристическое число соответствующего однородного стационарного уравнения переноса с заменой ${{\sigma }_{c}} \mapsto {{\sigma }_{c}} + \frac{\lambda }{{\left| {\mathbf{v}} \right|}}$. Эти условия, в частности, имеют место для односкоростного процесса переноса частиц в ограниченной среде с достаточно быстро убывающей по времени плотностью источника.

Теорема 1. Если выполняются соотношения

соотношение (2) и условия леммы 1, то

Доказательство. Прямое интегрирование с учетом соотношений (2), (3) дает равенства

причем $\varepsilon (t)$, ${{\varepsilon }_{1}}(t) \to 0$ для $t \to + \infty $. Интегрируя функцию $J{\text{'}}(t)$ в пределах $(\tau , + \infty )$ при $\tau \to \infty $ для λ < 0, получаем $J{\text{'}}(t) = C\lambda {{e}^{{\lambda t}}}[1 + {{\varepsilon }_{1}}(t)]$, т.е. ${{C}_{1}} = \lambda C$. В случае λ > 0 тот же результат получается путем введения дополнительного поглощения с коэффициентом $\sigma _{c}^{{(0)}} > \frac{\lambda }{{\left| {\mathbf{v}} \right|}}$. Это завершает доказательство теоремы 1.

4. Определяемую леммой 1 и теоремой 1 оценку метода Монте-Карло можно рандомизировать с целью изучения флуктуаций значения λ для процесса переноса частиц в случайной среде [2]. Будем полагать, что σ(r) – случайное поле, причем отношения $\frac{{{{\sigma }_{s}}}}{\sigma }$, $\frac{{{{\sigma }_{f}}}}{\sigma }$, а также индикатрисы рассеяния  и деления фиксированы. При этом J(t) := J(t, σ), J'(t) := J'(t, σ) и $\lambda (\sigma ) \approx \frac{{J{\text{'}}(t,\sigma )}}{{J(t,\sigma )}}$ соответственно теореме 1. В дальнейшем аргумент t будет опускаться, т.е. $J(t,\sigma ) \mapsto J(\sigma )$.

Практически важными являются величины Eλ(σ) и Dλ(σ). Логически простейший (прямой) способ их оценки методом Монте-Карло состоит в реализации достаточно точных оценок величины $\frac{{J{\text{'}}(t,\sigma )}}{{J(t,\sigma )}}$ для выборки σ достаточно большого объема. Однако такой способ может быть слишком трудоемким для реалистических моделей среды и процесса переноса [5]. Поэтому в настоящей работе для построения рандомизированного алгоритма используется соотношение

где J₀ ≈ EJ(σ).

Простейшая (элементарная) несмещенная рандомизированная оценка величины λ_n на основе леммы 1 строится путем реализации n + 1 условно независимых траекторий частиц для выбранного значения σ: ${{\lambda }_{n}} = {\text{E}}{{\tilde {\lambda }}_{n}}$, где

Для осреднения произведений из (5) целесообразно использовать несовпадающие сочетания по s элементов из группы ${{\Omega }_{1}}, \ldots ,{{\Omega }_{{i - 1}}},{{\Omega }_{{i + 1}}}, \ldots ,{{\Omega }_{{n + 1}}}$ после подстановки ${{\Omega }_{i}} \to {{\Omega }_{{n + 1}}}$ в (5), в связи с тем, что в качестве выделенной (n + 1)-й траектории можно рассматривать любую из совокупности $\{ {{\Omega }_{i}}\} $, i = 1, 2, …, n + 1. Соответствующие различным сочетаниям произведения и суммы наиболее просто вычислять последовательным пересчетом соответственно возрастанию s. Получаемую таким образом статистическую оценку обозначим через $\tilde {\lambda }_{n}^{{(i)}}$. Пусть $\tilde {\lambda }_{n}^{*} = \frac{1}{{n + 1}}\sum\limits_{i = 1}^{n + 1} \tilde {\lambda }_{n}^{{(i)}}$.

Теорема 2. Если выполняются условия леммы 1, причем $\rho ({{K}_{p}}(\sigma )) < 1 - \varepsilon $ (ε > 0) $\forall \sigma $, то ${\text{E}}\tilde {\lambda }_{n}^{*} = {{\lambda }_{n}}$ и ${\text{D}}\tilde {\lambda }_{n}^{*} \leqslant {\text{D}}{{\tilde {\lambda }}_{n}} < + \infty $.

Доказательство cтроится путем прямого осреднения и оценок на основе сформулированных выше утверждений и условной независимости траекторий ${{\Omega }_{1}},\; \ldots ,\;{{\Omega }_{{n + 1}}}$.

Оценку $\tilde {\lambda }_{n}^{*}$ естественно называть комбинированной. Отметим, что величина ${\text{D}}\tilde {\lambda }_{n}^{*}$ может быть существенно меньше величины ${\text{D}}{{\tilde {\lambda }}_{n}}$ при ${\text{D}}\sigma \ll 1$ вследствие слабой коррелированности реализаций произведений в $\tilde {\lambda }_{n}^{*}$ для фиксированных значений s.

Аналогичная несмещенная элементарная оценка на основе леммы 1 для ${{\mu }_{n}} \approx {\text{E}}{{\left[ {\frac{{J{\text{'}}(\sigma )}}{{J(\sigma )}}} \right]}^{2}}$ строится по формуле

с использованием n + 2 условно независимых траекторий. Более эффективная комбинированная оценка $\tilde {\mu }_{n}^{*} = \frac{{\sum\limits_{i = 1}^{n + 1} \mu _{n}^{{(i)}}}}{{n + 1}}$ строится по аналогии с $\tilde {\lambda }_{n}^{*}$ для фиксированного значения $\xi {\text{'}}({{\Omega }_{{n + 2}}},\sigma )$.

Для численного тестирования был рассмотрен односкоростной процесс переноса частиц с изотропным рассеянием (в том числе и после акта деления) в шаре радиуса R = 7.72043 с параметрами: σ – случайная величина, равномерно распределенная в интервале (0.7; 1.3);

Как указано в [2], такой шар при σ ≡ 1 с достаточной степенью точности критичен. Соответствующая прецизионная оценка величины λ методом Монте-Карло подтверждается расчетами в улучшенном диффузионном приближении [6], которые дали значение λ_d (σ = 1) = 0.000000… С использованием этого приближения путем моделирования значений σ были получены следующие тестовые статистические оценки: ${\text{E}}{{\lambda }_{d}}$ = –0.00103, ${\text{D}}{{\lambda }_{d}} = 0.00022$. Среднеквадратическая погрешность имеет здесь порядок последнего десятичного разряда. Для построения эффективного алгоритма метода Монте-Карло на основе теоремы 2 в сформулированную модель было введено поглощение с постоянным неслучайным коэффициентом σ_c/${v}$, который приводит к замене $\lambda \mapsto \lambda - \frac{{{{\sigma }_{c}}}}{{v}} $ $\forall \sigma $. Отметим, что такой прием является универсальным и может существенно повысить эффективность весового метода, исключая необходимость ветвления моделируемых траекторий, как это видно из дальнейшего.

Используя уравнение переноса (см., например, [1, 2]), можно сделать замену ${{\sigma }_{f}} \mapsto {{\sigma }_{f}} + {{\sigma }_{c}}$, $\nu \mapsto \frac{{\nu {{\sigma }_{f}}}}{{({{\sigma }_{f}} + {{\sigma }_{c}})}}$. В численном эксперименте моделировался процесс переноса с константами: $\sigma _{s}^{*} = {{\sigma }_{s}}$, $\sigma _{f}^{*} = {{\sigma }_{f}} + {{\sigma }_{c}}$, $\nu * = 1$. Вспомогательные веса при этом определяются формулой Q_n = = ${{\left( {\frac{{\nu {{\sigma }_{f}}}}{{{{\sigma }_{f}} + {{\sigma }_{c}}}}} \right)}^{{{{n}_{1}}}}}$, где n₁ – число делений, предшествующих данному столкновению [1]. Было использовано значение σ_c = 0.059, для которого  $\frac{{\nu {{\sigma }_{f}}}}{{{{\sigma }_{f}} + {{\sigma }_{c}}}}$ < 1 и тем самым [1] ρ(K_p) < 1 $\forall \sigma $. Плотность распределения первых столкновений взята в виде

где f_d(r) = $\frac{{C\sin (\unicode{230} r)}}{r}$ – диффузионное приближение к пространственной характеристической функции для σ = 1, $\unicode{230} $ = 0.3739866 [6]. Полагали также h(r, v) ≡ $\frac{{\sin (\unicode{230} (\rho )r)}}{r}$, где (см. [6])

Расчеты показали, что использование таких функциональных параметров алгоритма существенно улучшает сходимость в (4) сравнительно с ${{f}_{d}}({\mathbf{r}}) \equiv {{\left( {\frac{{4\pi {{R}^{3}}}}{3}} \right)}^{{ - 1}}}$ и даже сравнительно с вариантом, в котором f(r, t) определяется формулой (6), а h(r, v) ≡ 1.

Нетрудно проверить, что для сформулированных выше характеристик вычислительной модели выполняются условия теорем 1, 2, и она тем самым позволяет построить оценку значений Eλ и Dλ. Для построения статистических оценок этих величин было реализовано 16 × 10⁹ случайных значений оптической плотности σ; моделирование траекторий осуществлялось параллельно на восьми процессорах Intel Core i7-3930K. При этом полагали n = 5, т.е. для оценки Eλ использовали шесть, а для оценки Eλ² – семь условно-независимых траекторий частиц. Результативные оценки определялись на основе анализа установления (с учетом вероятностной погрешности) двух старших значащих цифр при возрастании параметров t и n. Таким образом получены следующие результативные оценки: Eλ = –0.00104 ± 0.00001 и Dλ = = 0.00020 ± 0.00007 для t = 10, n = 4. В качестве погрешности здесь приведено среднеквадратическое отклонение. Эти оценки вполне соответствуют диффузионным, достаточная точность которых проверена лишь для σ = 1.

Отметим, что полученные результаты переносятся на случай многоскоростного процесса переноса с анизотропным рассеянием соответственно [2].