Równanie postaci
a1x1 + a2x2 + ... + anxn = b
nazywamy równaniem liniowym o n niewiadomych x1, x2, ..., xn. Prawa strona równania b to tzw. wyraz wolny.
Równanie liniowe postaci
a1x1 + a2x2 + ... + anxn = 0
nazywamy równaniem liniowym jednorodnym.
Rozpatrzmy dalej układ m równań liniowych o n niewiadomych:
m n mn m
m
n n
n n n
b x a x
a x a
b x a x
a x a
b x a x
a x a
= +
+ +
= +
+ +
= +
+ +
K
K K K K K K K
K K
2 2 1 1
2 2
2 22 1 21
1 2
12 1 11
Układ taki moŜemy zapisać w postaci macierzowej:
= ,
gdzie to macierz m × n współczynników przy niewiadomych
=
mn m
m
n n
a a
a
a a
a
a a
a
L M O M M
L L
2 1
2 22
21
1 12
11
,
to wektor n niewiadomych a – wektor m wyrazów wolnych (prawych stron równań):
=
xn
x x
M
2 1
,
=
bm
b b
M
2 1
.
Jeśli = , to jest to jednorodny układ równań liniowych.
Jak widać, liczba niewiadomych nie musi być równa liczbie równań. Jeśli m > n (więcej równań niŜ niewiadomych), to układ nazywamy nadokreślonym. Jeśli m < n, (mniej równań niŜ niewiadomych, to układ nazywamy niedookreślonym).
Macierzą uzupełnioną macierzy nazywamy macierz m × (n+1) powstałą przez dopisanie do macierzy kolumny wyrazów wolnych układu:
=
m mn m
m
n n
b b b
a a
a
a a
a
a a
a
M L
M O M M
L L
2 1
2 1
2 22
21
1 12
11
.
Istnienie rozwiązań układu moŜna określić badając rzędy macierzy A i U. Przypomnijmy (wykład 1.3), Ŝe rząd macierzy jest równy liczbie liniowo niezaleŜnych kolumn (lub wierszy).
Dla praktycznego zastosowania dogodniejsze będzie równowaŜne określenie rzędu macierzy:
Podmacierzą kwadratową danej macierzy nazywamy dowolną macierz kwadratową, jaką moŜna otrzymać wykreślając z naszej macierzy wiersze i kolumny (w szczególności macierz kwadratowa jest swoją podmacierzą powstałą przez wykreślenie zero kolumn i zero wierszy).
Rząd macierzy jest równy stopniowi największej kwadratowej podmacierzy o niezerowym wyznaczniku.
Prawdziwe jest następujące twierdzenie (Kroneckera-Capellego):
Niech A będzie macierzą współczynników równań układu m równań liniowych z n niewiadomymi a U macierzą uzupełnioną tego układu.
1. Jeśli rank(A) = rank(U) = n, to istnieje dokładnie jedno rozwiązanie tego układu (układ jest oznaczony)
2. Jeśli rank(A) = rank(U) = r < n, to układ ma nieskończenie wiele rozwiązań (jest nieoznaczony) dających się wyrazić za pomocą n – r parametrów.
3. Jeśli rank(A) < rank(U), to układ nie ma rozwiązań (jest sprzeczny).
ZauwaŜmy, Ŝe warunek z pkt. 2 oznacza, Ŝe w układzie da się znaleźć podmacierz stopnia co najwyŜej r o niezerowym wyznaczniku. n – r zmiennych, które nie weszły do wyznacznika (z kolumn wykreślonych przy tworzeniu podmacierzy) moŜna potraktować jako parametry słuŜące do określenia pozostałych r zmiennych.
Jeśli jakieś równania nie weszły do wyznacznika, z którego wyznaczono rząd A (zostały wykreślone przy tworzeniu podmacierzy), to moŜna je pominąć (bo są liniowo zaleŜne od pozostałych a więc wynikają z nich).
2.2 Jednorodny układ równań liniowych.
RozwaŜania dotyczące istnienia rozwiązań z pkt. 2.1 dotyczą zarówno jednorodnych jak i niejednorodnych układów równań liniowych. PoniewaŜ układy jednorodne często pojawiają się w problemach fizycznych, warto tu w sposób jawny wyartykułować wnioski z twierdzenia Kroneckera-Capellego.
Jednorodny układ m równań liniowych z n niewiadomymi ma postać:
0 0 0
2 2 1 1
2 2
22 1 21
1 2
12 1 11
= +
+ +
= +
+ +
= +
+ +
n mn m
m
n n
n n
x a x
a x a
x a x
a x a
x a x
a x a
K
K K K K K K K
K K
,
czyli
Ax = 0.
ZauwaŜmy, Ŝe rank(A) = rank(U), bo dołoŜenie do macierzy kolumny z samych zer nie moŜe zwiększyć jej rzędu. Wynika stąd wniosek, Ŝe na pewno istnieje jakieś rozwiązanie układu niejednorodnego.
Rzeczywiście:
Jednorodny układ równań liniowych z n niewiadomymi ma zawsze rozwiązanie x1 = x2 = ... = xn = 0.
Jest to tzw. rozwiązanie trywialne.
Jeśli rank(A) = n, to rozwiązanie trywialne jest jedynym, jeśli natomiast rank(A) < n, to istnieją nietrywialne rozwiązania układu.
Inny waŜny wniosek dotyczy układu z kwadratową macierzą A:
Jednorodny układ n równań z n niewiadomymi ma nietrywialne rozwiązanie wtedy i tylko wtedy, gdy det(A) = 0.
2.3 Wzory Cramera.
Zajmiemy się teraz układem n równań liniowych z n niewiadomymi:
Ax = b
(w szczególności moŜe być teŜ b = 0).
Oznaczmy
W = det(A) =
nn n
n
n n
a a
a
a a
a
a a
a
L M L M M
L L
2 1
2 22
21
1 12
11
nn n
n
n
a a
b
a a
b
a a
b
L M L M M
L L
2
12 22
2
1 12
1
W1 = ,
nn n
n
n n
a b
a
a b
a
a b
a
L M L M M
L L
1
2 2
21
1 1
11
W2 = , … Wn =
n n
n a b
a
b a
a
b a
a
L M L M M
L L
2 1
2 22
21
1 12
11
(czyli Wi to wyznacznik macierzy powstałej przez zastąpienie i-tej kolumny macierzy A kolumną wyrazów wolnych.)
Jeśli W ≠ 0, to układ ma dokładnie jedno rozwiązanie, które zadane jest wzorami Cramera:
W Wi
xi = . (tw. Cramera)
Jeśli W = 0 i przynajmniej jeden z wyznaczników Wi jest róŜny od zera, to układ jest sprzeczny.
Jeśli W = W1 = W2 = ... = Wn = 0, to równania są liniowo zaleŜne i układ moŜe być sprzeczny lub nieoznaczony. W tym drugim przypadku eliminując równania zaleŜne liniowo od pozostałych dostajemy układ niedookreślony.
Ze względu na koszt obliczania wyznaczników (por. 1.4) wzory Cramera nie nadają się do praktycznego zastosowania przy rozwiązywaniu układów równań liniowych z wyjątkiem małych wartości n.
ZauwaŜmy jeszcze, Ŝe gdy W (czyli det(A) ) ≠ 0, to istnieje A-1 i teoretycznie nasz układ Ax = b
moŜna rozwiązać uŜywając macierzy odwrotnej do A A-1Ax = A-1b
czyli
x = A-1b
PoniewaŜ jednak obliczenie macierzy odwrotnej z wykorzystaniem macierzy dopełnień algebraicznych wymaga obliczania wyznaczników, takŜe i ten sposób nie nadaje się do praktycznego zastosowania.
2.4 Eliminacja Gaussa. Rozkład LU.
Zajmijmy się niejednorodnym układem n równań liniowych z n niewiadomymi
n n nn n
n
n n
n n n
b x a x
a x a
b x a x
a x a
b x a x
a x a
= +
+ +
= +
+ +
= +
+ +
K
K K K K K K K
K K
2 2 1 1
2 2
2 22 1 21
1 2
12 1 11
czyli
Ax = b.
Jego macierz uzupełniona ma postać:
=
n nn n
n
n n
b b b
a a
a
a a
a
a a
a
M L
M O M M
L L
2 1
2 1
2 22
21
1 12
11
U .
(pionową linią dla przejrzystości oddzieliliśmy kolumnę wyrazów wolnych) Wykonując na układzie następujące operacje:
- zamiana dwu równań miejscami
- pomnoŜenie dowolnego równania przez niezerową stałą
- dodanie do dowolnego równania innego równania pomnoŜonego przez stałą otrzymujemy układ równowaŜny, tzn. mający to samo rozwiązanie, co wyjściowy.
Operacjom tym odpowiadają działania na wierszach macierzy uzupełnionej:
- zamiana dwu wierszy
- pomnoŜenie wiersza przez niezerową stałą
- dodanie do wiersza innego wiersza pomnoŜonego przez stałą.
MoŜemy w ten sposób próbować przekształcić nasz układ do układu równowaŜnego, dla którego łatwiej wyznaczyć rozwiązanie.
W tym celu od drugiego równania odejmijmy pierwsze pomnoŜone przez czynnik a21/a11, od równania trzeciego równanie pierwsze pomnoŜone przez a31/a11, itd. aŜ do równania n-tego, od którego odejmujemy pierwsze pomnoŜone przez an1/a11. W ten sposób przy pomocy równania pierwszego wyzerowaliśmy współczynnik przy x1 w kolejnych równaniach (czyli w pierwszej kolumnie macierzy U mamy zera z wyjątkiem pierwszego wiersza).
Następnie przy pomocy równania drugiego zerujemy współczynniki przy zmiennej x2 w równaniu trzecim i następnych. Procedurę tę powtarzamy dla kolejnych równań.
Współczynnik przy eliminowanej zmiennej to tzw. element główny.
W wyniku tej procedury zwanej eliminacją Gaussa sprowadzamy macierz układu równań do postaci trójkątnej górnej, tzn. wszystkie elementy pod główną przekątną są równe zero a niezerowe znajdują się na niej i powyŜej. (Tak będzie, gdy det(A) ≠ 0, czyli układ ma rozwiązanie; w przeciwnym wypadku pojawi się wiersz(e) z samych zer – patrz ćwiczenia).
Dla układu równań w tej postaci łatwo juŜ znaleźć rozwiązanie: zaczynamy od ostatniego równania, przy którym mamy niezerowy współczynnik tylko przy ostatniej zmiennej, co
pozwala łatwo ją wyznaczyć. Wartość tej zmiennej podstawiamy do równania przedostatniego, dzięki czemu moŜemy obliczyć przedostatnią zmienną. Kontynuując podstawienia obliczonych zmiennych wyznaczamy z kolejnych równań następne aŜ do równania pierwszego. Proces ten nazywamy podstawieniem wstecz.
Nietrudno zauwaŜyć, Ŝe tak samo łatwo moŜna rozwiązać układ równań liniowych z macierzą trójkątną dolną – zaczynamy wtedy od pierwszego równania (a procedura to podstawienie w przód).
Procedura eliminacji Gaussa pozwala rozwiązać lub stwierdzić brak rozwiązań nie tylko dla układu równań z macierzą kwadratową, ale dla dowolnego układu m równań liniowych z n niewiadomymi (patrz ćwiczenia).
W trakcie eliminacji moŜe się okazać, Ŝe współczynnik przy zmiennej, którą chcemy w danym kroku eliminować z kolejnych równań (czyli element główny) jest równy 0. Wtedy odszukujemy wśród następnych równań takie, przy którym współczynnik przy interesującej nas zmiennej jest róŜny od zera, zamieniamy równania miejscami i kontynuujemy (jeśli nie da się takiego znaleźć, to znaczy, Ŝe det(A) = 0).
W numerycznej realizacji tego algorytmu zamianę wierszy w macierzy (czyli przestawienie równań) wykonuje się z jeszcze innego powodu. OtóŜ analiza błędów pokazuje, Ŝe najkorzystniej dla dokładności rozwiązania jest wybierać w kaŜdym kroku największy co do wartości bezwzględnej element główny. Najczęściej odbywa się to właśnie przez zamianę równań (czyli tzw. częściowy wybór elementu głównego w kolumnie), ale moŜna wybierać element główny w wierszu (zamieniając kolumny, co odpowiada przenumerowaniu zmiennych) albo łączyć te obie strategie (pełny wybór elementu głównego).
ZauwaŜmy jeszcze, Ŝe współczynniki, przez jakie mnoŜyliśmy równania zaleŜą tylko od macierzy współczynników A a nie zaleŜą od wyrazów wolnych b. Oznacza to, Ŝe gdy mamy do rozwiązania kilka układów róŜniących się jedynie prawymi stronami, to macierz A przekształcamy tylko raz (oczywiście musimy przekształcić wektory b, ale jest to znacznie mniej kosztowne). Eliminacja Gaussa pozwala zatem efektywnie rozwiązywać układy o wspólnej macierzy współczynników.
Eliminacja Gaussa ściśle wiąŜe się z rozkładem macierzy na iloczyn macierzy trójkątnych.
Jeśli A jest macierzą nieosobliwą (jej wyznacznik jest róŜny od 0), to istnieje jej jednoznaczny rozkład na czyniki:
A = LU,
gdzie L jest macierzą trójkątną dolną o elementach diagonalnych równych 1 a U jest macierzą trójkątną górną (rozkład trójkątny). (Bez określenia elementów diagonalnych którejś macierzy rozkład nie byłby jednoznaczny).
Rozkład taki moŜna otrzymać przez eliminację Gaussa (macierz U to macierz układu równań po eliminacji a macierz L tworzą mnoŜniki eliminacji). ZauwaŜmy jednak, Ŝe jeśli w czasie eliminacji zamienialiśmy wiersze, to nie otrzymaliśmy rozkładu wyjściowej macierzy A, tylko rozkład macierzy, która w porównaniu do A ma spermutowane wiersze.
Eliminacja Gaussa prowadzi zatem do rozkładu PA = LU,
gdzie P to macierz permutacji. Iloczyn PA daje macierz z poprzestawianymi wierszami.
(Macierz P ma w kaŜdym wierszu i w kaŜdej kolumnie dokładnie jeden element równy 1, pozostałe są równe 0).
Jeśli znamy rozkład trójkątny macierzy A, to moŜemy wykorzystać go do rozwiązania układu Ax = b
dzięki temu, Ŝe łatwo rozwiązać układ z macierzą trójkątną górną (podstawienie wstecz) lub dolną (podstawienie w przód).
Mamy bowiem PA = LU i obliczamy
b’ = Pb
następnie wyznaczamy c rozwiązując układ
Lc = b’ (łatwe, bo L jest trójkątna dolna) a potem x rozwiązując układ
Ux = c (teŜ łatwe, bo U jest trójkątna górna)
Koszt eliminacji Gaussa jest rzędu n3 mnoŜeń, koszt rozwiązania układu równań przy znanym rozkładzie trójkątnym jest rzędu n2.
ZauwaŜmy jeszcze, Ŝe operacje, które wykonywaliśmy na macierzy A mogły najwyŜej zmienić znak jej wyznacznika (zamiana wierszy). Oznacza to, Ŝe wyznacznik iloczynu LU z dokładnością do znaku jest taki sam, jak wyznacznik A. PoniewaŜ wyznacznik iloczynu macierzy jest równy iloczynowi wyznaczników, wyznacznik macierzy trójkątnej jest iloczynem jej elementów diagonalnych a L ma na diagonali same jedynki, to z dokładnością do znaku wyznacznik A jest równy iloczynowi elementów diagonalnych macierzy U. Znak moglibyśmy ustalić na podstawie wyznacznika macierzy P (jest on równy ±1), ale zwykle nie tworzy się macierzy permutacji a jedynie zapamiętuje informacje o przestawieniach wierszy (transpozycje). Wystarczy potem policzyć, czy ich liczba była parzysta, czy nieparzysta.
Eliminacja Gaussa daje zatem moŜliwość efektywnego obliczania wyznaczników macierzy.
PokaŜmy jeszcze, Ŝe rozwiązywanie układów równań liniowych (eliminacja Gaussa) daje narzędzie do praktycznego wyznaczania macierzy odwrotnej.
RozwiąŜmy bowiem n układów równań:
Axi = bi
Wektor bi ma i-tą składową równą 1, pozostałe są równe 0. Ustawiając obok siebie wektory bi dostaniemy zatem macierz jednostkową I. Skoro wektory xi są rozwiązaniami powyŜszego układu równań, to ustawiając je obok siebie dostajemy macierz X, która spełnia równanie
AX = I
Ale to oznacza przecieŜ, Ŝe X = A-1. Zatem macierz odwrotną do A dostajemy ustawiając kolumnami rozwiązania układów równań liniowych z macierzą A i prawymi stronami, w których kolejna składowa jest równa 1 a reszta 0. Dzięki rozkładowi trójkątnemu wymaga to rzędu n3 mnoŜeń na rozkład + n · n2 mnoŜeń na rozwiązanie n układów równań.
Zapamiętajmy zatem wniosek: w praktyce numerycznej układu równań liniowych n × n nie rozwiązuje się obliczając wyznaczniki, ani nie wykorzystuje się do tego macierzy odwrotnej.
Przeciwnie, rozwiązywania układu równań uŜywa się do obliczenia wyznacznika lub znalezienia odwrotności macierzy.
2.5 Rozkład Cholesky’ego.
Omówiony wyŜej rozkład LU nie jest jedynym moŜliwym. Przy znajdowaniu rozkładu macierzy na iloczyn czynników trójkątnych często ułatwieniem bywa szczególna postać macierzy.
Rozkład A = LU macierzy symetrycznej dodatnio określonej (patrz 1.3) moŜna przedstawić w innej postaci wykorzystując zachodzącą równość
LT = D-1U, gdzie D = diag(di), di > 0 Wtedy
DLT = DD-1U DLT = U
i rozkład moŜna przedstawić w postaci symetrycznej A = LDLT
Przyjmując ˆ LD1/2
L= , gdzie D1/2 = diag
( )
dkmamy rozkład A = Lˆ LˆT (daszek ^ uŜyty tu został jedynie po to, by odróŜnić czynniki trójkątne dolne we wzorze wyŜej).
Symetryczny rozkład macierzy symetrycznej i dodatnio określonej postaci A = LLT
nazywamy rozkładem Cholesky’ego.
Dysponując rozkładem Cholesky’ego macierzy A moŜemy rozwiązać układ Ax = b
wyznaczając c z trójkątnego układu Lc = b
a następnie x rozwiązując układ trójkątny LTx = c
Wyznacznik macierzy A moŜna obliczyć jako kwadrat iloczynu diagonalnych elementów macierzy L
det(A) =
2
1
∏
= n
i
lii
II. WYBRANE ASPEKTY IMPLEMENTACJI ALGEBRY LINIOWEJ.
2.6 Ile operacji wymaga mnoŜenie macierzy?
Zastanówmy się, ile mnoŜeń wymaga pomnoŜenie dwu macierzy n × n. Ze wzoru
∑
=
= n
k
kj ik
ij a b
c
1
(*)
wynika, Ŝe dla obliczenia pojedynczego elementu macierzy wynikowej musimy wysumować iloczyny przechodząc przez wiersz pierwszej macierzy (kolumnę drugiej), czyli mamy n mnoŜeń. PoniewaŜ musimy obliczyć wszystkie n2 elementów, potrzebujemy wykonać n⋅n2 = n3 mnoŜeń. Organizacja obliczeń bezpośrednio odpowiadająca wzorowi na element macierzy oznacza więc koszt mnoŜenia macierzy rzędu n3.
Zobaczmy to na macierzach 2 × 2:
=
22 21
12 11 22 21
12 11 22
21 12 11
b b
b b a a
a a c
c c c
Potrzeba 23 = 8 mnoŜeń.
c11 = a11b11 + a12b21 c12 = a11b12 + a12b22 c21 = a21b11 + a22b21
c22 = a21b12 + a22b22
Jak się jednak okazuje, moŜna wykonać obliczenia wykonując jedynie 7 mnoŜeń, za to kosztem dodatkowych dodawań:
M1 = (a11 + a22)(b11+ b22) M2 = (a21 + a22)b11
M3 = a11(b12 – b22)
M4 = a22(b21 – b11) (**)
M5 = (a11 + a12)b22
M6 = (a21 – a11)(b11 + b12) M7 = (a12 – a22)(b21 + b22) wtedy
c11 = M1 + M4 – M5 + M7
c12 = M3 + M5
c21 = M2 + M4
c22 = M1 – M2 + M3 + M6
Redukcja liczby mnoŜeń o 1 nie wydaje się być warta zwiększenia liczby dodawań i komplikacji algorytmu, ale...
ZauwaŜmy, Ŝe jeśli macierz o wymiarze 2n × 2n podzielimy na bloki o wymiarze n × n, to wzór (*) pozostaje prawdziwy dla bloków macierzy:
=
22 21
12 11 22 21
12 11 22
21 12 11
B B
B B A A
A A C
C C
C .
MoŜemy zatem zastosować powyŜsze wzory do bloków macierzy.
Na tym opiera się algorytm Strassena. MnoŜąc macierze kwadratowe stopnia 2n (inne moŜemy uzupełnić kolumnami i wierszami z zer do potrzebnego rozmiaru) dzielimy je na
bloki stopnia 2n-1, dla których moŜemy zastosować wzory (**). W celu obliczenia potrzebnych iloczynów moŜemy te bloki znowu podzielić na bloki stopnia 2n-2 i zastosować (**). MoŜemy postępować tak rekurencyjnie, aŜ dojdziemy do mnoŜenia liczb (macierzy 1 × 1). W praktyce przy dostatecznie małych macierzach przechodzimy do klasycznego mnoŜenia macierzy. Istota redukcji czasu obliczeń zasadza się na tym, Ŝe zyskujemy jedno mnoŜenie z ośmiu na kaŜdym etapie. Przy dostatecznie duŜych macierzach zysk na mnoŜeniach przewaŜa koszty dodatkowych dodawań.
Algorytm Strassena potrzebuje nlog27≈ n2.807 mnoŜeń. Widać zysk w porównaniu do klasycznego n3.
Najszybszy znany obecnie algorytm Coppersmitha-Winograda wymaga n2.376 mnoŜeń.
2.7 Biblioteki algebry liniowej. BLAS i LAPACK.
Zagadnienia algebry liniowej powszechnie pojawiają się w zadaniach numerycznych w naukach ścisłych, często dla bardzo duŜych zagadnień. Niezwykle waŜna jest zatem ich efektywna implementacja w językach programowania. Tworzenie potrzebnych procedur od podstaw przez kaŜdego uŜytkownika byłoby stratą czasu i prowadziłoby do mało wydajnych kodów. Problem ten rozwiązują powszechnie dostępne biblioteki procedur numerycznych, pozwalające na budowę efektywnie działających programów przy minimalnej znajomości metod numerycznych.
Ze względów historycznych znaczna część tych bibliotek napisana jest w Fortranie (język ten powstał do programowania obliczeń naukowych, na co wskazuje jego nazwa: FORmula TRANslator).
Podstawowe cegiełki do rozwiązywania zagadnień algebry liniowej zawiera biblioteka BLAS (Basic Linear Algebra Subprograms) napisana w Fortranie. Procedury biblioteki BLAS podzielone są na 3 poziomy:
BLAS level 1: operacje wektorowe typu αx + y, takŜe obliczanie iloczynów skalarnych i norm
BLAS level 2: operacje macierz-wektor typu αAx + βy, takŜe inne procedury, np.
rozwiązywanie układu Tx = y z trójkątną macierzą T BLAS level 3: operacje macierz-macierz typu αAB + βC
W bibliotece BLAS stosuje się konwencję kodowania w nazwie procedury informacji o precyzji zmiennych uŜywanych w obliczeniach, typu macierzy i wykonywanej operacji.
Operacje zdefiniowane są dla czterech dokładności (odpowiednie typy zmiennych w Fortranie):
S – rzeczywiste, pojedyncza dokładność D – rzeczywiste, podwójna dokładność C – zespolone, pojedyncza dokładność Z – zespolone, podwójna dokładność Przykładowe typu macierzy to:
GE – ogólna SY – symetryczna HE – hermitowska TR – trójkątna
Zaś przykładowe kody operacji to:
DOT – iloczyn skalarny
AXPY – suma wektorów αx + y MV – mnoŜenie macierz-wektor Ax SV – rozwiązanie układu macierz-wektor MM – mnoŜenie macierzy
SM – rozwiązanie układu macierz-macierz Przykłady:
poziom 1
SSCAL, DSCAL, CSCAL, ZSCAL obliczają αx w odpowiedniej dokładności i dla danegu typu danych
SDOT i DDOT obliczają iloczyn skalarny dwu rzeczywistych wektorów odpowiednio w pojedynczej lub podwójnej dokładności
DNRM2 i ZNRM2 obliczają w podwójnej dokładności normę euklidesową odpowiednio wektora rzeczywistego i zespolonego
poziom 2
SGEMV, DGEMV obliczają iloczyn Ax dla ogólnej macierzy A
CHEMV, ZHEMV obliczają iloczyn macierz –wektor dla macierzy hermitowskiej STRSV, DTRSV rozwiązują dla liczb rzeczywistych układ z macierzą trójkątną Tx poziom 3
SGEMM, DGEMM, CGEMM, ZGEEM – mnoŜenie macierzy, w tym procedura mnoŜenia macierzy rzeczywistych w podwójnej dokładności powszechnie stosowana w testach mierzących szybkość komputerów.
STRSM, DTRSM, CTRSM, ZTRSM rozwiązywanie układu trójkątnego z wieloma prawymi stronami
Operacja wykonywana przez procedurę moŜe być bardziej rozbudowana, niŜ moŜna by było przypuszczać z nazwy, np. DGEMM wykonuje operację
αAB + βC
przy czym jedna z macierzy A i B lub obie moŜe być wcześniej poddana transpozycji.
Zwykłemu obliczaniu AB odpowiada wyłączenie transponowania oraz ustawienie α=1, β=0.
Kody źródłowe procedur w Fortranie 77 dostępne są na stronie Netlib http://www.netlib.org/blas/
Tym niemniej nie naleŜy ich kompilować samodzielnie. Pod róŜnymi systemami operacyjnymi dostępne są bowiem skompilowane wersje biblioteki zapewniające lepszą wydajność. Oprócz bibliotek dostępnych pod większość wersji Linuxa istnieją takŜe biblioteki zawierające BLASa zoptymalizowane przez producenta konkretnej architektury, np.
– Intel Math Kernel Library (procesory Intel) – AMC Core Math Library (AMD)
– implementacje firm HP, SGI czy Sun Istnieją takŜe interfejsy BLASa do języka C.
Rozwiązania podstawowych problemów algebry liniowej zawiera biblioteka LAPACK (Linear Algebra PACKage) stworzona w języku Fortran 77. Obecnie istnieje wersja w
Fortranie 95 oraz implementacje w C i C++. Procedury LAPACKa bazują na operacjach udostępnianych przez bibliotekę BLAS.
LAPACK stosuje podobną konwencję dotyczącą nazewnictwa, jak BLAS, tzn. pierwsza litera nazwy podaje typ zmiennej i dokładność, dwie następne typ macierzy a trzy kolejne wykonywaną operację.
Podprogramy biblioteki LAPACK pozwalają m.in. wykonywać faktoryzację macierzy, rozwiązywać układy równań liniowych czy rozwiązywać zagadnienie własne.
PoniewaŜ rozwiązanie zagadnienia zwykle rozkłada się na etapy często bardzo wygodnie jest korzystać z tzw. driver routines wywołujących odpowiednie podprogramy.
Na przykład procedura DGESV rozwiązuje w podwójnej dokładności układ równań liniowych z ogólną macierzą wywołując w tym celu podprogram DGETRF przeprowadzający rozkład LU a następnie uŜywającą DGETRS do rozwiązania układu trójkątnego.
Podobnie DPOSV jest driverem dla rozwiązania układu liniowego z macierzą symetryczną, dodatnio określoną wywołującym DPOTRF dla przeprowadzenia rozkładu Cholesky’ego a następnie DPOTRS dla rozwiązania układu z macierzą trójkątną.
Procedura DSYEV oblicza wartości i wektory własne symetrycznej, rzeczywistej macierzy uŜywając m. in. DSYTRD do sprowadzenia jej przekształceniami ortogonalnymi do postaci trójdiagonalnej i DSTEQR do znalezienia wartości i wektorów własnych trójdiagonalnej macierzy symetrycznej metodą QR.
Oczywiście zamiast korzystać z driver routines moŜna uŜyć tzw. expert routines (pozwalających określić więcej parametrów) lub budować rozwiązanie z podprogramów rozwiązujących poszczególne etapy.
Kody źródłowe procedur LAPACKA są dostępne na stronie Netlib:
http://www.netlib.org/lapack/
jednak (podobnie jak w przypadku BLASa) najlepiej skorzystać ze skompilowanych bibliotek dostępnych pod dany system operacyjny. Wspomniane wyŜej optymalizowane biblioteki typu Intel MKL czy ACML zawierają takŜe procedury LAPACKa.
Procedury LAPACKA wykorzystywane są w popularnych pakietach obliczeniowych jak R, Matlab, Octave, Scilab, itd.
Oprócz uŜycia gotowych skompilowanych bibliotek algebry liniowej istnieje moŜliwość dopasowania pakietu w celu optymalnego wykorzystania moŜliwości konkretnego komputera, np. przez optymalizację wykorzystania pamięci cache. Do automatycznego wygenerowania zoptymalizowanego BLASa (i niektórych procedur LAPACKa) słuŜy ATLAS (Automatically Tuned Linear Algebra Software):
http://math-atlas.sourceforge.net/
W zastosowaniach algebry liniowej w chemii kwantowej pojawia się często potrzeba znalezienia kilku wybranych wartości i wektorów własnych bardzo duŜej, rzadkiej (tzn.
mającej mało niezerowych elementów) macierzy. Do takich celów tworzone są specjalne biblioteki, np. ARPACK
http://www.caam.rice.edu/software/ARPACK/
III. RÓWNANIA NIELINIOWE.
2.8 Metoda bisekcji. Metoda Netwona. Metoda siecznych.
Zajmiemy się teraz znajdowaniem miejsc zerowych nieliniowej funkcji jednej zmiennej, czyli rozwiązywaniem równania:
f(x) = 0.
Pierwsza z metod to metoda bisekcji (połowienia przedziału). Opiera się na fakcie, iŜ jeśli funkcja f jest funkcją ciągłą w przedziale [a,b] i jej wartości na końcach przedziału mają róŜne znaki, to funkcja ta musi mieć zero w przedziale (a,b). Obliczamy zatem c= 2
(
a+b)
1
(środek przedziału) i sprawdzamy znak f(c). Wybieramy następnie ten z przedziałów [a,c]
oraz [c,b], na którego końcach wartości funkcji mają róŜne znaki. Otrzymujemy w ten sposób przedział dwa razy krótszy od poprzedniego, zawierający zero funkcji i całą procedurę moŜemy powtórzyć. Błąd, z jakim wyznaczamy miejsce zerowe funkcji określony jest przez długość przedziału. Proces kontynuujemy do momentu, aŜ długość przedziału zmaleje poniŜej pewnej dostatecznie małej wielkości, albo do chwili, gdy wartość funkcji w środku przedziału będzie dostatecznie bliska zeru.
Oczywiście, jeśli przegapimy nieciągłość funkcji, to metoda moŜe znaleźć się w kłopotach:
Druga metoda – metoda Newtona wykorzystuje wzór Taylora. Niech r będzie szukanym zerem funkcji, a x jego przybliŜeniem. Jeśli istnieją odpowiednie pochodne funkcji f, to mamy
0 = f(r) = f(x + h) = f(x) + hf’(x) + ...
gdzie h = r – x.
Jeśli h jest małe (czyli jesteśmy blisko zera), to moŜemy pominąć człony z wyŜszymi potęgami h i dostajemy
f(x) + hf’(x) = 0, skąd
) ( '
) (
x f
x h=− f .
Zatem z wyjściowego przybliŜenia pierwiastka x dostajemy lepsze, równe )
( '
) (
x f
x x− f .
[ ] a b
W metodzie Newtona zakładamy zatem startowe przybliŜenie miejsca zerowego x0 a następnie rekurencyjnie obliczamy lepsze przybliŜenia:
) ( '
) (
1
n n n
n f x
x x f
x + = − .
Metoda Newtona przybliŜa badaną funkcję liniową. Graficznie oznacza to, Ŝe rysujemy styczną do wykresu funkcji f w punkcie bliskim zeru r i znajdujemy punkt, w którym przecina ona oś x (dający nam lepsze przybliŜenie zera funkcji).
Dla pewnych funkcji metoda będzie miała powaŜne problemy ze zbieŜnością:
Jeśli jednak metoda działa, to zbieŜność jest szybka. Okazuje się, Ŝe jeśli r jest zerem pojedynczym f i f” jest ciągła, to istnieje takie otoczenie punktu r, Ŝe startując z tego otoczenia metoda zbieŜna jest do r.
Metoda Newtona wymaga obliczania pochodnej funkcji. MoŜna spróbować zastąpić pochodną f’(x) ilorazem róŜnicowym:
1 1) ( ) ) (
( '
−
−
−
≈ −
n n
n n
n x x
x f x x f
f .
r xn+1 xn
f(xn) f(x)
f(x1) r
f(x0)
x2
x1 x0
f(x)
Otrzymujemy wtedy metodę siecznych:
) ( ) ) ( (
1 1 1
−
+ − −
− −
=
n n
n n n
n
n f x f x
x x x
f x
x .
Do startu obliczeń potrzeba dwu punktów początkowych; znalezienie kaŜdej kolejnej wartości xn+1 wymaga obliczenia jednej nowej wartości funkcji.
Graficzna interpretacja metody siecznych róŜni się od metody Newtona zastąpieniem stycznej sieczną.
2.9 Układy równań nieliniowych.
Rozwiązanie układu równań nieliniowych bazuje na zastosowaniu wzoru Taylora. Na przykład, dla układu dwu funkcji dwu zmiennych
f1(x1,x2) = 0 f2(x1,x2) = 0
zakładamy, Ŝe para (x1,x2) jest przybliŜonym rozwiązaniem. W celu obliczenia poprawek h1 i h2 uŜywamy rozwinięcia Taylora funkcji dwu zmiennych zatrzymując tylko człony liniowe.
Mamy stąd
2 1 2 1 1 1 2 1 1 2 2 1 1
1( , ) ( , )
0 x
h f x h f x x f h x h x
f ∂
+ ∂
∂ + ∂
≈ + +
=
2 2 2 1 2 1 2 1 2 2 2 1 1
2( , ) ( , )
0 x
h f x h f x x f h x h x
f ∂
+ ∂
∂ + ∂
≈ + +
=
Pochodne cząstkowe są tu obliczane w x1, x2. Otrzymujemy stąd układ równań liniowych
−
=
∂
∂
∂
∂
∂
∂
∂
∂
) , (
) , ( /
/
/ /
2 1 2
2 1 1 2
1 2 2 1 2
2 1 1 1
x x f
x x f h
h x f x f
x f x f
Jego macierz to jakobian J czyli macierz pochodnych cząstkowych.
Jeśli macierz J jest nieosobliwa, to moŜemy rozwiązać ten układ wyznaczając poprawki h1 i h2 dostając lepsze przybliŜenie rozwiązania (x1 + h1, x2 + h2).
f(xn)
xn-1
r xn+1 xn
f(xn-1) f(x)
Procedurę tę stosujemy iteracyjnie:
+
=
+ +
) ( 2 ) ( 1 ) ( 2 ) ( 1 )
1 ( 2
) 1 ( 1
k k k
k k
k
h h x
x x
x
wyznaczając poprawki jako rozwiązania układu liniowego
−
=
) , (
) , (
) ( 2 ) ( 1 2
) ( 2 ) ( 1 1 )
( 2 ) ( 1
k k
k k k
k
x x f
x x f h
J h .
Analogicznie rozwiązujemy większe układy równań. Dla układu fi(x1, x2, ..., xn) = 0 (1 ≤ i ≤ n)
mamy
+
=
+ + +
) (
) ( 2 ) ( 1
) (
) ( 2 ) ( 1
) 1 (
) 1 ( 2
) 1 ( 1
k n
k k
k n
k k
k n
k k
h h h
x x x
x x x
M M
M ,
przy czym poprawki wyznaczamy rozwiązując układ n równań liniowych
−
=
) , , , (
) , , , (
) , , , (
) ( ) ( 2 ) ( 1
) ( ) ( 2 ) ( 1 2
) ( ) ( 2 ) ( 1 1
) (
) ( 2 ) ( 1
k n k
k n
k n k
k
k n k
k
k n
k k
x x
x f
x x
x f
x x
x f
h h h J
L M
L L
M
z macierzą jakobianu
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
=
n n n
n
n n
x f x
f x f
x f x
f x f
x f x
f x f J
/ /
/
/ /
/
/ /
/
2 1
2 2
2 1 2
1 2
1 1 1
L M M
M M
L L
.
Pochodne cząstkowe są obliczane w (x1(k),x(2k),K,xn(k)).
