Решение эллиптических уравнений несколькими методами
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Филиал федерального государственного автономного образовательного учреждения высшего профессионального образования «Казанский (Приволжский) федеральный университет» в г. Набережные Челны
ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ
И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
КАФЕДРА ПРИКЛАДНОЙ МАТЕМАТИКИ И ИНФОРМАТИКИ
Специальность: 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 | Пользовательское соглашение
ПРОФЕССИОНАЛЬНАЯ ПОМОЩЬ СТУДЕНТАМ