Решение эллиптических уравнений несколькими методами

Автор работы: Пользователь скрыл имя, 21 Октября 2014 в 10:05, курсовая работа

Краткое описание

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

Содержание

Аннотация 3
Введение 4
Раздел 1. Математическое описание алгоритмов и операций 6
Раздел 2. Библиотека функций 11
Раздел 3. Тестирование 12
Вывод 16
Заключение 17
Список использованной литературы 18
Приложения 19

Скачать полностью (203.82 Кб) Сколько стоит заказать работу?

Прикрепленные файлы: 1 файл

Полный отчет по курсовой работе.doc

— 459.50 Кб (Скачать документ)

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

Заключение

В ходе работы была изучена разностная схема «крест» для нахождения численного решения уравнения Лапласа (эллиптическое уравнение), задача Дирихле.

Также были усвоены и закреплены навыки:

Использования указателей в среде matlab.
Программирования в среде matlab.
Разработки численных методов для уравнений эллиптического типа.

Были созданы четыре метода реализующих разностную схему «крест» и один метод для «метода установления».

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

Самарский А.А. Введение в численные методы. Учебное пособие для вузов. 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;

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'); % Последний столбец матрицы

%Собственно вычисление... Применяется метод Гаусса-Зейделя

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;

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'); % Последний столбец матрицы

%Собственно вычисление... Применяется метод верхней релаксации

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;

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;

%Произвести заполнение 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;

%Собственно вычисление... Применяется метод Гаусса с исключением нулевых

%элементов

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;

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;

iter=iter+1;

end;

figure;

display('Всего потребовалось итераций: ');

display(iter);

surf(x,y,U); % Вывод поверхности!

1 Агошков В.И., Дубровский П.Б., Шутяев В.П. Методы решения задач математической физики

2 Турчак Л. И., Плотников П. В. Основы численных методов

3 Бахвалов Н.С., Жидков Н. П., Кобельков Г. М. Численные методы

4 Бахвалов Н.С., Жидков Н. П., Кобельков Г. М. Численные методы

5 Самарский А.А. Введение в численные методы

Информация о работе Решение эллиптических уравнений несколькими методами