Метод матричной прогонки
Перепишем разностную схему (10.22) так:
Um+1,n − 2Umn + Um−1,n + α(Um,n+1 − 2Umn + Um,n−1) = h2f(xm, yn),
m = 1, 2, . . . , M − 1, n = 1, 2, . . . , N − 1 U0n = φ(0, yn), UMn = φ(a, yn) Um0 = φ(xm, 0), UmN = φ(xm, b)
h2
l2 = α > 0.
Пусть M > N. Введем обозначение
|
Um,1 |
|
|
|
|
. . . |
|
Um = Um,N. . .−1 |
|
, m = 0, . . . , M
(10.35)
(10.36)
Предположим, что в формулах (10.35) n = 1, . . . , N − 1 и, используя (10.36), запишем систему уравнений (10.35) в векторной форме:
U0 = φ0, UM = φa. |
|
− |
|
(10.37) |
Um+1 + AUm + Um−1 |
= fm, m = 1, . . . , M |
|
1 |
|
где A — трехдиагональная матрица порядка N − 1 вида
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
|
−(2 a |
|
(2 + 2a) a . . . 0 |
|
|
|
0 |
|
|
|
|
0 |
|
|
||||||||||||||||
A = |
|
+ 2a) |
|
|
|
a |
|
|
|
|
0 . . . 0 |
|
|
|
0 |
|
|
|
|
0 |
|
, |
||||||||
. . . |
− . . . |
|
|
|
. . . . . . . . . |
|
|
|
. . . |
|
|
|
. . . |
|||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0 |
|
|
|
0 |
|
|
|
|
0 . . . a |
|
|
(2 + 2a) |
|
a |
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
|
|
|
|
|
|
|
|
|
|
|
0 |
|
|
|
0 |
|
|
|
|
0 . . . |
|
0 |
|
a |
|
|
|
|
(2 + 2a) |
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
h f(xm, y1) − aφ(xm, 0) |
|
|
|
|
|
|
|
|
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
h2f(xm, y2) |
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
f |
m |
= |
|
|
|
|
|
|
|
. . . |
|
|
|
|
|
|
|
|
|
, |
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
h2f(x |
m |
, y |
N−2 |
) |
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
h2f(x |
m |
, y |
N−1 |
) |
|
|
aφ(x |
m |
, b) |
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
φ(0, y |
) |
|
|
|
|
|
|
|
|
|
|
φ |
a, y |
) |
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
1 |
|
|
, φ = |
|
|
|
( |
|
|
1 |
|
|
. |
|
|
|
|||||||
|
|
φ = |
|
φ(0, y2) |
|
φ(a, y2) |
|
|
|
|
|
|||||||||||||||||||
|
|
0 |
|
|
|
. . . |
|
|
|
|
|
|
a |
|
|
|
|
|
|
. . . |
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
φ(0, y |
N−1 |
) |
|
|
|
|
|
|
|
φ(a, y |
N |
−1 |
) |
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Задача (10.37) аналогична задаче, которая решалась нами методом прогонки. Ее отличие лишь в том, что она имеет векторную форму. Приведем здесь алгоритм решения задачи (10.37), который называется
матричной прогонкой.
1.По формуле Rm+1 = −(A + Rm)−1, m = 1, . . . , M − 1, пологая R1 = 0, вычисляем матрицы Rm = (Rij(m)), m = 1, 2, . . . , M. Порядок этих матриц (N − 1)(N − 1). После этого положим вектор S1 = φ0, а затем по формуле Sm+1 = Rm+1(Sm − fm), m = 1, 2, . . . , M − 1 вычислим векторы:
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
|
|
S2(m) |
|
|
|
S1(m) |
|
Sm = |
|
|
|
|
. . . |
||
|
|
S |
(m) |
|
|
N−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
, m = 1, 2, . . . , M.
2.Предположим UM = φa. По формуле Um−1 = RmUm + SM , m = M, M − 1, . . . , 1, вычислим последовательно искомые значения решения задачи (10.33): UM , UM−1, . . . , U1, U0.
Итерационный метод решения разностной схемы для задачи Дирихле
Выразим из схемы (10.24) значения Umn:
1 |
|
2 |
|
|
|
||
|
|
|
m = 1, 2, . . . , M − 1 n = 1, 2, . . . , N − 1 |
||||
Umn = 2(1+a) [Um+1,n + Um−1,1 + a(Um,n+1 + Um,n + Um,n−1) − h f(xm, yn)] |
(10.38) |
||||||
|
Um0 = φ(xm, 0), |
UmN = φ(xm, b) m = 1, 2, . . . , M −21 |
|
(10.39) |
|||
|
U0n = φ(0, yn), UMn = φ(a, yn) n = 1, 2, . . . , N − 1 |
|
|
||||
|
|
|
a = hl2 . |
|
|||
|
|
|
|
|
|
|
|
В равенстве (10.38) значение Umn выражается через 4 соседних значения Uij по пятиточечному шаблону. В итерационном методе решения полагают значения Umn во всех внутренних точках области Дh равным некоторым начальным условиям
Umn(0) , m = 1, 2, . . . , M − 1, n = 1, 2, . . . , N − 1.
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
Затем с помощью формул (10.38), (10.39) вычисляют новые значения Umn(1) , затем Umn(2) и т.д. до тех пор, пока максимальное отклонение значений сеточных функций на предыдущей и текущей итерациях не станет меньше по модулю, чем некоторая заданная точность:
MAX max |
U(l) |
U(l−1) |
|
< ε. |
(10.40) |
|
mn − |
|
|
|
|
= m,n |
mn |
|
|
|
Уравнение параболического типа. Явные и неявные конечноразностные методы
Рассмотрим параболическое уравнение |
|
|
|
||
|
∂T |
∂2T |
|
(10.41) |
|
|
|
= C |
|
+ f(t, x), |
|
|
|
∂x2 |
|||
|
∂t |
|
|
||
прототипом которого есть уравнение теплопроводности.
Пример. Стержень длиной L с постоянным по длине сечением погружен в изолирующий материал так, что с окружающей средой взаимодействует только его левый торец. В начальный момент времени t0
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
L
стержень |
X |
|
изоляция
Рис. 10.6.
стержень имеет равномерное распределение температуры по всей длине T = 0, а его левый торец скачком приобретает температуру T = 1000C. Требуется определить, как будет изменяться во времени температура в точках стержня, расположенных на разных расстояниях от его левого торца. Зависимость температуры от времени и координаты описывается следующим уравнением в частных производных, при f(t, x) = 0:
∂T = a2 ∂2T , ∂t ∂x2
где 0, x, L, t > 0, a2 — коэффициент теплопроводности стержня, зависящий от его теплопроводности, удельной теплоемкости и плотности a2 = l/(cr).
Граничные условия:
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
T (0, t) = T0+ = 100, |
∂T (L, t) |
= 0. |
||
∂x |
|
|||
|
|
|||
Начальное условие:
T (x, 0) = T0 = 0.
Общность численного решения можно повысить, если привести переменные к безразмерному виду:
x = |
x |
, t = |
ta2 |
, T = |
T − T0 |
. |
|
L2 |
|
||||
|
L |
|
|
T0+ − T0 |
||
Тогда исходное уравнение будет иметь вид (звездочки сверху у переменных опущены):
∂T = ∂2T . ∂t ∂x2
Изменится форма граничных условий
|
|
|
|
|
T (0, t) = 1, |
∂T (1, t) |
= 0 |
|
|
|||
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
∂x |
|
|
|
|
|
и начального условия T (x, 0) = 0. |
|
|
|
|
|
|
|
|
|
|||
Представив частные производные в разностном виде, получим: |
|
|
||||||||||
∂2T |
∂T |
= |
|
Ti+1,j − 2Ti,j + Ti−1,j |
|
|
Ti,j+1 − Ti,j |
= 0 |
||||
|
∂x2 |
− |
∂t |
|
|
h2 |
|
− |
|
k |
|
|
Результатом будет разностное уравнение
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
|
Tjm+1 − Tjm |
= c |
Tjm+1 − 2Tjm + Tjm−1 |
, |
|
∂t |
|||||
|
|
∂x2 |
|
которое аппроксимирует дифференциальное уравнение в точке (tm, xj).
Здесь производная по x заменена центрально-разностной, а производная по t — правой разностной производной. Отметим, что вычислительная сетка в данном случае не обязательно должна состоять из квадратных ячеек, так как h может и не быть равным k. Приняв r = k/h2, получим для Ti,j+1:
Ti,j+1 = rTi+1,j + (1 + r)Ti,j + rTi−1,j.
Это конечно-разностное выражение справедливо для всех внутренних узлов и позволяет явным образом вычислить температуру в момент времени t + k через температуру в момент времени t. Такая постановка задачи позволяет обойтись без системы совместных уравнений и, следовательно, не требует применения итерационных методов. Полученное решение будет содержать погрешности порядка h2 и k2, поскольку таковы порядки членов, отброшенных при конечно-разностной аппроксимации. В результате рассмотрения данной задачи может создаться впечатление, что величину шага для переменных можно выбирать независимым образом. К сожалению, это не так. Ранее отмечалось, что решения некоторых разностных уравнений могут сильно отклоняться от точного значения, и в этой связи было введено понятие устойчивости. Полученные решения дифференциальных уравнений в частных производных могут быть расходящимися или колеблющимися. Применительно к рассматриваемой задаче устойчивость решения можно проверить, меняя величину k. При r < 0, 5 решение устойчиво, но могут существовать колебания. Хотя при r = 0, 5 и наблюдаются некоторые колебания решения, удобно выбрать это значение, поскольку тогда рекуррентное соотношение приобретает вид
Ti,j+1 = 0, 5(Ti+1,j + rTi−1,j).
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
Ясно, что если расстояние между узлами по горизонтали уже выбрано, то расстояние между ними по вертикали не может быть произвольным. Так, если h = 0, 2, а r = 0, 5, то k = 0, 02. На рис. 10.7 показана сетка, соответствующая этим параметрам.
y,(i) |
|
6,25 25,0 56,0 |
|
||||
|
|
||||||
0 |
|
100 |
|||||
|
|
|
|
|
|||
|
|
|
|
|
|||
0 |
|
|
|
|
|
100 |
|
|
|
|
|
|
|||
0 |
|
|
|
|
|
100 |
|
|
|
|
|
|
|||
0 |
|
|
|
|
|
100 |
|
|
|
|
|
|
|||
0 |
|
|
|
|
|
100 |
|
25 50 75 |
x,(i) |
||||||
|
|||||||
Рис. 10.7. Двухмерная сетка с граничными и начальными условиями
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
(tm+1,xj)
|
|
|
|
|
|
|
|
|
(tm,xj-1) |
(tm,xj) |
(tm,xj+1) |
||||
(tm+1,xj-1) |
(tm+1,xj) |
(tm+1,xj+1) |
|||||
|
|
|
|
|
|
|
|
(tm,xj)
Рис. 10.8.
Разностная схема — это совокупность разностных уравнений, которые аппроксимируют начальное дифференциальное уравнение во всех внутренних узлах сетки и начальные условия в граничных узлах сетки. Разностную схему по аналогии с дифференциальной задачей назовём разностной задачей.
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
В данном случае разностная схема имеет вид
|
Tjm+1 − Tjm |
= c |
Tjm+1 − 2Tjm + Tjm−1 |
+ ϕm; |
(10.42) |
|
|
||||||
|
∂t |
|
|
∂x2 |
j |
|
j = 1, 2, . . . , N − 1, |
|
m = 0, 1, . . . , k − 1, hN = 1, |
kτ = t; |
|||
T m = µ1(tm), T m = µ2(tm), m = 0, 1, . . . , k; |
|
|||||
0 |
|
|
N |
|
|
|
Tj0 = g(xj), |
j = 1, 2, . . . , N − 1. |
|
|
|||
Эта схема — система линейных алгебраических уравнений, количество которых равно количеству неизвестных. Решать такую систему следует по узлам. Узлом называется множество всех узлов сетки, которые имеют одну и ту же временную координату. Если решение Tjm, j = 0, 1, . . . , N на узле m уже найдено, то решение Tjm+1 на m + 1 узле ищется по явной форме:
Ujm+1 = Ujm + µ(Ujm+1 − 2Ujm + Ujm−1), j = 1, 2, . . . , N − 1, µ = ct/h2.
Эта схема называется явной разностной схемой.
Погрешность разностной схемы определяется как разность между решением исходящей задачи и
нашей задачи |
Tjm+1 − Tjm |
= c |
Tjm+1 − 2Tjm + Tjm−1 |
. |
|
∂t |
|||||
|
|
∂x2 |
|||
Пусть T (t, x) — точное решение уравнения (10.41), начальными и граничными условиями. Если подставить это точное решение в (10.42), то оно удовлетвориться не полностью, а с некоторой погрешностью, которая называется локальной ошибкой дискретизации. Таким образом, локальная ошибка дис-
кретизации в точке (t, x) равняется |
|
|
|
|
|
|
|
|
||
e = |
T (t + t, x) − (t, x) |
− |
c |
[T (t, x + x) |
− |
2T (t, x) + T (t, x |
− |
x)]. |
(10.43) |
|
t |
(Δx)2 |
|||||||||
|
|
|
|
|
||||||
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель
Допустим, что нам известно |
значение точного решения |
T (t, x) |
для некоторого |
t |
и всех |
x [0, 1] |
и мы |
||||||||||
m |
|
|
|
|
|
|
|
|
|||||||||
желаем использовать (10.42) при ϕj = 0 |
для получения приближённого значения в момент t+ |
t. Если |
|||||||||||||||
|
|
|
|
|
ˆ |
|
t, x), тогда |
|
|
|
|
|
|
|
|
||
обозначить это приближённое значение как T (t + |
|
|
|
|
|
|
|
|
|||||||||
|
ˆ |
− T (t, x) |
|
|
c |
|
|
|
|
|
|
|
|
|
|
|
|
o = |
T (t + t, x) |
|
|
[T (t, x + x) |
|
2T (t, x) + T (t, x |
|
|
x)]. |
|
|
|
|||||
t |
|
− (Δx)2 |
− |
− |
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
||||||
Вычитая это равенство из (10.43), получаем
ˆ −
T (t + t, x) T (t + t, x) = e t.
Таким образом, погрешность, которую представили на одном шаге по времени в разностной схеме (10.42), равняется локальной погрешности дискретизации, умноженной на t. Значение e в (10.43) легко оценить при помощи t и x. Действительно, если считать T функцией только от t, считая x фиксированным, то из превращения в ряд Тейлора
T (t + t, x) = T (t, x) + Tt0(t, x)Δt + 0[(Δt)2]
получим
T (t + t, x) − T (t, x)
t
Рассмотрим T как функцию от x, считая t фиксированным. Превращая в ряд Тейлора T (t, x + x),
T (t, x + x), имеем T (t, x + x) |
− |
2T (t, x) + T (t, x + x) = T 00 |
(t, x)(Δx)2 |
+ 0[(Δx)4 |
]. |
||
Значит, и |
|
|
x |
|
|
|
|
|
|
|
|
|
|
|
|
T (t, x + |
x) − 2T (t, x) + T (t, x + |
x) |
= Tx00(t, x) + 0[(Δx)2] |
|
|||
|
|
||||||
|
|
(Δx)2 |
|
|
|
|
|
•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель