Численный метод нахождения алгебраического решения ИСЛАУ, основанный на треугольном расщеплении Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 4, № 4, 1999

ЧИСЛЕННЫЙ МЕТОД НАХОЖДЕНИЯ АЛГЕБРАИЧЕСКОГО РЕШЕНИЯ ИСЛАУ, ОСНОВАННЫЙ НА ТРЕУГОЛЬНОМ РАСЩЕПЛЕНИИ

А. Ю. Карлюк

Новосибирский государственный университет, Россия

The paper investigates stationary single-step iteration methods for computing algebraic solutions to interval linear systems, which are based on the triangular decomposition of the matrix of this system. We derive conditions sufficient for the convergence of successive approximations as well as an estimate of the convergence speed.

Введение

Основным объектом нашего исследования является интервальная система линейных алгебраических уравнений (ИСЛАУ) в полной интервальной арифметике Каухера IR:

где С = (Су) и в = (ф).

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

Задача нахождения алгебраического решения возникает при использовании алгебраического подхода к оцениванию различных множеств решений интервальных линейных систем, состоящего в замене исходной задачи интервального оценивания на задачу нахождения алгебраического решения некоторой вспомогательной ИСЛАУ в полной интервальной арифметике — так называемой системы в дуализациях [6, 7].

С. П. Шарым в [6, 8] были предложены два общих подхода к построению стационарных одношаговых итерационных методов для вычисления алгебраического решения системы (1). В первом существенно используется погружение пространства интервальных

© А.Ю. Карлюк, 1999.

<

c11x1 + c12x2 + • • • + c1 nxn — di ,

c21x1 + c22x2 + • • • + c2nxn = d2 , ........

. . . Cj, x,, d, G Ж-

Сп1х1 + Сп2х2 + • • • + Сппхп )

V

Для краткости будем также пользоваться матричной формой записи:

Сх = а,

(1)

векторов Жга в линейное пространство в два раза большей размерности М2га. Во втором подходе предлагается привести исходную систему к виду

х = Тх, (2)

Т : Жга ^ Мга и, выбрав некоторое начальное приближение х(0) € Жга, реализовать итерационный процесс

х(к+1) = Тх(к) (3)

непосредственно в интервальном пространстве Жга.

Как известно, в интервальных арифметиках, как классической, так и полной арифметике Каухера, не выполняется свойство дистрибутивности умножения относительно сложения: в общем случае ах + Ьх = (а + Ь) х. В связи с этим приведение системы (1) к

виду (2), в отличие от вещественного случая, является не вполне тривиальной задачей,

так как в (2) присутствуют два вхождения переменной х, а в (1) — только одно.

Для приведения системы к указанному виду предлагается найти дистрибутивное расщепление матрицы С [8], т. е. представление С в виде суммы

С = О + н такое, что Сх = (С + Н)х = Ох + Нх.

Если для оператора умножения на матрицу О можно построить обратный оператор Т,то из уравнения

Ох = а © Нх,

эквивалентного (1), получим

х = Т(в © Нх),

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

x

(k+1) —

T(d е Hx(ll)). (4)

Подобные алгоритмы дают существенно более низкую скорость сходимости, чем, например, субдифференциальный метод Ньютона [8, 12], основанный на погружении в линейное пространство. К достоинству же их относится то, что они могут гарантировать единственность решения, чего нельзя достичь при использовании субдифференциального метода Ньютона.

1. Постановка задачи

Интервалы-элементы Ж будем обозначать жирными буквами латинского алфавита. В случаях, когда имеет смысл указать концы интервала, будем пользоваться обозначением [x, x]. Детальное описание полной интервальной арифметики (арифметики Каухера) можно найти, например, в [9].

Напомним, что правильной проекцией интервала называется интервал

f [x, x], если x < x,

pro x =

[ [x, x], иначе.

В Ж для любого интервала существует противоположный:

opp x — [—x, —x], и определена операция внутреннего вычитания:

x е y — x + opp y-

Модулем интервала называется величина

|x| — max{ |x|, |x| }• Расстояние между интервалами определяется по формуле

q (x y) — max{|x — y1, |x — y|} Близостью интервала к нулю называется величина

min{ |a|, |a| }, если pro a ^ 0, 0, иначе •

(5)

Приведем также несколько соотношений, имеющих место в Ж, которые нам понадобятся в дальнейшем [9]:

Целью данной работы является исследование сходимости итерационного процесса (4), получаемого в результате треугольного расщепления матрицы C, когда G и H берутся в виде нижней и верхней треугольных интервальных матриц, причем на главной диагонали матрицы G — (gij) стоят обратимые элементы (pro g,, ^ 0), а у матрицы H главная диагональ состоит из нулевых элементов.

В дальнейшем для краткости вместо записи (4) будем использовать также (3), считая

2. Псевдометрика и теорема Шредера

На пространстве интервальных векторов Жга можно ввести псевдорасстояние [1, 3]:

каждая компонента которого является расстоянием между соответствующими компонентами интервальных векторов в пространстве Ж и определяется формулой (5):

Для изучения итерационного процесса (4) воспользуемся общей теоремой Шредера

о неподвижной точке в псевдометрическом пространстве. Формулировку этой теоремы можно найти в [3]. Для наших целей достаточно следующего следствия.

Если существует линейный непрерывный положительный оператор Р : Мга ^ Мга такой, что

q (ab, ac) < |a| ■ q (b c) q (a + b, c + d) < q (a, c) + q (b, d), q (opp a, opp b) — q (a, b)

Tx — T(d е Hx)

q : IRn x IRn ^ R+,

q2

( q(xb У1) \ q(^ У2)

q(x y)

V qn / \ q(xn,yn) )

q (Tx, Ty) < Pq (x, y),

(6)

и ряд

ГО

£Р7 *

7=0

сходится при любом а € К, то итерационный процесс (3) сходится при любом начальном приближении х(0) € Жга к единственной неподвижной точке х*, причем имеет место оценка

д(х*,х(к)) < а - ай, где а = Иш ак, а ак определяются рекуррентно:

{ ао = 0 ,

1 ак+1 = Рак + д (х(0), х(1)) .

3. Нахождение оператора Липшица

Для того чтобы воспользоваться теоремой Шредера, нужно найти оператор Р : Мп ^ Мп, удовлетворяющий свойству (6).

Распишем выражение для т (а 0 Нх) покомпонентно, обозначив для краткости компоненты вектора (а 0 Нх) через £*:

( а1 0 ЕП=2 х \

т (а 0 Нх) = т

а2 0 Е?=3 Ь2І х

аз 0 ЕТ=4 Ьзі х

а*0 Е п=і+і Ьіі х

вГО—1 0 hn-1.ro х

в

т

( Г1 \ Г2 Гз

Гп-1

V п /

/ ёп ■ Г1

ё22 0 ё21 ёП Ъ

ёЗз1 0 ё31 ёи Ъ 0 ё32 ё22 Г2

ё44 0 ё41 ёи Ъ 0 ё42 ё22 ^ 0 ё43 ё33 ^

ёпП 0 ёп1 ёп1 Ґ1 0 ёп2 ё221 Г2 0 ёп3 ё331 Г3 0 ••• 0 ёп,п-1 ёп-1,п-1 Гп-1

Обозначим компоненты полученного вектора через к'

к

( к1 \ к2

у кп у

т (а 0 Нх).

Попытаемся покомпонентно оценить величину д(Тх, Тх'), х, х' Є

д(Тх, Тх')

( д1(Тх,Тх/) \ д2(Тх, Тх') ^ д„(Тх, Тх') !

д1(Тх Тх') = д(к^ к1) = д(ёц^ъ ёи^) < 1 ё і /1 ■ д(Г1, =

п п \ / п \ / п 4

кйіЧЕЧ-х-^Ч-х'- < кйіЧ £д(Чх-,Ч-х'-)) < ШЕ1 ■ д(х-,х'-)

чі=2 -=2

\І=2

чІ=2

д2 (Тх Тх') = д(^ к2) = д( ё221(Г2 0 ё21к1) , ё221(Г'2 0 ё21к^ ) <

< ^ д (Ъ 0 ё21 к1, ї'2 0 ё21 к'1) — |ё221| д ї + ОРР(ё21^), ї'2 + Орр (ё21к1_)) <

< |ё221^ д(^ Г'2)+ д( ОРР (ё21к1) , ОРР(ё21к1^ =

|ё221М д £ Ч-х- ^ h2jх'- + д(ё21кЪ ё21к1)

<

-=3 -=3

<!ё22^ —3 ^ 1 ■ д(х- , х'-) + |ё21| ■ д(кЪ кі^ = |ё221^£ ^2- 1 ■ д(х-, х'-) + |ё21| ■ д1(Тх Тх')

Производя аналогичные преобразования, далее получаем:

д3(Тх Тх') = д(k3, к3) = д( ё331(Г3 0 ё31 к1 0 ё32к2) , ё331(Г'3 0 ё31к1 0 ё32к2)) <

< |ё3311 ( X! 1 ^ д(х- , х'-) + |ё31| ■ д1 (Т^ Тх') + |ё32| ■ д2(Тх Тх')

\І=4

Тх') = д^ к') <

п

<|ё''11 ( X! |^д(х- , х'- Жёг^&ОГ^ + Тх') + . •• + 1 ё','-1^'- 1(Tx, Тх0

Ч-='+1

дп-і(Тх, Тх') = д(кп-1, кП-1) <

< |ёп-1,п-11 (|hn- 1,п1 ■ д(хпх'п) + |ёп-1,11 ■ д1(Тх Тх')+

+ |ёп -1,21 ■ д2(Тх, Тх') + • • • + |ёп- 1,п - 2 | ■ дп - 2 (Тх, Тх')),

п

п

п

п

п

дп(Тх, Тх') = д(кп, кп) <

< |ё-,п1 ( |ёп,1| ■ д1(Тх Тх') + 1 ёп,2 1 ■ д2(Тх Тх') + • • • + |ёп,п-1| ■ С/п-ОВ терминах матрицы С полученные соотношения для д- (Тх, Тх') будут выглядеть сле-

дующим образом:

('-1 п \

£ |с-1 ■ д(Tx,Тх') + £ !су1 ■ д(х-,х'-) I , і = ^ • •,п (7)

-=1 -='+1 /

Введем в рассмотрение следующие вещественные п х п-матрицы:

Б = diag( а1, • • •, ап) с элементами а- = Іс- , (8)

Г І С-1, если і > ^',

Ь = (/-) с элементами = < (9)

[ 0, если і < ^',

Г 0, если і > ^’,

Я = (г-) с элементами г- = < (10)

[ | С'-1, если і<^\

Тогда, с учетом новых обозначений, полученные соотношения (7) можно переписать в матричном виде:

д (Тх, Тх') < Б(Ь ■ д (Тх, Тх') + Я ■ д (х, х'))

Для краткости вектор д (х, х') будем обозначать символом р, а вектор д (Тх, Тх') — символом д . Оценим сверху правую часть последнего соотношения, воспользовавшись им же самим:

д < Б Яр + БЬд < БЯр + БЬ ■ Б(Яр + Ьд ) = БЯр + БЬБЯр + (БЬ)2д <

< БЯр + БЬБЯр + (БЬ)2Б(Яр + Ьд) = Б Яр + БЬБЯр + (БЬ)2БЯр + (БЬ)3д < • • •

• • • < БЯр + БЬБЯр + (БЬ)2БЯр + (БЬ)3БЯр + • • • + (БЬ)п-1БЯр + (БЬ)пд < • • •

Матрица БЬ — нижняя треугольная и имеет на главной диагонали нули, т. е. нильпотент-ная, (БЬ)п = 0, значит

д < (Ъя + БЬБЯ + (БЬ)2БЯ + (БЬ)3БЯ + • • • + (БЬ)п-1Бя)р •

Искомая оценка имеет вид

д (Тх, Тх') < Рд (х, х'),

Р = БЯ + БЬБЯ + (БЬ)2БЯ + (БЬ)3БЯ + • • • + (БЬ)п-1 БЯ^ (11)

4. Условия сходимости

Как известно, для выполнения условий теоремы Шредера достаточно, чтобы ряд

ГО

* *

7=0

сходился при любом начальном приближении а £ [3]. В свою очередь, необходимым и

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

Заметим, что собственные числа матрицы РР нулевые, и для нее имеет место теорема Неймана о разложении обратной матрицы в ряд

где I виде

I + РР + (РР)2 + ... + (РР)п—1 = £(РР)7 = (I - РР)-1,

7=0

единичная матрица. С учетом этого соотношения формула (11) перепишется в

V1/

р = (I - ррррр

Введем обозначения:

/1\ /вЛ

у вп у

^ ' р* , р* — элементы матрицы Р. 7=1

Очевидно,

й = Ре = (I - РР)-1РРе, (I - РР)й = РРе.

Матрица (I — РР) — нижнетреугольная: над главной диагональю у нее стоят нули, на главной диагонали — единицы, а под главной диагональю — элементы, противоположные по знаку соответствующим элементам матрицы РР. Учитывая, что (РР)* = аг/у, (РЛ)г7- = агГу, получим

/ П

«1 = Е а 1 г 1* , 7=2

П

«2 = а2^21 ' «1 + Е а2г27 ,

7=3

\ г—1 п

«г = Е аг1г7 ' «7 + Е аггг7 ,

7=1 7=г+1

П—1

вп /Ж • ві

І=1

го

п

Є

в

или, в терминах матрицы C,

г—1

(г—1 n N

icij I sj + i j=1 j = i+1 /

' 11 ■ I 2^ 1^7 *7 -I-

7=г

Если зг < 1 при любом г £ {1,... , п}, то (так как шах{зг} — норма оператора Р) выпол-

няется требуемое ограничение р (Р) < 1 [4].

Нетрудно понять, что выполнения неравенств

|cii 1|^|cij1 < 1> * £{l,---,n} (12)

j=i

достаточно для справедливости соотношений s < 1. Заметим, что

|cii 1| = I [ 1/cii, 1/cii] I = max{ 11/cii|, |1/c-| }

min{|cii|, |cii|} (cii)

С учетом последнего (12) перепишется в виде

(cii) > |cj|-

j = i

Это есть условие диагонального преобладания для интервальной матрицы [11].

5. Оценка сходимости

Теорема Шредера [3] дает следующую оценку сходимости стационарного итерационного метода:

q (x*, x(k)j < а — afc ,

где x* — неподвижная точка оператора T, а а = lim ак — предел последовательности,

k^-ro

определяемой рекуррентным правилом

ао = 0,

afc+1 = Pafc + q (x(0), x(1)) .

Мы можем адаптировать эти результаты на исследуемый интервальный итерационный процесс. Преобразуем выражения для ак:

а1 = q (x(0), x(1)) , а2 = Pa1 + а1 ,

аз = Ра2 + а1 = P (Ра1 + а1) + а1 = Р2а1 + Ра1 + а1 ,

k—1

ак = Рак—1 + а1 = Pk—1а1 + Pk—2а1 + ... + Ра1 + а1 = Е Pjа1 ,

j=0

ГО

а = lim ак = Pj а1 = (I — P) 1а1.

k^ro v '

j=0

г

1

1

Окончательно получаем оценку:

й—1

д (V, х(к)) < ^ (I - Р) 1 - £ Р^ Я (х(0),х(1)) .

6. Результаты

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

Сх = а (1)

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

п—1

Р = £(РР)7РР = (I - РР)—1РР

7=0

меньше единицы (Р, Р, Р определяются формулами (8) - (10)). При этом имеет место оценка

д (V, х(к)) < ^(! - Р) —1 - 7^ Р7^ д (х(0), х(1)) .

Для выполнения условия р (Р) < 1 достаточно выполнения следующих ограничений на матрицу С системы:

1 / г— 1 п \

«г = 1--Г I £ |Сг7 | «7 + £ |Сг7 И < 1 для всех г £ {1, . . . , п}.

(Сгг) \ 7=1 7=г+1 /

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

строгого диагонального преобладания:

(Сгг) > £ |Сг7 | для всех г £ {1, . . . , п}.

7=г

Список литературы

[1] Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. Мир, М., 1987.

[2] Калмыков С. А., Шокин Ю.И., Юлдашев З. Х. Методы интервального анализа. Наука, Новосибирск, 1986.

[3] Коллатц Л. Функциональный анализ и вычислительная математика. Мир, М., 1969.

[4] Ланкастер П. Теория матриц. Наука, М., 1982.

[5] ОРТЕГА Дж., РЕйНБОЛДТ В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. Мир, М., 1975.

[6] ШАРЫй С. П. Численное нахождение алгебраического решения интервальных линейных систем. В “Дискретная математика”. КГТУ, Красноярск, 1996, 126-142.

[7] ШАРЫй С. П. Алгебраический подход к анализу линейных статических систем с интервальной неопределенностью. Изв. РАН. Теория и системы управления, №3, 1997, 51-61.

[8] ШАРЫй С. П. Численное нахождение алгебраического решения интервальных линейных систем. Вычисл. технологии, 4, №3, 1999, 85-102.

[9] KAUCHER E. Interval analysis in the extended interval space IR. Comp. Suppl., 2, 1980, 33-49.

[10] KEARFOTT R. B. Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht, 1996.

[11] NEUMAIER A. Interval Methods for Systems of Equations. Cambridge University Press, 1990.

[12] SHARY S. P. Algebraic approach to the interval linear static identification, tolerance, and control problems, or One more application of Kaucher arithmetic. Reliable Comp., 2, No. 1, 1996, 3-33.

Поступила в редакцию 29 января 1999 г.