Rozwiązywanie układów równań liniowych
Rozważania zamieszczone poniżej dotyczyć będą jedynie takich przypadków
gdy ilość równań jest równa ilości niewiadomych. To oznacza, że jeśli istnieje
rozwiązanie, jest ono jednoznaczne.
Każdy układ równań liniowych można zapisać w postaci macierzowej
Ax = b
a11 a12 L a1n
ł
ęa a22 L a2nś
gdzie A jest macierzą kwadratową (n,n)
21
ę ś
A =
współczynników układu postaci:
ę M M O M ś
ęa an2 L annś
n1
b1
ł
ęb ś
gdzie b jest wektorem kolumnowym (n,1)
2
ę ś
b =
wyrazów wolnych układu
ę M ś
ęb ś
n
x1
ł
ęx ś
gdzie x jest rozwiązaniem układu,
2
ę ś
x =
wektorem kolumnowym (n,1)
ę M ś
ęx ś
n
Z kursu algebry wiadomo, że układ równań jest układem niezależnym, czyli ma
rozwiązanie, gdy wyznacznik macierzy A jest różny zera (detA ą 0). Rozwiązaniem
układu jest wektor
x = A-1b
gdzie A-1 jest macierzą odwrotną macierzy A czyli taką, że ich iloczyn jest
macierzą jednostkową A-1 A = A A-1 = I.
Rachunek macierzowy dostarcza zależności za pomocą których sprawdzić można
istnienie rozwiązania układu, czyli wyznaczyć wyznacznik macierzy oraz sposoby
wyznaczania rozwiązania czyli konstruowanie macierzy odwrotnej. Zależności te,
z punku widzenia metod numerycznych są nieprzydatne. Wymagają tak dużego
nakładu obliczeń że ich implementacja numeryczna jest bezcelowa. Dlatego
powstało wiele innych metod rozwiązywania układów równań liniowych o
znacznie mniejszym nakładzie obliczeń oraz łatwych do zaimplementowania w
postaci algorytmu numerycznego. Metody te dzielimy na metody dokładne, czyli
takie które przy nieskończenie wielkiej dokładności obliczeń dałyby wynik
dokładny, i metody iteracyjne, gdzie konstruowany jest ciąg przybliżeń
rozwiązania zbieżny do wyniku dokładnego. Metody dokładne, ze względu na
skończoną dokładność obliczeń, mogą dawać wyniki obarczone poważnym
błędem. Metody iteracyjne wymagają spełnienia przez macierz A wielu
dodatkowych warunków by proces był zbieżny i prowadził do rozwiązania.
Przykład.
Znana ze szkoły średniej metoda Cramera dla układu n równań wymaga
obliczenia n+1 wyznaczników
a11 a12 L a1n b1 a12 L a1n
a11 b1 L a1n
a21 a22 L a2n b2 a22 L a2n
a21 b2 L a2n
Wg = Wx =
. . .
Wx =
1
2
M M O M M M O M
M M O M
an1 an2 L ann bn an2 L ann
an1 bn L ann
Rozwiązanie układu jest równe
Wx Wx . . .
1 2
x1 = x2 =
Wg Wg
Obliczanie wyznacznika jest bardzo kosztowne. Przy rozwinięciu wg kolumny
lub wiersza wymaga n! mnożeń. Czyli rozwiązanie układu równań metodą
Cramera to (n+1)! mnożeń.
Dla pewnych typów macierzy np. macierzy diagonalnych i trójkątnych, obliczania
wyznacznika i wyznaczenie rozwiązania układu jest wyjątkowo łatwe.
a11 0 L 0 a11 0 L 0 a11 a12 L a1n
ł ł ł
ę ś ęa a22 L 0 ś ę
0 a22 L 0 0 a22 L a2nś
21
ę ś ę ś ę ś
ę M M O M ś ę M M O M ś ę M M O M ś
ę ęa an2 L annś ę
0 0 L annś 0 0 L annś
n1
trójkątna trójkątna
diagonalna
dolna górna
Dla wszystkich pokazanych typów macierzy wyznacznik jest iloczynem wyrazów
n
na głównej przekątnej
det A =
a
ii
i=1
natomiast rozwiązanie układu A x = b jest następujące:
bk
xk = k = 1,2,..., n
- macierz diagonalna
akk
k-1
ć
1
bk -
xk =
a xi k = 1,2,3,..., n
ki
- macierz trójkątna dolna
akk
Ł i=1 ł
n
ć
1
bk -
- macierz trójkątna górna
xk =
a xi k = n, n -1, n - 2,...,1
ki
akk
Ł i=k+1 ł
Norma macierzy definicja w przestrzeni rzeczywistej
Niech
y : Rm,n [0,+Ą] będzie przekształceniem spełniającym warunki
"A Rm,n y (A) = 0 A = 0,
"A Rm,n "u R y (u * A) = u y (A),
"A, B Rm,n y (A + B) Ły (A) +y (B)
Każde takie przekształcenie y nazywamy normą w Rm,n i oznaczamy
y (A) = A
Norma jest miarą wielkości macierzy i dlatego jest uznawana za
A - B
miarę odległości macierzy A i B
Norma nazywamy submultiplikatywną jeśli spełnia warunek
A* B Ł A B "A,B Rn,n
Wektory są szczególnymi przypadkami macierzy. Przykładem norm wektorów
są normy Schura zdefiniowane następująco
1
n
p
ć
x = xi p
dla 1 Ł p <
p
Ł i=1 ł
przykłady
n
x = xi
p = 1
1
i=1
n
2
x =
x
i
p = 2
2
i=1
p Ą
x = max xi
Ą
1ŁiŁn
Normy p-te macierzy indukowane przez normy p-te wektorów definiowane są
następująco
A * x
p
A = sup
p
xą0 x
p
Z definicji indukowanej normy p-tej macierzy wynikają następujące związki:
- spełniają warunek zgodności z normą wektora
"A Rm,n "x Rn A* x Ł A x
p p p
- normy te są submultiplikatywne
"A Rm,l "B Rl,n A* B Ł A B
p p p
Przykłady norm macierzy
- normy indukowane
m
A = max aij
(max. suma modułów w kolumnie)
1
j=1,2,..,n
i=1
n
A = max aij
(max. suma modułów w wierszu)
Ą
i=1,2,..,m
j=1
A = największa wartość własna macierzy (ATA )
2
- norma nie indukowana (Frobeniusza)
m m
2
A =
a
ij
E
i=1 j=1
Dla dowolnej normy macierzy kwadratowej zgodnej z normą wektora
zachodzi
li Ł A
gdzie li jest dowolną wartością własną macierzy A
r(A) = max li
Liczba nazywana jest promieniem spektralnym macierzy.
i=1,2,..,n
Stąd wynika, że
r(A) Ł A p = 1,2,Ą, E
p
Numeryczne rozwiązywanie układów równań liniowych, ze względu na skończoną
dokładność obliczeń, zawsze obarczone jest pewnym błędem. Jeśli macierz A i
wektor b zastąpimy zaburzonymi wartościami A+dA oraz b+db uzyskamy
rozwiązanie x+dx.
Rozpatrzmy błąd rozwiązania w najprostszym przypadku gdy :
db `" 0, dA = 0
Wówczas bład rozwiązania wynika z zależności
A(x + dx) = b + db dx = A-1db
Korzystając z własności norm można zapisać
dx Ł A-1 p db
p p
A-1 p
Liczba jest miarą wpływu zaburzenia na wielkość błędu rozwiązania
Najczęściej interesuje nas nie wartość bezwzględnego błędu rozwiązania
lecz oszacowanie jego wartości względnej czyli wpływu zaburzenia
db dx
p p
na błąd względny
b x
p p
Wykorzystując z własności norm mamy
A-1 p b
dx db
p
p p
Ł
x x b
p p p
A-1 p b
p
k =
Liczba nazywana wzmocnieniem zaburzenia ma taką
x
p
wadę że trudno jest ją wyznaczyć bez znajomości rozwiązania
A-1 p Ax
p
Ale rozpisując dostaniemy k = Ł A-1 p A = K
p
p
x
p
Liczba Kp (niezależna od rozwiązania ) nazywana jest wskaznikiem
uwarunkowania układu równań
METODY DOKAADNE
Metoda eliminacji Gaussa
Metoda ta polega na sprowadzeniu równania Ax = b z pełną macierzą kwadratową
~
~ ~
A
do równoważnej postaci Ax = b , gdzie jest górną macierzą trójkątną natomiast
~ przekształconą kolumną wyrazów wolnych. Proces ten, nazywany często
b
eliminacją zmiennych, przebiega następująco:
- pierwszy wiersz macierzy A oraz pierwszy element wektora b mnożymy przez
-a21/a11 i dodajemy do wiersza drugiego,
- pierwszy wiersz macierzy A oraz pierwszy element wektora b mnożymy przez
-a31/a11 i dodajemy do wiersza trzeciego itd.. W rezultacie dostaniemy
a11 a12 L a1n b1
ł ł
ę ś ę ś
a21
0 a22 - a12 a21 L a2n - a1n a21ś
ę ęb - b1 a11ś
2
a11 a11 ś
ę ę ś
M M O M M
ę ś ę ś
ę ębn - b1 an1ś
0 an2 - a12 a21 L ann - a1n an1 ś
ę ę
a11 a11 ś a11ś
~ ~z macierzą
Ax = b
Postępując analogicznie dla wierszy 2,3,...,n dostaniemy układ
trójkątną którego rozwiązanie otrzymamy stosując uprzednio podane wzory
1 1
Metoda eliminacji Gausa wymaga M = n3 + n2 - n mnożeń i dzieleń oraz
3 3
1 1 5
D = n3 + n2 - n
dodawań.
3 2 6
Rozwiązując układ równań metodą Cramera (wyznaczników) konieczne jest
wykonanie (n+1)! mnożeń. To oznacza że stosując metodę Cramera dla układu
20 równań maszyna wykona 2432902008176640000 mnożeń. Przyjmując moc
obliczeniową 109 mnożeń na sekundę, konieczny czas wykonania obliczeń to
2432902008 sekund czyli ponad 77 lat. Dla metody eliminacji Gaussa potrzeba
w tym wypadku wykonać 5910 operacji co zajmie około 610-6 sekundy
Metoda eliminacji Gaussa nie jest metodą niezawodną. W przypadku gdy
element a11 macierzy A jest zerowy dostaniemy dzielenie przez zero czyli
krytyczny błąd obliczeń. By tego uniknąć stosowane jest takie przesuwanie
wierszy lub kolumn by element przez który dzielimy był zawsze równy od zera
Metody te nazywane są metodami wyboru elementu podstawowego
Metoda z pełnym wyborem elementu podstawowego (pivoting j.ang) polega
na tym że przed każdym kolejnym kroku redukcji kolejnej niewiadomej
przeglądamy macierz układu pod kątem odszukania największego co do modułu
jej elementu. Staje się on tzw. elementem podstawowym, czyli kolejną redukcję
prowadzi się dla niewiadomej skojarzonej z numerem kolumny elementu
podstawowego. Przed redukcją następuje cykl przestawienia kolumny i wiersza
elementu podstawowego do właściwej pozycji.
1 -2 2 x1 1 0 4 1 x1 3 4 0 1 x2 3
ł ł ł ł ł ł ł ł ł
ę3 0 -1ś ęx ś ę2ś ę3 0 -1ś ęx ś ę2ś ę ę2ś
= = 0 3 -1ś ę x1 ś =
22
ę ś ę ś ę ś ę ś ę ś ę ś ę ś ę ś ę ś
ę0 4 1 ś ęx3ś ę3ś ę1 -2 2 ś ęx3ś ę1ś ę-2 1 2 ś ęx3ś ę1ś
redukcja x2
4 0 1 x2 3 4 0 1 x2 3
ł ł ł ł ł ł
ę0 3 -1ś ę ś ę ś ę0 3 -1ś ę ś ę ś
x1 = 2 x1 = 2
ę ś ę ś ę ś ę ś ę ś ę ś
ę0 1 2.5ś ęx3ś ę2.5ś ę0 1 2.5ś ęx3ś ę2.5ś
Uwagi:
- z przestawieniem wiersza wiąże się przestawienie także wiersza prawych stron.
- z przestawieniem kolumny wiąże się przestawienie wiersza niewiadomych
(numeru niewiadomej)
Metoda z częściowym wyborem elementu podstawowego polega na tym że
przed każdym kolejnym kroku redukcji kolejnej niewiadomej przeglądamy nie
całą macierz układu pod kątem odszukania największego co do modułu jej
elementu, ale tylko kolumnę odpowiadającą aktualnie redukowanej
niewiadomej. Niepotrzebne jest wtedy przestawianie kolumn i reorganizacja
wektora niewiadomych.
1 -2 2 x1 1 3 0 -1 x1 2
ł ł ł ł ł ł
ę3 0 -1ś ęx ś ę2ś ę1 -2 2 ś ęx ś ę1ś
= =
22
ę ś ę ś ę ś ę ś ę ś ę ś
ę0 4 1 ś ęx3ś ę3ś ę0 4 1 ś ęx3ś ę3ś
redukcja x1
3 0 -1 x1 2
ł ł ł
ę0 -2 2 ś ęx ś ę ś
11
=
2
33
ę ś ę ś ę ś
ę0 4 1 ś ęx3ś ę3ś
Uwagi:
- z przestawieniem wiersza wiąże się przestawienie także wiersza prawych
stron.
Metoda zawsze prowadzi do rozwiązania (gdy ono istnieje i jest jednoznaczne)
nie zostanie przerwana na dzieleniu przez zero!!
Przypadki układów równań dla których nie ma potrzeby stosowania
wyboru (częściowego) elementu podstawowego.
Macierz diagonalnie (słabo) dominująca: moduły na diagonali są niemniejsze
od sumy pozostałych elementów w tym samym wierszu
3 0 -1 x1 2
ł ł ł
n
ę0 -2 2 ś ęx ś ę1ś
aii ł aik i = 1,2...,n przykad : =
2
ę ś ę ś ę ś
k =1,k ąi
ę0 4 -4ś ęx3ś ę3ś
Macierz diagonalnie (silnie) dominująca: moduły na diagonali większe od
sumy pozostałych elementów w tym samym wierszu
Macierz diagonalnie (silnie) dominująca kolumnowo: wtedy gdy AT jest
(silnie) diagonalnie dominująca.
TW. Jeśli macierz jest nieosobliwa i diagonalnie kolumnowo dominująca, to
przy eliminacji Gauss a nie ma potrzeby przestawiania wierszy (czyli
stosowania wyboru elementu podstawowego).
Jeżeli macierz A jest symetryczna i dodatnio określona nie trzeba stosować
wyboru elementu podstawowego w algorytmie Gauss a.
Metoda eliminacji Jordana
a11 a12 L a1n
ł
óó
1 0 L a1n
ł
ęa a22 L a2nś
ę0 1 L a2nś
óó
21
ę ś
ę ś
ę M M O M ś
ęM M O M ś
ęa an2 L annś
ę0 0 L annś
óó
n1
krok 2 krok 3
krok 1
ó ó
1 a12 L a1n
1 0 L 0
ł
ł
ę0 a22 L a2nś
ę0 1 L 0ś
ó ó
ę ś
ę ś
ęM M O M ś
ęM M O Mś
ę0 an2 L annś
ę0 0 L 1ś
ó ó
1 1
M = n3 + n2
Metoda eliminacji Jordana wymaga mnożeń i dzieleń oraz
2 2
1 1
D = n3 -
dodawań, czyli około 1,5 razy więcej aniżeli eliminacja Gaussa
2 2
Rozkład na macierze trójkątne (metoda Doolittla a)
Jeśli macierz kwadratową A zapisać można jako iloczyn dwóch macierzy:
trójkątnej dolnej L z jedynkami na głównej przekątnej, oraz trójkątnej górnej U
a11 a12 L a1n 1 0 L 0 u11 u12 L u1n
ł ł ł
ęa a22 L a2nś ęl
1 L 0ś . ę 0 u22 L u2nś
21 21
ę ś ę ś ę ś
=
ę M M O M ś ę M M O Mś ę M M O M ś
ęa an2 L annś ęl ln2 L 1ś ę
0 0 L unnś
n1 n1
wówczas dla każdego i = 1, 2,..., n
i-1
uij = aij -
l ukj
ik
j = i, i+1,..., n
k=1
i-1
aji -
l uki
jk
k=1
j = i+1, i+2,..., n
l =
ji
uii
Rozwiązanie układu Ax = b dostaniemy rozwiązując układy Ly = b oraz Ux = y
(ilość operacji arytmetycznych taka sama jak w metodzie eliminacji Gaussa).
Jeśli chodzi pozostałe własności (np. wybór elementu podstawowego i sytuacje
w których nie jest on potrzebny sa takie same jak dla alg. Gauss a)
Rozkład na macierze trójkątne (inne metody rozkładu)
W szczególnych wypadkach (gdy macierz A układu jest symetryczna i/lub dodatnio
określona ) można zastosować inne rozkłady macierzy A na czynniki, co pozwala
na zmniejszenie liczby niezbędnych operacji arytmetycznych (czyli przyspieszenie
obliczeń).
Rozkład A=LDLT , jest jednoznaczny pod warunkiem że A jest symetryczna.
d11 00 1 l21 ln1
a11 a12 L a1n 1 0 L 0 ł
ł
ł ł
ę0 1 ln2 ś
ęa a22 L a2nś ęl
0 d22 0
1 L 0ś ęś ś
21 21
ęś ę
ę ś ę ś
=
ęś ęś
ę M M O M ś ę M M O Mś
ś
ęa an2 L annś ęl ln2 L 1ś ę 0 0 dnn ś ę0 0
1
n1 n1
.
Algorytm rozkładu przebiega następująco:
j-1
ć
= aij -ik
c ljk / d
lij j
k =1
Łł
d1 = a11 cij = d lij j = 1,2,...,i -1 wykonywać naprzemiennie dla i = 2,3,...,n
j
i-1
di = aii -
c lik
ik
k =1
Liczba działań 2 razy mniejsza iż przy eliminacji Gauss a, ponadto wystarczy
pamiętać połowę macierzy A (ze względu na symetrię).
Rozkład na macierze trójkątne (inne metody rozkładu)
Rozkład A=LLT , jest jednoznaczny pod warunkiem że A jest symetryczna i
dodatnio określona. Zwany rozkładem Cholesky ego lub Banachiewicza.
a11 a12 L a1n
ł l11 00 l11 l21 ln1
ł ł
ęa a22 L a2nś
ęl l22 ś ę
0 0 l22 ln2 ś
21
21
ę ś
ęś ęś
=
ę M M O M ś ęś ęś
ęa an2 L annś ęl ln2 lnn ś ę 0 0 lnn ś
n1 n1
Algorytm rozkładu przebiega następująco:
i-1
2
lii = aii -
l dla i = 1,2,...,n
ik
k =1
.
i-1
ć
lji =
l lik / lii dla j = i +1,i + 2,...,n
ji jk
a -
k=1
Łł
Liczba działań około 2 razy mniejsza iż przy eliminacji Gauss a, ponadto wystarczy
pamiętać połowę macierzy A (ze względu na symetrię). Dodatnia określoność A
gwarantuje że elementy diagonali L będą niezerowe i nie będzie kłopotów z
dzieleniem przez zero.
Metody dokładne dla szczególnego typu macierzy c.d.
Macierzą rzadką nazywamy taką macierz A układu równań, która zawiera
mniej niż 20% elementów niezerowych. Jeśli struktura (położenie tych elementów
niezerowych) jest odpowiednia, to można te informacje wykorzystać do pominięcia
zbędnych operacji (mnożenia przez zero) w algorytmach eliminacji. Ponadto,
odpowiedni sposób kodowania takiej macierzy w pamięci operacyjnej komputera
pozwala na znakomitą oszczędność. Wiele metod dyskretnych rozwiązywania
zagadnień transportu (układy równań różniczkowych cząstkowych rzędu 2)
prowadzi do sformułowania układu równań algebraicznych liniowych jako
końcowego etapu rozwiązania. W wielu stosowanych metodach uzyskujemy
macierze rzadkie o szczególnych własnościach. Np. zastosowanie metody
różnicowej skutkuje otrzymaniem układu o macierzy A mającej niezerowe
elementy na określonych przekątnych macierzy A:
- Dla zagadnienia 1 wymiarowego i gwiazdy 3 punktowej macierz A jest 3
przekątniowa
- Dla zagadnienia 2 wymiarowego i gwiazdy 5 punktowej macierz A jest 5
przekątniowa
- Dla zagadnienia 3 wymiarowego i gwiazdy 7 punktowej macierz A jest 7
przekątniowa
Metody dokładne dla szczególnego typu macierzy
Macierz 3 przekątniowa (algorytm THOMASA a)
b1 c1 0 0 0
ł
ęa b2 c2 00 ś
2
ęś
Jeżeli macierz A (po lewej stronie) jest diagonalnie
ę 0 a3 b3 c3 0 ś
dominująca (wierszowo) to metoda eliminacji
ęś
ęś
Gauss a jest niezawodna. Jeśli zastosuje się
ęś
00 an-1 bn-1 cn-1
algorytm rozkładu A na iloczyn LU to procedura
ęś
0 0 0 an bn
obliczeniowa przebiega następująco:
l = b1 y1 = d1 / l
ui = ci / l l = bi+1 - uiai+1 yi+1 = (di+1 - ai+1yi ) /l dla i = 1,2,...,n -1
xn = yn
xi = yi - uixi+1 dla i = n -1,n - 2,...,1
Liczba operacji arytmetycznych do wykonania jest rzędu 6n . Zaletą metody jest
również to, że przystosowaniu podwójnej precyzji obliczeń można rozwiązać ta
metoda układy z precyzją 10% nawet gdy wskaznik uwarunkowania jest nie większy
od 10 do potęgi 10.
Podobne algorytmy można stosować dla macierzy 5 przekatniowej.
Metoda najszybszego spadku (steepest descent)
Przy założeniu że macierz A układu równań jest symetryczna i dodatnio
określona, rozwiązanie układu Ax=b można poszukiwać metodą minimalizacji
funkcji
1
J : Rn R , J (y) = yTAy - yTb
2
dla której minimum=0 jest osiągane dla y które jest rozwiązaniem układu Ax=b
Metoda ta jest metodą iteracyjną. Startujemy z próbnego rozwiązania x(0) i
generujemy ciąg rozwiązań x(k) takich że w kolejnych krokach J(x(k+1)
Kierunek przejścia z rozwiązania x(k) do kolejnego rozwiązania jest określony przez
wektor p(k) ,zaś odległość na którą przemieszczamy się w tym kierunku w celu
poszukiwania lepszego przybliżenia jest regulowana przez pewien współczynnik a(k)
x(k+1) = x(k ) +a(k )p(k )
Wartość a(k) jest dobierana drogą minimalizacji J zgodnie z kryterium:
J (x(k+1)) = min J (x(k ) +ap(k ))
a
Dla funkcji kwadratowej J procedura minimalizacji jest oczywista i daje wartość
T
unikalną wartość a(k)
p(k ) r(k )
( )
a(k ) = gdzie r(k ) = b - Ax(k )
T
p(k ) Ap(k )
( )
Metoda najszybszego spadku (steepest descent) c.d.
W metodzie najszybszego spadku wektor kierunku poszukiwania p(k) = DJ = r(k)
bowiem łatwo pokazać, że gradient funkcji J ma postać
T
ł
śJ śJ śJ
ŃJ (x) = ŃJ (y) = Ax - b = r
ęśx , śx2 ,..., śxn ś
1
Metoda taka wydaje się być bardzo oczywista, ale posiada pewne wady. Obliczenie
a(k) wymaga operacji mnożenia Ap(k) , która jest najbardziej pracochłonna (o ile
macierz nie jest rzadka) , i również policzenie r(k) wymaga takiego samego nakładu
obliczeń. Drugiego iloczynu możemy uniknąć, obliczając r(k+1) z formuły
rekurencyjnej
r(k+1) = r(k ) -a(k )Ap(k )
x2
Interpretacja geometryczna metody dla 2
niewiadomych. Startując z pkt. 1 poszukujemy w
kierunku linii czerwonej minimum, osiągamy pkt.2
itd. Jeśli funkcja J reprezentowana czarnymi liniami
stałej jej wartości, ma głęboką i długą dolinę, to
2
proces dojścia do minimum może być długi. Ale,
3
gdyby linie stałej wartości J tworzyły okręgi
współśrodkowe, to moglibyśmy dojść do minimum
1
w 1 kroku!!! (gradient jest prostopadły do linii !)
x1
Metoda najszybszego spadku (steepest descent) c.d.
Spowolnienie zbieżności metody spowodowane złym uwarunkowaniem układu,
powodującym długą dolinę funkcji J, może być poprawione drogą przekształcenia
układu do nowej jego równoważnej formy przez wymnożenie układu przez macierz
odwrotną do macierzy M zwanej preconditioner matrix. Najprostsza formą tej
macierzy jest metoda macierzy Jacobieg o, M-1=D-1 , gdzie D jest diagonalą
macierzy A. W tym wypadku, macierz Jacobieg o wpłynie dodatnie na zbieżność
tylko gdy elementy D znacznie się różnią co do rzędu wielkości.
Algorytm najszybszego spadku przy zastosowaniu macierzy M-1
r Ź r - Ax
Na początku rozwiązanie próbne jest w wektorze
p Ź M-1r
x, wektor prawych stron b jest załadowany do
k Ź 0
wektora r.
Pętla po k do osiągnięcia zbieżnosci
q Ź Ap
a = pTr / pTq
x Ź x + ap
2
r Ź r -aq
p Ź M-1r
k Ź k +1
powrót
Metoda Gradientów Sprzężonych (CG conjugate gradients)
Metoda najszybszego spadku posiada ograniczoną
pamięć , bowiem wyznacza nowe rozwiązanie
r Ź r - Ax
x(k+1) na podstawie tylko poprzedniego x(k) . Metoda
s Ź M-1r
CG różni się zasadniczo tylko tym, że posiada
v Ź rTs
pamięć i wyznacza x(k+1) na podstawie wielu
p Ź s
poprzednich rozwiązań ( a w zasadzie poprzednich
k Ź 0
kierunków poszukiwań p).
Pętla po k do osiągnięcia zbieżnosci
q Ź Ap
Algorytm CG przy zastosowaniu macierzy M-1
m Ź pTq
a Ź v / m
Na początku rozwiązanie próbne jest w wektorze
x Ź x + ap
r Ź r -aq x, wektor prawych stron b jest załadowany do
wektora r.
s Ź M-1r
v+ Ź rTs
b Ź v+ / v
Koszt obliczeniowy CG jest niewiele większy od
2
p Ź s + bp
metody najszybszego spadku
v Ź v+
k Ź k +1
powrót
Metoda CG i najszybszego spadku szybkość zbieżności
Porównanie szybkości zbieżności metody najszybszego spadku i CG dokonano dla
zagadnienia różniczkowego
ś2T ś2T
+ + f = 0 z w.b. T = g na brzegu
śx2 śy2
co po zastosowaniu różnic skończonych prowadziło do układu równań (stały
parametr siatki h)
-Ti-1, j - Ti, j-1 + 4Ti, j - Ti, j+1 - Ti+1, j = h2 fi, j i, j = 1,...,m n = m2
Testowano obie metody z/bez użycia macierzy preconditioning M. Macierz M
była typu SSOR:
w 11
ć
M = D - ED-1 ć D - ET
2 - w w w
Ł ł Ł ł
gdzie D diagonala A, -E dolna trójkątna część macierzy A, parametr w musi
być z zakresu (0,2)., w obliczeniach użyto w=1.9
Uwaga: nie ma potrzeby obliczania M-1, w to miejsce stosuje się jedną iterację
SSOR dla Ax=r startując z x=0.
2
Naszybszego CG z
Naszybszego
spadku z CG macierzą
spadku
macierzą M M
liczba iteracji 4010 58 118 28
METODY ITEARACYJNE
Tw. Ciąg określony wzorem x(i+1) = Mx(i) + w, gdzie M jest macierzą kwadratową
a w wektorem, jest zbieżny do jedynego punktu granicznego, dla dowolnego
wektora x(0), wtedy i tylko wtedy, gdy r(M) < 1 (promień spektralny macierzy
jest mniejszy od 1)
~
~ ~+
x
W punkcie będącym rozwiązaniem układu Ax = b musi zachodzić x = Mx w
Przyjmując, że wektor w zapisać można jako w = Nb, gdzie N jest pewną
macierzą kwadratową dostaniemy
A-1b = MA-1b + Nb czyli A-1 = MA-1 + N
Mnożąc obustronnie przez macierz A otrzymamy
M = I NA gdzie I jest macierzą jednostkową
Powyższa zależność pozwala na konstruowanie całej rodziny metod iteracyjnych
postaci
x(i+1) = (I NA)x(i) + Nb
umożliwiających wyznaczenie rozwiązania układu jeśli tylko r(I - NA) < 1
Metoda Jacobiego
Przyjmijmy, że macierz A układu równań Ax = b można rozłożyć następująco
A = L + D + U
0 a12 L a1n
a11 0 L 0
ł
0 0 L 0 ł
ł
ę0 0 L a2nś
ę ś
ęa
0 a22 L 0
0 L 0ś
21
ę ś
ę ś
ę ś D =
U =
L =
ęM M O M ś
ę M M O M ś
ę M M O Mś
ę0 0 L 0 ś
ę
ęa an2 L 0ś
0 an2 L annś
n1
Przyjmując N = D-1 mamy M = - D-1(L + U)
czyli
x(k+1) = D-1(- L - U)xk + D-1b
Zastosowanie metody Jacobiego wymaga takiego przestawienia wierszy lub
kolumn by współczynniki na przekątnej głównej były różne od zera, bowiem
tylko wtedy można wyznaczyć macierz odwrotną D-1. Nie gwarantuje to jednak
zbieżności metody czyli spełnienia warunku r(D-1(L + U)) < 1
Metoda Gaussa-Seidla
Zakładając jak w metodzie Jacobiego, że A = L + D + U oraz przyjmując
N = (D + L)-1 otrzymamy M = (D + L)-1U co daje wzór iteracyjny postaci
x(k +1) = D-1(- Lx(k +1) - Ux(k ))+ D-1b
Należy zwrócić uwagę, że po prawej stronie równania występuje człon Lx(k+1)
który może budzić wątpliwości co do słuszności wzoru. Okazuje się jednak, że
przy obliczaniu pierwszej współrzędnej wektora x(k+1) wykorzystywane są jedynie
współrzędne wektora x(k). Licząc kolejne współrzędne wykorzystujemy tylko
uprzednio wyznaczone współrzędne wektora x(k+1) oraz elementy wektora x(k).
Podobnie jak w metodzie Jacobiego konieczne jest przestawienia wierszy lub
kolumn by współczynniki na przekątnej głównej były różne od zera. Nie
gwarantuje to jednak zbieżności metody. Jeśli jednak potrafimy uzasadnić że
metoda Jacobiego jest zbieżna dla danej macierzy A, to zbieżna jest również
metoda Gaussa-Seidla i to znacznie szybciej.
Metody relaksacyjne
Rozwiązując układ równań metodą iteracyjną dążymy do minimalizacji
wektora residuum
r = b - A~
x
~
gdzie jest kolejnym przybliżeniem rozwiązania
x
Metoda Gaussa-Seidla należy do rodziny metod relaksacyjnych bowiem
przy wyznaczaniu k-tego przybliżenia rozwiązania uwzględniamy
przybliżenie k oraz k-1. Analizując szczegółowo tą metodę dochodzimy do
wniosku, że elementy wektora kolejnych przybliżeń rozwiązania dane są
zależnością
(k
rii )
xi(k ) = xi(k -1) +
aii
przy czym poszczególne elementy składnika wektora residuum
wynoszą
i-1 n
(k ) -1) -1)
rii ) = bi -
a x(jk - a x(jk - aiixi(k
ij ij
j=1 j=i+1
Rozważania powyższe pozwalają to stworzyć całą rodzinę metod relaksacyjnych
przy założeniu że kolejne poprawki przybliżeń
(k
rii )
aii
brane będą z mnożnikiem innym niż jeden, czyli
(k
rii )
xi(k ) = xi(k -1) +w
aii
Jeśli 0 < w <1 metodę nazywamy podrelaksacyjną, gdy w > 1 jest to metoda
nadrelaksacyjna.
W zapisie macierzowym wyniki kolejnych iteracji podane są wzorem
-1 -1
x(k ) = (D -wL) [(1-w)D +wU]x(k-1) +w(D -wL) b
Obliczanie wyznacznika
Obliczanie wyznacznika przy stosowaniu reguły rozwijania względem kolumny
lub wiersza jest bardzo kosztowne. Dlatego stosowane są inne metody np.
wykonanie rozkładu na macierze trójkątne LU oraz wykorzystanie zależności
detA = detLU = detLdetU = detU
Wyznaczenie wyznacznika macierzy U nie stanowi problemu gdyż jest to iloczyn
współczynników na głównej przekątnej
Odwracanie macierzy
Macierz odwrotną A-1 można obliczyć poprzez n-krotne rozwiązanie układu
równań
Ax(i) = e(i) i =1,2, ..., n
gdzie e(i) = [0, ..., 1, ...,0]T (jedynka na i-tej pozycji, wersor przestrzeni Rn )
Wyznaczone rozwiązania x(i) są kolejnymi kolumnami szukanej macierzy
odwrotnej A-1.
Numeryczne rozwiązywanie układów równań nieliniowych
Układ równań nieliniowych w ogólnej postaci zapisać można następująco
f1(x1, x2,...xn)
ć
f2(x1, x2,...xn)
F(x) =
...........................
fn(x1, x2,...xn)
Ł ł
Rozwiązanie układu to znalezienie takiego wektora x = [x1, x2,..., xn]T dla
którego zachodzi F(x) = 0.
Metoda Newtona
Metoda jest uogólnieniem znanej metody stycznych dla funkcji n zmiennych
śf1 śf1 śf1
-1
ł
x(i+1) = x(i) -[F'(x(i))] F(x(i))
ęśx1 śx2 L
śxn ś
ęśf śf2
śf2 ś
2
ę ś
L
[F'(x)]-1 F'(x) =
gdzie jest odwrotnością macierzy Jacobiego
ę śxn ś
śx1 śx2
ę ś
M M O M
ęśfn śfn
śfn ś
ęśx śx2 L
śxn ś
1
Przykład: Przy pomocy metody Newtona znalezć rozwiązanie układu równań
x2 + y2 - 5 = 0
x + y -1 = 0
Rozpisując zgodnie z wymaganiami metody Newtona dostaniemy
2 2
ł
x1 + x2 - 5
2x1 2x2
ł
F(x) =
ę ś F'(x) = J =
ę ś
x1 + x2 -1
1 1
1
ł
x(0) =
Przyjmijmy punkt startowy
ę0ś
wówczas
1
0ł
2 0
ł
- 4
ł ę ś
-1
2
F'(x0 ) =
F(x0 ) =
[F '(x0 )] =
ę1 1ś
ę ś ę ś
1
0
ę 1
- ś
2
1
0ł- 4ł 3 ł
1
ę ś
ł
2
x(1) = - =
ę ś
ę0ś ę ś ę
1
ę 1
- ś 0 - 2ś
2
2,00005
ł 2
2,2 2,01176 ł
ł ł
x(4) =
x(5) =
x(2) = x(3) =
ę ś
ę ś
ę ś ę ś
-1,00005
-1
-1,2 -1,01176
Wyszukiwarka
Podobne podstrony:
Wyklad 2 3 MACIERZE WYZNACZNIK UKLADY ROWNAN
uklady rownan liniowych
MN MiBM zaoczne wyklad 1 uklady rownan
Układy równań zadania
Macierze i układy równań przykłady
uklady rownan
C 02 Uklady równan
uklady rownan
4 uklady rownan liniowych
układy równań sprawozdanie7
t5 uklady rownan liniowych
BOiE układy równań liniowych
Uklady rownan 2
wykład 11 układy równań liniowych
4 Układy równań
więcej podobnych podstron