Решение эллиптических уравнений несколькими методами
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Филиал федерального государственного
автономного образовательного учреждения высшего профессионального образования
«Казанский (Приволжский) федеральный университет» в г. Набережные Челны
ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ
И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
КАФЕДРА ПРИКЛАДНОЙ МАТЕМАТИКИ И
ИНФОРМАТИКИ
Специальность: 010501.65 - Прикладная
математика и информатика
КУРСОВАЯ
РАБОТА
VIII СЕМЕСТР
ТЕМА:
«Решение эллиптических уравнений несколькими методами»
Дисциплина: численные методы
Выполнила
студентка Гатина Г.И.
группа 4606 курс 4
Научный руководитель
Марданшин Р.Г.
к. ф.-м. н., доцент
Набережные Челны 2009
Оглавление
Аннотация
Введение
Раздел 1. Математическое описание алгоритмов и операций
Раздел 2. Библиотека функций
Раздел 3. Тестирование
Вывод
Заключение
Список использованной литературы
Приложения
численное решение алгоритм эллиптическое
уравнение
Аннотация
В работе сначала приводятся основные понятия и математическое толкование
разностной схемы для уравнения Лапласа, далее приводятся разработанные в ходе
исследований методы. В третьем разделе описываются работы методов и выявляется
устойчивость различных разностных схем. Далее делается вывод о целесообразности
применении тех или иных схем и листинги разработанных методов.
Введение
Численное решение прикладных задач всегда интересовало математиков.
Крупнейшие представители прошлого сочетали в своих исследованиях изучение
явлений природы, получение их математического описания, как иногда говорят,
математической модели явления, и его исследование. Анализ усложненных моделей
потребовал создание специальных, как правило, численных или асимптотических
методов решения задач. Названия некоторых из таких методов - методы Ньютона,
Эйлера, Лобачевского, Гаусса, Чебышева, Эрмита, Крылова - свидетельствуют о
том, что их разработкой занимались крупнейшие ученые своего времени.
Настоящее время характерно резким расширением приложений математики, во
многим связанным с созданием и развитием средств вычислительной техники. В
результате появления ЭВМ (электронно-вычислительных машин, или как часто
говорят, компьютеров) с программным управлением менее чем за 50 лет скорость
выполнения арифметических операций возросла от 0.1 операции в секунду при
ручном расчете до 1012 операций на современных серийных ЭВМ, т.е.
примерно в 1013 раз.
В настоящее время разработка методов и алгоритмов решения задачи Коши для
обыкновенных дифференциальных уравнений продвинута настолько, что зачастую
исследователь, имеющий дело с этой задачей, не занимается выбором метода ее
решении, а просто обращается к стандартной программе.
В случае с уравнений с частными производными число принципиально
различных постановок задач существенно больше. В курсе уравнений с частными
производными обычно рассматривается незначительная часть таких постановок,
главным образом связанных с постоянными коэффициентами. При этом существует
очень малое количество задач, решаемых в явном виде. Многообразие постановок в
теории уравнений с частными производными связано с многообразием окружающего
нас мира.
Среди всех типов уравнений математической физики эллиптические уравнения
с точки зрения вычислителей стоят особняком. С одной стороны, имеется хорошо
развитая теория решения эллиптических уравнений и систем. Достаточно легко
доказываются теоремы об устойчивости разностных схем для эллиптических
уравнений. Цель работы: разработать сеточный метод, позволяющих решать
задачу Дирихле методом разностных схем на примере уравнения Лапласа. В качестве
среды разработки был выбран пакет matlab 6.5.
Раздел 1. Математическое описание
алгоритмов и операций
В данном разделе дается математическое толкование работы основных функций
и процедур библиотеки.
Рассмотрим сначала некоторые необходимые понятия из теории сеток:
Пусть
имеется пространство , где -
функция непрерывного аргумента . На
отрезке введем конечное множество точек , которое назовем сеткой. Точки , будем называть узлами сетки . Множество без
узлов и будем
обозначать . Если расстояние между
соседними узлами постоянно (не зависит от i), для всех , то
сетку называют равномерной (с шагом h), в противном
случае - неравномерной. Вместо функции ,
определенной для всех , будем рассматривать сеточную функцию , целочисленного аргумента или узла сетки , а заменим
конечномерным (размерностью N+1) пространством сеточных
функций. Очевидно, что сеточную функцию можно
рассматривать как вектор [1].
Можно
также провести дискретизацию и пространства функций многих переменных, когда - точка p-мерного евклидова пространства
(p>1). Так на плоскости можно ввести сетку ,
как множество точек (узлов) пересечения перпендикулярных прямых , , , где - шаги
сетки по направлениям и соответственно.
Сетка , очевидно равномерна по каждому из переменных в
отдельности. Вместо функции будем
рассматривать сеточную функцию
Многие стационарные физически задачи (исследование задач, связанных с
магнитными, гравитационными явлениями и теплопроводности) сводятся к решению
уравнения Лапласа вида
(1)
Решение для этого уравнения будем искать для некоторой ограниченной
области G изменения независимых переменных x, y. Границей области G является замкнутая линия L. Для
полной формулировки краевой задачи кроме уравнения Лапласа нужно задать
граничное условие на границе L.
Примем его в виде
(*)
Разностные схемы
В
области введем равномерную сетку
с шагами h по x и y. Заменяя производную по x и y
разностными выражениями
в уравнении (1) получаем:
(2)
В граничных условиях заменяем непрерывную функцию дискретной на области
прямоугольника:
С помощью данных уравнений можно записать систему линейных алгебраических
уравнений (2) в виде (схема «Крест»)
(3)
Система
(3) содержит 5 неизвестных. И образует систему линейных уравнений , где A - блочно-трехдиагональная матрица вида
, где
а
матрицы - диагональные матрицы с элементами на диагонали
равными единице. и размерности
.
Правая
часть системы имеет вид .
Для
решения такой системы можно применить метод Гаусса с исключением
несодержательных операций, т.е. обнуление элементов . В данном случае общее число операций составит
величину
Еще
одним из наиболее распространенных методом решения этой системы является
итерационный метод. Каждое уравнение запишем в виде, разрешенном относительно
значения в центральном узле:
(4)
Данную схему можно решить с помощью итерационных методов:
1. Гаусса-Зейделя,
(5)
2. Якоби,
3. Верхней релаксации.
(7)
В
алгоритме предусмотрен выбор начальных значений .
Итерационный процесс контролируется максимальным отклонением M
значений сеточной функции в узлах для последовательных итераций. Если его
величина достигнет некоторого допустимого значения точности , итерации прекращаются.
Существуют
и другие разностные схемы для решения уравнения Лапласа. Представив уравнение
(1) в виде (8), где t - переменная времени, то
исходная задача сводится к разностной схеме для уравнения теплопроводности, где
граничные условия совпадают с условиями .Условие (9)можно выбирать практически в произвольном виде,
согласованном с граничным условием.
Имеем
, отсюда можно найти явное выражение для значения
сеточной функции на (k+1)-м слое:
(10)
Граничные условия принимают вид
(11)
Вышеописанный метод называется «метод установления».
Процесс
численного решения представляет собой итерационный процесс решения задачи (8) с
условием (9) и (*) состоит в переходе от
произвольного значения (9) к искомому стационарному решению. Счет ведется до
выхода решения в стационарный режим. Естественно, ограничиваются решением, при
некотором достаточно большом , если
значения слоев совпадают с заданной степенью точности.
Условие
устойчивости вышеописанной схемы определяется неравенством (12)
Раздел 2. Библиотека функций
В работе разработаны следующие функции для решения
неоднородного одномерного уравнения теплопроводности:
1. ElipYa(X,N,fi,E) -реализует схему «Крест» при помощи
метода Якоби,
2. Elip1(X,N,fi,E)- решает схему «Крест» при помощи
метода Гаусса-Зейделя.
3. ElipR(X,N,fi,t,E) - решает схему «Крест» при помощи метода верхней
релаксации.
4. ElipGauss(X,N,fi,E) - решает схему
«Крест» методом Гаусса.
5. ElipT1(X,N,tau,fi,E) - сводит решаемую задачу к задаче теплопроводности и
реализует явную схему.
где x - длина стороны квадрата сетки, N - количество узлов сетки по x, fi - граничное условие, E - точность вычислений. В методе верхней релаксации t - параметр, а в явной схеме tau - шаг по времени.
Каждая
функция возвращает массивы, составляющие сетку: и и значения сеточной функции . В качестве начальных приближений используется
функция rand(I,J), которая возвращает матрицу случайных чисел
размерностью, указанной в скобках. Элементы матрицы находятся в диапазоне .
Раздел
3. Тестирование
Примем следующие параметры разработанных функций:
X=10
длина стороны квадрата;
N=15 -
кол-во узлов сетки;
fi=cos(x)+cos(y)- функция граничных условий;
E
=0.001 - точность;
t=1.2-
параметр для метода верхней релаксации;
tau=0.1
- шаг по времени для «метода установления».
Вводим начальные данные в функцию ElipYa - схема «крест» метод Якоби
Рисунок
1
Всего итераций 167.
Теперь с теми же параметрами произведем расчеты с помощью метода
Гаусса-Зейделя (Elip1).
Рисунок
2
Всего итераций 69.
Протестируем метод верхней релаксации (ElipT1)
Всего итераций 53.
Проверим метод Гаусса для решения схемы «крест»
Результаты совпадают с предыдущими схемами.
А теперь проверим явную схему - «метод установления»
Рисунок
3
Всего итераций 89.
Уменьшим шаг по сетке до h=0.59.
Т.е. X=10, а N=17.
Рисунок
4
Схема
проявила неустойчивость. Очевидно, что при данном раскладе условие (12) не выполнилось.
Исходя
из вышеизложенного, приходится при достаточно малом шаге выбирать очень мелкий шаг по времени.
Вывод
Запишем результаты в таблице
Метод
|
Количество итераций
|
Якоби
|
167
|
Гаусса-Зейделя
|
69
|
Метод верхней релаксации
|
53
|
Явный метод установления
|
89
|
Для метода Гаусса не приходится учитывать количество итераций.
1. Установлено
что явная схема имеет существенный недостаток: накладываются ограничения на
шаги и по
сетке. Чего лишены неявные схемы.
. Итерационные
методы идеально подходят для сеточной схемы «крест», метод верхней релаксации
оказался самым быстросходящимся.
. Для
метода Гаусса приходится хранить огромную матрицы в памяти ЭВМ, что при
достаточно больших N будет накладно.
Заключение
В ходе работы была изучена разностная схема «крест» для нахождения
численного решения уравнения Лапласа (эллиптическое уравнение), задача Дирихле.
Также были усвоены и закреплены навыки:
2. Программирования в среде matlab.
. Разработки численных методов для уравнений эллиптического типа.
Были созданы четыре метода реализующих разностную схему «крест» и один
метод для «метода установления».
Список использованной литературы
1.Самарский А.А. Введение в численные методы. Учебное
пособие для вузов. 3-е изд., стер. -СПб.: Издательство «Лань», 2005. - 288 с.:
ил. - (Учебники для вузов. Специальная литература).
.Бахвалов Н.С., Жидков Н. П., Кобельков Г. М.
Численные методы - 3-е изд, доп., и перераб. - М.: БИНОМ. Лаборатория знаний,
2004. - 636 с., илл.
.Агошков В.И., Дубровский П.Б., Шутяев В.П. Методы
решения задач математической физики / Под ред. Г. И. Марчука: Учеб. пособие. -
М: ФИЗМАТЛИТ, 2002.-320с.
.Мэтьюз Джон Г., Финкс Куртис Д Численные методы
использования MATLAB,3-е издание. /Пер с англ. - М. Издательский до «Вильямс»,
2001 - 720с.
.Турчак Л. И., Плотников П. В. Основы численных
методов: Учеб. пособие. - 2-е изд., перераб. И доп. - М.: ФИЗМАТЛИТ, 2005. -304
с.
.http://www.intuit.ru/department/calculate/nmdiffeq/6/
- электронная книга
.http://www.intuit.ru/department/calculate/vnmdiffeq/12
- видеолекция
Приложения
ElipYa.m
|
function ElipYa(X,N,fi,E)
%сеточный метод итерационный, процесс Якоби % X - Сторона
прямоугольника % Y - Сторона прямоугольника % N - Количество
узлов по х % K - Количество узлов по y % fi
- функция граничного условия % E - точность h=X/N;
%Равномерная сетка по пространственным переменным x=[0:h:h*N];
% Разбиваем сетку по переменной х y=[0:h:h*N];
% Разбиваем сетку по переменной у U=zeros(N+1,N+1);
% матрица значений сеточной функции F=zeros(N+1,N+1);
% матрица значений сеточной функции f(x,y)
% "дискретизация начальных условий" и начальное заполнение матрицы
значений сеточной функции U(1,:)=feval(fi,x(1),y); % Первая строка матрицы U(N+1,:)=feval(fi,x(N+1),y);
% Последняя строка матрицы U(:,1)=feval(fi,y(1),x'); % Первый столбец матрицы U(:,N+1)=feval(fi,y(N+1),x');
% Последний столбец матрицы U(2:N,2:N)=rand(N-1,N-1); % Для решения нам необходимо начальное
приближение. Осуществляем случайное начальное приближение! %Собственно
вычисление... Применяется метод Якоби M=E+0.01;
d=0; iter=0; while (M>E) % Пока не достигнем требуемой точности U0=U;
for i=2:N % произведем расчеты Xb=U(i,:); for j=2:N
v=0.25*(U0(i+1,j)+U0(i-1,j)+U0(i,j+1)+U0(i,j-1)); U(i,j)=v; end;
d=norm(U(i,:)-Xb); if(M>d)
% Если условия выполняется, то найдено новое значение остигнутой точности! M=d; end; end; iter=iter+1; end; figure; display('Всего потребовалось итераций: '); display(iter); surf(x,y,U); % Вывод поверхности!
|
Elip1.m
|
function
[x,y,U]=Elip1(X,N,fi,E) %сеточный метод, итерационный процесс Гаусса-Зейделя
% X - Сторона прямоугольника % Y - Сторона прямоугольника % N - Количество
узлов по х % K - Количество узлов по y % fxy - функция правой части % fi -
функция граничного условия % E - точность h=X/N; %Равномерная сетка по
пространственным переменным x=[0:h:h*N]; % Разбиваем сетку по переменной х
y=[0:h:h*N]; % Разбиваем сетку по переменной у U=zeros(N+1,N+1); % матрица
значений сеточной функции F=zeros(N+1,N+1); % матрица значений сеточной
функции f(x,y) % "дискретизация начальных условий" и начальное
заполнение матрицы значений сеточной функции U(1,:)=feval(fi,x(1),y); %
Первая строка матрицы U(N+1,:)=feval(fi,x(N+1),y); % Последняя строка матрицы
U(:,1)=feval(fi,y(1),x'); % Первый столбец матрицы
U(:,N+1)=feval(fi,y(N+1),x'); % Последний столбец матрицы
U(2:N,2:N)=rand(N-1,N-1); % Для решения нам необходимо начальное приближение.
Осуществляем случайное начальное приближение! %Собственно вычисление...
Применяется метод Гаусса-Зейделя M=E+0.01; d=0; iter=0; while (M>E) % Пока
не достигнем требуемой точности for i=2:N % произведем расчеты Xb=U(i,:); for
j=2:N v=0.25*(U(i+1,j)+U(i-1,j)+U(i,j+1)+U(i,j-1)); U(i,j)=v; end;
d=norm(U(i,:)-Xb); if(M>d) % Если условия выполняется, то найдено новое
значение остигнутой точности! M=d; end; end; iter=iter+1; end; %После
завершения расчетов построим поверхность и вернем значения сетки и сеточной
функции, опредеоенной на ней!. figure; display('Всего потребовалось итераций:
'); display(iter); surf(x,y,U); % Вывод поверхности!
|
ElipR.m
|
function
[x,y,U]=ElipR(X,N,fi,t,E) %сеточный метод, итерационный процесс верхней
релаксации % X - Сторона прямоугольника % Y - Сторона прямоугольника % N -
Количество узлов по х % K - Количество узлов по y % fi - функция граничного
условия % E - точность h=X/N; %Равномерная сетка по пространственным
переменным x=[0:h:h*N]; % Разбиваем сетку по переменной х y=[0:h:h*N]; %
Разбиваем сетку по переменной у U=zeros(N+1,N+1); % матрица значений сеточной
функции F=zeros(N+1,N+1); % матрица значений сеточной функции f(x,y) %
"дискретизация начальных условий" и начальное заполнение матрицы
значений сеточной функции U(1,:)=feval(fi,x(1),y); % Первая строка матрицы
U(N+1,:)=feval(fi,x(N+1),y); % Последняя строка матрицы
U(:,1)=feval(fi,y(1),x'); % Первый столбец матрицы
U(:,N+1)=feval(fi,y(N+1),x'); % Последний столбец матрицы U(2:N,2:N)=rand(N-1,N-1);
% Для решения нам необходимо начальное приближение. Осуществляем случайное
начальное приближение! %Собственно вычисление... Применяется метод верхней
релаксации M=E+0.01; d=0; iter=0; while (M>E) % Пока не достигнем
требуемой точности for i=2:N % произведем расчеты Xb=U(i,:); for j=2:N
v=(t/4)*(U(i+1,j)+U(i,j+1))-t*(1-1/t)*U(i,j)+(t/4)*(U(i-1,j)+U(i,j-1));
U(i,j)=v; end; d=norm(U(i,:)-Xb); if(M>d) % Если условия выполняется,
то найдено новое значение остигнутой точности! M=d; end; end;
iter=iter+1; end; figure; display('Всего потребовалось итераций: ');
display(iter); surf(x,y,U); % Вывод поверхности!
|
ElipGauss.m
|
function
[x,y,U]=ElipGauss(X,N,fi) % X - Сторона прямоугольника % N - Количество узлов
по х % K - Количество узлов по y % fi - функция граничного условия % E -
точность h=X/N; %Равномерная сетка по пространственным переменным
x=[0:h:h*N]; % Разбиваем сетку по переменной х y=[0:h:h*N]; % Разбиваем сетку
по переменной у U=zeros(N+1,N+1); % матрица значений сеточной функции
F=zeros(N+1,N+1); % матрица значений сеточной функции f(x,y)
Aii=-4*eye(N-1)+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1); %главная блочная
диагональ Aij=eye(N-1,N-1); % нижняя и верхняя диагонали
A=zeros((N-1)*(N-1),(N-1)*(N-1)); B=zeros((N-1)*(N-1),1); % правые части
pos11=0; pos12=0; pos21=0; pos22=0; % "дискретизация начальных
условий" и начальное заполнение матрицы значений сеточной функции
U(1,:)=feval(fi,x(1),y); % Первая строка матрицы U(N+1,:)=feval(fi,x(N+1),y);
% Последняя строка матрицы U(:,1)=feval(fi,y(1),x'); % Первый столбец матрицы
U(:,N+1)=feval(fi,y(N+1),x'); % Последний столбец матрицы
%U(2:N,2:N)=rand(N-1,N-1); % Для решения нам необходимо начальное
приближение. Осуществляем случайное начальное приближение! for i=1:N-1 %
Заполнение блочной матрицы - основной матрицы системы if(i==1) pos11=1;
pos12=N-1; else pos11=pos12+1; pos12=pos12+N-1; end; for j=1:N-1
if(j==1) pos21=1; pos22=N-1; else pos21=pos22+1; pos22=pos22+N-1; end;
if(i==j) A(pos11:pos12,pos21:pos22)=Aii; else if(i==j-1 | i==j+1)
A(pos11:pos12,pos21:pos22)=Aij; end; end; end; %Произвести заполнение B
if(pos11==1 & pos12==N-1) B(pos11:pos12,1)=-U(1,2:N)';
B(pos11,1)=B(pos11,1)-U(2,1); B(pos12,1)=B(pos12,1)-U(2,N+1); else if
(pos11==(N-1)*(N-1)-(N-2) & pos12==(N-1)*(N-1))
B(pos11:pos12,1)=-U(N+1,2:N)'; B(pos11,1)=B(pos11,1)-U(N,1);
B(pos12,1)=B(pos12,1)-U(N,N+1); else B(pos11,1)=B(pos11,1)-U(i+1,1);
B(pos12,1)=B(pos12,1)-U(i+1,N+1); end; end; end; %Собственно вычисление...
Применяется метод Гаусса с исключением нулевых %элементов d=0; iter=0; X=A\B;
for i=2:N if(i==2) pos11=1; pos12=N-1; else pos11=pos12+1; pos12=pos12+N-1;
end; U(i,2:N)=X(pos11:pos12,1)' end; figure; surf(x,y,U);
|
ElipT1.m
|
function
[x,y,U]=ElipT1(X,N,tau,fi,E) % X - Сторона квадрата % N - Количество узлов %
tau - шаг по времени % fi - функция граничного условия % E - точность
%Равномерная сетка по пространственным переменным h=X/N; x=[0:h:h*N]; %
Разбиваем сетку по переменной х y=[0:h:h*N]; % Разбиваем сетку по переменной
у U=zeros(N+1,N+1); % матрица значений сеточной функции U0=zeros(N+1, N+1); %
матрица, содержащая предыдущий слой решения по времени F=zeros(N+1,N+1); %
матрица значений сеточной функции f(x,y) % "дискретизация начальных
условий" и начальное заполнение матрицы значений сеточной функции
U(1,:)=feval(fi,x(1),y); % Первая строка матрицы U(N+1,:)=feval(fi,x(N+1),y);
% Последняя строка матрицы U(:,1)=feval(fi,y(1),x'); % Первый столбец матрицы
U(:,N+1)=feval(fi,y(N+1),x'); % Последний столбец матрицы % Заполнение
матрицы U начальными значениями функции psi(x,y) - % произвольная функция для
начального слоя for i=2:N for j=2:N U(i,j)=rand(1,1);% задаем некое случайное
начальное приближение end; end; M=E+0.1; sigma=tau/h^2; iter=0; while(M>E
& iter<=500) % процесс стремится к заданной точности U0=U; %
запоминаем предыдущее приближение for i=2:N for j=2:N
U(i,j)=sigma*(U0(i+1,j)+U0(i-1,j)+U0(i,j+1)+U0(i,j-1))+(1-4*sigma)*U0(i,j);%-tau*Fi(i,j);
end; d=norm(U(i,:)-U0(i,:)); if(M>d) M=d; end; end; iter=iter+1; end;
figure; display('Всего потребовалось итераций: '); display(iter);
surf(x,y,U); % Вывод поверхности!
|