Решение задачи Дирихле методом случайного блуждания на плоскости

В работе рассматриваются способы решения задачи Дирихле на плоскости и система (программа) решения задачи Дирихле на плоскости с помощью случайного блуждания. Программа написана на C++. Приводится математическое объяснение правильности реализации. Рассматривается решение для треугольника, круга двух квадратов разного размера. Система создана для демонстрации студентам математических методов.

Введение

В двухмерном пространстве оно записывается следующим образом:

В трехмерном пространстве уравнение Лапласа выглядит так:

Пусть на плоскости дана некоторая область M с кусочно-гладкой границей G. Внутри области и на ее границе определена функция u(x,y). В области G эта функция имеет двойные частные производные по x и по y. Оператор Лапласа ∆ от этой функции внутри M равен 0.

Обозначим функцию на границе области M через f(x,y). Требуется найти u(x,y) внутри области, если известны значения f(x,y) на границе. Такая задача называется первой краевой задачей или внутренней задачей Дирихле.

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

Для решения задач математической физики часто используется метод конечных разностей (Тихонов, Самарский, 2004; Самарский, 1977; Самарский, 1971; Демидович и др., 2010; Демидович, 1967; Тихонов, 2002; Бахвалов, 2023; Андреев, 2020; Яковлев, 2008).

Он состоит в том, что производные заменяются разностными схемами. Метод конечных разностей является сеточным методом.

В данной статье рассматриваются примеры решения задач математической физики на плоскости не только традиционными методами, а также с помощью траекторий блуждания (Демидович и др., 2010; Демидович 1967).

Решение задачи Дирихле с помощью конечных разностей

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

Обычно при решении с помощью конечных разностей используется метод сеток (сетки также используются в основном разделе, где для решения задачи применяется теория вероятностей). При этом производятся следующие действия (Демидович и др., 2010; Демидович 1967):

• В области M создается сетка Mh из клеток шириной и высотой h.

• Для узлов созданной сетки пишется соответствующее конечно-разностное уравнение.

• С помощью граничных условий задаются значения решения в узлах границы области Mh.

На плоскости для уравнения (1) получается система конечно-разностных уравнений с большим числом неизвестных. Решение этой системы дает значение функции в узлах системы. Для сетки при перемещении по ней (по узлам) можно задать соотношения: xi= x0 + ih, yj= y0 + jh (i,j = 0, ±1, ±2, …)

При этом надо стремиться, чтобы узлы сетки или принадлежали области M, или отстояли от нее меньше чем на h. Чтобы получить конечно-разностное уравнение в уравнении (1) производные заменяют отношениями конечных разностей (Демидович, 1967, с. 257):

Если правые части просуммировать, прировнять 0 и перенести функцию u (x, y) в правую часть, то получим соотношение:

Обозначим ui,j = u (xi, yj). Тогда для каждой внутренней точки (xi, yj) согласно (2) можно написать:

ui,j = ¼ (ui-1,j + ui+1,j + ui,j-1 + ui,j+1 ) (3)

Узлы сетки называются соседними, если находятся на расстоянии h друг от друга.

Если узел расположен внутри области M и все четыре соседних с ним узла также расположены в сетке, то он называется внутренним. Иначе он считается граничным.

Граничным узлом первого рода считается такой узел, у которого есть соседний внутренний узел. В противном случае это граничный узел второго рода.

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

Полученная система (3) является неоднородной линейной системой. Число неизвестных равно числу уравнений. Она имеет единственное решение (Демидович, 1967, с. 263):

Случайное блуждание

Предположим на плоскости имеется ограниченная область M в виде прямоугольника. В точке А этой области есть частица, которая перемещается на растояние h с вероятностью ¼ налево, направо, вверх или вниз. С вероятностью 1 (при большом числе шагов) эта частица достигнет границы M (см. рис. 2). На рисунке изображены 2 траектории (одна красная, а другая - зеленая). Обе начинаются в одной точке А (x0, y0).

Использование случайного блуждания для решения задачи Дирихле

На плоскости дана некоторая область M с кусочно-гладкой границей G. Внутри области и на ее границе определена функция u(x,y). Оператор Лапласа от этой функции внутри M равен 0 (см. формулу (1)). ∆u = 0

Кратко суть решения состоит в следующем. Строится такая же сетка, которая была рассмотрена выше в разделе «Решение задачи Дирихле с помощью конечных разностей». Чтобы найти значение функции u(x,y) в некотором внутреннем узле (x,y) области M, с помощью датчика случайных чисел порождается большое число траекторий, начинающихся в этом узле и оканчивающихся в граничном узле сетки. В результате получается случайная величина из значений функции u(x,y) в конечных узлах сетки. Математическое ожидание для конечного значения траекторий берется в качестве значения для фунции u(x,y).

Ниже доказывается справедливость этого утверждения. Данный алгоритм составлен на основе теории изложенной в (Демидович, 1967).

Алгоритм ссылается на приведенный выше раздел "Решение задачи Дирихле с помощью конечных разностей".

Узлы на границе области назовем граничными. Предполагается, что если при блуждании частица попадает на границу, то частица останавливается и траектория заканчивается в этом узле. Узлы (и точки) расположенные в области M и на ее границе будем обозначать Mij(xi,yj), Mij, Mpq .

Пусть P(i,j;p,q) – вероятность того, что частица из узла Mij попадает в узел Mpq. После попадания частицы на границу в узел (p,q) передвижение ее прекращается. Достижение границы происходит с вероятностью 1. Поэтому можно записать:

Суммирование производится по всем узлам границы.

При этом, если узел (i,j) находится на границе то:

1. P(i,j;p,q) = 1 , если i = p и j = q.

2. P(i,j;p,q) = 0 , если (i,j) и (p,q) различные граничные узлы.

Рассмотрим величины Vij , определяемые по формуле:

Здесь узел Mpq пробегает всю границу G. Введем функцию f(p,q) = u(p,q) на границе G. Функцию f(p,q) можно считать случайной величчиной, когда создаются различные траектории. Тогда согласно (5) Vi,j будет являться математическим ожиданием f(p,q) для траекторий начинающихся в узле (i,j) и заканчивающихся на границе G. Частица из узла (i,j) после первого шага с вероятностью ¼ перемещается налево, направо, вверх или вниз. Случайное блуждание распадается на 4 новых случайных блуждания. Используем формулу полной вероятности:

Умножим левую и правую части (6) на граничные значентя f(p,q) и просуммируем по всем значениям на границе. Получим:

Также:

Vp,q = f(p,q), так как P(p,q;p,q)= 1. (8)

Решить задачу Дирихле можно также методом конечных разностей. Это рассмотрено в вышеприведенном разделе "Решение задачи Дирихле с помощью конечных разностей". Повторим некоторые из этих рассуждений.

Для этого дифференциальное уравнение (1) в каждой внутренней точке x(i,j) (или узле) можно заменить соотношением:

Если правые части просуммировать, прировнять 0 и перенести функцию u(x, y) в правую часть, то получим соотношение:

Это же соотношение (для оценки точности) можно получить с помощью разложения в ряд Тейлора.

Заменяя u(xi,j) на ui,j получаем:

up,q = f(p,q) (10)

Формулы (7), (8) совпадают с (9), (10) с точностью до обозначений. Поэтому ui,j можно рассматривать как математические ожидания Vi,j . Эти уравнения аналогичны рассмотренным в (2) и (3), которые определчют неоднородную систему. Если выполнить большое число N равномерных случайных блужданий (траекторий) исходящих из определенного узла Mij то функцию ui,j можно получить так:

Суммирование производится по траекториям. Каждая траектория заканчивается на границе. Аргументы для функции f находятся на границе сетки.

Здесь различные p, q зависят от траектории и повторяются много раз при большом N. Несколько траекторий может заканчиваться в одном узле границы.

Опираясь на изложенную теорию, написана программа в системе программирования C++Bulder, которая решает задачу Дирихле с помощью случайного блуждания. В качестве функции u(x,y) выброна функция из (Демидович, 1967, с. 272): u(x,y) = x * y / 4

Программная реализация решения задачи Дирихле

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

Ниже рассматриваются основные особенности программы написанной на C++ (Архангельский, 2008; Архангельский, 2010; Стивен Гильберт 2003) и решающей задачу Дирихле с помощью случайного блуждания на плоскости.

В программе задача решается для четырех различных областей:

• -квадрат 9х9 (h =0,5)

• -треyгольник из квадрата 9х9 (h =0,5)

• -квадрат 17х17 (h =0,5)

• -круг из квадрата 17х17 (h =0,5)

Напомним вид уравнения Лапласа:

Для демонстрации решения взята функция u= x*y / 4.

У всех областей кусочно-гладкая граница. Сетка строится с помощью переменных (таблиц) типа TStringGrid. Ширина и высота ячейки равна h. Считается, что правая сторона таблицы и нижняя сторона таблицы не содержат узлов. Вычисления и обработка данных для всех областей производится похожим образом. В ячейках сетки часто располагаются значения функции u(x,y) для внутренних или граничных узлов. Это поясняется ниже. Предполагается, что узел, соответствующий значению в ячейке, расположен в левом верхнем углу ячейки.

Ниже на рис. 3 представлен начальный вид отображаемый в программе при расчете области «квадрат 9х9».

Начало координат находится в левом нижнем углу пересечения осей координат. Поддерживается адресация по номерам ячеек. При рассмотрении траекторий удобно пользоваться номерами ячеек. Левый верхний угол ячейки (0,0) соответствует началу координат. Ячейка (1,1) соответствует координате (h,h).

Задача решается примерно одинаково для различных областей, но зависит от их формы и величины шага h. В настоящее время программа работает с функцией: u= x*y / 4.

Рассмотрим подробнее случай с квадратом 9х9.

Граничные узлы находятся в левом верхнем углу ячеек с числами. Числа дают значения функции u(x,y) в граничных узлах. Внутренние узлы находятся в левом верхнем углу ячеек с буквой В.

Обеспечивается возможность расчета значения функции u(x,y) для любого выбранного внутреннего узла (ячейки). Для этого нужно номер соответсвующей ячейки задать в полях "x:" , "y:" и нажать кнопку "go1" (смотрите рис. 5). Задан номер узла (ячейки): x=2, y=2. Получено значение 0,24891 (левая таблица), а точное значение подсчитанное по формуле u= x*y / 4, равно 0,25 (правая таблица). Количество генераций для датчика случайных чисел при создании траекторий задано равным 60000.

При нажатии "all" производится расчет значений для всех внутренних узлов. При этом описанное выше действие для одного узла повторяется для всех внутренних узлов (смотрите рис. 5). 

Как происходит расчет значения функции u(x,y) для одного внутреннего узла

Во всех граничных узлах (левая таблица) определяется значение u(x,y). Нахождение значения u(x,y) для любого внутреннего узла (x,y) для левой таблицы производится следующим образом. С помощью датчика случайных чисел генерируются случайные траектории исходящие из узла (x,y). С вероятностью 1 траектория достигает границы области. После чего генерация этой траектории заканчивается и начинается новая траектория из узла (x,y). Граничное значение u(x,y) траекторий рассматривается как случайная величина. Математическое ожидание найденных граничных значений принимается в качестве результирующего значения для узла в соответствии с формулой (11).

Обработка полученных данных

При нажатии кнопки "разность" отображается таблица с разностями между полученными значениями и аналитическими для каждого узла, а при нажатии кнопки "траектории" отображается таблица с количеством траекторий для каждого внутреннего узла:

На рис. 8 отображается выдача статистистической информации после нажатия кнопки "статистика". Каждая строка соответствует одному прогону (кнопка "all") .

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

В поле max(разность) указывается максимальная разность для ячеек между подсчитанным и правильным значениями.

В поле "ср кв разн" указывается средний квадрат разности для ячеек.

В поле "ср абс разн" указывается средняя абсолютная разности для ячеек.

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

Выполнение вычисления для фигуры в виде треугольника

Аналогичные вычисления можно получить для фигуры типа "Треугольник" (см рис. 9).

Граничные узлы расположены по границам треугольника. Выдача статистистической информации для области "треугольник" отображена на рис. 10.

 

 

Выполнение вычислений для круга

При проведении расчетов для области "Квадрат Круг" появляются 2 таблицы, изображенные на рис. 11. Числа в ячейках левой таблицы задают значения для граничных узлов. В остальные ячейки заносится символ "*". При генерации траектории из внутреннего узла выход на ячейку, содержащую число означает выход на граничный узел и генерация этой траектории заканчивается.

На рис. 12 отображены расчитанные значения для круга.

Вышеприведенные вычисления (см. рис. 12) проводились при генерации 60000 случайных значений для каждого узла. На рис. 13 показаны статистические данные при вычислении серии : сначала при генерации 60000 случайных значений для каждого узла, затем 6000 и в конце 120000 случайных значений для каждого узла. С увеличением числа генераций случайных значений (а значит и траекторий) точность увеличивается.

Выполнение вычислений для квадрата 17х17

Проведены аналогичные вычисления для области "Квадрат17x17" . Такие же как для области “квадрат 9х9”. Рузультаты приведены на рис. 14.

Заключение

В статье рассмотрена программная система, которая решает задачу Дирихле на плоскости с помощью генерации случайных траекторий. Система создана для демонстрации студентам математических методов статистики. Приводится решение для треугольника, круга и двух квадратов разного размера. Программа реализована на C++. Программа позволяет генерировать различное число случайных траекторий, измерять отличие полученных результатов от аналитических величин. Проведено сравнение полученных результатов с аналитическими значениями. Результаты показывают, что при увеличении числа траекторий ошибка (максимальная величина отклонения от аналитических значений) уменьшается. Планируется сравнить полученные результаты с аналогичными данными получаемыми при решении подобных задач разностными методами.

1. Андреев, В.Б. (2020). Численные методы: учебное пособие (2-е изд.). Издательство МАКС Пресс.
   Andreev, V.B. (2020). Chyslennye metody: uchebnoe posobie [Numerical Methods: Textbook] (2th ed.). Moscow MAKS Press. (In Russ.).                        
2. Архангельский, А.Я. (2008). Компоненты C++Builder. Справочное пособие. Издательство: Бином.
   Arhangel’skii, A.Ja. (2008). Komponenty C++Builder. Spravochnoe posobie. [C++Builder. Components Reference Guide.]. Binom. (In Russ.).
3. Архангельский, А.Я. (2010). Программирование C++Builder (7-е изд.).  Издательство: Бином.
   Arhangel’skii, A.Ja. (2010). Programmirovanie C++Builder [C++Builder Programming] (7th ed.). Binom. (In Russ.).
4. Бахвалов, Н.С., Жидков, Н.П., Кобельков, Г.П. (2023). Численные методы (11-е изд.). Издательство Лаборатория знаний.
   Bakhvalov, N.S., Zhidkov N.P., Kobelkov G.P. (2023). Chyslennye metody [Numerical Methods] (11th ed.). Moscow Laboratory of Knowledge. (In Russ.).
5. Демидович, Б.П. (1967). Численные методы анализа: приближение функций, дифференциальные и интегральные уравнения (3-е изд.).  Издательство М. Наука.
   Demidovich, B. P. (1967).  Chyslennye metody analiza: priblizhenia funkcii, differencialnye i integralnye uravnenia [Numerical Methods of Analysis: Approximation of Functions, Differential and Integral Equations] (3th ed.).  Moscow Nauka. (In Russ.).
6. Демидович, Б.П., Марон И.А., Шувалова Э.З. (2010). Численные методы анализа: приближение функций, дифференциальные и интегральные уравнения (5-е изд.). Издательство Санкт-Петербург.
   Demidovich, B.P., Maron I.A., Shuvalova E.Z. (2010).  Chyslennye metody analiza: priblizhenia funkcii, differencialnye i integralnye uravnenia [Numerical Methods of Analysis: Approximation of Functions, Differential and Integral Equations] (5th ed.).  Saint-Petersburg. (In Russ.).
7. Самарский, А.А. (1971) Введение в теорию разностных схем. Издательство M. Наука.
   Samarskii, A.A. (1971) Vvedenie v teoriu raznostnyh shem [Introduction to the Theory of Difference Schemes]. Moscow Nauka. (In Russ.).
8. Самарский, А.А. (1977). Теория разностных схем: изд-во. Издательство M. Наука.
   Samarskii, A.A. (1977). Theoriia raznostnyh shem [Theory of Difference Schemes]. Moscow Nauka. (In Russ.).
9. Стивен Гильберт, Билл Маккарти (2003). Самоучитель VISUAL C++ 6 в примерах. Издательство М. DiaSoft.
   Stephen Gilbert, Bill McCarthy (2003). Samouchitel VISUAL C++ 6 v primerah. [VISUAL C++ 6 Tutorial in examples].  M. DiaSoft (In Russ.).
10. Тихонов, А.Н. (2002). Дифференциальные уравнения (4-е изд.). Издательство М. ФИЗМАТЛИТ.
    Tikhonov, A.N. Differentialnye uravneniya [Differential Equations] (4th ed.).  Moscow FIZMATLIT. (In Russ.).
11. Тихонов, А.Н., Самарский, А.А. (2004). Уравнения математической физики (7-е изд.). Издательство Московского университета.
    Tikhonov, A.N., Samarskii, A.A. (2004). Uravneniya matematicheskoi fiziki [Equations of Mathematical Physics] (7th ed.). Moscow State University. (In Russ.).
12. Яковлев В.Б. (2008) Вычислительная математика: Учебное пособие. Издательство М. МИЭТ.
    Yakovlev V.B. (2008) Vychislitelnaia matematika: uchebnoe posobie [Computational Mathematics: A Textbook]. Moscow MIET. (In Russ.).