Mariusz PYRZ
SIMR (PW), Instytut Pojazdów
Metody numeryczne w mechanice
Układy algebraicznych równań liniowych
3.
Mariusz PYRZ
SIMR (PW), Instytut Pojazdów
Metody numeryczne w mechanice
Układy algebraicznych równań liniowych
3.
Wst
ę
p
Rozwiązywanie duŜych układów równań liniowych (oraz nieliniowych) jest
nieuniknione w analizie numerycznej zagadnień opisanych za pomocą równań
róŜniczkowych cząstkowych (np. za pomocą metody elementów skończonych,
metodą róŜnic skończonych, ...).
Rozpatrywane układy mogą liczyć wiele tysięcy równań - niezbędne jest zatem
dysponowanie efektywnym programem obliczeniowym do ich rozwiązania.
WaŜna jest ponadto znajomość głównych zasad z których korzystają takie
2
WaŜna jest ponadto znajomość głównych zasad z których korzystają takie
procedury.
M.Pyrz Metody numeryczne w mechanice – Układy równa
ń
liniowych 10.2011
Rachunek macierzowy – główne pojęcia
Rząd macierzy – liczba wektorów kolumnowych (albo wierszowych) liniowo niezaleŜnych.
Macierz:
♦
kwadratowa A=(a
ij
) i=1,…n, j=1,…n
♦
diagonalna (a
ij
=0 dla i≠j, a
ii
=1)
♦
symetryczna (a
ij
=a
ji
)
♦
antysymetryczna (a
ij
=-a
ji
)
♦
hermita (a
ij
=a
*
ji
– elementy połoŜone symetrycznie względem przekątnej są zespolone
sprzęŜone)
♦
trójkątna dolna L (elementy zerowe nad przekątną)
Rysunek 1
3
♦
trójkątna dolna L (elementy zerowe nad przekątną)
♦
trójkątna górna U (elementy zerowe pod przekątną)
♦
sprzęŜona A
*
z A : a
*
ij
=â
ij
gdzie â
ij
jest liczbą zespoloną sprzęŜoną z a
ij
.
♦
transponowana
Wyznacznik macierzy
(kwadratowej) rzędu n:
det
det
det
det
det
det(
)
det
det(
)
det
det
1
det(
)
det
n
gdy
istnieje
α
α
−
−
=
=
=
=
=
=
T
*
1
1
I
1
A
A
A
A
A
A
AB
A
B
A
A
A
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Rysunek 1
Rachunek macierzowy – główne pojęcia
Macierz odwrotna A
–1
macierzy kwadratowej A istnieje wtedy i tylko wtedy gdy det(A) ≠ 0
(A jest wtedy macierzą regularną).
JeŜeli det(A) = 0 to macierz A jest osobliwa.
Macierze podobne mają taki sam wyznacznik.
A’ jest
podobna
do A jeŜeli A’=S
-1
AS , gdyŜ det(S
-1
AS)=det(S
-1
) det(A) det (S).
Równanie Ax=b ma jednoznaczne rozwiązanie jeŜeli det(A) ≠ 0 (x=A
-1
b).
Macierz diagonalna
4
Macierz diagonalna
a
ij
=0 dla i≠j, wyznacznik
Macierze kwadratowe trójkątne
dolna a
ij
=0 dla i<j, górna a
ij
=0 dla i>j
Wyznacznik macierzy kwadratowej trójkątnej
Iloczyn dwóch macierzy trójkątnych (dolnych)
Niech L
1
, L
2
będą macierzami trójkątnymi dolnymi rzędu n.
Ich iloczyn L
3
=L
1
L
2
jest macierzą trójkątną dolną.
det A
=
=
∏
a
ii
i
n
1
det
,
A
≠ ⇔
≠ ∀
≤ ≤
0
0
1
a
i
i
n
ii
det A
=
=
∏
a
ii
i
n
1
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Rysunek 2
Układ m algebraicznych równań liniowych
• Układ m równań z n niewiadomymi: x
1
, x
2
, ..., x
n
11 1
1
1
21 1
2
2
...
...
n
n
n
n
a x
a x
b
a x
a x
b
+ +
=
+ +
=
5
JeŜeli m > n to układ równań nadokreślonych
JeŜeli m < n t o układ równań niedookreślonych
W dalszych rozwaŜaniach rozpatrujemy przypadek m = n
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
21 1
2
2
1 1
...
...
n
n
m
mn
n
m
a x
a x
b
+ +
=
Układ n algebraicznych równań liniowych
• W postaci macierzowej
11
12
1
1
1
21
22
2
2
2
...
...
...
...
...
...
...
...
n
n
a
a
a
x
b
a
a
a
x
b
=
Ax = b
6
Układ n liniowych równań algebraicznych o n niewiadomych
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
1
2
...
...
...
...
...
...
...
n
n
nn
n
n
a
a
a
x
b
=
Metody rozwiązywania
Problem:
Rozwiązać Ax=b,
macierz A=(a
ij
) (1 ≤ i ≤ n, 1 ≤ j ≤ n) jest znana,
wektor b=(b
i
) (1 ≤ i ≤ n) jest podany,
wektor x=(x
i
) (1 ≤ i ≤ n) jest poszukiwany (nieznany).
Zakładamy, Ŝe A ( o elementach R
n
lub C
n
) jest odwracalna tzn. det(A)≠0.
7
Zakładamy, Ŝe A ( o elementach R
n
lub C
n
) jest odwracalna tzn. det(A)≠0.
Metody bezpośrednie – generują rozwiązanie x w skończonej liczbie operacji
Metoda Gaussa, Jordana, Choleskiego, …
Metody iteracyjne - tworzą ciąg wektorów {x
1
,x
2
, …} dąŜący do dokładnego
rozwiązania x (uzyskuje się rozwiązanie przybliŜone, ale metody te są skuteczne
w przypadku macierzy duŜego rozmiaru, macierzy rzadkich tj. mających wiele
elementów zerowych, …)
Metoda Jacobi’ego, Gauss-Seidel’a, metody gradientowe, ...
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metody bezpośrednie
Układ o macierzy diagonalnej
A=(a
ij
), a
ij
=0 gdy i ≠ j
Rozwiązanie :
liczba potrzebnych operacji N
PRZEKATNA
= n dzieleń.
1
i
i
ii
b
x
dla
i
n
a
=
≤ ≤
8
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przykład 1
Układ o macierzy trójkątnej dolnej
(a
ij
=0 dla i<j)
Rozwiązanie Ax=b
jest równoznaczne rozwiązaniu
Rozwiazywanie "schodząc"
Metody bezpośrednie
a x
b
i
n
a x
a x
b
i
n
ij
j
i
j
i
ij
j
ii
i
i
j
i
=
≤ ≤
+
=
≤ ≤
=
=
−
∑
∑
1
1
1
1
1
i
=
−
L
M
O
P
=
−
∑
1
1
Rozwiazywanie "schodząc"
liczba niezbędnych operacji: patrz przypadek macierzy trójkątnej górnej
9
x
a
b
a x
i
n
i
ii
i
ij
j
j
=
−
L
NM
O
QP
=
=
∑
1
1 2
1
, ,...,
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Rozwiązanie układu o macierzy trójkątnej
10
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przyklad 2
Metody bezpośrednie
Układ o macierzy trójkątnej górnej
(a
ij
=0 dla i>j)
Rozwiązywanie "wchodząc"
a x
b
i
n
a x
a x
b
i
n
ij
j
i
j i
n
ij
j
ii
i
i
j i
n
=
≤ ≤
+
=
≤ ≤
=
= +
∑
∑
1
1
1
x
a
b
a x
i
n n
i
ii
i
ij
j
j n
n
=
−
L
NM
O
QP
=
−
= +
∑
1
1
1
1
,
,...,
11
liczba niezbędnych
działań:
N
TRIAN
= n+n(n-1)/2+ n(n-1)/2=n
2
działań elementarnych ( + - * / )
a
ii
j n
N
Q
= +
1
(
1)
1 2
(
1)
2
(
1)
1 2
(
1)
2
n
dzielen
n n
n
dodawan
n n
n
mnozen
−
+ + + − =
−
+ + + − =
⋯
⋯
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Aby rozwiązać układ równań liniowych Ax=b (A jest macierzą rzędu n, kwadratową,
regularną), dokonujemy jego przekształcenia (w n-1 etapach) do równowaŜnego
układu równań Ux=c (gdzie U jest macierzą trójkątną górną) – co moŜna nazwać
„tringularyzacją”.
Otrzymany
układ
równań
rozwiązuje
się
procedurami
przystosowanymi do macierzy trójkątnych.
12
Metoda Gaussa odpowiada dekompozycji (rozkładowi)
A=LU (wtedy c=L
-1
b)
A
x =
b
=
x
c
U
A
=
U
L
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Uwaga 1:
JeŜeli musimy rozwiązać wiele układów równań postaci Ax
1
=b
1
, Ax
2
=b
2
… to
korzystne jest najpierw przeprowadzenie rozkładu A=LU a następnie kolejne
rozwiazywania układów równań: Lz
1
=b
1
z
1
, Ux
1
=z
1
x
1
, Lz
2
=b
2
z
2
,
Ux
2
=z
2
x
2
, … w których macierze L, U są stałe dla wszystkich układów równań.
13
Uwaga 2:
JeŜeli A jest symetryczna to poszukuje się raczej rozkładu typu LDL
T
gdzie D jest
macierzą diagonalną (uzyskując zmniejszenie liczby działań i niewiadomych) –
oferuje nam to na przykład metoda Cholesky’ego (wymagająca n
3
/6 działań, w
przeciwieństwie do n
3
/3 w rozkładzie LU).
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Ma na celu zastąpienie układu wyjściowego przez układ o macierzy trójkątnej,
posiadający takie samo rozwiązanie.
Numer etapu będzie oznaczony jako górny indeks w nawiasach, np. układ do
rozwiązania będzie zapisany A
(1)
=A, b
(1)
=b
A
(1)
x=b
(1)
.
Dzięki zamianie (permutacji) linii (lub kolumn) macierzy A
(1)
moŜemy załoŜyć,
Ŝe a
11
(1)
≠
0.
14
Dla linii 2
≤
i ≤ n :
mnoŜymy pierwsza linie
macierzy A
(1)
przez
r
i1
= -a
i1
(1)
/a
11
(1)
i dodajemy wynik
do
i-tego
równania
a
a
r a
j
n
i
n
a
a
j
n
a
i
n
b
b
r b
i
n
b
b
ij
ij
i1
j
j
j
i
i
i
i
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
,
2
1
1
1
1
2
1
1
1
2
2
1
1 1
1
1
2
1
1
2
2
1
0
2
2
=
+
≤ ≤
≤ ≤
=
≤ ≤
=
≤ ≤
=
+
≤ ≤
=
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Równanie A
(1)
x=b
(1)
jest równowaŜne A
(2)
x=b
(2)
gdzie A
(2)
=M
(1)
A
(1)
, b
(2)
=M
(1)
b
(1)
.
A
b
(2)
(2)
=
L
N
M
M
M
M
O
Q
P
P
P
P
=
L
N
M
M
M
M
O
Q
P
P
P
P
a
a
a
a
a
a
a
b
b
b
n
n
n
nn
n
11
1
12
1
1
1
22
2
2
2
2
2
2
1
1
2
2
2
0
0
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
15
Uwaga :
det(A
(2)
)=a
(1)
11
Minor(a
11
(1)
) ≠ 0
co upowaŜnia do przejścia
do następnego etapu.
Minor – wyznacznik macierzy kwadratowej powstałej z danej
macierzy przez skreślenie pewnej liczby wierszy i kolumn
M
(1)
=
L
N
M
M
M
M
M
M
O
Q
P
P
P
P
P
P
1
0
0
0
1
0
0
0
1
0
0
0
1
21
31
1
⋯
⋯
⋯
⋯
⋯ ⋯ ⋯ ⋯
⋯
r
r
r
n
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Uwaga: operacje na wierszach macierzy A i wektora b są identyczne. Aby
uprościć zapis moŜna zbudować macierz zwana rozszerzoną, w której
wektor b staje się (n+1) kolumną, a składowa b
i
będzie oznaczona a
i,n+1
.
A
b
( )
.
(.)
(.)
(.)
(.)
(.)
(.)
(.)
(.)
(.)
(.)
(.)
,
(.)
(.)
(.)
(.)
,
(.)
=
L
M
M
M
O
P
P
P
=
L
M
M
M
O
P
P
P
+
+
a
a
a
b
a
a
a
b
a
a
a
a
a
a
a
a
n
n
n
n
n
n
11
12
1
1
21
22
2
2
11
12
1
1
1
21
22
2
2
1
⋯
⋯
⋯
⋯
16
Macierz rozszerzona A
(2)
jest zdefiniowana jako
A
b
( )
.
(.)
(.)
(.)
(.)
,
(.)
(.)
(.)
,
(.)
=
N
M
M
M
Q
P
P
P
=
N
M
M
MM
Q
P
P
PP
+
+
a
a
a
b
a
a
a
b
a
a
a
a
a
a
a
a
n
n
n
nn
n
n
n
n
n
nn
n n
21
22
2
2
1
2
21
22
2
2
1
1
2
1
⋯
⋮
⋮
⋯
⋮
⋮
⋯
⋯
⋮
⋮
⋯
⋮
⋮
⋯
a
a
r a
i
n
j
n
a
a
j
n
a
i
n
ij
ij
i1
j
j
j
i
( )
( )
( )
( )
( )
( )
,
;
2
1
1
1
1
2
1
1
1
2
2
1
1
1
1
0
2
=
+
≤ ≤
≤ ≤ +
=
≤ ≤ +
=
≤ ≤
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
17
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przykład 3
Metoda eliminacji Gaussa
Na k-tym etapie
A
(k)
x=b
(k)
A
b
(k)
(k)
=
L
N
M
M
M
M
M
M
M
M
M
O
Q
P
P
P
P
P
P
P
P
P
=
L
N
M
M
M
M
M
M
M
M
M
O
Q
P
P
P
P
P
P
P
P
P
−
−
−
−
−
−
−
a
a
a
a
a
a
a
a
a
b
b
b
b
n
n
k
k
k
k
n
k
kk
k
kn
k
k
k
k
k
k
k
k
11
1
12
1
1
1
22
2
2
2
1
1
1
1
1
1
1
2
2
1
1
0
0
( )
( )
( )
( )
( )
,
(
)
,
(
)
( )
( )
( )
( )
( )
( )
(
)
( )
( )
⋯
⋯
⋯
⋯
⋯
⋯
⋮
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
18
Przejście od etapu (k) (A
(k)
x=b
(k)
) na etap (k+1) (A
(k+1)
x=b
(k+1)
)
dokonywane jest przez pomnoŜenie k-tej linii A
(k)
przez r
ik
= - a
ik
(k)
/a
kk
(k)
i dodanie
wyniku do i-tej linii dla numerów k+1 ≤ i ≤ n.
Macierz rozszerzona A
(k+1)
przybiera postać
NM
QP
NM
QP
a
a
b
nk
k
nn
k
n
k
0
0
0
( )
( )
( )
⋯
(
1)
( )
( )
(
1)
( )
1
,
1
1
1
1, 1
;
k
k
k
ij
ij
ik
kj
k
k
ij
ij
a
a
r a
k
i
n
k
j
n
a
a
j
n
i
k
pozostale elementy zerowe
+
+
=
+
+ ≤ ≤
+ ≤ ≤ +
=
≤ ≤ +
≤ ≤
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Układ A
(k+1)
x=b
(k+1)
moŜna otrzymać przez M
(k)
A
(k)
x=M
(k)
b
(k)
gdzie
M
(k)
=
L
M
M
M
M
M
M
M
M
O
P
P
P
P
P
P
P
P
+
1
0
0
0
0
1
0
0
0
1
0
1
0
0
0
1
0
1
⋯
⋯ ⋯
⋯
⋯ ⋯
⋯
⋯
⋯ ⋯
⋯
⋯ ⋯
⋯ ⋯ ⋯
⋯
⋯
⋯
⋯
⋯ ⋯
⋯
⋯
r
k
k
,
19
Po ostatnim etapie (n-1) uzyskujemy uklad o macierzy trójkątnej górnej
A
(n)
x=b
(n)
gdzie
A
(n)
=M
(n-1)
M
(n-1)
…M
(1)
A
(1)
b
(n)
=M
(n-1)
M
(n-1)
…M
(1)
b
(1)
Podstawiając U=A
(n)
oraz L=[M
(n-1)
M
(n-1)
…M
(1)
]
-1
moŜna zapisać macierz A w postaci A=LU.
N
M
M
M
Q
P
P
P
1
0
0
0
0
1
⋯
⋯ ⋯
⋯
⋯
⋯
⋯
r
nk
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa - wyznacznik
Liczby a
kk
(k)
nazywamy
elementami podstawowymi.
We wszystkich przekształceniach załoŜyliśmy, Ŝe elementy podstawowe są
róŜne od 0. Algorytm powinien zatem zawierać procedurę sprawdzania
niezerowej wartości elementu podstawowego (uzupełnioną o ewentualną
zamianę linii).
20
Obliczanie wyznacznika
Metoda Gaussa pozwala na łatwe obliczenie wyznacznika macierzy A.
PoniewaŜ det(M)=1, mamy :
( )
1
det
det
( 1)
n
p
k
kk
k
a
=
=
= −
∏
(n)
A
A
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Wybór elementów podstawowych
Ze względu na stabilność algorytmu, mała wartość elementu podstawowego jest
odradzana i wykorzystać moŜna wtedy następujące strategie:
Dobór elementu podstawowego „w kolumnie”
Na k-tym etapie jako element podstawowy a
ik
(k)
(k ≤ i ≤ n) wybieramy
Wybieramy zatem jako linię elementu podstawowego A
(k)
tę, która zawiera element o
maksymalnej wartości bezwzględnej w kolumnie k i dokonujemy permutacji linii numer k
oraz linii w której znaleziono element maksymalny.
a
a
ik
k
k i n
ik
k
( )
( )
sup
=
≤ ≤
21
oraz linii w której znaleziono element maksymalny.
Dobór elementu podstawowego „w podmacierzy”
Na k-tym etapie jako element podstawowy obieramy wartość
Jako element podstawowy wybierzemy element o maksymalnej wartości bezwzględnej w
podmacierzy «nieprzekształconej» (n-k+1) x (n-k+1) i dokonujemy odpowiedniej permutacji
k-tej kolumny i tej, w której znaleziono maksimum oraz linii k i odpowiedniej linii w której
znaleziono maksimum.
Uwaga:
jeśli dokonamy permutacji kolumn w macierzy A, to musimy dokonać
odpowiadającej jej permutacji elementów wektora rozwiązania x !
a
a
ij
k
k i n
k
j n
ij
k
( )
( )
sup
=
≤ ≤
≤ ≤
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Gaussa – liczba działań
Przekształcenie do macierzy trójkątnej
dla stałego k:
(n-k) dzieleń aby obliczyć r
ik
(n-k)
2
dodawań i (n-k)
2
mnoŜeń aby obliczyć a
ij
(k+1)
(n-k) dodawań i (n-k) mnoŜeń aby obliczyć b
i
(k+1)
a
a
r a
k
i
n
k
j
n
b
b
r b
k
i
n
ij
k
ij
k
ik
kj
k
i
k
i
k
ik
k
k
(
)
( )
( )
(
)
( )
( )
,
+
+
=
+
+ ≤ ≤
+ ≤ ≤
=
+
+ ≤ ≤
1
1
1
1
1
22
(n-k) dodawań i (n-k) mnoŜeń aby obliczyć b
i
(k+1)
rozwiązanie « wchodząc »
n dzieleń, n(n-1)/2 dodawań, n(n-1)/2 mnoŜeń
(
)
(
)
(
)
(
)
(
)(
)
(
)
(
)
n
k
n n
n
k
n
k
q
q
n n
n
n n
n n
k
n
q
n
q
n
k
n
−
=
−
−
+ −
=
+
=
−
− + − =
−
=
−
=
−
=
−
=
−
∑
∑
∑
∑
1
2
1 2
1
6
1
3
1
3
1
1
2
2
1
1
1
1
1
1
2
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Metoda Gaussa wymaga ogółem
2
3
2
2
3
2
(
1)
(
1)
2
2
(
1)
(
1)
2
3
5
2
3
6
(
1)
(
1)
2
3
5
n n
n n
n
dzielen
n n
n n
n
n
n
dodawan
n n
n n
n
n
n
mnozen
−
+
+
=
−
−
+
−
+
=
−
−
+
−
+
=
23
działań
(tj rzędu 2n
3
/3 dla odpowiednio duŜego n).
Przykład: dla n=10 : N
GAUSS
=805
Metoda współczynników Cramera: (n+1)! mnoŜeń
N
n
n
n
GAUSS
=
+
−
4
9
7
6
3
2
(
1)
(
1)
2
3
5
2
3
6
n n
n n
n
n
n
mnozen
−
−
+
−
+
=
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana
Polega na przekształceniu początkowego układu równań Ax=b do układu
o macierzy diagonalnej SAx=Sb gdzie SA jest macierzą jednostkową
A
x =
b
=
x
d
1
1
1
1
1
1
0
0
24
Procedura transformacji do macierzy jednostkowej jest dokonywana w
n etapach, z których kaŜdy składa się z fazy normalizacji (dzielenie
wszystkich elementów linii przez wybrany element) po której następuje
faza redukcji (kaŜda znormalizowana wybrana linia jest mnoŜona
przez element z odpowiedniej kolumny innej linii a wynik odejmowany
jest od tej “innej” linii).
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana - algorytm
Etap 1:
normalizacja
redukcja
Etap 2:
normalizacja
redukcja
a
a
a
j
n
j
j
1
2
1
1
11
1
1
1
( )
( )
( )
=
≤ ≤ +
a
a
a a
j
n
i
n
ij
ij
i
j
( )
( )
( )
( )
,
2
1
1
1
1
2
1
1
2
=
−
≤ ≤ +
≤ ≤
a
a
a
j
n
j
j
2
3
2
2
22
2
2
1
( )
( )
( )
=
≤ ≤ +
a
a
a
a
j
n
i
n i
( )
( )
( )
( )
,
,
3
2
2
3
2
1
1
2
=
−
≤ ≤ +
≤ ≤
≠
25
redukcja
Etap k:
normalizacja
redukcja
Końcowa postać macierzy A to macierz jednostkowa I.
W kolumnie (n+1) macierzy poszerzonej otrzymujemy rozwiązanie x
i
=
a
i,n+1
(n)
gdzie 1 ≤ i ≤ n.
a
a
a
a
j
n
i
n i
ij
ij
i
j
( )
( )
( )
( )
,
,
3
2
2
2
2
3
2
1
1
2
=
−
≤ ≤ +
≤ ≤
≠
a
a
a
k
j
n
kj
k
kj
k
kk
k
(
)
( )
( )
+
=
≤ ≤ +
1
1
a
a
a
a
k
j
n
i
n i
k
ij
k
ij
k
ik
k
kj
k
(
)
( )
( )
(
)
,
,
+
+
=
−
≤ ≤ +
≤ ≤
≠
1
1
1
1
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana
26
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przykład 4
Metoda eliminacji Gaussa-Jordana – liczba działań
2
2
(
1)
2
(
1)
2
(
1)
n n
dzielen
n n
dodawan
n n
mnozen
+
−
−
27
Ogółem N
JORDAN
=n
3
+n
2
/2-n/2.
Dla odpowiednio duŜego n, liczba działań jest rzędu n
3
.
Przykład: dla n=10 : N
JORDAN
=1045.
(
1)
2
n n
mnozen
−
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Obliczanie macierzy odwrotnej
Aby otrzymać wektory kolumnowe x
i
macierzy odwrotnej A
-1
wystarczy
rozwiązać n układów równań Ax
i
=e
i
gdzie e
i
jest i
tym
wektorem bazowym
e
1
=[1 0 0 … 0], e
2
=[0 1 0 … 0], … e
n
=[0 0 0 … 1].
28
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przykład 5
Obliczanie macierzy odwrotnej
Macierz A
-1
moŜna uzyskać po jednej transformacji do macierzy diagonalnej
stosując macierz rozszerzoną, składającą się z n linii i 2n kolumn
[ A | e
1
| e
2
| … |e
n
].
Kolumny macierzy odwrotnej A
-1
zajmują miejsca kolumn wstępnie
zajmowanych przez odpowiednie wektory bazowe e
i
.
29
Algorytm:
a
a
a
k
j
n
a
a
a
a
k
j
n
i
n i
n
kj
k
kj
k
kk
k
ij
k
ij
k
ik
k
kj
k
(
)
( )
( )
(
)
( )
( )
(
)
,
+
+
+
=
≤ ≤
=
−
≤ ≤
≤ ≤
≠
1
1
1
2
2
1
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Uklady o macierzach symetrycznych
Macierz symetryczna moŜe być zawsze przedstawiona w postaci LU lub
LDR gdzie
L – macierz trójkątna z jedynkami na przekątnej,
U – macierz trójkątna górna,
U=DR gdzie D jest macierzą diagonalną
R jest macierzą trójkątną górną z jedynkami na przekątnej.
30
R jest macierzą trójkątną górną z jedynkami na przekątnej.
JeŜeli macierz regularna, symetryczna A posiada rozkład LDR
to wtedy A= LDL
T
.
Dowód:
PoniewaŜ A jest symetryczna A=LDR=(LDR)
T
=R
T
DL
T
.
Jednoznaczność rozkładu daje L=R
T
, R=L
T
.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Choleskiego
Rozkład regularny Choleskiego macierzy A przybiera postać A=BB
T
gdzie
B jest macierzą trójkątną dolną regularna.
JeŜeli współczynniki na przekątnej macierzy B są dodatnie, mówimy o
dodatnim rozkładzie Choleskiego.
Macierz A posiada regularny rozkład Choleskiego wtedy i tylko wtedy gdy
jest ona dodatnio określona (wówczas posiada ona rozkład dodatni
31
jest ona dodatnio określona (wówczas posiada ona rozkład dodatni
jednoznaczny).
Definicja :
Macierz A jest dodatnio określona jeśli iloczyn
skalarny (Ax,x) > 0, dla kaŜdego x z R
(C)
-{0}.
JeŜeli macierz A jest symetryczna i dodatnio określona, to posiada
jednoznaczny rozkład A=LDL
T
(elementy na przekątnej są dodatnie).
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Choleskiego - algorytm
Rozwiązywanie Ax=b gdzie A=BB
T
sprowadza się do rozwiązania dwóch
układów
By=b
B
T
x=y.
W praktyce zakładamy Ŝe B = (b
ij
) jest macierzą trójkątną dolną.
Równanie A=BB
T
pozwala zapisać dla i
≤
j
a
b b
ij
ik
jk
j
=
∑
32
gdyŜ A jest symetryczna (b
jk
=0 dla k>j)
Stad
co pozwala na wyznaczenia pierwszej kolumny macierzy B.
ij
ik
jk
k
=
∑
1
b
a
b
a
b
i
n
i
i
11
11
1
1
11
2
=
=
=
,...,
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Cholesky’ego - algorytm
Postępując podobnie, konstruujemy macierz B „kolumna po kolumnie”.
Kolumna j macierzy jest obliczana w następujący sposób:
a
b
b
b
b
a
b
j
n
jj
jk
k
j
jj
jk
k
j
jj
jj
jk
k
j
=
= +
=
−
=
=
=
−
=
−
∑
∑
∑
2
1
2
2
1
1
2
1
1
2,...,
a
b b
b b
b b
b
b
a
b b
i
j
n
j
n
ij
ik
jk
k
j
ij
jj
ik
jk
k
j
ij
jj
ij
ik
jk
k
j
=
=
+
=
−
F
HG
I
KJ
= +
=
=
=
−
=
−
∑
∑
∑
1
1
1
1
1
1
1
2
,...,
,...,
33
MoŜna wykazać, Ŝe w linii j element b
jj
jest dodatni.
Oznaczmy przez A
j
i B
j
macierze rzędu j utworzone przez j pierwszych linii i
kolumn macierzy A i B (odpowiednio).
PoniewaŜ przez rekurencje b
jj
2
jest dodatnie dla 1
≤
i
≤
j-1, wówczas mamy
b
jj
2
>0 i
A
B B
A
B
j
j
j
T
j
j
=
=
=
>
=
∏
det(
)
det(
)
2
2
1
0
b
ii
i
j
det( )
A
=
=
∏
b
ii
i
n
2
1
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Choleskiego
–
liczba działań
Obliczenie b
jj
wymaga :
(j-1) dodawań+ (j-1) mnoŜeń+ 1 pierwiastek kwadratowy
Obliczenie b
ij
wymaga :
(j-1) dodawań+ (j-1) mnoŜeń+ 1 dzielenia
Na j-tym etapie :
(m-j+1)(j-1) dodawań, (m-j+1)(j-1) mnoŜeń, (m-1) dzieleń, 1 pierwiastek kwadratowy
(
1)
(
)
2
n
n
pierwiastkow kwdratowych
n n
n
j
dzielen
−
− =
∑
34
W sumie rozkład wymaga:
PoniewaŜ rozwiązanie układów By=b, B
T
x=y wymaga :
2n dzieleń, n(n-1) dodawa
ń
, n(n-1) mnoŜeń
W sumie działań elementarnych.
1
2
2
2
2
2
(
1)
(
1)(
1)
6
(
1)
(
1)(
1)
6
j
n
j
n
j
n n
n
j
j
dodawan
n n
n
j
j
mnozen
=
=
=
−
− +
− =
−
− +
− =
∑
∑
∑
N
n
n
n
CHOLESKY
=
+
−
2
15
6
3
2
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Bilans metod bezpośrednich
metoda
Macierz A
Liczba działa
ń
n=10
n=100
Diagonalna
n
10
100
Trójk
ą
tna
n
2
100
10000
Gauss
Dowolna
(4n
3
+9n
2
-7n)/6
805
681550
Jordan
Dowolna
n
3
+n
2
/2-n/2
1045
1004950
Choleski
Symetryczna
(2n
3
+15n
2
-n)/6
585
358317
35
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metody iteracyjne
Metody iteracyjne konstruują ciąg x
(i)
, i=1,2,… zbieŜny do rozwiązania
dokładnego x*. Otrzymywane jest w ten sposób rozwiązanie przybliŜone.
Metody iteracyjne są bardzo skuteczne w przypadku układów duŜego
rozmiaru (lub w przypadku macierzy rzadkich – mających wiele elementów
zerowych).
Zasada: Niech dany będzie ciąg wektorów x
(1)
, x
(2)
, …, x
(...)
zdefiniowany jako
x
(0)
– wektor zadany,
36
x
(0)
– wektor zadany,
x
(k+1)
= Bx
(k)
+c, k=1,2,…
gdzie B jest macierzą kwadratową (regularną) a c jest wektorem.
Ciąg ten będzie wykorzystywany do rozwiązania w sposób iteracyjny układu
Ax=b. Macierz B i wektor c będą zdefiniowane przez taki rozkład A aby ciąg
x
(k)
zbiegał się do rozwiązania układu wyjściowego.
RóŜne rozkłady macierzy A prowadza do róŜnych metod iteracyjnych.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
ZbieŜność metod iteracyjnych
Metoda iteracyjna jest zbieŜna jeŜeli
dla dowolnego x
(0)
.
Twierdzenie: Metoda iteracyjna jest zbieŜna jeŜeli B
k
0 dla k
∞
.
Dowód :
x
(k+1)
=Bx
(k)
+c
x*=Bx*+c (w granicy nieskończonej)
x
(k+1)
-x*= B(x
(k)
-x*) gdzie
e
(k+1)
=Be
(k)
(dla dowolnego k )
e
k
k
( )
=
− →
→∞
x
x
(k)
*
0
37
x
(k)
-x*= B
k
(x
(0)
-x*) gdzie
e
(k)
=B
k
e
(0)
ZbieŜność będzie miała miejsce wtedy i tylko wtedy gdy
Twierdzenie 2: Warunkiem
koniecznym i wystarczającym Ŝeby metoda iteracyjna była
zbieŜna jest aby promień spektralny macierzy B był mniejszy od jedności
ρ
(B)<1.
Warunkiem koniecznym i wystarczającym Ŝeby metoda iteracyjna była zbieŜna jest ||B||<1
dla dowolnej normy macierzowej.
lim
0
lim
0
k
k
czyli
→∞
→∞
∀
−
=
=
(0)
k
(0)
k
x
B (x
x*)
B
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metody rozkładu
JeŜeli zapiszemy macierz A wyjściowego układu równań w postaci A=M–N
gdzie M – macierz ”łatwo” odwracalna (to jest diagonalna lub trójkątna) to
moŜna dokonać rozkładu
Mx
(k+1)
=Nx
(k)
+b
gdzie x
(0)
dowolne natomiast
x
(k+1)
=M
-1
Nx
(k)
+M
-1
b=Bx
(k)
+c.
38
x
(k+1)
=M
-1
Nx
(k)
+M
-1
b=Bx
(k)
+c.
RóŜne rodzaje rozkładu macierzy A prowadzą do metody iteracyjnych:
metoda Jacobi’ego
metoda Gauss’a-Seidel’a
metoda kolejnych relaksacji
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Jacobi’ego
Niech A=(a
ij
) będzie macierzą rzędu n taka, Ŝe a
ii
≠
0 (i=1,…,n).
Przyjmijmy
A=D-(L+U)
gdzie D=(d
ij
)= a
ij
δ
ij
(diagonalna)
L=(l
ij
)= -a
ij
dla i>j, 0 w przeciwnym razie („trójkątna” dolna*)
U=(u
ij
)= -a
ij
dla i<j, 0 w przeciwnym razie („trójkątna” górna *)
39
U=(u
ij
)= -a
ij
dla i<j, 0 w przeciwnym razie („trójkątna” górna *)
*) zera na przekątnej
Taki rozkładowi moŜe towarzyszyć metoda iteracyjna
Dx
(k+1)
=(L+U)x
(k)
+b
lub w postaci n równań skalarnych
a x
a x
b
i
n
ii
i
k
ij
j
k
i
j
j i
n
(
)
( )
+
=
≠
= −
+
≤ ≤
∑
1
1
1
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Jacobi’ego - algorytm
Ciąg x
(k)
o wyrazach składowych
x
1
(k)
, x
2
(k)
,…, x
n
(k)
jest zdefiniowany jako:
test zatrzymania
(0)
(
1)
( )
1
ln
1
i
n
k
k
i
i
ij
j
j
ii
j i
x
dowo
y
x
b
a x
a
+
=
≠
=
−
∑
x
x
x
x
k
k
i
i
k
i
k
(
)
( )
(
)
( )
sup
+
∞
+
−
=
−
≤
1
1
ε
40
test zatrzymania
Macierz J=D
-1
(L+U) jest nazywana macierzą Jacobi’ego stowarzyszona z A.
ZbieŜność metody:
Warunkiem wystraczającym aby metoda Jacobi’ego była zbieŜna jest aby A
miała elementy diagonalne dominujące
.
Uwaga : moŜna dokonać permutacji linii aby zapewnić w/w wymagania.
i
i
i
∞
1
1
n
ii
ij
j
j i
a
a
dla
i
n
=
≠
>
≤ ≤
∑
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda Jacobi’ego
41
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Przykład 6
Metoda Gauss’a-Seidel’a
ZałóŜmy, Ŝe a
ii
≠
0 dla kaŜdego i, wówczas (D-L)
-1
istnieje.
Przyjmujemy
A=(D-L)-U
(D-L)x
(k+1)
=Ux
(k)
+b
Algorytm
Ciąg x
(k)
jest zdefiniowany jako
(0)
1
(
1)
(
1)
( )
1
i
i
n
k
k
k
x
dowolne
a x
a x
a x
b
i
n
−
+
+
+
= −
+
≤ ≤
∑
∑
42
Ciąg x jest zdefiniowany jako
test stopu : podobnie jak w metodzie Jacobi’ego
ZbieŜność: Warunek konieczny i wystarczający
ρ
(B)<1 jest trudny do
sprawdzenia w praktyce – zadowalamy sie warunkiem koniecznym.
Twierdzenie : jeŜeli macierz A ma elementy diagonalne dominujące, wówczas
metoda Gauss’a-Seidel’a jest zbieŜna.
Macierz G=(D-L)
-1
U jest nazywana macierzą Gaussa –Seidela stowarzyszoną z A.
(
1)
(
1)
( )
1
1
1
k
k
k
ii
i
ij
j
ij
j
i
j
j i
a x
a x
a x
b
i
n
+
+
=
= +
+
= −
+
≤ ≤
∑
∑
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda kolejnych relaksacji
W algorytmach poprzednio opisanych metod wprowadza się parametr
rzeczywisty
ω≠
0 aby polepszyć zbieŜność (SOR – Successive Over
Relaxation).
Np. w metodzie Gauss’a–Seidel’a
Obliczamy za pomocą algorytmu „klasycznego” metody wartość
i
n
−
L
M
O
P
1
1
43
i budujemy nowe przybliŜenie
x
a
a x
a x
b
i
n
i
k
ii
ij
j
k
j
i
ij
j
k
j i
n
i
(
)
(
)
( )
+
+
=
−
= +
= −
+
−
L
NM
O
QP
≤ ≤
∑
∑
1
1
1
1
1
1
1
x
x
x
x
i
k
i
k
i
k
i
k
(
)
( )
(
)
( )
+
+
=
+
−
1
1
ω
x
x
a
b
a x
a x
i
k
i
k
ii
i
ij
j
k
j
i
ij
j
k
j i
n
(
)
( )
(
)
( )
+
+
=
−
=
=
+
−
−
L
NM
O
QP
∑
∑
1
1
1
1
ω
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Metoda kolejnych relaksacji
Metoda odpowiada rozkładowi
Macierz R
ω
ω
ω
ω
=[D-
ω
L]
-1
[(1-
ω
)D+
ω
U] jest nazywana macierzą kolejnych relaksacji
stowarzyszoną z A.
ZbieŜność metody:
1
1
ω
ω
ω
ω
D
L x
D
U x
b
(k 1)
(k)
−
= −
+
L
NM
O
QP
+
+
44
ZbieŜność metody:
Twierdzenie :
Dla kaŜdej macierzy A, promień spektralny R
ω
ω
ω
ω
spełnia warunek
ρ
( R
ω
ω
ω
ω
)>|
ω
-1|.
Twierdzenie :
Warunkiem koniecznym zbieŜności metody kolejnych
relaksacji jest 1<
ω
<2.
Uwaga : przyjmując
ω
=1 otrzymujemy metodę Gauss’a-Seidel’a.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Uwarunkowanie macierzy
Rozwiązanie układu równań
Ax=b
na komputerze polega w rzeczywistości na
rozwiazywaniu układu przybliŜonego
(
A
+
dA
)
x
=
b
+
db
gdzie
dA
– macierz zaburzeń współczynników
A
,
db
– wektor zaburzeń
b
.
45
Nawet jeśli elementy macierzy są dokładne, ich reprezentacja binarna w
pamięci komputera pociąga za sobą błędy zaokrągleń.
Mówimy Ŝe problem jest źle uwarunkowany, jeŜeli małe zaburzenia danych
wejściowych powodują stosunkowo duŜe zaburzenia rozwiązania.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Uwarunkowanie macierzy
Przykład:
ale
2
6
8 0
2
6 00001
8 00001
1 0
1 0
x
y
x
y
x
y
+
=
+
=
→
=
=
RS
T
.
.
.
.
.
46
Rozwiazywanie numeryczne takiego źle uwarunkowanego problemu nie ma
sensu jeśli błędy w definiowaniu zadania (czy błędy przybliŜenia) są rzędu
szóstej cyfry po przecinku.
2
6
8 0
2
5 99999
8 00002
x
y
x
y
+
=
+
=
→
RS
T
.
.
.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
10 0
2 0
x
y
=
= −
.
.
Uwarunkowanie wstępne macierzy
Niech || . || oznacza normę macierzową. Uwarunkowaniem regularnej
macierzy A nazywamy liczbę
cond(A)=|| A || || A
-1
|| .
Macierz jest „dobrze uwarunkowana” jeśli jej uwarunkowanie nie jest duŜo
większe od jedności (im bardziej cond(A) jest bliskie 1, tym szybciej zbiega
47
większe od jedności (im bardziej cond(A) jest bliskie 1, tym szybciej zbiega
się algorytm rozwiązania).
JeŜeli powyŜszy wskaźnik jest duŜy (zwłaszcza dla duŜych układów równań),
wówczas błędy zaokrąglenia, które się kumulują podczas obliczeń, mogą
całkowicie zaburzyć wynik.
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Uwarunkowanie wstępne macierzy
Uwarunkowanie macierzy A polega na zastąpieniu rozwiazywania Ax=b
rozwiazywaniem równowaŜnego układu równań C
-1
Ax=C
-1
b gdzie C
-1
musi
być dobrana tak, aby cond(C
-1
A) był znacznie mniejszy niŜ cond(A).
Teoretycznie najlepszym wyborem jest C
-1
=A
-1
, w praktyce trzeba znaleźć C
-1
jak najbliŜsze A
-1
bez prowadzenia jednak zbyt kosztownych obliczeń w celu
48
jak najbliŜsze A bez prowadzenia jednak zbyt kosztownych obliczeń w celu
wyznaczenia C
-1
.
Istnieją róŜne metody uwarunkowania wstępnego macierzy A :
SSOR Evans’a,
rozkład niepełny Choleskiego,
metody IC(n), …
M.Pyrz Metody numeryczne w mechanice – Układy równań liniowych 10.2011
Normy wektorowe i macierzowe
normy wektorowe
1/ 2
2
1
2
1
1
1/
1
1
1
sup
n
n
i
i
i
i
p
n
p
i
i
p
i n
i
x
norma
x
norma euklidesowa
x
norma nieskonczona
x
=
=
∞
≤ ≤
=
=
=
=
=
∑
∑
∑
x
x
x
x
49
Normy macierzowe (stowarzyszone z normami wektorowymi)
A
Ax
Ax
x
A
Ax
x
A
Ax
x
A A
A
Ax
x
x C
x
x
x
x
*
x
n
=
=
=
=
L
NM
O
QP
=
=
=
=
L
NM
O
QP
∈
=
≠
≠
≤ ≤
=
≠
∞
≠
∞
∞
≤ ≤
=
∑
∑
sup
sup
sup
sup
sup
(
)
sup
sup
/
1
0
1
0
1
1
1
1
2
0
2
2
1 2
0
1
1
j n
ij
i
n
i n
ij
j
n
a
a
ρ
M.Pyrz Metody numeryczne w mechanice – Układy równań nieliniowych 10.2011