\documentclass[12pt,reqno]{report}
\usepackage[russian]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{graphicx}
\usepackage{graphicx,color}

\definecolor{red}{rgb}{0.8, .0, .0}
\definecolor{blue}{rgb}{.0, 0.0, 0.9}
\def\blue{\color{blue}}
\def\red{\color{red}}
\textheight 24cm \textwidth 16.5cm \topmargin -2cm \oddsidemargin
0.5cm

\def\sgn{\mathrm{sgn}}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}
\newtheorem{problem}{Задача}
\newtheorem{remark}{Замечание}
\renewcommand{\abstractname}{*}

\renewcommand{\thesection}{\arabic{section}.}
\renewcommand{\thedefinition}{\arabic{definition}.}
\renewcommand{\thelemma}{\arabic{lemma}.}
\renewcommand{\thetheorem}{\arabic{theorem}.}
\renewcommand{\thecorollary}{\arabic{corollary}.}
\renewcommand{\theremark}{\arabic{remark}.}
\renewcommand{\theequation}{\arabic{equation}}


\usepackage{float}
\begin{document}
{\bf УДК} 517.958 \vskip0.2cm


\begin{center}
{\bf \large  Численное решение начально-краевой задачи  для
уравнения переноса излучения  с  френелевскими условиями
сопряжения на границе раздела сред }\footnote[1]{Работа выполнена
при финансовой поддержке Российского научного фонда (проект №
14-11-00079).}
\end{center}

\section{Численные методы}
Для решения операторного уравнения (31) справедливо представление в виде ряда Неймана (\ref{neumann_series}).
\begin{equation}\label{neumann_series}
I = \sum_{i=0}^\infty\ (\mathcal{PB} + \mathcal{ES})^i (\mathcal P h + \mathcal E J)
\end{equation}
Рекуррентная запись этого представления (\ref{neumann_recur}) является базой для построения численных алгоритмов вычисления искомой функции $I$.
\begin{equation}\label{neumann_recur}
\begin{array}{cc}
I_0 = \mathcal P h + \mathcal E J
\\
\\
I_N = (\mathcal{PB} + \mathcal{ES}) I_{N-1} + I_0
\\
\\
I = \displaystyle{\lim_{N \to \infty}I_N}
\end{array}
\end{equation}
Представление искомой функции в таком виде имеет наглядное физическое представление. Каждый член последовательности $I_N$ является аппроксимацией функции $I$ с учётом не более, чем $n$ кратного рассеяния. Применение оператора $\mathcal{PB} + \mathcal{ES}$ к некоторой функции плотности потока позволяет оценить плотность этого потока после его однократного взаимодействия со средой. Следовательно, применение её к предыдущему элементу последовательности $I_{N-1}$ даст нам вклад излучения, испытавшего от $1$ до $n$ актов взаимодействия со средой, слагаемое $I_0$ учтёт вклад излучения, пришедшего напрямую от источников.

Будем использовать схему Монте-Карло для вычисления результатов многократного применения интегрального оператора. Рассмотрим три модификации этого метода.
\subsection{Сравнение трёх модификаций метода Монте-Карло}
Введём обозначения и для самих граничных точек: пространственную координату обозначим за $z_r(r,\omega,t) = r - d(r,-\omega, t)\omega$, а временную $z_t(r,\omega, t) = t - d(r,-\omega, t) / v_i$
\begin{itemize}
\item{С ветвлением --} для вычисления значения $I_N$ в некоторой точке фазового пространства $x$ будем вычислять значения $I_{N-1}$ в трёх других точках. Пространственная и временная координата двух точек будет совпадать и являться ближайшей граничной точкой в направлении $-\omega$, в которой излучение, претерпев отражение или преломление, поменяло своё направление на $\omega$. Третья точка будет выбираться случайно в направлении $-\omega$ на расстоянии $\xi$, с плотностью распределения $\displaystyle{\frac{\exp(-\mu_i \xi)}{(1-\exp(-\mu_i d(r, -\omega, t)))}} $ и с направлением $\zeta$, определяемым индикатрисой рассеяния. В этой точке будет вычислен вклад рассеянного излучения и излучения, дошедшего напрямую от внутренних источников.
\begin{multline}
I_N(r,\omega,t) = \exp(-\mu d(r,-\omega,t)) (h(z_r(r,\omega,t), \omega, z_t(r,\omega,t)) +
\mathcal{B}I_{n-1}^+ (z_r(r,\omega,t), \omega, z_t(r,\omega,t))) + \\ +
\frac{1 - \exp(-\mu d(r,-\omega,t))}{\mu_i} (J(r - \xi\omega, \omega, t - \xi /c_i) + \sigma_i I_{N-1}(r - \xi\omega, \zeta, t - \xi / c_i))
\end{multline}

Отметим что при данном подходе вклад отражённого и преломлённого излучения на каждом шагу вычисляется точно, за исключением отбрасывания вклада излучения, которое провзаимодействовало со средой более чем $N$ раз.
\item{Без ветвления --} при вычислении значения $I_N$ в точке $x$ будем вычислять только одно значение соответствующее либо отражённому, либо преломлённому, либо рассеянному излучению, с вероятностями $p_{tr} = \exp(-\mu d(r,-\omega,t)) T(z_r(r,\omega,t), \omega)$, $p_{re} = \exp(-\mu d(r,-\omega,t)) R(z_r(r,\omega,t), \omega)$ и $p_{sc} = 1 - \exp(-\mu d(r,-\omega,t))$, соответственно.
\begin{equation}
I_N(r,\omega,t) =
\left\{
\begin{array}{cc}
    h(z_r(r,\omega,t), \omega, z_t(r,\omega,t)) + I^+_{N-1} (z_r(r,\omega,t), \omega_{tr}, z_t(r, \omega, t)), &  p_{tr}\\ \\

    h(z_r(r,\omega,t), \omega, z_t(r,\omega,t)) + I^+_{N-1} (z_r(r,\omega,t), \omega_{re}, z_t(r, \omega, t)), & p_{re}\\ \\

    \displaystyle{\frac{1}{\mu_i} J(r - \xi\omega, \omega, t - \xi /v_i)
    + \frac{\sigma_i}{\mu_i}I_{N-1}(r - \xi\omega, \zeta, t - \xi / v_i))}, & p_{sc}
\end{array}
\right.
\end{equation}

\item{С частичным ветвлением --} данный способ является смесью первых двух и заключается в случайном выборе между вычислением вклада излучения с границы и вклада рассеянного излучения. В случае выбора излучения с границы, вычисляется вклад и отражённого и преломлённого излучения, что позволяет вычислить его точно, за исключением ошибок, которые будут допущены на следующих шагах и при отбрасывании остаточного члена ряда Неймана.
\begin{equation}
I_N(r,\omega,t) =
\left\{
\begin{array}{cc}
    h(z_r(r,\omega,t), \omega, z_t(r,\omega,t)) +
\mathcal{B}I_{n-1}^+ (z_r(r,\omega,t), \omega, z_t(r,\omega,t)), & p_{re} + p_{tr} \\  \\

    \displaystyle{\frac{1}{\mu_i} J(r - \xi\omega, \omega, t - \xi /c_i)
    + \frac{\sigma_i}{\mu_i} I_{N-1}(r - \xi\omega, \zeta, t - \xi / c_i))}, & p_{sc}
\end{array}
\right.
\end{equation}

\end{itemize}
Для определения границ эффективности каждого метода, проведём ряд вычислительных экспериментов, варьируя различные параметры среды. В каждом тесте граничные условия будут задаваться на сферической поверхности радиуса $r_0 = 1$ с центром в $c_0 = (0, 0, 0)$ по следующему правилу $h(z, \omega, t) = \exp(2 - |z - (0, 0, 1)| - |\omega - (0, 0, -1)|)$. Измерять плотность потока излучения будем в $(0, 0, 0.7)$ в направлении $(0, 0, 1)$. Внутрь среды будут помещены два вложенных друг в друга сферических включения с радиусами $r_1 = 0.2$ и $r_2 = 0.05$ и центрами в $c_1 = (0, 0.02, 0)$ и $c_2 = (0, 0.01, 0)$.

Коэффициенты рассеяния для основной среды и большого включения равны $\mu_0 = \mu_1 = 1.5$ во втором и третьем тесте, и $\mu_0 = \mu_1 = 0.5 $ --- в первом, в малом включении $\mu_3 = 7.5$. Кроме того во всей среде будет присутствовать 10\% поглощение, т.е. $\sigma_i = 0.9 \mu_i$. Коэффициенты преломления в первом и втором  тесте равны $k_1 = 1, k_2 = 3, k_3 = 10$, а в третьем --- $k_2 = 1.4$.

 Для вычисления искомой величины будем также усреднять значения случайной величины $I_N$ по выборке некоторого фиксированного размера $M$.
\begin{equation}
    \overline{I_N} = \frac{1}{M}\sum_{i = 0}^M I_{N,i}
\end{equation}
 При таком подходе адекватной оценкой эффективности алгоритма является его трудоёмкость, равная произведению среднего времени вычисления одной случайной величины на её дисперсию  $C_{I_N} = D(I_N) T(I_N)$.

 \begin{table}[H]

 \begin{tabular}{l||c|c|c}
 Метод & Среда а) & Среда б) & Среда в) \\
 \hline \hline
С ветвлением &
$ 0.093 \cdot 597.2 = 55.82 $ &
$ 0.208 \cdot 579.0 = 120.6 $ &
$ 0.227 \cdot 624.0 = 141.6 $
 \\
С частичным &
$ 1.574 \cdot 82.80 = 130.3 $ &
$ 0.861 \cdot 60.10 = 51.76 $ &
$ 1.334 \cdot 60.80 = 81.14 $
\\
Без ветвления &
$ 8.463 \cdot 25.10 = 212.4 $  &
$ 2.380 \cdot   25.27 =  60.18 $ &
$ 2.510 \cdot   26.00 =  65.27 $
\end{tabular}

\caption{Результаты сравнительного тестирования ветвящегося, неветвящегося и гибридного методов.}
\label{test_res}
\end{table}

Как видно из результатов тестирования, представленных в табл. \ref{test_res} применение ветвящегося метода наиболее эффективно для слабо рассеивающих сред. Выбор же между гибридным и неветвящимся методом следует делать в зависимости от разности долей отражённого и преломлённого излучения на внутренних границах. Так, например, во втором тесте доля отражённого излучения при прохождении границы большего включения под прямым углом составляет всего 7\%, а в третьем тесте она уже равна 27\%. Таким образом, гибридный метод даёт выигрыш при более равном разделении излучения на отражённое и преломленное.
\subsection{Визуализация среды}
Проведём тестирование численного алгоритма на предмет получения визуализации среды. Сравнение полученных результатов с естественными представлениями о распространении излучения в рассеивающих средах поможет проверить корректность работы метода. Чтобы получить изображение некоторой среды $G$ в изометрической проекции, будем вычислять значения функции плотности потока излучения в узлах прямоугольной сетки на плоскости $\Pi \cap G \ne \emptyset$, в направлении $\omega \perp \Pi$. В работе (ссылка на предыдущую работу) была показана стабилизация решения нестационарного уравнения к стационарному, при фиксированных граничных условиях. Рассмотрим случай, в котором граничные условия также являются нестационарными, т.е. зависят от временной координаты $t$. Изображения, построенные для значений из некоторой возрастающей последовательности $t = t_0, t_1, \hdots, t_n$, будут описывать процесс последовательного распространения излучения в среде.

Как и в предыдущем тесте, границей исследуемой среды будет единичная сфера. Входной поток задаётся выражением: $h(z, \omega, t) = 0.5 + \exp(4 - 3 |z - (0, 1, 0)|^2 - 3 |\omega - (0, -1, 0)|^2 - 10 (t - 10) ^ 2)$. Таким образом источник излучения будет достигать максимума в точке $(0, 1, 0)$, в направлении $(0, -1, 0)$, в момент времени $t = 10$ и экспоненциально убывать при удалении от этой фазовой точки. Экраном --- областью, в точках которой строится изображение, является квадрат с центром в (0.3, 0.1, 0.7), развёрнутый к началу координат.

Сама сцена разделяется на подобласти тремя плоскостями, параллельными координатным осям. Горизонтальная плоскость проходит через начало координат, среда ниже этой плоскости имеет следующие характеристики: $\mu_1 = 8, \sigma_1 = 7, k_1 = 1.1$. Две горизонтальные плоскости помещены на расстояние $0.2$ от начала координат, характеристики среды за плоскостью слева от экрана: $\mu_2 = 10, \sigma_2 = 9, k_2 = 10$ --- вдали от экрана: $\mu_3 = 2, \sigma_3 = 1, k_3 = 2$. Оставшаяся часть заполнена веществом, слабо взаимодействующим с излучением и имеющим единичный коэффициент преломления $\mu_0 = \sigma_0 = 0.0001, k_0 = 1$.

На горизонтальной плоскости находятся 3 шара. Шар в углу сильно рассеивающий, коэффициенты среды в нём в точности повторяют коэффициенты среды ниже горизонтальной плоскости. Справа от него находится зеркальный шар, среда в нём аналогична среде за левой плоскостью. В центре, ближе к экрану помещён прозрачный шар, коэффициенты среды внутри него равны: $\mu_4 = 1, \sigma_4 = 0.5, k_4 = 1.05$. На \ref{visual} показана визуализация среды для 6 моментов времени начиная с $t = 10.8$, с фиксированным шагом $\Delta t = 0.3$.

\begin{figure}[H]
  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0000.png}} t = 10.8 \\
  \end{minipage}
  \hfill
  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0001.png}} t = 11.1 \\
  \end{minipage}
  \hfill
  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0002.png}} t = 11.4 \\
  \end{minipage}

\medskip

  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0003.png}} t = 11.7 \\
  \end{minipage}
  \hfill
  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0004.png}} t= 12.0 \\
  \end{minipage}
  \hfill
  \begin{minipage}[h]{0.328\linewidth}
    \center{\includegraphics[width=1\linewidth]{0005.png}} t = 12.3 \\
  \end{minipage}
  \caption{Визуализация среды}
  \label{visual}
\end{figure}

\end{document}