CC BY
23
УДК 519.632.4
ОБ АСИМПТОТИЧЕСКОЙ СХОДИМОСТИ НЕЯВНОГО ИТЕРАЦИОННОГО ПОЛИНЕЙНОГО
РЕКУРРЕНТНОГО МЕТОДА РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ
А. А. Фомин, Л. Н. Фомина
ON ASYMTOTIC CONVERGENCE OF IMPLICIT ITERATION LINE-BY-LINE RECURRENCE METHOD FOR SOLVING A DIFFERENCE ELLIPTICAL EQUATIONS
А. A. Fomin, L. N. Fomina
Анализируется результирующая матрично-векторная форма записи алгоритма неявного итерационного полинейного рекуррентного метода решения систем разностных эллиптических уравнений с пятидиагональными матрицами положительного типа. Путем оценки норм входящих в нее матричных операторов доказывается сходимость метода в асимптотическом случае, когда шаг разностной сетки стремится к нулю, а параметр компенсации близок к оптимальному значению.
Matrix-vector form of recording the algorithm of implicit iteration line-by-line recurrence method is analyzed for solving five-diagonal matrix systems of difference elliptical equations with positive type matrixes. By an estimation of norms of matrix operators entering into it convergence of the method in an asymptotic case is proved, when a difference grid step aspires to zero and compensation parameter is near to the optimum value.
Ключевые слова: разностные эллиптические уравнения, итерационный метод решения, сходимость метода.
Keywords: a difference elliptic equations, iteration method, convergence of the method.
Разностные аппроксимации дифференциальных постановок краевых задач тепло- и массопереноса приводят, как правило, к системам линейных алгебраических уравнений (СЛАУ) вида
АФ = I, (1)
где А - матрица системы, ф - искомый вектор решения, I - вектор правой части системы. Причем при выполнении простых и вполне разумных правил разностной аппроксимации получаемые матрицы СЛАУ имеют положительный тип, а именно: они неразложимы, имеют строчное диагональное преобладание, элементы главной диагонали строго положительны, а внедиагональные элементы неположительны [1, 2]. В случае использования регулярной сетки при разностной аппроксимации двумерных краевых задач матрица А зачастую имеет пятидиагональную структуру.
до сих пор не получили в литературе достаточно полного отражения. В этом плане следует отметить работу [6], в которой путем матрично-векторных преобразований системы показывается, что если в процессе итераций вектор приближенного решения -~(к)
Ф сходится к некоторому предельному значению
-м -м
Ф , то этот предельный вектор Ф есть также решение исходной системы (1). Цель настоящего исследования - доказать сходимость неявного итерационного полинейного рекуррентного метода в асимптотическом случае, то есть когда шаг разностной сетки 11 ^ 0, а параметр компенсации метода 9 близок к своему оптимальному значению 9ор1.
В общем случае для двумерной задачи пятиточечное разностное уравнение имеет следующий вид:
ар.. = aE..Ф!+1j
+aN.. +1
aS:®j
(2)
где 1 < i < n, 1 < j < m,
Рис. 1
Для решения подобных СЛАУ был разработан неявный итерационный полинейный рекуррентный метод, который при решении систем с различными матрицами, в том числе характеризующимися большими числами обусловленности, показал высокую практическую эффективность [3 - 5]. Однако вопросы теоретического обоснования этого метода
причем ар > аЕ + а№ + а3 + ам. Здесь п, т -
количество узлов сеточного разбиения области по координатам х, у соответственно. Пусть матрица А системы из таких разностных уравнений будет матрицей положительного типа. Ее структура приведена на рис. 1. Общее число неизвестных и, следовательно, размерность матрицы А равно числу N = пхт. Для удобства матрица разбита на отдельные квадратные клетки, каждая из которых на рис. 1 поименована как Аи. Если количество клеток, для определенности, равно пхп, тогда размерность каждой клетки - тхт. Здесь и далее для наглядности и не в ущерб общности изложения материала выбрана разностная сетка 5х5. Видно, что клетки А11 являются трехдиагональными, клетки А11+1, А11-1 - диагональными, а все остальные - нулевыми.
В [4, 6] показано, что суть рассматриваемого метода - это последовательность эквивалентных и
приближенных преобразований исходной системы уравнений, которая приводит ее в конечном счете к системе с четырехдиагональной матрицей. Ее структура представлена на рис. 2. Здесь черные крестики означают элементы, которые остались неизменными от исходной матрицы А, а белые крестики - видоизмененные элементы. Пусть при этом сами уравнения преобразованной системы в общем случае будут иметь следующий вид:
д(п) =
Е (п) - Е (п) и(п) Е(0) Е(р) И
(п ) (п)
+Е(р)н )
о
(п — 1 ш — 1)
(+)
(п)
в ) =
Е (п) — Е (п) и (п) Е(0) Е(р) И
д(п—1)
(п — 1 2) (-)
(п—1) ,
о
(п) (п)
+Е(р))и )
в
_(п—1 ш—1)
V)
в
в
(п—12)
(-)
(5)
^ Ф(к +1) = а Ф(к +1) ар. — аЕ..^г + 1 3
а8.. . Ф .-1
(к + 1)
(4)
(п) (п) (п)
где матричные коэффициенты д , в , £ и М( ) вычисляются по формулам:
£(п) =
(п) (п) (п)
Т71 т? ТТ
■Е / п \ -Е / \ и
(0) (Р)
£
(п—1)
(п ) (п)
+Е((р)) н )
__(п — 1 ш — 1) __(п — 1 2)
1(+) +1 (—)
М'" =
Еп) — Еп) и(п) Е(0) Е(р)И
1 (п—13) ш (п—13)
П М+ + ПМ—>
3=ш-1 3=2
Суть первичных матриц элементарных преобра-
(п) (п)
+Ер))и)
(.3) (. *) (. в)
зований М(+) , В(+) , Ь(+) .
(3 (. в) (. в)
М(—), в(—) , х(—) ^
и(0, Е)г
п(* )
+ Ф(*+11)
N . . 11 + 1
(3)
где, в отличие от (2), уже присутствует номер итерации к как следствие выполненных приближенных преобразований. При этом свободный член Ъ3-
является линейной функцией Ф(к), в противном случае метод будет прямым, что, в общем случае, невозможно. Очевидна цель проводимых преобразований - получить в конечной системе разложимую матрицу (рис. 2), откуда сразу же следует последо-вательно-поклеточная разрешимость преобразованной системы трехточечными скалярными прогонками.
-(к+1) -(к+1) -(к)
Пусть ДФ = Ф — Ф - приращение решения на к+1 итерации. Тогда в матрично-векторной форме записи преобразованная система уравнений согласно [6] имеет вид:
(п) -*(к+1) (п) -*(к+1)
д ф + ев дф
г(п) Л(к) аА(п) 1
= — £ Ф + М /,
^(0), Е(р) обсуждается в [6]. Здесь же, в первую очередь, следует обратить внимание на рекуррентную структуру этих формул относительно ко-
_(п) Л(п) (п) _ .(п)
эффициентов д , В , £ и М . Если по-
следовательно подставлять в них рекуррентные вы-
(п—1) ^(п—1)
ражения вида (5) сначала для д , в ,
Ап—1) а/п—1) ^(п—2) го(п—2) г(п—2)
£ и М через д , В , £
а/п—2) ^(п—2)
и ‘М соответственно, затем для д ,
.З(п—2) ,>—2) а/п—2) Ап—3)
в , £ и М через д
г(п—3) ш(п—3)
£ и М
В
(п—3)
и так далее, то есть выполнить обратное по отношению к работе [6] преобразование системы уравнений, то в итоге получатся следующие явные выражения для коэффициентов уравнения (4):
Л(п) Л(п) „ .(п) (п)
д + ев = М А + Л^. \
л(п) Ап) (п)
ев — £ = Лп ,
м (п) = П М? .
Здесь
п—1
Л=х
I=2
I+1
П МИ)
л (I) л (п)
ЛИ +Ли
(6)
(7)
(8)
(9)
-(.)
(.)
Ми = Е(0) +
(0)
Е (.) и (!) Е (р )н
1
П
3 = ш-1
М
_(. — 13) (+)
ш (.—13)
П М (—)
3 = 2
— Е
(10)
дН = EP)H(i) j e
л (і-i j)
П M(+)
j=m
(i-1 s-1) (i-1 s-1)
3(+) - L(+)
m-2
ь£
s =1
А (і-i j)
nM(-) j=1
(i-1 s +1) (i-1 s +1)
3(-) - L(-)
(ГГ)
При выводе (6) - (7) учтено, что
G
(n)
(n)
8B ) =
G
(n)
L
(n)
_(n) (n)
6Ф ) - L
Следовательно, система (4) теперь может быть записана в виде:
M(n) A
д
(n)
Mk+1)
Ф
(n) ~*(k)
(n) -
— д' Ф +M f
(Г2)
или
-(k+1) -(k)
Ф = sQ Ф +
M(n) A
д
(n)
M '■1 f,
M(n) A
д
(n)
-1
A (n)
дЕ - матрица пе-
рехода. Если теперь показать, что при определенных условиях спектральный радиус матрицы перехода
р(Se) < 1, то согласно теореме 1.9 [2] сходимость
метода будет доказана.
Для наглядности изложения материала имеет смысл привести примеры не тривиальных (ненулевых и неединичных) клеток матрицы A и матриц первичных элементарных преобразований для алгоритма LR1. Поскольку структура матриц преобразований идентично повторяется от линии к линии i = const, для определенности их содержание и структуры можно рассмотреть для сеточной линии i = 1.
Первые две клетки первого ряда матрицы A имеют вид:
Aii =
ap p її aN N її О О О
-as 12 ap p 12 -% i2 О О
О -as s 13 ap p 13 -aN 13 О
О О -as s 14 ap p 14 -% i4
О О О -as s 15 ap p 15
' - aE О 11 О О О
О a e E12 О О О
О О a e E13 О О
О О О a e E14 О
О О О О a e E15 )
A12 =
остальные клетки матрицы А выписываются подобным образом согласно уравнению (2) и схеме на рис. 1.
Первые диагональные клетки матриц первого
(11) (1т)
шага эквивалентных преобразований М+) и М( ) имеют вид:
M
(11)
11(+)
M
(1m) ii(-)
(12) (lm-1)
M+ и M(-)
M
(12)
11(+)
M
(lm-1) ii(-)
1 О О О
as sl2 1 О О
ap p11
О О 1 О
О О О 1
О О О О
1 О О О О
О 1 О О О
О О 1 О О
О О О + a Nim -
ap p1 m
О О О О 1
га эквивалентных преобра
1 О О О О
О 1 О О О
О as s13 1 О О
a p p12
О О О 1 О
О О О О 1
1 О О О О'
О 1 О О О
О О 1 aN+m-2 О
4-1
О О О 1 О
О О О О 1
остальные диагональные клетки матриц M,
.(11)
(+)
(1 2) (1 m)
M+) , M-)
M
(1 m-1)
- единичные матрицы, а внедиагональные клетки - нулевые матрицы. Здесь
a XT a N 1 a i^0 “N+i 1 ap p11
4» — ap -1 p+ m-1 aN1m a -1 sim
чае: ap — p+ j ap - p+j - as s+j aN+j-i
4 = ap p+j aN+j as sij+i
/ ар , а в общем слу-
/ %-1 ;
/ Ч+1.
Нетривиальные клетки матриц приближенных преобразований, отражающие механизм компенсации второго шага имеют вид:
Д12)
"12(+)
ь
(1т—1) 12(-)
в
(12)
12(+)
в
(1т—1) 12(—)
— ПГ2 О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
О
— 2п О О
О
О
О
О
О
О О
О О
О О
О О
О О
О О
-пт—1
О
О
Здесь п“2 =
П1т —1 =
в общем случае
О О О'
О О О
1?2 ПГ2 О О
О О О
О О О
О О О'
О О О
П1т—1 2 1 —1 О
О О О
О О О
/ аРц К аБ / Б13 аР ), Р12 /
/ аР1т К ая Я1т—2 / Ч—
Пц = а5Е„
у М1,—1
= Умв1.--------------;
3 13 ?л..
а — ___13 +1 I а
^БЕ, ~ 1 Е ,
13 а „ ' 13—1
29а
БЕ„
1мЕи =
1 (3+1 + 291Ш13+1
&е — а
Я13-
13—1 а г
7е1 3- = % 97№Е13-+1
'13—1
7Р13+1
Все остальные клетки матриц Ь
(12)
(+)
ь
(1т—1)
(-)
(12) (1т—1)
В(+) и В( ) - нулевые. Алгоритмы ЬЯ2 и
[4,5] отличаются от ЬЯ1 только содержанием (но не
сутью!) матриц в(+) и В^^') - их элементами также
являются составляющие соответствующих компенсационных формул.
(2)
Матрица Н имеет на главной диагонали единичные клетки, а первая клетка второго ряда следующая:
Н (2) = Н 21 —
аШ21 /
О
21 Х11
О
О
+
О
О
О аШш / (аР12 + УР12 аР12 )
О О Ч3 / (
О О
О О
(2)
все остальные клетки матрицы Н - нулевые.
И наконец, матрицы «разделения» рядов клеток:
(2)
Е(О) - почти единичная, у которой только диаго-
нальная клетка Е
(2)
22(О)
нулевая;
(2)
Е(р) - почти нулевая, у которой только диаго-
(2)
нальная клетка Е22(р) - единичная матрица. По-
нятно, что Е,
(2)
(О)
22(р)
(2)
Е(р) — Е - единичная матрица.
а№,т / аР
Здесь следует обратить внимание на то, что все коэффициенты и их комплексы по построению положительны [4 - 6], иными словами, знаки элементов рассмотренных клеток матриц представлены явно.
Для оценки величины оптимального значения параметра компенсации 9ор необходимо проанализировать формулу вычисления свободных членов Ъ3-
преобразованной системы уравнений. Она имеет вид [4]:
г+13'
г+1}
г+13
Р.
1рц+1
р.- 1 + аБЕ ( Ф* . 2 —29 Ф* . 1 + 9 Ф* .
^%3—1 1 БЕц—1\ г+13—2 г+13—11 г+13
к,+1+-1яе„+1 (+^29 ф;+ц+1+9 ф;+ц)
Р
а
Смысл входящих в (13) коэффициентов подробно изложен в [4], здесь же следует обратить внимание на две круглые скобки, отражающие явную зависимость Ъ3- от Ф*, поскольку целью оптимизации параметра 9 является минимизация выражений в этих круглых скобках.
В случае линейной экстраполяции, приращение решения в точке (1+1, >2)
ДФг+13—2 — Ф*+13—2 — Ф*+13—2 определяется экстраполяционным выражением 2ДФ%+13-—1 — ДФ%+1з. В силу своей приближенности, оно не позволяет точно определить приращение ДФ%+13—2. Поэтому вместо него реально определяется приближенное значение величины ДФ%+1з—2:
ДФ % +13 —2 — 2ДФ г+13—1 — ДФ%+13 . (14)
Для оценки привнесенной погрешности, функции ДФ%+13 и ДФ%+13—1 раскладываются в ряд Тейлора в окрестности точки (1+1, >2):
ДФ%+13
+ЛДФ-
AФ і+lj-1 = ^ t+1j-2
h2
+1 j-2 + ~ AФ Г+1 j-2 + О(h3 ),
(l5)
ДФг+13 — ДФг+13—2 +
+2ЬДФ ^+13—2 + 2Ь2 ДФ ^+13—2 + О (Ь3).
Комбинация (15) и (16) дает
ДФ г+13—2 — (2ДФ г+13—1 — ДФ г+13 ) +
+Ь2ДФ г"+13-—2 + О (Ь3).
Введение множителя 9 призвано устранить разницу между (14) и (17). Следовательно, оптимальное значение 9ор4 параметра компенсации должно быть
(16)
(17)
таким, чтобы AФ,
+1j-2
^opt A<®s+1j-2
или
разностной аппроксимации краевой задачи в области О, покрытой регулярной сеткой с шагом И. При этом А - ленточная пятидиагональная матрица положительного типа. Тогда при условии, что решение исходной краевой задачи в О есть функция, дважды непрерывно дифференцируемая и ограниченная вместе со своими производными, итерационный метод (12) сходится при И 0 и параметре компенсации 9 = 9ор4.
Доказательство. В работе [6] было показано, что используемая в методе последовательность линейных эквивалентно-приближенных преобразований системы уравнений сохраняет свойство невырожденности матрицы со строчным диагональным преобладанием на каждом этапе преобразований. Следовательно, итоговая матрица в (12)
A
(n )
является невырожденной и для
нее существует обратная матрица, которая, в силу конечномерности рассматриваемых пространств векторов Ф , порожденных сеточным разбиением О, является ограниченной [7]. Иными словами, существует константа МО:
M(n А
A
(n)
< Mо .
Нетрудно показать, что М0 не зависит от сеточного разбиения области О. Действительно, при любом сеточном разбиении О структура матрицы
M’ А
A
(n )
остается четырехдиагональной с
(2ДФ г+13—1 — ДФ г+13 ) + Ь2ДФ Г+13—2 +
+О(Ь3 ) = 9ор, (2Дфг+13—1 — ДФг+13 ),
откуда следует, что
9ор* — 1 +[ДФ'+13—2 / ( 2ДФг+13—1 —ДФг+13 )]Ь + О*3).
Поскольку стоит задача только оценки порядка величины 9 в зависимости от И, то окончательно можно записать
9ор; — 1 + О (Ь2). (18)
Поскольку по определению компенсация не может быть более чем полной, то есть 9 < 1 [2], то
из (18) следует, что 9ор1 е [ 1 — О(Ь2 ),1].
Проведенные рассуждения позволяют сформулировать следующую теорему о сходимости неявного итерационного полинейного рекуррентного метода в асимптотическом случае.
Теорема. Пусть система линейных алгебраических уравнений (1) корректно определена путем
ненулевой главной диагональю. Такую систему уравнений всегда можно привести к виду, при котором ее коэффициенты при неизвестных на главной диагонали будут иметь порядок 0(1) + 0(ИЧ), где
q = 1,2,3,_, путем умножения каждого уравнения
системы на И в соответствующей степени. При этом коэффициенты на побочных диагоналях не могут иметь более низкий порядок по И (то есть, не возможен случай ч < 0), иначе нарушится свойство строчного диагонального преобладания. Тогда оценка
первой нормы матрицы
A
(n )
будет сле-
дующей [S]:
An)
(n)
aE,, ) >
> max (aP ) = —0 + O (hq) > —- > 0,
г,з V P ' 0 ' ' h ^0 2
где mo - const, не зависящая от h. В этом случае
по теореме 1 главы 5 [7]
-1
(n) (n)
M А + Aej
< —= Mo.
С другой стороны, рассмотренная ранее струк-
(і)
тура и содержание матриц M(+) , M(
Лі)
Ч-)
H , E,
(О)
1
р( ^
Е(р) говорят о том, что элементами данных матриц
являются либо нули, либо единицы, либо положительные отношения меньше единицы, не зависящие от 9 и И в отрицательной степени (см. рассуждения выше), причем в одну строку располагается не более двух ненулевых элементов. Поэтому сами матрицы и их линейные комбинации являются ограниченными. Тогда, как следует из (11), оценка норм мат-д (г)
риц Дн сводится к оценке норм матриц
( * *) (і в) ?(+) - 1(+)
(і в) (і в)
V) - V)
элементы которых есть нули либо составляющие формулы компенсации.
По определению нормы матрицы [7, 8] приме-
нительно, для примера, к матрице
( і в) (і в)
вв(+) - ь(+)
можно записать:
(і в) (і в)"
вв(+) - V)
8ир { |фг
ІІФІИ
= ^ир
||ф| 1=1
+18 + 1 — 29Фі + 18
(і в) (і в)
^+) - Ь(+)
Ф
вФ
+1в-1
I}-
/
. Подставляя в по-
лученное соотношение 9 = 9ор из (18) и учитывая, что 9ор < 1, можно получить следующую оценку:
( і в) ( і в)
в(+) - х(+) „^р { {
11Ф1 и
+1в+1 2вФі+1в
= п^Р І Пі,
||ф| Іи11
= ||Т {
||ф||и1
(+1в+1 - 2фі+1. +
-(+1в - фі+1в-1
вФ
ф
+1в-1
1}
+1в-1
О (к2
И-2Ф+Ь - фі +Ь-1 )\° и2 ) =
(к2 )} = о (к2 ).
фі+1в + 1 - 2фі+1в
ф
+1в-1
посколь-
можно сразу записать:
(-)
(-)
= О(к2).
сации. Отсюда следует, что
Д
иЖ(г3) 1,ж(г3)
ограниченности матриц М(+) , М(
= О (к2) в силу
(і)
(-)
Н
(п)
бинаций соответствующих ненулевых рядов клеток из матриц дН и ограниченных матричных коэффициентов в виде произведений матриц мН ), то оцен-
Д(п)
ка нормы результирующей матрицы ДЕ не изменится. То есть при 9 = 9
ор1
Д
(п)
= О(к2).
Следовательно, оценка спектрального радиуса матрицы перехода $е имеет вид:
р(е ) = р
М(п) А + Д<п)
<
<
Ап) (п)
М А + ДЕ
-1
Д
(п)
М(п )А + Д
(п)
Д
-1
< (п)
Д
(п)
<
(19)
< М0 ■ О (к2).
Здесь учтено, что вектор Ф является покоординатно «гладким», то есть
ку искомое решение задачи непрерывно и ограничено вместе со своими производными. По аналогии
(г*) (г*)
— ь
Анализ структуры формулы (11) говорит о том,
д(%) 14 Т?(г)
что в матрице Дн : 1) из-за множителя Ер) только
1-я строка клеток содержит ненулевые клетки; 2) каждая строка ненулевой клетки содержит не более одного набора составляющих формулы компен-
В итоге, поскольку из (9) следует, что ДЕ есть суммарный набор по рядам клеток линейных ком-
Понятно, что соответствующим подбором малого И можно всегда добиться выполнения неравенства р (Б9) < 1, откуда следует доказательство теоремы.
Замечание 1. Очевидно, что метод сходится и
при 9 е [ 9оР1 — О(Ь2 ),11 .
Замечание 2. Если рассмотреть алгоритм ЬЯ2 [4] (используются квадратичные экстраполяционные формулы компенсации в отличие от линейных для ЬЯ1 (14)), то по аналогии нетрудно убедиться, что оценка спектрального радиуса матрицы перехода
будет р (Бе )< еопзЬ ■ О (к3), а сам метод будет
сходиться при 9 е [ 9ор1 — О (Ь3) ,1 ].
Замечание 3. Из оценки (19) следует, что при И ^ 0 спектральный радиус р () становится
сколь угодно малым и на точной арифметике итерационный метод (12) стремится к прямому.
Литература
1. Патанкар, С. Численные методы решения задач теплообмена и динамики жидкости [Текст] / С. Патанкар. - М.: Энергоатомиздат,1984. - 152 с.
2. Ильин, В. П. Методы неполной факторизации для решения алгебраических систем [Текст] / В. П. Ильин. - М.: Физматлит, 1995. - 288 с.
3. Фомин, А. А. Сравнение эффективности высокоскоростных методов решения разностных эллиптических СЛАУ [Текст] / А. А. Фомин, Л. Н. Фомина // Вестник Томского государственного университета. Математика и механика. - 2009. -№ 2. - С. 71 - 77.
4. Фомина, Л. Н. Использование полинейного рекуррентного метода с переменным параметром компенсации для решения разностных эллиптических уравнений [Текст] / Л. Н. Фомина // Вычисли-
тельные технологии. - ИВТ СО РАН. - 2009. - Т. 14. - № 4. - С. 108 - 120.
5. Фомин, А. А. Об одном варианте полиней-ного рекуррентного метода решения разностных эллиптических уравнений [Текст] / А. А. Фомин, Л. Н. Фомина // Вестник Томского государственного университета. Математика и механика. - 2010. -№ 2. - С. 20 - 27.
6. Фомин, А. А. Обоснование корректности неявного итерационного полинейного рекуррентного
метода решения разностных эллиптических уравнений [Текст] / А. А. Фомин, Л. Н. Фомина // Вестник Кемеровского государственного университета. -2011. - № і. - С. 39 - 45.
7. Самарский, А. А. Методы решения сеточных уравнений [Текст] / Е. С. Николаев, А. А. Самарский. - М.: Наука, 1978. - 590 с.
8. Самарский, А. А. Численные методы [Текст] / А. А. Самарский, А. В. Гулин. - М.: Наука, 1989. -432 с.
