Численное решение уравнения теплопроводности методом конечных разностей

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

Северо-Восточный федеральный университет имени М.К. Аммосова

Институт математики и информатики

Кафедра информационных технологий

Курсовая работа на тему

Численное решение уравнения теплопроводности методом конечных разностей

Студента 3 курса ИМИ

группы ФИИТ-11

Климова Сергея Сергеевича

Якутск 2014



ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ

1. УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ

1.1 Вывод уравнения теплопроводности

1.2 Условия однозначности

2. МЕТОД КОНЕЧНЫХ РАЗНОСТЕЙ

2.1 Формулировка метода из рядов Тейлора

2.2 Построение разностной схемы

3. РЕШЕНИЕ СЛАУ

3.1 Градиентные методы

3.1.1 Метод наискорейшего спуска

3.1.2 Метод сопряженных градиентов

4. РЕАЛИЗАЦИЯ

4.1 Функциональные требования

4.2 Параллелизация вычислений

4.3 Результаты

ЗАКЛЮЧЕНИЕ

СПИСОК ЛИТЕРАТУРЫ



ВВЕДЕНИЕ

Основным преимуществом численных методов является возможность замены дорогостоящего или трудновыполнимого физического эксперимента, а также возможность моделирования процессов, не поддающихся аналитическому решению. Необходимость в численном моделировании процесса теплопроводности возникает во многих отраслях современной техники. Одним из наиболее простых численных методов решения уравнения теплопроводности является метод конечных разностей.

Объект исследования - начально-краевая задача для нестационарного уравнения теплопроводности. Предмет исследования - решение начально-краевой задачи для нестационарного уравнения теплопроводности методом конечных разностей.

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

Для достижения указанной цели необходимо решить следующие задачи:

1. Изучение теории построения разностных схем.

2. Изучение методов решения СЛАУ.

3. Параллелизация вычислений.



1. УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ

Уравнение теплопроводности описывает распространение тепла в заданной области пространства в зависимости от времени.

Является параболическим дифференциальным уравнением в частных производных. В двумерном случае уравнение имеет вид:

		(1.1)

1.1 Вывод уравнения теплопроводности

Представим однородное тело и вычленим из него элементарный объем со сторонами , , (рисунок 1).

Рисунок 1. Контрольный объем в прямоугольной системе координат

Входящие потоки тепла, расположенные перпендикулярно к поверхностям обозначим как , , . Потоки на противоположных поверхностях выразим из рядов Тейлора:

		(1.1.1)

Внутри тела так же могут быть внутренние источники тепла, если и стоки, если :

Изменение внутренней энергии:

		(1.1.2)
		(1.1.3)
		(1.1.4)

Подставим уравнения (1.1.1), (1.1.3) и (1.1.4) в уравнение (1.1.2) и получим:

		(1.1.5)

Подставим уравнения (1.1.1) в получившееся уравнение (1.1.5):

		(1.1.6)

Потоки тепла выразим из закона Фурье:

		(1.1.7)

Подставив их в уравнение (1.1.6), получим уравнение теплопроводности в общем виде для трехмерного пространства:

Введем коэффициент температуропроводности:

и опустим внутренние источники тепла. Получим уравнение теплопроводности в трехмерном пространстве без внутренних источников тепла:



1.2 Условия однозначности

Уравнение (1.1) описывает процесс в общем виде. Для ее применения к конкретной задаче необходимы дополнительные условия, называемые условиями однозначности. Данные условия включают в себя геометрические(форма и размеры тела), физические (физические свойства тела), временные(начальное распределение температуры) и граничные условия(описывают процесс теплообмена с окружающей средой).

Граничные условия можно разделить на три основных рода [5]:

1. Граничные условия Дирихле: задано значение функции на границе.

В случае задачи теплопроводности задают значения температуры на поверхности тела.

2. Граничные условия Неймана: задана нормальная производная функции на границе.

Задают плотность теплового потока на поверхности тела.

3. Граничные условия Робена: задана линейная комбинация значения функции и ее производной на границе.

Описывают теплообмен между поверхностью тела и окружающей средой по закону Ньютона-Рихмана.

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



2. МЕТОД КОНЕЧНЫХ РАЗНОСТЕЙ

2.1 Формулировка метода из рядов Тейлора

Суть метода конечных разностей может быть выражена через определение производной функции u в точке :

		(2.1.1)

Данная формула может стать подходящей заменой производной в случае, если функция u является непрерывной и является достаточно малой, но конечной величиной.

Представим разложение в ряд Тейлора функции u в окрестности точки в правом и левом направлениях:

		(2.1.2)
		(2.1.3)

Данные выражения формируют основу для конечно-разностной аппроксимации производной первого порядка в окрестности точки После перестановок в формулах (2.1.2) и (2.1.3) получим правую и левую конечно-разностные аппроксимации производной первого порядка соответственно:

		(2.1.4)
		(2.1.5)

где обозначает ошибку аппроксимации, показывающую разность между производной и её конечно-разностным представлением.

Для правой конечной разности ошибка имеет вид:

		(2.1.6)

Получим центральную конечную разность путем вычитания (2.1.5) из (2.1.4):

		(2.1.7)

Полученная аппроксимация имеет ошибку:



2.2 Построение разностной схемы

Для построения разностной схемы необходимо:

1. Совершить переход из области непрерывного изменения аргумента в её дискретный аналог.

2. Заменить производные их конечно-разностными аналогами.

Для замены области непрерывного изменения аргумента дискретным аналогом его изменения необходимо выбрать в данной области конечное множество точек - сетку. Тогда приближенное решение искать нужно будет только в узлах этой сетки. Функцию, определенную в узлах сетки будем называть сеточной функцией [1].

Пусть функция зависит от одной переменной . Разобьем отрезок на равных частей. Тогда шаг сетки . Узлами сетки будут точки деления . Получили сетку, состоящую из множества всех её узлов.

Тогда аппроксимацию производной первого порядка в узле можно произвести с помощью:

Правой конечной разности:

		(2.2.1)

Левой конечной разности:

		(2.2.2)

Центральной конечной разности:



Произведем аппроксимацию производной второго порядка

Для этого представим производную второго порядка функции как производную первого порядка некоторой функции :

Аппроксимируем производную функции в точке правой конечной разностью (2.2.1):

Аппроксимируем производные функции левой конечной разностью (2.2.2):

		(2.2.3)



В зависимости от выбора способа аппроксимации производной по времени можно получить два основных вида разностных схем:

1. Явная схема: значение узла на новом временном слое зависит только от значений узлов на предыдущем слое, то есть значение может быть вычислено явно из предыдущего слоя(рисунок 2). Данная разностная схема является устойчивой только при следующих условиях [6]:

Одномерный случай:

Двумерный случай (:

Рисунок 2. Явная разностная схема

2. Неявная схема: значение узла на новом слое зависит и от соседних узлов на новом слое, и от значения на предыдущем слое (рисунок 3). Данная схема всегда является устойчивой [6].



Рисунок 3. Неявная разностная схема

Аппроксимируем частную производную функции u по переменной t в точке используя правую конечную разность (2.2.1):

		(2.2.4)

Аппроксимируем частную производную второго порядка функции по переменной в точке с помощью (2.2.3):

		(2.2.5)

Аналогично аппроксимируем частную производную второго порядка функции по переменной в точке и получаем:

		(2.2.6)

Подставим получившиеся выражения (2.2.4), (2.2.5) и (2.2.6) в уравнение теплопроводности (1.1) и получим:



Элементы с n-ым шагом по времени перенесем вправо, элементы с (n+1)-ым шагом по времени перенесем влево:

Для удобства будем считать, что и введем параметр:

		(2.2.7)

При этом для учета граничных условий Дирихле, значения внешних узлов, граничащих с внутренними узлами необходимо перенести в правую часть.

Запишем полученную разностную схему (2.2.7) в виде СЛАУ:

		(2.2.8)
	(2.2.9)

Координаты узлов из двумерных переведем в одномерные в лексикографическом порядке сверху вниз, слева направо:

		(2.2.10)
		(2.2.11)



3. РЕШЕНИЕ СЛАУ

Для решения СЛАУ с большим количеством переменных и разряженной симметричной положительно-определенной матрицей одним из наиболее популярных методов является метод сопряженных градиентов[7; с. 2]. Перед методом сопряженных градиентов будет рассмотрен более простой, но медленный метод наискорейшего спуска.

3.1 Градиентные методы

Чтобы найти решение СЛАУ(2.2.8) с положительно-определенной матрицей (2.2.9), достаточно найти набор значений , при которых соответствующая квадратичная форма достигает своего минимального значения.

Ввиду положительной определенности матрицы A, график квадратичной формы имеет вид параболоида с ветвями образующих парабол вверх (рисунок 4).

Рисунок 4. График квадратичной формы с положительно-определенной матрицей

Квадратичная форма имеет вид:



		(3.1.1)

где - вектор правой части СЛАУ

- вектор неизвестных

- симметричная, положительно-определенная матрица

Симметричная матрица является положительно-определенной, если

		(3.1.2)

Для соблюдения условия (3.1.2) необходимо и достаточно, чтобы все главные миноры матрицы были положительны [2].

Для того чтобы показать что минимизирующий также служит решением СЛАУ возьмем градиент квадратичной формы (3.1.1) и приравняем его нулю:

Ввиду симметричности матрицы A ( получим:

В градиентных методах направления спуска выбираются исходя от градиента функции в текущей точке. Градиентные методы отличаются друг от друга способом выбора направления спуска и длиной шага. В данной работе будут рассмотрены два градиентных метода: метод наискорейшего спуска и метод сопряженных градиентов.

3.1.1 Метод наискорейшего спуска

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

Рисунок 5. Контурный график квадратичной формы с изображенными линиями спуска [7]

Направление спусков выбирается в направлении наискорейшего убывания функции в данной точке (то есть по направлению антиградиента)

Каждая следующая точка находится как:

где

- точка, из которой выполняется спуск

-длина шага

- антиградиент в точке

		(3.1.1.1)

Из уравнения (3.1.1.1) видно, что вектор также является невязкой. Его можно использовать в условии остановки последовательности спусков, когда относительная невязка становится меньше заданного числа.

Длина шага выбирается из учета того, что спускаться нужно только до тех пор, пока функция на линии, соответствующей выбранному направлению (3.1.1.1) убывает. Представим параболу, образованную пересечением параболоида и вертикальной плоскости (рисунок 6).

Рисунок 6. Пересечение графика квадратичной формы с вертикальной плоскостью, образующее параболу g(x) на которой необходимо найти минимум

Нужно найти такую , чтобы достигала на ней своего минимума.

		(3.1.1.2)

Из формулы (3.1.1.2) следует, что направление каждого последующего спуска должно быть ортогонально направлению предыдущего спуска.

Путем подстановок и преобразований получим длину шага [8; с. 6]:

В итоге, метод наискорейшего спуска выглядит следующим образом:

		(3.1.1.3)
		(3.1.1.4)

Приведенный выше алгоритм требует два умножения матрицы на вектор для каждой итерации. Чтобы избавиться от одного матрично-векторного произведения (3.1.1.3) в уравнении (3.1.1.4) обе части умножим на ( и прибавим :

Получили следующий алгоритм:

Вычислительная сложность данного алгоритма в основном определена векторно-матричным произведением. Умножение вектора на матрицу требует операций, где - количество ненулевых элементов матрицы. Предположим, что мы хотим произвести достаточно итераций, чтобы уменьшить норму ошибки

в раз:

		(3.1.1.5)

Тогда максимальное количество итераций необходимое для достижения условия (3.1.1.5) будет [8]:

где - число обусловленности.

Временная сложность метода наискорейшего спуска: . Пространственная сложность: Недостатком метода наискорейшего спуска является медленная сходимость, если квадратичная функция имеет «овраг»: спуск может пойти мелкими «зигзагами» вдоль оврага, потратив на это большое количество шагов.

3.1.2 Метод сопряженных градиентов

В отличие от метода наискорейшего спуска, где спуск часто производится по тем же направлениям, которые использовались ранее, в методе сопряженных градиентов спуск производится по сопряженным направлениям без их повторного использования. Данный метод позволяет минимизировать квадратичную функцию за шагов (если не учитывать ошибки округления) [4].

Множество векторов - направлений спуска называются сопряженными по отношению к матрице (или - ортогональными), если

Особенность метода сопряженных градиентов состоит в способе генерации сопряженных векторов: для вычисления вектора нужно знать только предыдущий вектор .

Каждое направление спуска вычисляется как линейная комбинация направления наискорейшего спуска и направления спуска на предыдущем шаге:

где находится из условия сопряженности :



В итоге получили алгоритм метода сопряженных градиентов:

Аналогично методу наискорейшего спуска, вычислительная сложность определяется умножением вектора на матрицу. Временная сложность метода сопряженных градиентов: . Пространственная сложность:

Максимальное количество итераций необходимое для достижения условия (3.1.1.5)[8]:



4. РЕАЛИЗАЦИЯ

Использование в алгоритмах методов наискорейшего спуска и сопряженных градиентов матричных операций, которые легко поддаются распараллеливанию, позволяет добиться ускорения вычислений путем их распараллеливания. Для организации передачи данных в параллельных вычислениях используются два основных подхода: использование разделяемой памяти(shared memory) и передача сообщений. В данной работе используется модель передачи сообщений, ввиду своей универсальности и умению выражать параллельные алгоритмы [3]. Спецификация MPI (Message Passing Interface) является одной из наиболее распространенных. Для реализации параллельных алгоритмов использовался язык программирования C и реализация интерфейса MPI MPICH2.

4.1 Функциональные требования

Программа должна обеспечивать возможность выполнения следующих функций:

-Главный процесс должен считывать из файла коэффициент температуропроводности, длину прямоугольной пластины, ширину прямоугольной области, длину по времени, количество узлов по длине, количество узлов по ширине, точность вычислений, максимальное количество итераций, температуру на верхней границе области, температуру на нижней границе области, температуру на левой границе области, температуру на правой границе области, совокупность температур на узлах в начальный момент времени.

-При допустимых входных данных программа должна уметь распределённым образом решить краевую задачу для нестационарного уравнения теплопроводности методом конечных разностей на прямоугольной области.

-Главный процесс должен выводить совокупность температур на узлах по строкам в текстовый файл.

4.2 Параллелизация вычислений

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

Вместо хранения матрицы A будем использовать функцию, возвращающую её элемент:

double A(int i, int j, int xys, double s)

{

if (i == j)

return 1 + 4 * s;

if (i + 1 == j || i - 1 == j || i + xys == j || i - xys == j)

return -s;

return 0;

}

Параллелизация скалярного произведения (рисунок 7):

С помощью MPI_Scatter разбиваем оба вектора на куски и рассылаем их соответствующим процессам.

Каждый процесс вычисляет свою часть скалярного произведения.

С помощью MPI_Reduce выполняем редукцию суммированием, получая результат скалярного произведения в нулевом процессе.

Рисунок 7. Параллелизация скалярного произведения

I. Параллелизация умножения матрицы на вектор (рисунок 8):

1. Так как значения элементов матрицы мы берем из функции, которая доступна всем процессам, то рассылать нужно только вектор с помощью MPI_Broadcast.

2. Каждый процесс вычисляет свой кусок результирующего вектора

3. Собираем куски вектора в один на нулевом процессе с помощью MPI_Gather.

Рисунок 8. Параллелизация умножения матрицы на вектор

II. Параллелизация суммирования векторов (рисунок 9):

1. С помощью MPI_Scatter разбиваем оба вектора на куски и рассылаем их соответствующим процессам.

2. Каждый процесс вычисляет свой кусок результирующего вектора

3. Собираем куски вектора в один на нулевом процессе с помощью MPI_Gather.

Рисунок 9. Параллелизация суммирования векторов

4.3 Результаты

Исследуем работу программы на следующей задаче:

Найти распределение температур на области в момент времени Для тестирования использовалась система на основе Intel Core i5-3470 @ 3.20 ГГц. Количество узлов сетки: 64 x 64.

Время, затраченное на решение СЛАУ и ускорение, соответствующее количеству процессов представлено в таблице 1:



Таблица 1. Время и ускорение решения СЛАУ

Метод	Кол-во итераций	1 процесс	2 процесса	4 процесса
		Время (с)	Время (с)	Ускорение	Время (с)	Ускорение
SD	2563	81.263	43.867	1.852	27.051	3.004
CG	85	2.805	1.560	1.798	0.999	2.807

Время, затраченное на одну итерацию, представлено в таблице 2:

Таблица 2. Время (в секундах), затраченное на одну итерацию:

Метод	1 процесс	2 процесса	4 процесса
SD	0.031706	0.017115	0.010554
CG	0.033	0.018352	0.011752

График ускорения от количества процессов представлен на рисунке 11:

Рисунок 11. Зависимость ускорения решения СЛАУ от количества процессов

Результаты работы программы для решения описанной задачи в различные моменты времени представлены на рисунке 10:



t=0	t=8
t=32	t=256

Рисунок 10. Распределение температуры в различные моменты времени



ЗАКЛЮЧЕНИЕ

В данной курсовой работе для достижения поставленной цели в ходе изучения теории построения разностных схем был выбран полностью неявный метод ввиду свободы выбора шага по времени. Был кратко рассмотрен метод сопряженных градиентов как наиболее популярный при решении разряженных симметричных положительно-определенных СЛАУ и наряду с ним был рассмотрен простой для понимания, но медленный метод наискорейшего спуска. Для реализации параллельных вычислений было решено использовать интерфейс MPI.

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



СПИСОК ЛИТЕРАТУРЫ

1. Самарский А. А. Введение в теорию разностных схем. - М.: Наука, 1971.

2. Gilbert G. T. Positive definite matrices and Sylvester's criterion //American Mathematical Monthly. - 1991. - С. 44-46.

3. Gropp W., Lusk E., Skjellum A. Using MPI: portable parallel programming with the message-passing interface. - MIT press, 1999. - Т. 1.

4. Nocedal J., Wright S. J. Conjugate gradient methods. - Springer New York, 2006. - С. 101-134.

5. Ozisik M. N. Boundary value problems of heat conduction. - Courier Dover Publications, 2013.

6. Ozisik N. Finite difference methods in heat transfer. - CRC press, 1994.

7. Refsnжs R. H. A brief introduction to the conjugate gradient method - 2009.

8. Shewchuk J. R. An introduction to the conjugate gradient method without the agonizing pain. - 1994.

Размещено на Allbest.ru

...