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

Глава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;
               }
           }
       }
   }

}