7. Układy równań liniowych
Wtej części rozważymy zagadnienie rozwiązywania zadań ogólnej postaci:
8
a11x1 + a12x2 + ::: + a1nxn = b1
>
>
<
a21x1 + a22x2 + ::: + a2nxn = b2
(115)
>
>
:
an1x1 + an2x2 + ::: + annxn = bn
zwanych w dalszej części wykładu URL, czyli układów n równań liniowych z n niewiadomymi o współczynnikach
aij i bi rzeczywistych. Odpowiednio oznaczając macierz współczynników układu oraz wektory wyrazów wolnych i
niewiadomych możemy krótko zapisać układ w postaci:
Ax = b:
Najważniejszym faktem, który wykorzystamy przy konstruowaniu metod rozwiązywania układów postaci (115) jest
równoważność układów. Przypomnijmy, że dwa URL są równoważne, jeśli jeden z nich jest otrzymany z drugiego za
pomocą skónczonego ciagu następujących operacji elementarnych: zamiana dwóch równań miejscami, pomnożenie
równania przez liczbę różną od zera oraz dodanie do równania iloczynu innego równania przez liczbe różną od zera.
7.1. Układy równań liniowych z macierzami o specjalnej strukturze
Rozpoczniemy od metod rozwiązywania układów równań, które mogą bżc łatwo rozwiazane, tj. takich, których
macierz współczynników ma pewna uproszczoną strukturę.
Najpierw, przypuścmy, że macierz A jest diagonalna, tzn. wszystkie niezerowe elementy macierzy leżą na głównej
przekątnej:
2 3 2 3 2 3
a11 0 ::: 0 x1 b1
6 7 6
0 a22 ::: 0 x2 7 6 b2 7
6 7 6 7 6 7
= (116)
4 5 4 5 4 5
::: ::: ::: ::: ::: :::
0 0 ::: ann xn bn
W tym przypadku rozwiązanie opisujemy krótko za pomocą wzoru:
bi
xi = ; i =1; 2; :::; n: (117)
aii
Oczywiście, układ (116) można numerycznie rozwiązać, o ile aii =0 dla wszystkich i.
6
Pzypuścmy teraz, że macierz A jest trójkątna dolna, tzn. jej wszystkie elementy niezerowe leżąnai poniżej głównej
przekątnej:
2 3 2 3 2 3
a11 0 ::: 0 x1 b1
6 7 6
a21 a22 ::: 0 x2 7 6 b2 7
6 7 6 7 6 7
= (118)
4 5 4 5 4 5
::: ::: ::: ::: ::: :::
an1 an2 ::: ann xn bn
Aby rozwiązać powŹzszy układ, założymy znów, że aii = 0 dla wszystkich i. Formalny algorytm nazywany jest
6
postępowaniem wprzód i jest nastepującej postaci:
PiĄ1
bi Ą aijxj
j=1
xi = ; i =1; 2; :::; n: (119)
aii
W podobny sposób rozwiązuje się układy z macierza trójkątną górną:
2 3 2 3 2 3
a11 a12 ::: a1n x1 b1
6
0 a22 ::: a2n 7 6 x2 7 = 6 b2 7
6 7 6 7 6 7
(120)
4 5 4 5 4 5
::: ::: ::: ::: ::: :::
0 0 ::: ann xn bn
Tym razem, formalny algorytm jest nazywany postępowaniem wstecz:
Pn
bi Ą aijxj
j=i+1
xi = ; i = n; n Ą 1; :::; 1: (121)
aii
UWAGA: Za pomocą powyższych metod można rozwiązywać również inne układy równań, o ile dadzą się trans-
formować do postaci (116), (118) lub (120) za pomocą permutacji równań.
1
7.2. Rozkłady typu LU
Przypuśćmy teraz, że macierz A układu (115) może byćrozłożona na iloczyn macierzy trójkątnej dolnej L i trójkąt-
nej górnej U, tzn. A = LU. Wówczas, aby rozwiązać układ Ax = b, wystarcza rozwiązać dwa układy z macierzami o
specjalnej strukturze:
Lz = b względem z
Ux = z względem x.
Pozostaje zatem tylko znalez
ć algorytm rozkładu macierzy A na iloczyn dwóch macierzy
2 3 2 3
l11 0 ::: 0 u11 u12 ::: u1n
6 7 6
l21 l22 ::: 0 0 u22 ::: u2n 7
6 7 6 7
L = oraz U = : (122)
4 5 4 5
::: ::: ::: ::: ::: ::: ::: :::
ln1 ln2 ::: lnn 0 0 ::: unn
Okazuje się, że L i U nie są jednoznacznie określone za pomocą równania A = LU, tzn. dlakażdego i możemy znalez
ć
niezerowąwartość albo lii, albo uii (ale nie obu). Wyróżnia się dwa szczególne przypadki rozkładu LU, a mianowicie:
(i) lii =1 dla i =1; 2; :::; n, wówczas L będzie macierzą trójkątną dolną z 1 na głównej przekątnej, (ii) uii =1 dla
i = 1; 2; :::; n, wówczas U będzie macierzą trójkątną górną z 1 na głównej przekątnej. Aby skonstruować rozkład,
wystarczy wykonać mnożenie macierzy L ó U i wykorzystać fakt, że L i U są macierzami trójkątnymi postaci (122),
wzwiązku z tym wiele składników w sumach stanowiących elementy obliczonej w ten sposób macierzy A zeruje się.
Macierze L i U buduje się w ten sposób, że w każdym kroku oblicza się nowe wartości wiersza macierzy U i kolumny
macierzy L. Algorytmy rozkładu LU możemy opisać nastepująco:
(i) gdy L zawiera 1 na głównej przekątnej - tzw. rozkład Doolittle a
dla k = 1; 2; :::; n
kĄ1
X
ukj = akj Ą lksusj; j = k; k +1; :::; n
s=1
PkĄ1
aik Ą lisusk
s=1
lik = ; i = k +1; k +2; :::; n (123)
ukk
(ii) gdy U zawiera 1 na głównej przekątnej - tzw. rozkład Crouta
dla k = 1; 2; :::; n
kĄ1
X
lik = aik Ą lisusk; i = k; k +1; :::; n
s=1
PkĄ1
akj Ą lksusj
s=1
ukj = ; j = k +1; k +2; :::; n (124)
lkk
LEMAT 68 (Tw. Sylvestera). Jeśli wszystkie n minory główne macierzy A są dodatnie, to macierz A posiada
rozkład LU.
DOWÓD pomijamy.
Kolejny, użyteczny w pewnych sytuacjach, rozkład macierzy A oparty jest na następującym:
TWIERDZENIE 69. Jeśli A jest rzeczywistą, symetryczną i dodatnio określoną macierzą, to ma jednoznaczny
rozkład A = LLT , gdzie L jest macierzątrójkątną dolnąz wartościami dodatnimi na głównej przekątnej.
DOWÓD pomijamy, ale przypomnijmy, że A jest symetryczna i dodatnio określona, gdy A = AT i xT Ax > 0 dla
dowolnego wektora x =0.
6
Algorytm zw. rozkładem Cholesky ego (albo Banachiewicza) jest szczególnym przypadkiem rozkładu LU (wówczas
U = LT ) i jest postaci:
dla k = 1; 2; :::; n
v
u
kĄ1
u
takk Ą X lks
2
lkk =
s=1
PkĄ1
aik Ą lislks
s=1
lik = ; i = k +1; k +2; :::; n (125)
lkk
2
Pl 2 2
UWAGA: Lem.68 gwarantuje, że lkk > 0, wzór na lkk daje również, że dla j � k zachodzi akk = lks � lkj, z
s=1
p
którego wynika, że jlkjj � akk dla wszystkich j. Oznacza to, że dowolny element macierzy L jest ograniczony przez
pierwiastek kwadratowy odpowiadajacego mu elementu diagonalnego macierzy A i implikuje, że elementy macierzy L
nie będą duże względem elementów macierzy A:
UWAGA: Podobną metodą rozwiązywania URL do rozkładu LU stanowi tzw. metoda LDLT , w której wyjściowy
układ rozkłada się na iloczyn trzech macierzy, gdzie L jest trójkątna dolna, D - diagonalna i LT - oczywiście trójkątna
górna. Szczegóły patrz np. [3].
7.3. Metody eliminacji
W tej części zajmiemy się tradycyjną formą metod eliminacji oraz ich mody
�kacjami, przede wszystkim metodą
eliminacji Gaussa. Główna idea tej metody sprowadza się do wykonywania operacji elementarnych na macierzy układu
A tak, by w końcowym efekcie uzyskać macierz trójkątną górną. W k-tym kroku eliminacji dokonujemy  wyzerowa-
nia wszystkich elementów macierzy leżących pod główna przekątną w k-tej kolumnie za pomocą operacji dodania
do danego wiersza k-tego wiersza pomnożonego przez odpowiednią liczbę różną od zera. Zatem po k Ą 1-ym kroku
macierz układu wygląda następująco:
2 3
a(k) ::: a(k) a(k) ::: a(k)
11 1kĄ1 1k 1n
6 7
::: ::: ::: ::: ::: :::
6 7
6
0 ::: a(k) a(k) ::: a(k) 7
6 7
kĄ1kĄ1 kĄ1k kĄ1n
6 7
6 7
0 ::: 0 a(k) ::: a(k)
kk kn
6 7
6
0 ::: 0 a(k) ::: a(k) 7
6 k+1k k+1n 7
4 5
::: ::: ::: ::: ::: :::
0 ::: 0 a(k) ::: a(k)
nn
nk
Algorytm, który opisuje proces eliminacji Gaussa jest nastepujący:
dla k = 1; 2; :::; n Ą 1
i = k +1; k +2; :::; n
aik
z = ; aij = aij Ą zakj; j = k +1; k +2; :::; n +1 (126)
akk
przy założeniu, że wszystkie elementy głównej przekątnej są różne od zera oraz oznaczeniu bi = ain+1. Uzyskany w
tensposóbukład można rozwiązać metodą (121).
Inną wersję eliminacji stanowi metoda Gaussa-Jordana, w której eliminujemy nie tylko elementy leżące poniżej
głównej przekątnej, ale równieżpowyżej, tak że w końcowym efekcie otrzymujemy macierz diagonalną. Algorytm jest
analogiczny, jak w przypadku eliminacji Gaussa z jedną różnicą
dla k = 1; 2; :::; nĄ 1
i = 1; :::; k Ą 1; k +1; k +2; :::; n
aik
z = ; aij = aij Ą zakj; j = k +1; :::; n +1 (127)
akk
Niestety metody eliminacji są niezadowalające dla pewnych, wydawałoby się, prostych układów równań.
PRZYKA
AD: Rozważmy układ:
� �� � � �
" 1 x1 1
;
1 1 x2 = 2
gdzie " jest pewnąmałą liczbąróżną od zera. Po zastosowaniu algorytmu Gaussa otrzymamy poniższy układ z macierza
trójkątną górną:
� �� � � �
" 1 x1 1
;
1 1
0 1 Ą x2 = 2 Ą
" "
którego rozwiązaniem jest para liczb:
1
2 Ą 2" Ą 1 1 Ą x2
"
x2 = = ź 1 oraz x1 = ź 0:
1
1 Ą " Ą 1 "
"
(w obliczeniach numerycznych, jeśli " jest dostatecznie małą liczbą, to wartości 2 Ą "Ą1 oraz 1 Ą "Ą1 będą obliczone
3
jako "Ą1). Oczywiście rozwiązaniem dokładnym jest para
1 1 Ą 2"
x1 = ź 1 oraz x2 = ź 1;
1 Ą " 1 Ą "
zatem znalezione rozwiązanie jest dokładne dla x2, ale ekstremalnie niedokładne dla x1. Wszystko oczywiście zależy
od wielkości " i używanej precyzji. Przy czym, tego typu niestabilność metody nie jest wcale spowodowana małą
wielkością współczynnika a11, ale jego względną małą wielkością w stosunku do pozostałych elementów w wierszu.
Z drugiej strony problem znika, gdy w rozważanym układzie zamienimy miejscami równania. Wowczas eliminacja
Gaussa prowadzi do układu:
� �� � � �
1 1 x1 2
;
0 1 Ą " x2 = 1 Ą 2"
którego rozwiązanie numeryczne jest już zadowalające (dokładne). Ostatecznie podsumowując, oczywiście taka niesta-
bilność metody eliminacji Gaussa jest spowodowana błędami zaokrąglenia operacji zmiennopozycyjnych.
Teraz krótko opiszemy algorytm zw. eliminacjąGaussaz częściowym wyborem elementu głównego (lub teżmetodą
Gaussa-Crouta). Jak wynikało z ostatniego przykładu, jednym ze sposobów ominięcia problemu niestabilności elimi-
nacji Gaussa jest zamiana miejscami wierszy, co w tym wypadku spowodowało przeniesienie  kłopotliwego elementu
poniżej głównej przekątnej, gdzie został wskutek procesu eliminacji  wyzerowany . Mody
�kacja klasycznej procedury
polega na znalezieniu najlepszego wiersza, za pomocą którego przeprowadzi się eliminację. W k-tym kroku, przed
dokonaniem eliminacji, szuka się wiersza z elementemnajw/
ekszym co do wartości bezwzględnej w k-tej kolumnie
(czyli tzw. elementu głównego) i zamienia się ten wiersz miejscami z k-tym wierszem. Można stosować rowni
z
metodę eliminacji z pełnym wyborem elementu głównego. Wowczas szuka się elementu głównego w całej pozostałej
do przekształcenia podmacierzy i dokonuje się odpowiedniej zamiany i wierszy, i kolumn.
Istniejąukłady równań liniowych, które dająsięrozwiązaćmetodą eliminacji Gaussa bez wyboru elementu głównego.
Jedna z takich klas tworzą układy z macierzami diagonalnie dominującymi, tzn. takimi, dla których spełniony jest
warunek
n
X
jaiij > jaijj dla 1 � i � n: (128)
j=1;j=i
6
UWAGA: (i) Eliminacja Gaussa bez wyboru elementu głównego zachowuje diagonalnie dominujący charakter
macierzy.
(ii) W poprzedniej części stwierdziliśmy, że rozkład Cholesky ego jest szczególnym przypadkiem rozkładu LU oraz
p
spełniona jest nierówność jlijj � aii. Zatem wybór elementu głównego w tej metodzie nie jest potrzebny, gdyż z
zacytowanej nierówności wynika, że elementy macierzy L są w porównaniu z elementami macierzy A zdecydowanie
mniejsze.
(iii) Algorytm eliminacji Gaussa jest również szczególnym przypadkiem rozkładu LU. Mianowicie macierzą U jest
macierz uzyskana po eliminacji, natomiast macierz L składa się z mnożników użytych w trakcie eliminacji. Również
wersja z wyborem elementu głównego jest metoda rozkładu: jeśli przez P oznaczymy macierz permutacji powstałą z
tablicy permutacji opisujących przestawianie wierszy, to PA = LU.
W pewnych zastosowaniach (interpolacja splajnami sześciennymi, równania różniczkowe, itd.) pojawia się często
uklad równań, w którym macierz układu A jest trójdiagonalna, tzn. poniższej postaci:
d1 c1 0 ::: 0 0
a1 d2 c2 ::: 0 0
0 a2 d3 ::: 0 0
(129)
::: ::: ::: ::: ::: :::
0 0 0 ::: dnĄ1 cnĄ1
0 0 0 ::: anĄ1 dn
W przypadku zastosowania eliminacji Gaussa, uwzględniając, że wiele elementów macierzy (129) jest równa 0, można
uzyskać się stosunkowo proste przekształcenia postaci:
dla i = 2; 3; :::; n
aiĄ1
miĄ1 =
diĄ1
di = di Ą miĄ1ciĄ1; bi = bi Ą miĄ1biĄ1 (130)
4
aby sprowadzić macierz układu stała się trójkątnie górna (w rzeczywistości dwudiagonalna). Postępowanie wstecz jest
również prostsze:
bn
xn =
dn
bi Ą cixi+1
xi = ; i = n Ą 1; nĄ 2; :::; 1 (131)
di
Alternatywną wersję powyższej metody stanowi transformacja macierzy układu do postaci trójkątnej dolnej (oczy-
wiście dwudiagonalnej) za pomocą następujących przekształceń:
dla i = n Ą 1; nĄ 2; :::; 1
ci
mi =
di+1
di = di Ą mi1ai; bi = bi Ą mibi+1 (132)
i zastosowanie wówczas postępowania wprzód:
b1
x1 =
d1
bi Ą aiĄ1xiĄ1
xi = ; i =2; 3; :::; n (133)
di
Na koniec zajmijmy się kosztem rozważanych metod. Koszt obu wersji rozkładu LU i metody eliminacji Gaussa
1 1 1 1 5
jest taki sam. Trzeba wykonaćdokładnie n3 + n2 Ą n mnożeńoraz n3 + n2 Ą n dodawań, w eliminacji Gaussa-
3 3 3 2 6
1 1 1 1 1
Jordana: n3 + n2 Ą n mnożeń oraz n3 Ą n dodawań,w metodzie Cholesky ego: n operacji pierwiastkowania,
2 2 2 2 2
1 3 1 1 7
n3 + n2 + n mnożeńoraz n3 +n2 Ą n dodawań. W przypadku macierzy trojdiagonalnej wystarczy tylko 5nĄ4
6 2 3 6 6
mnożeńoraz 3n Ą 3 dodawań.
7.4. Analiza błędu metod dokładnych
Aby przedyskutować błędy w problemach numerycznych zawierających wektory i macierze, wykorzystuje się po-
jęcie normy. W analizie numerycznej najczęściej używa się następujących norm wektorów:
� !1=2
n
X
kxk2 = x2
i
i=1
kxk1 = max jxij
i=1;:::;n
n
X
kxk1 = jxij
i=1
nazywanych odpowiednio normą euklidesową, normą `1 (maksymalną) i normą `1 oraz indukowanych norm macierzy:
n
X
kAk1 = max jaijj
i=1;:::;n
j=1
n
X
kAk1 = max jaijj
j=1;:::;n
i=1
UWAGA: Zależności pomiędzy powyższymi normami wektorów przedstawia nierówność:
kxk1 �kxk2 �kxk1
Niech URL ma standardowąpostać Ax = b i przypuścmy, że macierz A jest odwracalna (tj. det A =0). Rozważmy
6
następuującą kwest/
e: jeśli macierz AĄ1 jest przybliżona przez pewną macierz B, to rozwiązanie x = AĄ1b jest
przybliżone przez wektor x = Bb. Obliczmy błąd bezwzględny takiego przybliżenia:
~
kx Ą xk = kx Ą Bbk = kx Ą BAxk = k(I Ą BA)xk �kI Ą BAkkxk ;
~
a jeśli interesuje nas błąd względny, to
kx Ą xk
~
�kI Ą BAk :
kxk
Powyższa nierówność daje oczywiście oszacowanie błędu względnego z góry.
5
~
Teraz zbadajmy błąd w sytuacji, gdy wekter b jest przybliżony pewnym wektorem b oraz x i x spełniają równania
~
~
Ax = b i A~ = b :
x
�
�AĄ1b b� � Ą ~ � � � � � Ą ~�
�AĄ1(b b)� �AĄ1� �b b�
kx Ą xk = Ą AĄ1~� = � � �
~ � � �
i wielkość względna
� � � �
�b Ą ~� �b Ą ~�
� �
b b
� � � �
� �
� �AĄ1� kAxk �AĄ1�
�
kx Ą xk ��AĄ1� �b Ą ~� = � � � kAkkxk ;
~ � �
b
kbk kbk
skąd
� �
�b Ą ~�
b
� �
�
kx Ą xk
~
�AĄ1� :
�
�kAk (134)
kxk kbk
DEFINICJA 70. Liczba
�
�AĄ1�
�
�(A) =kAk (135)
jest nazywana wskaz
nikiem uwarunkowania macierzy A (jego wielkość zależy od wybranej normy).
Nierówność (134) oznacza, że błąd względny nie przekracza więcej niż �(A) razy błędu względnego b. Wid�c
również, że jeśli wskaz
nik uwarunkowania jest mały, to niewielkie zakłócenia wektora b prowadzą domałych zakłóceń
rozwiązania x. Przy czym, zawsze �(A) � 1.
PRZYKA
AD: Zbadajmy wskaz
nik uwarunkowania macierzy
� � � �
1
1 1 + " 1 Ą1 Ą "
A = , wówczas AĄ1 = :
1 Ą " 1 Ą1+" 1
"2
�
�AĄ1� = 2+" Zatem z de
�nicji �(A) = Ą2+"ó2. J�sli
�
Jeśli użyjemy normy maksymalnej, to kAk1 =2 +" oraz
"2 "
1
" � 0:1, to �(A) � 441, zaś jeśli " � 0:01, to �(A) � 40401. W takim przypadku, małe względne zaburzenie wektora
b może spowodować ponad 40000 razy większe zaburzenie rozwiązania układu Ax = b.
Jeśli rozwiązujemy numerycznie URL, oczywiście otrzymujemy rozwiązanie przybliżone x zamiast dokładnego x.
~
Różnica między nimi jest wektorem błędu
e = x Ą x
~
Można zbadać jakość przybliżenia, znajdując tzw. wektor residuum
r = b Ą A~
x:
Zależność Ae = r między wektorem błędu i wektorem residuum jest bardzo ważna.
~, ~
Zauważmy, że x jest rozwiązaniem dokładnym układu A~ = b który ma zakłóconą prawą stronę b = b Ą r.
~ x
Spróbujemy zbadać jaki jest związek między błędami względnymi x i b. Poniższe twierdzenie pokazuje, że wskaz
nik
uwarunkowania odgrywa tu ważną rolę:
TWIERDZENIE 71.
1 krk kek krk
� � �(A)
�(A) kbk kxk kbk
DOWÓD: Nierówność po prawej stronie może być zapisana jako
�
�AĄ1�
�
kekkbk �kAk krkkxk
i jest prawdziwa, gdyż
�
�AĄ1r� kAxk � � �
� �AĄ1� krkkAkkxk
kekkbk =
(w rzeczywistości nierówność po prawej stronie to (134)). Nierówność polewej stronie może być zapisana jako
�
�AĄ1�
�
krkkxk �kAk kbkkek
i wynika z następującej
�
�AĄ1b� �kAkkek � �
� �AĄ1� kbk :
krkkxk = kAek
Ą
Macierz z dużym wskaz
nikiem uwarunkowania jest nazywana z
le uwarunkowaną. Dla z
le uwarunkowanej macierzy
A rozwiązanie układu Ax = b jest bardzo wrażliwe na małe zmiany wektora b. Innymi słowy, aby osiągnąć duża pre-
6
cyzję przy obliczeniu wketora x, powinniśmy wymagać znacznie wyższej precyzji wektora b. J�sli wskaz
nik uwarunk-
owania jest niezbyt duży, to macierz nazywamy dobrze uwarunkowaną.
TWIERDZENIE 72. Jeśli A jest macierzą rozmiaru n Ł n taką, że kAk < 1 (gdzie kók jest dowolną normą
macierzową), to macierz I Ą A jest odwracalna oraz
1
X
(I Ą A)Ą1 = Ak: (136)
k=0
DOWÓD pomijamy, ale zauważmy, że szeregi postaci jak po prawej stronie (136) są nazywane szeregami Neu-
manna. PowŹzsze twierdzenie ma istotne teoretyczne i praktyczne konsekwencje związane z uwarunkowaniem macierzy.
Zauważmy, że z równania (136) otrzymujemy
1 1
� �
�(I Ą A)Ą1� � X � � � X kAkk = 1 :
�Ak�
� �
1 ĄkAk
k=0 k=0
7.5. Iteracyjne poprawianie rozwiązania
Jeśli x(0) jest przybliżonym rozwiązaniem równania Ax = b, torozwiązanie dokładne jest dane przez
ł �
x = x(0) + AĄ1 b Ą Ax(0) = x(0) + e(0);
Ą ó
gdzoe e(0) = AĄ1 b Ą Ax(0) jest wektorem błędu. Wektorem residuum odpowiadającym przybliżeniu x(0) jest
r(0) = b Ą Ax(0). Oba wektory można obliczyć bez kłopotu, jeśli zauważymy, że wektor e(0) może być otrzymany
jako rozwiązanie równania
Ae(0) = r(0);
(wtedy nie musimy szukać AĄ1). Powyższe rozważania prowadzą do procedury nuemrycznej zw. iteracyjnym popraw-
ianiem rozwiązania, którą szczegółowo opiszemy w postaci algorytmu:
1. Rozw/
azać dowolną metodą dokładną wyjściowy URL Ax = b. Oznaczmy przez x(0) uzyskane w ten sposób
rozwiązanie.
2. Obliczyć r(0) = b Ą Ax(0).
3. Rozwiązać URL Ae(0) = r(0), aby otrzymać wektor e(0).
4. Obliczyć x(1) = x(0) + e(0).
5. Aby otrzymać lepsze przybliżenia x(2), x(3), ... powtórzyć kroki 2-4.
UWAGA: (i) Zauważmy, że w kroku 2 rozwiązujemy układ równań z tą samą macierzą współczynników, co układ
wyjściowy, tylko z innym wektorem wyrazów wolnych. Zatem użycie dowolnej metody rozkładu umożliwia szybkie
znalezienie e(0) bez konieczniości stosowania pełnej procedury rozwiązywania URL.
(ii) Idealnie byłoby otrzymać r(k) jako wektor zerowy, co jest możliwe, ale mało prawdopodobne, zatem jako kryterium
końca używa się zwykle nierówności
� �
(k)
�e �
� �
� 0:1:
(k+1)
�x �
Przeanalizujmy powyższy algorytm teoretycznie. Proces iteracyjny może być zapisany jako
ł �
x(k+1) = x(k) + B b Ą Ax(k) dla k � 0; (137)
gdzie B oznacza przybliżenie AĄ1.
TWIERDZENIE 73. Jeśli kI Ą BAk < 1, to metoda iteracynego porawiania rozwiązania opisana przez (137)
generuje ciąg wektorów zbieżny do rozwiązania układu Ax = b.
DOWÓD: Wykorzystując (137), otrzymamy
ł �
x(k+1) Ą x = x(k) Ą x + B b Ą Ax(k) =(I Ą BA)(x(k) Ą x);
skąd wynika, że
� �
�x(k+1) Ą x� �kI Ą BAk � Ą x� �kI Ą BAk2 � Ą x� �kI Ą BAkk � Ą x� :
�x(k) � �x(kĄ1) � �x(0) �
� � � � � � � �
� �
�x(k+1) Ą x� zbiega do 0 przy k !1 dla dowolnego x(0).Ą
Ponieważ kI Ą BAk < 1, tobłąd
7
7.6. Wstępne uwarunkowywanie macierzy URL
Przy rozwiązywaniu układów równań liniowych w szczególnie kłopotliwych sytuacjach stosuje się pewne techniki,
które mogą zminimalizować problemy wynikające ze złego uwarunkowania macierzy:
1. wstępne uwarunkowywanie przez równoważenie wierszy:
równoważenie wierszy jest procesem dzielenia każdego wiersza macierzy współczynników przez maksymalny, co do
wartości bezwzględnej, element w tym wierszu, tj. pomnożenie i-tego wiersza przez
1
ri = dla i =1; :::; n.
maxj=1;:::;n jaijj
Wowczas nowe elementy aij spełniają warunek maxj=1;:::;n j^ijj =1 dla i =1; :::; n. W praktyce mnożnik ri wybiera
^ a
się jako liczbę binarną, tj. postaci 2m najbliższą odwrotności powyższego maksimum. Dokonuje się tego, aby uniknąć
wprowadzania dodatkowych błędów zaokrąglenia. Oczywiście przez ri trzeba również pomnożyć wartości bi. Zatem
w rzeczywistości rozwiązuje się układ równań (RA)x =(Rb), gdzie R = diag(ri), czyli
n
X
(riaij)xj = ribi dla i =1; :::; n:
j=1
2. wstępne uwarunkowywanie przez równoważenie kolumn:
równoważenie kolumn jest procesem podobnym do poprzedniego, tzn. monoży się j-tą kolumnę przez
1
cj = dla j =1; :::; n.
maxi=1;:::;n jaijj
Tutaj nowy układ równań jest postaci
� ś
n
X
xj
(cjaij) = bi dla i =1; :::; n
cj
j=1
lub inaczej (AC)(CĄ1x) =b, gdzie C = diag(cj).
3. pełny wybór elementu głównego - opisany w rozdz. 7.3. Metody eliminacji.
4. wstępne uwarunkowywanie podczas procesu eliminacji - jeden z dwóch sposobów wstępnego uwarunkowywania
stosowany w trakcie realizacji procesu eliminacji, tj. przed wykonaniem każdego kroku.
5. iteracyjne poprawianie rozwiązania - technika dyskutowana w poprzedniem rozdziale.
PRZYKA
AD: Porównajmy wstępne uwarunkowywanie macierzy przez równoważenie wierszy i kolumn dla macierzy
� �
1 R
A = , gdzie R =0:
6
2 0
Równoważenie wierszy
� � � �
1 R RĄ1 1
!
2 0 1 0
a równoważenie kolumn
� � � �
1
1 R 1
2
! :
2 0 1 0
Wskaz
niki uwarunkowania powyższych macierzy wynoszą (dla kók1 i kók1):
Ą ó2
1 9
�(A) = (R +1) , �(AZW ) = 1+RĄ1 , �(AZK) = ;
2 4
zatem np. dla R =108 mamy
9
�(A) =5 ó 107, �(AZW ) ź 1, �(AZK) = :
4
7.7. Metody iteracyjne rozwiązywania układów równań liniowych
Metody rozkładu i eliminacji rozwiązywania URL określa sięjako bezpośrednie lub dokładne. Umożliwiają znalezie-
nie w skończonej i znanej z góry liczbie kroków znalezienie rozwiązania dokładnego, jeśli nie pojawiałby się błąd
zaokrąglenia związany z operacjami zmiennoprzecinkowymi. W przeciwieństwie do nich, metody przybliżone generują
ciąg, który zbiega do rozwiązania w nieznanej z góry liczbie kroków i dlatego nazywane są metodami iteracyjnymi.
Obliczenia zatrzymuje się poosiągnięciu pewnej założonej dokładności.
8
Metody iteracyjne mają przewagę nad dokładnymi przy rozwiązywaniu dużych układów liniowych, jeśli chodzi
o szybkość zbi
zności oraz wykorzystanie pamięci maszyny. Są one zwłaszcza bardzo efektywne dla układów z
macierzami, które zawierają dużo zer.
Rozważymy teraz metody iteracyjne w ogólnym sensie. Załóżmy, że macierz A jest nieosobliwa. Wówczas układ
Ax = b ma rozwiązanie dla dowolnego wektora b. Wyjściowy problem można zapisać wrównoważnej postaci
Qx =(Q Ą A)x + b (138)
dla pewnej macierzy Q nazywanej macierzą rozszczepienia. PowŹzszy zapis sugeruje proces iteracyjny, zde
�niowany
następująco:
Qx(k) =(Q Ą A)x(kĄ1) + b dla k � 1: (139)
Wektor początkowy x(0) może być dowolny. Naszym celem jest wybrać takie Q, by:
� �
(i) wyrazy ciągu x(k) były łatwe do obliczania,
� �
(ii) ciąg x(k) był szybko zbieżny do rozwiązania.
� �
Zauważmy, że jeśli ciąg x(k) jezt zbieżny, powiedzmy do wektora x, to x jest rozwiązaniem, bowiem przechodząc
do granicy w równaniu (139), otrzymamy równanie (138), które jest równoważne wyjściowemu URL.
Załóżmy teraz dodatkowo, że również macierz Q jest nieosobliwa, wówczas równanie (139) może być rozwiązane
względem wektora x(k):
Ą ó
x(k) = I Ą QĄ1A x(kĄ1) + QĄ1b: (140)
Taka postać jest wygodna do analizy teoretycznej, choć numerycznie otrzymuje się x(k) bez użycia QĄ1.
Rozwiązanie x spełnia równanie
Ą ó
x = I Ą QĄ1A x + QĄ1b;
Ą ó
czyli x jest punktem stałym odwzorowania x ! I Ą QĄ1A x + QĄ1b. Odejmując stronami powyższe równanie od
(140), otrzymamy
ł �
Ą ó
x(k) Ą x = I Ą QĄ1A x(kĄ1) Ą x :
Teraz, używając pewnej normy wektorowej i indukowanej przez nią macierzowej, będziemy mieli
� �
�
�x(k) Ą x� � � Ą QĄ1A� � Ą x� ;
�x(kĄ1) �
�I
� � � �
i powtarzając wielokrotnie szacowanie z góry, otrzymamy nierówność
� �
�k �x(0) �
�x(k) Ą x� � � Ą QĄ1A� � Ą x� :
�I
� � � �
� �
�I
Wówczas, jeśli Ą QĄ1A� < 1, to
� �
�x(k) �
lim Ą x� =0
�
k!1
� �
�I
dla dowolnego x(0). Zauważmy, że założenie Ą QĄ1A� < 1 implikuje na podstawie tw. 72 o szeregu Neumanna,
odwracalność macierzy QĄ1A i skąd mamy
�A, �
�I
TWIERDZENIE 74. Jeśli Ą QĄ1A� < 1 dla pewnej indukowanej normy macierzowej, to ciąg generowany
przez wzór (139) jezt do rozwiazania układu Ax = b dla dowolnego wektora początkowego x(0).
�zbieżny � � �
(k)
�I �x
Jesli norma ą = Ą QĄ1A� < 1, to bezpiecznie można zatrzymać proces iteracyjny, gdy norma Ą x(kĄ1)�
jest mała. Można udowodnić, że
� � �
�x(k) Ą x� � ą � Ą x(kĄ1)� :
�x(k)
� � � �
1 Ą ą
Wynika z tego, że przerywanie iteracji po uzyskaniu poprawki mniejszej od dopuszczalnego błędu jest poprawne tylko
wtedy, gdy ą � 1=2.
7.7.1. Ogólna analiza
Rozważmy proces iteracyjny zde
�niowany równaniem
x(k) = Gx(kĄ1) + c; (141)
gdzie G jest ustaloną macierzą wymiaru n Ł n, c jest ustalonym wektorem z Rn. Zauważmy, że taka iteracja jest
bardziej ogólna niż opisana wzorem (140), mianowicie G = I Ą QĄ1A oraz c = QĄ1b są szczególnym przypadkiem
powyższej. Spróbujmy znalez
ć warunki konieczny i dostateczny dotyczące G, które zapewnią zbieżność iteracji (141)
9
dla dowolnego wektora poczatkowego.
OKREŚLENIE: Wartości własne macierzy A są liczbami zespolonymi �, dla których macierz A Ą �I nie jest
odwracalna. Zatem są to pierwiastki równania charakterystycznego macierzy A postaci:
det (A Ą �I) =0:
Promień spektralny macierzy A de
�niuje się jako
� (A) =max fj�j : det (A Ą �I) =0g
czyli jako największącodowartości bezwzglednej wartośćwłasną macierzy lub inaczej jako promień najmniejszej kuli
o środku 0 w przestrzeni zespolonej zawierającej wszystkie wartości własne.
A
atwo pokazać, że wartościami własnymi macierzy trójkątnej są elementy na głównej przekątnej.
TWIERDZENIE: Promieńspektralnyjakofunkcja spełnia równanie
� (A) =inf kAk ;
kók
w którym kres dolny jest brany po wszystkich indukowanych normach macierzowych.
Powyższy fakt pozwoli nam udowodnić twerdzenie o warunku koniecznym i dostatecznym zbieżnośći metod itera-
cyjnych:
TWIERDZENIE 75. Warunkiem koniecznym i dostatecznym na to, by formuła iteracyjna (141) generowała ciąg
zbieżny do (I Ą G)Ą1 c dla dowolnego wektora początkowego x(0), jest nierówność
� (G) < 1.
WNIOSEK 76. Formuła iteracyjna (139), tj. Qx(k) =(Q Ą A)x(kĄ1) + b generuje ciąg zbieżny do rozwiązania
układu Ax = b dla dowolnego x(0), wtedy i tylko wtedy, gdy
Ą ó
� I Ą QĄ1A < 1:
7.7.2. Metoda Richardsona
Najprostszą z możliwych sytuacji powyższego typu jest metoda Richardsona, w której macierzą Q jest macierz
jednostkowa I, tj.
x(k) =(I Ą A) x(kĄ1) + b = x(kĄ1) + r(kĄ1); (142)
gdzie r(kĄ1) jest wektorem residuum, zde
�niowanym jako r(kĄ1) = b Ą Ax(kĄ1). Zgodnie z tw.75, iteracja Richard-
sona jest skuteczna, jeśli kI Ą Ak < 1 (np. dla macierzy diagonalnie dominujących). Wzór (142) można zapisać w
wygodniejszej postaci:
n
X
x(k) = x(kĄ1) + bi Ą aijx(kĄ1) dla i =1; :::; n.(143)
i i j
j=1
7.7.3. Metoda Jacobiego
W tej metodzie macierz Q jest macierzą diagonalną, której wyrazy niezerowe pokrywają się z wyrazami macierzy
A, tzn. aii. W tym przypadku dowolnym elementem macierzy QĄ1A jest aij=aii, a diagonalne elementy są 1, skąd
Ż Ż
n
X
� �
�I Ą QĄ1A�1 = max Żaij Ż :
Ż Ż
Ż Ż
i=1;:::;n
aii
j=1;j=i
6
Ztej równości i wzoru (128) wynika następujące:
TWIERDZENIE 77. Jeśli macierz A jest diagonalnie dominująca, to ciąg generowany przez iterację Jacobiego jest
zbieżny do rozwiązania układu Ax = b dla dowolnego wektora początkowego.
Iterację Jacobiego podaje następujący wzór:
Pn
bi Ą aijx(kĄ1)
j=1;j=i j
6
x(k) = dla i =1; :::; n: (144)
i
aii
UWAGA: Ten algorytm oraz inne, póz
niejsze mogą być zrealizowane mniejszym kosztem przez wykonanie wszys-
tkich dzieleń, zanim zacznie się proces iteracyjny, czyli
bi aij
bi = oraz aij = dla odpowiednich i oraz j;
aii aii
o ile nie potrzebujemy już dalej macierzy A ani wektora b. Zatem w rzeczywistości sprowadza się to do rozwiązania
10
przeskalowanego układu DĄ1Ax = DĄ1b, gdzie D = diag(aii).
7.7.4. Metoda Gaussa-Seidela
W metodzie Jacobiego najpierw są obliczane wszystkie nowe współrzędne wektora x, a potem dokonywana jest za-
miana. W iteracji Gaussa-Seidela, którąmożna traktowaćjakomody
�kację Jacobiego, przy obliczaniu i-tej współrzęd-
nej wektora x, wykorzystuje sięjuż obliczone wcześniej i Ą 1 współrzędnych, co daje większą efektywność tej metody,
ale uniemożliwa dokonywanie obliczeń równoległych.
Macierz Q w metodzie Gaussa-Seidela jest dolnątrójkątna częściąmacierzy A, włącznie z głównąprzekątną, zatem
algorytm można zapisać następująco:
PiĄ1 Pn
bi Ą aijx(k) Ą aijx(kĄ1)
j=1 j ;j=i+1 j
x(k) = (145)
i
aii
TWIERDZENIE 78. Jeśli macierz A jest diagonalnie dominująca, to metoda Gaussa-Seidela jest zbieżna dla
dowolnego wektora początkowego.
7.7.5. Metoda nadrelaksacji
Następnym ważnym przykładem metod iteracyjnych jest metoda nedrelaksacji zw. krótko SOR. Najpierw zapoz-
namy się z twierdzeniem zawierającym warunki zbieżności tej metody:
TWIERDZENIE 79. Przypuśćmy, że macierz rozszczepienia Q jest wybrana jako !Ą1D Ą C z parametrem rzeczy-
wistym !, gdzie D jest dowolną dodatnio określoną macierzą Hermite a i C jest dowolną macierzą spełniającą rów-
nanie C + Cń = D Ą A. Jeśli A jest dodatnio określoną macierzą Hermite a, Q jest nieosobliwa, 0 < ! < 2, to
metoda nadrelaksacji jest zbieżna dla dowolnego wektora poczatkowego.
Pozostaje teraz jeszcze przypomnieć, że macierz A jest nazywana macierzą Hermite a, jeśli A = Ań oraz że Ań oz-
nacza sprzężoną transpozycję macierzy A, czyli Ań =(aji) w ogólnym przypadku macierzy o elementach zespolonych.
ą
Powyższe twierdzenie nie precyzuje, jak wybierać macierze D i C. Zwykle przyjmuje s/
e, że D jest macierzą
diagonalną z elementami macierzy A, zaś C dolną trójkątną częścią macierzy A bez wyrazów głównej przekątnej.
Jeśli zauważymy, że wzór (145) opisujący metodę Gaussa-Seidela można zapisać w postaci:
(kĄ1)
x(k) = x(kĄ1) + ri ; gdzie
i i
PiĄ1 Pn
bi Ą aijx(k) Ą aijx(kĄ1)
(kĄ1) j=1 j ;j=i j
ri = ;
aii
to łatwo zobaczyć, że metoda SOR jest pewnym jej uogólnieniem, które pozwala często przyspieszyć zbieżność:
(kĄ1)
x(k) = x(kĄ1) + !ri ; (146)
i i
(kĄ1)
gdzie ri jest określone tak jak powyżej.
7.7.6. Macierze iteracyjne
Przypuścmy, że macierz A jest podzielona w następujący sposób
A = D Ą CL Ą CU ;
gdzie D =diag(A), CL jest ujemną ściśle dolnątrójkątnączęścią A i CU jest ujemną ściśle górnątrójkątnączęścią A.
Dla wszystkich przedstawionych metod iteracyjnych, opiszemy kluczowe macierze i ogólny proces iteracyjny:
1. Richardson
�
Q = I
G = I Ą A
x(k) =(I Ą A)x(kĄ1) + b
2. Jacobi
�
Q = D
G = DĄ1(CL + CU)
Dx(k) =(CL + CU)x(kĄ1) + b
11
3. Gauss-Seidel
�
Q = D Ą CL
G =(D Ą CL)Ą1CU
(D Ą CL)x(k) = CUx(kĄ1) + b
4. SOR
�
Q = !Ą1(D Ą !CL)
G =(D Ą !CL)Ą1 (!CU +(1Ą !)D)
ł �
(D Ą !CL)x(k) = ! CU x(kĄ1) + b +(1Ą !)Dx(kĄ1)
7.8. Zastosowanie rozwiązywania URL
Metody rozwiązywania układów równan liniowych mogą byc stosowane również do innych celów zwiazanych z
algebrą macierzy, a mianowicie: obliczania wyznacznika macierzy oraz odwracania macierzy.
7.8.1. Obliczanie wyznacznika macierzy
Operacje elementrane wykonywane na macierzy (wymienione na początku rozdziału 7) nie zmieniająwartości wyz-
nacznika, z jednym wyjątkiem: zamiana miejscami dwóch wierszy lub kolumn powoduje zmianę znaku wyznacznika.
Zatem do obliczania wyznacznika można zastosować:
1. dowolną metodę eliminacji rozwiązywania URL, która pozwala sprowadzić macierz do postaci diagonalnej lub
trójkątnej, a następnie wzór
n
Y
det(A) =(Ą1)p aii;
i=1
gdzie p oznacza liczbę zamian wierszy lub kolumn, jeśli były one wykonywane podczas eliminacji,
2. dowolną metodę rozkładu, która pozwala rozłożyć macierz A na iloczyn dwóch macierzy trójkątnych LU (gdyż
zgodnie z tw. Cauchy ego wyznacznik iloczynu macierzy jest równy iloczynowi wyznaczników poszczególnych czyn-
ników) i wzór
n n n
Y Y Y
2
det(A) = lii, det(A) = uii lub det(A) = lii
i=1 i=1 i=1
w zależności odmetodyrozkładu (Crouta, Doolitle a lub Cholesky ego).
7.8.2. Odwracanie macierzy
Niech A bedzie macierzą nieosobliwą i niech y(1); y(2); :::; y(n) oznaczają kolumny macierzy odwrotnej AĄ1.
Wówczas
h i h i
I = AAĄ1 = A y(1); y(2); :::; y(n) = Ay(1); Ay(2); :::; Ay(n) ;
czyli Ay(j); j = 1; :::; n stanowi kolumnę macierzy jednostkowej, czyli wektor e(j), którego j-ta współrzędna jest
równa 1, a pozostałe- 0. W ten sposób otrzymujemy zbiór n-układów równań liniowych z ta sama macierząwspółczyn-
ników i różnymi prawymi stronami:
Ay(1) = e(1); :::; Ay(n) = e(n):
Innymi słowy; kolumny macierzy odwrotnej AĄ1 są rozwiązaniami układów równań liniowych z macierzą współczyn-
ników A i z prawymi stronami równymi kolumnom macierzy jednostkowej I. Zatem szczególnie poleca się tumetody,
które pozwalają prowazdić obliczenia równoległe, tzn. metody eliminacji, zwłaszcza jeżeli zamiast jednego wektora b
użyjemy macierzy jednostkowej.
UWAGA: Można udowodnić, że macierze odwrotne do trójkątnych również są trójkątne tego samego typu, zatem
znając rozkład macierzy A, łatwo znalez
ć macierz odwrotną AĄ1, gdyż
AĄ1 =(LU)Ą1 = UĄ1LĄ1:
Odwrotności macierzy trójkątnych oblicza się stosunkowo łatwo: współrzędne wektora y(j); j = 1; :::; n można
obliczyć za pomocą podstawiania wprzód lub wstecz (w zależności odtypu trójkątności ).
12
7.9. Układy liniowe nadokreślone
Pojęcie układu liniowego nadokreślonego pojawia się, gdy zamierzamy dopasować do pewnych danych liniowy
model matematyczny i wykonuje się więcej pomiarów niż jest niewiadomych. Wówczas w układzie równań liniowych
Ax = b macierz A jest prostokątna o rozmiarach m Ł n, gdzie m � n i wektor b ma m składowych. Należy
znalez
ć wektor x o n składowych, nazywany rozwiązaniem średniokwadratowym, taki, że Ax jest najlepszym przy-
bliżeniem wektora b. Na ogół nie można go spełnić dokładnie i istnieje wiele możiwych sposobów określenia najlep-
szego rozwiązania. Jednym ze sposobów jest metoda najmniejszych kwadratów, prowadzacą do stosunkowo prostych
obliczeń, w której szuka się wektora x minimalizującego normę euklidesową residuum, tzn. wielkość
p
krk2 = rT r;
gdzie r = b Ą Ax.
W pewnych zastosowaniach bardziej naturalne mogą bżc inne normy, ale może to spowodować, że znalezienie
wektora x jest znacznie trudniejsze. Często rozwiązanie średniokwadratowe dobrze przybliża również rozwiązania w
innych normach.
7.10.1. Równania normalne
Istotną własność rozwiązania średniokwadratowego podaje następujące:
TWIERDZENIE 80. Niech A będzie macierzą rzeczywistą o rozmiarach m Ł n, a b - wektorem o m składowych.
Jeśli x spełnia równanie
AT (b Ą Ax) =0; (147)
to jest rozwiązaniem średniokwadratowym.
Z równania (147) wynika, że dla każdego wektora z zachodzi równość
(Az)T (b Ą Ax) =0;
która oznacza,że residuum b Ą Ax dla rozwiązania średniokwadratowego x jest ortogonalne do wszystkich wektorów
z przestrzeni R(A) rozpiętej na kolumnach macierzy A. Geometrycznie interpretując, rozwiązanie średniokwadratowe
rozkłada wektor b na dwie składowe: b = Ax + r, z których Ax jest rzutem ortogonalnym wektora b na R(A), z�s r
jest ortogonalne do R(A).
Z równania (147) wynika również, że rozwiązanie średniokwadratowe spełnia tzw. równania normalne
(AT A)x = AT b: (148)
Macierz AT A powyższego układu jest symetryczna, dodatnio określona i rozmiaru n Ł n. Zatem układ (148) jest
układem n równań liniowych z n niewiadomymi składowymi wektora x. Jesli macierz AT A jest nieosobliwa, to
można go rozwiązać za pomocą najmniej kosztownych metod dokładnych, np. metodą Cholesky ego lub symteryczną
eliminacją Gaussa (i tak najbardziej kosztowne jest utworzenie równań normalnych).
Warunek jednoznaczności rozwiązania podaje poniższe:
TWIERDZENIE 81. Macierz AT A jest nieosobliwa wtedy i tylko wtedy, gdy kolumny macierzy A są liniowo
niezależne.
W przeciwnym razie, jeśli kolumny macierzy A są liniowo zależne, to istnieje wiele rozwiązań średniokwadra-
towych x. Ze względu na jednoznaczność rzutu ortogonalnego wszystkie te rozwiązania daja to samo residuum r.
PRZYKA
AD 1: Rozważmy nastepujące zagadnienie: oszacować wysokość trzech punktów A, B i C nad poziomem
morza, mając dane różnice wysokości pomiędzy nimi oraz pomiędzy nimi a punktami D, E i F leżącymi na poziomie
morza. Każdy pomiar dał związek liniowy dla wysokości xA, xB i xC, comożna wyrazić w postaci macierzowej
2 3 2 3
1 0 0 1
6 7 2 3 6 7
0 1 0
6 7
xA 6 2 7
6 7 6 7
0 0 1 3
6 7 4 5 6 7
xB = ;
6 7 6 7
Ą1 1 0 1
6 7
xC 6 2 7
4 5 4 5
0 Ą1 1
Ą1 0 1 1
skąd mamy układ równań normalnych
2 3 2 3 2 3
3 Ą1 Ą1 xA Ą1
4 5 4 4 5
Ą1 3 Ą1 xB 5 = 1 :
Ą1 Ą1 3 xC 6
13
Za pomocą symetrycznej eliminacji Gaussa otrzymamy układ trójkątny
2 3 2 3 2 3
3 Ą1 Ą1 xA Ą1
8 2
4 4 5
0 Ą4 5 4 xB 5 = ;
3 3 3
0 0 2 xC 6
Ł ńT
5 7 1 3
którego rozwiązaniem są xA = , xB = i xC =3. Odpowiednie residuum Ą1; Ą1; 0; ; ; Ą3 jest ortogonalne
4 4 4 4 2 4 4
do wszystkich kolumn macierzy A.
Wpraktyce często spotyka się przypadek, gdy kolumny macierzy A są prawie liniowo zależne i wówczas równani
normalne nie dają najlepszego efektu.
PRZYKA
AD 2: Rozważmy układ Ax = b dla
2 3 2 3
1 1 1 1
6 7 6 7
" 0 0 0
6 7 6 7
A = oraz b = ;
4 5 4 5
0 " 0 0
0 0 " 0
gdzie j"j ż1. Może to odpowiadać sytuacji, gdy sumę x1 + x2 + x3 można zmierzyć z w/
ekszą dokadnóscią niż
poszczególne składniki. Otrzymamy
2 3 2 3
1+"2 1 1 1
4 5 4 5
AT A = 1 1 + "2 1 oraz AT b = 1 ;
1 1 1 + "2 1
h iT
1 1 1
skąd x = ; ; . J�sli " =10Ą4 i używamy precyzji 7- lub 8-cyfrowej, to suma 1+"2 =1:00000001
1+"2 1+"2 1+"2
zostanie zaokrąglona do 1 i macierz AT A jest osobliwa.
UZUPEA
NIENIE: Symetryczna eliminacja Gaussa jest szczególnym przypadkiem klasycznej eliminacji Gaussa,
gdy macierz ukłądu jest symetryczna. Jeśli eliminację wykonuje się bez przestawiania wierszy i kolumn, to a(k) = a(k)
ij ji
dla k � i; j � n, tzn. przekształcane elementy tworzą macierze symetryczne stopnia n +1Ą k dla k =2; 3; :::; n.
Wystarczy zatem obliczać tylko elementy leżące w A(k) na głównej przekątnej i nad nią. Algorytm jest realizowany
według wzorów:
dla k = 1; 2; :::; n Ą 1
i = k +1; k +2; :::; n
aki
z = ;
akk
aij = aij Ą zakj; j = i; i +1; :::; n
Liczba operacji w tym przypadku zmniejsza się w przybliżeniu dwukrotnie.
7.10.2. Metoda ortogonalizacji
Wtej części zajmiemy się metodą, która omija trudności pojawiające sie przy wykorzystaniu równań normalnych,
ponieważ w ogóle ich nie konstruuje. Metoda opiera sięnarozkładzie QR, który jest charakteryzowane przez poniższe:
TWIERDZENIE 82. Niech A będzie macierzą prostokątną rozmiaru m Ł n, gdzie m � n, której kolumny są
liniowo niezależne. Istnieje wtedy jedyny rozkład A = QR na czynniki: macierz Q o rozmiarach m Ł n taką, że
QT Q = D, gdzie D jest macierzą diagonalną o wszystkich elementach dodatnich i macierz trójkątna górną R z
elementami rkk =1; k =1; :::; n.
Zilustrujmy zastosowanie rozkładu QR do zagadnienia poszukiwania rozwiązania średniokwadratowego. Jego włas-
ność (147) można teraz zapisać w postaci:
RT QT (b Ą Ax) =0:
Ponieważ macierz R jest nieosobliwa i QT A = QT QR = DR, topowyższy układ jest równoważny układowi
Rx = y, gdzie y = DĄ1QT b:
Zatem jeśli znane będą macierze Q i R, to rozwiązanie średniokwadratowe otrzymamy rozwiązując układ z macierzą
trójkątną. Pozostaje tylko znalez
ć odpowiedni algorytm rozkładu, który pozwoli obliczyć Q, R i y. Można to zrobić
np. za pomocą zmody
�kowanej metody Grama-Schimdta.
14
7.11. Macierze rzadkie
7.11.1. Macierze liczbowe
DEFINICJA 83. Macierz A nazywa się wstęgową, j�sli wjej wszystkie niezerowe elementy znajdują s/
e wewnątrz
obszaru ograniczonego przekątnymi równoległymi do głownej przekątnej. W ten sposób aij =0, jeśli i Ą j > Ż lub
j Ą i >�i ak;kĄŻ =0 i ak;k+� =0 chociaż dla jednej wartości k. Ż + � +1 jest szerokościąwstęgi. Wstęgą macierzy
6 6
A nazywa się zbiór elementów, dla których i Ą j � Ż lub j Ą i � �. Innymi słowami: dla dowolnego i-tego wiersza do
wstęgi należą wszystkie elementy o indeksach kolumnowych od i Ą Ż do i + �, razemŻ + � +1elementów.
UWAGA: Jeśli macierz jest symetryczna, to wystarczy oczywiście przechowywać jej półwstęgę. Np. dolna półw-
stęga składa się z elementów znajdujących się wgórnej części wstęgi, tj. takich że 0 W schemacie diagonalnym dla przechowywania symetrycznej macierzy wstęgowej A rozmiaru nŁn wykorzystuje
się tablicę o rozmiarze n Ł (Ż +1). Główną przekątną przechowuje się w ostatniej kolumnie, a dolne przekątne - w
pozostałych kolumnach z przesunięciem o jedną pozycję wdół przy każdym przejściu w lewo.
PRZYKA
AD 1. Roz�azmy symetryczną macierz wstęgową z szerokóscią wstęgi 5, dla przexchowania której
potrzebna jest tablica rozmiaru 7 Ł 3:
2 3 2 3
1 1
6 7 6 7
2 8 9 0 2
6 7 6 7
6 7 6 7
8 3 0 8 3
6 7 6 7
6 7 6 7
A = 9 4 10 AN(i; j) = 9 0 4
6 7 6 7
6 7 6 7
10 5 11 12 0 10 5
6 7 6 7
4 5 4 5
11 6 0 11 6
12 7 112 0 7
W przypadku macierzy niesymetrycznych konieczna jest oczywiście tablica o rozmiarach n Ł (Ż + � +1), wtedy
dolną półwstęgę i główna przekątną przechowuje się tak, jak poprzednio, natomiast górne przekątne - w prawej częci
tablicy z przesunięciem o jedną pozycję w górę przykażdym przejściu w prawo.
UWAGA: (i) Stosowanie schematu diagonalnego jest opłacalne, jeśli max fŻ; �g żn.
(ii) Schemat diagonalny jest schematem łatwego dostępu, w tym sensie, że istnieje proste wzajemnie jednoznaczne
odwzorowanie między położeniem elementu w macierzy A i w tablicy AN: aij jest przechowywany w AN(i; j Ą i +
Ż +1) (dla macierzy symetrycznej)
Macierze wstęgowe majątąważnąwłasność, że szerokość wstęgi zależy od porządku, w jakim położone sąwiersze
i kolumny. Dlatego można szukać przestawienia wierszy i kolumn, prowadzącego do zmniejszenia szerokości wstęgi.
Mniejsza szerokość wstęgi oznacza oczywiście zmniejszenie rozmiaru tablicy AN, co oznacza mniejsze wymagania
pamięci i zmniejsze kosztu obliczeń. Symetria macierzy bedzie zachowana, jeśli dla wierszy i kolumn wykonane będzie
takie same przestawienie.
W przypadku rozwiązywania układu równań liniowych z macierzą wstęgową metodami eliminacji bez wyboru ele-
mentu głównego dowolny nowy element niezerowy, jeśli się pojawi, to tylko wewnątrz wstęgi. Również przy obliczaniu
wartości i wektorów własnych nie potrzeba dodatkowej pamięci, oprócz tablicy AN.
Macierze wstęgowe dużego rozmiaru mogą mi�c szeroką wstęgę i dużą liczbę zer wewnątrz niej. Dla takich
macierzy schemat diagonalny może być nieekonomiczny. W takim przypadku, bardziej efektywny jest schemat prze-
chowywania macierzy symetrycznych nazywany pro
�lowym lub zmiennej wstęgi.
Dla każdego i-tego wiersza macierzy symetrycznej A przyjmijmy
Żi = i Ą jmin(i);
gdzie jmin(i) jest najmniejszym indeksem kolumnowym i-tego wiersza, dla którego aij =0. W ten sposób, pierwszy
6
niezerowy element i-tego wiersza znajduje sięnapozycji Żi na lewo od głównej przekątnej. Zatem półszerokośćwstęgi
Ż = max Żi:
i=1;:::;n
Powłoka macierzy A jest zbiorem elementów aij, dla których 0 wszystkie elementy z indeksami kolumnowymi od jmin(i) do i Ą 1, razem Żi elementów (elementów przekątniowych
nie zalicza się dopowłoki). Pro
�l macierzy A określa się jako liczbę elementówwpowłoce
n
X
prof(A) = Żi:
i=1
15
Wszystkie elementy powłoki, uporządkowane według wierszy, przechowuje się, włączając zera, w jednowymi-
arowej tablicy AN. Element przekątniowy danego wiersza umieszcza się za jego elementami. Rozmiar tablicy AN
jest równa sumie pro
�lu i stopnia macierzy A. Konieczna jest jescze jedna jednowymiarowa tablica IA wskaz
ników,
której elementy wskazują położenie elementów przekątniowych w tablicy AN (rozmiar tablicy IA jest oczywiście
równy stopniowi macierzy A). Zatem, dla i >1 elementy i-tego wiersza znajdują się w tablicy AN na pozycjach od
IA(i Ą 1) + 1 do IA(i). Innymi słowy
aij = AN(IA(i) Ą i + j):
PRZYKA
AD 2: Rozważmy macierz wstęgową z przykładu 1. Wówczas
Ł ń
AN = 1 2 8 3 9 0 4 10 5 11 6 12 0 7
Ł ń
IA = 1 2 4 7 9 11 14
UWAGA: (i) Możliwe są jeszcze inne warianty schematu pro
�lowego: przechowywanie powłoki kolumnami, prze-
chowywanie elementów głownej przekątnej w oddzielnej tablicy.
(ii) Jak w przypadku schematu diagonalnego, można osiągnąć w/
ekszą efektywność przestawiając wiersze i kolumny,
wówczas zmniejszy się pro
�l macierzy. Istnieją bardzo skuteczne algorytmy minimalizacji pro
�lu dla macierzy rzad-
kich, np. algorytm odwrotny Cuthill-McKee wykorzystujący grafy.
(iii) Nowe elemnty niezerowe, w przypadku stosowania metod eliminacji, mogą pojawić się tylko wewnątrz powłoki,
gdzie jest już dla nich zarezerwowane miejsce.
(iv) Schemat pro
�lowy jest bardzo przydatny przy rozwiązywaniu problemów statycznych, czyli takich, w których
macierz nie zmienia swojej postaci.
Teraz zostanie omówiony tzw. schemat rozrzedzony wierszowy, który ma najmniejsze wymagania pamięciowe i
dlatego jest bardzo przydatny dla wielu obliczeń związanych z macierzami rzadkimi, niekoniecznie symetrycznymi.
Wartości niezerowych elementów macierzy i odpowiednie indeksy kolumnowe przechowuje się według wierszy w
dwóch jednowymiarowych tablicach (AN i JA). Dodatkowo, niezbędna jest jeszcze jedna jednowymiarowa tablica
IA wskaz
ników, okreslających pozycje w tablicach AN i JA, na których znajdują s/
e pierwsze elementy danego
wiersza. Dodatkowy element tablicy IA zawiera wskaz
nik do pierwszej wolnej pozycji w JA i AN.
PRZYKA
AD 3: Rozważmy macierz
2 3
0 0 1 2 0 0 0 4 0 0
4 5
A = 0 0 0 0 0 0 0 0 0 0 :
0 0 0 0 0 3 0 5 0 0
Przechowamy ją wnastępujący sposób:
Ł ń
AN = 1 2 4 3 5
Ł ń
JA = 3 4 8 6 8
Ł ń
IA = 1 4 4 6
Opis pierwszego wiersza zaczyna się na pozycji IA(1) = 1 tablic AN i JA, drugiego wiersza od IA(2) = 4, z�s
trzeciego od IA(3) = 4. IA(4) = 6 podaje pierwszą wolną pozycję.
W ogólnym przypadku, opis i-tego wiersza macierzy A znajduje się w tablicach JA i AN na pozycjach od IA(i)
do IA(i +1) Ą 1, z wyłaczeniem sytuacji, gdy IA(i) =IA(i +1), która oznacza, że i-ty wiersz jest pusty (tzn. nie
zawiera elementów niezerowych).
Powyższą wersję nazywa się pełną, jesli przedstawia całą macierz A, i uporządkowaną, j�sli elementy każdego
wiersza przechowuje się odpowiednio ze wzrostem indeksów kolumnowych, tzw. RR(C)O (Row-wise Representation
Complete and Ordered):
Inne warianty schematu rozrzedzonego:
1. RR(C)U (wierszowy, nieuporządkowany) - uporządkowanie wierszy jest zachowane, ale wewnątrz każdego
wiersza elementy są przechowywane w dowolnym porządku. W przypadku macierzy z przykładu 3:
Ł ń
AN = 4 1 2 5 3
Ł ń
JA = 8 3 4 8 6
Ł ń
IA = 1 4 4 6
Nieuporządkowanie może być przydatne w sytuacji, gdy w macierzy pojawiają s/
e nowe elementy niezerowe. Nie
wystepuje wtedy zwiększenie kosztu związane z potrzebąporządkowania elementów w obrebie opisu każedgo wiersza.
16
2. CR(C)O (kolumnowy, uporządkowany) - przechowuje sie wtedy wartości oraz indeksy wierszowe według kol-
umn. Dla przykładu 3 wówczas
Ł ń
AN = 1 2 3 4 5
Ł ń
JA = 1 1 3 1 3
Ł ń
IA = 1 1 1 2 3 3 4 4 6 6 6
Zauważmy, że jest to jednoczesnie schemat RR(C)O dla macierzy AT .
Jesli macierz rzadka A jest symetryczna, to oczywiście wystarczy przechowywaćtylkojej główną przekątną i jeden
z trójkątów. Mianowicie:
PRZYKA
AD 4: Rozważmy macierz
2 3
12
6 7
3 4
6 7
B = ;
4 5
5
2 4 6
wówczas przechowując przekątną i górny trójkąt w wariancie RR(DU)O (wierszowy, przekątniowy i górny, uporzad-
kowany), otrzymamy
Ł ń
BN = 1 2 3 4 5 6
Ł ń
JB = 1 4 2 4 3 4
Ł ń
IB = 1 3 5 6 7
UWAGA: Oczywiście eleemnty głownej przekątnej w przypadku macierzy symetrycznych można przechowywaćw
oddzielnej tablicy i w wielu wypadkach jest to bardziej efektywne, choć czasami może komplikować prostyalgorytm.
PRZYKA
AD 5: Macierz z przykładu 4 w wariancie z oddzielnie przechowywaną głowną przekątną (RR(U)O -
wierszowy, górny, uporządkowany):
Ł ń
BD = 1 3 5 6
Ł ń
BN = 2 4
Ł ń
JB = 4 4
Ł ń
IB = 1 2 3 3 3
Zajmiemy się teraz zwartym schematem przechowywania macierzy rzadkich trójkątnych, który jest w pewnym
sensie wariantem schematu rozrzedzonego czyli tzw. schematem Shermana. Omówimy wersję dla macierzy trójkątnej
górnej U.
Aby uniknąc powtórnego zapisu jednakowych ciągów indeksów kolunowych, w zwarty sposób przechowuje się te
indeksy wykorzystując jeszcze jedną jednowymiarową tablicę wskaz
ników. Elementy głównej przekątnej przechowuje
się w tablicy UD, z�s elementy trójkąta według wierszy w tablicy UN. Tablica IU zawiera wskaz
niki dla UN, jak
poprzednio. Indeksy kolumnowe, odpowiadające znajdującym się w UN eleemntom niezerowym, przechowuje się w
tablicy JU w zwartej formie, wykorzystująć to, że pewne zbiory wierszy mają podobną strukturę (a dokładniej wiersze
o indeksach nie większych od j mają takąsamą strukturę na prawo od kolumny j). W ten sposób między UN i JU nie
istnieje proste odwzorowanie. Indeksy kolumnowe wybiera się z JU za pomocą dodatkowej tablicy wskaz
ników IJ.
Pewne wiersze mogą być opisane w JU tym samym zestawem indeksów kolumnowych.
PRZYKA
AD 6: Rozważmy macierz trójkątna górną
2 3
1 11 12 13
6 7
214 15
6 7
6 7
316 17
6 7
6 7
418
6 7
6 7
5 19 20
6 7
U = ;
6 7
6 21 22 23
6 7
6 7
7 24 25
6 7
6 7
8 26
6 7
4 5
9 27
10
17
która będzie za pomocą schematu Shermana przechowywana następująco:
Ł ń
UD = 1 2 3 4 5 6 7 8 9 10
Ł ń
UN = 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
Ł ń
IU = 1 4 6 8 9 11 14 16 17 18 18
Ł ń
JU = 3 7 10 6 9 8 9 10 8 9 10
Ł ń
IJ = 1 4 2 5 2 6 9 10 11
Zauważmy, że np. wiersze 1, 3 i 5 mają jednakową strukturę po prawo od kolumny 5. Wiersz 1 opisują indeksy 3, 7,
10 i jednocześnie indeksy 7 i 10 opisują wiersze 3 i 5. W ten sposób: IJ(1) = 1, IJ(3) = IJ(5) = 2. J�sli chcemy
wybrać pełny wiersz, np. pierwszy, to postępujemy następująco: wartości elementów znajdują się w UN na pozycjach
od IJ(1) = 1 do IJ(2) Ą 1 =3, zaś ich indeksy kolumnowe znajdują się w JU na tych samych pozycjach.
W ogólnym przypadku, opis i-tego wiersza macierzy U znajduje się w tablicach UJ i UN na pozycjach od IJ(i)
do IJ(i +1)Ą 1, z wyłączeniem sytuacji, gdy IJ(i) >IJ(i +1)Ą 1, która oznacza, że i-ty wiersz jest opisany w UJ
i UN na pozycjach od IJ(i) do pozycji, na której w JU kończy się wzrastanie indeksów.
7.11.2 Macierze blokowe
Idea rozbijania wielkich macierzy na podmacierze (bloki) powstaje ze względu na fakt, że można, posługiwać s/
e
nimi tak, jak byłyby elementami macierzy. Również w technologii macierzy rzadkich ma to duże znaczenie, gdyż wiele
algorytmów skonstruowanych dla macierzy liczbowych mozna zastosować dla przypadku macierzy blokowych. Przez
przechowywanie macierzy blokowych rozumiemy przechowywanie zbioru jej podmacierzy.
Omówimy najpierw schemat niejawny przeznaczony przede wszystkim dla macierzy symetrycznych z kwadra-
towymi blokami przekątniowymi. Bloki przekątniowe rozpatruje się tak, jak gdyby tworzyły jedną macierz i prze-
chowuje się za pomocą schematu pro
�lowego (wartości elementów - w tablicy AN i wskaz
niki do eleemntów diago-
nalnych każdego wiersza - w tablicy IA). Bloki podprzekątniowe ropzatruje się tak, jak gdyby tworzyły druga macierz
i przechowuje się za pomocą schematu rozrzedzonego wierszowego (elementy niezerowe - w tablicy AP, odpowiednie
indeksy kolumnowe - w tablicy JA, wskaz
niki do początku wierszy - w tablicy IP ). Samo rozbicie określa się nu-
merami pierwszych wierszy każdego bloku, przechowywanymi w jednowymiarowej tablicy LP , gdzierównieżdołacza
się dodatkowy skłądnik oznaczający pierwsza wolna pozycję w AP i JA.
UWAGA: Schemat niejawny jest bardzo przydatny dla metod eliminacji.
PRZYKA
AD 7:Rozważmy symetryczną macierz blokową:
2 3
1 3 17
6 7
2 4 6 21
6 7
6 7
3 4 5 18 24
6 7
6 7
6 7 19 20
6 7
6 7
17 8 10 23
6 7
A = ;
6 7
18 19 9 22
6 7
6 7
20 10 11 25
6 7
6 7
21 22 12 13 15
6 7
4 5
23 13 14
24 25 15 16
której wierszei kolumnysazgrupowanezapomocaindeksów 1Ą4, 5Ą7 i 8Ą10. W rezultacie A zostaje rozdzielona na
9 bloków: przekątniowe są kwadratowe, gdyż wiersze i kolumny zostały pogrupowane w ten sam sposób. Zastosowanie
schematu niejwabego daje następujące tablice:
Ł ń
AN = 1 2 3 4 5 6 0 7 8 9 10 0 11 12 13 14 15 0 16
Ł ń
IA = 1 2 5 8 9 10 13 14 16 19
Ł ń
AP = 17 18 19 20 21 22 23 24 25
Ł ń
JA = 1 3 4 4 2 6 5 3 7
Ł ń
IP = 1 1 1 1 1 2 4 5 7 8 10
Ł ń
LP = 1 5 8 11
18
Inny schemat, przeznaczony zwłaszcza dla macierzy, które mają wiele bloków zerowych, jest nazywany hiper-
macierzowym. Używa się w nim dowolnego zwykłego schematu, różnica polega tylko na tym, że zamiast wartości
liczbowych elementów przechowuje się informację o położeniu i rozmiarach bloków. Natomiast same bloki prze-
chowywane są w zależności od ich formatu. Tego typu metody pozwalają na rozwiązywanie wielkich zadań przy
niewielkich wymaganiach pamięci, a zwłaszcza na wygodne wykonywanie operacji na macierzach takich, jak do-
dawanie, mnożenie przez liczbę, przestawianie, itp.
19