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

УДК 5І9.632.4

ОБОСНОВАНИЕ КОРРЕКТНОСТИ НЕЯВНОГО ИТЕРАЦИОННОГО ПОЛИНЕЙНОГО

РЕКУРРЕНТНОГО МЕТОДА РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ

А. А.Фомин, Л. Н. Фомина

ARGUMENTATION OF CORRECTNESS 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 proposed for solving five-diagonal matrix systems of linear algebraic equations with a positive type matrixes. Structure of a resultant operator converting original system to the one with an about upside triangular matrix is deduced by means of stepwise stable transformations. Argumentation of correctness of the method is given.

Ключевые слова: разностные эллиптические уравнения, итерационный метод решения, корректность метода.

Keywords: a difference elliptic equations, iteration method, correctness of the method.

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

АФ = / .

Здесь А - матрица системы, Ф - искомый вектор решения, f - вектор правой части системы. Подобные матрицы обладают высоким порядком и ленточной структурой [1]. Для многомерных задач, полученные таким образом СЛАУ разрешаются, как правило, итерационными методами.

В работах [2, 3] излагаются различные алгоритмы и рассматриваются результаты тестовых расчетов неявного итерационного полинейного рекуррентного метода для случаев, когда А - матрица положительного типа. Данный метод демонстрирует свою высокую эффективность, однако в силу того, что в общем случае записать его в канонической форме [1]

рк+1

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

полнении условия Ф

k+1

_ Фк

О ), полу-

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

ҐХ X XXX XXX XXX X X X X X X X Г \ \

X X X X X XX XXX XXX XXX XX X X X X X 1 )

X X X X X X X XXX XXX XXX XX X X X X X

( л X X X X X XX XXX XXX XXX XX X X X X X

V ч J X X X X X X X XXX XXX XXX XX)

Рис. 1.

Пусть для двумерных задач пятиточечное разностное уравнение имеет следующий вид:

ЛР..Фг] = ЛЕ.. Фг+1] + Фг-1] +

г] г]

+ЛМ.. Фг]+1 + ЛБ.. Фг]-1 + Ьг],

г] .]

1 < . < п, 1 < ] < т,

где

причем

щ + а3 + ам . Здесь п, т - количество

узлов сеточного разбиения расчетной области по координатам х, у соответственно. Тогда матрица А имеет пятидиагональную структуру (рис. 1). Общее число неизвестных и, следовательно, размерность матрицы А равно числу N = пхт. Как известно, такая матрица может быть разбита на отдельные квадратные клетки, количество которых, для определенности, будет пхп, тогда размерность каждой клетки - тхт. Здесь и далее для наглядности и не в ущерб общности изложения материала выбрана разностная сетка 5х5. Понятно, что в рамках рассматриваемого случая только клетки на диагонали и около диагональные будут отличны от нуля, все остальные - нулевые. Поклеточная структура матрицы А и состав двух клеток (для примера) имеют вид:

А =

А11 А12 А13 А14 А15

А21 А22 А23 А24 А25

А А А А А

31 32 33 34 35

А А А А А

41 42 43 44 45

,А51 А52 А53 А54 А55 у

А11 =

ар — ам 0 0

р11 л 11

—а а& — аА

Ь12 р12 л

0 — аЬ

ар13 ал

0 'ї 0 0

0 0 — аЬ

ар аN

р14 л14

А12 =

0 0 0 — аЬ

Ь15

Г—аР 0 0 0 0 'ї

-Ъ11

0 — ар 0 0 0

Ъ12

0 0 — аЪ 0 0

Ъ13

0 0 0 — ар 0

Ъ14

0 0 0 0 — а„

Преобразования СЛАУ начинаются вдоль направления по координате у для ] от 1 до т на линии I = 1. Первые три уравнения для узлов (1,1), (1,2) и (1,3) имеют следующий вид

аР11 Ф11 = аЕ11 Ф21 + аМп Ф12 + Ь11;

Ф12 = аЕг,ф 22 + аМ,„ ф13 + %_ Ф11 + Ь12 ;

аР,3 Ф13 = аЕ,3 Ф 23

аЫ 13 Ф14

aS,, Ф12

Уравнение (3) умножается на отношение

/

и подставляется в правую часть уравне-

ния (4) вместо третьего слагаемого. После приведения подобных получается соотношение

(аР12 — aN1laS12 / аР11 )Ф12 = %12 ф13 + аЕ12 Ф22 + +(аЕ11аЯ12 / аР11 )Ф21 + Ь11аЯ12 / аР11 + Ь12'

Переобозначения коэффициентов позволяют переписать последнее уравнение в виде:

аР12 Ф12 = aN12 Ф13 + аЕ12 Ф 22 + aSE12 Ф 21 + в12>

где аР = аР — аы а„ / аР ,

Р12 Р12 ^11 °12 Р11

aN — aN , аЕ — аЕ ,

7У12 7У12 12 12

аЭЕ12 = аЕ11аЯ12 / аР11 , в12 = Ь11аЯ12 / аР11 + Ь12 .

Подобное преобразование проводится и для уравнения с центральным узлом (1,3). Для этого в уравнении (5) исключается слагаемое Ф:2, путем подстановки вместо него правой части (6) умноженной на отношение а3 / аР . После приведения

подобных получается соотношение

(аР13 — а^2 а£13 / аР12 )Ф13 = aN 13 Ф14 +

+аЕ13Ф23 + аЕ12аБ13 / аР12Ф22 +

+а

. О Т7І Ы/О

^12 ‘-'ІЗ

/ аР12 Ф21 + в12а^„ / аР„ +

В матричном виде эти же самые первые два шага преобразований полинейного рекуррентного метода (вывод уравнений (6) и (7)) представляют собой воздействие на исходную систему (1) двух не-

вырожденных операторов М и М [2]:

При этом только диагональные клетки матриц

операторов М(+) и М++^ отличны от нуля, из них

первые клетки в верхнем ряду имеют следующий вид:

М.

11

М.

11

1 0 0 0 0

1 сС (3 2 СгГ в 0 0 0

= 0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

' 1 0 0 0 0'

0 1 0 0 0

= 0 а„ /аР Ь13 Р12 1 0 0

0 0 0 1 0

0 0 0 0 1

или, более кратко, структуру:

М-

11

М.

11

X 0 0 0 0 XX 0 0 0 0 0 X 0 0 0 0 0 X 0 0 0 0 0 X X 0 0 0 0 0 X 0 0 0 0 XX 0 0 0 0 0 X 0 0 0 0 0 X

остальные диагональные клетки - единичные мат-

рицы. Понятно, что матрицы

м(+) и м(+) представляют собой частный случай так называемых элементарных нижних треугольных матриц, с помощью которых матрица исходной СЛАУ преобразуется к треугольному виду [4]. Такие матрицы не вырождены, а само преобразование будет устойчивым в силу свойства строчного диагонального преобладания исходной матрицы А. Следовательно, клетки произведения имеют следующую структуру:

Чз-

М(Vм(+) А

(12) (11) М ((+)М (+) А

11

12

XX 0 0 0 0 XX 0 0 0 0 XX 0 0 0 XXX 0 0 0 XX

X 0 0 0 0 XX 0 0 0 XXX 00 0 0 0 X 0 0 0 0 0 X

Произведение М( ) М( ) А можно представить

в виде суммы С(+) + Ь(+) , причем матрица

г

Ь(+) состоит из нулей, кроме второй клетки в верхнем ряду, которая имеет всего лишь один ненулевой элемент - коэффициент «внешаблонного» узла [2]:

' 0 0 0 0 0'

0 0 0 0 0 Х0 0 0 0 0 0 0 0 0 0 0 0 0 0

12

О

а матрица О(+) содержит все остальные элементы

матрицы М(+) М(+) А . Используя данное представление, система (8) перепишется в виде:

(12)- (12)-*

О Ф + ь Ф =

О

П М+

3=2

I .

При помощи одного из способов аппроксимации значения искомой функции во «внешаблонном» узле производится замена члена уравнения, содер-

Ь Ф ^ Ь Ф + ВБ ДФ

жащего его:

где

-^

ДФ = ((к+1) — Ф(к))- приращение решения. При этом алгоритм принимает итерационный

характер. Структура матрицы Б++^ отражает собственно механизм аппроксимации значения искомой функции во «внешаблонном» узле (2,1) третьего уравнения системы, то есть ее ненулевые элементы являются коэффициентами приближенных преобра-

зований уравнений системы. Так как Б(+) умножа-

Мк+1)

ется на приращение решения ДФ , то при сходящемся решении данное произведение стремится к нулю.

Вторая клетка в первом ряду матрицы Б(+) в

случае линейной экстраполяции [2] имеет вид:

В

12

' 0 0

0 0

„ аБЕ,

0 2

а

0 0

0 0

Е12С*13

0 0

или, более кратко, структуру

В

0 0 0 0 0 0 0 0 0 0 0 XX0 0 0 0 0 0 0 0 0 0 0 0

Все остальные клетки, включая диагональные

БЦ(+) , являются нулевыми матрицами. Очевидно,

что то же самое можно сказать и в случае квадратичной экстраполяции [2] с той лишь разницей, что вид второй клетки в первом ряду будет несколько иным. А именно:

В

12

0

0

0 -3-

0 0 0

0 0 0

о г?! с о *Е12 *13 (а о & сь о о *Е12 *13 аБЕг^

а р р12 0 3 ар р12 0 ар рг 0

0 0 0

или, более кратко, структуру

В

12

0 0 0 0 0 0 0 0 0 0 0 XXX 0 0 0 0 0 0 0 0 0 0 0

Тогда понятно, что комбинация 9В( ) , добав

,

ленная к матрице , не изменит содержимое

диагональных клеток последней. Соответственно

элементы матрицы ВБ(+) и последующих подобных матриц в процессе преобразований с помощью элементарных нижних треугольных матриц не будут влиять на элементы главной и двух прилегающих побочных диагоналей итоговой матрицы сис-

не изменяет свойство вырожденности (не вырож-денности) преобразованной СЛАУ.

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

темы. Следовательно, добавление матрицы 9В

о; ф + 9 в; дф

(12) 2(к) +ь(+)Ф “

П М

3=2

12

Последующее эквивалентное преобразование (11) также выражается в виде воздействия невыро-

жденного оператора М(+) , построенного по аналогии (9), на обе части (11), а именно

О Ф

О

_(1 т—1) -^(к) +ь(+) Ф =

_

вВ(+) дф

1

П М(^

3=т-1

I.

<>Ч>Ф +

+м,18,ьт,Ф14 =

^М

вМ, ’ В, ’ дф +

1

П М+

3=3

I.

При этом структуры первой и второй клеток в

„

первом ряду матрицы О(+) и вторых клеток в

_

первом ряду матриц Ь и Б( ) будут сле-

Как и на предыдущем шаге алгоритма пусть

д#

м(|) О,, = О, + Ь, ^ , где матрица Ь(, ч имеет

дующими:

О

структуру аналогичную ь(+) , то есть

11

XX0 0 0 X 0 000

0 XX 0 0 (1 т —1) X X 000

00 X X0 , О12(+) = 0 XX00

000 X X 12(+) 00 XX 0

000 0 X 00 0 XX

г

12

Замена Ь(+) Ф

0 0000 0 0000

0 0000

0 0000 . _(1 т—1) 0 0000

0X 000 0 0000 Ь12(+) = X 0000 0X 000

00X 00

• Ь Ф

Ь

І) дф

позволяет записать следующее соотношение:

О Ф

О

В,

М

дф

-

___

В 12

!(13) + V) +

-

Ф =

П м

3 = 3

В

_

12

или, при условии

В

= В

М

00 0 00 00 0 00 0 XX 00 0 0 XX 0 000 XX

0 0 0 0 0 0 0 0 0 0 0 XXX 0 0 0 XXX 0 0 0 XX

- линейная экстраполяция,

квадратичная экстрапо-

т<13) = г

Ь

О Ф

О

-(13) Чк) +Ь(+))ф =

вВ(+) ДФ

1

П М™

3=3

ляция.

Аналогичным образом рассматривается проход по локальному направлению на линии / = 1 в сторону уменьшения индекса / Выполнение подобной цепочки преобразований по ] ^ 2 приводит, в итоге, к следующему матричному уравнению

I.

В общем случае для произвольного ] = 3 , при увеличении индекса ], на линии / = 1 будет иметь место соотношение

-

О Ф

О

—(12) “*(к) +ь(—)Ф =

-(12) Чк +1) вВ( дф

П М

3=2

I,

О Ф

О

-(17) Чк) +ь(+))ф =

вВ+ дф I! м+

3 = 7

Л#

где М(—) О(—) = О(—) + Ь ,

В

I,

= В

л/и )п

где М(+) О(+) = О(+) + ь

+м—/ в((—;

(12) —(13) М(—) Ь(-) .

В 7 = В 7

В

) = т(и) Ь(+) = Ь(+) '

(17 ) -(17—1) +М(+) В((+) =

М(+) Ц+) ) .

Т<12) = г

Ь

А в общем случае при обратном проходе для произвольного ] = 3 имеют место соотношения

Окончание прохода по локальному направлению ] = 1, т — 1 на линии / = 1 приводит к итоговому матричному уравнению:

М

В

7<1и) ь

= ь

= О 7 + Т{17) О

ЬМ—/ В((—) }

М (17) Ь7+1) М(—) ь(—) .

Структуры первой и второй клеток в первом

„

ряду матрицы О’( ) и вторых клеток в первом ряду

?(12) д(12) Л

матриц Ь(-) и Б(-) будут, соответственно, следующими:

G

ii

X оооо XX ооо

XX ооо , g(1(!\ = о XX о о

о XX оо о о XX о

оо XX о ії(-) ооо XX

ооо XX ооо 0 X

L 12) Li:

B

ії

B

ії

XX0 оо оX X оо о оXXо о оо оо о оо оо

XX о о о XX/о о о XXX о о о о о о о о о о о

оо X оо оооX о оооо X ооооо ооооо

линеиная экстраполяция,

- квадратичная экстраполяция.

Здесь еще раз следует обратить внимание на существо переходов

r(i j ]2 L Ф ^

L(+r

r(l j) 2

и L )Ф

Jij) ^

L Ф

L(+)Ф

Jij) -r(k) L Ф

L

0B((+) ДФ

(lj) -(k+l) 0B((-j) ДФ ,

r

в которых матрицы L(+) и L(-) содержат коэффи

циенты при неизвестном во «внешаблонном» узле, а

,

,

матрицы B(+) и B(-) выражают собой механизм

аппроксимации этого неизвестного через неизвестные в узлах основного шаблона. В качестве таких механизмов ранее были рассмотрены так называемые линеиная и квадратичная экстраполяции, хотя понятно, что они не единственные и могут быть другие способы выражения неизвестного во «внешаблонном» узле через неизвестные в узлах основного шаблона [3]. Как будет показано далее, с точки зрения обоснования корректности метода, основным является тот момент, что каков бы не был в общем случае механизм аппроксимации (компенсации), все коэффициенты этого механизма располагаются только в матрицах В. Поэтому в последующих рассуждениях механизм аппроксимации более не детализируется.

В силу совпадения верхних наддиагоналей кле-

и нижних поддиагоналей клеток

ток Aii и Gii

Aii и Gii(;

n(im-i) + -(i:) разность G , + G,

- A

ii

пред-

_) 11

ставляет собой диагональную клетку с положительными элементами. Тогда понятно, что вычитание из

-

позволяет получить систему с «зеркальными» шаблонами [2] на линии I = 1

3

„

_

'L

-

_

-0B

ДФ

-

П M^+n M-

j=m -l j =1

Рис. 2.

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

,

'11

,

Т11

-

G „\ - a11 ] + в[в\

a" = G

тогда A1 = [G ства нулю диагональных клеток матриц B

^(1() л 1

Gll(-) - All ] в силу равен-

B

причем, как уже отмечалось ранее, клетка

,

А11 - диагональная матрица, следовательно к ней

просто найти обратную матрицу [аЦ ] . Далее

вводится нижняя треугольная матрица преобразова-

ния («деления») Н , у которой первая клетка вто-

г(2) , Г,(1)-1-1

рого ряда

H11 =- ЛКї ]-

все диагональные

клетки - единичные матрицы, а остальные клетки -

нулевые матрицы. Умножение (17) слева на Н и приводит к уравнению вида:

Н 1О + Н О— — Н 1А

+в

Н В

н (2) в;—2;

дф

Ф

-

Н V) + Н Ь

-

Ф =

18)

Н () П МЇЇ

Н

I,

3=т-1 3=2

причем структура матрицы уравнения

(1), записанное для Ф , с уравнением (18) таким образом, чтобы в (1) второй ряд клеток был заменен на второй ряд из уравнения (18). Для этого вводятся

р

две диагональные матрицы Е(0) и Е(р) , причем у

р

Е(0) почти все диагональные клетки - единичные

матрицы, кроме Е22(0), которая является нулевой, а

у Е(р) все наоборот - все диагональные клетки ну-

левые, кроме Е22(р), которая является единичной.

Понятно, что Е

Е(р) = Е - тождественный

оператор. Полученная при этом сумма:

Н О + Н О— — Н А

-

+в

+

Н В++ + Н в

-(к+1) дф +

(2)Л1т— 1) (2) -(12)'

Н і

-

н П М+) + н ПМ—) — н

3 =т-1 3=2

или

О

_

В

_

^ +Ь

—(к+1) дф

О12) 0—)

—(к+1) ф

Чк)

Ф =

ЕР) Н‘

е0)) —еР)н12 1+е» н2 II М1?+пМ

3=т-1

3=2

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

Я -

+е

О

А

в = е(Р; ну>

В

В

у~

Ь = Е

Е (2) _ Е (2) Н (2) Е(0) Е(р) Н

(2) (2) Ч)Н

П М

3=т-1

П Мїї

3=2

уравнение (19) окончательно можно переписать в виде

Я Ф

(2) -(к+1) вв дф

В заключении выкладок, связанных с линией I = 1 следует заметить, что все проведенные преобразования исходной системы представляют собой линейные комбинации уравнений, в которых в качестве сомножителей всегда используются коэффициенты по модулю меньше единицы в силу изна-

чального свойства построчного диагонального преобладания матрицы СЛАУ. Следовательно, проводимые преобразования: 1) устойчивы; 2) сохраняют свойство построчного диагонального преобладания; 3) сохраняют противоположность знаков диагональных и внедиагональных элементов. Иными словами матрица системы

в в

(рис. 4) продолжает

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

матрицах Я , Ь только в матрице в

,

а приближенные -

Повторение проведенных эквивалентноприближенных преобразований по всем расчетным

линиям і = 2, п — 1 (рядам клеток матрицы СЛАУ) приводит к тому, что исходная система преобразуется к виду:

Я иФ

в )дф

і Ап) л(к) пл(п) 7 +Ь Ф = м I,

в которой эквивалентные преобразования содержат-

Лп) Лп) п/(п) ся в матрицах у , Ь , М , а приближенные

^->

только в матрице В , причем матрица щ(п) (в) = у(п) + 0(п) по свой структуре четырехдиагональная и положительного типа (рис. 5), откуда следует последовательно-поклеточная разрешимость полученной системы.

В случае сходимости итерационного процесса -

имеет место стремление ДФ ^ 0 , в результате которого влияние члена, содержащего приближенные преобразования, становится несущественным.

Т-Г Лп) Г

Поскольку операторы у , Ь и М - суть

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

вид оператора воздействия на первоначальную систему уравнений:

(п) (п) (п

4) + 4!Н

1

П М

3 =2

— Е

Е

Е(0) + Е(р) Н

1

П М

3=т-1 3 =2

- сЯ , ^(2^(2) 1

... X - Е(0) + Е(„) Н

—Е

І0)

-ч*) -

АФ — I

\р)

=0.

Т-Г

П М

=т-1 3=2

— Е

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

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

Щ(п)(е)(к+1 — Фк ) = ^(Фк — /),

»(п)/ \

где щ (В) - удобно разрешаемый оператор с четырехдиагональной почти верхнетреугольной матрицей,

^ - оператор невырожденных эквивалентных преобразований (произведение фигурных скобок в (22)),

а В - итерационный параметр компенсации.

|(п) / \

При этом матрицы операторов щ (В)и ^ выписываются явным образом. Нетрудно видеть, что в этом случае вопрос о корректности метода разрешается естественным образом.

Литература

1. Ильин, В. П. Методы неполной факторизации для решения алгебраических систем [Текст] / В. П. Ильин. - М.: Физматлит. - 1995. - 288 с.

2. Фомина, Л. Н. Использование полинейного рекуррентного метода с переменным параметром компенсации для решения разностных эллиптических уравнений [Текст] / Л. Н. Фомина // Вычислительные технологии. - ИВТ СО РАН. - 2009. -Т. 14. - № 4. - С. 108 - 120.

3. Фомин, А. А. Об одном варианте полинейно-го рекуррентного метода решения разностных эллиптических уравнений [Текст] / А. А. Фомин, Л. Н. Фомина // Вестник Томского государственного университета. Математика и механика. - 2010. -№ 2. - С. 20 - 27.

4. Самарский, А. А. Численные методы [Текст] / А. А. Самарский, А. В. Гулин. - М.: Наука. -1989.