|
|
| Linia 1: |
Linia 1: |
|
| |
|
| =Układy równań liniowych=
| |
|
| |
| Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń
| |
| zmiennoprzecinkowych, przy czym wiele znanych, <strong>matematycznie równoważnych</strong> metod
| |
| rozwiązywania takich zadań, ma <strong>diametralnie różne własności numeryczne</strong>.
| |
| Bardzo ważną klasą takich zadań jest rozwiązywanie układów równań liniowych
| |
|
| |
| <center><math>\displaystyle
| |
| Ax = b,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle A</math> jest nieosobliwą macierzą <math>\displaystyle N\times N</math>, a dany wektor prawej strony <math>\displaystyle b\in R^N</math>.
| |
|
| |
| W
| |
| praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
| |
| specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
| |
|
| |
| Rozwiązywanie układów równań liniowych jest sercem wielu innych algorytmów
| |
| numerycznych, dlatego nie dziwi, że szacuje się, że około 75 procent czasu
| |
| obliczeniowego superkomputerów jest wykorzystywanych właśnie na rozwiązywanie
| |
| takich zadań.
| |
|
| |
| Okazuje się, że kilka znanych w matematyce sposobów rozwiązywania układów równań
| |
| liniowych, takie jak:
| |
| * metoda wyznacznikowa (wzory Cramera)
| |
| * obliczenie macierzy <math>\displaystyle A^{-1}</math> i następnie <math>\displaystyle x = A^{-1}b</math>
| |
|
| |
| <strong>nie nadają się</strong> do numerycznego rozwiązywania takich zadań.
| |
|
| |
| O tym, jak <strong>skutecznie</strong> rozwiązywać takie zadania, jakie jest ich
| |
| uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z
| |
| dobrych metod jest metoda eliminacji Gaussa.
| |
|
| |
| ==Proste układy równań==
| |
|
| |
| Niektóre układy równań można bardzo łatwo rozwiązać. Zgodnie z zasadą
| |
|
| |
| <blockquote style="background-color:#fefeee">
| |
| Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
| |
| </blockquote>
| |
|
| |
| w dalszej kolejności pokażemy, jak dowolny układ równań sprowadzić do sekwencji
| |
| dwóch (czasem trzech) łatwych do rozwiązania układów równań. Ale... jakie układy
| |
| równań są "łatwe"?
| |
|
| |
| ====Układy z macierzą trójkątną====
| |
|
| |
| Rozważmy układ z macierzą
| |
| trójkątną <math>\displaystyle A</math>. Będą nas szczególnie interesować macierze
| |
| <strong>trójkątne górne</strong>, dla których <math>\displaystyle a_{i,j}=0</math> gdy <math>\displaystyle i>j</math>, oraz
| |
| macierze <strong>trójkątne dolne</strong> z jedynkami na przekątnej, tzn.
| |
| <math>\displaystyle a_{i,j}=0</math>, <math>\displaystyle i<j</math>, oraz <math>\displaystyle a_{i,i}=1</math>. Macierze pierwszego rodzaju
| |
| będziemy oznaczać przez <math>\displaystyle U</math>, a drugiego rodzaju przez <math>\displaystyle L</math>.
| |
|
| |
| <center><math>\displaystyle L = \begin{pmatrix}
| |
| 1 & & & & & \\
| |
| * & 1 & & & & \\
| |
| * & * & 1 & & & \\
| |
| * & * & * & 1 & & \\
| |
| \vdots & \vdots & \vdots & \ddots & \ddots & \\
| |
| * & * & * & \cdots & * & 1
| |
| \end{pmatrix} ,
| |
| \qquad
| |
| U = \begin{pmatrix}
| |
| * & * & * & * & \cdots & * \\
| |
| & * & * & * & \cdots & * \\
| |
| & & * & * & \cdots & * \\
| |
| & & & * & \ddots & \vdots \\
| |
| & & & & \ddots & * \\
| |
| & & & & & * \end{pmatrix}
| |
| </math></center>
| |
|
| |
| Układ z nieosobliwą macierzą trójkątną górną
| |
|
| |
| <center><math>\displaystyle
| |
| U\, x\;=\; c,
| |
| </math></center>
| |
|
| |
| <math>\displaystyle U=(u_{i,j})</math>, <math>\displaystyle c=(c_j)</math>, można rozwiązać stosując algorytm:
| |
|
| |
| {{algorytm|Podstawienie w tył||
| |
| <pre>
| |
|
| |
| <math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
| |
| for (i = N-1; i >= 1; i--)
| |
| <math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
| |
| u_{i,j}x_j^*\right)/u_{i,i}</math>;
| |
| </pre>}}
| |
|
| |
| (Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy
| |
| implikuje, że <math>\displaystyle u_{i,i}\ne 0</math>, <math>\displaystyle \forall i</math>.) Podobnie, układ
| |
| <math>\displaystyle L x= c</math> rozwiązujemy algorytmem:
| |
|
| |
| {{algorytm|Podstawienie w przód||
| |
| <pre>
| |
|
| |
| <math>\displaystyle x_1 = c_1</math>;
| |
| for (i=2; i <= N; i++)
| |
| <math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;
| |
| </pre>}}
| |
|
| |
| Oba algorytmy wymagają rzędu <math>\displaystyle N^2/2</math> mnożeń lub dzieleń i
| |
| <math>\displaystyle N^2/2</math> dodawań lub odejmowań, a więc łącznie <math>\displaystyle O(N^2)</math>
| |
| działań arytmetycznych.
| |
|
| |
| ====Układy z macierzą ortogonalną====
| |
|
| |
| Równie tanio można rozwiązać układ równań
| |
|
| |
| <center><math>\displaystyle
| |
| Q x = b,
| |
| </math></center>
| |
|
| |
| gdy <math>\displaystyle Q</math> jest macierzą ortogonalną, to znaczy <math>\displaystyle Q^TQ = I</math>. Rzeczywiście, z
| |
| ortogonalności mamy natychmiast, że
| |
|
| |
| <center><math>\displaystyle
| |
| x = Q^T b
| |
| </math></center>
| |
|
| |
| i w konsekwencji <math>\displaystyle x</math> można wyznaczyć kosztem takim jak koszt mnożenia macierzy
| |
| przez wektor, czyli <math>\displaystyle O(N^2)</math> operacji.
| |
|
| |
| Podobnie, gdy <math>\displaystyle Q\in C^{N\times N}</math> jest unitarna, to znaczy <math>\displaystyle Q^HQ = I</math>,
| |
| rozwiązaniem układu równań jest
| |
|
| |
| <center><math>\displaystyle
| |
| x = Q^H b.
| |
| </math></center>
| |
|
| |
| ==Metoda eliminacji Gaussa==
| |
|
| |
| [[grafika:Gauss.jpg|thumb|right||Carl Friedrich Gauss<br> [[Biografia Gauss|Zobacz biografię]]]]
| |
|
| |
| W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym
| |
| rozwiązywania układu równań
| |
|
| |
| <center><math>\displaystyle Ax=b</math></center>
| |
|
| |
| okazuje się popularna
| |
| eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm
| |
| ten wyrazimy w języku tzw. <strong>rozkładu LU</strong> macierzy, to znaczy,
| |
| sprowadzającego zadanie do znalezienia
| |
| macierzy trójkątnej dolnej <math>\displaystyle L</math> (z jedynkami na diagonali) oraz trójkątnej górnej
| |
| <math>\displaystyle U</math> takich, że
| |
| <center><math>\displaystyle
| |
| A = LU,
| |
| </math></center>
| |
|
| |
| a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:
| |
|
| |
| {{algorytm|Rozwiązanie układu z wykorzystaniem rozkładu LU||
| |
| <pre>
| |
|
| |
| Znajdź rozkład <math>\displaystyle A=LU</math>;
| |
| Rozwiąż <math>\displaystyle Ly = b</math> przez podstawienie w przód;
| |
| Rozwiąż <math>\displaystyle Ux = y</math> przez podstawienie w tył;
| |
| </pre>}}
| |
|
| |
| Przypuśćmy, że taki rozkład <math>\displaystyle A=LU</math> istnieje. Wówczas, zapisując macierze w postaci
| |
| blokowej, eksponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy,
| |
| mamy
| |
|
| |
| <center><math>\displaystyle
| |
| \begin{pmatrix}
| |
| a_{11} & a_{12}^T\\
| |
| a_{21} & A_{22}
| |
| \end{pmatrix}
| |
| =
| |
| \begin{pmatrix}
| |
| 1 & 0^T\\
| |
| l_{21} & L_{22}
| |
| \end{pmatrix}
| |
| \begin{pmatrix}
| |
| u_{11} & u_{12}^T\\
| |
| 0 & U_{22},
| |
| \end{pmatrix}
| |
| </math></center>
| |
|
| |
| skąd (mnożąc blokowo macierz <math>\displaystyle L</math> przez <math>\displaystyle U</math>) wynika, że
| |
| * <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest
| |
| kopią pierwszego wiersza <math>\displaystyle A</math>,
| |
| * <math>\displaystyle l_{21} = a_{21}/u_{11}</math>, więc pierwsza kolumna <math>\displaystyle L</math> powstaje przez
| |
| podzielenie wszystkich elementów wektora <math>\displaystyle a_{21}</math> przez element na diagonali
| |
| <math>\displaystyle a_{11}</math>,
| |
| * <math>\displaystyle A_{22} - l_{21}u_{12}^T = L_{22}U_{22}</math>, a więc znalezienie podmacierzy
| |
| <math>\displaystyle L_{22}</math> oraz <math>\displaystyle U_{22}</math> sprowadza się do znalezienia rozkładu LU zmodyfikowanego
| |
| bloku <math>\displaystyle A_{22}</math> macierzy <math>\displaystyle A</math>,
| |
| wymiaru <math>\displaystyle (N-1)\times (N-1)</math>.
| |
|
| |
| Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie
| |
| rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
| |
| klasycznych pętli. Jest to ważne w praktyce numerycznej, gdyż rekurencja
| |
| kosztuje: zarówno pamięć, jak i czas.
| |
|
| |
| Ponadto zauważmy, że opisany algorytm możemy wykonać w miejscu, nadpisując
| |
| elementy <math>\displaystyle A</math> elementami macierzy <math>\displaystyle U</math> i <math>\displaystyle L</math> (jedynek z diagonali <math>\displaystyle L</math> nie musimy
| |
| pamiętać, bo wiemy ''a priori'', że tam są).
| |
|
| |
| {{algorytm|Rozkład LU metodą eliminacji Gaussa||
| |
| <pre>
| |
|
| |
| for k=1:N-1
| |
| if <math>\displaystyle a_{kk}</math> == 0
| |
| STOP;
| |
| end
| |
| for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| |
| <math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
| |
| end
| |
| for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
| |
| for i=k+1:N
| |
| <math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
| |
| end
| |
| end
| |
| end
| |
| </pre>}}
| |
|
| |
| Łatwo przekonać się, że <math>\displaystyle k</math>-ty obrót zewnętrznej pętli (tzn. <math>\displaystyle k</math>-ty krok
| |
| algorytmu rozkładu LU) kosztuje rzędu <math>\displaystyle 2(N-k)^2</math> operacji arytmetycznych, skąd łączny koszt tego
| |
| algorytmu rozkładu LU wynosi około <math>\displaystyle \frac{4}{3}N^3</math>.
| |
|
| |
| Jeśli więc wykorzystać rozkład LU do rozwiązywania układu równań <math>\displaystyle Ax=b</math>, to mamy
| |
| następujące zestawienie kosztów:
| |
| * Koszt znalezienia rozkładu <math>\displaystyle A=LU</math>: <math>\displaystyle O(N^3)</math>;
| |
| * Koszt rozwiązania układu <math>\displaystyle Ly=b</math>: <math>\displaystyle O(N^2)</math>;
| |
| * Koszt rozwiązania układu <math>\displaystyle Ux=y</math>: <math>\displaystyle O(N^2)</math>.
| |
|
| |
| Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania
| |
| wynosi już tylko <math>\displaystyle O(N^2)</math>.
| |
|
| |
| {{uwaga|Złożoność obliczeniowa zadania rozwiązania układu równań liniowych||
| |
|
| |
| Z powyższego wynika, że łączny koszt rozwiązania równania liniowego wynosi
| |
| <math>\displaystyle O(N^3)</math>. Można zastanawiać się, jaka jest <strong>najmniejsza możliwa</strong> liczba
| |
| operacji zmiennoprzecinkowych potrzebnych do rozwiązania układu równań
| |
| liniowych.
| |
|
| |
| Można pokazać, że minimalny koszt rozwiązania układu <math>\displaystyle N</math> równań
| |
| liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch
| |
| macierzy <math>\displaystyle N\times N</math>. Tymczasem znany jest całkiem prosty algorytm rekurencyjny,
| |
| wyznaczający iloczyn dwóch macierzy kosztem <math>\displaystyle 4.7\cdot N^{log_27} \approx 4.7
| |
| \cdot N^{2.807}</math> (algorytm Strassena). Bardziej skomplikowany (i praktycznie
| |
| nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt
| |
| <math>\displaystyle O(N^{2.376})</math>. Tak więc równania liniowe daje się (teoretycznie) rozwiązać
| |
| kosztem <math>\displaystyle O(N^{2.376})</math>.
| |
|
| |
| Jednak w praktyce nawet prosty algorytm Strassena zazwyczaj nie jest stosowany.
| |
| Wynika to stąd, że ma trochę gorsze własności numeryczne oraz, co istotniejsze, wymaga dużo
| |
| dodatkowej pamięci na przechowywanie pośrednich wyników.
| |
| }}
| |
|
| |
| ==Wybór elementu głównego==
| |
|
| |
| Opisany powyżej algorytm rozkładu LU niestety czasem może się załamać: gdy
| |
| napotka w czasie działania zerowy element na diagonali zmodyfikowanej
| |
| podmacierzy, np. chociaż macierz
| |
|
| |
| <center><math>\displaystyle
| |
| A = \begin{pmatrix} 0 & 1\\ 1 & 0
| |
| \end{pmatrix}
| |
| </math></center>
| |
|
| |
| jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od
| |
| razu zetknie się z dzieleniem przez <math>\displaystyle a_{11}=0</math>. Ale wystarczy zamienić ze sobą
| |
| kolejnością wiersze macierzy <math>\displaystyle A</math> (to znaczy, w układzie równań, zamienić ze sobą
| |
| miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez
| |
| problemu.
| |
|
| |
| W praktyce obliczeniowej, aby uzyskać algorytm o [[|Dodaj link: możliwie dobrych własnościach
| |
| numerycznych]], wykorzystujemy tzw. strategię <strong>wyboru elementu głównego w
| |
| kolumnie</strong>.
| |
| Polega to na tym, że zanim wykonamy <math>\displaystyle k</math>-ty krok algorytmu rozkładu LU,
| |
| * szukamy w pierwszej kolumnie podmacierzy <math>\displaystyle A(k:N,k:N)</math> szukamy elementu o
| |
| największym module (taki element, na mocy założenia nieosobliwości macierzy,
| |
| jest niezerowy) --- to jest właśnie element główny
| |
| * zamieniamy ze sobą wiersz <math>\displaystyle A(k,1:N)</math> z wierszem, w którym
| |
| znajduje się element główny
| |
| * zapamiętujemy dokonaną permutację, bo potem --- gdy już przyjdzie do
| |
| rozwiązywania układu równań --- będziemy musieli dokonać
| |
| analogicznej permutacji wektora prawej strony
| |
|
| |
| Wynikiem takiego algorytmu jest więc rozkład
| |
|
| |
| <center><math>\displaystyle PA = LU,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle P</math> jest pewną (zerojedynkową) macierzą permutacji (tzn. macierzą
| |
| identyczności z przepermutowanymi wierszami).
| |
|
| |
| Oprócz wyboru elementu głównego w kolumnie, stosuje się czasem inne strategie,
| |
| m.in. wybór w kolumnie oraz tzw. <strong>wybór pełny</strong>, gdy elementu głównego szukamy w
| |
| <strong>całej</strong> podmacierzy <math>\displaystyle A(k:N,k:N)</math>, co znacznie zwiększa liczbę porównań
| |
| niezbędnych do wskazania elementu głównego, ale też trochę poprawia własności
| |
| numeryczne takiego algorytmu.
| |
|
| |
| W praktyce całą informację o wykonanych permutacjach przechowujemy nie w
| |
| zerojedynkowej macierzy j.w. ale w jednym wektorze.
| |
|
| |
| {{algorytm|Rozkład LU z wyborem elementu głównego w kolumnie||
| |
| <pre>
| |
|
| |
| P = 1:N; /* tu zapamiętujemy wykonane permutacje */
| |
| for k=1:N-1
| |
|
| |
| w wektorze A(k:N,k) znajdź element główny <math>\displaystyle a_{pk}</math>;
| |
| zamień ze sobą wiersze A(k,1:N) i A(p,1:N);
| |
| P(k) = p; P(p) = k;
| |
| if <math>\displaystyle a_{kk}</math>
| |
| STOP: macierz osobliwa!
| |
| end
| |
|
| |
| /* kontunuuj tak jak w algorytmie bez wyboru */
| |
|
| |
| for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| |
| <math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
| |
| end
| |
| for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
| |
| for i=k+1:N
| |
| <math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
| |
| end
| |
| end
| |
| end
| |
| </pre>}}
| |
|
| |
| <div class="thumb tright"><div><flash>file=Macierz.swf</flash><div.thumbcaption>Eliminacja Gaussa z wyborem elementu głównego</div></div></div>
| |
|
| |
| Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
| |
| już rozkład LU wykonany z wyborem w kolumnie.
| |
|
| |
| {{algorytm|Rozwiązywanie układu równań metodą eliminacji Gaussa z wyborem elementu głównego w kolumnie||
| |
| <pre>
| |
|
| |
| znajdź rozkład <math>\displaystyle PA = LU</math>;
| |
| rozwiąż względem <math>\displaystyle y</math> układ z macierzą górną trójkątną <math>\displaystyle Uy = Pb</math>;
| |
| rozwiąż względem <math>\displaystyle x</math> układ z macierzą dolną trójkątną <math>\displaystyle Lx = y</math>;
| |
| </pre>}}
| |
|
| |
| Dla niektórych ważnych klas macierzy wiadomo, że rozkład LU jest wykonalny <strong>bez wyznaczania elementu głównego</strong>, co istotnie może zmniejszyć całkowity czas
| |
| działania algorytmu. Jest tak m.in. dla macierzy
| |
| * symetrycznych, dodatnio określonych: <math>\displaystyle A=A^T</math> oraz <math>\displaystyle x^TAx > 0</math>, <math>\displaystyle \forall x
| |
| \neq 0</math>,
| |
| * silnie diagonalnie dominujących: macierz <math>\displaystyle A</math> (lub <math>\displaystyle A^T</math>) spełnia
| |
|
| |
| <center><math>\displaystyle |a_{ii}| > \sum_{j\neq i} |a_{ij}|, \qquad \forall i.
| |
| </math></center>
| |
