Журнал вычислительной математики и математической физики, 2022, T. 62, № 2, стр. 249-254

ВВЕДЕНИЕ

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

При построении статистического алгоритма решение краевой задачи записывается в виде математического ожидания некоторой случайной величины $\xi $, т.е. случайная величина $\xi $ является несмещенной оценкой решения задачи. Обычно несмещенные оценки для решения краевых задач строятся на траекториях марковского случайного процесса с дискретным временем (случайного блуждания). Само случайное блуждание происходит в замкнутой ограниченной области $D \subset {{\mathbb{R}}^{m}}$, в которой решается краевая задача. Вероятность перехода $P(x,A)$ из точки $x \in D$ в множество $A \subset D$ является субстохастическим ядром. Исходная краевая задача при этом заменяется эквивалентным ей интегральным уравнением в пространстве непрерывных функций $C(D)$:

где правая часть $F(x)$ уравнения (0.1) определяется через правую часть дифференциального уравнения и граничные условия. Пусть $x(0) = x,x(1),x(2), \ldots $ – цепь Маркова с переходными вероятностями $P(x,B)$. Вообще говоря, это обрывающаяся цепь. На траекториях цепи стандартным образом строится последовательность $\xi (n)$, $n = 0,1, \ldots $, несмещенных оценок для $u(x)$. Пусть $N$ – момент обрыва цепи ($N = \infty $, если цепь не обрывается), $N \wedge n$ – минимум из $N$ и $n$, $\eta (i)$ – несмещенная оценка для $F(x(i))$, тогда

Из марковского свойства цепи и уравнения (0.1) следует

Лемма 1. Пусть $F(x) \geqslant 0$, $u(x)$ – ограниченное решение уравнения (0.1) и оценки $\eta (i)$ неотрицательны. Тогда последовательность несмещенных оценок (0.2) образует мартингал. Если он равномерно интегрируемый, то он сходится с вероятностью 1 к случайной величине $\xi (\infty )$, имеющей конечное математическое ожидание. Для любого марковского момента $\tau $ случайная величина $\xi (\tau )$ является несмещенной оценкой $u(x).$

Поведение траекторий процесса характеризует

Лемма 2. Пусть координатные функции ${{x}_{i}}$, $i = 1,2, \ldots m$, являются супергармоническими для ядра $P(x,A)$, тогда цепь $x(n)$ сходится с вероятностью $1$ к некоторому случайному вектору $x(\infty )$.

Доказательство лемм можно найти в [3]. При этом мы пользуемся терминологией и результами из книги П. Мейера [4]. Вектор $x(\infty )$, как правило, распределен на границе области $D$. В качестве марковского момента $\tau $ обычно выбирается момент первого попадания траектории блуждания в $\varepsilon $-окрестность границы области $D$.

Интегральное уравнение (0.1) обычно получают по формуле Грина, используя функцию Грина или функцию Леви для семейства областей $T(x) \subset D,\quad x \in D$. При этом носителем меры $P(x,A)$ будет либо граница области, либо сама область $T(x).$ Вид носителя меры $P(x,A)$ определяет название блуждания и название статистического метода решения задачи: блуждание по сферам, блуждание по шарам, блуждание по эллипсоидам, блуждание по полусферам, блуждание по границе.

Таким образом, в статье рассматриваются стохастические алгоритмы вычисления значения решения краевой задачи $u(x)$ в точке $x$ ограниченной области по известным значениям этой функции в граничных точках.

1. ПОСТАНОВКА ЗАДАЧИ И ТЕОРЕМЫ О СРЕДНЕМ ЗНАЧЕНИИ

Пусть $D$ – ограниченная область в ${{\mathbb{R}}^{m}}$, $m \geqslant 3$, с границей $\Gamma $, которая состоит из двух подобластей ${{D}_{1}}$, ${{D}_{2}}$ и их общей кусочно-гладкой границы $\gamma $. Для функций $u \in {{C}^{2}}({{D}_{1}} \cup {{D}_{2}}) \cap C(\bar {D})$ определим оператор $M$ равенством

где матрица $\left\| {{{a}_{{i,k}}}(x)} \right\|_{{i,k = 1}}^{m}$ – симметричная, а ее элементы ${{a}_{{i,k}}} \in {{C}^{1}}({{D}_{1}} \cup {{D}_{2}}),$ причем их производные первого порядка непрерывны вплоть до границы в каждой подобласти. Предположим также, что матрица коэффициентов – положительно-определенная, т.е. для некоторой постоянной $\nu > 0$ при всех $x \in D$ и всех $z \in {{\mathbb{R}}^{m}}$ справедливо неравенство Таким образом, $M$ – эллиптический оператор.

Рассмотрим краевую задачу

где $n(x) = ({{n}_{1}}(x),...,{{n}_{m}}(x))$ – внешняя нормаль к границе области ${{D}_{1}}$. На внешней границе области $D$ функция $u(x)$ известна и равна $\varphi (x)$.

Практически важным примером поставленной задачи является задача для уравнения

с условиями согласования нормальных производных

Функции ${{\varepsilon }_{1}}(y)$ и ${{\varepsilon }_{2}}(y)$ предполагаются ограниченными снизу постоянной $\nu > 0$ и непрерывно дифференцируемыми в ${{\bar {D}}_{1}}$ и ${{\bar {D}}_{2}}$ соответственно.

В случае, когда функции ${{\varepsilon }_{1}}(y)$ и ${{\varepsilon }_{2}}(y)$ являются различными положительными постоянными ${{\varepsilon }_{1}}$ и ${{\varepsilon }_{2}}$ и граница раздела сред плоская, несмещенные оценки решения задачи (1.4), (1.5) построены в [3]. Аналогичные, но смещенные оценки, в областях с произвольной границей построены в работах Н.А. Симонова (см. [5]$,$ [6])$.$

Для построения несмещенных и малосмещенных оценок поставленной краевой задачи (1.4), (1.5) с переменными коэффициентами нам потребуется представление решения $u(x)$ для точек $x \in \gamma $ в интегральной форме. Для этого мы воспользуемся формулами Грина для подобласти $T$ с кусочно-гладкой границей, расположенной в ${{\bar {D}}_{1}}$ или ${{\bar {D}}_{2}}$:

справедливой для $u,{v} \in {{C}^{2}}(T) \cap C(\bar {T}).$ Нормаль в этой формуле внешняя по отношению к области.

Пусть граница $\gamma $ гладкая в окрестности $x \in \gamma .$ Рассмотрим шар $T(x) \subset D$ с центром в точке $x$ и радиуса $R = R(x).$ Рассмотрим два его полушара ${{T}_{1}}(x) = T(x) \cap {{D}_{1}}$ и ${{T}_{2}}(x) = T(x) \cap {{D}_{2}}$ и шар ${{T}_{\delta }}$ с центром в точке $x$ и радиуса $\delta .$ В качестве функции ${v}(y)$ рассмотрим функцию Леви

где $r = {\text{|}}y - x{\kern 1pt} {\text{|}}$ – расстояние между точками $x$ и $y.$ Применяя формулу Грина к области $T = {{T}_{1}}{{\backslash }}{{T}_{\delta }}$, при $\delta \to 0$ получаем вторую формулу Грина в которой ${{\sigma }_{m}}$ обозначает “площадь поверхности” единичной сферы в ${{\mathbb{R}}^{m}}.$ Аналогично получаем вторую формулу Грина для области ${{T}_{2}}$

Заметим, что

равен нулю на сферической части границы, сама функция ${v}(y)$ – тоже. Пусть $\gamma (x) = T(x) \cap \gamma $ – часть границы $\gamma ,$ лежащая в шаре. Упрощая формулы, получаем равенства Нормали в этих формулах являются противоположными векторами. При сложении этих формул, используя условие на нормальные производные решения, получаем

Для плоской границы $\gamma $ нормальная производная в этих формулах равна нулю. После несложных вычислений находим значение оператора

и получаем формулу среднего значения

Отметим, что осредняющее ядро в формуле (1.8) неотрицательно при выполнении условия

которое справедливо при достаточно малых $R$ в силу сделанных предположений о параметрах краевой задачи. Кроме того, оно стохастическое, так как после подстановки в формулу (1.8) решения краевой задачи $u(x) \equiv 1$ получим, что интеграл от ядра равен единице.

Для точек, не лежащих на границе $\Gamma \cup \gamma $, верны формулы среднего значения в шаре, а именно, справедлива следующая лемма.

Лемма 3. Если точка $x \in {{D}_{1}}$ и $T(x) \subset {{D}_{1}},$ то

Аналогично, если $x \in {{D}_{2}}$ и $T(x) \subset {{D}_{2}}$, то

Доказательство леммы непосредственно следует из формулы (1.6) и равенства нулю на границе шара как самой функции ${v}(y),$ так и ее градиента.

2. СТОХАСТИЧЕСКИЙ АЛГОРИТМ В СЛУЧАЕ ПЛОСКОЙ ГРАНИЦЫ РАЗДЕЛА ОБЛАСТИ

В случае переменных ${{\varepsilon }_{1}}$ и ${{\varepsilon }_{2}}$ для построения малосмещенной оценки решения краевой задачи предлагается следующий алгоритм.

1. Задаем призвольное $\delta > 0,$ которое будет определять порядок погрешности приближенного решения задачи.

2. Определяем две $\delta $-окрестности границы ${{\gamma }_{{1,\delta }}}$ и ${{\gamma }_{{2,\delta }}}$ как множества точек из ${{D}_{1}}$ и ${{D}_{2}}$ соответственно, отстоящих от границы $\gamma $ не более, чем на $\delta .$ Также определяем ${{\Gamma }_{\delta }}$ – окрестность внешней границы.

3. Для каждой точки $x$ из произвольной $\delta $-окрестности границы определяем точку $x{\kern 1pt} *$ как ближайшую к $x$ точку границы.

4. Фиксируем начальную точку $x$, для которой будет строиться статистическая оценка $\xi (x)$ для значения $u(x)$ и полагаем ${{x}_{0}}: = x$.

5. Если ${{x}_{0}} \in {{\Gamma }_{\delta }},$ то полагаем $\xi (x): = \varphi (x{\kern 1pt} *)$ и останавливаемся.

6. Если ${{x}_{0}} \in {{\gamma }_{{1,\delta }}}$ или ${{x}_{0}} \in {{\gamma }_{{2,\delta }}},$ то ${{x}_{0}}: = x_{0}^{ * }$ и переходим из точки ${{x}_{0}}$ в точку ${{x}_{1}},$ распределенную внутри области $T({{x}_{0}})$ с плотностью, которая является ядром в формуле (1.8). Полагаем ${{x}_{0}}: = {{x}_{1}},$ переходим к пункту 5.

7. Переходим из точки ${{x}_{0}}$ в точку ${{x}_{1}},$ распределенную внутри области $T({{x}_{0}})$ с плотностью, которая является ядром в формуле (1.9) или формуле (1.10). Полагаем ${{x}_{0}}: = {{x}_{1}},$ переходим к пункту 5.

Назовем построенное случайное блуждание блужданием по шарам. Все переходы в этом блуждании легко моделируются методом отбора. Процедура моделирования подробно описана в [2].

Работоспособность алгоритма обеспечивается следующими его свойствами.

– Блуждание по шарам в области ${{D}_{1}}$ (в области ${{D}_{2}}$) достигает ее $\delta $ границы за конечное число шагов. При этом среднее число шагов конечно. (Эти свойства доказаны для любого блуждания по эллипсоидам, в том числе и для блуждания по шарам в [2].)

– Среднее число посещений границы блужданием по шарам конечно. Этот факт является очевидным следствием строго марковского свойства блуждания и ограниченности снизу вероятности события, состоящего в том, что процесс попадет впервые в ${{\Gamma }_{\delta }}$ раньше, чем в ${{\gamma }_{\delta }}.$ (Асимптотика этих средних при $\delta \to 0$ конечно важна, но ее удается определить крайне редко.)

Замечание 1. Как известно (см. [3]), в случае постоянных коэффициентов ${{\varepsilon }_{1}}$ и ${{\varepsilon }_{2}}$ для обеспечения выхода блуждания на плоскую границу раздела $\gamma $ можно использовать обычную функцию Грина для оператора Лапласа в полусфере. В алгоритме блуждания по полусферам существенно меньше смещение статистической оценки, поскольку оно определяется лишь заменой точного значения функции в последней точке блуждания на значение в ближайшей к ней точке границы $\Gamma $.

3. СТОХАСТИЧЕСКИЙ АЛГОРИТМ В СЛУЧАЕ ГЛАДКОЙ ГРАНИЦЫ РАЗДЕЛА ОБЛАСТИ

В случае гладкой границы $\gamma $ вместо формулы (1.8) приходится использовать формулу (1.7), которая не имеет вероятностного смысла в связи с присутствием интеграла

Чтобы его убрать, нужно его оценить сверху величиной которая должна быть меньше $\delta \tfrac{{{{\sigma }_{m}}}}{2}R(m - 2)({{\varepsilon }_{1}}(x) + {{\varepsilon }_{2}}(x)).$ Фактически получается неравенство на величину телесного угла, под которым граница $\gamma (x)$ видна из точки $x.$ Выбор радиуса шара теперь зависит от кривизны поверхности в окрестности точки $x.$ После выбрасывания интеграла получаем формулу применение которой позволяет построить статистическую оценку решения $u(x).$

В случае постоянных ${{\varepsilon }_{1}}$ и ${{\varepsilon }_{2}}$ можно использовать функцию Леви

Тогда получится процесс блуждания по сферам.

Замечание 2. Используя формулу (1.6), можно совсем легко построить и обосновать для уравнения (1.4) алгоритм решения задачи со смешанными граничными условиями в областях с кусочно-гладкой границей, на одной части которой задано условие Неймана, а на другой – уcловие Дирихле. Такие алгоритмы рассматривались для уравнений Пуассона и Лапласа в [3], [5]–[7].

4. ЗАДАЧА ДЛЯ ПРОИЗВОЛЬНОГО ОПЕРАТОРА

Вернемся теперь к задаче для произвольного самосопряженного эллиптического оператора (1.1). Введем необходимые обозначения. Матрицу коэффициентов оператора будем обозначать через $A(x),$ а ее определитель – ${\text{|}}A(x){\kern 1pt} {\text{|}}$. Рассмотрим квадратичную форму ${{\sigma }^{2}}(y,x) = (y - x){\kern 1pt} {\text{'}}{{A}^{{ - 1}}}(x)(y - x),$ где $(y - x){\kern 1pt} {\text{'}}$ обозначает строку, транспонированную к столбцу $(y - x)$. Для точек $x \in {{D}_{1}} \cup {{D}_{2}}$ определим функцию $\mathcal{L}(y,x)$ равенством

Рассматривая в качестве множеств $T(x)$ области $T(x) = \{ y \in D:{{\sigma }^{2}}(y,x) < {{R}^{2}}\} ,$ можно определить процесс блуждания по эллипсоидам и на его траекториях уже описанным методом построить статистические оценки для решения задачи (1.2), (1.3).

Авторы выражают благодарность рецензенту за замечания, способствовавшие улучшению работы.