1
UKŁADY RÓWNAŃ LINIOWYCH
1. Wstęp
Rozwiązanie układu równań liniowych o postaci:
a x
11 1
+ a x
12
2
+ L + a x
1
= b
n
n
1
a x
21 1
+ a x
22
2
+ L + a x
2
= b
n
n
2
[1]
M
M
O
M
M
a x
1 1
+ a x
2
2
+ L + a x
= b
n
n
nn
n
n
polega na znalezieniu wartości zmiennych x ( i = ,
1 ,
2 ..., n) jeśli znane są wartości
i
współczynników a i b ( i = ,
1 ,
2 ..., n,
j = ,
1 ,
2 ..., n) .
ij
i
Układ [1] można zapisać za pomocą sumy:
n
∑ a x = b ( i = ,1 ,2..., n)
[2]
ij
j
i
j=1
lub w postaci macierzowej:
AX=B
[3]
gdzie:
a
a
L
a
b
x
11
12
n
1
1
1
a
a
L
a
b
x
A = 21
22
2 n B = 2 X = 2
[4]
M
M
O
M
M
M
a
a
L
a
b
x
n 1
n 2
nn
n
n
Jeśli macierz A jest nieosobliwa (det A ≠ 0) to układ [1] posiada tylko jedno rozwiązanie. W całym rozdziale zakładamy, że macierz A i wektor B są rzeczywiste.
Istnieje wiele metod znajdowania rozwiązania układu [1]. Metody te można podzielić na dwie grupy:
− metody dokładne (wzory Cramera, metoda Gaussa, metoda Gaussa-Jordana, metoda wyboru elementu podstawowego, metoda Choleskiego czyli rozkładu na macierze L i U),
− metody iteracyjne (metoda Jacobiego zwana też metodą iteracji prostej, metoda Seidela, metoda relaksacji).
2. Wzory Cramera
Jeśli macierz A jest nieosobliwa to układ [1] jest oznaczony i można wykazać, że rozwiązania dokładne ma postać:
det A
x
i
=
( i = ,
1 ,
2 ..., n)
[5]
i
det A
2
gdzie:
a
L
a
b
a
L
a
11
,
1 i−1
1
,
1 i+1
n
1
A
i = M
O
M
M
M
O
M
L
L
a
a
b
a
a
n 1
n, i−1
n
n i
, +1
nn
jest macierzą utworzoną z macierzy A, w której w miejsce i-tej kolumny wstawiono wektor B. Metoda oparta na wzorach Cramera wymaga obliczenia wartości n+1 wyznaczników. Dla dużych układów równań metoda ta jest mało efektywna ze względu na dużą złożoność obliczeniową oraz istotny wpływ błędów operacji arytmetycznych.
3. Metoda eliminacji Gaussa
Metoda eliminacji Gaussa polega na sprowadzeniu układu równań AX=B do układu o postaci A(n)X=B(n) gdzie A(n) jest macierzą trójkątną górną, a następnie rozwiązaniu tego trójkątnego układu równań.
Niech dany będzie układ równań:
)
1
(
( )
1
)
1
(
( )
1
)
1
(
a x
+ a x
+ a x
+ L + a x
= b
11
1
12
2
13
3
1 n
n
1
)
1
(
( )
1
)
1
(
( )
1
)
1
(
a x
+ a x
+ a x
+ L + a x
= b
21
1
22
2
23
3
2 n
n
2
( )
1
( )
1
)
1
(
( )
1
)
1
(
a
x
+ a x
+ a x
+ L + a x
= b
[6]
31
1
32
2
33
3
3 n
n
3
M
M
M
O
M
M
)
1
(
( )
1
)
1
(
( )
1
)
1
(
a x
+ a x
+ a x
+ L + a x
= b
1
n
1
n 2
2
n 3
3
nn
n
n
zakładając że ( )
1
a
≠ 0 pomnóżmy pierwsze równanie przez ( )1
)
1
(
a
a
i odejmijmy od i-tego
11
i 1
11
równania ( i = ,
2 ,
3 ..., n) , otrzymamy wtedy
1
( )
( )
1
1
( )
1
( )
( )
1
a x
+ a x
+ a x
+ L + a x
= b
11
1
12
2
13
3
1 n
n
1
( 2)
(2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
22
2
23
3
2 n
n
2
( 2)
(2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
[7]
32
2
33
3
3 n
n
3
M
M
O
M
M
( 2)
(2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
n 2
2
n 3
3
nn
n
n
)
1
(
( )
1
a a
( )
1
( )
1
a b
gdzie
( 2)
( )
1
i 1
1
a
= a
j
−
, (2)
)
1
(
i 1
1
b
= b −
( i, j = ,
2 ,
3 ..., n) .
ij
ij
)
1
(
a
i
i
( )
1
a
11
11
Następnie pomnóżmy drugie równanie przez (2)
( 2)
a
a
( (2)
a
≠ 0 ) i odejmijmy od i-
i 2
22
22
tego równania ( i = ,
3 ,
4 ..., n) , otrzymamy więc
( )
1
)
1
(
)
1
(
)
1
(
)
1
(
a
x
+ a x
+ a x
+ L + a x
= b
11
1
12
2
13
3
1 n
n
1
( 2)
(2)
( 2)
( 2)
a
x
+ a x
+ L + a x
= b
22
2
23
3
2 n
n
2
(3)
(3)
(3)
a
x
+ L + a x
= b
[8]
33
3
3 n
n
3
M
O
M
M
(3)
(3)
(3)
a
x
+ L + a x
= b
n 3
3
nn
n
n
3
( 2)
( 2)
a
a
( 2)
( 2)
a
b
gdzie
(3)
(2)
i 2
2
a
= a
j
−
, (3)
( 2)
i 2
2
b
= b −
( i, j = ,
3 ,
4 ..., n) .
ij
ij
( 2)
a
i
i
( 2)
a
22
22
Kontynuując takie postępowanie otrzymamy układ trójkątny
( )
1
)
1
(
)
1
(
)
1
(
)
1
(
a
x
+ a x
+ a x
+ L + a x
= b
11
1
12
2
13
3
1 n
n
1
( 2)
(2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
22
2
23
3
2 n
n
2
(3)
(3)
(3)
a
x
+ L + a x
= b
[9]
33
3
3 n
n
3
O
M
M
( n )
( n )
a
x
= b
nn
n
n
( k )
( k )
a
a
( k )
( k )
+
a
b
gdzie
( k + )
1
( k )
ik
kj
a
= a −
, ( k )
1
( k )
ik
k
b
= b −
, ( k )
a
kk
≠ 0
ij
ij
( k )
a
i
i
( k )
a
kk
kk
( k = ,
1 ,
2 ..., n −1 i, j = k + ,
1 k + ,
2 ..., n) .
Rozwiązanie tego układu równań jest proste. Rozwiązujemy go od „tyłu” tzn.
obliczamy kolejno x , x
, x
,..., x
n
n 1
−
wg wzorów:
n−2
1
n
( i )
( i )
b
− ∑ a x
( n )
i
ij
j
b
= +
j i 1
n
x =
,
x =
( i = n − ,
1 n − ,
2 ..., ,
2 )
1 .
[10]
n
( n )
i
( i )
a
a
nn
ii
4. Metoda wyboru elementu podstawowego
W trakcie przekształcania układu równań do układu z macierzą trójkątną górną musi być spełniony warunek
( k )
a
( k = ,
1 ,
2 ..., n) .
[11]
kk
≠ 0
Również przypadek gdy moduł współczynnika
( k )
a
ma bardzo małą wartość może
kk
prowadzić do znacznych błędów zaokrąglania spowodowanych operacją dzielenia przez bardzo małą wartość.
Aby uniknąć dzielenia przez zero lub przez małe liczby stosuje się metodę wyboru elementu podstawowego. Algorytm jest następujący:
1. W pierwszym kroku wyszukujemy element o maksymalnym module wśród wszystkich współczynników ( )
1
a
ij
( )
1
( )
1
max a
= a .
ij
rs
i, j= .
1 .. n
Niech nim będzie element leżący w r-tym wierszu i s-tej kolumnie. Zamieniamy r-ty wiersz z pierwszym a następnie s-tą kolumnę z pierwszą. Dalej dokonujemy redukcji macierzy wg wzoru [7].
2. W drugim kroku wyszukujemy elementu o maksymalnym module wśród wszystkich współczynników (2)
a
ij
( 2)
( 2)
max a
= a
.
ij
pq
i, j= .
2 .. n
Zamieniamy odpowiednio p-ty wiersz z drugim wierszem i q-tą kolumnę z drugą itd.
4
Jeśli w którymś kroku znaleziony element o maksymalnym module jest równy zeru to obliczenia przerywamy. Oznacza to bowiem, że det A = 0 , czyli układ jest albo sprzeczny albo nieoznaczony.
Należy również zapamiętywać wszystkie kolejne zamiany kolumn, ponieważ powodują one zmianę kolejności niewiadomych x w wektorze X.
i
5. Metoda Gaussa-Jordana
Metoda ta stanowi pewną modyfikację metody Gaussa. Otóż przekształca się układ równań [6] do układu w którym macierz współczynników jest macierzą jednostkową AX=B → EX=B(n)
wg następującego algorytmu. Dzielimy obustronnie pierwszy wiersz układu równań [6] przez współczynnik ( )
1
a
a następnie mnożymy przekształcony pierwszy wiersz przez współczynnik
11
( )
1
a
i odejmujemy od i-tego wiersza ( i = ,
2 ,
3 ..., n) . Otrzymamy
i 1
(2)
( 2)
(2)
( 2)
x
+ a x
+ a x
+ L + a x
= b
1
12
2
13
3
1 n
n
1
(2)
( 2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
22
2
23
3
2 n
n
2
(2)
( 2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
32
2
33
3
3 n
n
3
M
M
O
M
M
(2)
( 2)
(2)
( 2)
a
x
+ a x
+ L + a x
= b
n 2
2
n 3
3
nn
n
n
( )
1
a
( )
1
b
gdzie
( 2)
1
a
j
=
( j = ,
2 ,
3 ..., n)
(2)
1
b
=
( )
1
a
≠ 0
1 j
( )
1
a
1
( )
1
a
11
11
11
( 2)
( )
1
)
1
(
( 2)
a
= a − a a , (2)
)
1
(
( )
1
(2)
b
= b − a b ( i, j = ,
2 ,
3 ..., n) .
ij
ij
i 1
1 j
i
i
i 1
1
W otrzymanym układzie równań dzielimy drugi wiersz przez współczynnik ( 2)
a
a
22
następnie mnożymy przekształcony drugi wiersz przez współczynnik (2)
a
i odejmujemy od
i 2
i-tego wiersza ( i = ,
1 ,
3 ..., n) . Otrzymamy
(3)
(3)
(3)
x
+
a
x
+ L + a x
= b
1
13
3
1 n
n
1
(3)
(3)
(3)
x
+ a x
+ L + a x
= b
2
23
3
2 n
n
2
(3)
(3)
(3)
a
x
+ L + a x
= b
33
3
3 n
n
3
M
O
M
M
(3)
(3)
(3)
a
x
+ L + a x
= b
n 3
3
nn
n
n
( 2)
a
( 2)
b
gdzie
(3)
2
a
j
=
( j = ,
3 ,
4 ..., n)
(3)
2
b
=
(2)
a
≠ 0
2 j
( 2)
a
2
( 2)
a
22
22
22
(3)
(2)
( 2)
(3)
a
= a − a a , (3)
( 2)
( 2)
(3)
b
= b − a b ( i = ,
1 ,
3 ..., n,
j = ,
3 ,
4 ..., n ) .
ij
ij
i 2
2 j
i
i
i 2
2
Kontynuując obliczenia po „k” krokach otrzymamy
5
( k + )
1
( k + )
1
( k + )
1
x
+
a
x
L
a
x
b
k
k
+
+
n
n
=
1
,
1
1
+
1
+
1
1
( k + )
1
( k + )
1
( k + )
1
x
+
a
x
L
a
x
b
k
k
+
+
n
n
=
2
2,
1
+
1
+
2
2
O
M
M
=
M
[12]
( k + )
1
( k + )
1
( k + )
1
x
a
x
L
a
x
b
k
+
k k
k
+
+
kn
n
=
,
1
+
1
+
k
M
M
M
( k + )
1
( k + )
1
( k + )
1
a
x
L
a
x
b
n k
k
+
+
nn
n
=
,
1
+
1
+
n
( k )
a
( k )
+
b
gdzie
( k + )
1
kj
a
=
( j = k + ,
1 k + ,
2 ..., n)
( k
)
1
k
b
=
( k )
a
kk
≠ 0
kj
( k )
a
k
( k )
a
kk
kk
( k + )
1
( k )
( k )
( k + )
1
a
a
a
a
, ( k+ )
1
( k )
( k )
( k + )
1
b
b
a
b
( i = ,
1 ,
2 ..., n, i ≠ k j = k + ,
1 ..., n ) .
i
= i −
ij
= ij − ik kj
ik
k
Po „n-1” krokach otrzymamy (realizujemy obliczenia wg wzoru [12] dla k = , 1 ,
2 ..., n −1)
( n )
( n )
x
+
a
x
= b
1
1 n
n
1
( n )
( n )
x
+
a
x
= b
2
2 n
n
2
O
M
M
( n )
( n )
x −
+ a
x
−
= b
n 1
n
,
1 n
n
n 1
−
( n )
( n )
a
x
= b
nn
n
n
W ostatnim n-tym kroku wystarczy podzielić ostatnie równanie przez współczynnik ( n )
a
a następnie wyrugować zmienną x z równań 1,2 do n-1, czyli rozwiązaniem jest: nn
n
( n )
bn
x =
( n)
a
[13]
nn
≠ 0
n
( n )
ann
( n )
( n )
x
.
i = bi
− a x
i
in
n
= ,
1 ,
2 ..., n −1
W trakcie eliminacji Gaussa-Jordana współczynniki ( k )
a
( k = ,
1 ,
2 ..., n) . Aby
kk
≠ 0
warunek ten był spełniony należy również stosować metodę wyboru elementu podstawowego.
6. Metoda rozkładu na macierze L i U
Rozpatrzmy układ równań opisany równaniem [3] AX=B.
Twierdzenie 1. Jeśli wszystkie wyznaczniki podmacierzy Aii ( i = 1,2,..., n- 1) utworzonych z
„i” pierwszych kolumn i wierszy macierzy A licząc odpowiednio z lewej do prawej strony i z góry na dół, są różne od zera ( det A
) to istnieje jednoznaczny rozkład macierzy A na
ii ≠ 0
dwie macierze trójkątne L i U, odpowiednio dolną z jedynkami na przekątnej głównej i górną o postaci:
A=LU
[14]
1
0
L
0
u
u
L
u
11
12
1 n
l
1
L
0
u
L
u
21
0
22
2 n
gdzie L =
U =
.
[15]
M
M
O
M
M
M
O
M
L
l
l
1
L
0
0
u
n 1
n 2
nn
Algorytm rozwiązania układu równań AX=B jest więc następujący. Wstawiając równanie [14] do [3] otrzymamy LUX=B, które możemy rozbić na dwa równania LY=B i
6
UX=Y. Z pierwszego równania wyznaczamy wektor Y a z drugiego wektor X. Zaletą takiego rozbicia jest fakt, że obydwa wektory otrzymujemy natychmiast, bowiem rozwiązujemy dwa układy równań z macierzami trójkątnymi:
y 1 = b
1
i−1
[16]
y
,
2 ,
3 ...,
i = bi − ∑ l
y
i
ik
k
=
n
k =1
yn
xn =
unn
n
[17]
yi
∑
−
u x
ik
k
xi =
k = i+1
i = n − ,
1 n − ,
2 ..., 1
uii
Elementy macierzy L i U można wyznaczyć dwoma sposobami.
1.
Istnieje ścisły związek między rozkładem macierzy A na macierze L i U a metodą eliminacji Gaussa. Można wykazać, że elementy kolejnych kolumn macierzy L są równe współczynnikom przez które mnożone są w kolejnych krokach wiersze układu równań celem dokonania eliminacji niewiadomych w odpowiednich kolumnach
( k )
aik
l =
k = ,
1 ,
2 ..., n
i = k + ,
1 ..., n .
[18]
ik
( k )
akk
Natomiast macierz U jest równa macierzy trójkątnej uzyskanej w eliminacji Gaussa, czyli
( k )
u = a
k = ,
1 ,
2 ..., n
i = k, k + ,
1 ..., n .
[19]
ki
ki
2.
Algorytm Doolittle’a. Mnożąc macierz L i U a następnie porównując z macierzą A odpowiednio pierwszy wiersz i pierwszą kolumnę, drugi wiersz i drugą kolumnę itd otrzymamy
u
,
1 ,
2 ...,
1 j = a
j
1 j
=
n
a j 1
[20]
l
,
2 ,
3 ...,
j 1 =
j =
n
u 11
dla i = ,
2 ,
3 ..., n
i−1
uij = aij − ∑ l u
j
ik
kj
= i, i + ,
1 ..., n
k =1
i−1
[21]
a ji − ∑
l u
jk
ki
l ji =
k =1
j = i + ,
1 i + ,
2 ..., .
n
uii
Rozkład macierzy na macierze L i U może być wykorzystany do rozwiązania innych problemów matematycznych, takich jak obliczanie wartości wyznacznika dowolnej macierzy lub wyznaczanie macierzy odwrotnej.
1.
Obliczanie wartości wyznacznika.
7
Niech A=LU. Wiemy że det A = det L*det U. Zauważmy, że wartość wyznacznika macierzy trójkątnej jest równa iloczynowi wartości elementów na przekątnej głównej n
czyli det A = ∏ u .
[22]
ii
i=1
2.
Wyznaczanie macierzy odwrotnej.
Zauważmy, że AA-1=E. Jeśli jako niewiadome przyjmiemy kolejne kolumny macierzy A-1 a jako wektor wyrazów wolnych odpowiednią kolumnę macierzy jednostkowej E, to problem sprowadzi się do rozwiązania n układów równań o n niewiadomych każdy co można zapisać następująco: niech X(k) oznacza k-tą kolumnę macierzy A-1 a E(k) odpowiednio k-tą kolumnę macierzy jednostkowej. Należy rozwiązać n układów równań AX(k) = E(k) k = 1, 2, ..., n. Korzystając z rozkładu macierzy A na macierze L i U otrzymamy LUX(k) = E(k) . Ostatecznie musimy rozwiązać 2n układów równań z macierzami trójkątnymi LY(k) = E(k) oraz UX(k) = Y(k) k = 1, 2, ..., n.
7. Metoda Jacobiego (iteracji prostej)
Metoda ta zaliczana jest do metod iteracyjnych. Układ równań [1] przekształcamy następująco. Dzielimy i-ty wiersz przez współczynnik a
( a ≠ 0) i = ,
1 ,
2 ..., n i
ii
ii
pozostawiamy zmienną x po lewej stronie równania a pozostałe wyrazy przenosimy na i
prawą stronę równania. Otrzymamy nowy układ równań
x
L
1
=
α x
12
2
+ α x
13
3
+
+
α x
1
+ β
n
n
1
x
L
2
= α x
21 1
+
α x
23
3
+
+
α x
2
+ β
n
n
2
[23]
M
M
M
M
O
M
M
x
= α x
L
1 1
+ α x
2
2
+ α x
3
3
+
+ α
x
, −1
−1
+ β
n
n
n
n
n n
n
n
0
i = j
b
gdzie α =
i
− a
β =
=
ij
,
i, j
,
1 ,
2 ..., n
[24]
ij
i ≠
i
j
a
a
ii
ii
Układ równań można przedstawić w zapisie macierzowym
X = αX + β
[25]
gdzie elementy macierzy α i wektora β określają wzory [24].
Układ [25] rozwiązuje się metodą kolejnych przybliżeń. Za zerowe X(0) przybliżenie przyjmujemy wektor β (X(0) = β) i kolejno wyliczamy: pierwsze przybliżenie X(1) = αX(0) + β, drugie przybliżenie X(2) = αX(1) + β itd. Ogólnie (k+1)-sze przybliżenie wyznacza się wg wzoru
X(k+1) = αX(k) + β k = ,
1 ,
2 ...
[26]
Wzór [26] można zapisać w postaci skalarnej
x(0) = β
i
i
n
i =
n
k =
( k + )
1
( k
,
1 ,
2 ...,
,
0 ,
1 ,
2 ...
[27]
x
= ∑α x ) + β
i
ij
j
i
j= ,
1 j≠ i
8
Twierdzenie 2. Jeśli ciąg przybliżeń X(0) , X(1) , X(2) , ..., X(k) , ... ma granicę ( k
lim X ) = X
g
k →∞
to stanowi ona rozwiązanie układu [25] czyli układu [1].
Z twierdzenia 2 wynika, że metoda ta nie jest metodą dokładną, posiada tzw. błąd metody.
Warunkiem wystarczającym zbieżności procesu iteracyjnego jest, by dowolna z norm macierzy α była mniejsza od jedności co można zapisać
α < 1 lub α < 1 lub α
< 1.
[28]
I
II
III
Z [24] i definicji normy [30] wynika, że elementy leżące na przekątnej głównej macierzy A powinny być silnie dominujące nad pozostałymi elementami.
Definicja 1. Norma macierzy A (oznaczana jako A ) jest liczbą, stanowiącą w pewnym sensie miarę macierzy, spełniającą następujące warunki:
1)
A ≥ 0 , A = 0 ⇔ A = ,
0
2)
cA = c A , c ∈ R ,
[29]
3)
A + B ≤ A + B ,
4)
A B ≤ A B .
Najczęściej stosuje się trzy następujące normy macierzy A:
n
A
= max∑ a ,
ij
I
i
j=1
n
A
= max ∑ a ,
[30]
ij
II
j
i=1
n
n
A
= ∑∑ a 2 .
ij
III
i=1 j=1
Jeśli macierzą jest wektor kolumnowy X to odpowiednie normy wektora X mają postać X
= max x ,
i
I
i
n
X
= ∑ x ,
[31]
i
II
i 1
=
n
X
= ∑ x 2 .
i
III
i=1
8. Metoda Seidela
Metoda ta stanowi modyfikację metody iteracji prostej. Jest więc również metodą niedokładną. Polega na tym, że przy obliczaniu (k+1)-szego przybliżenia niewiadomej ( k + )
1
x
wykorzystuje się już obliczone (k+1)-sze przybliżenia zmiennych ( k+ ) 1
( k + )
1
x
, ..., x
oraz k-
i
1
i 1
−
te przybliżenia pozostałych zmiennych
( k )
( k )
x
, ..., x
i 1
+
. Czyli przekształcamy układ równań
n
AX=B do postaci X=αX+β gdzie α i β określone są wg wzorów [24]. Natomiast wzór [27]
ulega przekształceniu do postaci
9
x(0) = β
i
i
n
k = ,
0 ,
1 ,
2 ...
( k + )
1
( k
x
[32]
1
= ∑α x )
1
+ β
j
j
1
i = ,
2 ,
3 ..., n
j=2
i−1
n
( k + )
1
( k + )
1
( k
x
= ∑α x
+ ∑α x ) + β
i
ij
j
ij
j
i
j=1
j= i+1
Warunki wystarczające zbieżności procesu iteracyjnego w metodzie Seidela są identyczne jak w metodzie Jacobiego (wzór [28]). Zazwyczaj proces Seidela jest szybciej zbieżny niż proces iteracji prostej (Jacobiego). Często jest zbieżny podczas gdy proces Jacobiego jest rozbieżny.
9. Stop obliczeń iteracyjnych
Jako warunek zatrzymania obliczeń iteracyjnych najczęściej przyjmuje się by jedna ze zmodyfikowanych norm różnicy wektorów k-tego i (k+1)-szego przybliżenia była mniejsza od przyjętej dokładności ε lub gdy liczba wykonanych iteracji osiągnie maksymalną przyjętą wartość. Jako zmodyfikowane normy wektora proponuje się
*
( k + )
1
X
− ( k)
X
=
( k + )
1
max x
x
,
i
− ( k)
i
≤ ε
I
i
n
( k + )
1
( k )
∑ xi
− xi
*
( k + )
1
( k )
i 1
X
− X
= =
≤ ε ,
[33]
II
n
2
n
( k + )
1
( k )
∑ xi
− xi
*
( k + )
1
( k )
i 1
X
− X
=
=
≤ ε .
III
n
To że
( k + )
1
X
− ( k)
X
≤ ε nie zawsze oznacza, że X
X
gdzie X jest
d −
( k + )
1
≤ ε
d
wartością dokładną.
10. Metoda relaksacji
Przekształćmy układ równań [23] przenosząc wyrazy stojące po lewej stronie znaku równości na stronę prawą. Otrzymamy
− x
+ α x
+ L + α x
β
n
n
+
= 0
1
12
2
1
1
α x −
x
+ L + α x
β
n
n
+
= 0
21 1
2
2
2
M
M
O
M
M
M
α
[34]
x
α x
L
α x
β
s
+
s
+
+
sn
n
+
s
= 0
1 1
2
2
M
M
O
M
M
M
α x
α x
L
x
β
n
+
n
+
−
n
+
n
= 0
1 1
2
2
Używając zapisu macierzowego mamy
– X + αX + β = 0
[35]
gdzie elementy macierzy α {α } oraz wektora β { β } określają wzory [24].
ij
i
Wektor residuów definiujemy jako
R = − X + α X + β
[36]
10
Algorytm metody relaksacji
1. Krok pierwszy
Niech X(0) będzie zerowym przybliżeniem rozwiązania układu [35]. Wstawiając X(0) do równania [36] otrzymamy wektor residuów
R(0) = – X(0) + αX(0) + β
Szukamy największej bezwzględnej wartości residuum
max r (0) = r(0)
i = ,
1 ,
2 ..., n
i
s
i
i określamy wektor pierwszej poprawki
( )
1
δ X następująco
r(0) i
)
1
(
s
= s
δ x
,
1 ,
2 ..., .
i
=
i =
n
0
i ≠ s
Wektor pierwszego przybliżenia określamy następująco
)
1
(
(0)
)
1
(
X
= X +δ X .
Wyliczmy wektor residuów dla tego przybliżenia
( )
1
( )
1
( )
1
(0)
( )
1
(0)
)
1
(
(0)
)
1
(
)
1
(
R
= − X +α X + β = − X −δ X + X
α
+αδ X + β = R −δ X +αδ X .
Zauważmy że
( )
1
(0)
( )
1
( )
1
r
r
δ x
α δ x
s
= s − s + ss s = 0
( )
1
(0)
( )
1
)
1
(
(0)
(0)
r
= r −δ x +α δ x = r +α r i = , 1 ,
2 ..., n
i ≠ s
i
i
i
is
s
i
is s
czyli maksymalne bezwzględne residuum zostało wyzerowane a pozostałe residua zwiększyły swoją wartość o
(0)
α r .
is s
2. K-ty krok obliczeń
Określamy wektor k-tej poprawki
( k )
δ X następująco
( k−
r
)
1
i
( k )
p
= p
δ x
,
1 ,
2 ...,
i
=
i =
n
0
i ≠ p
gdzie max r( k− )
1
= r( k− )1 i = ,
2
,
1
..., n .
i
p
i
Wektor k-tego przybliżenia określamy następująco
( k )
( k − )
1
( k )
X
= X
+δ X
i wyliczmy wektor residuów dla tego przybliżenia
( k )
( k )
( k )
( k − )
1
( k )
( k − )
1
( k )
( k − )
1
( k )
( k )
R
= − X
+ X
α
+ β = − X
−δ X + X
α
+αδ X + β = R
−δ X +αδ X .
Zauważmy że
( k )
( k
)
1
( k )
( k )
r
r
δ x
α δ x
p
= −
p
−
p
+ pp p = 0
( k )
( k − )
1
( k )
( k )
( k − )
1
( k − )
1
r
r
δ x
α δ x
r
α r
i = ,
1 ,
2 ..., n
i ≠ p .
i
= i
− i + ip p = i
+ ip p
Tak więc metoda relaksacji polega na zerowaniu w kolejnych krokach obliczeniowych maksymalnej bezwzględnej wartości residuum. Cel ten osiągamy przez nadanie odpowiadającemu temu residuum składnikowi odpowiedniego wektora poprawek wartości maksymalnego residuum. Pozostałe składniki wektora poprawek mają wartość zero.
Obliczenia przerywamy, gdy wszystkie residua danej iteracji są mniejsze od przyjętej dokładności ε, czyli
max r( k ) ≤ ε
i = ,
1 ,
2 ..., n .
i
i
11. Program ćwiczenia
1. Napisać program rozwiązujący układ równań metodą Gaussa oraz rozwiązać zadane przez prowadzącego układy równań.
11
2. Napisać program rozwiązujący układ równań metodą Gaussa-Jordana oraz rozwiązać zadane przez prowadzącego układy równań.
3. Napisać program rozwiązujący układ równań metodą Choleskiego oraz rozwiązać zadane przez prowadzącego układy równań.
4. Napisać program rozwiązujący układ równań metodą iteracji prostej (Jacobiego) oraz rozwiązać zadane przez prowadzącego układy równań.
5. Napisać program rozwiązujący układ równań metodą Seidela oraz rozwiązać zadane przez prowadzącego układy równań.
6. Napisać program rozwiązujący układ równań metodą relaksacji oraz rozwiązać zadane przez prowadzącego układy równań.
7. Dla metod dokładnych sprawdzić czy macierz współczynników jest macierzą źle uwarunkowaną.
8. Dla metod iteracyjnych sprawdzić warunki zbieżności.
9. Dla metod iteracyjnych porównać dokładność rozwiązania wykonując tę samą liczbę iteracji i przyjmując ten sam wektor zerowego przybliżenia.
10. Napisać sprawozdanie zawierające:
a) tekst programu i wyniki przeprowadzonych obliczeń,
b) opis przeprowadzonych badań,
c) analizę uzyskanych wyników,
d) wnioski.
12. Zestaw pytań kontrolnych
1. Podaj warunki jakie musi spełniać układ równań liniowych oznaczony, nieoznaczony i sprzeczny.
2. Podaj wzory Cramera, jakie są wady tej metody.
3. Przedstaw algorytm metody eliminacji Gaussa.
4. Podaj rozwiązanie układu równań z macierzą trójkątną górną.
5. Na czym polega metoda wyboru elementu podstawowego i po co się ją stosuje?
6. Podaj algorytm metody Choleskiego oraz warunki rozkładu na macierze L i U.
7. Co to są macierze L i U? Jak wyznaczyć wyznacznik macierzy A stosując do tego celu rozkład na macierze L i U?
8. Jak w oparciu o rozkład macierzy A na macierze L i U wyznaczyć macierz A-1 ?
9. Przedstaw algorytm metody iteracji prostej (Jacobiego).
10. Co stanowi rozwiązanie dokładne w metodzie Jacobiego i jakie są warunki zbieżności tej metody?
11. Przedstaw algorytm metody Seidela.
12. Co stanowi rozwiązanie dokładne w metodzie Seidela i jakie są warunki zbieżności tej metody?
13. Co to jest norma macierzy i jakie znasz przykłady norm? Co to jest macierz żle uwarunkowana?
14. Podaj warunki stopu obliczeń iteracyjnych.