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

Глава1

Введение

О чём эта работа?

Эта работа представляет собой теорию и практику численного решения задачи Дирихле для уравнения Лапласа.

Из чего же состоит эта работа?

Во-первых, это само введение, которое вы сейчас и читаете. Оно призвано кратко объяснить суть работы.

Во-вторых, это теоретическая часть. Здесь я определяю основные термины и рассказываю, что в принципе делается в третьей части — в практической части.

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

Теперь я хотел бы сказать пару слов о тексте самой работы и перейти к изложению самой сути.

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

Я, как автор данного реферата, заявляю о лицензировании работы под лицензией CC BY-SA. Это означает, что вы можете свободно распространять и/или изменять части работы так, как вам вздумается при соблюдении двух условий: ссылаться на авторскую работу и распространять продукт под той же лицензией, под которой вы его получили, то есть под CC BY-SA.

В тексте я буду использовать фрагменты, взятые из открытых источников, в том числе и из русской Википедии, что прекрасно согласуется с условиями выбранной лицензии. Такие ссылки я буду обозначать буквой «w» и иным шрифтом, что будет означать ссылку на статью, одноимённую слову на момент написания реферата, то есть на 27.04.2013.

Приступим.

Теория

При исследовании стационарных процессов вроде явлений теплопроводности, диффузии и т. д. обычно получаются дифференциальные уравнения второго порядка, частным случаем которых является уравнение Лапласа.

Так как в большинстве случаев нас интересуют численные оценки явления, то мы рассматриваем решение уравнение Лапласа на какой-то области, а следовательно, задаём некие граничные условия системы; а тогда исходная задача называется задачей Дирихле для уравнения Лапласа.

Примечание:

Уравнение Лапласа — дифференциальное уравнение в частных производных:

\Delta u(x_{1},x_{2},\dots ,x_{n})=0,

а в трёхмерном пространстве уравнение Лапласа называется уравнением Гельмгольца.

Задача Дирихле — задача отыскания в области

D

евклидова пространства гармонической функции, которая на границе

{\mathfrak {U}}

области

D

совпадает с наперёд заданной непрерывной функцией

f

.

В целом ряде случаев дифференциальные уравнения не имеют полного аналитического решения, то есть вы не можете просто сесть и найти точное решение на бумаге, но можете «в лоб» высчитать нужную вам функцию на компьютере.

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

Да простят меня математики за следующую фразу, но вещественных чисел «больше», чем рациональных, а если говорить точнее, то у множеств ${\mathfrak {Q}}$ и ${\mathfrak {R}}$ разные мощности. Но так как компьютер оперирует только с числами с конечным количеством знаков после запятой, то он неизбежно сводит ${\mathfrak {Q}}$ к ${\mathfrak {R}}$ , а, следовательно, теряет все знаки после запятой начиная с энного, что приводит к погрешности.

Второе, на что стоит обратить внимание — округление входящих данных самим пользователем. Никто не будет, например, прописывать огромное количество знаков числа $0,(3)$ ; обычно ограничиваются чем-то вроде « $0,333333333333$ », что меньше $0,(3)$ .

Третье — не на всех компьютерах, находящихся в массовой продаже можно сделать достаточное число итераций из-за ограничений в размере жёсткого диска, оперативной памяти или — и это особенно важно — из-за ограничений в скорости работы процессора, а, значит, для того, чтобы получить решение уравнения раньше заката, нужно увеличивать и без того большие ошибки компьютера малой точностью вычислений.

Практическая часть

Итак, задача. Решить методом сеток (я потом объясню, что это такое) двумерное уравнение Лапласа внутри области Ԇ, ограниченной двумя кривыми

{\mathfrak {U_{1}}}

:

x^{2}+y^{2}-x=16(x^{2}+y^{2});

{\mathfrak {U_{2}}}

:

x^{2}+y^{2}=1;

с граничными условиями

\partial \mid u({\mathfrak {U_{1}}})=0;

\partial \mid u({\mathfrak {U_{2}}})=\pi .

Я решил использовать декартову систему координат из-за простоты реализации. Я знаю, что можно решать и в сферических, но корректного решения мне получить в них не удалось.

Для решения данной задачи Дирихле используется конечно-разностный метод, в котором (опуская промежуточные выкладки) уравнение Лапласа $\Delta u(x,y)=0$ записывается в виде $u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)=4u(x,y)$ . То есть значение функции в данной точке выражается через среднее значение её соседей по сетке.

Для решение задачи методом сеток необходимо покрыть данную нам область ${\mathfrak {U}}$ сеткой из прямых линий, при этом точки пересечений линий — узлы — могут лежать как в области ${\mathfrak {U}}$ , так и вне неё.

Сначала вычислим значение функции в точках на границе. Ну, то есть мы-то с вами их знаем, а компьютер — ещё нет; слово «вычислить», с вашего позволения, я иногда буду использовать в смысле «записать в память компьютера». После этого заполним массив значений функции внутри области некоей константой, пусть средним арифметическим всех точек границы, конкретное значение константы не так важно, как кажется на первый взгляд; важно, чтобы оно не сильно отличалось от порядка значений на границе. Назовём это распределение значений функции 0-системой.

А теперь начинаем, собственно говоря, считать нашу функцию. Берём каждую и-житую точку, считаем в ней значение функции как среднее арифметическое её соседей. Делаем так для всех точек на области и получаем 1-систему.

А теперь повторим эти подсчёты для области ${\mathfrak {U}}$ ещё несколько раз до достижения нужной точности, получив в конце концов n-систему, где n — число повторов.

Сходимость метода пропорциональна числу узлов сетки.

Практическая задача

Поставленная задача решалась с помощью программы, написанной на C#. Количество итераций — 1500. Из-за ограничений в вычислительных мощностях у меня дома мне не удалось просчитать весь график целиком, но лишь часть его — один сектор. Графики строились с помощью gnuplot.

Вывод

В данной работе была дана теория одного из численных метод решения уравнений в частных производных эллиптического типа – метод сеток, а в частности метод конечных разностей. Была написана программа по реализации этого метода, получены и проанализированы результаты решения задачи Дирихле уравнения Лапласа в заданной области.

Список литературы

Положий Г.П. — Уравнения математической физики;
Корн Г., Корн Т. — Справочник по математике для научных работников и инженеров;
ru.wikipedia.org.

Приложение

namespace DirichletProblem {

   public class Point
   {
       public double[] x = new double[2000];
       public double[] y = new double[2000];
   }
   public partial class Form1 : Form
   {
       public Form1()
       {
           InitializeComponent();
       }
       private void makeComputing(object sender, EventArgs e)
       {
           const Int64 sizeOfUArray = 2001;
           Point points = new Point();
           double[,] u = new double[sizeOfUArray, sizeOfUArray];
           double[,] uTemp = new double[sizeOfUArray, sizeOfUArray];
           int i, j;
           int n = 1000; // итерации
           int k = -n;
           int sizeOfPointsArray;
           double h = 0.05; // шаг
           int iterationCount = 0;            
           for (i = 0; i < 2 * n; i++) // заполняем узлы сетку
           {
               points.x[i] = k * h;
               points.y[i] = k * h;
               k++;
           }
           sizeOfPointsArray = points.x.Count(); 
           for (i = 0; i < sizeOfPointsArray; i++) // отбираем точки, лежащие внутри области и на ее границе; считаем U
           {
               for (j = 0; j < sizeOfPointsArray; j++)
               {
                   double x = points.x[j];
                   double y = points.y[i];
                   if (isInAreaIncludeBorders(x, y))
                   {
                       u[j, i] = borders(x, y);
                   }
               }
           }
           for (int g = 1; g <= 20; g++)
           {
               iterationCount++;
               for (i = 1; i < sizeOfPointsArray; i++) // отбираем точки, лежащие внутри области; считаем U для внутренних точек
               {
                   for (j = 1; j < sizeOfPointsArray; j++)
                   {
                       double x = points.x[j];
                       double y = points.y[i];
                       if (isInAreaExcludeBorders(x, y))
                       {
                           double u1, u2, u3, u4;
                           u1 = u[j + 1, i];
                           u2 = u[j - 1, i];
                           u3 = u[j, i + 1];
                           u4 = u[j, i - 1];
                           u[j, i] = (u1 + u2 + u3 + u4) / 4;
                       }
                   }
               }
           }
           TextWriter tw = new StreamWriter("data.txt");
           for (i = 0; i < sizeOfPointsArray; i++)
           {
               for (j = 0; j < sizeOfPointsArray; j++)
               {
                   if (u[j, i] != 0)
                       tw.WriteLine(points.x[j].ToString() + " " + points.y[i].ToString() + " " + u[j, i].ToString());
               }
           }
           tw.Close();
       }
       private double borders(double x, double y)
       {
           return Math.Sqrt(y);
       }
       private bool isInAreaIncludeBorders(double x, double y)
       {
           if (((y >= -Math.Sqrt(8 + x - x * x - 4 * Math.Sqrt(4 + x))) || (y <= Math.Sqrt(8 + x - x * x - 4 * Math.Sqrt(4 + x))) || (y >= -Math.Sqrt(8 + x - x * x + 4 * Math.Sqrt(4 + x))) || (y <= Math.Sqrt(8 + x - x * x + 4 * Math.Sqrt(4 + x)))) && ((y >= Math.Sqrt(1 - x * x)) || (y >= Math.Sqrt(1 - x * x))))
           {
               return true;
           }
           else
           {
               return false; 
           }
       }
       private bool isInAreaExcludeBorders(double x, double y)
       {
           if (((y > -Math.Sqrt(8 + x - x * x - 4 * Math.Sqrt(4 + x))) || (y < Math.Sqrt(8 + x - x * x - 4 * Math.Sqrt(4 + x))) || (y > -Math.Sqrt(8 + x - x * x + 4 * Math.Sqrt(4 + x))) || (y < Math.Sqrt(8 + x - x * x + 4 * Math.Sqrt(4 + x)))) && ((y > Math.Sqrt(1 - x * x)) || (y > Math.Sqrt(1 - x * x)))) // 
           {
               return true;
           }
           else
           {
               return false;
           }
       }
       private void clearSquareArray(double[,] array)
       {
           int size = array.GetLength(0);
           for (int i = 0; i < size; i++)
           {
               for (int j = 0; j < size; j++)
               {
                   array[i, j] = 0;
               }
           }
       }
   }

}