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

 

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ

Филиал федерального государственного автономного образовательного учреждения высшего профессионального образования «Казанский (Приволжский) федеральный университет» в г. Набережные Челны

ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ
И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
КАФЕДРА ПРИКЛАДНОЙ МАТЕМАТИКИ И ИНФОРМАТИКИ

Специальность: 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.Якоби,


(6)


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 будет накладно.


Заключение


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

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

1.Использования указателей в среде matlab.

2.Программирования в среде matlab.

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

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


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


1.Самарский А.А. Введение в численные методы. Учебное пособие для вузов. 3-е изд., стер. -СПб.: Издательство «Лань», 2005. - 288 с.: ил. - (Учебники для вузов. Специальная литература).

.Бахвалов Н.С., Жидков Н. П., Кобельков Г. М. Численные методы - 3-е изд, доп., и перераб. - М.: БИНОМ. Лаборатория знаний, 2004. - 636 с., илл.

.Агошков В.И., Дубровский П.Б., Шутяев В.П. Методы решения задач математической физики / Под ред. Г. И. Марчука: Учеб. пособие. - М: ФИЗМАТЛИТ, 2002.-320с.

.Мэтьюз Джон Г., Финкс Куртис Д Численные методы использования MATLAB,3-е издание. /Пер с англ. - М. Издательский до «Вильямс», 2001 - 720с.

.Турчак Л. И., Плотников П. В. Основы численных методов: Учеб. пособие. - 2-е изд., перераб. И доп. - М.: ФИЗМАТЛИТ, 2005. -304 с.

.#"center">Приложения


ElipYa.mfunction 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.mfunction [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.mfunction [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.mfunction [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.mfunction [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); % Вывод поверхности!


ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ Филиал федерального государственного автономного образовательного учреждения высшего профессионального образования «

Больше работ по теме:

КОНТАКТНЫЙ EMAIL: [email protected]

Скачать реферат © 2017 | Пользовательское соглашение

Скачать      Реферат

ПРОФЕССИОНАЛЬНАЯ ПОМОЩЬ СТУДЕНТАМ