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

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

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

Якутск - 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 Результаты

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



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

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

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

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

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

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

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



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

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

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

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

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



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

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

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

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

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

Граничные условия можно разделить на три основных рода: (10 Ozisik Boundary value problems of heat conduction)

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

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

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

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

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

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



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

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

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

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

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

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



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

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

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

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

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

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

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

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

Для замены области непрерывного изменения аргумента дискретным аналогом его изменения необходимо выбрать в данной области конечное множество точек - сетку. Тогда приближенное решение искать нужно будет только в узлах этой сетки. Функцию, определенную в узлах сетки будем называть сеточной функцией.( Самарский А. А. Введение в теорию разностных схем. - М. : Наука, 1971.)

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

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

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

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

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

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

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

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



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

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

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

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

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



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

2. Неявная схема: значение узла на новом слое зависит и от соседних узлов на новом слое, и от значения на предыдущем слое(рисунок 3). Данная схема всегда является устойчивой, (1 Ozisik Finite difference methods in heat transfer) поэтому использовать ее.

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

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

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



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

Подставим получившиеся выражения в уравнение теплопроводности и получим:

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

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

Введем параметр:

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



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



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

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

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

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

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

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

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

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

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

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

Для соблюдения условия (3.1.2) необходимо и достаточно, чтобы все главные миноры матрицы были положительны. (Gilbert G. T. Positive definite matrices and Sylvester's criterion //American Mathematical Monthly. - 1991. - С. 44-46.)

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

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

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

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

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



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

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

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

где

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

-длина шага

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

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

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

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

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

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

Теперь путем подстановок и преобразований получим длину шага (Shewchuk J. R. An introduction to the conjugate gradient method without the agonizing pain. - 1994.)



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

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

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



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

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

В отличие от метода наискорейшего спуска, где спуск часто производится по тем же направлениям, которые использовались ранее, в методе сопряженных градиентов спуск производится по сопряженным направлениям без их повторного использования. Данный метод позволяет минимизировать квадратичную функцию за шагов (если не учитывать ошибки округления).( Nocedal J., Wright S. J. Conjugate gradient methods. - Springer New York, 2006. - С. 101-134.)

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

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

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

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



:

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



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

4.1 Организация

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2 Результаты

Intel Core i5-3470 @ 3.20 GHz

Метод	Кол-во итераций	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

Количество узлов сетки: 64 x 64



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

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

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

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

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

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

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

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

8. Strikwerda J. C. Finite difference schemes and partial differential equations. - Siam, 2004.

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

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

...