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
Wyszukiwarka
Podobne podstrony:
MNM 8 2014 id 304166 Nieznany
MNM 1 2014
MNM 4 2014
MNM 2 2014
MNM-6-2014
MNM 7 2014 id 304165 Nieznany
MNM 9 2014 id 304167 Nieznany
MNM 8 2014 id 304166 Nieznany
więcej podobnych podstron