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
Wa\
na jest ponadto znajomość głównych zasad z których korzystają takie
procedury.
2
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:
f& kwadratowa A=(aij) i=1,& n, j=1,& n f& diagonalna (aij=0 dla i`"j, aii=1)
f& symetryczna (aij=aji) f& antysymetryczna (aij=-aji)
f& hermita (aij=a*ji  elementy poło\
one symetrycznie względem przekątnej są zespolone
sprzÄ™\
one)
f& trójkątna dolna L (elementy zerowe nad przekątną)
f& trójkątna dolna L (elementy zerowe nad przekątną)
Rysunek 1
Rysunek 1
f& trójkątna górna U (elementy zerowe pod przekątną)
f& sprzÄ™\
ona A* z A : a*ij=âij gdzie âij jest liczbÄ… zespolonÄ… sprzÄ™\
onÄ… z aij.
f& transponowana
det I = 1 det AT = det A det A* = det A
Wyznacznik macierzy n
det(Ä…A) = Ä… det A det(AB) = det A det B
(kwadratowej) rzędu n:
1
det(A-1) = gdy A-1 istnieje
det A
3
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
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-1AS , gdy\
det(S-1AS)=det(S-1) det(A) det (S).
Równanie Ax=b ma jednoznaczne rozwiązanie je\
eli det(A) `" 0 (x=A-1b).
Macierz diagonalna
Macierz diagonalna
n
det A =
aij=0 dla i`"j, wyznacznik "a
ii
det A `" 0 Ô! aii `" 0 "i, 1 d" i d" n
i=1
Macierze kwadratowe trójkątne
dolna aij=0 dla ij
n
det A =
"a
ii
Wyznacznik macierzy kwadratowej trójkątnej
i=1
Iloczyn dwóch macierzy trójkątnych (dolnych)
Niech L1, L2 będą macierzami trójkątnymi dolnymi rzędu n.
Rysunek 2
Ich iloczyn L3=L1L2 jest macierzą trójkątną dolną.
4
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Układ m algebraicznych równań liniowych
" Układ m równań z n niewiadomymi: x1, x2, ..., xn
a11x1 +...+ a1nxn = b1
a21x1 +...+ a2nxn = b2
21 1 2n n 2
...
am1x1 +...+ amnxn = bm
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
5
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Układ n algebraicznych równań liniowych
Ax = b
" W postaci macierzowej
a11 a12 ... a1n x1 b1
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ïÅ‚a
a22 ... a2n śł ïÅ‚x2 śł ïÅ‚b2 śł
21
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
=
=
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
... ... ... ... ... ...
... ... ... ... ... ...
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
an2 ... ann ûÅ‚ ðÅ‚xn ûÅ‚ ðÅ‚bn ûÅ‚
n1
ðÅ‚a
Układ n liniowych równań algebraicznych o n niewiadomych
6
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metody rozwiÄ…zywania
Problem: Rozwiązać Ax=b,
macierz A=(aij) (1 d" i d" n, 1 d" j d" n) jest znana,
wektor b=(bi) (1 d" i d" n) jest podany,
wektor x=(xi) (1 d" i d" n) jest poszukiwany (nieznany).
Zakładamy, \
e A ( o elementach Rn lub Cn) jest odwracalna tzn. det(A)`"0.
Zakładamy, \
e A ( o elementach Rn lub Cn) 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 {x1,x2, & } 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, ...
7
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metody bezpośrednie
Układ o macierzy diagonalnej A=(aij), aij=0 gdy i `" j
bi
RozwiÄ…zanie :
xi = dla 1 d" i d" n
aii
liczba potrzebnych operacji NPRZEKATNA= n dzieleń.
Przykład 1
8
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metody bezpośrednie
Układ o macierzy trójkątnej dolnej (aij=0 dla ii
1 d" i d" n
"a xj = bi
ij
j=1
RozwiÄ…zanie Ax=b
i-1
jest równoznaczne rozwiązaniu
1 d" i d" n
"a xj + aii xi = bi
ij
j=1
i-1
L O
L O
1
1
Rozwiazywanie "schodzÄ…c"
Rozwiazywanie "schodzÄ…c"
xi = bi - =
= -
M P
M P
" ij
"a xj i = 1,2,...,n
aii
j=1
N Q
liczba niezbędnych operacji: patrz przypadek macierzy trójkątnej górnej
9
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Rozwiązanie układu o macierzy trójkątnej
Przyklad 2
10
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metody bezpośrednie
Układ o macierzy trójkątnej górnej (aij=0 dla i>j)
n
"a xj = bi 1 d" i d" n
ij
j=i
n
"a xj + aii xi = bi 1 d" i d" n
ij
j=i+1
n
L O
1
xi = bi -
RozwiÄ…zywanie "wchodzÄ…c"
M P
"a xj i = n,n -1,...,1
ij
aii
aii
j=n+1
j=n+1
N Q
N Q
n dzielen
liczba niezbędnych
n(n -1)
1+ 2 +ï"+ (n -1) = dodawan
działań:
2
n(n -1)
1+ 2 +ï"+ (n -1) = mnozen
2
NTRIAN= n+n(n-1)/2+ n(n-1)/2=n2 działań elementarnych ( + - * / )
11
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.
A x = b U x = c
Metoda Gaussa odpowiada dekompozycji (rozkładowi) A=LU (wtedy c=L-1b)
U
A L
=
12
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 Ax1=b1, Ax2=b2 & to
korzystne jest najpierw przeprowadzenie rozkładu A=LU a następnie kolejne
rozwiazywania układów równań: Lz1=b1 z1 , Ux1=z1 x1, Lz2=b2 z2 ,
Ux2=z2 x2, & w których macierze L, U są stałe dla wszystkich układów równań.
Uwaga 2:
Je\
eli A jest symetryczna to poszukuje się raczej rozkładu typu LDLT gdzie D jest
macierzą diagonalną (uzyskując zmniejszenie liczby działań i niewiadomych) 
oferuje nam to na przykład metoda Cholesky ego (wymagająca n3/6 działań, w
przeciwieństwie do n3/3 w rozkładzie LU).
13
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 a11(1) `" 0.
( ( (
Dla linii 2 d" i d" n :
aij2) = aij1) + ri1a11) 2 d" j d" n 2 d" i d" n
j
mno\
ymy pierwsza linie
( ( 2)
a12) = a11) 1 d" j d" n, ai(1 = 0 2 d" i d" n
macierzy A(1) przez j j
(
ri1= -ai1(1)/a11(1)
bi(2) = bi(1) + ri1b11) 2 d" i d" n
i dodajemy wynik
(2) (
b1 = b11)
do i-tego równania
14
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
(1) (1 (1 (
L L O
a12) a1n) b11)
11
Równanie A(1)x=b(1)
Ma0 a22) ï" a22)O Mb P
P
(2 ( (2)
ï"
n 2
jest równowa\
ne A(2)x=b(2)
M P M P
A(2) = b(2) =
M M P
ï" ï" ï"
gdzie A(2)=M(1)A(1), b(2)=M(1)b(1).
Mï" an2) ï" ann)P Mb P
P
( (2 (2)
0 ï"
N 2 Q N n Q
L
Mr1 0 0 ï" 0O
P
1 0 ï" 0
Uwaga :
21
Mr P
M
M(1) = 0 1 det(A(2))=a(1)11Minor(a11(1)) `" 0
31
Mï" ï" ï" ï" 0P
P
ï" ï" co upowa\
nia do przejścia
Mr P
M P
0 0 ï" 1
do następnego etapu.
N n1 Q
Minor  wyznacznik macierzy kwadratowej powstałej z danej
macierzy przez skreślenie pewnej liczby wierszy i kolumn
15
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 bi będzie oznaczona ai,n+1.
(.) (.) (.) (.)
(.) (.) (.)
L a1n a1,n+1
a1n (.)
11
11
Ma a12 ï" a2n b1 O L a12 ï" a2n a2,n+1O
(.) (.) (.) (.)
(.) (.) (.) (.)
b2
21
21
21
21
Ma P
Maî" a22 ï" P
Ma a22 ï" P
Maî" a22 ï" a2n b2 P = Ma a22 ï" a2n a2,n+1P
P
A(.) b =
A(.) b = =
M P
M
M P
î" î" î"
î" î" î"
Ma P
Ma an2 ï" ann bn P M an2 ï" ann an,n+1P
P
(.) (.) (.) (.)
(.) (.) (.) (.)
ï"
ï"
M P
n1
N n1 Q
N Q
Macierz rozszerzona A(2) jest zdefiniowana jako
( ( (
aij2) = aij1) + ri1a11) 2 d" i d" n, 1 d" j d" n +1
j
( ( 2)
a12) = a11) 1 d" j d" n +1; ai(1 = 0 2 d" i d" n
j j
16
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Przykład 3
17
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
(1 (1 (1 (
L O L O
a11) a12) ï" ï" ï" a1n) b11)
M P M P
(2 ( (
0 a22) ï" ï" ï" a22) b22)
n
M P M P
Na k-tym etapie M P M P
î" ï" ï" ï"
M P M P
( ( (k
A(k)x=b(k) A(k) = b(k) =
akk -1)-1 ï" ï" akk -1) bk -1)
M -1,k -1,nP M -1 P
(k (k (k
M P M P
0 akk ) ï" akn ) bk )
M P M P
ï" ï" ï"
M P M P
(k (k (
(k ) (k ) (k )
M P M P
M P M P
0 0 0 ank ) ï" ann ) bnk )
ï"
N Q N Q
N Q N Q
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 rik= - aik(k)/akk(k) i dodanie
wyniku do i-tej linii dla numerów k+1 d" i d" n.
Macierz rozszerzona A(k+1) przybiera postać
( ( (k
aijk +1) = aijk ) + rikakj ) k +1 d" i d" n, k +1 d" j d" n +1
( (
aijk +1) = aijk ) 1 d" j d" n +1, 1 d" i d" k; pozostale elementy zerowe
18
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
L
M1 0 0 ï" ï" ï" 0O
0
Mï" 1 0 ï" ï" ï" 0P
P
M 0 1
Mï" ï" ï" ï" ï" ï" 0P
P
M(k) = 1 ï" ï" ï"
M P
M
+1,k
Mï" ï" ï" ï" ï" 1 0P
Mï" ï" 0 rkï" ï" 0 0P
P
ï" 1 0
M
Mï" ï" ï" P
M0 0 ï" rnk ï" 0 1P
P
N Q
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.
19
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa - wyznacznik
Liczby akk(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).
Obliczanie wyznacznika
Metoda Gaussa pozwala na Å‚atwe obliczenie wyznacznika macierzy A.
Poniewa\
det(M)=1, mamy :
n
(k )
det A = det A(n) = (-1)p
"a
kk
k =1
20
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
( (
aikk ) = aikk )
sup
Na k-tym etapie jako element podstawowy aik(k) (k d" i d" n) wybieramy
k d"id"n
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.
oraz linii w której znaleziono element maksymalny.
( (
Dobór elementu podstawowego  w podmacierzy
aijk ) = aijk )
sup
k d"id"n
Na k-tym etapie jako element podstawowy obieramy wartość
k d" jd"n
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 !
21
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
( ( (
aijk +1) = aijk ) + rikakjk ) k +1 d" i d" n, k +1 d" j d" n
dla stałego k:
(
bi(k +1) = bi(k ) + rikbk k ) k +1 d" i d" n
(n-k) dzieleń aby obliczyć rik
(n-k)2 dodawań i (n-k)2 mno\
eń aby obliczyć aij(k+1)
(n-k) dodawań i (n-k) mno\
eń aby obliczyć bi(k+1)
(n-k) dodawań i (n-k) mno\
eń aby obliczyć bi(k+1)
n-1
"(n - k) = n(n2-1)
k =1
n-1 n-1 n-1
2
(n - k)2 + (n - k) =
" "q + "q = n(n -1)(2n -1) + n(n3-1) = n(n23-1)
6
k =1 q=1 q=1
rozwiÄ…zanie « wchodzÄ…c
n dzieleń, n(n-1)/2 dodawań, n(n-1)/2 mno\
eń
22
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa
Metoda Gaussa wymaga ogółem
n(n -1) n(n +1)
n + = dzielen
2 2
n(n -1) n(n2 -1) 2n3 + 3n2 - 5n
+ = dodawan
2 3 6
n(n -1) n(n -1) 2n + 3n - 5n
n(n -1) n(n2 -1) 2n3 + 3n2 - 5n
+ = mnozen
+ = mnozen
2 3 6
4n3 + 9n2 - 7n
działań
NGAUSS =
6
(tj rzędu 2n3/3 dla odpowiednio du\
ego n).
Przykład: dla n=10 : NGAUSS=805
Metoda współczynników Cramera: (n+1)! mno\
eń
23
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Ä…
1 0
1
1
A x = b x = d
1
1
0
1
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).
24
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana - algorytm
(
a11)
j
(
Etap 1: normalizacja
a12) = 1 d" j d" n +1
j
(1
a11)
( ( 1) (
redukcja
aij2) = aij1) - ai(1 a12) 1 d" j d" n +1, 2 d" i d" n
j
(
a22)
j
(
Etap 2: normalizacja
a23) = 2 d" j d" n +1
j
(2
a22)
( ( 2) (
redukcja
redukcja
a(3) = a(2) - a(2)a23) 2 d" j d" n +1, 1 d" i d" n, i `" 2
aij3) = aij2) - ai(2 a(3) 2 d" j d" n +1, 1 d" i d" n, i `" 2
j
(
akjk )
(
Etap k: normalizacja
akjk +1) = k d" j d" n +1
(k
akk )
( ( ( (
redukcja
aijk +1) = aijk ) - aikk )akjk +1) k d" j d" n +1, 1 d" i d" n, i `" k
Końcowa postać macierzy A to macierz jednostkowa I.
W kolumnie (n+1) macierzy poszerzonej otrzymujemy rozwiÄ…zanie xi =
ai,n+1(n) gdzie 1 d" i d" n.
25
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana
Przykład 4
26
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda eliminacji Gaussa-Jordana  liczba działań
n(n +1)
dzielen
2
n(n2 -1)
dodawan
2
n(n -1)
n(n2 -1)
mnozen
mnozen
2
Ogółem NJORDAN=n3+n2/2-n/2.
Dla odpowiednio du\
ego n, liczba działań jest rzędu n3.
Przykład: dla n=10 : NJORDAN =1045.
27
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Obliczanie macierzy odwrotnej
Aby otrzymać wektory kolumnowe xi macierzy odwrotnej A-1 wystarczy
rozwiązać n układów równań Axi=ei gdzie ei jest itym wektorem bazowym
e1=[1 0 0 & 0], e2=[0 1 0 & 0], & en=[0 0 0 & 1].
Przykład 5
28
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
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 | e1 | e2 | & |en ].
Kolumny macierzy odwrotnej A-1 zajmują miejsca kolumn wstępnie
zajmowanych przez odpowiednie wektory bazowe ei .
(
akjk )
(
Algorytm:
akjk +1) = k d" j d" 2n
(k
akk )
( ( ( (
aijk +1) = aijk ) - aikk )akjk +1) k d" j d" 2n 1 d" i d" n, i `" n
29
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.
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= LDLT.
Dowód: Poniewa\
A jest symetryczna A=LDR=(LDR)T=RTDLT.
Jednoznaczność rozkładu daje L=RT, R=LT.
30
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Choleskiego
Rozkład regularny Choleskiego macierzy A przybiera postać A=BBT 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
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=LDLT (elementy na przekątnej są dodatnie).
31
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Choleskiego - algorytm
Rozwiązywanie Ax=b gdzie A=BBT sprowadza się do rozwiązania dwóch
układów By=b BTx=y.
W praktyce zakładamy \
e B = (bij) jest macierzą trójkątną dolną.
j
Równanie A=BBT pozwala zapisać dla i d" j
aij =
"
"b bjk
ij ik jk
ik
k =1
gdy\
A jest symetryczna (bjk=0 dla k>j)
Stad
b11 = a11
ai1
bi1 = i = 2,...,n
b11
co pozwala na wyznaczenia pierwszej kolumny macierzy B.
32
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:
j j-1
j j-1
2 2
aij =
a =
"b bjk =bijbjj +"b bjk
ik ik
jj "b =b2 + "b
jk jj jk
k =1 k =1
k =1 k =1
j-1
j-1
i = j +1,...,n
Fa -"b bjkI
1
2
bij =
bjj = ajj -
G J
"b j = 2,...,n
jk
j = 2,...,n
bjj ij k =1 ik
H K
k =1
Mo\
na wykazać, \
e w linii j element bjj jest dodatni.
Oznaczmy przez Aj i Bj macierze rzędu j utworzone przez j pierwszych linii i
kolumn macierzy A i B (odpowiednio).
j
2
Aj = BjBT det(Aj) = det(Bj)2 =
j "b > 0
ii
i=1
Poniewa\
przez rekurencje bjj2 jest dodatnie dla 1 d" i d" j-1, wówczas mamy
n
bjj2 >0 i
2
det(A) =
"b
ii
i=1
33
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Choleskiego  liczba działań
Obliczenie bjj wymaga : (j-1) dodawań+ (j-1) mno\
eń+ 1 pierwiastek kwadratowy
Obliczenie bij 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
n pierwiastkowkwdratowych
n
"
"(n - j) = n(n2-1) dzielen
2
j=1
W sumie rozkład wymaga:
n
"(n - j +1)( j -1) = n(n2 -1) dodawan
6
j=2
n
"(n - j +1)( j -1) = n(n2 -1) mnozen
6
j=2
Poniewa\
rozwiązanie układów By=b, BTx=y wymaga :
2n dzieleń, n(n-1) dodawań, n(n-1) mno\
eń
2n3 +15n2 - n
W sumie działań elementarnych.
NCHOLESKY =
6
34
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 n2 100 10000
Gauss Dowolna (4n3+9n2-7n)/6 805 681550
Jordan Dowolna n3+n2/2-n/2 1045 1004950
Choleski Symetryczna (2n3+15n2-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,
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.
36
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Zbie\
ność metod iteracyjnych
Metoda iteracyjna jest zbie\
na je\
eli
e(k ) =kx(k) - x* 0
dla dowolnego x(0) .
"
Twierdzenie: Metoda iteracyjna jest zbie\
na je\
eli Bk 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 )
x(k)-x*= Bk(x(0)-x*) gdzie e(k)=Bke(0)
Zbie\
ność będzie miała miejsce wtedy i tylko wtedy gdy
"x(0) lim Bk (x(0) - x*) = 0 czyli lim Bk = 0
k" k"
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.
37
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-1Nx(k)+M-1b=Bx(k)+c.
x(k+1)=M-1Nx(k)+M-1b=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
38
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Jacobi ego
Niech A=(aij) będzie macierzą rzędu n taka, \
e aii`"0 (i=1,& ,n).
Przyjmijmy A=D-(L+U)
gdzie D=(dij)= aij´ij (diagonalna)
L=(lij)= -aij dla i>j, 0 w przeciwnym razie ( trójkątna dolna*)
U=(uij)= -aij dla iU=(uij)= -aij dla i*) 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
n
)
aii xi(k +1) = -
"a x(k + bi 1 d" i d" n
ij j
j=1
j`"i
39
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Jacobi ego - algorytm
xi(0) dowo ln y
Ciąg x(k) o wyrazach składowych
ëÅ‚ öÅ‚
n
1
)
ìÅ‚b -
x1(k), x2(k),& , xn(k) jest zdefiniowany jako:
xi(k +1) =
"a x(jk ÷Å‚
aii ìÅ‚ i j=1 ij ÷Å‚
ìÅ‚ ÷Å‚
j`"i
íÅ‚ Å‚Å‚
x( k +1) - x(k ) " = sup xi( k +1) - xi(k ) d" µ
i i
test zatrzymania
test zatrzymania
"
i
i
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
n
aii > aij dla1 d" i d" n
miała elementy diagonalne dominujące .
"
j=1
j`"i
Uwaga : mo\
na dokonać permutacji linii aby zapewnić w/w wymagania.
40
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Jacobi ego
Przykład 6
41
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda Gauss a-Seidel a
Załó\
my, \
e aii`"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
xi(0) dowolne
i-1 n
CiÄ…g x(k) jest zdefiniowany jako
CiÄ…g x jest zdefiniowany jako
+1)
+1)
a x(k +1) + a x(k ) + b 1 d" i d" n
aiixi(k +1) + aij x(jk ) + bi 1 d" i d" n
"a x(k = - "
"a x(jk = - "
ij
j=1 j=i+1
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)-1U jest nazywana macierzÄ… Gaussa  Seidela stowarzyszonÄ… z A.
42
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-1 n
i-1 n
L O
L" O
(k +1)
1
1
)
x = - aij x(k +1) +
i
M P
M P
"a x(jk - bi 1 d" i d" n
ij
aii j=1 j
j=i+1
N Q
i budujemy nowe przybli\
enie
(k +1)
xi(k +1) = xi(k ) + É x
i - xi(k )
i-1 n
L O
É
+1) )
xi(k +1) = xi(k ) + bi -
M P
"a x(k -"a x(k
ij j ij j
aii
j=1 j=i
N Q
43
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Metoda kolejnych relaksacji
Metoda odpowiada rozkładowi
1
L
(k)
D - É L x(k+1) =
M1- É D + UOx + b
P
É É
N Q
Macierz RÉ=[D-ÉL]-1[(1-É)D+ÉU] jest nazywana macierzÄ… kolejnych relaksacji
É
É
É
stowarzyszonÄ… z A.
Zbie\
ność metody:
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.
44
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.
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 z
le uwarunkowany, je\
eli małe zaburzenia danych
wejściowych powodują stosunkowo du\
e zaburzenia rozwiÄ…zania.
45
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Uwarunkowanie macierzy
Przykład:
2x + 6y = 8.0 x = 1.0
R

S
2x + 6.00001y = 8.00001 y = 1.0
T
ale
2x + 6y
R x = 10.0

S2x + 5.99999y= 8.0
= 8.00002
y = -2.0
T
Rozwiazywanie numeryczne takiego z
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.
46
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
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
większe od jedności (im bardziej cond(A) jest bliskie 1, tym szybciej zbiega
siÄ™ algorytm rozwiÄ…zania).
Je\
eli powy\
szy wskaz
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.
47
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-1Ax=C-1b gdzie C-1 musi
być dobrana tak, aby cond(C-1A) był znacznie mniejszy ni\
cond(A).
Teoretycznie najlepszym wyborem jest C-1=A-1, w praktyce trzeba znalez
ć C-1
jak najbli\
sze A bez prowadzenia jednak zbyt kosztownych obliczeń w celu
jak najbli\
sze A-1 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), &
48
M.Pyrz Metody numeryczne w mechanice  Układy równań liniowych 10.2011
Normy wektorowe i macierzowe
normy wektorowe
1/ 2
n n
îÅ‚ Å‚Å‚
x = xi norma 1 x = xi2 śł norma euklidesowa
" "
ïÅ‚
1 2
i=1 ðÅ‚ i=1 ûÅ‚
1/ p
n
îÅ‚ Å‚Å‚
x = sup xi norma nieskonczona x = xip śł
"
ïÅ‚
" p
1d"id"n
ðÅ‚ i=1 ûÅ‚
Normy macierzowe (stowarzyszone z normami wektorowymi)
n
Ax
Ax
L" aij O
1
A = sup Ax = sup A = sup = sup
M P
1
x x
x `"0 x `"0 1d" jd"n
x"Cn
Ni=1 Q
1
x =1
n
Ax Ax L" aij O
1/2
2 "
A = sup = Á(A*A) A = sup = sup
M P
2 "
x x
x `"0 x `"0 1d"id"n
j=1
N Q
2 "
49
M.Pyrz Metody numeryczne w mechanice  Układy równań nieliniowych 10.2011