Wykład 7: Podstawowe metody numeryczne dla
Wykład 7: Podstawowe metody numeryczne dla
liniowych układów algebraicznych
liniowych układów algebraicznych
Rozważmy układ n równań liniowych z niewiadomymi x1, x2,..., xn
ìð
E1 : a1,1x1 +ð a1,2x2 +ð...+ð a1,n-ð1xn-ð1 +ð a1,nxn =ð b1
ïð
E2 : a2,1x1 +ð a2,2x2 +ð...+ð a2,n-ð1xn-ð1 +ð a2,nxn =ð b2
ïð
ïð
íð
ïð
En-ð1 : an-ð1,1x1 +ð an-ð11,2x2 +ð...+ð an-ð1,n-ð1xn-ð1 +ð an-ð1,nxn =ð bn-ð1
ïð
ïð
En : an,1x1 +ð an,2x2 +ð...+ð an,n-ð1xn-ð1 +ð an,nxn =ð bn
îð
W notacji macierzowo wektorowej powyższy układ równań ma postać
Ax =ðb
a1,1 a1,2 ... a1,n-ð1 a1,n
éðÅ‚ð
Ä™ð
a2,1 a2,2 ... a2,n-ð1 a2,n Å›ð x1 b1
éð Å‚ð éð Å‚ð
Ä™ðÅ›ð
Ä™ð Å›ð Ä™ð Å›ð
... ... ... ... ...
gdzie A =ð Ä™ðÅ›ð , x =ð , b =ð
Ä™ð Å›ð Ä™ð Å›ð
Ä™ða
an-ð1,2 ... an-ð1,n-ð1 an-ð1,n Å›ð
Ä™ð Å›ð Ä™ð Å›ð
n-ð1,1
ëðxn ûð ëðbn ûð
Ä™ðÅ›ð
an,1 an,2 ... an,n-ð1 an,n ûð
Ä™ðÅ›ð
ëð
Zakładamy, że rozwiązanie istnieje i jest jednoznaczne, tj. wyznacznik det(A) jest niezerowy.
Omówimy ogólną metodę rozwiązania takiego układu znaną pod nazwą Metody Eliminacji
Gaussa (MEG). Podstawowa koncepcja tej metody jest prosta: manipulujÄ…c w odpowiedni
sposób równaniami sprowadzimy wyjściowy układ do równoważnego (tj. posiadającego to
samo rozwiązanie) układu z macierzą górną trójkątną, a następnie w łatwy sposób
rozwiążemy ten układ.
Rozważmy pierwszy krok procedury eliminacji. W kroku tym eliminujemy niewiadomą x1 z
równań E2, E3,.., En. W tym celu równanie E1 jest mnożone przez odpowiednio dobrane
liczby, a następnie jest odejmowane od kolejnych równań układu. Transformacja ta wygląda
następująco:
E1 : a1,1x1 +ð a1,2x2 +ð ... +ð a1,nxn =ð b1
ìð
E1 : a1,1x1 +ð a1,2x2 +ð ... +ð a1,nxn =ð b1
ìð
ïð
ïðE : a2,1x1 +ð a2,2x2 +ð ... +ð a2,nxn =ð b2
(1) (1) (1) (1)
a2,2x2 +ð ... +ð a2,nxn =ð b2
ïð ïðE :
2 2
Þð
íð íð
ïð ïð
(1) (1) (1) (1)
ïðEn : an,1x1 +ð an,2x2 +ð ... +ð an,nxn =ð bn ïðEn :
an,2x2 +ð ... +ð an,nxn =ð bn
îð
îð
ai,1
ai(1) =ð ai, j -ð a1, jli,1 , li,1 =ð
, j
a1,1
gdzie
bi(1) =ð bi -ð b1li,1 , i, j =ð 2,...,n
Macierzowo-wektorowa postać układu liniowego otrzymanego w efekcie pierwszego etapu
eliminacji to A(1)x =ð b(1), gdzie
a1,1 a1,2 a1,n
b1
éðÅ‚ð
éð Å‚ð
(1) (1)
Ä™ð (1)
Ä™ðb Å›ð
0 a2,2 a2,n Å›ð
2
Ä™ðÅ›ð
Ä™ð Å›ð
A(1) =ð b(1) =ð
Ä™ðÅ›ð
Ä™ð Å›ð
Ä™ð (1) (1)
Ä™ðb(1) Å›ð
0 an,2 an,n Å›ð
ëð n ûð
ëðûð
Procedura eliminacji kolejnych niewiadomych jest kontynuowana analogicznie aż do
osiągnięcia docelowej formy układu. W szczególności, po (k-1)-szym kroku układ równań ma
następującą formę
E1 : a1,1x1 +ð a1,2x2 +ð ... +ð a1,k xk +ð ... +ð a1,nxn =ð b1
ìð
Górny indeks informuje ile razy
(1) (1) (1) (1) (1)
ïð
E2 : a2,2x2 +ð ... +ð a2,k xk +ð ... +ð a2,nxn =ð b2
modyfikowana była dana wielkość do
ïð
tego momentu obliczeń. W
ïð
następnym, tj. k-tym kroku, równanie
ïð
( ( ( -ð1) (
Ekk-ð1) : akkk-ð1)xk +ð ... +ð akkn xn =ð bkk-ð1) (
íð
,,
Ekk-ð1) bÄ™dzie wykorzystane do
( -ð ( ( ( -ð
ïð
Ekk11) : akk-ð1)xk +ð ... +ð akk-ð1)xn =ð bkk11)
+ð +ð1,k +ð1,n +ð
eliminacji niewiadomej xk z równań
ïð
( -ð (
Ekk11),..., Enk-ð1).
ïð
+ð
( ( ( -ð1) (
ïð
Enk-ð1) : ankk-ð1)xk +ð ... +ð ankn xn =ð bnk-ð1)
, ,
îð
Formuły transformacji współczynników macierzy układu oraz elementów wektora prawych
stron układu w kroku k-tym mają postać
ai(,k-ð1)
( -ð1) k
ai(,kj) =ð ai(,kj-ð1) -ð akkj li,k , li,k =ð
,
(
akkk-ð1)
,
(
bi(k ) =ð bi(k ) -ð bkk-ð1)li,k , i, j =ð k +ð1,...,n
Macierzowo-wektorowa postać ukÅ‚adu po k krokach eliminacji to A(k )x =ð b(k )
, gdzie
éðÅ‚ð
a1,1 a1,k+ð1 a1,n
Ä™ðÅ›ð
*ð
Ä™ðÅ›ð
(1) ( ( ) (
( ) ( )
Ä™ðÅ›ð
, b(k ) =ð [b1,b2 ,..,bkk-ð1),bkk1,..,bnk ) ]T
A(k ) =ð *ð *ð akk1,k+ð1 akk1,n
+ð
+ð +ð
Ä™ðÅ›ð
*ð *ð
Ä™ðÅ›ð
( (
Ä™ðûð
*ð *ð ankk)+ð1 ankn) Å›ð
, ,
ëð
Procedura eliminacji kończy się po n-1 krokach otrzymaniem układu równań o następującej
strukturze
ìð
E1 : a1,1x1 +ð a1,2x2 +ð ... +ð a1,k xk +ð ... +ð a1,nxn =ð b1
ïð
(1) (1) (1) (1) (1)
a2,2x2 +ð ... +ð a2,k xk +ð ... +ð a2,nxn =ð b2
ïðE :
2
ïð
íð
ïðE(n-ð2) ( -ð ( -ð (n-ð2)
: ann1,2)-ð1xn-ð1 +ð ann1,2)xn =ð bn-ð1
n-ð1 -ð n -ð n
ïð
( ( -ð1) (
ïðEnn-ð1) :
annn xn =ð bnn-ð1)
,
îð
Ć
W zapisie macierzowo-wektorowym mamy Ux =ð b, gdzie
a1,1 a1,2 a1,n-ð1 a1,n
éðÅ‚ð
(1) (1) (1)
Ä™ð
*ð a2,2 a2,n-ð1 a2,n Å›ð
Ä™ðÅ›ð
(1) ( -ð2) (
Ć
*ð *ð
U ºð A(n-ð1) =ð Ä™ðÅ›ð , b ºð b(n-ð1) =ð [b1,b2 ,..,bnn1 ,bnn-ð1) ]T
-ð
Ä™ð ( -ð ( -ð
*ð *ð *ð ann1,2)-ð1 ann1,2) Å›ð
-ð n -ð n
Ä™ðÅ›ð
( -ð1)
*ð *ð *ð *ð annn ûð
Ä™ðÅ›ð
,
ëð
Otrzymany układ z macierzą górną trójkątną U można rozwiązać łatwo  od końca .
Istotnie, ostatnie równanie układu zawiera tylko ostatnią niewiadomą, którą można
natychmiast obliczyć. Wówczas przedostatnie równanie jest efektywnie równaniem również z
jedną niewiadomą (czyli xn-1) ponieważ xn jest już znana. Ogólnie, możemy zapisać
następujący algorytm rozwiązania układu z macierzą U następująco
( ( -ð1)
xn =ð bnn-ð1) / annn
,
n
éðÅ‚ð
1
(i-ð1) -ð1)
xi =ð -ð ai(,ij xj Å›ð , i =ð n -ð1,n -ð 2,...,2,1
åð
ai(,ii-ð1) Ä™ðbi
j=ði+ð1
ëðûð
Zauważmy, że pętla realizowana jest w kierunku malejącego indeksu (wstecz). Pseudokod
kompletnego algorytmu eliminacji Gaussa przedstawiamy na następnej stronie.
Metoda Eliminacji Gaussa(wariant podstawowy )
% Etap 2 -ð Rozwiazanie
% Etap 1 -ð Transformacja do ukladu trojkatnego
x(n) Źð b(n) / a(n, n)
for k =ð 1 : n -ð 1
for j =ð n -ð 1 : 1 : -ð1
for j =ð k +ð 1: n
x( j) Źð b( j)
l( j, k ) Źð a( j, k ) / a(k, k )
®ð for k =ð j +ð 1: n
a( j, k +ð 1 : n) Źð a( j, k +ð 1 : n) -ð l( j, k ) ×ð a(k, k +ð 1 : n)
x( j) Źð x( j) -ð a( j, k ) ×ð x(k )
b( j) Źð b( j) -ð l( j, k ) ×ð b(k)
end
end
x( j) Źð x( j) / a( j, j)
end
end
UWAGA: algorytm w opisanej formie może zawieść nawet jeśli został zastosowany do
macierzy nieosobliwej (det(A) `" 0). O tym dalej &
Faktoryzacja LU
Pokażemy, że ważnym  produktem ubocznym eliminacji Gaussa rozkład wyjściowej
macierzy układu do postaci iloczynu dwóch macierzy: dolnej i górnej trójkątnej.
Rozkład taki nazywany jest faktoryzacją (przedstawieniem w postaci iloczynu) LU.
(1) (1) (1)
Rozważmy ponownie 1-szy krok eliminacji. Wprowadziliśmy uprzednio liczby l2,1 ,l3,1 ,...,ln,1 .
Użyjemy ich teraz do skonstruowania specjalnej macierzy L(1) , a mianowicie
1 *ð *ð *ð
éðÅ‚ð
(1)
Ä™ð
Ä™ð-ðl2,1 1 *ð *ðÅ›ð
Å›ð
L(ð1)ð =ð
*ð *ð
Ä™ðÅ›ð
Ä™ð (1) Å›ð
n,1
ëð-ðl *ð *ð 1ûð
Symbol  gwiazdka oznacza elementy zerowe. Zauważmy teraz (ćwiczenie dla Czytelnika!),
że pierwszy etap eliminacji transformuje oryginalny układ równań do układu z macierzą
A(1) =ð L(1)A!
l(kk) , j =ð k +ð1,..,n za pomocÄ…
Ogólnie, w k-tym etapie eliminacji zdefiniowaliśmy liczby ,
j,
których możemy skonstruować macierz L(k ) postaci
1 *ð *ð *ð *ð *ð
éðÅ‚ð
Ä™ð*ð
*ð *ð *ð *ðÅ›ð
Ä™ðÅ›ð
*ð *ð 1 *ð *ð *ð
Ä™ðÅ›ð
L(k ) =ð
Ä™ð*ð *ð -ðlkk1,k 1 *ð *ðÅ›ð
( )
+ð
Ä™ðÅ›ð
*ð *ð *ð *ð
Ä™ðÅ›ð
Ä™ð*ð *ð -ðlnkk) *ð *ð 1Å›ð
(
,
ëð ûð
Ponownie, macierz współczynników układu po k-tym kroku eliminacji może być obliczona ze
wzoru A(k ) =ð L(k )A(k-ð1).
Z powyższego wynika natychmiast, że docelowa macierz górna trójkątna U może być
przedstawiona jako iloczyn oryginalnej macierzy A oraz macierzy typu L z kolejnych kroków
eliminacji, a mianowicie
U ºð A(n-ð1) =ð L(n-ð1)A(n-ð2) =ð L(n-ð1)L(n-ð2)...L(2)L(1) A =ð LA
L
Kapitalną własnością macierzy L jest to, że macierz do niej odwrotna wyraża się zaskakująco
prosto, a mianowicie
1 *ð *ð *ð *ð *ð
éðÅ‚ð
(1)
Ä™ðl
1 *ð *ð *ð *ðÅ›ð
2,1
Ä™ðÅ›ð
(1) (2
l3,1 l3,2) *ð *ð *ð
Ä™ðÅ›ð
L ºð L-ð1 =ð
Ä™ð
1 *ð *ðÅ›ð
Ä™ðÅ›ð
( -ð
lnn1,2) 1 *ð
Ä™ðÅ›ð
-ð n-ð2
Ä™ðl ln2)
(1) ( ( (
lnnn-ð22) lnnn-ð1) 1Å›ð
n,1 ,2 , -ð , -ð1
ëðûð
Dowód tego (nieoczywistego) faktu można znalez
ć w dobrych podręcznikach numerycznej
algebry liniowej np. w znakomitej książce Trefethena i Bau pt.  Numerical Linear Algebra .
Dla ambitnych godna polecenia jest również monografia  Analiza Numeryczna Kincaida i
Cheneya (WNT, Warszawa 2006).
Otrzymaliśmy zatem następującą formułę faktoryzacyjną dla macierzy (nieosobliwej!) A
A =ðLU
Pseudokod podstawowej wersji faktoryzacji LU wygląda następująco:
Faktoryzacja LU -ð wariant podstawowy
for k =ð 1: n -ð1 do
ìð
ïð
for j =ð k +ð 1: n do
ïð
a( j, k) Źð a( j, k) / a(k, k)
ïð
íð
a( j, k +ð1: n) Źð a( j, k +ð1: n) -ð a( j, k) ×ð a(k, k +ð1: n)
ïð
ïð
end
ïð
îðend
Uwaga :
1. Oryginalna zawartosc A zostala zniszczona!
2. Glowna diagonala i gorny trojkat A zawiera
niezerowe elementy macierzy U
3. Dolny trojkat A zawiera niezerowe elementy
macierzy L. Nie ma potrzeby przechowywania
elementow diagonalnych L (sÄ… to jedynki).
Co można powiedzieć o koszcie numerycznym MEG i faktoryzacji LU? Analizując
pseudokody tych algorytmów można wyciągnąć następujące wnioski:
1. W etapie eliminacji mamy trzy zagnieżdżone pętle. Wynika z tego, że całkowita liczba
operacji zmiennoprzecinkowych jest proporcjonalna do n3, gdzie n jest rozmiarem
układu. Dokładniejsza analiza pokazuje, że dla dużych wartości n dominujący koszt
wynika z operacji mnożenia i dodawania/odejmowania oraz, że liczba tych operacji jest
2
bliska n3.
3
2. W etapie rozwiązania układu z macierzą górną trójkątną mamy jedynie dwie
zagnieżdżone pętle, wobec czego koszt numeryczny będzie (z grubsza) proporcjonalny
n2. Wnioskujemy, że dla układu o dużym rozmiarze koszt numeryczny etapu
rozwiązania jest znikomy w porównaniu z kosztem sprowadzenia tego układu do
postaci trójkątnej.
3. Koszt numeryczny faktoryzacji LU (zastosowanej do ogólnej macierzy kwadratowej
macierzy) jest proporcjonalny do trzeciej potęgi jest rozmiaru.
Kiedy opłaca się zastosować faktoryzację LU zamiast  zwyczajnej MEG? Rozważmy
sytuację, w której chcemy wyznaczyć rozwiązania szeregu m układów liniowych o tej samej
macierzy A, ale różnych wektorach prawych stron
Ax( j) =ð b( j) , j =ð 1,..,m
Używając MEG, całkowity koszt numeryczny rozwiązania tych układów będzie
proporcjonalny ( z grubsza) do m×ðn3. Taka estymacja kosztu wynika z faktu, że każdorazowo
przeprowadzana jest kompletna procedura eliminacji (kosztem n3) pracujÄ…cÄ… na tej samej
macierzy. Niepotrzebnego powtarzania tych samych operacji można jednak uniknąć stosując
faktoryzację LU. Idea polega na jednorazowym obliczeniu faktorów U i L (kosztem n3), a
następnie rozwiązaniu m układów równań kosztem proporcjonalnym jedynie do n2 (oba
faktory są trójkątne). Oto jak to działa:
1) Wyznacz faktory L i U macierzy A (koszt n3)
2) dla j =ð 1,..,m wykonaj
ìð
Ly =ð b( j)
LUx( j) =ð b( j) Þð
íð
( j)
îðUx =ð y
Pseudokod procedury rozwiązania układu liniowego przy pomocy faktorów L i U
Rozwiazanie ukladu (LU)x = b
% Ux =ð y
% Ly =ð b
x(n) Źð y(n) / a(n, n)
ìð
y(1) Źð b(1)
ìð
ïð
for j =ð n -ð 1 : 1 : -ð1
ïð
ïð
for j =ð 2: n
ïð
x( j) Źð y( j)
ïð
y( j) Źð b( j)
ïð
ïð
Þð for k =ð j +ð 1: n
ïð ïð
for k =ð 1 : j -ð 1
íð íð
x( j) Źð x( j) -ð a( j, k ) ×ð x(k )
ïð ïð
y( j) Źð y( j) -ð a( j, k ) ×ð y(k )
ïð ïð
end
end
ïð ïð
x( j) Źð x( j) / a( j, j)
ïðend ïð
îð
ïðend
îð
Trudności z podstawowym wariantem MEG. Wybór elementu głównego.
Przedstawiony wyżej podstawowy wariant MEG może zawieść nawet gdy macierz A jest
(
nieosobliwa. Istotnie, jeÅ›li w k-1-szym etapie eliminacji z równania Ekk-ð1) zniknęła również
(
niewiadoma xk , to współczynnik akkk-ð1) jest równy zeru i w kroku k-tym nie mam możliwoÅ›ci
,
obliczenia żadnej z liczb l(k) (mamy dzielenia przez zero)!
W skrajnym przypadku, pierwsze z równań może w ogóle nie zawierać niewiadomej x1 i
proces eliminacji nie może być rozpoczęty! Oto banalny przykład nieodpowiedniej macierzy
0 1
éðÅ‚ð
A =ð
Ä™ð1 1Å›ð
ëðûð
To jednak nie koniec problemów. Rozważmy mianowicie układ równań z następującą
macierzą współczynników
éðÅ‚ð
10-ð10 1
A =ð
Ä™ð
1 1Å›ð
ëðûð
1 0
éðÅ‚ð
10-ð10 1
éðÅ‚ð
Jej faktory L i U to (sprawdzić  ćwiczenie!): L =ð
Ä™ðÅ›ð
10
Ä™ð10 1Å›ð , U =ð 0 1-ð1010
ëðûð
ëðûð
Jak wiemy, arytmetyka komputera nie jest dokładna. Zobaczymy co się stanie jeśli do
rozwiązania układu liniowego z powyższą macierzą zastosujemy komputer, którego
arytmetyka jest dziesiętna i ma dokładność 7-miu miejsc znaczących (w mantysie). W uwagi
na prostotę układu obliczenia takiego hipotetycznego komputera możemy z łatwością
zasymulować & ręcznie!
Zacznijmy od tego, że element u22 faktora U będzie miał reprezentację przybliżoną, a
mianowicie
1010 -ð1 =ð 1×ð1010 -ð 0.0000000001×ð1010 =ð (1-ð 0.0000000 001) ×ð1010 @ð1010
tylko 7 cyfr jest
uwzglednionych
Zatem, dokładny faktor U będzie w naszym komputerze zastąpiony faktorem przybliżonym
postaci
éðÅ‚ð
10-ð10 1
U =ð
Ä™ð
0 -ð1010 Å›ð
ëðûð
Zauważmy, że faktor L bÄ™dzie natomiast reprezentowany Å›ciÅ›le (L ºð L).
Jeśli pomnożymy faktory macierzy A tak, jak widzi je nasz komputer to otrzymamy
éðÅ‚ð
10-ð10 1
A =ð LU =ð Ä…ð A
Ä™ð
1 0Å›ð
ëðûð
Jak widać, różnica pomiędzy A and A jest niebagatelna!
Zobaczmy jaki wpływ na rozwiązanie układu liniowego mają błędy zaokrągleń w naszym
przykładzie. Rozważmy układ równań postaci
1
éð Å‚ð
Ax =ð b , b =ð
Ä™ð0Å›ð
ëð ûð
Jego dokładne rozwiązanie to
-ð1
éðÅ‚ð
-ð(1-ð10-ð10)-ð1
éð Å‚ð
x =ðð
Ä™ð
1
(1-ð10-ð10)-ð1 Å›ð Ä™ð Å›ð
ëð ûð
ëðûð
Zastosujemy metodę faktoryzacji LU uwzględniając fakt przybliżonej reprezentacji faktora U
i błędy zaokrągleń popełnione w trakcie obliczeń. Mamy wówczas
1 1 0 y1 1 1
éð Å‚ð éð Å‚ð éð Å‚ð éð Å‚ð éð Å‚ð
L Ux =ð Þð =ð Þð y =ð
10 10
Ä™ð0Å›ð Ä™ð10 1Å›ð Ä™ð Å›ð Ä™ð0Å›ð Ä™ð Å›ð
y2
y
ëð ûð ëð ûð ëð ûð ëð ûð ëð-ð10 ûð
!!!
x1 10
éðÅ‚ð
10-ð10 1 éð Å‚ð éð Å‚ð éð Å‚ð
=ð Þð x =ð Ä…ð x
Ä™ð
Ä™ð1Å›ð
0 -ð1010 Å›ð Ä™ðx 2 Å›ð Ä™ð 10 Å›ð
ëð ûð ëð-ð10 ûð ëð ûð
ëðûð
Ponownie, otrzymane rozwiązanie znacząco odbiega od rozwiązania dokładnego! Widać
(1)
również, że  winnym zaobserwowanego efektu jest wielki co do wartości mnożnik l2,1 ,
równy aż 1010!
Okazuje się, że prostym lekarstwem na opisany problem jest zmiana kolejności równań,
czyli następująca transformacja
1 1 1 0
éðÅ‚ð
10-ð10 1
éð Å‚ð éð Å‚ð éð Å‚ð
Ć
A =ð , b =ð Þð Â =ð , b =ð
Ä™ð
-ð10
Ä™ð0Å›ð Ä™ð10 Ä™ð1Å›ð
1Å›ð
1 1Å›ð
zmiana
ëð ûð ëð ûð ëð ûð
ëðûð
kolejnosci
rownan
Poniższy rachunek pokazuje, że tym razem metoda faktoryzacji LU działa znakomicie!
1 0 1 1 1 1
éð Å‚ð éð Å‚ð éð Å‚ð
Ć
L =ð , Û =ð @ð
-ð10
Ä™ð10 Ä™ð0 1-ð10-ð10 Å›ð Ä™ð0 1Å›ð =ð Û
1Å›ð
ëð ûð ëð ûð ëð ûð
0 1 0 y1 0 0
éð Å‚ð éð Å‚ð éð Å‚ð éð Å‚ð éð Å‚ð
Ć Ć
LUx =ð Þð =ð Þð w
=ð
-ð10
Ä™ð1Å›ð Ä™ð10 Ä™ð1Å›ð
1Å›ð Ä™ð y2 Å›ð Ä™ð1Å›ð
ëð ûð ëð ûð ëð ûð ëð ûð ëð ûð
1 1 x1 0 -ð1
éð Å‚ð éð Å‚ð éð Å‚ð éð Å‚ð
Ć
=ð Þð x =ð ð x
Ä™ð0 1Å›ð Ä™ðx Å›ð Ä™ð1Å›ð Ä™ð Å›ð
1
ëð ûð ëð 2 ûð ëð ûð ëð ûð
(1)
Zauważmy, że tym razem mnożnik l2,1 =ð 10-ð10, tj. jest bardzo maÅ‚y i nie powoduje
 katastrofalnego błędu obcięcia . Ogólnie mówiąc, należy unikać dużych wartości
mnożników l!
Ogólna strategia stosowana w celu otrzymania algorytmu eliminacji odpornego na opisane
wyżej efekty nosi nazwę Metody Eliminacji Gaussa z wyborem elementu głównego.
Wybór elementu głównego w k-tym kroku eliminacji polega na:
(k
·ð znalezieniu takiego równania EEG-ð1) wÅ›ród równaÅ„ E(k -ð1), j Îð{k,k +ð1,..,n}, w którym
j
współczynnik przy niewiadomej xk jest największy co do modułu,
(k (
·ð przestawieniu w ukÅ‚adzie równania EEG-ð1) i równania Ekk -ð1) .
W ten sposób wartości bezwzględne wszystkich mnożników l nigdy nie przekraczają
jedności i katastrofalna utrata dokładności nie ma miejsca (chyba, że macierz A jest fatalnie
uwarunkowana, o czym przy innej okazji & ).
Pseudokod Metody Eliminacji Gaussa z wyborem elementu głównego
for k =ð 1: n -ð1
ìð
ïð
Znajdz i e" k taki, ze ai,k jest maksymalny x(n) Źð b(n) / a(n, n)
ìð
ïð
ïð
ïð for j =ð n -ð1: 1: -ð1
a(k, k : n) «ð a(i, k : n) przestaw wiersze i -ð ty z k -ð tym
ïð
ïð
x( j) Źð b( j)
b(k) «ð b(i) przestaw elementy wektora b ïð
ïð
ïð
ïð
for k =ð j +ð 1: n
for j =ð k +ð 1: n
ïð
Þð
íðíð
x( j) Źð x( j) -ð a( j, k) ×ð x(k)
l( j, k) Źð a( j, k) / a(k, k)
ïðïð
ïðïð
end
a( j, k +ð1: n) Źð a( j, k +ð1: n) -ð l( j, k) ×ð a(k, k +ð1: n)
ïðïð
x( j) Źð x( j) / a( j, j)
b( j) Źð b( j) -ð l( j, k) ×ð b(k)
ïðïð
ïðîðend
ïð
end
ïð
îðend
Wybór elementu głównego a faktoryzacja LU
Zrozumienie w jaki sposób wybór elementu głównego wpływa na postać faktoryzacji LU
wymaga wprowadzenia tzw. macierzy permutującej. Ogólnie, macierz permutująca
otrzymywana jest w wyniku przestawienia wierszy (lub - równoważnie  kolumn) macierzy
jednostkowej I. Na przykład:
1 0 0 0 0 1
éðÅ‚ð éðÅ‚ð
Ä™ð0 1 0Å›ð Þð P :=ð Ä™ð0 1 0Å›ð
I =ð
Ä™ðÅ›ð Ä™ðÅ›ð
Ä™ðÅ›ð Ä™ðÅ›ð
ëð0 0 1ûð ëð1 0 0ûð
Co się stanie jeśli dana macierz A zostanie pomnożona przez macierz permutującą? Oto
odpowiedz
:
0 0 1 1 2 3 7 8 9
éð Å‚ð éð Å‚ð éð Å‚ð
Ä™ð0 1 0Å›ð Ä™ð4 5 6Å›ð =ð Ä™ð4 5 6Å›ð przestawienie wierszy 1 i 3 w macierzy A
PA =ð -ð
Ä™ð Å›ð Ä™ð Å›ð Ä™ð Å›ð
Ä™ð Å›ð Å›ð Å›ð
ëð1 0 0ûð Ä™ð 8 9ûð Ä™ð 2 4ûð
ëð7 ëð1
1 2 3 0 0 1 3 2 1
éð Å‚ð éð Å‚ð éð Å‚ð
Ä™ð4 5 6Å›ð Ä™ð0 1 0Å›ð =ð Ä™ð6 5 4Å›ð przestawienie kolumn 1 i 3 w macierzy A
AP =ð -ð
Ä™ð Å›ð Ä™ð Å›ð Ä™ð Å›ð
Ä™ð Å›ð Å›ð Å›ð
ëð7 8 9ûð Ä™ð 0 0ûð Ä™ð 8 7ûð
ëð1 ëð9
Rozważmy teraz proces eliminacji z wyborem EG. Ponieważ wynikające z wyboru EG
przestawienie równań układu poprzedza eliminacje kolejnej niewiadomej, zatem
transformacja macierzy A(k-ð1) do macierzy A(k ) może być zapisana w postaci nastÄ™pujÄ…cej
A(0) ºð A A(k ) =ð L(k )P(k )A(k-ð1) , k =ð 1,2,..,n -ð1
,
W konsekwencji
U ºð A(n-ð1) =ð L(n-ð1)P(n-ð2)A(n-ð2) =ð L(n-ð1)P(n-ð1)L(n-ð2)P(n-ð2)...L(1)P(1)A
Fenomenalna własność macierzy permutujących P i faktorów L polega na tym, że ich
naprzemienny iloczyn da się zapisać następująco
(n-ð1)
éðL(n-ð1)L(n-ð2)...L(1) Å‚ð
éðûð
L(n-ð1)P(n-ð1)L(n-ð2)P(n-ð2)...L(1)P(1) =ð×ð
ëðP P(n-ð2)...P(1) Å‚ð
ëðûð
gdzie oznaczyliśmy
L(k ) =ð P(n-ð1)...P(k+ð1)L(k )[P(k+ð1) ]-ð1...[P(n-ð1) ]-ð1
-ð1
éðL(n-ð1)L(n-ð2)...L(1) Å‚ð
Co wiÄ™cej, macierz odwrotna L =ð istnieje i jest dolna trójkÄ…tna!
ëðûð
WprowadzajÄ…c zatem globalnÄ… macierz permutujÄ…cÄ… P =ð P(n-ð1)P(n-ð2)...P(1) otrzymujemy
następującą postać faktoryzacji LU
PA =ðLU
Faktoryzacja LU z wyborem EG może być wykorzystana do rozwiązania liniowego układu
równań z następujący sposób
ìð
Ly =ð Pb
Ax =ð b Þð PAx =ð Pb Þð LUx =ð Pb Þð
íð
Ux =ð y
îð
Jak widać, jedyna komplikacja polega na tym, że w programie komputerowym należy
zapamiętać historię przestawień kolejności równań dokonanych w procesie eliminacji. W
praktycznej implementacji załatwia się do za pomocą dodatkowego wektora zawierającego
odpowiednio poprzestawiane liczby naturalne od 1 do n = dim(A).
Inne użyteczne zastosowania faktoryzacji LU:
1. Obliczanie wyznacznika macierzy A
PA =ð LU Þð det(PA) =ð det(LU)
ßð
det(P) det(A) =ð det(L) det(U)
1(even permut.) =ð 1( product of
or -ð1(odd perm.) diagonal elem.)
ßð
det(A) =ð Ä…ðdet(U) =ð u11u22...unn
2. Obliczanie macierzy odwrotnej (na ogół należy go unikać!)
AA-ð1 =ð I Þð Axi =ð ei , ei Źð i -ð ta kolumna macierzy I , i =ð 1,2,..,n
Ostatecznie A-ð1 =ð [x1 | x2 |... | xn ]
Uwaga: Powyższe układy mogą być rozwiązane łącznym kosztem proporcjonalnym do n3.
Naiwne zastosowanie eliminacji Gaussa dałoby koszt całkowity proporcjonalny do n4.
UWAGI KOC
COWE
·ð IstniejÄ… inne użyteczne warianty faktoryzacji macierzy. W szczególnoÅ›ci każda macierz może być
przedstawiona w formie iloczynu macierzy ortogonalnej Q i górnej trójkątnej R (faktoryzacja QR).
Dla macierzy symetrycznych i dodatnio okreÅ›lonych (A =ð AT , xTAx >ð 0 dla każdego x Ä…ð 0)
istnieje faktoryzacja Choleskiego (A =ð LLT ). Metody te omawiane sÄ… m.in. w trakcie kursu
 Metody Numeryczne prowadzonym na WMEiL.
·ð Metoda eliminacji Gaussa i faktoryzacja LU należą do tzw.  metod dokÅ‚adnych , tj. algorytmów,
które dają dokładne (przynajmniej w teorii) rozwiązanie układu po skończonej, z góry określonej
liczbie operacji arytmetycznych. Dla układów o bardzo dużych rozmiarach koszt numeryczny i
wymagania co do pamięci komputera są niepraktyczne i/lub nierealistyczne. Wielkie układy
równań (o wymiarach rzędu n 104i więcej) rozwiązuje się metodami przybliżonymi, tj.
metodami w których rozwiązanie osiągane jest na drodze kolejnych iteracji. Niektóre z tych
metod omawiane są w kursie  Metody Numeryczne , a także w innych prowadzonych na Wydziale
MEiL kursach poświęconych zastosowaniom metod komputerowych w rozmaitych działach
techniki.