« 4. Zagadnienie własne IV  Spis treści  6. Metody iteracyjne Kryłowa »  

Zagadnienia

5. Proste metody iteracyjne rozwiązywania układów równań liniowych

Zajmiemy się przybliżonymi metodami rozwiązywania układów równań liniowych

A⁢x=b

z nieosobliwą rzeczywistą macierzą A. Jak już sugerowaliśmy we wstępie, wraz z coraz większymi modelami pojawiającymi się w praktyce obliczeniowej, częściej zachodzi potrzeba rozwiązywania zadań algebry liniowej o bardzo wielkiej liczbie niewiadomych. Gdy A jest dużego rozmiaru, metody bezpośrednie (wykorzystujące na przykład rozkład LU metodą eliminacji Gaussa lub inne rozkłady macierzy, np. QR) mogą być zbyt kosztowne lub mniej wygodne w użyciu. Przykładowo, jeśli niewiadomych jest kilka milionów, a macierz nie ma żadnej specjalnej struktury, to korzystając z eliminacji Gaussa musielibyśmy wykonać rzędu 1019 działań arytmetycznych (i zarezerwować rzędu 1013 bajtów — czyli 10000 GB — pamięci na czynniki rozkładu). Zakładając optymistycznie, że nasz komputer jest w stanie wykonać 1010 działań na sekundę (czyli jest w stanie realnie osiągnąć wydajność 10 Gflopów/s), dawałoby to 109 sekund, czyli około 32 lat… Rozsądną alternatywą może być wtedy rozwiązanie układu w sposób przybliżony, ale za to znacznie tańszy lub wygodniejszy. Jednym ze sposobów przybliżonego rozwiązywania układów równań liniowych są metody iteracyjne, szczególnie użyteczne wówczas, gdy A jest macierzą rozrzedzoną.

5.1. Macierze rozrzedzone

Macierze, w których bardzo wiele elementów jest zerowych, nazywamy macierzami rozrzedzonymi lub, potocznie, rzadkimi. Dla odróżnienia, macierze, które nie są rzadkie, nazwiemy gęstymi — przykładem takiej macierzy jest macierz, której wszystkie elementy są niezerowe. Intuicyjnie możemy spodziewać się, że praktyczne zadania liniowe wielkiego wymiaru będą prowadziły właśnie do macierzy rozrzedzonej, gdyż związki pomiędzy niewiadomymi w pojedynczym równaniu zadanego układu nie będą raczej dotyczyć wszystkich, tylko wybranej (nielicznej) grupy innych niewiadomych.

Wykorzystanie rozrzedzenia macierzy nie tylko może doprowadzić do algorytmów istotnie szybszych od ich analogonów dla macierzy gęstych (jak pamiętamy, standardowe algorytmy oparte na rozkładzie macierzy gęstej, np. LU, potrzebują O⁢N3 działań arytmetycznych), ale wręcz może być jedynym realnym sposobem na to, by niektóre zadania w ogóle stały się rozwiązywalne przy obecnym stanie techniki obliczeniowej!

Gdy A jest rozrzedzona, mnożenie takiej macierzy przez wektor jest bardzo tanie (koszt jest proporcjonalny do liczby niezerowych elementów macierzy). Dlatego, konstruując metodę przybliżonego rozwiązywania układu, warto oprzeć się na iteracji, których głównym składnikiem jest operacja mnożenia przez A lub jej część.

5.1.1. Przykłady macierzy rozrzedzonych

Jak zobaczymy poniżej, rzeczywiście łatwo można spotkać realne zadania matematyki stosowanej, w których macierz wymiaru N ma tylko O⁢N niezerowych elementów. Omówimy tu dwie klasy takich zadań: dyskretyzacje równań różniczkowych oraz łańcuchy Markowa z dyskretną przestrzenią stanów.

Równania różniczkowe cząstkowe

Jednym ze szczególnie ważnych źródeł układów równań z macierzami rozrzedzonymi są równania różniczkowe cząstkowe (pochodzące np. z modeli pogody, naprężeń w konstrukcji samochodu, przenikania kosmetyków do głębszych warstw skóry, itp.).

Przykład 5.1 (Dyskretyzacja jednowymiarowego laplasjanu)

Rozważmy modelowe eliptyczne równanie różniczkowe z jednorodnym warunkiem brzegowym Dirichleta

-u′′⁢x=f⁢x⁢⁢∀x∈0,1,(5.1)
u⁢0=u⁢1=0,(5.2)

w którym u jest szukanym rozwiązaniem, a f — zadaną funkcją. Aby znaleźć jego przybliżone rozwiązanie, możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym,

-ui-1-2⁢ui+ui+1h2=f⁢xi⁢⁢∀i=1,…,N,(5.3)
u0=uN+1=0,(5.4)

gdzie ui ma odpowiadać wartości rozwiązania w węźle dyskretyzacji xi, ui≈u⁢xi, oraz xi=i⋅h, przy czym h=1/N+1. Odpowiedzią na pytanie, jak dobrze (i czy w ogóle) takie rozwiązanie dyskretne aproksymuje rozwiązanie dokładne, zajmujemy się na wykładzie z Numerycznych równań różniczkowych; tutaj zobaczymy, do jakiego zagadnienia algebraicznego prowadzi równanie różnicowe (5.3).

Jest to oczywiście równanie liniowe na współczynniki UN=u1,…,uNT,

TN⁢⁢UN=h2⁢FN,

gdzie

TN=2-1-12-1⋱⋱⋱-12-1-12N×N(5.5)

oraz FN=f⁢x1,…,f⁢xNT. Ponieważ interesują nas dobre przybliżenia (to znaczy przypadek, gdy h jest bardzo małe), znaczy to, że N będzie bardzo duże. Nasza macierz jest rzeczywiście macierzą rozrzedzoną, bo ma nieco mniej niż 3⁢N niezerowych elementów. Jej współczynnik wypełnienia (ang. density), czyli stosunek liczby niezerowych do liczby wszystkich elementów macierzy, jest rzędu 3/N i maleje ze wzrostem N.

Otrzymana powyżej macierz, aczkolwiek ma bardzo wiele cech charakterystycznych dla dyskretyzacji ,,poważniejszych” równań różniczkowych, wydaje się jednak zbyt trywialna, by traktować ją metodą inną niż eliminacja Gaussa (uwzględniając jej pasmową strukturę): rzeczywiście, jest to macierz trójdiagonalna (a przy tym: symetryczna i dodatnio określona), zatem właściwie zrealizowany algorytm eliminacji Gaussa (lub lepiej: rozkładu L⁢D⁢LT) pozwoli nam wyznaczyć rozwiązanie układu TN⁢UN=h2⁢FN kosztem O⁢N, a więc — z dokładnością do stałej — optymalnym.

Jednak, gdy przejdziemy do wyższych wymiarów i być może dodatkowo regularną siatkę węzłów dyskretyzacji zastąpimy na przykład nieregularną siatką elementu skończonego, uzyskamy macierze o znacznie bardziej skomplikowanej strukturze, które już nie tak łatwo poddają się eliminacji Gaussa.

Przykład 5.2 (Dyskretyzacja wielowymiarowego laplasjanu)

Przez analogię z poprzednim przykładem, rozważmy równanie Poissona w obszarze Ω⊂Rd, z jednorodnym warunkiem Dirichleta:

-Δ⁢u⁢x=f⁢x⁢⁢∀x∈Ω,(5.6)
u|∂⁡Ω=0,(5.7)

w którym u:Ω→R jest szukanym rozwiązaniem, a f:Ω→R — zadaną funkcją. Δ oznacza operator Laplace'a,

Δ≡∂2∂⁡x12+…+∂2∂⁡xd2.

W dalszym ciągu przyjmiemy, dla ominięcia dodatkowych trudności, że Ω jest kostką jednostkową: Ω=0,1d.

Aby znaleźć przybliżone rozwiązanie, znów możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym. Wybierzemy tu najmniej efektywną, ale za to koncepcyjnie najprostszą metodę dyskretyzacji, opartą na równomiernej siatce węzłów w Ω.

Przykładowo, jeśli d=2, to weźmiemy siatkę węzłów pi⁢j=xi,yj∈Ω, gdzie xi=i⁢hx,yj=j⁢hy, przy czym hx=a/nx+1 jest krokiem siatki w kierunku x i analogicznie hy=b/ny+1, zob. rysunek LABEL:regulargrid.pstex_t.

Rys. 5.1. Regularna siatka na Ω.

Zastępując pochodne ilorazami różnicowymi [8]

Δ⁢u⁢xi,yj≈Δh⁢ui⁢j=1h2⁢4⁢ui,j-ui+1,j-ui-1,j-ui,j+1-ui,j-1,(5.8)

dostajemy równanie różnicowe

-ui-1,j-2⁢ui,j+ui+1,jhx2-ui,j-1-2⁢ui,j+ui,j+1hy2=fi⁢j⁢⁢∀i=1,…,nx,⁢j=1,…,ny.(5.9)

W pozostałych węzłach siatki rozwiązanie przyjmuje wartość równą zero na mocy warunku brzegowego Dirichleta. Więcej na temat metod dyskretyzacji takich i podobnych równań różniczkowych można dowiedzieć się na wykładzie Numeryczne równania różniczkowe.

W dalszym ciągu, dla uproszczenia zapisu będziemy zakładać, że nx=ny=n i w konsekwencji hx=hy=h. W takim układzie mamy do wyznaczenia N=n2 niewiadomych:

UN=u1,1,u2,1,u3,1,…,uNx,1︸niewiadome pierwszej kolumny,u1,2,u2,2,u3,2,…,uNx,2︸drugiej kolumny,…,u1,Ny,u2,Ny,u3,Ny,…,uNx,Ny︸ostatniej kolumnyT.

które spełniają liniowy układ równań o postaci macierzowej

PN⁢⁢UN=h2⁢FN,

gdzie PN jest tym razem pięciodiagonalną macierzą o blokowej strukturze,

PN=Tn+2⁢I-I-ITn+2⁢I-I⋱⋱⋱-ITn+2⁢I-I-ITn+2⁢IN×N.(5.10)

W kolejnych diagonalnych blokach PN znajdują się macierze trójdiagonalne Tn+2⁢I, gdzie Tn jest znaną już nam macierzą dyskretyzacji jednowymiarowego laplasjanu. Używając więc symbolu Kroneckera dla iloczynu tensorowego macierzy mamy przy naszej numeracji niewiadomych

PN=I⊗Tn+Tn⊗I.

W przypadku tej macierzy, prosta strategia polegająca na eliminacji Gaussa (z wykorzystaniem dodatkowo faktu, że PN jest pasmowa, o szerokości pasma 2⁢n), dałaby nam szansę rozwiązać to równanie kosztem O⁢n4, co jest, przy dużych wielkościach n, wartością trudną do zaakceptowania.

Prowadząc analogicznie dyskretyzację na jednorodnej siatce trójwymiarowego zadania Poissona na kostce jednostkowej (d=3), dostalibyśmy N=n3 oraz zadanie z macierzą siedmiodiagonalną (o szerokości pasma 2⁢n2):

SN=I⊗I⊗Tn+I⊗Tn⊗I+Tn⊗I⊗I.(5.11)

Nawet dla średnich wartości n, układ ten trudno rozwiązać metodą bezpośrednią, opartą na prostym rozkładzie pasmowej macierzy (tutaj na przykład: rozkładzie Cholesky'ego), zarówno ze względu na koszt obliczeniowy, jak i koszt pamięciowy. Rzeczywiście, koszt rozkładu Cholesky'ego uwzględniającego to, że w czynnikach rozkładu poza pasmem nie pojawią się elementy niezerowe wyniesie O⁢n4 w przypadku dwuwymiarowego równania Poissona i O⁢n7 w przypadku trójwymiarowego równania Poissona (dlaczego?). To oczywiście mniej niż pesymistyczne O⁢N3=O⁢n9, ale gdy n jest duże, może być wystarczająco zachęcające do tego, by poszukać tańszych rozwiązań, dających sensowne rozwiązanie w sensownym czasie.

Przykład 5.3 (ewolucyjne równanie dyfuzji)

Rozważmy równanie

ut⁢x,t-Δx⁢u⁢x,t=f⁢x,tdla ⁢x,t∈Ω×0,T,
u⁢x,t=0dla ⁢x,t∈∂⁡Ω×0,T.

Dla uproszczenia przyjmiemy, że Ω jest kwadratem jednostkowym. Wprowadzając dyskretyzację tego równania po zmiennej przestrzennej jak poprzednio, a po zmiennej czasowej — niejawny schemat Eulera (patrz rozdział o schematach różnicowych Numerycznych równań różniczkowych), otrzymujemy, po uwzględnieniu warunku brzegowego, liniowy układ równań algebraicznych na przybliżenie UNk rozwiązania w chwili tk,

UNk-UNk-1τ+PN⁢UNk=fNk,

gdzie PN jest określone jak w poprzednim przykładzie oraz

fNk=f⁢xi⁢j,tk

i węzły w Ω numerujemy zgodnie z regułą podaną powyżej. Grupując niewiadome po lewej stronie, dostajemy równanie

I+τ⁢PN⁢UNk=UNk-1+τ⁢fNk,

a więc takie, w którym macierz układu jest dalej rozrzedzona, ale — dla dostatecznie małych τ — jest ,,małym zaburzeniem” macierzy jednostkowej.

Gdy siatka nie jest regularna1Zazwyczaj w zaawansowanych metodach dyskretyzacji równań różniczkowych stosuje się właśnie takie siatki, często — adaptacyjnie dopasowane do przebiegu rozwiązania., macierze dyskretyzacji także tracą regularną strukturę (w powyższych przykładach była ona odwzwierciedlona w tym, że niezerowe elementy macierzy układały się wzdłuż nielicznych diagonali). Macierze bez regularnej struktury są jeszcze trudniejsze dla metod bezpośrednich ze względu na rosnące kłopoty w uniknięciu wypełnienia czynników rozkładu macierzy.

Przykład 5.4 (Macierz z kolekcji Boeinga)

Spójrzmy na macierz sztywności dla modelu silnika lotniczego, wygenerowaną swego czasu w zakładach Boeinga i pochodzącą z dyskretyzacji pewnego równania różniczkowego cząstkowego metodą elementu skończonego. Przykład pochodzi z kolekcji Tima Davisa. Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych).

Opis zastępczy
Rys. 5.2. Struktura niezerowych elementów macierzy bcsstk38.

Jej współczynnik wypełnienia (to znaczy, stosunek liczby niezerowych do wszystkich elementów macierzy) wynosi jedynie 0.006, a więc macierz jest bardzo rozrzedzona: w każdym wierszu są średnio tylko 44 niezerowe elementy.

We wspomnianej kolekcji znajdzuje się jeszcze więcej macierzy rozrzedzonych, pochodzących z najrozmaitszych realnych modeli, nie tylko opartych na równaniach różniczkowych. Choć sporo z nich — tak jak tutaj przedstawiona — będzie symetrycznych, nie jest to regułą; inne układy równań ze zgromadzonego zbioru nie muszą mieć żadnej z pożądanych przez numeryka cech: regularnej/pasmowej struktury, symetrii, ani dodatniej określoności. Jedną z nich przedstawiamy poniżej.

Symulacje działania układów elektronicznych

Współczesne układy elektroniczne mogą zawierać miliony części (tranzystorów, rezystorów, kondensatorów), a ich projektowania jest procesem bardzo złożonym i kosztownym. Co więcej, w miarę postępów opracowywania, coraz trudniej o dokonanie w nich poważniejszych zmian projektu. Z tego powodu konieczna jest komputerowa symulacja ich działania. W ogólności prowadzi to do ogromnego układu równań różniczkowo–algebraicznych, ale także bada się modele (prawa Kirchhoffa!), sprowadzające się do rozwiązania układu równań liniowych wielkiego rozmiaru (N na poziomie kilku milionów).

Wizualizacja grafu macierzy Circuit5M o rozmiarze powyżej pięciu milionów

(Źródło: The University of Florida Sparse Matrix Collection)

Duże łańcuchy Markowa

Modele wielostanowych systemów kolejkowych (np. routera obsługującego wiele komputerów) także prowadzą do gigantycznych układów równań z macierzami rozrzedzonymi o specyficznej strukturze. Jeśli taki układ może z danego stanu przejść tylko do niewielu innych, macierz reprezentująca łańcuch Markowa będzie rozrzedzona.

Jednym z zadań, w którym rozważa się bodaj największy ,,naturalnie” pojawiający się łańcuch Markowa jest zadanie wyznaczania tzw. PageRank — wektora wartościującego strony internetowe w wyszukiwarce Google (zob. artykuł ,,Miara ważności” K. Diksa w Delcie 8/2008): macierz ta jest oryginalnie bardzo rozrzedzona, całkowicie nieregularna i niesymetryczna, a jej rozmiar jest równy liczbie indeksowanych stron WWW (a więc mniej więcej rzędu 109).

Macierz jest bardzo rzadka, bez wyraźnie regularnej struktury
Rys. 5.3. Struktura niezerowych elementów macierzy PageRank dla portalu http://www.mimuw.edu.pl.
Ćwiczenie 5.1

Wykaż, że zadanie TN⁢x=b można rozwiązać, korzystając z eleminacji Gaussa bez wyboru elementu głównego kosztem O⁢N.

Ćwiczenie 5.2

Niech macierz A będzie macierzą pasmową rozmiaru N×N, o szerokości pasma 2⁢k:

ai⁢j=0,⁢ gdy ⁢i-j>k.

Wykaż, że rozkład LU (bez wyboru elementu głównego) macierzy A można wyznaczyć kosztem O⁢N⁢k2 działań arytmetycznych. Jak zmieni się odpowiedź, gdy trzeba będzie dokonać wyboru elementu głównego w kolumnie?

Ćwiczenie 5.3

Zbadaj koszt rozkładu LU macierzy PN, w którym wykorzystuje się jedynie informację o tym, że PN jest macierzą pasmową.

Rozwiązanie: 

Mamy N=n2 i szerokość pasma jest równa n. Zatem na mocy poprzedniego zadania, koszt jest rzędu O⁢n4.

5.2. Macierze rozrzedzone: implementacja

Zanim przystąpimy do omawiania metod rozwiązywania układów równań liniowych z macierzami rozrzedzonymi, warto zapoznać się ze sposobami reprezentacji (formatami) macierzy rozrzedzonych. Ponieważ macierze rozrzedzone mają dużo zerowych elementów, ogólną zasadą jest zapamiętywanie tylko tych różnych od zera, w połączeniu z informacją o ich lokalizacji w macierzy. Spośród wielu struktur danych wygodnych dla przechowywania macierzy rozrzedzonych, opisanych m.in. w monografii [13], największą popularnością wśród numeryków cieszą się dwa: format współrzędnych oraz format spakowanych kolumn (lub wierszy)2Mniej popularny format diagonalny spotkamy m.in. w reprezentacji macierzy pasmowych w LAPACKu.. Liczbę niezerowych elementów macierzy A rozmiaru N×M będziemy oznaczali, przez analogię do funkcji MATLABa, nnz⁢A.

Format współrzędnych

Jest to najprostszy sposób reprezentacji macierzy rozrzedzonych. Do zapamiętania macierzy A rozmiaru N×M, liczącej nnz⁢A niezerowych elementów, wykorzystujemy trzy wektory: I, J — oba typu int — oraz V, typu double, wszystkie o długości nnz⁢A, przy czym zachodzi

AI⁢k,J⁢k=V⁢k,⁢∀k=0,…,nnz⁢A-1,

zob. rysunek 5.4.

Schemat rozmieszczenia elementów macierzy w wektorach
Rys. 5.4. Format współrzędnych reprezentacji macierzy rozrzedzonych.
Przykład 5.5 (Macierz jednowymiarowego laplasjanu w formacie współrzędnych)

Zdefiniujemy niezbędne struktury danych i następnie wypełnimy je tak, by reprezentować macierz TN jednowymiarowego laplasjanu (5.5) w formacie współrzędnych. Macierz TN jest kwadratowa rozmiaru N×N, o co najwyżej 3⁢N niezerowych elementach.

Poniższa struktura: 

typedef struct {
int I;
int J;
double V;
} SpElem;
będzie odpowiadała jednemu elementowi reprezentowanej macierzy. Cała macierz TN zmieści się zatem w tablicy 3⁢N elementów typu SpElem, na które przydzielamy miejsce:

SpElem *T;

/* alokujemy miejsce na 3N elementów */
T = (SpElem *) calloc(3*N, sizeof(SpElem));

Zwróćmy uwagę na to, że do przydzielenia pamięci użyliśmy tym razem funkcji calloc(), a nie malloc(), jak zwykle. To dla większego bezpieczeństwa, bo calloc() wypełnia zerami przydzieloną pamięć.

Teraz możemy wypełnić tablicę. Jak widzimy z dalszego ciągu kodu, zmienną T_NNZ inkrementujemy od początkowej wartości 0 za każdym razem, gdy do macierzy T zapisujemy kolejny element — w ten sposób będziemy dokładnie wiedzieli, ile w rzeczywistości (potencjalnie) niezerowych elementów wpisaliśmy do zaalokowanej tablicy.

int nnz = 0; /* bieżąca liczba wpisanych do tablicy elementów */

/* wypełniamy macierz */
for(i = 1; i <= N-1; i++)
{
	/* element diagonalny */
	T[nnz].I = i;
	T[nnz].J = i;
	T[nnz].V = 2;
	nnz++;
	/* element naddiagonalny w i-tym wierszu */
	T[nnz].I = i;
	T[nnz].J = i+1;
	T[nnz].V = -1;
	nnz++;
	/* element poddiagonalny w i-tej kolumnie */
	T[nnz].I = i+1;
	T[nnz].J = i;
	T[nnz].V = -1;
	nnz++;
}
/* ostatni element, (N,N) */
T[nnz].I = N;
T[nnz].J = N;
T[nnz].V = 2;
nnz++;

Nie zawsze taka forma implementacji formatu współrzędnych będzie możliwa: często wybór implementacji będzie narzucony przez biblioteki, jakich zechcemy używać.

Format spakowanych kolumn (lub wierszy)

Format współrzędnych nie narzucał żadnego uporządkowania elementów macierzy — można było je umieszczać w dowolnej kolejności. Z drugiej strony, narzucenie sensownego porządku mogłoby wspomóc realizację wybranych istotnych operacji na macierzy, na przykład, aby wygodnie było realizować działanie (prawostronnego) mnożenia macierzy przez wektor, wygodnie byłoby przechowywać elementy macierzy wierszami. Tak właśnie jest zorganizowany format spakowanych wierszy (CSR, ang. Compressed Sparse Row). Analogicznie jest zdefiniowany format spakowanych kolumn (CSC, Compressed Sparse Column), którym zajmiemy się bliżej.

Podobnie jak w przypadku formatu współrzędnych, macierz w formacie CSC jest przechowywana w postaci trzech wektorów. V jest wektorem typu double o długości nnz⁢A, zawierającym kolejne niezerowe elementy macierzy A wpisywane kolumnami, I jest wektorem typu int, także o długości nnz⁢A, zawierającym numery wierszy macierzy A, odpowiadających elementom z V. Natomiast zamiast tablicy J, jak to było w formacie współrzędnych, mamy krótszy wektor typu int, P, o długości M+1, zawierający na j-1 miejscu indeks pozycji w V, od której w V rozpoczynają się elementy j-tej kolumny macierzy A (zob. rysunek 5.5).

Schemat rozmieszczenia elementów macierzy w wektorach
Rys. 5.5. Format CSC reprezentacji macierzy rozrzedzonych.

W ten sposób, wartości kolejnych elementów j-tej kolumny macierzy A znajdują się w tablicy V na pozycjach o indeksach P[j-1]…P[j]-1. Konsekwentnie, zawsze P[0]=0, a ostatni element, P[M], jest równy nnz⁢A. Dzięki temu, takie operacje, jak na przykład wydobycie jednej kolumny z macierzy A, lub mnożenie przez wektor, są łatwiejsze do implementacji w schludnej, bezwariantowej pętli [13].

Ostatecznie mamy więc zależność

AI⁢k,j=V⁢k,⁢j=1,…,M,⁢k=P⁢j-1,…,P⁢j-1.

Z tego rodzaju formatu reprezentacji macierzy rzadkich (z drobnymi modyfikacjami) korzystają np. pakiety Octave, MATLAB i UMFPACK.

Przykład 5.6 (Procedura mnożenia macierzy w formacie CSC przez wektor)

Mając macierz w formacie CSC, zadaną tablicami I[NNZ], P[M+1], V[NNZ], możemy napisać zgrabną procedurę mnożenia jej przez wektor x[M]. Wynik będziemy zapisywać w wektorze y[N]:

for(i = 0; i < N; i++)
	y[i] = 0.0; /* zerujemy wektor wynikowy */

for(j = 0; j < M; j++)
	for(k = P[j]; k < P[j+1]; k++)
		y[I[k]] += V[k]*x[j];
Ćwiczenie 5.4

Zapisz w pseudokodzie procedurę mnożenia macierzy rozrzedzonej przez wektor, gdy jest ona reprezentowana

  1. w formacie współrzędnych,

  2. w formacie CSR.

5.3. Metody stacjonarne rozwiązywania układów równań liniowych

Teraz zajmiemy się prostymi i historycznie najstarszymi metodami iteracyjnego rozwiązywania układu równań

A⁢x=b,

gdzie nieosobliwa macierz A jest (zapewne wielkiego) rozmiaru N. Są one bodaj najprostsze w analizie i implementacji, ale — jak można się domyślić — w praktyce najmniej efektywne. Z drugiej strony, stanowią one ważny składnik jednej z najszybszych metod rozwiązywania niektórych trudnych układów równań (por. rozdział 8.3).

Rozważane przez nas w tym rozdziale metody opierają się na rozkładzie macierzy A na część ,,łatwo odwracalną”, M, i ,,resztę”, Z=M-A. Dokładniej, równanie A⁢x=b można zapisać jako zadanie punktu stałego

M⁢x=Z⁢x+b,

a jeśli M jest nieosobliwa, to równoważnie:

x=M-1⁢Z⁢x+b.

Ponieważ jest to zadanie znajdowania punktu stałego postaci x=Φ⁢x, można spróbować zastosować doń metodę iteracji prostej Banacha:

xk+1=M-1⁢Z⁢xk+b.(5.12)

Takie metody nazywamy stacjonarnymi metodami iteracyjnymi.

Ćwiczenie 5.5

Wykaż, że każdą stacjonarną metodę iteracyjną można zapisać w wygodniejszej (dlaczego?) do praktycznej implementacji postaci

xk+1=xk+M-1⁢rk,

gdzie rk=b-A⁢xk.

Rozwiązanie: 
xk+1=M-1⁢Z⁢xk+b=M-1⁢M-A⁢xk+b=xk+M-1⁢b-A⁢xk.

Jest to wygodniejsza postać, bo wymaga operowania wektorem residuum rk, na którym zapewne będziemy opierać kryterium stopu iteracji (więc i tak wtedy byśmy go wyznaczali)). Nie wymaga też jawnego Z.

Metody stacjonarne można zapisać w ogólnej postaci

xk+1⁢=⁢B⁢xk⁢+⁢c,(5.13)

gdzie macierz B oraz wektor c są tak dobrane, by punkt stały x* równania (5.13) był rozwiązaniem równania wyjściowego A⁢x=b.

Dla stacjonarnej metody iteracyjnej, B=M-1⁢Z=I-M-1⁢A oraz c=M-1⁢b. Macierz B nazywamy macierzą iteracji, gdyż zachodzi

xk+1-x*⁢=⁢B⁢xk-x*=…=Bk+1⁢x0-x*.

Jasne jest, że z twierdzenia Banacha o punkcie stałym wynika od razu następujący warunek zbieżności:

Stwierdzenie 5.1 (warunek wystarczający zbieżności metody stacjonarnej)

Jeśli B<1, to metoda (5.13) jest zbieżna co najmniej liniowo do x* dla dowolnego x0∈RN. Tutaj ∥⋅∥ oznacza dowolną normę macierzową indukowaną przez normę wektorową.

Dowód

Przy naszych założeniach, mamy do czynienia z kontrakcją Φ⁢x=B⁢x⁢+⁢c, która odwzorowuje kulę K=x∈RN:x-x*≤x0-x* w siebie:

Φ⁢x-x*=Φ⁢x-Φ⁢x*=B⁢x-B⁢x*≤B⁢x-x*<x-x*.

Ponieważ, z założenia, stała Lipschitza dla Φ wynosi B<1, to z twierdzenia Banacha o kontrakcji (por. twierdzenie 9.1) wynika teza stwierdzenia.

∎
Wniosek 5.1

Jeśli I-M-1⁢A<1, to ciąg określony przez stacjonarną metodę iteracyjną jest zbieżny do x*=A-1⁢b oraz

xk+1-x*≤I-M-1⁢A⁢xk-x*.
Dowód

Aby stwierdzić zbieżność, wystarczy w poprzednim stwierdzeniu podstawić B=M-1⁢M-A=I-M-1⁢A oraz x=xk. Oszacowanie błędu jest konsekwencją własności kontrakcji.

∎

Warunek konieczny i dostateczny zbieżności tej iteracji dla dowolnego wektora startowego x0 podaje poniższe twierdzenie. Podkreślmy, że chodzi nam tu o zbieżność metody do rozwiązania układu, gdy startujemy z dowolnego przybliżenia początkowego x0.

Twierdzenie 5.1 (warunek konieczny i dostateczny zbieżności metody stacjonarnej)

Niech B=I-M-1⁢A będzie macierzą iteracji metody (5.13). Ciąg xk określony tą metodą jest zbieżny do rozwiązania x* dla dowolnego x0∈RN wtedy i tylko wtedy, gdy

ρ⁢B=max⁡λ:λ⁢ jest wartością własną ⁢B<1.
Dowód

Na mocy (5.13) otrzymujemy równanie błędu ek=xk-x*,

ek=B⁢ek-1=…=Bk⁢e0.

Konieczność warunku jest oczywista: jeśli metoda jest zbieżna dla każdego x0, to w szczególności dla takiego, że e0 jest równe wektorowi własnemu macierzy B odpowiadającego zadanej (dowolnej) wartości własnej λ. W konsekwencji,

ek=Bk⁢e0=λk⁢e0→0,

a to oznacza, że λ<1.

To, że warunek wymieniany w tezie twierdzenia jest także wystarczający dla zbieżności metody, wynika z faktu (zob. [11]), że dla dowolnego ϵ>0 istnieje norma macierzowa indukowana przez normę wektorową ∥⋅∥ϵ taka, że

Bϵ≤ρ⁢B+ϵ.

Skoro więc z założenia ρ⁢B<1, to można dobrać takie ϵ, by Bϵ<1 i w konsekwencji

ekϵ≤Bϵk⁢e0ϵ,

skąd xk→x*.

∎
Ćwiczenie 5.6

W powyższym dowodzie, jeśli λ jest zespolona, to odpowiadający jej wektor własny też jest zespolony, a więc — formalnie — nie może być uznany za kontrprzykład zbieżności, bo ta ma być ,,dla każdego x0∈RN”, a nie — ,,dla każdego x0∈CN”. Uzupełnij tę lukę pokazując, że jeśli nawet jeśli λ jest zespolona, to gdy λ≥1 można wskazać taki rzeczywisty wektor x0, że zbieżność nie będzie zachodzić.

Zaletą stacjonarnych metod iteracyjnych jest ich prostota powodująca, że są one wyjątkow wdzięczne do szybkiego zaprogramowania. Zobaczymy to na przykładzie kilku klasycznych metod stacjonarnych: Jacobiego, Gaussa–Seidela i SOR. Wszystkie będą bazować na podziale macierzy A na trzy części: diagonalną, ściśle dolną trójkątną i ściśle górną trójkątną:

A=L+D+U,(5.14)

gdzie

L=0a210a31a320⋮⋮⋱0aN⁢1aN⁢2⋯⋯0,⁢D=a11a22⋱⋱aN⁢N,⁢U=0a12a13⋯a1⁢N0a23⋯a2⁢N0⋱⋮0⋮0

5.3.1. Metoda Jacobiego

Biorąc w (5.12) M=D, gdzie D jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy A (zob. (5.14)), otrzymujemy (o ile na przekątnej macierzy A nie mamy zera) metodę iteracyjną

xk+1⁢=D-1⁢b-L+U⁢xk,

zwaną metodą Jacobiego.

Rozpisując ją po współrzędnych, dostajemy układ rozszczepionych równań (numer iteracji wyjątkowo zaznaczamy w postaci górnego indeksu):

xik+1=1ai⁢i⁢bi-∑j≠iai⁢j⁢xjk,

co znaczy dokładnie tyle, że w i-tym równaniu wyjściowego układu przyjmujemy za współrzędne x wartości z poprzedniej iteracji i na tej podstawie wyznaczamy wartość xi.

Widzimy więc, że metoda rzeczywiście jest banalna w implementacji, a dodatkowo jest w pełni równoległa: każdą współrzędną nowego przybliżenia możemy wyznaczyć niezależnie od pozostałych.

Twierdzenie 5.2 (O zbieżności metody Jacobiego)

W metodzie Jacobiego warunek dostateczny zbieżności, D-1⁢L+U∞<1, jest spełniony np. wtedy, gdy macierz A ma dominującą przekątną, tzn. gdy

ai⁢i⁢>⁢∑j≠iai⁢j,⁢∀i=1,…,N.(5.15)
Dowód

Rzeczywiście, ponieważ wyraz i,j macierzy M-1⁢Z wynosi 0 dla i=j oraz ai⁢j/ai⁢i dla i≠j, to

M-1⁢Z∞=max1≤i≤N⁡∑j=1,j≠iNai⁢j/ai⁢i=max1≤i≤N⁡∑j=1Nai⁢j/ai⁢i⁢- 1⁢< 1,

przy czym ostatnia nierówność wynika z warunku diagonalnej dominacji.

∎

Niestety, w wielu ważnych wypadkach metoda Jacobiego, choć zbieżna, będzie zbieżna zbyt wolno, by nas zadowolić.

Przykład 5.7 (Macierz laplasjanu)

Macierz N×N, zwana macierzą jednowymiarowego laplasjanu

TN=2-1-12⋱⋱⋱-1-12

pojawia się w bardzo wielu zastosowaniach (patrz przykład 5.1), również jako sztucznie wprowadzane podzadanie w niektórych algorytmach numerycznych. Ta macierz jest macierzą taśmową, symetryczną i dodatnio określoną, więc układ równań z tą macierzą można bez trudu rozwiązać metodami bezpośrednimi, kosztem O⁢N. Jak jednak za chwilę się przekonamy, układ równań z macierzą TN będzie wyjątkowo trudny dla klasycznej metody stacjonarnej.

Stosując do TN metodę Jacobiego mamy M=2⁢I oraz Z=TN-M. Obliczając normę macierzy iteracji Jacobiego dostajemy M-1⁢Z∞=1, co — jak wiemy z twierdzenia 5.1 — nie rozstrzyga jeszcze o jej zbieżności lub niezbieżności, ale już może stanowić dla nas poważne ostrzeżenie.

Potrzebna będzie bardziej subtelna analiza. Okazuje się, że są znane wzory na wartości własne λj macierzy TN (por. ćwiczenie 5.8):

λj=2⁢1-cos⁡j⁢πN+1=4⁢sin2⁡j⁢π2⁢N+1,

(a więc, 0<λj<4) dla 1≤j≤N. W konsekwencji, wartościami własnymi M-1⁢Z=12⁢Z=12⁢TN-I są liczby μi=12⁢λi-1. Ponieważ -1<μi<1, znaczy to, że metoda Jacobiego jest zbieżna dla macierzy TN.

Z drugiej strony, nie dajmy się zwieść optymizmowi matematyka (,,nie martw się, jest zbieżny…”): nietrudno sprawdzić, że ρ⁢M-1⁢Z=cos⁡π/N+1=1-π22⁢N+12+O⁢N-4, co oznacza, że metoda Jacobiego — choć zbieżna! — dla dużych N staje się zbieżna tak wolno, że w praktyce bezużyteczna: na przykład, gdy N=100, ρ⁢M-1⁢Z≈0.9995, zatem potrzebowalibyśmy około 4600 iteracji, by mieć pewność, że początkowy błąd zredukowaliśmy zaledwie… 10 razy.

Dopiero w rozdziale 8.3 przekonamy się, jak — przez niebanalne wykorzystanie głębszych informacji o samym zadaniu — można wykorzystać metodę stacjonarną do konstrukcji metody iteracyjnej o optymalnym koszcie, którą będziemy mogli stosować także do macierzy dwu- i trójwymiarowego laplasjanu.

Ćwiczenie 5.7

Wykaż, że wartości własne macierzy dwuwymiarowego laplasjanu PN — zob. (5.10) — są postaci

λj⁢k=4⁢sin2⁡j⁢π2⁢nx+1+sin2⁡k⁢π2⁢ny+1,

dla j=1,…,nx,k=1,…,ny. Podobnie pokaż, że, wartości własne macierzy trójwymiarowego laplasjanu SN — zob. (5.11) — są postaci

λj⁢k⁢l=4⁢sin2⁡j⁢π2⁢nx+1+sin2⁡k⁢π2⁢ny+1+sin2⁡l⁢π2⁢nz+1,

dla j=1,…,nx,k=1,…,ny,l=1,…,nz.

Rozwiązanie: 

Niech TN=V⁢Λ⁢VT będzie rozkładem macierzy jednowymiarowego laplasjanu takim, że Λ jest macierzą diagonalną zawierającą wartości własne, a V taka, że VT⁢V=I zawiera wektory własne TN. Ponieważ macierz PN=I⊗TN+TN⊗I to konsekwentnie

PN=⁢….

Stąd wynika teza. Analogicznie postępujemy dla macierzy trójwymiarowego laplasjanu.

Ćwiczenie 5.8

Wykaż, że jeśli b≠0 oraz c≠0, to rzeczywista trójdiagonalna macierz Toeplitza3Macierzą Toeplitza nazywamy macierz o stałych współczynnikach na diagonalach.

T=dbcdb⋱⋱⋱cdbcdN×N(5.16)

ma jednokrotne wartości własne, równe

λk=d+2⁢c⁢α⁢cos⁡k⁢θ,⁢k=1,…,N,(5.17)

gdzie

α=bc,⁢θ=πN+1.

Zwróć uwagę na to, że gdy b⁢c<0, wartości własne są zespolone. Wektor własny vk macierzy T odpowiadający λk ma współrzędne dane wzorem

vkj=α-j⁢sin⁡j⁢k⁢θ,⁢j=1,…,N.(5.18)
Wskazówka: 

Dla X=diag⁡α,α2,…,αN zachodzi

T=d⁢I+c⁢α⁢X-1⁢TN0⁢X=X-1⁢d⁢I+c⁢α⁢T0⁢X,

gdzie

T0=01101⋱⋱⋱10110N×N.

Znaczy to, że T jest podobna do macierzy d⁢I+c⁢α⁢T0 — a więc ma te same wartości własne. Nietrudno sprawdzić wprost (ze wzorów trygonometrycznych), że T0 ma wektory i wartości własne takie, jak przewiduje wzór.

Przykład 5.8

Niech B będzie losową kwadratową macierzą N×N, która ma średnio 10 elementów w wierszu (zatem jej współczynnik wypełnienia wynosi 10/N) i niech

A=B+p⁢I.

Dobierając do B wartość p możemy sterować stopniem diagonalnej dominacji macierzy A i, ostatecznie, zagwarantować sobie, że A jest diagonalnie dominująca.

W poniższym eksperymencie komputerowym możemy naocznie przekonać się, jak stopień dominacji diagonali wpływa na szybkość zbieżności metody Jacobiego. Dla porównania zobaczymy, jak z naszą macierzą radzi sobie metoda bezpośrednia: oparta na rozkładzie LU macierzy A i starająca się wykorzystać w inteligentny sposób jej rozrzedzenie.


Ćwiczenie 5.9

  • Przekonaj się, czy stopień diagonalnej dominacji (wartość parametru p w skrypcie) ma istotny wpływ na szybkość zbieżności.

  • Domyślnie testy prowadzimy dla macierzy rozmiaru 2500. Sprawdź, jak zmieni się czas wykonania skryptu, gdy N będzie jeszcze większe: która metoda: iteracyjna, czy bezpośrednia zagwarantuje nam sensowny wynik w krótszym czasie?

  • Przekonaj się, że gdy rozmiar macierzy jest niewielki, na przykład N=100, nie ma większego sensu stosować metody iteracyjnej. Przy jakim poziomie tolerancji tol metoda iteracyjna staje się konkurencyjna dla bezpośredniej (pod względem czasu działania)?

Ćwiczenie 5.10

  • Uzupełnij powyższy skrypt o wyznaczenie teoretycznej wartości współczynnika redukcji błędu w normie ∥⋅∥∞ i porównaj ją z faktyczną szybkością redukcji błędu.

  • Sprawdź, że symetria A nie ma większego wpływu na szybkość zbieżności.

Przykład 5.9

Zdarzają się macierze — niestety, nie są to sztucznie generowane, akademickie przypadki — które są patologicznie trudne dla stacjonarnej metody iteracyjnej, a bardzo łatwe dla metody bezpośredniej. Jedną taką macierz już znamy: jest to macierz jednowymiarowego laplasjanu TN.

W naszym teście, dla zadania z macierzą A=TN+p⁢I i różnych wartości parametru p≥0 zbadamy szybkość zbieżności metody Jacobiego i porównamy czas jej działania z metodą bezpośrednią (w przypadku macierzy trójdiagonalnej jest ona praktycznie optymalna).


Tym razem zmiana p wyraźnie wpływa na szybkość zbieżności, a dla p=0 zbieżność jest — jak już zdążyliśmy się przekonać także analitycznie — jest bardzo wolna.

Ćwiczenie 5.11

  • Wyjaśnij przyczynę szybkiej zbieżności metody dla dużych p.

  • Domyślnie eksperymenty prowadzimy dla macierzy rozmiaru 25. Sprawdź, jak zmieni się czas wykonania skryptu, gdy N będzie większe: która metoda zagwarantuje nam wynik w krótszym czasie? Jak odpowiedź zależy od p?

  • Jeśli zależałoby nam na obliczeniu dokładnego rozwiązania układu, powinniśmy znacząco zmniejszyć warunek stopu metody, tol, na przykład do poziomu 10-16. Jak teraz wypada porównanie metody Jacobiego i bezpośredniej, gdy p=0?

Ćwiczenie 5.12

Rozważmy macierz

A=1ααα1ααα1.

Wykaż, że

  • A>0⇔-1<2⁢α<2,

  • metoda Jacobiego dla tej macierzy jest zbieżna, wtedy i tylko wtedy, gdy -1<2⁢α<1.

Czy rezultat o zbieżności metody Jacobiego można uogólnić na przypadek, gdy A rozmiaru N×N ma na diagonali same jedynki, a poza nią wszystkie jej elementy są równe α?

5.3.2. Metoda Gaussa–Seidela

Heurystyka tej metody opiera się na zmodyfikowaniu metody Jacobiego tak, by w każdym momencie iteracji korzystać z najbardziej ,,aktualnych” współrzędnych przybliżenia rozwiązania x. Rzeczywiście, przecież wykonując jeden krok metody Jacobiego, czyli rozwiązując kolejno równania skalarne względem xik+1 dla i=1,…,N:

xik+1=1ai⁢i⁢bi-∑j≠iai⁢j⁢xjk,

nietrudno zauważyć, że w części sumy, dla j<i, moglibyśmy odwoływać się — zamiast do ,,starych” xjk — do ,,dokładniejszych”, świeżo wyznaczonych, wartości xjk+1, tzn. ostatecznie wyznaczać

xik+1=1ai⁢i⁢bi-∑j<iai⁢j⁢xjk+1-∑j>iai⁢j⁢xjk.

W języku rozkładu macierzy A=M-Z i iteracji xk+1=M-1⁢Z⁢xk+b mielibyśmy więc M=L+D (dolny trójkąt macierzy A z diagonalą) oraz Z=-U (ściśle górny trójkąt A) i konsekwentnie zapis macierzowy iteracji

xk+1=L+D-1⁢b-U⁢xk.
Twierdzenie 5.3 (O zbieżności metody Gaussa–Seidela)

Jeśli macierz A jest diagonalnie dominująca, to metoda Gaussa–Seidela jest zbieżna do x* dla dowolnego wektora startowego x0.

Inny wariant tej metody dostalibyśmy, biorąc za M górny trójkąt macierzy A.

Uwaga 5.1

Obie metody, Jacobiego i (zwłaszcza) Gaussa–Seidela stosuje się także czasem w prostych algorytmach rozwiązywania układów równań nieliniowych: ich zaletą jest to, że głównym składnikiem iteracji jest rozwiązywanie skalarnego równania nieliniowego na każdym kroku metody.

Metoda Gaussa–Seidela jest w wielu przypadkach rzeczywiście szybciej zbieżna od metody Jacobiego, np. tak jest w przypadku macierzy jednowymiarowego laplasjanu, patrz wniosek 5.2. Wciąż jednak, dodajmy, dla tego zadania jej zbieżność jest zbyt wolna, by ją stosować jako samodzielną metodę.

Przykład 5.10

Kontynuując przykład 5.8, porównamy szybkość zbieżności metody Gaussa–Seidela i Jacobiego na tym samym zadaniu.


Ćwiczenie 5.13

Porównaj szybkość zbieżności metod Jacobiego i Gaussa–Seidela na macierzy TN jednowymiarowego laplasjanu.

Wskazówka: 

Możesz przeprowadzić porównanie, prowadząc dobrze zaplanowane eksperymenty numeryczne. W MATLABie łatwo utworzysz macierz TN poleceniem 

e = ones(N,1);
TN = spdiags([-e, 2*e, -e], [-1,0,1], N, N);
Jeśli interesuje Cię wynik teoretyczny, czytaj dalej.

5.3.3. Metoda SOR

Zbieżność metody Gaussa–Seidela można przyspieszyć, wprowadzając parametr relaksacji ω i kolejne współrzędne nowego przybliżenia xk+1 wyznaczać, kombinując ze sobą poprzednie przybliżenie xik oraz współrzędną nowego przybliżenia x~ik+1, uzyskanego metodą Gaussa–Seidela:

xik+1=1-ω⁢xik+ω⁢x~ik+1.

Gdyby ω=1, dostalibyśmy z powrotem metodę Gaussa–Seidela. Ponieważ zwykle wybiera się ω>1, powstałą metodę nazywa się metodą nadrelaksacji (ang. successive overrelaxation, SOR), która jest na swój sposób szczytowym osiągnięciem wśród metod stacjonarnych. Rozkład macierzy A=M-Z odpowiadający metodzie SOR z parametrem ω jest zadany przez

M=1ω⁢D+L,⁢Z=1ω-1⁢D-U.
Twierdzenie 5.4 (lemat Kahana o dopuszczalnych wartościach parametru relaksacji SOR)

Jeśli A ma niezerową diagonalę, a metoda SOR jest zbieżna, to musi być 0<ω<2.

Dowód

Wystarczy zbadać promień spektralny macierzy iteracji metody SOR, która jest równa M-1⁢Z=D+ω⁢L-1⁢1-ω⁢D-ω⁢U. Ponieważ macierz D+ω⁢L jest macierzą trójkątną o diagonali D, to det⁡D+ω⁢L-1=det⁡D-1 i w konsekwencji

det⁡M-1⁢Z=det⁡D+ω⁢L-1⁢det⁡1-ω⁢D-ω⁢U=det⁡D-1⁢det⁡1-ω⁢D=det⁡1-ω⁢I=1-ωN.

Ze względu na to, że wyznacznik macierzy jest równy produktowi jej wartości własnych, dochodzimy do wniosku, że ρ⁢M-1⁢Z≥1-ω, co kończy dowód.

∎

Jako ciekawostkę odnotujmy poniższe

Twierdzenie 5.5 (Ostrowskiego i Reicha)

Jeśli A=AT>0, to SOR jest zbieżna dla dowolnego ω∈0,2.

Dowód tego twierdzenia można znaleźć na przykład w [15, rozdział 8.3].

Dla klasy macierzy obejmującej niektóre spotykane w praktycznych zastosowaniach, można wskazać brzemienny w konsekwencje związek promienia spektralnego macierzy iteracji metody SOR i metody Jacobiego.

Definicja 5.1

Rzeczywistą macierz kwadratową A rozmiaru N, o niezerowych elementach na diagonali, nazywamy macierzą zgodnie uporządkowaną, jeśli wartości własne macierzy

Jα=α⁢L~+1α⁢U~,

gdzie L~=D-1⁢L, U~=D-1⁢U, są niezależne od α∈C∖0.

Zauważmy, że J-1=-L~-U~ jest macierzą iteracji metody Jacobiego dla macierzy A.

Szczególnym przypadkiem macierzy zgodnie uporządkowanych są macierze postaci

A=D1U2L1D2,(5.19)

w których D1,D2 są nieosobliwymi macierzami diagonalnymi.

Rzeczywiście, dla macierzy postaci (5.19) mamy

A=L+D+U,

gdzie

L=00L10,⁢D=D1D2,⁢U=0U200

i w konsekwencji

Jα=1α⁢D1-1⁢U2α⁢D2-1⁢L1=Iα⁢I⁢D1-1⁢U2D2-1⁢L1︸J1⁢Iα⁢I-1.

Zatem macierz Jα jest podobna do J1 i w konsekwencji ma te same wartości własne, co J1 (niezależnie od α≠0), a więc — jest zgodnie uporządkowana.

Ćwiczenie 5.14

Wykaż, że jeśli A jest macierzą trójdiagonalną (taką, jak na przykład TN) lub blokowo trójdiagonalną (taką, jak na przykład PN) to A można poddać takiej permutacji P wierszy i kolumn, że uzyskana macierz P⁢A⁢PT jest postaci (5.19)4W rzeczywistości, można wykazać więcej, że nawet bez tej permutacji macierze TN i PN są zgodnie uporządkowane..

Rozwiązanie: 

Rzeczywiście, w przypadku macierzy trójdiagonalnej TN wystarczy wziąć najpierw zmienne o indeksach nieparzystych, a potem — o indeksach parzystych: P⁢1,2,…,NT=1,3,…,2,4,…T. Dla macierzy PN wybieramy podobnie (tzw. red-black ordering).

Twierdzenie 5.6 (Vargi o związku pomiędzy zbieżnością SOR i metody Jacobiego)

Niech A będzie macierzą zgodnie uporządkowaną i niech ω>0. Wtedy jeśli μ jest wartością własną macierzy iteracji metody Jacobiego, to -μ też nią jest. Ponadto, jeśli λ jest różną od zera wartością własną macierzy iteracji metody SOR z parametrem relaksacji ω oraz zachodzi

λ+ω-12=λ⁢⁢ω⁢⁢μ2

to μ jest wartością własną macierzy iteracji Jacobiego. Na odwrót, jeśli μ jest wartością własną macierzy iteracji Jacobiego, to λ zadane powyższym wzorem jest wartością własną macierzy iteracji SOR.

Dowód

Macierz iteracji metody Jacobiego to po prostu

-L~+U~=J⁢-1=J⁢1,

(ostatnia równość na mocy założenia). Stąd oczywiście wynika część pierwsza tezy.

Macierz iteracji metody SOR ma postać

Gω=D+ω⁢L-1⁢1-ω⁢D-ω⁢U=I-ω⁢L~-1⁢1-ω⁢I+ω⁢U~.

Niech teraz λ będzie wartością własną macierzy iteracji metody SOR,

det⁡λ⁢I-Gω=0.

Ponieważ det⁡I-ω⁢L~=1, mamy równoważnie

0=det((I-ωL~)(λI-Gω))=det(λ(I-ωL~)-(1-ω)I-ωU~)
=det⁡λ+ω-1⁢I-λ⁢ω⁢L~-ω⁢U~.

Zatem λ jest wartością własną Gω wtedy i tylko wtedy, gdy λ+ω-1ω jest wartością własną macierzy λ⁢L~+U~.

W konsekwencji, jeśli tylko λ≠0, to λ+ω-1ω⁢λ jest wartością własną macierzy J⁢λ. Na mocy założenia, J⁢λ ma te same wartości własne, co J⁢-1, co kończy dowód drugiej części twierdzenia. Dowód części trzeciej pozostawiamy jako ćwiczenie.

∎
Wniosek 5.2 (O zbieżności metody Gaussa–Seidela)

Jeśli A jest zgodnie uporządkowana i promień spektralny macierzy iteracji metody Jacobiego dla A jest równy μ<1, to promień spektralny macierzy iteracji metody Gaussa–Seidela jest równy μ2. Znaczy to, że metoda Gaussa–Seidela jest dwukrotnie szybsza od metody Jacobiego.

Dowód

Rzeczywiście, na mocy twierdzenia Vargi, każdej niezerowej wartości własnej macierzy SOR odpowiada wartość własna macierzy iteracji Jacobiego μ, spełniająca równanie

λ+ω-12=λ⁢ω⁢μ2.

Dla metody Gaussa–Seidela, ω=1, i w konsekwencji λ=μ2.

∎

Korzystając z twierdzenia Vargi można udowodnić

Twierdzenie 5.7 (Twierdzenie Younga o optymalnym parametrze SOR)

Niech macierz A będzie symetryczną macierzą zgodnie uporządkowaną i niech μ<1 będzie równe promieniowi spektralnemu macierzy iteracji metody Jacobiego. Wtedy optymalna wartość parametru relaksacji ω, gwarantująca najmniejszy promień spektralny macierzy iteracji metody SOR, jest równa

ωo⁢p⁢t=21+1-μ2.

Promień spektralny macierzy iteracji SOR z tym parametrem wynosi wtedy

μ1+1-μ22.

Dowód powyższego twierdzenia można znaleźć na przykład w [15, rozdział 8.3].

Jakkolwiek piękne teoretycznie, nawet twierdzenia Younga i Vargi są dość bezlitosne dla metody SOR w wielu praktycznych zastosowaniach.

Przykład 5.11

Dla macierzy jednowymiarowego laplasjanu rozmiaru N, oznaczając ξ=π/N+1, mamy μ=cos⁡ξ i w konsekwencji ωo⁢p⁢t=2/1+sin⁡ξ. Zatem, dla N=100, gdy μ≈0.9995, promień spektralny macierzy iteracji SOR z optymalnym parametrem relaksacji jest równy aż 1-sin⁡ξ/1+sin⁡ξ≈0.9396, a więc jest praktycznie tak samo blisko jedności, jak w przypadku metody Jacobiego!

Przykład 5.12

Kontynuując przykład 5.9, porównamy szybkość zbieżności metody SOR, Gaussa–Seidela (czyli SOR z parametrem ω=1) i Jacobiego na tym samym, ekstremalnie trudnym dla metod iteracyjnych zadaniu. Na początek weźmiemy N=20 i zbadamy zależność szybkości zbieżności SOR od wartości ω. Korzystając z wyniku poprzedniego przykładu, będziemy mogli także uruchomić SOR z optymalnym parametrem.


Zwróć uwagę na dramatyczną przewagę SOR z optymalnym parametrem nad metodą Gaussa–Seidela. Zbadaj, jak zmienią się wyniki, gdy wyraźnie zwiększysz N, np. do tysiąca.

Ćwiczenie 5.15

Gdy macierz iteracji B metody stacjonarnej

xk+1=B⁢xk+d

jest daleka od normalnej, numeryczna realizacja iteracji w arytmetyce zmiennoprzecinkowej ograniczonej precyzji może napotkać na problemy: w początkowych iteracjach residuum może znacząco rosnąć, co w efekcie może doprowadzić do praktycznej utraty zbieżności.

Zbadaj to na przykładzie–zabawce (por. [2, 10.2.3]),

B=λ10λ,

gdy λ∈0.51.0. Wyznacz promień spektralny tej macierzy, jej normę, a także wykres normy xk w zależności od k.

Rozwiązanie: 

Dla λ≈1 (wyraźnie większych od 1), macierz iteracji B, choć o promieniu spektralnym mniejszym niż 1, ma tę właściwość, że norma jej kolejnych potęg początkowo rośnie. 

lam = 0.9995; B = [lam 1; 0 lam];
resid=[]; x = [1;1];
for i = 1:100000,
	x = B*x; resid(i) = norm(x);
end;
semilogy(resid);

Ćwiczenie 5.16

Metoda Richardsona z parametrem τ∈R jest określona wzorem

xk+1=xk+τ⁢b-A⁢xk.

Niech macierz A będzie symetryczna.

  1. Sprawdź, że jest to stacjonarna metoda iteracyjna oparta na rozszczepieniu A=M-Z, gdzie M=1τ⁢⁢I.

  2. Wykaż, że jeśli A jest nieokreślona, dla żadnego wyboru τ nie jest prawdą, że metoda Richardsona jest zbieżna do x* dla dowolnego x0.

  3. Wykaż, że jeśli A jest dodatnio określona, to metoda Richardsona jest zbieżna dla 0<τ<2/λmax⁢A, a dla

    τopt=2λmax⁢A+λmin⁢A,

    promień spektralny macierzy iteracji M-1⁢Z jest najmniejszy możliwy i równy

    ρ⁢M-1⁢Z=M-1⁢Z2=κ-1κ+1,⁢ gdzie ⁢κ=λmax⁢Aλmin⁢A.
Ćwiczenie 5.17

Niech Asym=A+AT2 i niech istnieją stałe 0<λ≤Λ takie, że dla każdego x∈RN

λ⁢xT⁢x≤xT⁢Asym⁢x⁢⁢oraz⁢⁢A⁢x22≤Λ⁢xT⁢Asym⁢x.

Wykaż, że wtedy metoda Richardsona jest zbieżna w normie euklidesowej dla 0<τ<2/Λ. Ponadto, jeśli τ=1/Λ, to współczynnik zbieżności ρ w tej normie można oszacować przez 1-λΛ.

Ćwiczenie 5.18 (Hackbusch)

Niech A,M będą macierzami symetrycznymi i dodatnio określonymi. Dla dwóch macierzy symetrycznych X,Y będziemy pisali X>Y, jeśli X-Y>0, czyli, innymi słowy, gdy xT⁢X⁢x⁢>⁢xT⁢Y⁢x dla każdego x≠0. Wykaż, że:

  1. Jeśli

    2⁢M>A,

    to iteracja

    xk+1=xk+M-1⁢b-A⁢xk

    jest zbieżna, a dokładniej, macierz iteracji, B=I-M-1⁢A, spełnia

    ρ⁢B=BA=BM<1.

  2. Jeśli dla pewnych 0<λ≤Λ zachodzi

    λ⁢W≤A≤Λ⁢W,

    to wszystkie wartości własne B=I-M-1⁢A są rzeczywiste oraz 1-Λ≤λ⁢B≤1-λ.

5.4. Metody blokowe

Pewnym remedium na powolną zbieżność metod stacjonarnych jest stosowanie metod blokowych, w których iterację konstruuje się tak jak dotychczas, ale zamiast na elementach, działa się na całych blokach macierzy (naturalnie, jeśli trzeba, musimy założyć ich odwracalność). Dla poprawienia zbieżności często stosuje się wariant z zakładką, zob. rysunek 5.6. W praktyce takie metody mogą być wyraźnie skuteczniejsze od dotychczas rozważanych: rzeczywiście, w skrajnym przypadku, pojedynczego ,,bloku” — rozmiaru N — wszystkie wyżej omówione metody stacjonarne redukują się do metody bezpośredniej…

Macierz A jest podzielona na kilka, częściowo nakładających się na siebie, bloków. Te bloki następnie są --- każdy z osobna --- odwracane.
Rys. 5.6. Schemat ideowy blokowej metody Jacobiego z zakładką: zamiast standardowego D-1⁢r=∑i=1Nei⁢di⁢i-1⁢eiT⁢r wyznaczamy ∑i=1KRi⁢Ai⁢i-1⁢RiT⁢r, gdzie Ri jest macierzą operatora obcięcia r do współrzędnych odpowiadających Ai⁢i.

Ciekawym i nietrywialnym uogólnieniem metod z zakładką jest addytywna metoda Schwarza: jeden z przedstawicieli klasy tzw. metod dekompozycji obszaru, będących bardzo efektywnymi i naturalnie równoległymi operatorami ściskającymi (mówimy o nich w rozdziale 7) dla dyskretyzacji równań różniczkowych cząstkowych, [14].

5.5. Metody projekcji

Dotychczas omówione metody opierają się na czysto formalnym, mechanicznym, podziale macierzy, bezpośrednio związanym z jej strukturą. Tymczasem można zdefiniować metody iteracyjne o podobnym koszcie iteracji, ale w których kolejne przybliżenia rozwiązania będziemy wybierać tak, by w jakimś sensie zminimalizować błąd.

Metodę iteracyjną będziemy definiować w formie aktualizacji poprzedniego przybliżenia,

xk+1=xk+δk.

,,Idealna poprawka” δk, oznaczmy ją δ*, powinna spełniać więc równanie

A⁢δ*=rk,

bo wtedy dostalibyśmy xk+1=x*. Ma ona jednak tę wadę, że do jej wyznaczenia musielibyśmy dokładnie rozwiązać nasz układ równań (co prawda z nieco inną prawą stroną, ale to bez znaczenia!) — a więc wpadlibyśmy z deszczu pod rynnę! Jednak, gdyby można było tanio rozwiązać ten układ równań w przybliżeniu, to może w ten sposób udałoby się nam określić zbieżną metodę?

Aby określić przybliżone rozwiązanie równania poprawki, najpierw przeformułujemy je nieco. Niech dwie macierze: U i V mają tę własność, że kolumny każdej z nich tworzą bazę RN. Wtedy rozwiązanie dokładne równania idealnej poprawki, δ*, możemy reprezentować w bazie V, δ*=V⁢a* dla pewnego a*∈RN. Co więcej, z nieosobliwości macierzy U i V wynika, że

A⁢δ*=rk⇔UT⁢A⁢V⁢a*=UT⁢rk.

To daje nam pomysł na definicję przybliżonej poprawki idealnej: niech Uk i Vk będą macierzami pełnego rzędu, tego samego zadanego z góry rozmiaru N×dk. Wtedy poprawkę δk określimy jako wektor

δk=Vk⁢ak

taki, że ak∈Rdk spełnia równanie

UkT⁢A⁢Vk⁢ak=UkT⁢rk.(5.20)

Jest to właśnie metoda projekcji.

Ponieważ macierz zredukowana Ak=UkT⁢A⁢Vk jest kwadratowa rozmiaru dk, to δk będzie tanie do wyznaczenia, gdy dk jest niewielkie. Nazwa metody bierze się stąd, że z powyższej definicji wynika, że wektor residualny równania poprawki, rk-A⁢δk jest prostopadły do podprzestrzeni rozpiętej przez kolumny Uk:

UkT⁢rk-A⁢δk=0.

Z drugiej strony, δk będzie dobrze określone tylko wtedy, gdy macierz zredukowana Ak będzie nieosobliwa — co nie zawsze musi być prawdą, nawet jeśli A jest nieosobliwa (pomyśl o macierzy 0111 i Uk=Vk=10). Aby zagwarantować sobie odwracalność macierzy zredukowanej, zwykle wybiera się macierz pełnego rzędu Vk i w zależności od własności macierzy A dobiera się macierz Uk:

  • Jeśli A=AT>0, to kładziemy Uk=Vk. Wtedy macierz zredukowana Ak=VkT⁢A⁢Vk jest symetryczna i dodatnio określona.

  • Jeśli A jest tylko nieosobliwa, to kładziemy Uk=A⁢Vk. Macierz zredukowana jest postaci Ak=VkT⁢AT⁢A⁢Vk, a więc jest macierzą zredukowaną poprzedniego rodzaju, ale dla macierzy równań normalnych, AT⁢A.

Jak możemy się domyślić, metody projekcji są metodami minimalizacji, co potwierdza poniższy dwuczęściowy lemat:

Lemat 5.1

Niech Vk będzie zadaną macierzą N×dk, pełnego rzędu. Oznaczmy przez Vk podprzestrzeń rozpiętą przez kolumny Vk.

  • Jeśli A=AT>0 oraz Uk=Vk, to xk+1 określone metodą projekcji (5.20) spełnia

    x*-xk+1A≤x*-xA⁢⁢∀x∈xk+Vk.

  • Jeśli A jest nieosobliwa i Uk=A⁢Vk, to xk+1 określone metodą projekcji (5.20) spełnia

    b-A⁢xk+12≤b-A⁢x2⁢⁢∀x∈xk+Vk.

Dowód zostawiamy jako ćwiczenie.

Ćwiczenie 5.19

Udowodnij lemat 5.1.

Wskazówka: 

Patrz [13] lub dowód twierdzenia 6.1.

5.5.1. Metoda najszybszego spadku

Jednym z bardziej prominentnych przykładów metody projekcji jest metoda najszybszego spadku, działająca w przypadku, gdy macierz A jest symetryczna i dodatnio określona. W tej metodzie wybieramy Uk=Vk=rk=b-A⁢xk. Ponieważ wymiar przestrzeni rozpiętej przez kolumny Uk jest równy dk=1, równanie poprawki upraszcza się do jednego równania skalarnego na ak∈R,

ak=rkT⁢rkrkT⁢A⁢rk

i w konsekwencji xk+1=xk+ak⁢rk.

Nazwa metody wywodzi się stąd, że wektor poprawki w tej metodzie jest proporcjonalny do residuum, które z kolei jest kierunkiem gradientu funkcjonału ϕ⁢x=x*-xA2 w punkcie xk:

∇⁡ϕ⁢xk=b-A⁢xk=rk.
Twierdzenie 5.8 (o zbieżności metody najszybszego spadku)

W metodzie najszybszego spadku,

x*-xk+1A≤κ-1κ+1⁢x*-xkA,

gdzie κ=cond2⁡A=λmax⁢A/λmin⁢A.

Dowód

Łatwo wykazać (por. [13, twierdzenie 5.2]), że jeśli rk≠0,

x*-xk+1A2≤1-rkT⁢rkrkT⁢A⁢rk⁢rkT⁢rkrkT⁢A-1⁢rk⁢x*-xkA2.

Teza wynika z lematu, którego elegancki dowód, pochodzący od Braessa, można znaleźć w [2]:

Lemat 5.2 (Kantorowicza)

Niech A=AT>0. Wtedy dla dowolnego x≠0,

xT⁢A⁢xxT⁢x⋅xT⁢A-1⁢xxT⁢x≤14⁢κ+1κ2.
∎
Ćwiczenie 5.20

Wykaż, że w metodzie najszybszego spadku zachodzi rk+1T⁢rk=0.

Rozwiązanie: 

Ponieważ xk+1=xk+ak⁢rk, to

rk+1=rk-ak⁢A⁢rk,

zatem

rkT⁢rk+1=rkT⁢rk-ak⁢rkT⁢A⁢rk=rk22-ak⁢rkT⁢A⁢rk=0

z definicji ak.

Przykład 5.13

Kontynuujemy przykład 5.8. Chcąc porównywać trzy metody: Jacobiego, Gaussa–Seidela oraz metodę najszybszego spadku, musimy zadbać o to, żeby były spełnione warunki zbieżności tej ostatniej — a więc, aby macierz A była symetryczna i dodatnio określona. Dlatego, tym razem położymy, dla dodatniego p,

A=BT⁢B+p⁢I.

Choć A wydaje się dosyć gęsta, w rzeczywistości wciąż jest macierzą rzadką. Zwróćmy uwagę nie tylko na samą szybkość zbieżności — mierzoną liczbą iteracji — ale także na efektywność metody: ile czasu zajmuje wyznaczenie przybliżenia z zadaną dokładnością. Metoda najszybszego spadku niewątpliwie ma najtańszą iterację (a najdroższą — metoda Gaussa–Seidela), co ostatecznie przekłada się na większą efektywność metody najszybszego spadku (pomimo mniejszej szybkości zbieżności).

Ćwiczenie 5.21

Wyjaśnij, dlaczego w powyższym przykładzie, dla małych wartości p, metoda Jacobiego czasem nie jest zbieżna.

Rozwiązanie: 

Bo macierz A nie musi być diagonalnie dominująca.

Przykład 5.14

Dla macierzy jednowymiarowego laplasjanu, mamy

κ=λmax⁢Aλmin⁢A=sin⁡π⁢N/2⁢N+1sin⁡π/2⁢N+12=O⁢N2,

zatem dla N=100 mamy współczynnik redukcji błędu na poziomie κ-1/κ+1≈0.9995. Znaczy to dokładnie tyle, że metoda nie nadaje się do rozwiązywania takich zadań, gdy N jest duże. Szybkość zbieżności metody najszybszego spadku jest w naszym przykładzie porównywalna z szybkością metody Jacobiego i gorsza od metody SOR z optymalnym parametrem relaksacji.

5.5.2. Metoda najmniejszego residuum

Gdy o macierzy A wiemy jedynie, że jest nieosobliwa, możemy zastosować metodę najmniejszego residuum. W tej metodzie wybieramy Vk=rk=b-A⁢xk oraz Uk=A⁢Vk=A⁢rk. Równanie poprawki znów upraszcza się do jednego równania skalarnego na ak∈R,

ak=rkT⁢AT⁢rkrkT⁢AT⁢A⁢rk

i w konsekwencji xk+1=xk+ak⁢rk.

Twierdzenie 5.9

Załóżmy, że macierz Asym=A+AT/2 jest dodatnio określona i oznaczmy μ=λmin⁢Asym>0. Wtedy

rk+12≤1-μ2A221/2⁢rk2.

Asym nazywana jest częścią symetryczną macierzy A.

Ćwiczenie 5.22

Przeprowadź dowód twierdzenia o zbieżności metody najmniejszego residuum.

Wskazówka: 

Zob. dowód twierdzenia 5.3 w [13].

Ćwiczenie 5.23

Przypuśćmy, że umiemy tanio rozwiązywać układy równań z (nieosobliwą) macierzą A rozmiaru N×N. Niech będzie dana macierz

A~=AvuTδ,

gdzie u,v∈RN oraz δ∈R.

  • Wskaż warunek konieczny i dostateczny na to, by macierz A~ była nieosobliwa.

  • Podaj możliwie tani algorytm rozwiązywania układu równań z macierzą A~.

Wskazówka: 

A~ jest nieosobliwa wtedy i tylko wtedy, gdy δ-uT⁢A-1⁢v≠0.

Ćwiczenie 5.24 (wzór Shermana–Morrisona)

Przypuśćmy, że umiemy tanio rozwiązywać układy równań z (nieosobliwą) macierzą A rozmiaru N×N. Niech będzie dana macierz

A~=A-v⁢uT.

gdzie u,v∈RN.

  • Wskaż warunek konieczny i dostateczny na to, by macierz A~ była nieosobliwa.

  • Podaj możliwie tani algorytm rozwiązywania układu równań z macierzą A~.

Wskazówka: 

Zauważ, że A~⁢x=b wtedy i tylko wtedy, gdy xt jest rozwiązaniem układu

AvuT1⁢xt=b0.

Patrz także stwierdzenie 11.1.

Quiz 5.1

  1. diagonalna


    Dobrze Oczywiście, macierz diagonalną możemy rozwiązać bezpośrednio kosztem N flopów — a więc optymalnym. Wiele metod jednak sprowadziłoby się do metody bezpośredniej — na przykład metoda Jacobiego…

    Źle Zastanów się, ile kosztowałoby rozwiązanie takiego układu równań metodą bezpośrednią, wykorzystującą specjalną postać macierzy

  2. trójdiagonalna


    Dobrze Rozkład LU macierzy trójdiagonalnej możemy wyznaczyć bezpośrednio kosztem O⁢N flopów, podobnie z rozwiązaniem samego układu (macierz L jest dwudiagonalna, a macierz U — co najwyżej trójdiagonalna). Stałe przy tym nie są zbyt duże

    Źle Zastanów się, ile kosztowałoby rozwiązanie takiego układu równań metodą bezpośrednią, wykorzystującą specjalną postać macierzy

  3. pełna, tzn. bez zerowych elementów


    Dobrze Ponieważ macierz jest pełna, to jedno mnożenie wektora przez A kosztuje O⁢N2 flopów, więc metoda iteracyjna będzie miała sens, gdy satysfakcjonujące przybliżenie dostaniemy po k≤N iteracjach

    Źle Rozkład LU macierzy pełnej możemy wyznaczyć bezpośrednio kosztem O⁢N3 flopów, co zwykle jest ponad siły współczesnego komputera, gdy N jest bardzo duże

Ćwiczenie 5.25

Czy metodę iteracyjną warto stosować do rozwiązywania układu A⁢x=b z macierzą N×N, w przypadku, gdy N jest bardzo duże oraz A jest rzadka?

Rozwiązanie: 

Tutaj, jak wynika z naszych dotychczasowych rozważań, odpowiedź jest zniuansowana. Jedno mnożenie wektora przez A kosztuje zapewne O⁢N flopów, więc metoda iteracyjna będzie miała sens, gdy satysfakcjonujące przybliżenie dostaniemy po k≪N iteracjach (a to, jak wiemy, nie zawsze jest gwarantowane). Z drugiej strony, może istnieć skuteczny reordering macierzy (czyli zmiana uporządkowania niewiadomych lub równań), pozwalający tanio wyznaczyć rozwiązanie metodą bezpośrednią [13].

« 4. Zagadnienie własne IV  Spis treści  6. Metody iteracyjne Kryłowa »  

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

strona główna | webmaster | o portalu | pomoc

© Wydział Matematyki, Informatyki i Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.