Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 507, № 1, стр. 81-85

1. Работа посвящена построению и оптимизации численно-статистических проекционных оценок одномерных функциональных характеристик $\varphi (z),0 \leqslant z \leqslant H < + \infty ,$ решения интегральных уравнений вида

или $\Phi = K\Phi + f.$ Здесь $x \in {{R}_{n}},$ $K \in [{{L}_{1}} \to {{L}_{1}}],$ $f \in {{L}_{1}}(x)$ и спектральный радиус $\rho ({{K}_{1}}) < 1$, где K₁ – интегральный оператор с ядром ${\text{|}}k(x{\kern 1pt} ',x){\text{|}}$. Функция $\varphi (z)$ получается интегрированием решения $\Phi (x)$ по всем переменным, кроме $z \in {{R}_{1}}$ в некоторой системе координат. Для простоты изложения далее предполагается, что z – одна из координат фазовой точки: $z = z(x)$.

Методом Монте-Карло уравнения типа (1) решаются с помощью моделирования специальных цепей Маркова [1]. Важной частью такого метода являются оценки функциональных зависимостей, в основном разные типы ядерных и статистических проекционных оценок (см., например, [1]). Решая конкретную задачу, необходимо проводить сравнительный анализ таких оценок, предварительно их оптимизируя. Следует отметить, что использовать статистические проекционные оценки в методе Монте-Карло впервые предложил Н.Н Ченцов [2], разработавший общую методику их оптимизации, которая требует существенной доработки и детализации в конкретных задачах. Простая приближенная оптимизация специальных статистических проекционных оценок для решения задач из областей эффективного применения метода Монте-Карло (переноса излучения с учетом поляризации, теории разреженных газов) использована в [3, 4]. Она может быть существенно улучшена с использованием излагаемых далее результатов.

В настоящей работе предлагается детализированная универсальная методика оптимизации для практически важного класса статистических проекционных оценок с использованием комбинации аналитических и соответствующих численных расчетов.

С целью построения функции $\varphi (z) \in {{L}_{2}}(0,H)$ используется ее разложение по полиномам Лежандра $\{ {{\psi }_{i}}(z)\} $ (см. [5]), ортонормированным в заданном интервале $(0,H)$:

В процессе расчетов значения полиномов вычисляются рекуррентно [6]:

Использование полиномов Лежандра связано с тем, что реализация ортогональных разложений с адаптированным весом $p \approx \varphi $ может быть арифметически весьма затруднительной, как это показало решение рассматриваемой далее типовой тестовой задачи теории переноса с экспоненциальной асимптотикой $\varphi (z)$ при $H \to \infty $.

Имея достаточно хорошее приближение $\varphi \, \approx \,{{\varphi }_{0}}$, можно использовать оценку $\varphi _{m}^{{(1)}} = {{\varphi }_{0}}{{\left( {\frac{\varphi }{{{{\varphi }_{0}}}}} \right)}_{m}}$, для которой ${\text{||}}\varphi - \varphi _{m}^{{(1)}}{\text{||}} \leqslant \max {{\varphi }_{0}}\left\| {\frac{\varphi }{{{{\varphi }_{0}}}} - {{{\left( {\frac{\varphi }{{{{\varphi }_{0}}}}} \right)}}_{m}}} \right\|$.

Следует отметить, что численно-статистическая оптимизация оценки $\varphi _{m}^{{(1)}}$ требует дополнительного исследования в каждом конкретном случае. В рассматриваемой далее тестовой задаче использование в качестве ${{\varphi }_{0}}$ асимптотики решения при $H \to \infty $ не дало улучшения статистической оценки, сравнительно с ${{\tilde {\varphi }}_{m}}$, вследствие увеличения дисперсий оценок коэффициентов разложения и недостаточного приближения $\varphi (z)$ функцией ${{\varphi }_{0}}(z)$для малых $z$.

2. Рандомизация оценки (2) получается путем вычисления линейных функционалов (φ, ψ_i) = = $\int \varphi (z){{\psi }_{i}}(z)dz$ методом Монте-Карло с использованием так называемой “оценки по столкновениям” (см., например, [1]) ${{\xi }_{i}} = \sum\nolimits_{k = 0}^{{{N}_{c}}} {{Q}_{k}}{{\psi }_{i}}(z({{x}_{k}})).$ Здесь ${{x}_{0}},{{x}_{1}},...,{{x}_{{{{N}_{c}}}}}$ – вспомогательная обрывающаяся с вероятностью 1 цепь Маркова, определяемая переходной плотностью $p(x{\kern 1pt} ',x)$ и плотностью начального состояния f₀(x), Q₀ = $f({{x}_{0}}){\text{/}}{{f}_{0}}({{x}_{0}}),$ ${{Q}_{k}} = {{Q}_{{k - 1}}}k({{x}_{{k - 1}}},{{x}_{k}}){\text{/}}p({{x}_{{k - 1}}},{{x}_{k}})$, причем выполняются условия несмещенности [1]. Предполагается, что $\rho ({{K}_{p}}) < 1$, где K_p – интегральный оператор с ядром ${{k}^{2}}(x{\kern 1pt} ',{{x}^{2}}){\text{/}}p(x{\kern 1pt} ',x)$. При сформулированных выше условиях имеем ${\text{E}}{{\xi }_{i}}\, = \,(\varphi ,{{\psi }_{i}})\, = \,{{a}_{i}},$ ${\text{D}}{{\xi }_{i}} < + \infty $.

Рандомизированная оценка после реализации $N$ траекторий строится следующим образом:

Здесь ${{\tilde {a}}_{i}} = \frac{1}{N}\sum \xi _{i}^{{(k)}}$, где ${{\{ \xi _{i}^{{(k)}}\} }_{{k = 1,...,N}}}$ – выборочные значения “оценки по столкновениям” ${{\xi }_{i}}$. Выполняются соотношения ${\text{E}}{{\tilde {\varphi }}_{m}} = {{\varphi }_{m}},$ ${\text{D}}{{\tilde {\varphi }}_{m}}\, < \, + {\kern 1pt} \infty $.

По аналогии с [2] имеем

где ${{L}_{2}}(m) = {\text{||}}\varphi - {{\varphi }_{m}}{\text{|}}{{{\text{|}}}^{2}}$.

3. Перейдем теперь к решению задачи минимизации квадрата среднеквадратической погрешности L(m). С этой целью сформулируем утверждение, которое является простым следствием теоремы 4.10 из [5].

Лемма 1. Если для $r \geqslant 1$ выполняется соотношение $\frac{{{{d}^{r}}\varphi (z)}}{{{{d}^{r}}z}} < c < + \infty $, при $0 \leqslant z \leqslant H,$ то

Расчеты, проведенные, в частности, при выполнении работ [3, 4] и настоящей работы, показали, что изменением величины ${\text{D}}{{\tilde {a}}_{i}}$, сравнительно с $a_{i}^{2}$, можно пренебречь, т.е. полагать ${{L}_{1}}(m) \approx {{C}_{1}}m = (d{\text{/}}n)m$.

Было также замечено, что последовательности $a_{i}^{2}$ для четных и нечетных номеров в асимптотике могут существенно различаться и поэтому целесообразно с целью минимизации $L(m)$ рассматривать ${{L}_{2}}(m)$ для нечетных m в виде

Это свойство коэффициентов было проверено на примере разложения функции φ(z) = = $\exp ( - z{\text{/}}5.4)$ по полиномам Лежандра в интервале $0 < z < 10$.

Лемма 1 дает основание предполагать, что в случае достаточно гладкой $\varphi (z)$ выполняется соотношение

Достаточно просто доказывается следующее утверждение.

Лемма 2. Если

то выполняются соотношения

Лемма 2 показывает, что при выполнении соотношения (5) асимптотически при $N \to \infty $ вследствие соотношения ${{C}_{1}} = d{\text{/}}N$ имеем

Для решения рассматриваемой задачи минимизации $L(m)$ необходимо с помощью предварительных расчетов оценивать величины ${{C}_{2}}$ и $k$ (${{C}_{1}}$ – это среднее значение слабо меняющейся величины ${\text{D}}{{\tilde {a}}_{i}}$). Использование соотношения (7) для оценки $k$ неэффективно из-за сильного возрастания относительной статистической погрешности оценок ${{\tilde {a}}_{i}}$. Поэтому для исследования асимптотики величины ${{L}_{2}}(m)$ используется тот факт, что соотношение ${{L}_{2}}(m) \asymp {{C}_{2}}{\text{/}}{{m}^{k}}$ соответствует соотношению $b_{i}^{2} \asymp \frac{{C_{2}^{{(1)}}}}{{{{i}^{{k + 1}}}}}$, причем $b_{i}^{2}{\text{/}}b_{{i + 1}}^{2} = ((i + 1){\text{/}}i{{)}^{{k + 1}}}$.

Определив интервал $({{i}_{1}},{{i}_{2}})$ слабого изменения отношения $b_{i}^{2}{\text{/}}b_{{i + 1}}^{2}$, можно оценить подходящее значение $k + 1$ путем осреднения отношений $\ln (b_{i}^{2}{\text{/}}b_{{i + 1}}^{2}){\text{/}}\ln ((i + 1){\text{/}}i)$, т.е. по формуле

Заметим, что значение k в (4) практически ограничивается возможным нарастанием производных; это было замечено и при решении рассматриваемой далее тестовой задачи. Оценив, как указано выше, значение k, далее можно оценить коэффициент ${{C}_{2}}$, используя те же значения $\tilde {b}_{i}^{2}$ по формуле

Замечательным свойством этой формулы является то, что она не требует знания ${{a}_{i}}$ для $i\, > \,2{{i}_{2}}$ + 2.

4. В качестве тестовой рассматривалась задача об оценке плотности $\varphi (z)$ столкновений частицы в слое $0 \leqslant z \leqslant H = 10$ рассеивающего и поглощающего вещества для источника столкновений с плотностью f(z₁, z₂, z; ω) = e^–zδ(z₁ – 0)δ(z₂ – 0)$\delta (\omega - {{\omega }_{0}}),$ $z > 0,$ где ${{\omega }_{0}} = (0,0,1)$ – направление скорости частицы, вызывающей начальное столкновение.

Параметры среды: коэффициент ослабления σ = 1, вероятность рассеяния ${{\sigma }_{s}}{\text{/}}\sigma = 0.9$, вероятность поглощения ${{\sigma }_{c}}{\text{/}}\sigma = 0.1$, индикатриса рассеяния Хеньи–Гринстейна:

Средний косинус угла рассеяния ${{\mu }_{0}} = 0.9$. При этом значение $\mu $ моделируется по формуле:

где ${\text{rand}}$ – случайное число, равномерно распределенное в $(0,1)$.

Отметим, что при $\sigma \equiv 1$ осредненная по ${{z}_{1}},{{z}_{2}}$ плотность столкновений φ(z) = $\int \int \int \Phi ({{z}_{1}},{{z}_{2}},z;\omega )d{{z}_{1}}d{{z}_{2}}d\omega ,$ где $\Phi $ – интенсивность излучения. Таким образом, фактически рассматривается задача, близкая к проблеме Милна [7]. Для данных параметров довольно высокую точность имеет транспортное приближение [7], которое дает для $\varphi (z)$ следующую асимптотическую оценку ${{\varphi }_{{as}}}(z) \asymp {{e}^{{ - \lambda z}}},\lambda \approx \frac{1}{{5.4}}$.

Для построения оптимальной согласно сказанному выше оценки ${{\tilde {\varphi }}_{m}}(z)$ в рамках сформулированной физико-вероятностной модели было реализовано 10⁹ траекторий частиц-квантов излучения. В результате расчетов были получены значения слагаемых ${{k}_{i}} + 1$ в сумме (8), приведенные для ${{i}_{1}} \leqslant i \leqslant {{i}_{2}} = 11$ в табл. 1. Эта таблица показывает практически достаточную устойчивость оценки параметра степенного изменения величины $b_{i}^{2}$ в интервале $3 \leqslant i \leqslant 11$. Из формулы (8) получаем $\tilde {k} = 4.97 \approx 5$.

По формуле (9) здесь получается значение ${{C}_{2}} \approx 0.094$, а осреднение величин ${\text{D}}{{\tilde {a}}_{i}}$ в интервале $0 \leqslant i \leqslant 30$ дает ${{C}_{1}} \approx {{10}^{{ - 9}}}$. С использоваеним этих оценок по формулам (6), (7) получено: ${{m}_{{opt}}} \approx 28,{\text{E}}{{(\varphi - {{\tilde {\varphi }}_{{28}}})}^{2}} \approx \tilde {L}(28)$ = 3.3 × 10^–8, т.е. среднеквадратическая погрешность оценивается величиной ${{\tilde {\delta }}_{2}} = \sqrt {\tilde {L}(28)} \approx 1.8$ × 10^–4.

Известно [1], что в данной задаче $\varphi (z)$ равно среднему числу пересечений частицей уровня $z$ с весом $1{\text{/|}}{{\omega }_{z}}{\text{|}}$. При ограничении ${\text{|}}{{\omega }_{z}}{\text{|}} > \varepsilon $ дисперсия соответствующей статистической “локальной оценки” ${{\tilde {\varphi }}_{l}}(z)$ конечна, а относительное смещение не превосходит $\varepsilon $. В расчетах было использовано $\varepsilon = {{10}^{{ - 6}}}$, а относительная статистическая погрешность при $N = {{10}^{9}}$ имеет порядок величины 0.001%. В результате расчетов было получено значение $\tilde {\delta }_{2}^{{(l)}} = {{\tilde {\delta }}_{2}}({{\tilde {\varphi }}_{{28}}} - {{\tilde {\varphi }}_{l}}) \approx 2.2$ × 10^–4, а минимальное значение $\delta _{2}^{{(l)}} = 2.0$ × 10^–4 реализуется при m = 36.

Полученное отличие $\tilde {\delta }_{2}^{{(l)}}$ от ${{\tilde {\delta }}_{2}}$ можно объяснить асимптотическим характером леммы 2, а также статистической погрешностью оценки ${{\tilde {\varphi }}_{l}}$. С помощью прямого дифференцирования интегрального уравнения вида (1) “локальная оценка” была распространена на производные функции $\varphi (z)$. Уже первая производная для больших значений $z$ оказалась на порядок больше $\varphi (z)$; этим объясняется полученное по формуле (8) значение $k = 5$ (хотя $\varphi (z)$, по-видимому, бесконечно дифференцируема при $0 < z < H)$).

Как уже было указано в конце пункта 1 текста статьи, разложение функции $\varphi {\text{/}}{{\varphi }_{0}}$ здесь не улучшает оценку из-за недостаточного приближения функции $\varphi (z)$ экспонентой $\exp ( - z{\text{/}}5.4)$. Попытка построить оценку ${{\tilde {\varphi }}_{m}}(z)$ на основе разложения по полиномам, ортогональным с весом $\exp ( - z{\text{/}}5.4)$, оказалась неудачной из-за необходимости расчетов со слишком большим количеством арифметических разрядов.

Использование разложения по полиномам Лагерра, ортонормированным на $(0, + \infty )$, здесь заведомо неэффективно из-за разрыва функции $\varphi (z)$ в точке $z = H = 10$. Однако соответствующий алгоритм был улучшен путем моделирования траекторий в полупространстве $z \geqslant 0$ с игнорированием вклада в оценку от столкновений на траекториях, до этого покинувших слой $0 \leqslant z \leqslant H$. Таким образом, $\varphi (z)$ была продолжена непрерывно на $z \geqslant H$, но с разрывом 1-го рода производной для z = H. При этом была получена оценка ${{L}_{2}}(m) \approx \frac{{{{C}_{2}}}}{m}$ и соответственно (6) ${{m}_{{opt}}}$ получалось путем уравнивания ${{L}_{2}}(m)$ и ${{L}_{1}}(m)$. Для такого “регуляризованного” алгоритма в случае $N = {{10}^{9}}$ было получено значение ${{\tilde {\delta }}_{2}} \approx 3.6$ × 10^–3, в то время как для алгоритма, связанного с полиномами Лежандра, имеем ${{\tilde {\delta }}_{2}} \approx 1.8$ × 10^–4.

Большой объем проведенных расчетов показал, что при малой величине среднеквадратического отклонения практически всегда малым является и равномерное отклонение рассматриваемой статистической проекционной оценки от искомого решения. По-видимому, это связано с достаточно точной среднеквадратической оценкой производной от решения, что требует дополнительного довольно сложного исследования.