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

Лабораторная работа №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 для контроля знаний проиллюстрирована контрольно – обучающей программой.

Уравнение теплопроводности

Разностные методы решения дифференциальных

Уравнений в частных производных.

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

Уравнение теплопроводности.

Рассмотрим задачу нахождения непрерывной на замкнутом прямоугольнике D=<0£x£l, 0£y£Y> функции u(x, y), удовлетворяющей уравнению теплопроводности

, (*)

если заданы начальные условия:

,

где f(x, y), j(x), p(y), s(y) –заданные, функции (n раз дифференц.) такие что

.

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

Идея метода заключается в следующем.

Разобьём отрезки [0, l] оси х и [0; T] оси у соответственно на n и n₁ равных частей и введём обозначения: . Через точки деления проведём прямые , параллельные соответствующим осям. В результате область D разобьётся на прямоугольники с вершинами (х_i, y_i), где x_i=(i-1)h, i=1,n+1, у_i=(j-1)×t, i=1,n₁+1.

Множество вершин прямоугольников называется сеткой, а отдельные вершины – узлами сетки. Узлы, имеющие одинаковый индекс j, образуют j слой. Числа h и t называют шагами сетки соответственно по переменным х и у.

По определению частная производная равна

Если рассматривать функцию только в узлах сетки, то частную производную можно записать в форме

где узел соответствует точке .

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

Такое выражение называется левой конечной разностью. Можно получить центральную конечную разность, найдя среднее этих выражений.

Теперь получим выражения для вторых производных.

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

Заменяя производные, входящие в уравнение (*) разностными отношениями, получим конечно-разностные уравнения

(**)

. (***)

Эти уравнения аппроксимируют исходное дифференциальное уравнение в узле сетки (х_i, у_i) с погрешностью порядка О(h 2 +t).

Для получения первого уравнения была использована конфигурация узлов (1), а для второго (2).

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

Тогда во внутренних точках сетки решение можно искать в явном виде по схеме (из уравнения (**)