Stabilność algorytmu eliminacji dla szczególnych klas macierzy*Część I

And r ze j Kie ł ba s iń s k i

Warszawa

Stabilność algorytmu eliminacji dla szczególnych klas macierzy*

Część I

(Praca wpłynęła do Redakcji 1987.11.11)

0. Wprowadzenie. W licznych zadaniach numerycznych albo w algoryt- mach ich rozwiązywania spotykamy macierze specjalnych klas, np.:

— trójdiagonalne symetryczne,

— symetryczne, dodatnio lub ujemnie określone,

— diagonalnie dominujące.

Niekiedy należy rozwiązać układ równań liniowych z taką macierzą, niekiedy tylko skonstruować jej rozkład trójkątny.

Oczywiście, znane są świetne algorytmy wykonujące takie zadania, por.

np. [1], [2]. Są to warianty algorytmu eliminacji z wyborem elementu głównego lub bez wyboru. Są również dobrze znane własności numeryczne tych algorytmów, gdy stosujemy je do ogólnych macierzy. Znacznie trudniej znaleźć w literaturze adekwatny opis tych własności, gdy algorytm jest stosowany do macierzy jednej z wyżej wymienionych klas. Prowadzi to niekiedy do paradoksalnej sytuacji, że posługujemy się rozkładem trójkątnym macierzy A, wyznaczonym numerycznie w arytmetyce fl z precyzją względną v, mimo że nie jest spełniony typowy warunek

(0.1) v K (n )-cond(A) < 1,

gwarantujący konstruowalność rozkładu oraz jego przydatność do wyznacza- nia sensownych rozwiązań układu liniowego. (K(n) oznacza tu chaaraktery- stykę kumulacji błędów użytego algorytmu.)

Warunek (0.1) nie jest warunkiem koniecznym, a więc sporadyczny sukces takiego postępowania nie powinien nikogo dziwić, nawet jeśli wielkość v • K (n) • cond (A) jest nieco większa od jedności. Tymczasem nie trudno trafić na szerokie klasy zadań, w których wielkość v • K (n) • cond (A) jest bardzo duża, a mimo to konstruowany rozkład trójkątny „udaje się” systematycznie.

Warto więc wyjaśnić, jakie to specyficzne własności macierzy szczególnych klas pozwalają konstruować numerycznie i stosować rozkład trójkątny. Przy

* CPBP 01.01.04.

okazji wskażemy adekwatne charakterystyki uzyskanej dokładności rozkładu lub wyznaczonego rozwiązania układu liniowego. W analizie tej pomijamy zjawiska nadmiaru i niedomiaru arytmetyki numerycznej.

1. Macierze trójdiagonalne symetryczne dodatnio określone

1.1 Oznaczenia i podstawowe własności. Symetryczną macierz trójdiago- nalną będziemy w tym rozdziale oznaczali przez

hi a2 a2 b2 <*3 T = td(af; bt; ai + l) =

a n -

1 h„_i a„

an bn_

uzupełniając to oznaczenie w miarę potrzeby podaniem zakresu wskaźnika:

T = td (at; bt; ui + 1)r=1.

Załóżmy, że macierz T jest dodatnio określona (a więc ht- > 0, i = l,...,n ). Wszystkie opisywane w tym rozdziale własności dotyczą również macierzy ujemnie określonych (z ewentualnymi zmianami kierunku niektórych nierówności).

Obok macierzy T będziemy rozważali macierze |T| i T:

|T| = td(|af|; bt; \ai+l\) = B + |A|, (1.1) T = td (- |a f|; bt; - |

+ 1|) = B -|A |,

B = diag(bi), A = T -B = td ^ ; 0; ai + l).

Obie te macierze są ortogonalnie-diagonalnie-podobne do macierzy T,

|T| = SiTSi, T = S2TS2, ( 1 . 2 )

Sq = diag(S^, ..., = ±1, q = 1, 2, i = l, ..., n, mają więc te same wartości własne, co macierz T. Z (1.1) wynika ponadto równość

f-1 = B - 1/2( I -C )-1B- 1'2 = B 12 £ CkB-1/2,

k= 0

gdzie C = B~1/2-|A|-B~1/2. Otrzymujemy stąd zależności I T - ^ T " 1, |||T -1|||p < | | f - 1||p,

lim ilp = ||T||p, pe{l, 2,

F}.

Zależności, dotyczące którejkolwiek z tych czterech norm, będziemy zapisy-

wali w dalszym ciągu tej pracy bez podawania wskaźnika normy.

Macierz dodatnio określona T rozkłada się na iloczyn czynników trój- kątnych, przy czym czynniki te są dwudiagonalne:

(1.4) T = td (af; bt; a( + 1) = td (ef; 1; 0)-td(0; df; gf).

Elementy ef, df, gf są określone jednoznacznie. Równość (1.4) jest równo- ważna układowi równości:

(1.5) bv = df, af = ef dfL u bt = ef gf_ i+ df, a{ = gf_ t dla i = 2, ..., n.

Twi er dz eni e

1. Macierz symetryczna trójdiagonalna

jest dodatnio okreś- lona wtedy i tylko wtedy, gdy ma rozkład postaci (1.4), przy czym spełnione są nierówności

(1.6) df > 0, i = 1, ..., n.

Dowód. Niech będzie spełniona równość (1.4) (a więc i (1.5)). Możemy wówczas napisać

td (0; df; gf) = td(0; df; ai + J = diag(df)-td(0; 1; ef+1), a więc

T = L D LT, (1.7)

D = diag (df), L = td (ef; 1; 0).

Macierz L jest nieosobliwa, więc macierz T przystaje do macierzy D. Jest zatem dodatnio określona wtedy i tylko wtedy, gdy spełnione są warunki

( ¹ . ⁶ ). .

Warunek (1.6) jest efektywnym kryterium sprawdzania, czy symetrycz- na trójdiagonalna macierz T jest dodatnio określona. Jest ponadto istot- nym elementem następującego algorytmu wyznaczania wielkości {df; ef), i =

..., n, określających czynniki rozkładu (1.4) lub (1.7). (Wobec równości d* = blt gf =ai+l wielkości tych nie potrzeba specjalnie wyznaczać.)

Al g o r y tm

1. Dla danej macierzy/ T

td(a,; bt; ai

wyznaczamy wielkości {d,}"= l5 {et}"= 2 algorytmem

(1.8) dt : = bx; for i : = 2 to^ n do^ [e,: = a jd ^ i; di'.= br—ei *a,]. ■ Z zależności (1.5) przy założeniu (1.6) uzyskujemy dowód indukcyjny następującego twierdzenia.

Twi er dz eni e 2.

Algorytm

jest wykonalny w dokładnej arytmetyce dla dodatnio określonej macierzy T = td (a,; bt; ai + i), przy czym obliczone dt, c, są elementami czynników rozkładu

(1.9) di = df, £?; = ef. m

1.2. Teoria zaburzeń. Rozważamy małe względne zaburzenia <5T ele- mentów macierzy T = td(af; ai + 1),

(1.10) |<5T| ^ rj\T\, ÓT = ÓTt , rj^O.

Udowodnimy następujące twierdzenie:

Twier dz enie

3. Jeśli macierz T

td (a{

ai +

jest dodatnio określona, a macierz <5T spełnia warunki (1.10), to warunek

(1.11) rj ||T~1 • |T| || < 1,

(gdzie T = td (— łaj; b{; — |ai + 1|)) wystarcza, aby macierz T + <5T była dodatnio określona.

Dowód. Gdyby macierz T-f<5T nie była dodatnio określona, wtedy dla pewnej liczby 6, 0 < 6 ^ 1, macierz T + 0-<5T = T(I + T -1 -0-<5T) byłaby osobliwa. Jest to możliwe tylko wtedy, gdy zachodzi nierówność

(1.12) 1 ^ ||T~1 -0-<5T||.

Z (1.3) i (1.10) wynikają jednak nierówności

||T - 1 • 6 ■ <5T|| ^ ||T - 1 • ÓT\\ ^ ||T" 1 • \ST\ || ^ r] ||T“ 1 • |T| ||.

Uzyskujemy zatem sprzeczność (1.12) z założeniem (1.11). ■ U waga 1. Ze względu na zależności (1.3) zachodzi nierówność

||T~1 • |T| || ^ ||T~ 1|| • ||T|| = cond(T).

Zatem warunek (1.11) jest słabszy niż typowy warunek postaci:

ri • cond (T) < 1.

Przykładem macierzy, dla których zachodzi silna nierówność

||T~1 • |T| || <ścond(T),

mogą być tzw. „macierze ustopniowane”. Na przykład dla macierzy 104 -1 0 0

-1 0 102 -1

0 -1 1

mamy

T 1

989900

99 10 10

10 104 104

10 104 999900

Zachodzi więc

IIT-1 -|T| || w 3.3 <ścond(T) * 104. ■

Uwaga 2. Warunek (1.11) nie jest oczywiście konieczny, aby macierz T + <5T była dodatnio określona. Istnieją jednak macierze T rozważanej klasy i macierze zaburzeń <5T, spełniające warunek (1.10), dla których złamanie warunku (1.11) wiąże się z nieokreślonością macierzy T + <5T, np.

T = T(e) = 1 1- e

1- e 1 <5T = <$T(e)= — >7 |T (e)|.

Dla 0 < e < l zachodzi równość ||T“ 1 • |T| ||2 = (2 — e)/e. Zatem dla

£ > £ * := 2rj/(l + ri) warunek (1.11) jest spełniony. Dla e = e* warunek ten jest złamany. Równocześnie macierz T(e*) + <5T(e*) jest osobliwa ■

Udowodnimy jeszcze lemat, pozwalający na korzystanie z warunku (1.11) w dowodach indukcyjnych.

L

1. Rozważmy dwie macierze dodatnio określone (n > 1) T = td f a ; bi; at + OJL i , T' = td (ą ; bt; at + 0?= /.

Spełniają one nierówność

(1.13) ||(T')~1 • |Tj || < ||T_ 1 • |T| ||.

Dowód. Rozważmy dla ee(0, 1) macierze

T' 0

T(e) = ean

0 sa„ bn

Elementy obu macierzy

W(e) = (T (e)r\ |T (e)|

są nieujemne i nie maleją przy rosnącym e. Równocześnie zachodzi równość

W(0)|T(0)| (T')“ 1 • |Tj 0

0 1

Stąd

||( T r 1 • |Tj II < ||W(0) • |T(0)1 II < IIW(1)|T(1)1 II = ||T -1 • |T| II.

1.3. Analiza błędów zaokrągleń algorytmu rozkładu trójkątnego. Udowod- nimy następujące

Twier dz en ie 4.

Algorytm

realizowany w arytmetyce zmiennopozycyjnej (fi) ze względną precyzją v, v < 610 —3, jest wykonalny dla symetrycznej, dodatnio określonej, trójdiagonalnej macierzy T = td(a,; bt; ai + 1)eR n,n, spełnia- jącej warunek

(1.14) 0.76-v||(T)-1-|T|||p < l

dla dowolnej z norm p = 1, 2, oo. Obliczone wartości {d,}f=1, {e,-}”= 2 spełniają przy tym zależności

(1.15) di > 0, i = 1, . . n,

(1.16) f

td (dr, br, ai + l)

td (er, 1; 0)-td(0; dr, ai+1), gdzie

|a*l, Wil l£il < p:=0.76v,

|£,-| ^ p : = 0.26v, (i = 1, - n). m

Uwaga 1. Równość (1.16) wyraża rozkład macierzy f na iloczyn czynni- ków trójkątnych. Algorytm 1 daje więc „mało zaburzone czynniki rozkładu trójkątnego mało zaburzonej macierzy T”. Liczby p, p szacujące zaburzenia względne danych i wyników są mniejsze niż precyzja arytmetyki fi. ■

Uwaga 2. Dla wielu macierzy T wskazanej klasy algorytm 1 może być wykonalny w arytmetyce fi z zachowaniem zależności (1.15)-(1.17), mimo że warunek (1.14) nie jest spełniony. Wystarczy np., aby był spełniony warunek (1.15) a będą spełnione również warunki (1.16)-(1.17) (por. dowód twierdze- nia 4). Jednakże w pełnej klasie zadań dopiero warunek (1.14) gwarantuje wykonalność algorytmu, por. uwagę 2 do twierdzenia 3. ■

Dowód. Twierdzenie jest oczywiście spełnione dla n = 1. Przyjmijmy założenie indukcyjne, że twierdzenie jest prawdziwe dla n = k dla pewnego k > 1. Niech macierz T = td(u,; bt; ai + l)i=j spełnia warunek (1.14). Na mocy lematu 1 również podmacierz T' = td (af; br, ai + x)f=1 macierzy T spełnia ten warunek. Z założenia indukcyjnego algorytm 1 jest więc wykonalny dla macierzy T', czyli instrukcje (1.8) wyznaczają wielkości

spełniające warunki df = Oi( 1 + a,), ,1.17, * - M1 + W

e{ = e,-(l+ £,•), d{ = dj( 1 + 3()

(1.18) di> 0, i = 1, ..., k,

oraz rownosci

(1.19) bi — ij i 2, k,

bi = iidi + di,

przy zależnościach (1.17). Ze względu na (1.18) w arytmetyce fi jest wykonal- ny następny krok algorytmu 1 dla macierzy T:

ek+1 : = ak + iA4 ? dk+l : = bk +1 — ek + i *ak + 1.

Istnieją zatem takie liczby X,

,

że zachodzą zależności

ak +1 (1 + ^) = ek +1 dk, bk+l = ek

ak + ± (1 — @)4-dk + i (1 —

), W, lei <

|

| < v/(l-v) « v.

(

1

20

Wprowadźmy małe zaburzenia ak + 1, ek+1 równościami

(l+ a k + 1)2 = ( l+ i) ( l + <5k) ( l- e ) 1/4, |ak + i| ^ 0.76v,

(1 +£k +

)2 = (1 +*)-■1 (1 +Sk)~1 (1 -

)114, \sk + l\ ^ 0.76v, oraz, zgodnie z (1.17), wielkości

(1.22) ak+ i = ak+ x (1 + ak + x), ek + l = ek+l (1

1).

Z (1.17), (1.20), (1.21) wynikają równości

(1.23) dk + i — &k+1' dk,

(1.24) bk +1 = ck+1 dk+! (1 — g)3/4 + dk +1 (1 — t).

Drugą z nich przepiszmy w postaci

(1-25) &k +1 = $k+1 dk+1 + dk +1, gdzie

bk +1 — bk+i(l+P\ dk+i = dk + 1(l+<5),

(1.26) ( i +0) = (1 - e ) “3/4, (1+ Ą = (i -T)(1 - e ) - 3/4,

|/?|s$0.76v, |5| < 1.76v.

Zależności (1.19), (1.23), (1.25) oznaczają rozkład trójkątny macierzy

T = T

0 dk +1

dk+1 0

bk+1

spełniającej oszacowanie, por. (1.21) i (1.26),

|T —f | < 0.76v |T|.

Wobec (1.14) i twierdzenia 3 macierz T jest dodatnio określona. Zatem, por.

twierdzenie 1, spełniony jest warunek dk + 1 > 0, czyli — wobec (1.26)

(1.27) *4+i>0.

Pozwala to nam inaczej, niż w (1.26) określić zaburzenia wielkości bk+1 i dk+l. (Nie możemy przyjąć zaburzeń fik+l = fi, ók + l = S, gdyż S nie spełnia oszacowania (1.17).) Określamy element dk+1 oraz małe zaburzenia 3, ę, <5fc + 1 równościami

(1.28)

1 -3 = ( l- e ) 3/4, |3| ^ 0.755v, 1 — (p = (1 —

)3/4, \ę\ < 0.755v, l + <4 + i = ( 1-

)1/4, |Ąk + 1|<0.26vf

*4+i = dk+1 (i+<4+1)- Z (1.24) wynika wówczas równość

(1.29) bk+1 + 4 +14 +1 $+*4+1 = 4 +1 *4+1+*4+1 • Wprowadzamy w końcu zaburzenie fik + i-

(1.30) fik +1 = (4 +1 *4 +1 $ + *4 +1 (p)/bk +1.

Zgodnie z (1.17), wykorzystując oszacowania (1.28), zapisujemy (1.31) bk + i — bk+i (1 +fik + i), \fik + il ^ 0.76v, oraz, z (1.29),

(1-^2) bk+1 = 4 +1 *4 +1+*4 +1 •

Równości (1.19), (1.23), (1.32) oraz zależności (1.18), (1.27), oszacowania (1.17) dla i ^ k oraz (1.21) i (1.31) oznaczają spełnienie twierdzenia 4 dla n = k+ 1.

Tym samym uzyskaliśmy pełny dowód twierdzenia 4. ■

U waga 3. Zastępując założenie v ^ 610 —3 założeniem silniejszym, np.

v

^ 10-6, v < 1 0 -14 itd., możemy zastąpić stałe kumulacji 0.76, 0.26 z twierdzenia 1 stałymi bliższymi odpowiednio 0.75, 0.25. W dalszej części tej pracy będziemy pomijali w oszacowaniach udział składników, proporcjonal- nych do wyższych potęg v. ■

Konstrukcja zaburzeń (danych lub wyników), równoważnych błędom

zaokrągleń w algorytmie, nie jest oczywiście jednoznaczna. Różne takie

konstrukcje mogą być przydatne do różnych celów. W następującym twier-

dzeniu przedstawiamy przykładowo kilka alternatywnych sformułowań, cha- rakteryzujących jakość uzyskanych rozkładów (1.4) lub (1.7) macierzy T.

Twier dz en ie 5.

Przy założeniach twierdzenia

obliczone elementy

\dt J spełniają następujące zależności:

(1.33) T + . T = (L + ^ L,(R + JR) (rozkład LR),

|(L + dL)(D + dD)(L + dL)r (rozkład LDLl]

gdzie

L = td (et; 1; 0), R = td(0; dt; ai + l), D = diag (di), dT = td (ot cii I ^ Pi ;ai+1 o) +1),

(1-34) dL = td («;£,.; 0; 0),

dR = td(0; k k ; at + 1 ai+1), dD = diag(kk)-

Zachodzą przy tym następujące warianty oszacowań Rozkład LR:

(a) dT = dTr, k |, \Pi\, |£f| ^ 0.75v, |,5i\ ^ 0.25v.

(b) ai + i= 0 , k |, \pi\ ^ v, dL =

(c) dT = dTr, k l < v, \Pi\ ^ 2v, dL = 0 = dR.

Rozkład LDLT:

(d) JT = J T r, K| < v, |/J,| ^ 2v, JL = 0 = dR.

(e)

Dowód. Przypadek (a) odpowiada twierdzeniu 4, został więc już udowod- niony. Przepiszmy zależności (1.20) w postaci

a«(l + k) = eidi- 1, bi = e1ai( l - 0i) + ^(l-T ;), (1.35)

w . y , w ^ v.

Kładąc a ,• = ai + 1 =0, /łf = (e..^ & + */,-T;)/k, otrzymujemy przy- padek (b).

Kładąc af = af = A*, & = (et at Ą + dt

,)/^- , gdzie 1 - Ą = (1 - &)/( 1 + a,), otrzymujemy przypadek (c).

Przepiszmy (1.35) w postaci

a.(l+k) = eid,_1, ^ = e l d i - i ^ - ^ + diil-Ti), (1.36)

w , y>

< v,

odpowiadającej rozkładowi LDLr . Przyjmując zaburzenia takie, jak w dowo- dzie przypadku (c), otrzymujemy z (1.36) przypadek (d).

Konstrukcja zaburzeń dla przypadku (e) jest bardziej złożona. Połóżmy najpierw

et = £?; (1 + £;), di-i = di-i (1 + $ - 1), (1.37)

1+

= 1+ Ą -

= ( l+ ż tr 1/2, m, ^ 0.5v.

Z tymi oznaczeniami równości (1.36) przyjmują postać (1.38) cii = eidi- 1, bi = M 1 + /?,) = ef j + dit gdzie

Pi =(e}di-1rii + diyi)/bi,

(1.39) 1 -»/,• = (1 - ft)(l + ^ )1/2, 1 -7i = (1 - m + l t + i)112, hU ly.-l. \Pt\ < i.5v.

Równości (1.38) oznaczają rozkład typu LDL7:

(1.40) td(ą; bi; ai + l) = L D Lr, gdzie

L = td (et; 1; 0), D = diag (di).

Mnożąc obustronnie obie strony równości (1.40) przez macierz S =

= diag((1 +Pi)~1/4) otrzymujemy, zgodnie z (1.33), rozkład LDLT macierzy T + dT = S • td (a i ;bi;ai+i)- S, z czynnikami

L + dL = SLS" \ D + dD = SDS.

Zachodzą wówczas równości

1 + «, = (1 + P i ) ~ il4( 1 + A _ i)" 1/4, 1 + A = (1 + A)1/2, i + Ą = : ( i + A r

1

4

i+Ą = (i+ ^ ) ( i + A r i/2.

Łatwo sprawdzamy, że tak określone zaburzenia spełniają oszacowania, wskazane dla przypadku (e). ■

Uwaga. W dowodach twierdzeń 4 i 5 korzystaliśmy kilkakrotnie z tego,

że wyrażenia „ą a,” oraz „d,” są tego samego znaku. Pozwalało to nam

uzyskać dobre oszacowania konstruowanych zaburzeń elementu bt. Jeśli

algorytm 1 da się wykonać w arytmetyce fl dla symetrycznej macierzy

T = td (Oi't bil Oi + i), niekoniecznie określonej, to wyrażenia „ ą a ” oraz „d,”

mogą się różnić znakiem. Obliczone elementy {ej, {ftj spełniają jednak wówczas nadal zależności (1.35). Wprowadzając wielkości

ft = « .(!+ a,), ft = ft (1 + ft), ft = et (1 + £ź),

1 +s, = (1 +«,)/(1 + Aj, |a,|, |c,t $ 1.5v, Iftl v, otrzymujemy z (1.35) równości

f t= f t- f t- i, bi = efdi- 1+di.

Oznaczają one rozkład typu LDLr :

f = td(ft; ft; di + l) = L • D • Lr,

L = td(ft; 1; 0), D = diag(ft).

Rozkład ten może być wykorzystany do obliczenia wyznacznika macierzy T lub do określenia liczby ujemnych wartości własnych tej macierzy. Nie może być jednak na ogół wykorzystany do rozwiązania układu Tx = f, gdyż mogłaby w takim obliczeniu wystąpić poważna utrata dokładności, jeśli (różniące się znakiem) wyrażenia „ft”, „e.a,” są co do modułu znacznie większe od odpowiedniego elementu ft. Taki algorytm wyznaczania rozwią- zania x (przy nieokreślonej macierzy T) jest na ogół niestabilny.

Zauważmy jeszcze, źe macierz T jest bliska macierzy T,

|f — T| < 1.5v|T|,

ale iloczyn L • D • Lr, gdzie L = td (e{; 1; 0), na ogół nie przybliża T, nawet w normie, tzn. wyrażenie ||L * D -Lr — T||/||T|| nie może być ograniczone z góry w klasie wszystkich symetrycznych macierzy T (dla których jest wykonalny algorytm 1). ■

1.4. Rozwiązywanie układu równań liniowych. Rozkład trójkątny (1.4) poz- wala rozwiązać układ równań liniowych

(1.41) Tx = f

dla dowolnego wektora f = (fi)eR n, przez rozwiązanie dwu układów o macierzach dwudiagonalnych

Lh = f, Rx = h, gdzie

L = td(e{; 1; 0), R = td (0; ft; ai + i).

W istoHe wielkości eit d{, ft grają tu tylko rolę pomocniczą. Interesuje nas jedynie przejście od danych początkowych {af, ft,^} do wyników |x j. Roz- ważmy więc algorytm, realizujący to zadanie.

8 — Matematyka Stosowana t. 31

Al g o r y t m

2. Dla danej macierzy dodatnio określonej T =

= td(flj; bil ai+1)i=i i danego wektora f= (fi)e R n wyznaczamy składowe {x.-}r=i wektora x e R n, będącego rozwiązaniem układu Tx = f:

K t : (rozkład LR macierzy T i rozwiązanie układu Lh = f):

di := bx; ht := fx

for i 2 to n do et : = ajd{ _!; dt : = bt: - et * at;

hi : = f - e i *hi. l

K 2: (rozwiązanie układu Rx = h): xn:=hn/dn;

for i : = n — 1 downto 1 do^ x{: = (/i, — ai + 1* xt + . ■ Udowodnimy następujące twierdzenie.

Twier dz enie 6.

Algorytm

realizowany

arytmetyce

ze względną precyzją v jest wykonalny dla dodatnio określonej macierzy T =

= td(a{; bt\ ai + l) eRn,n, spełniającej warunek (1.14), oraz dowolnego wektora f eRn. Obliczone rozwiązanie ?T=(xf) spełnia układ równań

(1.42)

(T + <5T)(x + <5x) = t + 6 f , przy czym zachodzą nierówności

(1.43)

|<5T|^4.5v|T|, |<5f|^|f|-1.75v,

\óx\ |x| -2.75v. ■

Dowód. Błędy zaokrągleń w algorytmie 2 wyrazimy równościami

0,(1 + ^) = bi = Ąf l i ( l - f t ) + di(l-T ,-),

f ( 1 + (pi) — ht A et /if'- i (1 + ot) (1 + (pt),

diXi(l

£i)

1xi

3i),

w, ifti, itii,

kii,

v,

2v.

Zdefiniujemy następujące wielkości zaburzone:

Xi

161

2v,

at

i-,

fi = f ( l + (Pi)> l<Pil

et = ei(l+£i),

N

2v, di

At)

3v,

A = + IAI ^ 6V,

(1.45)

gdzie 1 -rji = (1 -&)(1 +£,•)“ 1 (1 + £i) Wykorzystując te oznaczenia, przepi- szemy (1.44) w postaci

Cli = €[ df _ j , bi €i Cli "I” j

(1.46) fi = hi + eihi-i,

di x, hi u,- + j Xj +1.

Równości te oznaczają, że wektor x = (x.) jest rozwiązaniem układu równań

(1.47) td(a£; Ą+1)x = f,

gdzie f = (fi). Równość tę przekształcimy do postaci (1.42), kładąc T 4- <5T = S • td (a,-; bt; a,- + i) • S,

(1.48) x"+<5xf = S_ 1 * x, i + ó f = S f, gdzie

S = d iag ((l+ /y -1,s).

Łatwo sprawdzamy, że zdefiniowane w (1.48) zaburzenia spełniają oszacowa- nia (1.43). ■

Uwaga. Wykorzystanie rozkładu (1.7) macierzy T do rozwiązania układu (1.41) sprowadza się do pewnej modyfikacji algorytmu 2, np. przez następują- cą zmianę kroku K 2:

K'2: xn:= h jd n; for i : = n - 1 downto 1 do^ 1*:= hi/di — ei+ lxi+1'].

Odpowiadający temu algorytmowi wariant twierdzenia 6 różni się jedynie oszacowaniami zaburzeń. Zamiast (1.43) mielibyśmy

|ÓT|^5.5v|T|, |<5f| ^ |f| • 2v, (1.43')

|£x| < |x|*3.5v. ■

1.5. Wnioski. Twierdzenia z paragrafów 1.3 i 1.4 wykazują, że rozkład trójkątny macierzy symetrycznej, dodatnio lub ujemnie określonej, trójdiago- nalnej T może być wykonany w arytmetyce fi (z dostatecznie wysoką precyzją) algorytmem eliminacji w sposób numerycznie poprawny.

To samo dotyczy rozwiązań układów liniowych z macierzą T.

Skonstruowane zaburzenia danych i wyników, równoważne błędom zaok- rągleń, są małe względem indywidualnych elementów. Kumulacja błędów nie przekracza kilku jednostek, nie zależy od wymiaru zadania.

Algorytm eliminacji zapewnia więc dla zadań tego typu zasadniczo

najwyższą możliwą jakość obliczeń.

1.6. Algorytm rozkładu trójkątnego ze skierowanym zaokrąglaniem.

W paragrafach 1.3-1.5 rozważaliśmy numeryczną realizację algorytmów w standardowej arytmetyce fl. Zakładaliśmy, że jest to arytmetyka z popraw- nym zaokrąglaniem wyniku każdego działania, tzn. dla dowolnych liczb a, b tej arytmetyki, dla których dokładny wynik x = a 0 b działania O (0 e {+,

—, *, /}) jest dobrze reprezentowalny w fl, wynik obliczony c, c : =a Ob , spełnia równość

(1.49) c = rd(x) = x -(1 +e), |e| ^ v.

Wytworzony błąd względny e = e{x) może być przy tym dodatni lub ujemny. Operację rd nazywamy zaokrągleniem. Jeśli przyporządkowuje ona liczbie x najbliższą liczbę arytmetyki fl, to nazywamy ją zaokrągleniem symetrycznym.

Zwykłe zaokrąglenie rd może być zastąpione, por. [3]:

— zaokrągleniem „skierowanym w górę”

(1.50) c = rdg(x) = x-(l +e), 0 ^ £ ^ v, lub

— zaokrągleniem „skierowanym w dół”

(1.51) c = rdd(x) = x*(l + £), - v ^ e ^ O .

Precyzja v zaokrągleń skierowanych w przypadku realizacji sprzętowej wiąże się z precyzją v zaokrąglania symetrycznego (1.49) zależnością v = 2v, a przy realizacji programistycznej zazwyczaj zależnością v = 4v (koszt typowej realizacji programistycznej zaokrąglenia skierowanego jest równy kosztowi jednego dodatkowego mnożenia w fl).

Zaokrąglenia skierowane mogą być wykorzystane w następujących dwu wersjach rozkładu trójkątnego macierzy trójdiagonalnej, symetrycznej, dodat- nio określonej T = td (ą ; bt; a{ + i)^ x.

Al g o r y t m lg.

dt : = h i; for i : = 2 to_ n do et rdd {ajd^i) wt : = rdd (et * a,)

di:

rdg

—

dx : = b i; for i: — 2 to_ n do ei : = rdg(ai/di_1) wi : = rdg(ei *ai)

di: = rdd (h,- — w,)

Udowodnimy następujące twierdzenia:

Twi er dz en ie 4g.

Algorytm

jest wykonalny dla każdej macierzy dodatnio określonej T = td (a,; ai + 1)"=i- Obliczone elementy {ej, {dj spełniają zale- żności

(1.52) di > O, i = 1, n,

(1.53) td(a,(l + af); M l+ ft); ai + 1(l+ a i + 1)) =

= tdfejl+e,); 1; 0)td(0; d jl + <5f); af+i (1 + ai + J), przy czym zaburzenia względne a,-, £,-, <5f spełniają nierówności

— 0.75v ^ a* < O, O < ^ ^ 0.75v, (1.54)

— 0.125v ^ £,• < 0.625v, -0.25v < Ą < 0. ■

(i) Algorytm ld jest wykonalny dla dodatnio określonej macierzy T = td (a{; £>,•; a, + 1)"=1, spełniającej warunek

(1.55) 0.75v||T_1 • |T| || < 1.

Obliczone elementy {ej, JdJ spełniają wówczas zależności (1.52), (1.53), przy czym zaburzenia względne cii, /?*, £f, <5f spełniają nierówności

O ^ a, < 0.75v,

— 0.75v ^ ^ ^ O, (1.56)

-0.625v^ą- <0.125v, O ^ Si ^ 0.25v.

(ii) Jeśli macierz T nie jest dodatnio określona, to algorytm ld wyznaczy dla pewnego i (i = 1, ..., n) wartość dt < 0. ■

Dowód twierdzenia 4g. Zaokrąglenia skierowane w algorytmie lg powodują, że spełnione są nierówności

|e,j < h M - J ,

(1.57) W i^ei-ai,

di > bi-Wi.

Prosty dowód indukcyjny wykazuje prawdziwość nierówności (1.58) k | < |eH, di>dT (i = 1, ..., n),

gdzie ef, df są elementami czynników rozkładu trójkątnego (1.4) macierzy T

(por. tw. 2). Zatem dla każdego i zachodzi d, > 0, więc algorytm lg jest wykonalny. Zależności (1.53), (1.54) otrzymujemy podobnie, jak w dowodzie twierdzenia 4. Wystarczy tylko zauważyć, że błędy w (1.20) spełniają (w przypadku algorytmu lg) nierówności

(1.59) — v ^ ^ 0, 0 < £ < v, 0 <

^

. Przy założeniu indukcyjnym

(1.60) — 0.25v < <$k ^ 0

otrzymujemy z (1.21), (1.28), (1.30) oszacowania (1.54). ■

Dowód twierdzenia 4d. Część (i) wynika z twierdzenia 4 oraz podob- nego szacowania konstruowanych zaburzeń, jak w dowodzie twierdzenia 4g.

Część (ii) wynika ze spostrzeżenia, że wyznaczone algorytmem ld elementy eit ^ muszą spełniać nierówności

kil > 131, di < di,

gdzie Ją}, {df} oznaczają wielkości, które byłyby wyznaczone algorytmem 1 w dokładnej arytmetyce. W przypadku macierzy T nie określonej dodatnio, dla pewnego i musiałaby zachodzić nierówność dt ^ 0. ■

Uwaga. Algorytmy lg czy ld są oczywiście numerycznie poprawne, w tym samym sensie co algorytm 1, ale z charakterystyką kumulacji dwukrot- nie (lub czterokrotnie) większą. Wybór jednego z tych algorytmów, zamiast algorytmu 1, zależy od konkretnego zadania, któremu konstruowany rozkład trójkątny macierzy T miałby służyć. ■

Prace cytowane *

[1] B. T. Sm ith, J. M. Boyle, Y. Ikebe, V. C. K lem a, C. B. M oler, Matrix Eigensystem Routines EISPACK Guide, Springer Verlag, New York 1970.

[2] J. D o n g arra , I. R. Bunch, C. B. M oler, G. W. S tew art, U N PA C K Users Guide, SIAM Publ. Philadelphia 1978.

[3] W. I. C ody i in., A proposed Radix- and World-length independent Standard for Floating- point Arithmetic, IEEE Magazine Micro, Aug. 1984, 86-100.