Метод сеток для уравнения теплопроводности

Лабораторная работа №7

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

Решить одномерное уравнение теплопроводности методом сеток.

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

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

, j = 1, . , n , ( 10 )

где , n — число узлов cетки по x. ( 11 )

1. В первой строке введем названия параметров: n , l , c , Dx, g, m, Dt, a, b, а под ними в соответствующих ячейках — числовые значения (n=10, l=1, c=1). Для Dx вводим соответствующую формулу =B2/A2 (Dx=l/n). Исходя из условий устойчивости явной схемы, выбираем m=0,5 и выражаем Dt через m из уравнения ( 11 ) (=$D$2*$D$2*0,5).Для параметров функций, задающих краевые и начальные условия, выбираем следующие значения: g=8, a=10000, b=250.

2. В столбце А, как мы уже делали в предыдущих работах, разместим вычисление значений x, соответствующих узлам сетки. В третьей строке разместим формулы, вычисляющие значения узлов по времени. В столбце Вразместим формулы, вычисляющие начальное распределение температуры по длине стержня =EXP($E$2*A4-$E$2*A4*A4), а в четвертой и четырнадцатой строках, начиная со столбца С, -формулы, вычисляющие значения температуры на концах стержня: =EXP(-$H$2*C$3*C$3+$I$2*C$3).

3. В ячейку С5вводим формулу, реализующую конечно-разностное уравнение ( 10 ) —=B5+$F$2*(B6-2*B5+B4). Распространяем эту формулу на всю область, ограниченную краевыми и начальными условиями.

4. Результирующая таблица и построенная с использованием данных из блока A3:J14 диаграмма представлены на рис. 24 и 25.

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

Рис. 25. Графическое изображение решения одномерного уравнения теплопроводности.

Литература

1. Фигурнов В. Э. IBM PC для пользователя. 6-е изд. — М.: ИНФРА-М, 1995.

2. Шиб Й. Windows 3.1 (русская версия ) : Пер. с нем. — М. : БИНОМ, 1995.

3. Николь Н., Альбрехт Р. Электронные таблицы Excel 5.0: Практ. пособ. / Пер. с нем. — М.: ЭКОМ.,1995.

4. Наймершайм Дж. Excel 5.0 for Windows: Учебное пособие / Пер. с англ. — М.: Междунар. отношения, 1995.

5. Осейко Н. Н. Excel 5.0 для пользователя: — К.: Торгово — издательское бюро ВНV, 1994.

6. Альтхаус М. Excel. Секреты и советы / Пер. с нем. М.: БИНОМ, 1995.

7. Основы работы с Excel 5.0 : Методические указания / ИГХТА. — Иваново, 1996.

8. Численные методы : Методические указания / ИХТИ. — Иваново, 1988.

9. Методические указания и задания к практическим занятиям по вычислительной математике : Методические указания / ИХТИ. — Иваново, 1988.

10. Моделирование сложных изотермических реакций, описываемых линейными дифференциальными уравнениями : Методические указания / ИХТИ. — Иваново, 1992.

Оглавление

Лабораторная работа № 1. Основные элементы работы в EXCEL . . . 3

Лабораторная работа № 2. Построение графиков и диаграмм . . . . . . 7

Лабораторная работа № 3. Вычисление определенных интегралов . . 15

Лабораторная работа № 4. Решение систем линейных уравнений . . . 17

Лабораторная работа № 5. Обработка экспериментальных данных . 20

Лабораторная работа № 6. Решение обыкновенных дифференциаль-

ных уравнений и систем обыкновенных дифференциальных уравне-

Лабораторная работа № 7. Решение уравнений в частных производ-

ных методом сеток . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Дата добавления: 2015-02-09 ; просмотров: 18 ; Нарушение авторских прав

Решение систем линейных алгебраических уравнений

Главная > Решение

7.1 Метод сеток для решения смешанной задачи для уравнения параболического типа (уравнения теплопроводности)

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

Рассмотрим смешанную задачу для однородного уравнения теплопроводности

, k =const>0.

Задано начальное условие

и заданы краевые условия первого рода

Требуется найти функцию u (x,t) , удовлетворяющую в области D (0 x a , 0 t T) условиям (7.5) и (7.6). Физически это можно представить как стержень, на концах которого поддерживается требуемый температурный режим, заданный условиями (7.6).

Рисунок 10 – Неявная схема

При проведении замены получим , т.е. k =1. Задача решается методом сеток : строим в области D равномерную сетку с шагом h по оси x и шагом  по t (см. рисунок 10).

Приближенное значение искомой функции в точке — обозначим через . Тогда ; ; i =0,1. n ; ;

j =0,1. m ; .
Заменим производные разностными отношениями

;

.

В результате получим неявную двухслойную схему с погрешностью O (  +h 2 )

.

Используя подстановку , выразим из этой схемы u i,j-1

,

где: u 0, j =  1 ( t j ) ; u n , j =  2 ( t j ) .

Получаем разностную схему, которой аппроксимируем уравнение (7.4). Эта схема (7.7) неявная, и выглядит так, как показано на рисунке 10. При построении схемы (7.7) получается система линейных уравнений с трехдиагональной матрицой. Решив ее любым способом (в частности, методом прогонки), получаем значения функции на определенных временных слоях. Так, на нулевом временном слое используем начальное условие U i,0 =f ( x i ), т.к. j =0. Эта неявная схема более устойчива для любых значений параметра >0.

Есть и явная схема (рисунок 11), но она устойчива только при , т.е. при .

Рисунок 11 — Явная схема

7.2 Решение задачи Дирихле для уравнения Лапласа методом сеток

Рассмотрим уравнение Лапласа

.

Уравнение (7.8) описывает распространение электромагнитных волн(полей). Будем рассматривать уравнение Лапласа в прямоугольной области с краевыми условиями

; ; ; ,

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

Введем обозначения u ij = u ( x i , y j ). Накладываем на прямоугольную область сетку ; i =0,1,…, n ; ; j =0,1,…, m . Тогда , .

Частные производные аппроксимируем по формулам

и заменим уравнение Лапласа конечно-разностным уравнением

Рисунок 12 – Схема “крест”

,

где: i =1,…, n -1, j =1. m -1 (т.е. для внутренних узлов).

Погрешность замены дифференциального уравнения разностным составляет величину О(). Выразим u i , j при h =l, и заменим систему

Получаем систему (7.10) линейных алгебраических уравнений, которые можно решить любым итерационным методом (Зейделя, простых итераций и т.д.). При этом построении системы использовалась схема типа “крест”(рисунок 12). Строим последовательность итераций по методу Гаусса-Зейделя

,

где s -текущая итерация.

Условие окончания итерационного процесса

.

Условие (7.11) ненадежно и на практике используют другой критерий

где .

Схема “крест “- явная устойчивая схема ( малое изменение входных данных ведет к малому изменению выходных данных).

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

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

Задача состоит в отыскании функции u ( x , t ) при t >0, удовлетворяющей уравнению гиперболического типа

,

где: 0 x a ; 0 t

и краевым условиям

Заменим С на с t и получим уравнение

и в дальнейшем будем считать С =1.

Для построения разностной схемы решение задачи (7.12)-(7.14) построим в области сетку ; i = 0,1,…, n ; ; ; j =0,1,…, m ;  m = T .

Аппроксимируем (7.12) разностными производными второго порядка точности относительно шага

.

Полагая  =  / h перепишем (7.15), выразив U i , j +1. Таким образом получим трехслойную разностную схему

,

где: i =1,…, n ; j =1,…, m . Задаем нулевые граничные условия  1 ( t ) =0,  2 ( t ) =0. Тогда в (7.16) , для всех j .

Схема (7.16) называется трехслойной, т.к. она связывает значения искомой функции на трех временных слоях j -1, j , j +1.

Численное решение задачи состоит в вычислении приближенных значений решения u ( x , t ) в узлах при i =1,…, n ; j = 1,…, m . Алгоритм решения основан на том, что решение на каждом следующем слое ( j = 2,3. n ) можно получить пересчетом решений с двух предыдущих слоев ( j = 0,1. n — 1) по формуле (7.16). При j =0 решение известно из начального условия . Для вычисления решения на первом слое ( j = 1) положим

,

тогда , i = 1,2,…, n . Теперь для вычисления решений на следующих слоях можно использовать формулу (7.16).

Описанная схема аппроксимирует задачу (7.12)-(7.14) с точностью O (  + h ). Невысокий порядок аппроксимации по  объясняется грубостью аппроксимации по формуле (7.17).

Схема будет устойчивой, если выполнено условие .

Лабораторная работа № 1

Решение систем линейных алгебраических уравнений

Входные параметры: n—целое положительное число, равное порядку n системы; а — массив из n х n действительных чисел, содержащий матрицу коэффициентов системы (а(1) = а 11 , а(2) = a 12… а(n) = а n 1 , а(n + 1) = а 12 , . а(n х n) = а nn ); b — массив из n действительных чисел, содержащий столбец свободных членов системы (b(1) = b 1 , b(2)=b 2, …b(n)=b n ) .

Выходные параметры: b—массив из n действительных чисел (он же входной); при выходе из программы содержит решение системы b(l) = x 1 , b(2) = x 2 , … b(n) = х n ; error—признак правильности решения (код ошибки): если ks = 0, то в массиве b содержится решение системы, если error= 1, исходная система не имеет единственного решения (определитель системы равен нулю).

Перед обращением к подпрограмме SIMQ необходимо:

1) описать массивы а и b. Если система содержит n уравнений, то массив а должен содержать n 2 элементов, а массив b – n элементов;

2) присвоить значение параметру n, который равен числу
уравнений системы;

3) присвоить элементам массивов а и b значения коэффициентов системы следующим образом: a(l) = a 11 , а(2) = а 21 , а(3) = а 31 ,…а(n) = а n1 а(n+1) = а 12 , а(n+2) = а 22 … а(n x n) = а nn . b(1) = b 1 , b(2)=b 2, …b(n)=b n

4) проверить соответствие фактических параметров по типу и порядку следования формальным параметрам подпрограммы SIMQ. Параметры а и b — величины вещественного типа, n и error — целого типа.

Задание. Используя программу SIMQ, решить заданную систему трех линейных уравнений. Схема алгоритма приведена на рисунке 13.

Порядок выполнения лабораторной работы:

1. Составить головную программу, содержащую обращение к SIMQ и печать результатов;

2. Произвести вычисления на ЭВМ.

Пример. Решить систему уравнений

Рисунок 13 – Схема алгоритма метода Гаусса

PROCEDURE SIMQ(Nn:Integer;Var Aa:TMatr;Var Bb:TVector;Var Ks:Integer);

Var Max,U,V : Real; I,J,K1,L : Integer;

For I:=1 To Nn Do Aa[i,Nn+1]:=Bb[i];

For I:=1 To Nn Do

For L:=I+1 To Nn Do If (Abs(Aa[l,i])>Max) Then

For J:=I To Nn+1 Do

Begin U:=Aa[i,j]; Aa[i,j]:=Aa[k1,j]; Aa[k1,j]:=U;

For J:=I To Nn+1 Do Aa[i,j]:=Aa[i,j]/V;

V:=Aa[l,i]; For J:=I+1 To Nn+1 Do Aa[l,j]:=Aa[l,j]-Aa[i,j]*V;

For I:=Nn-1 Downto 1 Do

For J:=I+1 To Nn Do Bb[i]:=Bb[i]-Aa[i,j]*Bb[j];

Вычисления по программе привели к следующим результатам:

X(1)= .100000E+01 Х(2)= .200000Е+01 Х(3)= .З00000Е + 01

признак выхода 0

Варианты заданий для решения систем линейных алгебраических уравнений методом Гаусса приведены в таблице 1.

Метод квадратных корней Холецкого

Входные параметры: n—целое положительное число, равное порядку n системы; а — массив из n х n действительных чисел, содержащий матрицу коэффициентов системы (а(1) = а 11 , а(2) = a 12… а(n) = а n 1 , а(n + 1) = а 12 , . а(n х n) = а nn ); b — массив из n действительных чисел, содержащий столбец свободных членов системы (b(1) = b 1 , b(2)=b 2, …b(n)=b n ) .

Выходные параметры: b—массив из n действительных чисел (он же входной); при выходе из программы содержит решение системы b(l) = x 1 , b(2) = x 2 , … b(n) = х n ; p—количество операций.

Схема алгоритма приведена на рисунке 14.

Пример. Решить систему уравнений

Procedure Holets(n:integer;a:TMatr;b:TVector;var x:TVector;var p:integer);

Var i,j,k:integer; a11:real;

For i:=1 To n Do Begin

If i<>1 Then Begin

If a[i,i]=0 Then Begin

p:=0; error:=2; MessageDlg(‘. ‘,mtError,[mbOk],0);

For j:=1 To i Do Begin

For k:=1 To j-1 Do Begin

For i:=1 To n Do Begin

For j:=1 To i-1 Do b[i]:=b[i]-a[i,j]*b[j];

If a[i,i]=0 Then Begin

p:=0; error:=2; MessageDlg(‘. ‘,mtError,[mbOk],0);

For i:=n DownTo 1 Do Begin

For j:=n DownTo i+1 Do b[i]:=b[i]-a[i,j]*b[j];

Вычисления по программе привели к следующим результатам:

X(1)= .100000E+01 Х(2)= .200000Е+01 Х(3)= .З00000Е + 01

Рисунок 14 — Схема алгоритма метода Холецкого

Тема лабораторной работы №1 для контроля знаний проиллюстрирована контрольно – обучающей программой.

Лекция 1: Исследование разностных схем для эволюционных уравнений на устойчивость и сходимость

1.1 Постановка некоторых задач для уравнений математической физики

Напомним некоторые корректные постановки задач для уравнений в частных производных , которые будут встречаться в дальнейшем.

Задача Коши для уравнения теплопроводности.

Найти функцию u(t , x) в области , удовлетворяющую уравнению

и начальным данным u(0, x) = u₀(x) , где u₀(x) — заданная функция .

Смешанная задача для уравнения теплопроводности.

Найти функцию u(t , x) в области , удовлетворяющую уравнению

начальным данным u(0, x) = u₀(x) и краевым условиям, записанным в общей форме

Смешанная задача для уравнения переноса.

Найти функцию u(t , x) в области , удовлетворяющую уравнению (для определенности положим c > 0 ):

начальным данным Коши

и левому краевому условию

Смешанная задача для волнового уравнения.

Найти функцию u(t , x) в области , удовлетворяющую уравнению

с начальными данными

и краевыми условиями .

Эллиптическая краевая задача (уравнение Пуассона).

Найти функцию u(x , y) в области , удовлетворяющую уравнению

и краевым условиям , , .

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

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

Для решения одномерной смешанной задачи для уравнений в частных производных параболического типа область определения искомой функции покрывается расчетной сеткой с узлами в точкахn, x_m> , n = 0, . , N , m = 0, . , M , , x_m = mh , , h = X/M , где — шаги сетки по времени и пространству соответственно. Приближенным решением задачи назовем сеточную функцию . Верхний индекс в такой форме записи сеточной функции традиционно указывает на номер слоя по времени, нижний (нижние) — на номер узла сетки по пространственной координате (рис. 1.1).

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

Явная разностная схема для приближенного решения уравнения теплопроводности во внутренних узлах сетки (не принадлежащим границе сеточной области) имеет вид

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

Эта схема аппроксимирует дифференциальное уравнение во внутренних точках (узлах) области интегрирования, т.е. при n = 1 , . , N — 1 , m = 1 , . , M — 1 . Проведем аппроксимацию начальных данных и краевых условий:

для определенности положим .

Расчет ведется по рекуррентной формуле на каждом временном слое от n = 1 до n = N от m = 1 до m = M — 1 во внутренних узлах; слой n = 0(t = t₀) соответствует начальным данным, лучи m = 0(x = x₀) и m = M(x = x_M) — левому и правому краевым условиям.

Запишем явную схему в виде

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

Для завершения расчета слоя t = t_{n + 1} необходимо вычислить и , для чего разрешаем левое и правое краевые условия относительно этих величин:

где и уже вычислены ранее. Реализация одного шага по времени занимает O(M) арифметических операций , всех слоев — O(NM) операций.

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

Простейшая неявная разностная схема имеет вид (для простоты положим a = 1 )