Решение эллиптических уравнений несколькими методами
2
1
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Филиал федерального государственного автономного образовательного учреждения высшего профессионального образования «Казанский (Приволжский) федеральный университет» в г. Набережные Челны
ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ
И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
КАФЕДРА ПРИКЛАДНОЙ МАТЕМАТИКИ И ИНФОРМАТИКИ
Специальность: 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. Установлено что явная схема имеет существенный недостаток: накладываются ограничения на шаги и по сетке. Чего лишены неявные схемы.
2. Итерационные методы идеально подходят для сеточной схемы «крест», метод верхней релаксации оказался самым быстросходящимся.
3. Для метода Гаусса приходится хранить огромную матрицы в памяти ЭВМ, что при достаточно больших N будет накладно.
Заключение
В ходе работы была изучена разностная схема «крест» для нахождения численного решения уравнения Лапласа (эллиптическое уравнение), задача Дирихле.
Также были усвоены и закреплены навыки:
1. Использования указателей в среде matlab.
2. Программирования в среде matlab.
3. Разработки численных методов для уравнений эллиптического типа.
Были созданы четыре метода реализующих разностную схему «крест» и один метод для «метода установления».
Список использованной литературы
1.Самарский А.А. Введение в численные методы. Учебное пособие для вузов. 3-е изд., стер. -СПб.: Издательство «Лань», 2005. — 288 с.: ил. — (Учебники для вузов. Специальная литература).
2.Бахвалов Н.С., Жидков Н. П., Кобельков Г. М. Численные методы — 3-е изд, доп., и перераб. — М.: БИНОМ. Лаборатория знаний, 2004. — 636 с., илл.
3.Агошков В.И., Дубровский П.Б., Шутяев В.П. Методы решения задач математической физики / Под ред. Г. И. Марчука: Учеб. пособие. — М: ФИЗМАТЛИТ, 2002.-320с.
4.Мэтьюз Джон Г., Финкс Куртис Д Численные методы использования MATLAB,3-е издание. /Пер с англ. — М. Издательский до «Вильямс», 2001 — 720с.
5.Турчак Л. И., Плотников П. В. Основы численных методов: Учеб. пособие. — 2-е изд., перераб. И доп. — М.: ФИЗМАТЛИТ, 2005. -304 с.
6.http://www.intuit.ru/department/calculate/nmdiffeq/6/ — электронная книга
7.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=AB; 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); % Вывод поверхности! |
Нужна похожая работа?
Оставь заявку на бесплатный расчёт