Wydział Podstawowych Problemów Techniki
Kierunek: Fizyka
Specjalność: Fizyka komputerowa
PRACA DYPLOMOWA
Środowisko programowe
do wyznaczania macierzy odwrotnej
do symetrycznej macierzy trójdiagonalnej
Praca magisterska
Aleksander Pasek
Promotor: dr hab. inż. Włodzimierz Salejda, prof. PWr
słowa kluczowe:
symetryczna macierz trójdiagonalna, macierz odwrotna,
programowanie, numeryczna algebra liniowa
krótkie streszczenie:
Wynikiem pracy jest program komputerowy do wyznaczania macierzy
odwrotnej do symetrycznej macierzy trójdiagonalnej. Służy on do
wykonywania obliczeń oraz badania dokładności i czasów wykonania
algorytmów numerycznych. Opracowano i zbadano 22 metody numeryczne,
z czego 6 to metody nowe, które oparte są o wzory analityczne.
Wrocław, marzec 2007
Pragnę złożyć serdeczne podziękowania opiekunowi naukowemu,
dr hab. inż. Włodzimierzowi Salejdzie, prof. PWr
za zainteresowanie tematem, cierpliwość i pomoc
przy realizacji pracy dyplomowej.
Spis treści
1. Wprowadzenie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2. Przegląd metod numerycznych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1. Metody bezpośrednie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1. Metoda wyznacznikowa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2. Metoda eliminacji Gaussa, rozkład LU . . . . . . . . . . . . . . . . . . . 7
2.1.3. Rozkład Cholesky ego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.4. Rozkład QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.5. Rozkład SVD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2. Metody iteracyjne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1. Metoda Jacobiego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.2. Metoda Richardsona . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.3. Metoda nadrelaksacji (SOR) . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.4. Metoda Czebyszewa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.5. Metoda gradientów sprzężonych (CG) . . . . . . . . . . . . . . . . . . . . 16
2.3. Metody oparte o wzory analityczne . . . . . . . . . . . . . . . . . . . . . . . . . 17
Metoda dla macierzy symetrycznej . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Metoda dla macierzy diagonalnie dominującej . . . . . . . . . . . . . . . . . . . . 20
3. Opis środowiska obliczeniowego (Program TIPS) . . . . . . . . . . . . . . . . . 22
3.1. y
ródła danych i algorytmów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2. Ocena dokładności obliczeń . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3. Interfejs użytkownika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1. Informacje podstawowe . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2. Zapis i odczyt opcji obliczeń . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.3. Zakładka  Macierze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.4. Zakładka  Algorytmy i wyniki . . . . . . . . . . . . . . . . . . . . . . . 31
3.3.5. Zakładka  Testy i wykresy . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3.6. Okno  Wykres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.7. Okno  Obliczenia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3.8. Korzystanie ze wzorów . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4. Przeprowadzone testy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.1. Porównanie czasów wykonania algorytmów . . . . . . . . . . . . . . . . . . . . . 45
4.2. Porównanie dokładności - macierz Poissona . . . . . . . . . . . . . . . . . . . . . 46
4.3. Porównanie dokładności - macierz Hilberta . . . . . . . . . . . . . . . . . . . . . 48
4.4. Porównanie dokładności - macierz losowa . . . . . . . . . . . . . . . . . . . . . . 52
5. Podsumowanie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
Spis ilustracji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Literatura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2
1. Wprowadzenie
Wyznaczanie macierzy odwrotnej jest jednym z podstawowych problemów algebry
numerycznej. Do tego celu stosuje się różne algorytmy numeryczne. Macierz odwrotna
przydaje się w problemach automatyki sterowania i przetwarzania sygnałów, gdyż po-
zwala szybko rozwiązywać układy równań dla zmieniających się wektorów wejściowych
podczas gdy macierz układu równań pozostaje stała. Wiele zjawisk opisuje się przez mo-
dele liniowe, a matematyczną reprezentacją operatorów liniowych są właśnie macierze.
Każdą macierz symetryczną można zredukować do symetrycznej macierzy trójdiagonal-
nej poprzez przedstawienie zagadnienia w innej bazie przestrzeni liniowej. Konieczność
wyznaczania macierzy odwrotnej do symetrycznej macierzy trójdiagonalnej pojawia się
zatem w wielu zagadnieniach fizyki, ekonomii oraz inżynierii.
W praktyce obliczeń numerycznych zwykle unika się bezpośredniego obliczania ma-
cierzy odwrotnej, zastępując je przez zagadnienie rozwiązywania układu równań, które
jest lepiej uwarunkowane [1, 2]. Istnieją jednak sytuacje, w których dokładna znajomość
macierzy odwrotnej jest bardzo przydatna [3, 4]. Na uwagę zasługuje na przykład zagad-
nienie rozwiązywania czasowego równania Schr�dingera [5], gdzie dla każdego kolejnego
kroku czasu należy rozwiązywać trójdiagonalny układ równań.
Znane jest wiele klasycznych metod numerycznych stosowanych do wyznaczania ma-
cierzy odwrotnej. Są to zarówno metody bezpośrednie, jak i metody iteracyjne. Zakoń-
czone sukcesem próby uzyskania rozwiązania analitycznego dla problemu odwracania
macierzy trójdiagonalnej pozwalają na sformułowanie metody numerycznej nowego ty-
pu. Metodę zwaną w tej pracy  metodą analityczną można zastosować do uzyskania
analitycznego rozwiązania zagadnienia odwracania symetrycznej macierzy trójdiagonal-
nej dla macierzy, których elementy są ustalone [6, 7]. Dla dowolnych elementów macierzy
ta sama metoda pozwala sformułować algorytm inny niż dotychczas stosowane dla tego
typu zagadnień.
Celem tej pracy jest opracowanie i napisanie programu komputerowego, który po-
zwoli zbadać własności numeryczne tej  metody analitycznej na tle własności innych
wybranych metod bezpośrednich i iteracyjnych.
W rozdziale 2 przedstawiono krótko każdą z metod numerycznych zaimplementowa-
nych w programie. Podrozdział 2.1 zawiera opis pięciu podstawowych metod bezpośred-
nich, czyli takich, które dają wynik po skończonej liczbie kroków algorytmu. W podroz-
dziale 2.2 natomiast znajduje się opis zastosowanych metod iteracyjnych. Opis metod
 analitycznych , wraz z ich częściowym wyprowadzeniem na podstawie prac [7, 8] za-
wiera podrozdział 2.3. Zdefiniowano w nim 6 wersji algorytmów numerycznych opartych
o te metody.
3
Rozdział 3 poświęcono na dokładniejszy opis środowiska obliczeniowego (zestawu
programów komputerowych). Opis ten podzielono na trzy części. W części pierwszej
(podrozdział 3.1) zawarty jest wykaz 22 programów obliczeniowych realizujących odwra-
canie symetrycznej macierzy trójdiagonalnej. Podrozdział 3.2 przedstawia zastosowane
w programie metody badania i porównywania dokładności obliczeń. Część trzecia, za-
warta w podrozdziale 3.3, to opis możliwości programu, wyglądu poszczególnych okien
i sposobu korzystania z dostępnych funkcji.
Przeprowadzono także szereg testów numerycznych, a ich wyniki przedstawione są
w rozdziale 4. Znajduje się tam porównanie czasów wykonania poszczególnych algoryt-
mów oraz porównanie dokładności wyników dla trzech przykładowych macierzy: Poisso-
na, Hilberta oraz macierzy losowej.
W podsumowaniu, znajdującym się w ostatnim rozdziale 5, zawarto wnioski z prze-
prowadzonych testów, ze szczególnym uwzględnieniem wniosków dotyczących metod
 analitycznych . Pracę zamyka spis ilustracji oraz lista pozycji cytowanej literatury.
Do pracy dołączono też płytę CD zawierającą oprogramowanie w wersji wykonywalnej
i z
ródłowej wraz z programem instalacyjnym.
2. Przegląd metod numerycznych
Rzeczywista symetryczna macierz trójdiagonalna stopnia N jest następującej postaci:
�ł �ł
d1 e1
�ł �ł
e1 d2 e2
�ł �ł
�ł �ł
. .
�ł �ł
. .
e2 . .
�ł �ł
TN = �ł �ł . (2.1)
. .
�ł �ł
. .
. .
�ł eN-2 �ł
�ł �ł
�ł
eN-2 dN-1 eN-1 łł
eN-1 dN
Macierz TN zawiera zawsze liczbę zer równą N2 - 3N + 2. Ze względu na symetrię
wystarczy 2N - 1 liczb aby zapisać macierz TN w pamięci operacyjnej komputera.
Macierz do niej odwrotna T-1 nie jest już w ogólności trójdiagonalna, czyli ma N2 nie-
N
zerowych elementów, a ze względu na symetrię można ją zapisać w pamięci przy pomocy
(N2 + N)/2 liczb. Algorytmy obliczające macierz odwrotną do macierzy trójdiagonalnej
powinny więc mieć złożoność obliczeniową nie mniejszą niż O(N2).
W kolejnych podrozdziałach przedstawione są wybrane metody bezpośrednie i itera-
cyjne stosowane do obliczania macierzy odwrotnej do symetrycznej macierzy trójdiago-
nalnej. Podrozdział 2.3 opisuje metodę analityczną, która pozwala wyprowadzać wzory
na poszczególne elementy macierzy T-1. Z metody tej wynika bezpośrednio algorytm
N
numeryczny, który został także zaimplementowany w programie TIPS. Nazwa progra-
mu TIPS to skrótowiec składający się z pierwszych liter nazwy angielskiej: Tridiagonal
Inverse Problem Solver.
2.1. Metody bezpośrednie
2.1.1. Metoda wyznacznikowa
Elementy macierzy odwrotnej do macierzy kwadratowej A są postaci [9, 10]

Mji
A-1 = (-1)i+j , (2.2)
ij
|A|
gdzie |A| jest wyznacznikiem macierzy A, a Mij jest minorem (N - 1) rzędu, czyli
wyznacznikiem macierzy A po skreśleniu z niej i-tego wiersza i j-tej kolumny.
5
W celu obliczenia wyznacznika macierzy (2.1), tworzymy ciąg wielomianów {pi(z)},
tzw. ciąg Sturma [11]
ńł
�ł
p0(z) = 1
�ł
p1(z) = d1 (2.3)
�ł
ół
pi(z) = (di - z)pi-1(z) - e2 pi-2(z), i = 2, ..., N.
i-1
Po (N - 1) krokach rekurencyjnych otrzymujemy wielomian pN (z), który jest wie-
lomianem charakterystycznym macierzy TN , a szukany wyznacznik |TN | jest wyrazem
wolnym tego wielomianu, czyli
|TN | = pN (0). (2.4)
Następnie obliczymy minor Mij macierzy TN . Jeśli skreślimy z macierzy trójdiago-
nalnej (2.1) jeden wiersz i jedną kolumnę, to otrzymamy macierz blokowo-diagonalną
(przykład na rysunku 2.1).
Rysunek 2.1. Minor M42 macierzy T6
Wyznacznik takiej macierzy jest iloczynem wyznaczników kwadratowych bloków sto-
jących na przekątnej jeśli wszystkie elementy macierzy ponad (albo poniżej) tych bloków
są zerami [10]. Na tej samej zasadzie wyznacznik macierzy trójkątnej jest iloczynem jej
elementów diagonalnych. W przypadku z rysunku 2.1 możemy zapisać
M42 = d1� e2e3� (d5d6 - e2).
5
Jeśli analogicznie do (2.3) zdefiniujemy ciąg {qj(z)} jako
ńł
�ł
qj(z) = (dj - z)qj+1(z) - e2qj+2(z), j = N - 1, ..., 1
�ł
j
qN (z) = dN , (2.5)
�ł
ół
qN+1(z) = 1
to możemy zapisać
4-1

M42 = p2-1(0) ek q4+1(0).
k=2
6
Ogólnie minor Mij ma postać
ńł �ł �ł
j-i
�ł

�ł
�ł
�ł �ł
�ł pi-1(0) ej-kłł qj+1(0) dla i < j,
�ł
�ł
�łk=1 �ł
Mij = (2.6)
i-j
�ł

�ł
�ł
�ł �ł
�ł pj-1(0) ej+k-1łł qi+1(0) dla i > j.
�ł
ół
k=1
Podstawiając (2.6) i (2.4) do wzoru (2.2) otrzymujemy następujące wyrażenie na element
macierzy odwrotnej do symetrycznej macierzy trójdiagonalnej:
ńł �ł �ł
j-i
�ł

�ł
�ł
�ł �ł
pi-1(0) ej-kłł qj+1(0) dla i < j,
�ł
�ł
�ł
�ł
k=1

(-1)i+j �ł
T-1 = pi-1(0) qj+1(0) dla i = j, (2.7)
N
�ł �ł
ij �ł
pN (0)
�ł
i-j
�ł
�ł
�ł
�ł
�ł
�ł pj-1(0) ej+k-1łł qi+1(0) dla i > j.
�ł
ół
k=1

Pamiętamy przy tym, że nasza macierz jest symetryczna, więc teoretycznie T-1 =
N
ij

T-1 , ale wyniki obliczeń numerycznych dla (i < j) i (i > j) mogą być różne.
N
ji
Spodziewamy się, że także równość pN (0) = q1(0) będzie spełniona tylko ze skończoną
dokładnością. Podczas bezpośredniego obliczania wyznaczników może wystąpić przepeł-
nienie [ang. floating point overflow]. W pozycji [1] na stronie 57 znajdziemy opis jak
można się przed tym zabezpieczyć.
2.1.2. Metoda eliminacji Gaussa, rozkład LU
Metoda eliminacji Gaussa jest metodą rozwiązywania układów równań liniowych
typu
AX = B, (2.8)
gdzie macierz A jest macierzą układu, X macierzą rozwiązań, a B macierzą wektorów
kolumnowych składających się z wyrazów wolnych odpowiednich równań. W metodzie
tej przeprowadza się na macierzy A układu i na wektorze (wektorach) prawej strony
(kolumny macierzy B) operacje, które nie zmieniają rozwiązania układu X. Operacje te
wykonuje się w celu sukcesywnego zerowania kolejnych elementów macierzy A poniżej
(powyżej) jej głównej przekątnej aby otrzymać macierz trójkątną górną (dolną). Wyko-
nuje się te operacje w takiej kolejności, żeby zera wygenerowane w poprzednich krokach
pozostały zerami. Operacje, o których tutaj mowa, to zastępowanie wierszy (równań),
przez nie trywialne kombinacje liniowe innych wierszy (równań). Kiedy już powstanie
macierz trójkątna górna (dolna) okazuje się, że współczynniki kolejnych kombinacji li-
niowych układają się również w macierz trójkątną, ale dolną (górną).
7
Eliminacja Gaussa może być zatem wykorzystana do otrzymania tzw. rozkładu trój-
kątnego macierzy
A = LU, (2.9)
gdzie macierz L jest trójkątna dolna [ang. lower], a macierz U  trójkątna górna [ang.
upper]. Warunkiem istnienia rozkładu (2.9) jest niezerowość wszystkich tzw. minorów
wiodących macierzy A [3]. Dla symetrycznej macierzy trójdiagonalnej można ten waru-
nek zapisać przy użyciu notacji z (2.41) jako W1,j = 0, dla j = 2, ..., N.
Kiedy znamy rozkład LU, rozwiązanie X równania (2.8) otrzymujemy w procesach
podstawiania, czyli rozwiązywania układów równań z macierzami trójkątnymi
LY = B, (2.10)
UX = Y. (2.11)
Standardowa eliminacja Gaussa nie oblicza rozkładu LU. Wykonuje ona jednocze-
śnie te same operacje na macierzach A i B, nie musi więc zapamiętywać macierzy L,
a macierz B przekształca się automatycznie w macierz Y i podstawianie (2.10) nie jest
już potrzebne. Podstawianie (2.11) trzeba jednak na koniec wykonać.
Faktycznie w procesie eliminacji Gaussa zerujemy tylko elementy poniżej albo po-
wyżej przekątnej macierzy A. Nic nie stoi na przeszkodzie, aby zerować jednocześnie
elementy poniżej i powyżej tej przekątnej. Modyfikacja taka prowadzi do algorytmu eli-
minacji Gaussa-Jordana. Nie trzeba już wykonywać ani (2.10), ani (2.11), bo macierz
B przekształca się automatycznie w macierz rozwiązań X.
Jeśli macierz B jest macierzą jednostkową B = I, to wcale nie trzeba jej pamiętać,
a wynik X = A-1 można zapisywać kolumna po kolumnie zamiast macierzy A (w tej
samej pamięci). W takim przypadku mówimy, że algorytm  pracuje w miejscu (tzn.
nadpisuje pamięć danych wynikami).
Metoda eliminacji Gaussa byłaby dobrą metodą gdyby nie to, że jest niestabilna
numerycznie ze względu na możliwość wystąpienia zerowych, lub bardzo małych elemen-
tów na przekątnej macierzy A [1, 2, 9]. Problem ten rozwiązuje się zauważając prosty
fakt: przekątna nie musi być przekątną. Przestawiając odpowiednie wiersze macierzy
A przenosimy problematyczne elementy diagonalne gdzie indziej. Przestawienie wierszy
uzyskuje się mnożąc macierz lewostronnie przez tzw. macierz permutacji P. Rozkład LU
nie jest już więc rozkładem macierzy A, ale macierzy PA
PA = LU (2.12)
Aby otrzymać wynik X, obliczamy najpierw Y przez podstawienie
LY = PB (2.13)
oraz stosujemy znane już podstawienie (2.11). Algorytm ten nazywany jest elimina-
cją Gaussa z częściowym wyborem elementu głównego albo rozkładem PLU.
Oczywiście macierzy P nie znamy a priori, gdyż powstaje ona podczas wykonywania
rozkładu. Rozkład PLU jest wykonalny dla macierzy A odwracalnej (nieosobliwej) [3].
8
W rzadkich przypadkach, kiedy częściowy wybór elementu głównego nie jest wystar-
czający proponuje się algorytm eliminacji Gaussa z pełnym wyborem elementu
głównego. Polega on na tym, że wybieramy element Aij do zamiany z elementem dia-
gonalnym Akk nie tylko spośród elementów i j = k (wybór częściowy), ale spośród
elementów i k oraz j k. Wiąże się to z koniecznością przestawiania nie tylko wierszy,
ale także kolumn macierzy A. Przestawienie kolumn formalnie uzyskuje się mnożąc ma-
�
cierz prawostronnie przez inną macierz permutacji P. Rozkład LU jest teraz rozkładem
�
macierzy PAP
�
PAP = LU. (2.14)
Aby otrzymać wynik X wykonujemy (2.13), a następnie
UZ = Y, (2.15)
�
X = PZ. (2.16)
Widzimy na podstawie (2.16), że przestawianie kolumn macierzy A powoduje przesta-
wienie wierszy w uzyskanym przez nas wyniku. Rozkład LU z pełnym wyborem elementu
głównego istnieje dla każdej macierzy kwadratowej [3].
Kiedy macierz układu jest trójdiagonalna symetryczna A = T i nie wykonujemy
zamiany wierszy ani kolumn, to jej rozkład LU ma następującą postać [12]:
�ł �ł �ł �ł �ł �ł
d1 e1 1 u1 e1
�ł �ł �ł �ł �ł �ł
. . . . .
�ł �ł �ł �ł �ł �ł
. .
. . . . .
e1 . . l1 .
�ł �ł �ł �ł �ł �ł
�ł �ł = �ł �ł �ł �ł . (2.17)
. . . . .
�ł
. . .
. . . . .
�ł eN-1 �ł �ł . . �ł �ł eN-1 �ł
łł �ł łł �ł łł
eN-1 dN lN-1 1 uN
Współczynniki {li} i {ui} dane są wzorem rekurencyjnym

li = ei/ui
u1 = d1; i = 1, ..., N - 1. (2.18)
ui+1 = di+1 - liei
Ciąg {ui} pojawia się często gdy rozpatrujemy macierz trójdiagonalną (patrz wzór
2.52). Macierz odwrotną obliczamy ze wzoru T-1 = U-1L-1. Można pokazać, że wyni-
kiem jest następujący wzór na elementy macierzy odwrotnej
�ł �ł �ł �ł
N k-1 k-1


1
�ł �ł
T-1 = (-1)i+j lpłł lqłł . (2.19)
ij
uk p=i q=j
k=max(i,j)
Widać stąd, że jeśli istnieje wyraz ciągu {ui} równy zero, to nie istnieje macierz odwrotna,
czyli macierz T jest osobliwa. Problemy numeryczne pojawiają się wtedy, gdy {ui} jest
bardzo małą liczbą, niekoniecznie dokładnie równą zero. Ten sam problem występuje
zawsze, gdy mamy do czynienia z ciągiem {ui}. Tutaj ratunkiem jest wybór elementu
głównego. Przy częściowym wyborze mamy do dyspozycji tylko jeden wiersz ((k+1)-szy),
który możemy zamienić z k-tym wierszem. Jeśli to zrobimy, macierz U będzie miała
9
niezerowe elementy także na drugiej naddiagonali. W ten sposób działa podprogram
DGTSV z biblioteki Lapack [13]. Pełnego wyboru elementu głównego zwykle się w tym
wypadku nie wykonuje, gdyż zamiana kolumn psuje trójdiagonalną strukturę macierzy,
a zapamiętanie nowej pozycji niezerowych elementów wymaga niepotrzebnie dodatkowej
pamięci. Jeśli wybór częściowy się nie powiedzie, to już wiemy, że macierz jest osobliwa
i kończymy obliczenia.
2.1.3. Rozkład Cholesky ego
Rozkład Cholesky ego jest podobny do rozkładu LU, z tą różnicą, że macierz górna
trójkątna jest transpozycją macierzy dolnej trójkątnej, czyli U = LT . Rozkład ten można
zastosować tylko do macierzy, które są dodatnio określone [ang. positive definite], tzn.
mają wszystkie wartości własne dodatnie [1, 2]. Macierze takie występują często w fizyce,
gdyż są one reprezentacjami operatorów hermitowskich o rzeczywistych wartościach wła-
snych (wartości te muszą być rzeczywiste, bo są wynikami pomiarów różnych wielkości
fizycznych).
Zapiszmy rozkład Cholesky ego LLT dla symetrycznej macierzy trójdiagonalnej:
�ł �ł �ł �ł �ł �ł
d1 e1 m1 m1 l1
�ł �ł �ł �ł �ł �ł
. . . . .
�ł �ł �ł �ł �ł �ł
. .
. . . . .
e1 . . l1 .
�ł �ł �ł �ł �ł �ł
=
�ł �ł �ł �ł �ł �ł .
. . . . .
�ł �ł �ł
. . .
. . . . .
�ł eN-1 �ł �ł . . łł �ł lN-1 �ł
łł �ł łł
eN-1 dN lN-1 mN mN
(2.20)
Bezpośrednio z (2.20) otrzymujemy związki rekurencyjne pozwalające obliczyć współ-
czynniki macierzy L:

li = ei/mi

m1 = d1; i = 1, ..., N - 1. (2.21)
2
mi+1 = di+1 - li
Rozkład ten nie wymaga przestawiania wierszy ani kolumn, a sam algorytm jest
poprawny numerycznie [1]. Możemy także uniknąć obliczania pierwiastka kwadratowego
� �T
i zamiast rozkładu LLT obliczyć rozkład LDL [1, 11], gdzie mi = 1 i D = diag(m2) =
�
i
diag(ui)
�ł
�
�ł
li = ei/ui
żł
u1 = d1; i = 1, ..., N - 1. (2.22)
e2
i
�ł
ui+1 = di+1 - �ł
ui
Ciągi {�} i {ui} z (2.22) to te same ciągi, co {li} i {ui} we wzorach (2.18). Rozkład
li
Cholesky ego jest zaimplementowany w podprogramie DPTSV biblioteki Lapack [13].
Zauważmy, że w (2.22) elementy ciągu {ui} mogą być ujemne, co sugerowałoby, że roz-
kład ten można wykonać także dla macierzy, które nie są dodatnio określone. Można, ale
algorytm nie jest wtedy numerycznie poprawny, co oznacza, że błąd może być dowolnie
duży nawet dla macierzy o małym wskaz
niku uwarunkowania [1].
10
Obliczamy X, macierz odwrotną do symetrycznej macierzy trójdiagonalnej przy po-
mocy rozkładu Cholesky ego przez rozwiązanie trójkątnych układów równań LZ = I
� �
oraz LT X = Z. Jeśli natomiast obliczymy rozkład T = LDLT , to macierz odwrotną
� �
otrzymamy rozwiązując równania macierzowe LDZ = I oraz LT X = Z.
2.1.4. Rozkład QR
Rozkład T = QR, zwany rozkładem ortogonalno-trójkątnym jest zawsze wykonalny
dla trójdiagonalnej macierzy symetrycznej [2, 11]. Macierz górna trójkątna R powstaje
w procesie tzw. dekompozycji, czyli po przeprowadzeniu (N -1) transformacji ortogonal-
nych macierzy T. Każda z tych transformacji na wyzerować jeden element poddiagonalny
macierzy T.
Transformacje, o których tutaj mowa to tzw. obroty Givensa. Mogą to być też inne
transformacje (np. odbicia Householdera albo zmodyfikowana ortogonalizacja Gramma-
-Schmidta), ale po ich zastosowaniu macierz R jest na ogół pełną macierzą trójkątną,
podczas gdy (N - 1) obrotów Givensa prowadzi do macierzy trójkątnej z tylko dwie-
ma naddiagonalami (2.25) (podobnie jak rozkład LU z częściowym wyborem elementu
głównego). Pojedynczy i-ty obrót Givensa ma postać macierzy
�ł �ł
1
�ł �ł
.
.
�ł �ł
.
�ł �ł
�ł �ł
ci si
�ł �ł
GT = �ł �ł . (2.23)
i
�ł -si ci �ł
�ł �ł
�ł �ł
.
.
�ł . łł
1
Kąt obrotu �i jest tak dobrany, aby wyzerować element poddiagonalny ei, co prowadzi
do związków
�i ei
cos �i = ci = , sin �i = si = , (2.24)
2 2
�i + e2 �i + e2
i i
gdzie �i to element di już zmieniony przez poprzedni obrót GT .
i-1
Proces dekompozycji GT � � � GT T = QT T = R prowadzi do macierzy trójkątnej
N-1 1
�ł �ł
r1 q1 v1
�ł �ł
.
�ł �ł
.
r2 q2 .
�ł �ł
�ł �ł
. . .
�ł �ł
. . .
. . .
�ł �ł
R = . (2.25)
�ł �ł
.
�ł
.
.
�ł qN-2 vN-2 �ł
�ł
�ł �ł
�ł
rN-1 qN-1 łł
rN
Macierz odwrotną otrzymujemy rozwiązując równanie macierzowe RX = QT względem
X, gdzie QT jest macierzą skumulowanych (tzn. pomnożonych przez siebie) transformacji
ortogonalnych.
11
Z postaci (2.23) wnioskujemy, że macierz QT jest prawie trójkątna dolna (tylko pierw-
sza naddiagonala jest niezerowa i zawiera wartości si). Jeśli w algorytmie dekompozycji
przedstawionym w [11] na stronie 64 zamienimy kolejność obliczania wielkości pi+1 i qi,
to w każdym kroku możemy wykonać o jedno mnożenie mniej. Po tych modyfikacjach
algorytm rozkładu QR symetrycznej macierzy trójdiagonalnej wraz z akumulacją trans-
formacji ortogonalnych przedstawia się następująco:
QT = I
{qi = vi = ei} i = 1, ..., N - 1
p = d1
ńł �ł

�ł 2 �ł
�ł �ł
ri = p2 + vi , c = p/ri, s = vi/ri
�ł �ł
�ł �ł

�ł �ł
�ł �ł
�ł �ł
QT = -s QT
�ł �ł
i+1,j i,j
�ł �ł
(2.26)
�ł j = 1, ..., i żł
QT = c QT
i,j i,j
i = 1, ..., N - 1
�ł �ł
�ł �ł
QT = s, QT = c
�ł �ł
i,i+1 i+1,i+1
�ł �ł
�ł �ł
�ł �ł
�ł p = cdi+1 - sqi, qi = cqi + sdi+1 �ł
�ł �ł
�ł �ł
ół
if i N - 2 then vi = sqi+1, qi+1 = cqi+1 �ł
rN = p.
2.1.5. Rozkład SVD
Rozkład SVD (Singular Value Decomposition) polega na zapisaniu dowolnej prosto-
kątnej macierzy rzeczywistej A w postaci iloczynu
A = UŁVT , (2.27)
gdzie macierze U i V są ortogonalne, a macierz Ł jest diagonalna. Niezerowe elementy
macierzy Ł to tzw. wartości szczególne macierzy A. Rozkład ten łatwo zastosować do
obliczenia macierzy odwrotnej do macierzy A, gdyż
A-1 = VŁ-1UT . (2.28)
Jeśli macierz A nie jest kwadratowa, to macierz (2.28) nazywamy macierzą pseudoodw-
rotną (znajduje ona zastosowanie w problemach minimalizacji funkcji metodą najmniej-
szych kwadratów).
Nowoczesne algorytmy służące do obliczania rozkładu SVD są dwuetapowe. Najpierw
przy pomocy transformacji ortogonalnych sprowadza się macierz A do postaci górnej
macierzy bidiagonalnej, a następnie oblicza się wartości i wektory szczególne tej macierzy.
W przypadku symetrycznej macierzy trójdiagonalnej zastosowanie tej metody napo-
tyka na trudności związane z koniecznością zapamiętania pełnych macierzy zawierają-
cych wektory szczególne oraz pomnożenia tych macierzy przez siebie. Z drugiej jednak
strony redukcja macierzy trójdiagonalnej do bidiagonalnej jest bardzo prosta i szybka
przy pomocy obrotów Givensa, a obliczenie macierzy odwrotnej do macierzy diagonalnej
Ł sprowadza się do obliczenia odwrotności jej elementów diagonalnych.
12
Proces redukcji dowolnej macierzy pasmowej do postaci górnej macierzy bidiagonal-
nej wykonuje procedura DGBBRD z biblioteki Lapack [13], a obliczenie rozkładu SVD
możemy wykonać jedną z dwóch metod (również z biblioteki Lapack):
" DBDSQR przeprowadza iterację metodą QR,
" DBDSDC stosuje algorytm typu  dziel i rządz
 .
Druga z podanych tu metod wymaga więcej pamięci i powinna działać szybciej dla
dużych macierzy. Powodem, dla którego należy zwrócić uwagę na metodę SVD jest jej
duża niezawodność [2]. Potrafi ona dawać bardzo dokładne wyniki nawet w przypadku
macierzy o wyznaczniku bliskim zeru (prawie osobliwych). Jeśli macierz będzie osobli-
wa, to pojawią się wartości szczególne równe zero, jednak liczba niezerowych wartości
szczególnych wciąż jest równa rzędowi macierzy wyjściowej.
2.2. Metody iteracyjne
Oprócz metod bezpośrednich (dających rozwiązanie po skończonej liczbie kroków)
istnieje cała rodzina metod i algorytmów iteracyjnych służących do rozwiązywania ukła-
dów równań liniowych. Iteracja polega na cyklicznym wykonywaniu tych samych operacji
w każdym kroku tak, aby wartości obliczone w poszczególnych krokach były zbieżne do
poszukiwanego rozwiązania. Aby rozpocząć iterację potrzebna jest znajomość tzw. wek-
tora początkowego x(0), który możemy wybrać dowolnie albo jako wcześniej obliczone
inną metodą rozwiązanie przybliżone (uwarunkowanie wstępne [ang. preconditioning]).
Iteracja w i-tym kroku ma postać mnożenia niezmienionej macierzy A przez wektor
i dodania innego wektora [1, 2, 3]:
x(i+1) = Mx(i) + w. (2.29)
Główny koszt wykonania metody iteracyjnej związany jest z mnożeniem macierzy
przez wektor, które jest czasochłonne. Jeśli jednak macierz jest macierzą rzadką, tzn.
liczba jej niezerowych elementów jest rzędu O(N), to mnożenie takie można wykonać
znacznie szybciej. Dzięki temu każda metoda iteracyjna potrafi skorzystać z faktu rzad-
kości macierzy w przeciwieństwie do metod bezpośrednich, które trzeba specjalnie za-
projektować, aby mogły skorzystać z tego samego faktu.
Liczba pojedynczych operacji zmiennoprzecinkowych, które trzeba wykonać w meto-
dzie iteracyjnej jest niekiedy znacznie większa niż w metodzie bezpośredniej zastosowanej
dla macierzy tych samych rozmiarów. Większe niebezpieczeństwo kumulacji pojedyn-
czych błędów zaokrągleń nie martwi nas już jednak tak bardzo jak w przypadku metod
bezpośrednich, ponieważ jeśli tylko metoda iteracyjna jest zbieżna, błędy wygenerowane
w jednym kroku mogą zostać poprawione w kolejnych krokach. Dzieje się tak dlatego, że
na każdym etapie iteracji istnieje stały i niezmienny czynnik: macierz iteracji M. Podczas
gdy algorytmy bezpośrednie manipulują wyrazami macierzy A, algorytmy iteracyjne raz
obliczoną na początku macierz M pozostawiają w niezmienionej postaci aż do końca obli-
czeń. Tak działają stacjonarne metody iteracyjne, ale istnieją też metody niestacjonarne,
w których zarówno macierz M, jak i wektor w mogą zmieniać się w kolejnych iteracjach
(zmiany te są zwykle niewielkie i mają na celu przyspieszenie zbieżności metody).
13
Warunkiem koniecznym zbieżności metody iteracyjnej postaci (2.29) jest prosta za-
leżność M < 1, gdzie � jest dowolną normą macierzy zgodną z normą wektorową.
Warunek ten jest konieczny, ale nie zawsze wystarczający [1].
Wystarczającym warunkiem jest aby promień spektralny macierzy iteracji M (naj-
większa co do modułu wartość własna M), spełniał [1, 3, 14]
�(M) = max |i| < 1. (2.30)
1 i N
Należy też zwrócić uwagę na fakt, że im mniejszy promień spektralny macierzy iteracji
M, tym szybsza zbieżność metody iteracyjnej [1, 3].
Aby zastosować metodę iteracyjną rozwiązywania N równań liniowych postaci Ax =
b do odwracania macierzy A, będziemy uruchamiać daną metodę N razy obliczając
kolejne kolumny macierzy odwrotnej przez podstawienie wektora ei = (0, ..., 0, ei =
1, 0, ..., 0)T za wektor b. Wyjątkiem jest metoda Richardsona z [14], zaimplementowana
w taki sposób, że przeprowadza jednocześnie iterację wszystkich (nawet więcej niż N)
kolumn macierzy B rozwiązując de facto równanie macierzowe AX = B. Przez podsta-
wienie B = I otrzymamy X H" A-1.
2.2.1. Metoda Jacobiego
Najprostszą z metod iteracyjnych jest metoda Jacobiego. Polega ona na rozkładzie
macierzy A na sumę macierzy wg wzoru
A = L + D + U, (2.31)
gdzie D jest macierzą diagonalną i zawiera elementy diagonalne macierzy A, a macierze
L i U są odpowiednio dolną i górną częścią macierzy A z zerami na głównej przekątnej.
Po podstawieniu (2.31) do układu równań Ax = b proste przekształcenia prowadzą
do wzoru [1, 14]
x = -D-1(L + U)x + D-1b. (2.32)
Porównując powyższy wzór z (2.29) mamy M = -D-1(L + U) oraz w = D-1b.
Zbieżność metody Jacobiego wymaga, aby macierz M była nieredukowalna i diagonalnie
słabo dominująca [1, 3]. Taka macierz jest też nieosobliwa.
2.2.2. Metoda Richardsona
Metoda ta również prowadzi do stałej macierzy iteracji M i stałego wektora w,
a oparta jest na następujących przekształceniach [1, 3, 14]:
Ax = b,
ąAx = ąb,
x + ąAx = x + ąb,
x = x - ąAx + ąb,
x = (I - ąA)x + ąb,
M = I - ąA,
w = ąb.
14
A
atwo pokazać [1], że jeśli macierz A jest symetryczna, to aby promień spektralny
macierzy M był mniejszy od jedności, wszystkie wartości własne macierzy A muszą być
tego samego znaku. Warunek ten spełnia symetryczna macierz dodatnio określona.
Parametr ą jest liczbą całkowitą dobraną tak, aby promień spektralny macierzy
M = I - ąA był jak najmniejszy. Najmniejszą wartość promienia spektralnego równą


max - min

�min(M) = (2.33)

max + min)
otrzymuje się dla
2
ą = . (2.34)
max + min
Żeby móc skorzystać z tej metody, trzeba w naszym przypadku znać ekstremalne war-
tości własne symetrycznej macierzy trójdiagonalnej T. Wiemy (np. z [9]), że wszystkie
wartości własne macierzy dodatnio określonej zawarte są w przedziale (0, T "), gdzie
� " to norma nieskończona macierzy, czyli maksymalna suma modułów w wierszu
(patrz (3.3)). Jeśli nie chcemy szukać wartości własnych możemy spróbować zastąpić
min przez 0, a max przez T ". Kiedy zbieżność metody wciąż nie jest zadowalają-
ca, możemy spróbować wyznaczyć największą i najmniejszą wartość własną macierzy
trójdiagonalnej przy pomocy stabilnej metody Martina-Deana opisanej w [11]. Okazuje
się bowiem, że liczba ujemnych wyrazów w ciągu (2.52) równa jest liczbie ujemnych
wartości własnych macierzy T. Wystarczy przesuwając widmo macierzy znalez
ć takie
przesunięcie, w otoczeniu którego liczba ta  przeskakuje o jeden i metoda bisekcji po-
zwoli znalez
ć przybliżenie jednej wybranej wartości własnej macierzy T (w tym wypadku
dwóch: maksymalnej i minimalnej).
2.2.3. Metoda nadrelaksacji (SOR)
Metoda ta korzysta z obu powyższych pomysłów i jest jakby połączeniem metody
Jacobiego z metodą Richardsona. W metodzie Jacobiego pomysłem był rozkład macierzy
na sumę (2.31) i przeniesienie jednego z powstałych w ten sposób czynników na drugą
stronę równania, a w metodzie Richardsona pomnożenie obu stron równania przez pewien
parametr i dodanie do obu stron składnika zawierającego wektor x. W metodzie SOR
stosujemy oba te przekształcenia [1, 3, 14, 15]:
Ax = b,
�(L + D + U)x = �b,
�Lx = -�Dx - �Ux + �b,
Dx + �Lx = Dx - �Dx - �Ux + �b,
(D + �L)x = ((1 - �)D - �U)x + �b,
x = (D + �L)-1((1 - �)D - �U)x + �(D + �L)-1b,
M = (D + �L)-1((1 - �)D - �U),
w = �(D + �L)-1b.
Właściwy wybór wartości parametru � jest trudny, ale niezwykle istotny, bo może
znacznie przyspieszyć zbieżność algorytmu [1, 15]. Wartość � musi spełniać 0 < � < 2,
15
przy czym dla � < 1 mówimy o podrelaksacji, a dla � > 1 o nadrelaksacji. Kiedy
parametr � jest równy jeden, to metoda SOR sprowadza się do metody Gaussa-Seidela,
która jest zbieżna dla tych samych warunków, co metoda Jacobiego (tylko jest zwykle
zbieżna szybciej).
Dla symetrycznej i dodatnio określonej macierzy trójdiagonalnej można wykazać, że
optymalna wartość parametru � wyraża się przez
2

� = , (2.35)
1 + 1 - 2
max
gdzie max jest największą wartością własną macierzy M z iteracji Jacobiego, czyli M =
-D-1(L + U).
2.2.4. Metoda Czebyszewa
Metoda Czebyszewa należy do niestacjonarnych metod iteracyjnych, czyli w każdym
kroku iteracji, jej parametry zależą od poprzednich kroków. Wykorzystuje ona relacje
rekurencyjne wielomianów Czebyszewa, aby uzależnić kolejne przybliżenie od dwóch po-
przednich przybliżeń obliczonych w dwóch poprzednich krokach iteracji. Iteracja w me-
todzie Czebyszewa ma postać [1]
r(i) = Ax(i), (2.36)
x(i+1) = x(i) + (p(i-1)(x(i) - x(i-1)) - r(i))/q(i). (2.37)
Współczynniki p(i) i q(i) oblicza się ze wzorów [1]
"
2
a + b b - a a + b + 2 ab
c = , d = , q" = ,
2 4 4
2d
p(-1) = 0, q(0) = c, p(0) = ,
c
"

a + b ab 2 ab d
q(1) = + , ł(1) = ,
4 a + b a + b q"
d p(i-1)ł(i-1)
p(i-1) = , ł(i) = , q(i) = q" - ł(i). (2.38)
q(i-1) q"
Parametry a i b to odpowiednio minimalna i maksymalna wartość własna macierzy
T, albo (jeśli nie są one znane) przyjmuje się, że a = 0, a b = T . Głównym kosz-
tem metody Czebyszewa jest mnożenie macierzy przez wektor przy obliczaniu wektora
residualnego r.
2.2.5. Metoda gradientów sprzężonych (CG)
Metoda gradientów sprzężonych, podobnie jak metoda Czebyszewa jest niestacjonar-
ną metodą iteracyjną. Można ją zastosować do symetrycznej macierzy dodatnio określo-
nej. Oparta jest na minimalizacji funkcji
1
f(x) = x � A � x - b � x (2.39)
2
16
poprzez doprowadzanie jej gradientu "f = A � x - b do zera. Minimalizacja tej funkcji
polega na poszukiwaniu tzw. kierunku najszybszego spadku i na sukcesywnym poru-
szaniu się w tym kierunku. W każdym kroku poszukuje się takiej wartości ąi, która
minimalizuje wyrażenie f(x(i+1)), gdzie x(i+1) = x(i) + ąip(i).
Algorytm dla metody gradientów sprzężonych wykonuje w każdym kroku iteracji
jedno mnożenie macierzy A przez wektor i dwa iloczyny skalarne wektorów. Ta pro-
stota jest przyczyną, dla której metoda CG jest taka popularna. Wystarczy pamiętać
trzy wektory: przybliżone rozwiązanie x, resztę r i kierunek poszukiwań p. A
atwo też
udowodnić zbieżność iteracji CG zakładając, że obliczenia wykonywane są z nieskoń-
czoną dokładnością [2]. Interpretacja wyników w arytmetyce zmiennoprzecinkowej jest
już trudniejsza i skuteczność algorytmu CG w tym przypadku jest zadziwiająco duża.
W teorii [2] jeśli macierz A ma N różnych wartości własnych, to iteracja CG (przy
nieskończonej precyzji obliczeń) osiąga minimum w co najwyżej N krokach. Mamy więc
metodę iteracyjną o cechach metody bezpośredniej, dla której zbieżność jest teoretycznie
pewna po skończonej liczbie iteracji! Jeśli mnożenie wektora przez macierz można wyko-
nać przy pomocy liczby operacji rzędu O(N), to sumaryczny koszt metody gradientów
sprzężonych powinien być rzędu O(N2).
Podstawową wersję algorytmu, którą opracowali Hestens i Stiefel, przytoczymy tutaj
za [2]. Tej samej wersji używa wykorzystana w programie TIPS procedura z [15]. Dla
wartości początkowych x0 = 0, r0 = b i p0 = r0, n-ty krok iteracji ma postać
rT rn-1
n-1
ąn = nowa wartość parametru ą
pT Apn-1
n-1
xn = xn-1 + ąnpn-1 kolejne przybliżenie
rn = rn-1 + ąnApn-1 reszta (2.40)
rT rn
n
�n = nowa wartość parametru �
rT rn-1
n-1
pn = rn + �npn-1 nowy kierunek poszukiwań
2.3. Metody oparte o wzory analityczne
Metoda dla macierzy symetrycznej
Wróćmy do metody wyznacznikowej z podrozdziału 2.1.1. Obie rekurencje (2.3)
i (2.5) obliczają wyznacznik macierzy TN , przy pomocy wyznaczników niższych stopni
(nazwijmy je minorami diagonalnymi) stosując rozwinięcie Laplace a [10] odpowied-
nio względem ostatniego i pierwszego wiersza. Minory diagonalne, to takie wyznaczniki,
których przekątna jest w całości zawarta w przekątnej macierzy TN .
17
Oznaczymy minory diagonalne rzędu j - i + 1 jako


di ei


ei di+1 ei+1


. . .

. . .
Wi,j = . (2.41)
. . .



ej-2 dj-1 ej-1


ej-1 dj
Wzór (2.7) dla i < j (symetria macierzy T-1) zapisujemy krócej jako
N
�ł �ł
j-1


(-1)i+j
T-1 = W1,i-1 �ł ekłł Wj+1,N . (2.42)
N
ij
W1,N
k=i
Stosunkowo łatwo wykazać indukcyjnie, że dla N 1 i j = 1, ..., N zachodzi
�ł �ł
j-1

�ł łł
e2 Wj+1,N = W1,j-1W2,N - W2,j-1W1,N , (2.43)
k
k=1
jeśli przyjmiemy, że dla i > j jest Wi,j = 2 - i + j. Gdy podstawimy Wj+1,N ze wzoru
(2.43) do (2.42) otrzymujemy

(-1)i+j W1,j-1W2,N - W W1,N
2,j-1
T-1 = W1,i-1 i-1 =
N
ij
W1,N
ek j-1 ek
k=1 k=1

W1,i-1W1,j-1 W2,N W2,j-1

= (-1)i+j - =
i-1
ek j-1 ek W1,N W1,j-1
k=1 k=1

(-1)i-1W1,i-1 (-1)j-1W1,j-1 (-1)jW2,j-1 (-1)N+1W2,N
= - . (2.44)
i-1
j-1
(-1)j-1W1,j-1 (-1)N W1,N
ek
ek
k=1
k=1
Definiujemy dwa nowe ciągi {pi} i {qj}:
(-1)iW1,i (-1)j+1W2,j
pi = , qj = . (2.45)
i
j
ek
ek
k=1
k=1
Nie są to te same ciągi, co (2.3) i (2.5), gdzie dla odróżnienia pisaliśmy zawsze argument
wielomianu z. Dzięki (2.45) wzór (2.44) sprowadza się do prostej postaci


qj-1 qN
T-1 = pi-1pj-1 - . (2.46)
N
ij
pj-1 pN
Wzór (2.46) pochodzi z pracy [8], jednak metoda jego wyprowadzenia przedstawiona
tutaj jest bardziej oczywista. Ciągi (2.45) są w [8] zdefiniowane rekurencyjnie
eipi = -ei-1pi-2 - dipi-1, p0 = 1, p1 = -d1/e1,
(2.47)
eiqi = -ei-1qi-2 - diqi-1, q0 = 0, q1 = 1/e1.
18
Korzyść, jaką uzyskujemy w stosunku do wzoru (2.7) polega na tym, że nie liczymy
już bezpośrednio wartości minorów, które mogą być bardzo duże i powodować problemy
numeryczne. Zamiast tego obliczamy rekurencyjnie według (2.47) elementy ciągów {pi}
i {qj}, które mają mniejszą wartość bezwzględną dzięki iloczynom w mianownikach wzo-
rów (2.45). Niestety te same iloczyny powodują, że wartości pN i qN potrzebne w (2.46)
nie są dobrze określone (nie istnieje w macierzy wyraz eN ). Dobrze określony jest nato-

miast ich iloraz qN /pN , ponieważ T-1 = -qN /pN . Wartości pN i qN nie otrzymamy
N
11

z ciągów (2.47), więc musimy obliczyć T-1 inną metodą i zamiast z (2.46) skorzystać
N
11
z następującego wzoru

T-1 = pi-1 qj-1 + pj-1 T-1 . (2.48)
N N
ij 11
Ciąg {qj(0)} (2.5) obliczany jest od końca (tzn. dla j = N -1, ..., 1) w przeciwieństwie
do ciągu {pi(0)} (2.3), który liczymy od początku (tzn. dla i = 2, ..., N). Natomiast ciągi
{pi} i {qj} (2.47) można obliczać jednocześnie, w jednej pętli, w miarę czytania kolejnych
wartości {di} i {ei}. Dzięki temu można ponownie wykorzystać pamięć przeznaczoną na
macierz TN .
Liniowe rekurencje drugiego rzędu (2.47) pozwalają wyznaczać kolejny wyraz każdego
z ciągów przy pomocy dwóch poprzednich wyrazów tego samego ciągu (metoda 1A). Tak
jednak być nie musi. W pracy [8] znajdujemy także łatwy do udowodnienia indukcyjnie
związek
ei(pi-1qi - piqi-1) = 1, dla i 1. (2.49)
Jeśli teraz najpierw obliczymy pi, to qi możemy wyznaczyć z (2.49) nie używając
jawnie wyrazu qi-2 (metoda 1B). Podobnie jeśli najpierw obliczymy qi, to pi możemy
wyznaczyć nie używając pi-2 (metoda 1C). Daje to dwie alternatywne w stosunku do
(2.47) metody obliczania ciągów {pi} i {qj}, które choć równoważne matematycznie,
mogą nie być równoważne numerycznie.
Do skorzystania z metod 1A, 1B i 1C opartych o wzór (2.48) musimy poznać wartość
pierwszego elementu diagonalnego macierzy odwrotnej. W pracy [8] podaje się odpo-
wiedni wzór w postaci ułamka łańcuchowego (ciągłego)

1
T-1 = . (2.50)
N
11
e2
1
d1 -
e2
2
d2 -
.
.
.
� � � - � � �
e2
N-2
e2
N-1
dN-1 -
dN
Obliczenie wzoru (2.50) sprowadza się do obliczenia odwrotności wyrazu v1 ciągu
(2.53).
19
Metoda dla macierzy diagonalnie dominującej
Wyprowadzona i udowodniona w pracy [7] metoda obliczania wyrazów macierzy
odwrotnej dotyczy macierzy trójdiagonalnej, która jest diagonalnie dominująca (tzn.
|di| > |ei-1| + |ei|), ale niekoniecznie symetryczna. Zastosowana jednak do macierzy
symetrycznej, po przyjęciu oznaczeń z (2.1), (2.3) i (2.5) ma następującą postać:
ńł

pi-1(0)
�ł
�ł
�ł -ei T-1 i < j
�ł
N
�ł
i+1,j
pi(0)
�ł
�ł
-1
�ł
�ł
pj-2(0) qj+2(0)
T-1 = . (2.51)
dj - e2 - e2 i = j
N j-1 j
ij �ł
pj-1(0) qj+1(0)
�ł
�ł
�ł

�ł
�ł qi+1(0)
�ł
�ł
-ei-1 T-1 i > j
ół
N
i-1,j
qi(0)
Ze wzoru (2.51) wynika, że jeśli najpierw obliczymy wyrazy diagonalne macierzy T-1,
N
to pozostałe wyrazy łatwo wyznaczyć, bo są one od siebie nawzajem zależne! Minory
z ciągów {pi(0)} i {qi(0)} wchodzą do wzoru (2.51) zawsze w postaci ilorazów dwóch
kolejnych minorów z tych samych ciągów.
Można zatem zdefiniować pomocnicze ciągi w następujący sposób:
pi(0) e2
i-1
u1 = d1, ui = = di - , i = 2, ..., N; (2.52)
pi-1(0) ui-1
e2
qj(0)
j
vN = dN , vj = = dj - , j = N - 1, ..., 1. (2.53)
qj+1(0) vj+1
Ciąg (2.52) jest znany i szeroko stosowany w zagadnieniach własnych z symetryczną
macierzą trójdiagonalną [11], ponieważ są to kolejne wyrazy diagonalne macierzy D
z rozkładu TN = LDLT (porównaj z (2.18) i (2.22)).
Dzięki (2.52) i (2.53) możemy uprościć (2.51) do postaci
ńł
ei
�ł
�ł - T-1 i < j
�ł
i+1,j
�ł ui N
�ł
-1
�ł
�ł
e2
j
T-1 = . (2.54)
uj - i = j
N
ij �ł
v
�ł
�ł

�ł
�ł ei-1 j+1
�ł
ół - T-1 i > j
i-1,j
vi N
Wzór (2.54) można obliczać na trzy sposoby. Metoda 2A polegać będzie na tym, że

najpierw obliczymy {vj}, a potem T-1 jednocześnie zapamiętując {ei/ui} w pamięci
N
ii
poprzednio zajętej przez {ei}. Wyrazy niediagonalne obliczymy ostatecznie ze wzoru
(2.54) dla i < j. Metoda 2B polegać będzie na tym, że najpierw obliczymy {ui}, a potem

T-1 jednocześnie zapamiętując {ej/vj+1} w pamięci poprzednio zajętej przez {ej}.
N
jj
Wyrazy niediagonalne obliczymy ostatecznie ze wzoru (2.54) dla i > j. Metoda 2C

polegać będzie na tym, że najpierw obliczymy {ui}, potem osobno {vj}, dalej T-1 ,
N
ii
20
a wyrazy niediagonalne macierzy odwrotnej obliczymy korzystając z obu wzorów (2.54)
odpowiednio dla i < j oraz i > j.
Problemem może wydawać się obliczenie ostatniego wyrazu diagonalnego macierzy
odwrotnej ze wzoru (2.54), gdyż znów mamy tam nie istniejący wyraz eN . Tutaj jednak,

w przeciwieństwie do (2.45), wystarczy przyjąć za [7], że eN = 0, a wtedy T-1 =
N
NN
u-1 = W1,N-1/W1,N . Widać przy tym podobieństwo do wzoru (2.44), z którego wynika
N
T-1 = W2,N /W1,N . Obydwa te wzory wynikają też w oczywisty sposób bezpośrednio
N
11
z (2.2).
3. Opis środowiska obliczeniowego (Program
TIPS)
3.1. y
ródła danych i algorytmów
Środowisko TIPS (Tridiagonal Inverse Problem Solver) służy do wyznaczania ma-
cierzy odwrotnej do symetrycznej macierzy trójdiagonalnej. Przeznaczone jest ono dla
R
wszystkich systemów operacyjnych rodziny Windows�, chociaż dostępne kody z
ródłowe
pozwalają na osobną kompilację środowiska w innym systemie operacyjnym.
Składa się z aplikacji głównej (tips.exe), programu losującego macierze trójdiago-
nalne (h matgen.exe), programu sprawdzającego poprawność wyniku (h tester.exe) oraz
szeregu osobnych programów obliczających macierz odwrotną.
TIPS jest aplikacją napisaną i skompilowaną przy pomocy pakietu Delphi 5, a po-
zostałe programy napisano w języku Fortran i skompilowano przy pomocy środowiska
Force w wersji 2.0.8 i dostępnego wraz z nim kompilatora G77 w wersji 0.5.25 z dodat-
kowymi programami pakietu binutils w wersji 2.15.91. Niektóre programy skompilowano
z fragmentami bibliotek LAPACK 3 oraz BLAS, które nie są do nich dołączone w postaci
z
ródłowej. Wykorzystano także bibliotekę Portable Network Graphics Delphi w wersji
1.564 napisaną przez Gustavo Daud ([e-mail]) i dostępną za darmo
pod adresem http://pngdelphi.sourceforge.net.
Wszystkie programy obliczeniowe wymienione są w pliku konfiguracyjnym tips.ini.
Plik ten użytkownik może dowolnie modyfikować, tzn. dodawać własne programy obli-
czeniowe pod warunkiem, że przestrzegają one reguł dotyczących nazewnictwa plików
wejściowych, wyjściowych oraz formatu i kolejności danych w tych plikach. Dokładny
opis tego interfejsu programistycznego, a także innych możliwości środowiska zawarty
jest osobno w pomocy do programu TIPS i nie jest tutaj przedstawiony szczegółowo.
Dane dla programu TIPS to macierze trójdiagonalne oraz dodatkowe parametry je-
śli wybrany algorytm numeryczny ich wymaga (np. tolerancje dla kryteriów zbieżno-
ści, maksymalna dozwolona liczba iteracji itp.). Macierze mogą być wylosowane przez
program TIPS (a dokładnie przez podprogram h matgen.exe), dostarczone z zewnątrz
w postaci plików tekstowych, albo obliczone przez TIPS na podstawie zadanego wzoru
(osobne wzory dla elementów diagonalnych i pozadiagonalnych).
Obliczenia wzorów w programie TIPS wykonywane są w rozszerzonej precyzji (80
bitów), a następnie zaokrąglane do podwójnej precyzji (64 bity). Aby zachować mak-
simum kompatybilności z innymi językami programowania wszelkie dane przekazywane
są między programami środowiska TIPS poprzez pliki tekstowe. Programy obliczeniowe
otrzymują dane podwójnej precyzji (17 dziesiętnych cyfr znaczących), wykonują oblicze-
22
nia w podwójnej precyzji i przekazują wyniki (tj. macierz odwrotną, czas obliczeń mie-
rzony w jednostkach 1/1000 sekundy oraz ewentualne błędy) do programu testującego,
który przeprowadza czasochłonne mnożenie macierzy, oblicza wybrane błędy numerycz-
ne (również w podwójnej precyzji) i umieszcza je w wielokolumnowym pliku tekstowym.
Plik ten może być potem czytany przez program TIPS, który rysuje dwuwymiarowe
wykresy wybranych testów (dowolne zależności między kolumnami pliku wynikowego).
Program TIPS (h matgen.exe) może generować symetryczne macierze trójdiagonalne
przy pomocy generatora pseudolosowego z [9]. Macierz tą można też uczynić diagonalnie
dominującą lub dodatnio określoną.
Algorytmy dla metod iteracyjnych zaczerpnięto z [15]. Wyjątki są dwa. Implementa-
cja metody Czebyszewa oparta jest w całości na [1]. Implementacja metody Richardsona
zawarta w [14], została dostosowana do konwencji używanych w pozostałych programach
napisanych w języku Fortran i może służyć także jako przykład takiego dostosowania.
W wielu programach obliczeniowych korzystam także z biblioteki Lapack [13], z której
pochodzą też w całości algorytmy obliczania rozkładów SVD (2 wersje), LDLT oraz
LU z częściowym wyborem elementu głównego. Eliminacja Gaussa-Jordana opiera się
na podprogramie gaussj z [5]. Rozkład QR jest nieco zmodyfikowanym algorytmem
dekompozycji z [11]. Pozostałe metody, tj. wyznacznikowa i różne metody  analityczne
(w 6 wersjach) oparte są na wzorach zawartych i częściowo wyprowadzonych w tej pra-
cy. Ich implementacja i przetestowanie były inspirowane wynikami prac [7, 8] oraz [6].
Kompletny spis dołączonych do programu TIPS programów obliczeniowych (ich nazwy
zaczynają się od p ) zawarty jest w tabeli na następnej stronie.
23
Program Metoda y
ródło/Wzory
p wyzn.exe wyznacznikowa (2.3, 2.5, 2.7)
p gauss.exe eliminacja Gaussa / rozkład DGTSV [13] ze zmiana-
LU mi
p dgtsv.exe eliminacja Gaussa z wyborem DGTSV [13]
elementu głównego / rozkład
LU
p gsjor.exe eliminacja Gaussa-Jordana oparta na gaussj [5]
p llt1.exe rozkład LLT (2.21)
p dptsv.exe rozkład LDLT DPTSV [13]
p qr1.exe rozkład QR (2.26), oparta na [11]
p svd1.exe rozkład SVD / metoda QR DGBBRD, DBDSQR
[13]
p svd2.exe rozkład SVD / metoda  dziel DGBBRD, DBDSDC
i rządz
 [13]
p alg1.exe  analityczna 1A (2.47), (2.48), (2.53), [8]
p alg2.exe  analityczna 1B (2.47), (2.48), (2.53),
(2.49), [8]
p alg3.exe  analityczna 1C (2.47), (2.48), (2.53),
(2.49), [8]
p alg4.exe  analityczna 2A (2.52), (2.53), (2.54), [7]
p alg5.exe  analityczna 2B (2.52), (2.53), (2.54), [7]
p alg6.exe  analityczna 2C (2.52), (2.53), (2.54), [7]
p jacobi.exe Jacobiego JacobiREVCOM[15]
p richar.exe Richardsona msSolveRichardson [14]
p gseid.exe Gaussa-Seidla SORREVCOM [15]
p sor.exe nadrelaksacji SOR SORREVCOM [15]
p cgrad1.exe gradientów sprzężonych (bez CGREVCOM [15]
uwarunkowania wstępnego)
p cgrad2.exe gradientów sprzężonych CGREVCOM [15]
(diagonalne uwarunkowanie
wstępne)
p cheby2.exe Czebyszewa (2.38) [1]
3.2. Ocena dokładności obliczeń
Numeryczne obliczenie macierzy odwrotnej do macierzy trójdiagonalnej T jest rów-
noważne rozwiązaniu układu równań macierzowych
Ć
TX = I, (3.1)
24
Ć Ć
gdzie X H" T-1. Dokładność wyznaczenia macierzy X możemy sprawdzać obliczając
reszty [4]
Ć
XT - I reszta lewostronna [ang. left residual],
(3.2)
Ć
TX - I reszta prawostronna [ang. right residual].
Z definicji macierzy odwrotnej wynika przemienność mnożenia A-1A = AA-1 = I
[10]. Ponieważ macierz T jest symetryczna, to spodziewamy się, że reszty (3.2) będą
Ć
jednakowe, jeśli tylko obliczona przybliżona macierz X także będzie symetryczna. Pa-
miętajmy jednak, że w trakcie mnożenia macierzy wykonujemy także dodawanie, a to
nie jest przemienne w sensie numerycznym i podczas dodawania bliskich siebie liczb
przeciwnych mogą wystąpić duże błędy obcięcia [ang. cancelation].
Ć
Nie możemy bezpośrednio obliczyć tzw. błędu bieżącego [ang. forward error] T-1-X,
bo nie znamy dokładnej wartości T-1. Czasem bada się algorytmy numeryczne wyko-
nując obliczenia w większej precyzji, a następnie zaokrągla się uzyskiwane wyniki sy-
mulując w ten sposób arytmetykę zmiennoprzecinkową o mniejszej precyzji. W takich
kontrolowanych warunkach uprawnione jest potraktowanie wyników dokładniejszych ja-
ko prawdziwe i oszacowanie przy ich pomocy błędu bieżącego [4].
W programie TIPS także możemy dysponować dokładniejszym rozwiązaniem korzy-
stając z procedury DGTRFS z biblioteki Lapack [13]. Procedura ta stosuje metodę itera-
cyjnego udokładniania [ang. refinement] rozwiązania, która nie tylko pozwala polepszyć
dokładność wyniku uzyskanego każdą inną metodą, ale także udostępnia oszacowania
błędu bieżącego FERR i wstecznego błędu komponentowego BERR. Opis definicji i do-
kładnego sposobu obliczania tych błędów znajdziemy w [13].
Reszty (3.2) są macierzami i porównywać je można w sensie dowolnej normy (przez
równoważność norm [1, 2, 13]) albo komponentowo [2, 4]. Dalej (podobnie jak w [13])
jako normy macierzy będziemy używali normy � " (tzw. normy nieskończonej, zwanej
też wierszową)
N

A " = max |Aij|. (3.3)
1 i N
j=1
Jest ona najmniejszą znaną normą macierzową, bo zachodzą związki [13]
" "
A " N A 1 oraz A " N A 2 N A F . (3.4)
Program TIPS może obliczać także pozostałe normy z (3.4):
N

" normę kolumnową (normę jeden) A 1 = max |Aij|,
1 j N
i=1
Ax
" normę spektralną (normę dwa) A 2 = max = �max(A),
x =0 x

" normę Frobeniusa A F = |Aij|2.
ij
Norma spektralna wyznaczana jest jako maksymalna wartość szczególna macierzy,
która jest szacowana przez pierwiastek kwadratowy z maksymalnej wartości własnej
symetrycznej, dodatnio określonej macierzy AT A. Specjalnie do znalezienia tej wartości
25
własnej zaimplementowano w programie (h tester.exe) metodę potęgową [1, 11, 14]. Ob-
liczenia te mogą być czasochłonne, dlatego w programie TIPS są domyślnie wyłączone.
Program TIPS oblicza lewo- i prawostronne błędy względne
Ć
XT - I "
LERR := , (3.5)
Ć
T " X "
Ć
TX - I "
P ERR := , (3.6)
Ć
T " X "
oraz błędy względne komponentowe (lewo- i prawostronne)
Ć Ć
CLERR := min{ : "ij|(XT - I)ij| |Tij||Xij|}, (3.7)
Ć Ć
CP ERR := min{ : "ij|(TX - I)ij| |Tij||Xij|}. (3.8)
Wzory (3.7) i (3.8) wymagają wyjaśnienia. Nazwa  komponentowy przyjęła się z ję-
zyka angielskiego. Komponent to po prostu element lub wyraz macierzy. Dla każdego
elementu i-tego wiersza i j-tej kolumny macierzy liczymy iloraz
|Rij|
,
Ć
|Tij||Xij|
a błąd względny komponentowy definiuje się jako największy z tych ilorazów. Powyż-
sza macierz R jest macierzą reszty lewo- albo prawostronnej (3.2). Problem polega na
Ć
tym, że nie policzymy tego ilorazu jeśli elementy Tij lub Xij są równe zero. Dlatego
właśnie zamiast ilorazów we wzorach (3.7) i (3.8) mamy iloczyny. Myląca jest funkcja
minimum występująca w tych wzorach, gdyż kwantyfikator " (dla każdego) niejako za-
mienia ją w maksimum. Jeśli dla pewnego elementu macierzy zdarzy się przypadkiem, że
Ć
Rij = Tij = Xij = 0, to jedyne możliwe wynosi zero, ale kwantyfikator " mówi nam, że
musiałoby tak być dla każdego elementu macierzy, więc takie przypadkowe zero zostanie
zignorowane (tzn. nie będzie stanowiło minimum). Błąd komponentowy jest narzędziem,
które nie tylko pozwala zignorować przypadkowe zera (w tym przypadku norma byłaby
po prostu mniejsza i wynik wydawałby się dokładniejszy niż jest w rzeczywistości), ale
celowo ignorować zera, które nie są przypadkowe. Dlatego błędy komponentowe (3.7)
i (3.8) mają przewagę nad błędami w sensie normy. Mogą one być duże podczas gdy błę-
dy w sensie normy pozostają bardzo małe. Błędy komponentowe pozwalają oszacować
maksymalne zaburzenie, jakie może zaistnieć po zastosowaniu danego algorytmu, a nie
ogólny wpływ tych zaburzeń na normę. W przypadku macierzy rzadkiej jaką jest na
przykład macierz trójdiagonalna, nie interesuje nas wpływ zaburzeń zerowych wyrazów
na wynik, bo w żadnej praktycznej sytuacji takie zaburzenie nie wystąpi (tych wyrazów
po prostu nie ma). Gdy obliczamy błędy w sensie normy, to poprzez uśrednianie (li-
czymy przecież sumy modułów) uwzględniamy takie zaburzenia, natomiast licząc błędy
komponentowe nie bierzemy niejawnie pod uwagę zaburzeń wyrazów macierzy równych
zero z definicji.
26
3.3. Interfejs użytkownika
3.3.1. Informacje podstawowe
Główne okno programu TIPS po uruchomieniu wygląda tak:
Rysunek 3.1. Wygląd okna głównego programu TIPS
Przed uruchomieniem obliczeń (F9) należy ustawić opcje obliczeń. Opcje te podzie-
lone są tematycznie na trzy zakładki: Macierze  opcje dotyczące macierzy i odpowie-
dzialne za dostarczenie danych dla algorytmów numerycznych, Algorytmy i wyniki
 opcje dotyczące algorytmów numerycznych i sposobu zapisania obliczonych macie-
rzy na dysku, Testy i wykresy  opcje dotyczące wyboru przeprowadzanych testów
i wykresów przedstawiających wyniki tych testów.
27
3.3.2. Zapis i odczyt opcji obliczeń
Zapis opcji odbywa się po naciśnięciu F5 albo przycisku  Zapisz w górnej części okna
TIPS. Opcje zapisywane są domyślnie do plików tekstowych z rozszerzeniem  .tips , któ-
re zostają związane z programem TIPS przy jego instalacji. Po podwójnym kliknięciu
R
pliku z rozszerzeniem  .tips w oknie Eksploratora Windows� otwiera się program TIPS
z opcjami już wczytanymi z tego pliku. Opcją jest wszystko to, co ustawiamy w oknie
głównym TIPS na trzech zakładkach. Opcjami obliczeń nie są więc wykresy i parametry
ich wyświetlania ani też parametry zapisu wyników do plików i nie są one zapisywa-
ne w plikach z rozszerzeniem  .tips . Dodatkowe parametry metod obliczeniowych są
zapisywane do plików  .tips . Ustawione opcje można przywrócić do ustawień domyśl-
nych (początkowych) przez naciśnięcie F4 albo przycisku  Nowy w górnej części okna
TIPS. Opcje początkowe nie pozwalają na natychmiastowe uruchomienie obliczeń, gdyż
konieczne minimum stanowi wybranie algorytmu, typu macierzy i sposobu dostępu do
danych (czy i jak są losowane, czy są wczytywane z pliku, czy są obliczane ze wzoru).
Wcześniej zapisane opcje można wczytać do programu TIPS naciskając F3 albo przycisk
 Otwórz w górnej części okna TIPS.
3.3.3. Zakładka  Macierze
Określanie rodzaju macierzy
Rodzaj macierzy określa w polu wyboru z rysunku 3.3. Macierz dowolna symetryczna
powstaje po prostu przez wylosowanie liczb z podanego rozkładu. Wybranie macierzy
diagonalnie dominującej powoduje, że najpierw odbywa się losowanie macierzy dowol-
nej symetrycznej z rozkładu jednostajnego (di " (dod, ddo), ei " (eod, edo)), a potem
wylosowane elementy są przeliczane według wzoru z rysunku 3.4.
Macierz dodatnio określona jest tworzona przez mnożenie macierzy bidiagonalnej
dolnej przez jej transpozycję (odwrócenie rozkładu LLT ). Liczone są wartości własne
danej macierzy przy pomocy procedury DSTEQR z biblioteki Lapack [13]. Jeśli oka-
że się, że wygenerowana macierz posiada jednak wartości własne ujemne, to widmo
tej macierzy jest przesuwane o podwojoną wartość bezwzględną najmniejszej (ujemnej)
wartości własnej.
28
Rysunek 3.2. Wygląd zakładki  Macierze
Rysunek 3.3. Rodzaje macierzy do wyboru
29
Rysunek 3.4. Wzory do przekształcenia macierzy w macierz diagonalnie dominującą
Wybranie numerów K macierzy
Program TIPS może wykonywać wiele obliczeń macierzy odwrotnych przy jednym
uruchomieniu. Każde takie obliczenie ma swój numer oznaczony przez liczbę K. Numeru
K można używać w polach edycyjnych we wzorach aby zmieniać wartości parametrów
w kolejnych obliczeniach. Określany jest numer K minimalny i maksymalny. Podczas
obliczeń numery K zmieniają się co jeden od najmniejszego, do największego.
Wybranie rozmiarów N macierzy
Przy każdym obliczeniu, razem z nowym K, można zmienić rozmiar macierzy. Wy-
starczy wpisać wzór, zawierający K. Rozmiary macierzy N mogą się zmieniać od 1 do
1000.
Określanie parametru ISEED
Parametr ISEED oznacza liczbę całkowitą, użytą do zainicjalizowania generatora
liczb pseudolosowych. Są dwie osobne wartości ISEED: jedna dla elementów diago-
nalnych, a druga dla elementów pozadiagonalnych macierzy trójdiagonalnej. Jeśli nie
podamy żadnej wartości, to generator zostanie zainicjalizowany liczbą wybraną (także
wylosowaną) przez program TIPS. Określenie parametru ISEED nie jest więc potrzebne
jeśli chcemy wykonywać kolejne obliczenia dla innych losowych macierzy. Jeśli jednak
podamy jakąś wartość dla ISEED, to będzie ona ustawiana przed każdym obliczeniem
(dla każdego K), więc generowane za każdym razem macierze będą takie same (chyba że
wzrośnie ich rozmiar, ale wtedy będą się różnić tylko ostatnimi elementami). Określenie
niezerowego parametru ISEED umożliwia też uruchomienie różnych algorytmów dla tej
samej macierzy losowej.
Parametry rozkładów
Jak widać na rysunku 3.2 dla rozkładu jednostajnego mamy dwa parametry od i do,
które wyznaczają przedział losowanych wartości. Rozkład Gaussa ma także dwa para-
metry: wartość średnią m i odchylenie standardowe od średniej �. Liczby losowane są
przy pomocy generatorów pseudolosowych z książki [9].
30
Macierze o elementach obliczonych ze wzoru
Zamiast losować macierze, można pozwolić programowi TIPS obliczyć ich elementy
z podanego przez nas wzoru. We wzorze tym możemy użyć czterech zmiennych: K, N, I
oraz RAND. Zmienna I przyjmuje wartości od 1 do N i zmienia się dla każdego elementu
diagonali (dla elementów nad- i poddiagonali maksymalna wartość to oczywiście N - 1).
Zmienna RAND przyjmuje wartość losową o rozkładzie jednostajnym na przedziale
(0, 1). Jej pseudolosowość także jest kontrolowana przez parametr ISEED. Podanie tej
samej wartości ISEED generuje te same ciągi liczb. Przy wprowadzaniu wzoru nie można
podać parametru ISEED. Żeby podać parametr ISEED, trzeba wybrać macierz losową,
wpisać wartość ISEED dla di albo ei i powrócić do wyboru opcji  wg wzoru .
Macierze wczytane z plików
Macierz może zostać wczytana z pliku tekstowego. Należy podać osobne pliki dla
elementów diagonalnych i pozadiagonalnych. Nazwa pliku może zmieniać się wraz z nu-
merem K. Sekwencja %K% wstawia do nazwy wartość K jako tekst. Jeśli chcemy, aby
wstawiona wartość K miała stałą liczbę cyfr (zera z lewej strony), to np. dla wartości
4-cyfrowej wpisujemy %000K% (trzy zera przed K). Pliki o podanych nazwach powinny
się znajdować w katalogu z programami obliczeniowymi (domyślnie to podkatalog  prg
katalogu, w którym zainstalowano program TIPS).
3.3.4. Zakładka  Algorytmy i wyniki
Wybranie algorytmu
Algorytm wybiera się w polu wyboru zatytułowanym  Algorytm . Na rysunku 3.6
przedstawiona jest lista dostępnych do wybrania algorytmów.
Dodatkowe parametry dla metody (algorytmu)
Jeśli wybrany algorytm potrzebuje dodatkowych parametrów, to pod oknem wybo-
ru algorytmu pojawia się arkusz z nazwami parametrów dodatkowych dla wybranego
algorytmu i ich domyślnymi wartościami. Wartości te można zmieniać. Zmiany te są za-
pisywane w plikach z rozszerzeniem  .tips po naciśnięciu przycisku  Zapisz albo klawisza
F5.
Przycisk i okno  Opcje wyników
Przycisk ten pokazuje okno  Opcje wyników z rysunku 3.7. Opcje te są opcjami
programu TIPS i nie są zapisywane razem z innymi opcjami w plikach z rozszerzeniem
 .tips . Dane (symetryczna macierz trójdiagonalna) znajdują się w plikach fort.8 i fort.9,
a wyniki (macierz odwrotna do symetrycznej macierzy trójdiagonalnej) w pliku fort.10.
Jeśli jednak obliczenia są powtarzane dla różnych wartości K, to dane z tych plików
są utracone, a dostępne są tylko te, które powstały dla największego (ostatniego) K.
Aby temu zapobiec i aby umożliwić zgromadzenie wszystkich danych i wyników z danej
sekwencji obliczeń (dla różnych K) na dysku, program h tester.exe zapisuje wyniki do
31
Rysunek 3.5. Zakładka  Algorytmy i wyniki
32
Rysunek 3.6. Lista dostępnych do wyboru algorytmów
Rysunek 3.7. Okno  Opcje wyników
33
plików fort.18, fort.19 i fort.21, które mogą być następnie przeniesione do dowolnego
innego katalogu pod dowolną inną nazwą. Po drodze przeprowadzane są transformacje
danych i wyników do wybranego formatu. Format liczb się nie zmienia, ale zmienia się
ich ułożenie w pliku. Wyboru można dokonać między formatami dla danych (rysunek
3.8) i dla wyników (rysunek 3.9). Jak widać można zrezygnować z zapisu całkowicie,
a nawet usuwać pliki o pasujących nazwach nie zastępując ich nowymi plikami.
Rysunek 3.8. Opcje zapisu pliku(-ów) z macierzą trójdiagonalną
Rysunek 3.9. Opcje zapisu pliku z macierzą odwrotną
Reguły nazewnictwa plików
Jeśli wybrano format zapisu danych w osobnych plikach, to w polu edycyjnym  Pliki
macierzy trójdiagonalnych należy podać dwie nazwy oddzielone przecinkiem. W na-
zwach można używać wartości zmiennych K i N zamienionych na tekst. Sekwencja %K%
wstawia do nazwy wartość zmiennej K, a sekwencja %N% wstawia wartość N. Jeśli
chcemy, aby wstawiona wartość miała stałą liczbę cyfr (zera z lewej strony), to np. dla
wartości 4-cyfrowej wpisujemy %000K% (trzy zera przed K). Tak samo można wstawić
zmienną N, czyli rozmiar macierzy. Pliki mogą się znajdować w innym katalogu, ale
katalog ten musi wcześniej istnieć (nie jest tworzony).
34
Rysunek 3.10. Zakładka  Testy i wykresy
35
3.3.5. Zakładka  Testy i wykresy
Wybór testów do przeprowadzenia
Po uzyskaniu wyniku jakim jest macierz odwrotna uruchamiany jest program testu-
jący h tester.exe. Wynikiem jego działania jest tzw. plik wynikowy. Każdy wiersz tego
pliku odpowiada kolejnemu obliczeniu (innemu K). W kolumnach znajdują się wybrane
testy. Testy FERR i BERR pochodzą z wywołania procedury DGTRFS z biblioteki
Lapack [13] po obliczeniu wyniku daną metodą. Procedura ta udokładnia otrzymane
rozwiązanie i jeśli wybierzemy co najmniej jeden z tych testów, to otrzymana macierz
odwrotna będzie już macierzą dokładniejszą, ale pozostałe testy są przeprowadzane wcze-
śniej i dotyczą macierzy mniej dokładnej, będącej bezpośrednim wynikiem wybranego
algorytmu obliczeniowego. Jeśli natomiast wybierzemy ostatnie dwa testy preFERR lub
preBERR, to samo udokładnianie wyniku zostanie uruchomione przed obliczaniem po-
zostałych testów, więc ich wyniki będą dotyczyły już udokładnionego rozwiązania. testy
FERR i BERR nie są przeprowadzane, jeśli wybrano testy preFERR lub preBERR. Aby
wybrać dany test, należy zaznaczyć pole wyboru z lewej strony jego nazwy. Po zazna-
czeniu testu nadawany mu jest numer, który pojawia się z przodu nazwy testu przed
kropką. Jest to numer kolumny w pliku wynikowym. W tej kolumnie zostaną zapisane
wyniki danego testu dla każdego kolejnego K. Testy nie wybrane nie są przeprowadzane,
chyba że ich przeprowadzenie jest konieczne do obliczenia wyników innych wybranych
testów. Dzięki temu wybranie mniejszej ilości testów spowoduje przyspieszenie testowa-
nia wyniku przez program TIPS. Jeśli obliczenia są już przeprowadzone lista  Kolumny
pliku wynikowego jest  wyszarzona (tzn. niedostępna) i nie można już zmieniać kolumn
pliku wynikowego. Ponieważ plik wynikowy istnieje, opcje obliczeń zostają z nim zwią-
zane. Powiązanie to, czyli nazwa pliku wynikowego, może być zapisane razem z opcjami
obliczeń w pliku  .tips .
Dodawanie nowych wykresów
Aby po obliczeniach zobaczyć wyniki wybranych testów, należy dodać wykresy. Doda-
nie wykresu polega na wybraniu numerów kolumn pliku wynikowego, w których znajdują
się wartości X (poziome) i Y (pionowe) dwuwymiarowego wykresu. Dodajemy wykre-
sy naciskając przycisk  Nowy wykres (rysunek 3.10). Pojawia się wtedy okno  Nowy
wykres (rysunek 3.11). Nowe wykresy (ich nazwy) pojawiają się na zakładce  Testy
i wykresy na liście  Wykresy . Podwójne kliknięcie na nazwie wykresy na tej liście
(albo zaznaczenie nazwy i kliknięcie przycisku  Pokaż wykres ) spowoduje pojawienie się
okna  Wykres zawierającego dany wykres. Wykresy można też usuwać. W tym celu po
zaznaczeniu wykresu do usunięcia w oknie  Wykresy naciskamy przycisk  Usuń wykres .
Nie trzeba usuwać wykresów aby zamykać okna z wykresami. Zmieniając kolumny X i Y
wykresu automatycznie wpisywane są domyślne podpisy osi wykresu i nazwa wykresu
umieszczona na górze okna  Wykres Nazwę tę i podpisy można zmienić teraz, albo póz
-
niej poprzez edycję parametrów wykresu. Kliknięcie przycisku  OK spowoduje dodanie
wykresu do listy  Wykresy na zakładce  Testy i wykresy . Można potem obejrzeć dany
wykres klikając podwójnie na jego nazwę na tej liście.
36
Rysunek 3.11. Okno  Nowy wykres
Dodawanie wielu wykresów jednocześnie
Ponieważ dodawanie wykresów jest dość czasochłonną czynnością można automa-
tycznie dodać wiele wykresów jednocześnie przez naciśnięcie przycisku  Wszystkie z tą
kolumną X (rysunek 3.11). Dodane zostaną wszystkie wykresy, które mają aktualnie
wybraną w tym oknie kolumnę X i każdą z kolumn Y (innych niż kolumna X). Nie są
wybierane kolumny, których nie ma w pliku wynikowym.
Rysunek 3.12. Okno  Edytuj wykres
37
Edytowanie parametrów wykresu
Po naciśnięciu przycisku  Edytuj wykres na zakładce  Testy i wykresy pokaże się
okno  Edytuj wykres (rysunek 3.12). Różni się ono od okna  Nowy wykres tylko tym,
że nie ma przycisku  Wszystkie z tą kolumną X . Przeznaczenie pozostałych elementów
okna jest takie samo. Naciśnięcie przycisku  OK nie dodaje jednak nowego okna, ale
zmienia parametry (może nawet kolumny X i Y) wykresu istniejącego i aktualnie wybra-
nego w oknie  Wykresy . Dodane wykresy i ich parametry nie są zapisywane w plikach
z rozszerzeniem  .tips i trzeba je dodawać od nowa po każdym wczytaniu opcji obliczeń
z takiego pliku.
Rysunek 3.13. Okno  Wykres
3.3.6. Okno  Wykres
Pokazanie wykresu
Aby pokazać wykres można kliknąć podwójnie lewym przyciskiem myszy na nazwie
wykresu w oknie  Wykresy , albo zaznaczyć nazwę pojedynczym kliknięciem i nacisnąć
przycisk  Pokaż wykres . Pokaże się wtedy okno wykresu podobne do tego z rysunku 3.13.
Pokazanie wykresu jest możliwe dopiero po obliczeniach (wtedy istnieje plik wynikowy
i jego kolumny faktycznie zawierają wyniki wybranych wcześniej testów), albo po wczy-
taniu pliku z opcjami obliczeń, kiedy ten plik był zapisany po tym jak obliczenia były
już wykonane (wtedy pozostało powiązanie opcji z plikiem wynikowym). Wykresy są
rysowane na podstawie danych z dwóch kolumn pliku wynikowego. Plik ten jest czytany
przy każdym wyświetlaniu okna z wykresem. Kolumny wybierane są w oknie  Edytuj
wykres (rysunek 3.12).
38
Wygląd wykresu
Rozmiary okna z wykresem można dowolnie zmieniać w jednym kierunku chwyta-
jąc i przeciągając myszką jedną z jego krawędzi, albo w dwóch kierunkach chwytając
i przeciągając jeden z rogów okna. Wykresy są  wykresami punktowymi . Poszczególne
punkty nigdy nie są łączone linią. Liczby w opisach osi nigdy nie mają więcej niż 4
cyfry (2 cyfry po kropce). Jeśli zapis liczby wymaga więcej cyfr, cała oś jest rysowana
w skali, przy czym mnożnik jest zapisany w nawiasie kwadratowym obok osi przy jej
końcu. Nazwa wykresu i podpisy osi mogą być zmienione w oknie  Edycja wykresu . Osie
wykresu nie przecinają się w punkcie (0, 0). Skala nie jest rysowana poniżej minimalnej
i powyżej maksymalnej wartości współrzędnych rysowanych punktów. Dzięki temu moż-
na rysować wykresy bardzo małych wartości (np. 1E-17) bez obawy, że nie zobaczymy
drobnych różnic w pobliżu zera. Jeśli rozrzut wartości jest duży, można je narysować
w skali logarytmicznej.
Rysunek 3.14. Menu podręczne do okna  Wykres
Menu podręczne
Menu z rysunku 3.14 pojawia się po kliknięciu prawym klawiszem myszy w obsza-
rze okna z wykresem. Dostępne tutaj opcje pozwalają manipulować opcjami wykresu
niedostępnymi nigdzie indziej w programie TIPS i nie zapisywanymi nigdzie (po za-
mknięciu okna z wykresem trzeba je ustawiać na nowo). Każdą z osi można ustawić
jako logarytmiczną, a jeśli istnieją punkty o współrzędnych równych zero, to można
je usunąć z wykresu. Do wykresu można też dodać (dorysować) czerwoną linię regresji.
Regresja liniowa pozostaje liniowa nawet jeśli co najmniej jedna z osi jest logarytmiczna.
Parametry regresji się wtedy zmieniają i pozwala to na wykorzystanie regresji zarówno
na wykresie zwykłym, jak i logarytmicznym. Na rysunku 3.15 pokazano przykładowy
wygląd linii regresji i okna z parametrami regresji dla tej linii. Okno to pokazuje się
również wskutek wybrania odpowiedniej pozycji z menu podręcznego.
Kiedy podejmiemy już decyzje czy osie mają być logarytmiczne i czy wyświetlić
linię regresji, to możemy zapisać obraz wykresu do pliku graficznego w formacie PNG
(Portable Network Graphics).
39
Rysunek 3.15. Okno  Wykres z linią regresji i okno  Parametry regresji
40
Skale logarytmiczne
Zarówno skala Y jak i skala X mogą być logarytmiczne. Oznacza to w tym wypad-
ku, że zamiast właściwych wartości rysowane są ich logarytmy dziesiętne. Opisy osi nie
zawierają jednak całych potęg 10, ale ułamkowe wartości właściwych logarytmów. Po-
nadto podpis osi, która jest logarytmiczna zostaje automatycznie poprzedzony tekstem
 log . Na osiach logarytmicznych nie jest rysowana skala, tzn. nie mają one mnożni-
ka w nawiasie kwadratowym oznaczającego potęgę 10, którą trzeba przemnożyć przez
wszystkie wartości na osi aby otrzymać właściwy wynik. Wartości wszystkich testów są
(powinny być) nieujemne, ale mogą tam wystąpić zera, szczególnie dla małych rozmiarów
N macierzy. Skala nie może być logarytmiczna gdy musi zawiewrać wartość zero (gdy
są punkty o współrzędnej zero). Takie punkty można jednak usunąć z wykresu (dane
w pliku wynikowym pozostają nie zmienione), a następnie narysować wykres logaryt-
miczny z pozostałych punktów. Usuniętych punktów nie można już przywrócić, chyba
że zamkniemy okno z wykresem i otworzymy je na nowo (przycisk  Pokaż wykres ),
co spowoduje ponowne odczytanie pliku wynikowego. Na przykładzie z rysunku 3.16
pokazano jak wykorzystać osie logarytmiczne i linię regresji do odczytania złożoności
czasowej algorytmu (tutaj jest to O(N2.24)).
3.3.7. Okno  Obliczenia
Okno z rysunku 3.17 pojawia się na czas wykonywania obliczeń. Po uruchomieniu
obliczeń okno główne zostaje zminimalizowane, a na środku ekranu pojawia się poniższe
okno z paskiem postępu, który przesuwa się K razy, raz po każdej macierzy w serii obli-
czeń. Jeśli podczas obliczeń naciśnięty zostanie przycisk zamknięcia okna, to obliczenia
zostaną wstrzymane i pojawi się pytanie z rysunku 3.18. Naciśnięcie  NIE spowoduje
wznowienie obliczeń i ich dokończenie, a odpowiedz
 TAK przerwie obliczenia i wyświe-
tli okno z potwierdzeniem przerwania obliczeń. Nie będzie można póz
niej dokończyć tak
przerwanych obliczeń (można zawsze uruchomić je od początku, albo dla zmienionych
opcji).
3.3.8. Korzystanie ze wzorów
Liczby
Wszystkie liczby zostaną zamienione na zmiennoprzecinkowe podczas obliczania wzo-
ru. Jeśli wynikiem wzoru ma być liczba całkowita, to zostanie on zaokrąglony przez
obcięcie części ułamkowej (czyli w stronę zera na osi liczbowej). Do oddzielenia części
ułamkowej używamy kropki, a nie przecinka. Można stosować w zapisie liczb notację
naukową, np. 1.345e-19. Wszystkie liczby we wzorach są zmiennoprzecinkowe o typie
rozszerzonym 80-bitowym. Podczas obliczania wyrażeń nie ma innego typu. Są jednak
dwa wyjątki: logarytm naturalny LN() i funkcja znaku SGN(). Argument tych funkcji
nie może być bliżej zera niż 1E-16. Jeśli będzie, spowoduje błąd.
41
Rysunek 3.16. Przykład wykorzystania osi logarytmicznych i linii regresji
Rysunek 3.17. Okno  Obliczenia
42
Rysunek 3.18. Wstrzymanie obliczeń i pytanie o potwierdzenie chęci ich przerwania
Zmienne i funkcje
We wzorach można używać nazw zmiennych oraz predefiniowanych nazw funkcji. Nie
można definiować własnych funkcji. Wielkość liter jest bez znaczenia. Poniżej przedsta-
wiona jest lista dostępnych funkcji.
" SQRT() pierwiastek kwadratowy
" ABS() wartość bezwzględna
" RAD() zamiana stopni na radiany (funkcje trygonometryczne przyjmują argumenty
w radianach)
" DEG() zamiana radianów na stopnie
" INV() odwrotność liczby
" EXP(x) liczba e do potęgi x
" LN() logarytm naturalny (podstawa jest liczbą e)
" SGN() znak liczby (-1,0 lub 1)
" CEIL() zaokrąglenie do liczby całkowitej w górę
" FLOOR() zaokrąglenie do liczby całkowitej w dół
" FRAC() część ułamkowa liczby
" MOD(a,b) dzielenie modulo
" POW(a,b), POWER(a,b) a do potęgi b
" MIN(a,b) mniejsza z liczb
" MAX(a,b) większa z liczb
" IF(a,b,c) funkcja IF, jeśli a = 0 to c inaczej b
" SIN() sinus
" COS() kosinus
" TG(), TAN() tangens
" CTG(), CTAN(), COTAN(), COT() kotangens
" SEC() sekans
" COSEC(), CSC() kosekans
" ARCSIN(), ASIN() arkus sinus
" ARCCOS(), ACOS() arkus kosinus
" ARCTG(), ATG(), ARCTAN(), ATAN() arkus tangens
" ARCCTG(), ACTG(), ACTAN(), ARCCOT(), ACOT() arkus kotangens
43
" ARCSEC(), ASEC() arkus sekans
" ARCCSC(), ACSC(), ARCCOSEC(), ACOSEC() arkus kosekans
" SH(), SINH() sinus hiperboliczny
" CH(), COSH() kosinus hiperboliczny
" TH(), TGH(), TANH() tangens hiperboliczny
" CTH(), CTGH(), CTANH(), COTH() kotangens hiperboliczny
" SECH() sekans hiperboliczny
" CSCH(), COSECH() kosekans hiperboliczny
Operatory
We wzorach oprócz liczb, zmiennych i funkcji można używać wielu operatorów ma-
tematycznych. Należy pamiętać, że operatory dwuargumentowe są lewostronnie łączne.
Oznacza to, że wyrażenia są obliczane od strony lewej do prawej. Wyrażeniu A+B+C
odpowiada (A+B)+C, a wyrażeniu A*B/C odpowiada (A*B)/C. Priorytet operatorów
jest zgodny z matematyką (wyjątkiem są operatory relacji). Poniżej znajduje się lista
dostępnych operatorów w kolejności od najbardziej wiążących:
+ - ! jednoargumentowe: plus, minus i negacja
(wyr) fun(wyr1,wyr2) nawias i wywołanie funkcji
daszek (Shift+6) potęgowanie
* / mod mnożenie, dzielenie i operacja modulo
+ - dodawanie i odejmowanie
Ą ż = Ą= ż= == Ąż != relacje
Operatory = i == działają jednakowo i oznaczają równość. Operatory Ąż i != dzia-
łają jednakowo i oznaczają różność. Lewostronna łączność nie dotyczy operatorów rela-
cji. Wynik relacji jest jedynką dla relacji prawdziwej albo zerem dla fałszywej. Gdyby
operatory relacji były lewostronnie łączne, to na przykład relacja (-1 < -10 < 1)
byłaby prawdziwa, a relacja (-1 < 0 < 1) fałszywa. Takie wielokrotne relacje są trak-
towane w specjalny sposób. Jeśli co najmniej jedna z relacji składowych jest fałszywa,
to cała relacja jest fałszywa. Na przykład relacja A < B < C jest traktowana jak
(A < B) " (B < C). Zwykły iloczyn działa na wyniki relacji jak logiczna koniunkcja
(AND), a zwykła suma relacji działa jak logiczna alternatywa (OR). Nie ma osobnych
operatorów AND i OR. Można jednak użyć zwykłego mnożenia i dodawania jeśli wymu-
simy konwersję liczby na wartość logiczną przez dwukrotny operator negacji, na przykład
!!A odpowiada funkcji if(A,1,0). Spotykaną w wielu językach programowania konstrukcję
(A OR B) można uzyskać tak (!!A + !!B). Podobnie (A AND B) zapisujemy jako (!!A
* !!B) albo też jako !!(A*B). Dla dodawania na ogół nie jest prawdą, że (!!A + !!B)
= !!(A+B). Funkcja if sama przeprowadza konwersję swojego pierwszego argumentu na
wartość logiczną, więc wyrażenia if((AĄB)*(CżD),1,2) oraz if(!!( (AĄB)*(CżD) ),1,2) dają
ten sam wynik. Funkcję if(A,B,C) można zapisać jako wyrażenie (!!A*B + !A*C).
4. Przeprowadzone testy
4.1. Porównanie czasów wykonania algorytmów
Algorytmy przetestowano na przykładzie macierzy Poissona (di = 2, ei = -1) o roz-
miarach od 100 do 1000. Jeśli zależność czasu obliczeń t od rozmiaru macierzy N jest
funkcją potęgową typu t(N) = ąN�, to po zlogarytmowaniu tego równania stronami ma-
my log t(N) = � log(N)+log(ą). Wykładnik potęgowy � zawarty w poniższych tabelach
został wyznaczony przy pomocy regresji liniowej tak, jak na przykładowym rysunku 3.16.
Czas wykonania metod iteracyjnych był mierzony aż do osiągnięcia dokładności 10-7.
Metody bezpośrednie
Metoda Wykładnik �
1. Wyznacznikowa 2.98
2. Eliminacja Gaussa 1 2.11
3. Eliminacja Gaussa 2 2.17
4. Eliminacja Gaussa-Jordana 3.41
5. Rozkład LLT 2.01
6. Rozkład LDLT 2.09
7. Rozkład QR 2.17
8. Rozkład SVD 1 2.89
9. Rozkład SVD 2 2.92
10.  Analityczna 1A 1.97
11.  Analityczna 1B 1.89
12.  Analityczna 1C 2.53
13.  Analityczna 2A 2.13
14.  Analityczna 2B 2.04
15.  Analityczna 2C 2.21
Metody iteracyjne
Metoda Wykładnik �
16. Jacobiego 3.58
17. Richardsona 3.19
18. Gaussa-Seidela 3.17
19. SOR 3.38
20. Gradientów sprzężonych 1 3.04
21. Gradientów sprzężonych 2 3.09
22. Czebyszewa 2.63
45
4.2. Porównanie dokładności - macierz Poissona
Norma (w sensie normy nieskończonej) macierzy odwrotnej do macierzy Poissona
rośnie jak O(N2). Powoduje to, że wskaz
nik uwarunkowania także rośnie jak O(N2)
i dla N = 1000 osiąga wartość 500000 (rysunek 4.1).
Rysunek 4.1. Wskaz
nik uwarunkowania macierzy Poissona
Dla metody wyznacznikowej nie obserwuje się dużych różnic między błędami lewo-
stronnymi i prawostronnymi. Zarówno błędy w sensie normy, jak i błędy komponentowe
są rzędu �mach dla macierzy Poissona o rozmiarach od N = 50 do N = 1000. Pozwala to
twierdzić, że metoda ta jest w tym przypadku numerycznie poprawna.
Obie metody eliminacji Gaussa (bez wyboru i z wyborem elementu głównego) zacho-
wują się podobnie, stąd wniosek, że wybór ten w przypadku macierzy Poissona nie jest
potrzebny. Występuje jednak znaczna różnica w charakterze zależności błędów (zarówno
komponentowych, jak i w sensie normy) od rozmiaru macierzy. Błędy lewostronne wy-
raz
nie rosną (rysunek 4.2), podczas gdy błędy prawostronne pozostają proporcjonalne
do �mach (rysunek 4.3). Ponieważ wraz z rozmiarem macierzy Poissona rośnie jej współ-
czynnik uwarunkowania, to nie wiemy czym dokładnie jest spowodowany wzrost błędu
lewostronnego.
Metoda Gaussa-Jordana nie wykazuje powyższych różnic między błędami lewo- i pra-
wostronnymi (wszystkie są rzędu �mach).
Błędy lewostronne metod korzystających z rozkładów LLT i LDLT wraz ze wzrostem
46
Rysunek 4.2. Błąd lewostronny metody Gaussa dla macierzy Poissona
Rysunek 4.3. Błąd prawostronny metody Gaussa dla macierzy Poissona
47
rozmiaru macierzy podobnie jak dla metody Gaussa, podczas gdy błędy prawostronne
oscylują wokół jednej wartości.
Ten sam wzrost błędów lewostronnych obserwujemy dla metody rozkładu QR, z tym
że błędy lewostronne są w tym przypadku prawie o 3 rzędy wielkości większe (prawo-
stronne pozostają rzędu �mach).
Obie metody wykorzystujące rozkład SV są do siebie bardzo podobne. Błędy le-
"D
wostronne rosną w nich proporcjonalnie do N, a błędy prawostronne są 4 razy większe
od lewostronnych (nie pozostają stałe, ale także rosną).
Wyniki trzech pierwszych wersji metody analitycznej (1A, 1B, 1C) dla macierzy
Poissona są dokładnie takie same (otrzymano identyczne liczby zmiennoprzecinkowe).
Błędy komponentowe lewostronne i prawostronne także są dokładnie takie same, ale
2.44
niestety dość duże, bo rzędu 10-9 i wykazują tendencję wzrostową typu N . Istnieją
jednak różnice w błędach względnych w sensie normy. Błędy lewostronne LERR są duże
(na poziomie 10-14), ale błędy prawostronne PERR są jeszcze większe (średnio 20 razy
większe).
Wyniki metod analitycznych 2A, 2B i 2C są do siebie podobne, ale nie są identyczne.
W przypadku metod 2A i 2B obserwuje się jednak identyczność błędów komponentowych
CLERR i CPERR dla tej samej metody. Oba te błędy, a także błędy w sensie normy
rosną liniowo wraz ze wzrostem N (rysunek 4.4). Błąd PERR osiąga 10-14 dla N = 1000,
ale błąd LERR jest od niego około 3.4 razy mniejszy. Błędy komponentowe, oprócz tego
że są identyczne dla mnożenia lewo- i prawostronnego, to są one tylko o rząd wielkości
większe i dla N = 1000 osiągają wartość 10-13.
Metoda 2C wymaga osobnego komentarza. Podczas gdy błędy lewostronne CLERR
i LERR są bardzo podobne do tych samych błędów w metodach 2A i 2B, błędy prawo-
stronne są znacznie mniejsze, rzędu �mach i nie wykazują liniowej tendencji wzrostowej!
Upodabnia to metodę  analityczną w wersji 2C do metod eliminacji Gaussa, jednak
"
tam błąd LERR rósł nie liniowo, ale jak N i dla N = 1000 był rzędu �mach (rysunek
4.2).
Błędy metod iteracyjnych można zbadać ustalając niemożliwą do osiągnięcia tole-
rancję dla maksymalnego błędu i wykonywać iteracje aż błąd przestanie maleć. Błąd
jednak nie zawsze maleje monotonicznie, a poza tym takie obliczenia byłyby bardzo
czasochłonne. Można też inaczej porównywać dokładność osiągniętą przez różne metody
iteracyjne po wykonaniu takiej samej liczby iteracji. Wtedy jednak oceniamy szybkość
zbieżności, a nie błąd metody.
4.3. Porównanie dokładności - macierz Hilberta
1 1
Macierz Hilberta ma postać: di = , ei = . Macierz ta nie jest dodatnio określona
2i-1 2i
i metod korzystających z rozkładów LLT i LDLT nie można zastosować. Wyniki metod
bezpośrednich porównano dla macierzy o rozmiarach od 5 do 100 z krokiem 5. Norma
macierzy odwrotnej do macierzy Hilberta przyjmuje wartości z jednej z trzech gałęzi, tak
jak to widać na rysunku 4.5. Gałąz
środkowa odpowiada rozmiarom macierzy podzielnym
przez 3, gałąz
górna zawiera macierze o rozmiarach N = 3m-1, a dolna dla N = 3m+1.
48
Rysunek 4.4. Liniowy wzrost błędu LERR metody 2C dla macierzy Poissona
Rysunek 4.5. Norma macierzy Hilberta
49
Dla metody wyznacznikowej błędy w sensie normy LERR i PERR maleją i pozostają
mniejsze od �mach gdy rozmiar N macierzy rośnie. Natomiast błędy względne kompo-
nentowe CLERR i CPERR są takie same i rosną jak N2.2 osiągając wartości rzędu 10-11
już dla N = 100.
Dla metod Gaussa błędy LERR i PERR również są bardzo małe, mniejsze niż �mach
i maleją gdy rozmiar macierzy rośnie od 5 do 100. Zawsze błąd prawostronny zdaje się być
nieznacznie mniejszy niż lewostronny. Jeśli nie wykonujemy wyboru elementu głównego,
to błędy CLERR są rzędu 10-11, a CPERR są około 8 razy mniejsze. Przy wyborze ele-
mentu głównego nie ma takiej jednoznacznej zależności, natomiast błąd LERR zmniejsza
się do poziomu błędu PERR. W metodzie Gaussa-Jordana błędy komponentowe nadal
są rzędu 10-11 dla N = 100, a błędy w normie nieskończonej są mniejsze od �mach.
Wciąż błędy komponentowe rosną co najmniej jak N2. Obserwujemy jednak, że w tej
metodzie błędy lewostronne są większe niż prawostronne (najwyżej 2 razy większe).
Metoda rozkładu QR jest porównywalna z eliminacją Gaussa, a obie metody rozkładu
SVD są o rząd wielkości gorsze.
Wyniki dla metod  analitycznych są do siebie bardzo podobne, ale istnieją jednak
pewne drobne różnice. Błędy CLERR i CPERR są jednakowe dla każdej metody z wy-
jątkiem 2C. Błędy PERR są porównywalne z LERR dla metod 1A i 2B, a są około 3 razy
większe dla metod 1B i 1C i 1.5 razy większe dla metody 2B. Metoda 2C jest wyjątkowa,
gdyż dla niej błędy PERR są 3 razy mniejsze niż LERR, co znów wyróżnia tą wersję
metody spośród pozostałych (rysunek 4.6).
Rysunek 4.6. Błędy PERR Metody 2C dla macierzy Hilberta
50
Rysunek 4.7. Prawostronne błędy komponentowe dla macierzy Hilberta
51
4.4. Porównanie dokładności - macierz losowa
Przeprowadzono także testy dla macierzy losowej, której elementy były losowane
z rozkładu jednostajnego na odcinku (0, 1). Obie metody SVD wykazują bardzo duży
błąd komponentowy i błąd w sensie normy o rząd wielkości większy niż �mach. Metoda
wyznacznikowa generuje macierz symetryczną, więc błędy komponentowe lewo- i pra-
wostronne są sobie równe. Błędy w sensie normy są mniejsze niż �mach. Metoda Gaussa
z częściowym wyborem elementu głównego jest dokładniejsza niż bez wyboru, ale po-
nieważ nie generuje macierzy symetrycznej, błąd LERR jest większy niż dla metody
wyznacznikowej (za to PERR jest nieco mniejszy). Metoda QR jest prawie o rząd wiel-
kości gorsza od metody Gaussa. Metody  analityczne 1A, 1B i 1C czasem nie dają
właściwego wyniku (nieskończoności), tak jak metoda rozkładu LLT nie dawała wyniku
dla macierzy, która nie jest dodatnio określona. Metoda 2B ma najmniejsze błędy kom-
ponentowe z metod 2A, 2B i 2C. Metoda 2A ma mniejsze błędy prawostronne podczas
gdy metoda 2B ma mniejsze błędy lewostronne. Metoda 2C ma nieco większe błędy
komponentowe od metody 2B, ale za to błędy w sensie normy są znacznie mniejsze
(PERR jest mniejszy o ponad 3 rzędy wielkości).
5. Podsumowanie
Na podstawie przeprowadzonych testów i obliczeń numerycznych można sformułować
następujące wnioski:
1. Metoda 2C jest w wielu przypadkach najdokładniejszą z badanych metod  analitycz-
nych , choć nie jest najszybsza. Najszybsza jest metoda 1B. Dla macierzy Poissona
błąd metody 2C rośnie liniowo wraz ze wzrostem rozmiaru macierzy N (rysunek 4.4).
Lepszy wynik uzyskujemy dla metody eliminacji Gaussa z częściowym wyborem ele-
"
mentu głównego (Eliminacja Gaussa 2), której błąd rośnie jak N (rysunek 4.2).
W najgorszym przypadku (macierz Hilberta) błąd komponentowy obu metod rośnie
w przybliżeniu jak N2 (rysunek 4.7).
2. Dla macierzy Hilberta obserwujemy wzrost błędów komponentowych wraz ze wzro-
stem rozmiaru macierzy (rysunek 4.7), podczas gdy błędy w sensie normy maleją
(rysunek 4.6). Być może wzrost błędów w sensie normy widoczny będzie dla macie-
rzy znacznie większych rozmiarów niż badane (100 dla macierzy Hilberta i 1000 dla
macierzy Poissona).
3. Dwie najdokładniejsze metody (tzn. Gaussa i 2C) nie generują symetrycznej ma-
cierzy odwrotnej. Jeśli dana metoda nie generuje symetrycznej macierzy odwrotnej,
to zwykle jeden z błędów (lewostronny lub prawostronny) jest znacznie mniejszy od
drugiego, ale nigdy nie są one porównywalne. Kiedy jakaś metoda generuje macierz
symetryczną, to najbardziej prawdopodobne jest, że wykorzystuje ona fakt syme-
tryczności wyniku po prostu kopiując wartości elementów xij znad (spod) diagonali
do odpowiadających im elementów xji (w ten sposób zostały zaimplementowane
wszystkie metody  analityczne oprócz 2C). Nie opłaca się więc tak formułować
algorytmów numerycznych odwracania macierzy, że wymuszają one symetryczność
wyniku przez zignorowanie połowy obliczeń.
4. Metod  analitycznych w wersjach 1A, 1B i 1C nie można stosować dla macierzy,
które nie są dodatnio określone. Obserwowano bardzo buże błędy, a nawet przekro-
czenie zakresu (nieskończoności) dla niektórych macierzy losowych. Z kolei metody
2A, 2B i 2C zachowują się poprawnie dla macierzy losowych. W teorii ich stabilność
jest jednak gwarantowana tylko dla macierzy diagonalnie dominujących [7].
5. Metody używające rozkładów macierzy, w których występują macierze pełne lub też
trójkątne pełne (SVD, QR) są znacznie mniej dokładne ze względu na większą liczbę
operacji zmiennoprzecinkowych, które muszą one wykonywać.
6. Metoda wyznacznikowa okazała się być stabilna numerycznie, co jest pewnym zasko-
czeniem. Dla macierzy Poissona jest nawet poprawna numerycznie. Złożoność tej me-
53
tody zbliża się jednak bardziej do złożoności metod iteracyjnych (O(N3)) niż innych
metod bezpośrednich, co raczej wyklucza możliwość jej praktycznego zastosowania.
7. Z metod iteracyjnych najszybsza okazała się metoda Czebyszewa. Może to wynikać
w faktu, że w implementacji tej metody obliczane są ekstremalne wartości własne ma-
cierzy trójdiagonalnej aby wybrać optymalny przedział (a, b). Z pozostałych metod
iteracyjnych najszybsza jest metoda gradientów sprzężonych, co nie jest zaskocze-
niem.
Program TIPS, który powstał jako główny rezultat tej pracy umożliwia skuteczne ba-
danie różnych metod numerycznych w zastosowaniu do zagadnienia odwracania syme-
trycznej macierzy trójdiagonalnej. Udostępnia do wyboru szereg testów (łącznie z normą
spektralną), które mogą służyć do analizowania i diagnozowania algorytmów zarówno
w warunkach praktycznych, jak i w pobliżu teoretycznych granic stosowalności danego
algorytmu. Możliwości opracowanych specjalnie dla programu TIPS modułów rysowania
wykresów i obliczania wzorów analitycznych znacznie wykraczają poza prezentowane
w tej pracy przykłady. Modułowa struktura środowiska obliczeniowego TIPS pozwala
na rozbudowę programu w przyszłości. A
atwo też połączyć jego funkcjonalność z innymi
środowiskami obliczeniowymi, gdyż wszystkie dane i wyniki zapisane są w przenośnym
formacie tekstowym z dokładnością 17 cyfr znaczących. Wszystkie te cechy powodują, że
program TIPS może być przydatnym narzędziem w analizie algorytmów odwracania ma-
cierzy (po pewnych modyfikacjach nie tylko trójdiagonalnych), a także może służyć jako
pomoc dydaktyczna w prezentacji wielu problemów związanych z analizą numeryczną.
54
Spis ilustracji
2.1. Minor M42 macierzy T6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.1. Wygląd okna głównego programu TIPS . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2. Wygląd zakładki  Macierze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3. Rodzaje macierzy do wyboru . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4. Wzory do przekształcenia macierzy w macierz diagonalnie dominującą . . . . . . . . 30
3.5. Zakładka  Algorytmy i wyniki . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.6. Lista dostępnych do wyboru algorytmów . . . . . . . . . . . . . . . . . . . . . . . . 33
3.7. Okno  Opcje wyników . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.8. Opcje zapisu pliku(-ów) z macierzą trójdiagonalną . . . . . . . . . . . . . . . . . . . 34
3.9. Opcje zapisu pliku z macierzą odwrotną . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.10. Zakładka  Testy i wykresy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.11. Okno  Nowy wykres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.12. Okno  Edytuj wykres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.13. Okno  Wykres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.14. Menu podręczne do okna  Wykres . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.15. Okno  Wykres z linią regresji i okno  Parametry regresji . . . . . . . . . . . . . . 40
3.16. Przykład wykorzystania osi logarytmicznych i linii regresji . . . . . . . . . . . . . . 42
3.17. Okno  Obliczenia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.18. Wstrzymanie obliczeń i pytanie o potwierdzenie chęci ich przerwania . . . . . . . . . 43
4.1. Wskaz
nik uwarunkowania macierzy Poissona . . . . . . . . . . . . . . . . . . . . . . 46
4.2. Błąd lewostronny metody Gaussa dla macierzy Poissona . . . . . . . . . . . . . . . . 47
4.3. Błąd prawostronny metody Gaussa dla macierzy Poissona . . . . . . . . . . . . . . . 47
4.4. Liniowy wzrost błędu LERR metody 2C dla macierzy Poissona . . . . . . . . . . . . 49
4.5. Norma macierzy Hilberta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.6. Błędy PERR Metody 2C dla macierzy Hilberta . . . . . . . . . . . . . . . . . . . . . 50
4.7. Prawostronne błędy komponentowe dla macierzy Hilberta . . . . . . . . . . . . . . . 51
55
Literatura
[1] M. Dryja, J. I M. Jankowscy, Przegląd metod i algorytmów numerycznych, cz. 2, WNT,
Warszawa 1988.
[2] L.L. Trefethen, D. Bau III, Numerical linear algebra, SIAM, Philadelphia 1997.
[3] A. Kiełbasiński, H. Schwetlick, Numeryczna algebra liniowa, WNT, Warszawa 1992.
[4] J. Du Croz, N.J. Higham, Stability of methods for matrix inversion, IMA, J. Numer. Anal.,
12 (1992), pp. 1-19.
[5] B.P. Flannery, W.H. Press, Numerical Recipes in FORTRAN 77. The Art of Scientific
Computing Volume 1, Volume 1 of Fortran Numerical Recipes, Cambridge University Press,
Cambridge 1993.
[6] Hu, G. Y.; O Connell, R. F., Analytical inversion of symmetric tridiagonal matrices, Journal
of Physics A: Mathematical and General, Volume 29, Issue 7, pp. 1511-1513 (1996).
[7] Huang, Y.; McColl, W. F., Analytical inversion of general tridiagonal matrices, Journal of
Physics A: Mathematical and General, Volume 30, Issue 22, pp. 7919-7933 (1997).
[8] Yamani, H. A.; Abdelmonem, M. S., COMMENT: The analytic inversion of any finite
symmetric tridiagonal matrix, Journal of Physics A: Mathematical and General, Volume
30, Issue 8, pp. 2889-2893 (1997).
[9] T. Pang, Metody obliczeniowe w fizyce. Fizyka i komputery, Wydawnictwo Naukowe PWN,
Warszawa 2001.
[10] T. Jurlewicz, Z. Skoczylas, Algebra liniowa 1. Definicje, twierdzenia, wzory, Wydanie piąte
powiększone, Oficyna Wydawnicza GiS, Wrocław 1998.
[11] W. Salejda, M.H. Tyc, M. Just, Algebraiczne metody rozwiązywania równania Schr�dingera,
Wydawnictwo Naukowe PWN, Warszawa 2002.
[12] N.J. Higham, Bounding the error in gaussian elimination for tridiagonal systems, SIAM, J.
Matrix Anal. Appl., Vol. 11, No. 4, pp. 521-530, (1990).
[13] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A.
Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK User s Guide, Third
Edition, SIAM, Philadelphia 1999.
[14] B. Baron, Algorytmy numeryczne w Delphi, Helion, Gliwice 2006.
[15] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R.
Pozo, C. Romine, H. Van der Vorst, Templates for the Solution of Linear Systems: Building
Blocks for Iterative Methods, Second Edition, SIAM, Philadelphia 1994.
56