МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ
БЕЛОРУССКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Факультет прикладной математики и информатики
Кафедра вычислительной математики


БАСАЛАЙ ДМИТРИЙ НИКОЛАЕВИЧ

БЕССЕТОЧНЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЯ ПУАССОНА

Отчет по преддипломной практике

студента 5 курса 5 группы

“Допустить к защите“ с предварительной оценкой ___ Руководитель практики _____________________ “___” ___________ 2010 г

Руководитель практики Никифоров Иван Васильевич, доцент кафедры Выч. мат., канд. физ.-мат. наук

Минск 2010

БЕЛОРУССКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Факультет прикладной математики и информатики

Кафедра вычислительной математики

“Утверждаю”

Заведующий кафедрой

_______________ П.А. Мандрик

“___”  _______________  2010 г.

ЗАДАНИЕ

ПО ПОДГОТОВКЕ ПРЕДДИПЛОМНОЙ ПРАКТИКИ

Студенту 5 курса Басалаю Дмитрию Николаевичу

1. 
   Тема работы БЕССЕТОЧНЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЯ ПУАССОНА
   
2. 
   Срок сдачи студентом законченной работы 20 февраля 2010 г.
   
3. 
   Исходные данные к работе
   

• 
  Методические материалы.
  
• 
  Размещенные в электронной библиотеке методические материалы, примеры из обучающих программных модулей.
  
• 
  Технические требования к электронным версиям отчетных документов, рекомендации составителям.
  

Библиографические описания источников, рекомендуемых студентам к ознакомлению при выполнении работы (для изучения предметной части задания, как правило, достаточно ознакомиться с любой из перечисленных в начале списка книг):

1. 
   J. Kuhnert, S. Tiwari. Grid free method for solving the Poisson equation. – Berichte des Fraunhofer ITWM, Nr. 25, 2001. – 20 с.
   
2. 
   T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl. Meshless methods. An overview and recent developments. 2001. – 77 с.
   
3. 
   К.Е. Афанасьев, И.В. Григорьева. Современное состояние бессеточных методов. Кемеровский государственный университет, 2007. – 7 с.
   
4. 
   [Электронный ресурс] http://eqworld.ipmnet.ru
   
5. 
   [Электронный ресурс] http://ru.wikipedia.org
   

4. 
   Перечень вопросов подлежащих разработке или краткое содержание работы
   

изучение литературы по решению уравнения Пуассона бессеточными методами, построение алгоритмов численного решения, проведение вычислительного эксперимента, анализ полученных результатов, оформление дипломной работы согласно ГОСТа.

5. 
   Перечень графического материала
   

• 
  Логотип БГУ для включения на слайды презентации.
  

6. 
   Дата выдачи задания «23» ноября 2009г.
   
7. 
   Календарный график работы на весь период (с указанием этапов работы и сроков их выполнения)
   

• 
  23.11.09 – 23.01.10 – ознакомление с предлагаемыми темами, выбор и согласование темы с руководителем, регистрация на кафедре Выч. мат., ознакомление с техническими требованиями к электронным отчетным документам и освоение правил, как их реализовать, изучение постановки задачи, основных теоретических вопросов;
  
• 
  23.01.10 – 20.02.10 – практическая реализация задач проекта, оформление результатов работы (отчета, презентации), подготовка доклада и отладка презентации на защиту;
  
• 
  20.02.10 – защита.
  

Руководитель ____________________ / ____________________ /

(подпись) (Ф.И.О.)
Задание принял к исполнению “___” ___________ 2009г.

______________________

(подпись студента)

Аннотация

Главной целью работы является изучение бессеточного метода решения уравнения Пуассона.

Анатацыя

Галоўнай мэтай дадзенай работы з’яўляецца вывучэнне бяссетачнага метада рашэння ураўнення Пуасона.

Annotation

The main objective of this work is learning of the grid free method for Poisson equation solving.

Реферат

Отчет по преддипломной практике, 37 с., 18 рисунков, 5 источников.

БЕССЕТОЧНЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЯ ПУАССОНА

Объект исследования – уравнение Пуассона, а также бессеточный метод его решения.

Цель работы – найти решение уравнения Пуассона, используя бессеточный метод.

Результаты работы – найдено решение уравнения Пуассона.

Содержание

Введение

Бессеточные методы были первоначально разработаны для решения модельных задач динамики жидкости. Это так называемые методы частиц. В рассматриваемых моделях жидкость заменяется дискретным набором точек или частиц, которые движутся вместе с ней. Изначально так называемый метод сглаженных частиц (SPH-метод) был разработан для решения астрофизических задач без границы.SPH-метод был распространен решение различных задач динамики жидкости.

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

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

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

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

Метод устойчив, и численное решение сходится к фиксированной точке при числе итераций, стремящемся к бесконечности. Он может быть применен как г грубому, так и более тонкому распределению точек. Скорость сходимости выше на более грубой сетке.

1Численный метод

1.1Аппроксимация функции и ее производных методом наименьших квадратов

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

Пусть скалярная функция, и ее значения в точках , где – общее число точек в . Мы решаем задачу аппроксимации функции и ее производных в точке , используя значения в соседних точках. Чтобы ограничить количество точек, мы используем весовую функцию с параметром . Константа подобна длине сглаживания в классическом методе сглаженных частиц. Весовая функция может быть произвольной. В наших вычислениях будем использовать весовую функцию Гаусса в следующей форме:

(1.1)

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

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

Рассмотрим разложение по формуле Тейлора функции в окрестности точки :

(1.2)

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

(1.3)

где

(1.4)

, , .

(1.5)

Символами , , обозаначены , , , соответственно, .

При эта система переопределена для десяти неизвестных .

Неизвестный вектор получается из взвешеннго метода наименьших квадратов минимизацией квадратичной формы

(1.6)

которая может быть записана в следующей форме

(1.7)

где

(1.8)

Минимизация формально дает результат

(1.9)

1.2Аппроксимация методом наименьших квадратов для уравнения Пуассона

Используя этот метод, мы не дискретизируем уравнение Пуассона непосредственно как в классических методах. Мы решим его, используя итерационный процесс. Рассмотрим следующее уравнение Пуассона

	(1.10)
	(1.11)
	(1.12)

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

В предыдущем пункте мы рассмотрели метод наименьших квадратов для приближения функции и ее производных в произвольной точке, используя ее значение в соседних точках. Теперь ситуация немного другая, так как значения в дискретных точках неизвестны априори. Это значит, что если вектор в формуле (1.9) неизвестен, мы не можем определить коэффициенты вектора . Поэтому мы задаем начальное значение функции во всех точках. Теперь рассмотрим задачу нахождения в произвольной точке по ее соседним точкам . Как было показано в предыдущем пункте, снова рассмотрим разложение в точке

(1.13)

, где известные начальные значения. Также применяем условие того, что точка является решением уравнения (1.10). Для этого добавляем следующее уравнение к упомянутым выше уравнениям:

(1.14)

Для граничных условий Неймана (1.12) мы имеем еще одно уравнение для точки , лежащей на границе

(1.15)

где , , – компоненты нормального вектора в граничной точке . Следовательно, мы имеем систему из уравнений с неизвестными и в общем случае .

Получим коэффициенты

(1.16)

для в точке минимизацией

(1.17)

Точно также минимизация дает

(1.18)

где матрицы и векторы отличаются от (1.4), (1.5) и имеют вид

(1.19)

(1.20)

,

(1.21)

Отметим, что можно использовать уравнение (1.14) в (1.13) для того, чтобы уменьшить количество неизвестных. Например, из (1.14) имеем

(1.22)

Подставляем в (1.13) и уменьшаем число неизвестных до для внутренних точек. Следовательно, вместо матицы нужно инвертировать матрицу . Это, конечно, не сильно уменьшает вычислительные затраты. В этом случае коэффициенты матрицы и вектора отличаются от предложенных выше. Но при этом в обоих случаях скорость сходимости остается одинаковой.

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

Завершаем итерации, если ошибка удовлетворяет следующему неравенству

(1.23)

и решение определяется как при стремящемся к бесконечности. Параметр – очень малая положительная константа, варьирующаяся в зависимости от задачи и значения . Параметр обычно выбирают согласно топологии. Мы рассмотрели постоянную . Скорость сходимости будет больше, если выбрать параметр большим, при этом увеличив расстояние между точками.

2Программа

2.1Общее описание

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

2.2Обзор основного функционала и работа с программой

Программа представляет собой веб-приложение (см. рисунок 2.1).

Рисунок 2.1 – Демонстрационное приложение

Для решения необходимо выбрать тип уравнения Пуассона:

Рисунок 2.2 – Выбор типа уравнения

размерность задачи (двумерная или трехмерная):

Рисунок 2.3 – Выбор размерности задачи

метод задания сетки:

Рисунок 2.4 – Выбор метода задания сетки

и наконец, ввести входные параметры метода:

Рисунок 2.5 – Входные параметры

Если все введенные данные корректны, то по окончанию рассчетов полученные результаты выведутся на экран (см. рисунок 2.6), будет построен график решения (см. рисунок 2.7):

Рисунок 2.6 – Решение

Рисунок 2.7 – График решения

Для анализа результатов выполнения программы с разными входными данными предусмотрена возможность просмотра статистики:

Рисунок 2.8 – Статистика

2.3Листинг

Ниже приведен листинг основной части приложения, отвечающей за решение уравнения. Используемая платформа разработки – Microsoft .NET Framework 4.0, язык программирования C#.

Отметим, что для нахождения обратной матрицы использовался метод LU-факторизации.

public class GridFreeMethodFor2D

{

#region Public properties

public Problem2D Problem { get; set; }
public bool AutoGeneratePoints { get; set; }
#endregion
#region Private properties
private List PrevU { get; set; }
#endregion
#region Constructors and Initializers
public GridFreeMethodFor2D(Problem2D problem): this(problem, true)

{

}

public GridFreeMethodFor2D(Problem2D problem, bool autoGeneratePoints)

{

AutoGeneratePoints = autoGeneratePoints;

Problem = problem;

Problem.Solution = new Solution2D();

Problem.ExactSolution = new Solution2D();
if (AutoGeneratePoints)

{

Problem.Solution.U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2);

}

}
public GridFreeMethodFor2D(int n,

double x1, double x2,

double y1, double y2,

double g, double epsilon, double h, double alpha)

{

Problem = new Problem2D

{

InputParameters = new InputParameters2D

{

N = n,
X1 = x1,

X2 = x2,

Y1 = y1,

Y2 = y2,
G = g,

Epsilon = epsilon,

H = h,

Alpha = alpha

},

Solution = new Solution2D

{

IterationsCount = 0

}

};
Problem.Solution.U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2);

}
public static double F(double x, double y)

{

return 1;

}
#endregion
#region Private methods
private Matrix GetM(Point3D p, IList neighboringPoints)

{

int n = neighboringPoints.Count;

var m = new Matrix(n + 1, 6);

for (int i = 0; i < n; i++)

{

var point = neighboringPoints[i];

m[i, 0] = 1;

m[i, 1] = point.X - p.X;

m[i, 2] = point.Y - p.Y;

m[i, 3] = 0.5 * (point.X - p.X) * (point.X - p.X);

m[i, 4] = (point.X - p.X) * (point.Y - p.Y);

m[i, 5] = 0.5 * (point.Y - p.Y) * (point.Y - p.Y);

}

m[n, 0] = Problem.Type == ProblemType.TypeThree ? -1 : 0;

m[n, 1] = 0; m[n, 2] = 0;

m[n, 3] = 1; m[n, 4] = 0; m[n, 5] = 1;

return m;

}
private Matrix GetW(Point3D p, IList neighboringPoints)

{

int n = neighboringPoints.Count;

var w = new Matrix(n + 1, n + 1);

for (int i = 0; i < n; i++)

{

var point = neighboringPoints[i];

double norm = (new Vector(point) - new Vector(p)).Norm();
if (norm / Problem.InputParameters.H <= 1)

{

w[i, i] = System.Math.Exp(-Problem.InputParameters.Alpha * norm * norm / (Problem.InputParameters.H * Problem.InputParameters.H));

}

else

{

w[i, i] = 0;

}

}

w[n, n] = 1;

return w;

}
private Vector GetB(Point3D p, IList neighboringPoints)

{

int n = neighboringPoints.Count;

var b = new Vector(n + 1);
for (int i = 0; i < n; i++)

{

var point = neighboringPoints[i];

b[i] = GetPoint(point, neighboringPoints).Z;

}

b[n] = F(p.X, p.Y);

return b;

}
public void Solve()

{

var startTime = DateTime.Now;

Problem.Solution.IterationsCount = 0;

do

{

PrevU = new List(Problem.InputParameters.N);

foreach (Point3D p in Problem.Solution.U)

{

PrevU.Add(new Point3D(p));

}
for (int i = 0; i < Problem.Solution.U.Count; i++)

{

Point3D point = Problem.Solution.U[i];

if (!point.IsBoundary)

{

var neighboringPoints = GetNeighboringPoints(point);

var m = GetM(point, neighboringPoints);

var w = GetW(point, neighboringPoints);

var b = GetB(point, neighboringPoints);
var mt = m.Transpose();

var tmp = mt * w;

var a = (tmp * m).Inverse() * tmp * b;

point.Z = a[0];

}

else

{

point.Z = Problem.InputParameters.G;

}

}

Problem.Solution.IterationsCount++;

} while (!IsIterationsStopped());
var endTime = DateTime.Now;
Problem.Solution.ExecutionTime = endTime - startTime;

}
private bool IsIterationsStopped()

{

double summ1 = 0;

double summ2 = 0;
for (int i = 0; i < Problem.Solution.U.Count; i++)

{

summ1 += System.Math.Abs(Problem.Solution.U[i].Z - PrevU[i].Z);

summ2 += System.Math.Abs(Problem.Solution.U[i].Z);

}

return summ1 / summ2 < Problem.InputParameters.Epsilon;

}
private Point3D GetPoint(Point3D point, IEnumerable points)

{

var list = points

.Where(p => p.X == point.X && p.Y == point.Y)

.ToList();

if (list.Count > 0)

{

return list[0];

}

return null;

}
private List GenerateRegularPoints(int count, double x1, double x2, double y1, double y2)

{

var points = new List(count);

int boundaryPointsCount = (int)System.Math.Pow(count, 1.0 / 2);
double hx = (x2 - x1) / boundaryPointsCount;

double hy = (y2 - y1) / boundaryPointsCount;

for (int i = 0; i <= boundaryPointsCount; i++)

{

for (int j = 0; j <= boundaryPointsCount; j++)

{

points.Add(new Point3D(x1 + i * hx, y1 + j * hy, 0,

i == 0 || j == 0 ||

i == boundaryPointsCount || j == boundaryPointsCount));

}

}
return points;

}
private IList GetNeighboringPoints(Point3D point)

{

var list = Problem.Solution.U

.Where(p => GetDistance(point, p) <= Problem.InputParameters.H &&

(p.X != point.X && p.Y != point.Y))

.ToList();
return list;

}
private static double GetDistance(Point3D a, Point3D b)

{

return System.Math.Sqrt((a.X - b.X) * (a.X - b.X) +

(a.Y - b.Y) * (a.Y - b.Y));

}
#endregion

}

3Тесты

Числовые эксперименты выполним в двумерном пространстве. Решим уравнение Пуассона на единичном квадрате, где известно аналитическое решение. Рассмотрим как регулярное, так и нерегулярное распределение точек. Регулярные точки сгенерированы с интервалом . Для нерегулярного случая границы области заменены набором точек с интервалом , а распределение внутренних точек выбрано произвольным неоднородным образом. Параметр в весовой функции выберем равным . В качестве начального приближения для возьмем во всех точках.

3.1Задача с граничными условиями Дирихле

Рассмотрим следующую задачу для уравнения Пуассона с граничными условиями Дирихле

(3.1)

Аналитическое решение этой задачи задается формулой:

(3.2)

Таблица 3.1 – Входные параметры.

Параметр	Значение
Количество точек	100
X1	0
X2	1
Y1	0
Y2	1
Значение функции на границе	0
	0.00001
	0.23
	0.00001

Результаты:

количество итераций – 28,

время выполнения – 00:00:00.7950455,

максимальная погрешность – 0.001502105.

Рисунок 3.1 – График решения

4Решения уравнения Пуассона при помощи MATLAB

MATLAB – это пакет прикладных программ для решения задач технических вычислений.

Partial Differential Equation (PDE) Toolbox из пакета MATLAB содержит средства для исследования и решения нестационарных дифференциальных уравнений второго порядка в частных производных. В пакете используется метод конечных элементов. Команды и графический интерфейс пакета могут быть использованы для математического моделирования PDE применительно к широкому классу инженерных и научных приложений, включая задачи сопротивления материалов, расчеты электромагнитных устройств, задачи тепломассопереноса и диффузии.

4.1Решение

Найдем решение задачи, рассмотренной в предыдуще пункте, используя PDE toolbox.

В командной строке MATLAB выполним команду pdetool. В результате откроется окно тулбокса:

Рисунок 4.1 –PDE toolbox

Зададим область :

Рисунок 4.2 – PDE toolbox: задание области

Зададим условия на границе:

Рисунок 4.3 – PDE toolbox: задание граничных условий

Проинициализируем сетку:

Рисунок 4.4 – PDE toolbox: инициализация сетки

Определим уравнение Пуассона:

Рисунок 4.5 – PDE toolbox: задание вида уравнения

Получим решение и построим его график:

Рисунок 4.6 – PDE toolbox: график решения

Полученные результаты можно экспортировать в файл.

4.2Сравнение результатов, полученных с помощью бессеточного метода и PDE toolbox

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

Рисунок 4.7 – Нерегулярная сетка

Получим решение, используя наше приложение:

Рисунок 4.8 – Сравнение полученных результатов

Рисунок 4.9 – Графики решений

Как видим, полученные решения отличаются лишь в 4 знаке после запятой, при этом в качестве входных параметров выбирались не самые оптимальные значения (для увеличения производительности), что говорит о высокой точности рассматриваемого метода.

Заключение

Итак, получена численная модель для решения уравнения Пуассона в бессеточной структуре. Рассматриваемый метод основан на методе наименьших квадратов. Он может применяться для любой эллиптической задачи, в особенности, если сетка имеет сложную нерегулярную структуру. Также было получено решение уравнения Пуассона с помощью пакета прикладных программ MATLAB. В сравнении с ним бессеточный метод показал отличный результат, что говорит о высокой степени эффективности его применения для подобных задач.

Список используемой литературы

1. 
   J. Kuhnert, S. Tiwari. Grid free method for solving the Poisson equation. – Berichte des Fraunhofer ITWM, Nr. 25, 2001. – 20 с.
   
2. 
   T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl. Meshless methods. An overview and recent developments. 2001. – 77 с.
   
3. 
   К.Е. Афанасьев, И.В. Григорьева. Современное состояние бессеточных методов. Кемеровский государственный университет, 2007. – 7 с.
   
4. 
   [Электронный ресурс] http://eqworld.ipmnet.ru
   
5. 
   [Электронный ресурс] http://ru.wikipedia.org