Metoda Gaussa-Seidla

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj

Metoda Gaussa-Seidlaiteracyjna metoda numeryczna rozwiązywania układów równań liniowych. Stosowana jest głównie do rozwiązywania ogromnych układów równań postaci \ A \mathbf{x} = \mathbf{b}, w których \ A jest macierzą przekątniowo dominującą. Równania tego typu, obejmujące tysiące a nawet miliony niewiadomych, występują powszechnie w numerycznych metodach rozwiązywania eliptycznych równań różniczkowych cząstkowych, np. równania Laplace'a. Nazwa metody upamiętnia niemieckich matematyków: Carla Friedricha Gaussa i Philippa Ludwiga von Seidla

Obecnie metoda Gaussa-Seidla ma charakter czysto akademicki. Dla małych układów równań dużo szybsze są metody bezpośrednie, np. metoda eliminacji Gaussa, natomiast dla ogromnych układów równań lepszą zbieżność zapewniają metody nadrelaksacyjne oraz wielosiatkowe (ang. multigrid).

Definicja[edytuj | edytuj kod]

Metoda Gaussa-Seidla jest metodą relaksacyjną, w której poszukiwanie rozwiązania rozpoczyna się od dowolnie wybranego rozwiązania próbnego \ x_0, po czym w kolejnych krokach, zwanych iteracjami, za pomocą prostego algorytmu zmienia się kolejno jego składowe, tak by coraz lepiej odpowiadały rzeczywistemu rozwiązaniu. Metoda Gaussa-Seidla bazuje na metodzie Jacobiego, w której krok iteracyjny zmieniono w ten sposób, by każda modyfikacja rozwiązania próbnego korzystała ze wszystkich aktualnie dostępnych przybliżonych składowych rozwiązania. Pozwala to zaoszczędzić połowę pamięci operacyjnej i w większości zastosowań praktycznych zmniejsza ok. dwukrotnie liczbę obliczeń niezbędnych do osiągnięcia zadanej dokładności rozwiązania.

Rozpatrzmy układ \ n równań liniowych z \ n niewiadomymi:

 A x = b.\,

gdzie

  • \ Amacierz układu (\ n\times n),
  • \ b – zadany wektor (\ n składowych),
  • \ x – poszukiwane rozwiązanie (wektor o \ n składowych).

Pojedynczą iterację metody Gaussa-Seidla można zapisać algebraicznie jako


x^{(k+1)}  = \left( {D + L} \right)^{-1} \left( {-U x^{(k)}  + b} \right),

gdzie \ D jest nieosobliwą macierzą diagonalną, a \ L i \ U są odpowiednio macierzą dolnotrójkątną i górnotrójkątną macierzy \ A (tzn. \ L oraz \ U mają zera na głównej przekątnej oraz \ A \equiv D+L+U), natomiast indeks \ k = 0,1,2,\ldots oznacza numer porządkowy iteracji.

Po rozpisaniu na składowe wzór ten przyjmuje postać używaną w implementacjach numerycznych:


x^{(k+1)}_i  = \frac{1}{a_{ii}} \left(b_i - \sum_{j=1}^{i-1}a_{ij}x^{(k+1)}_j-\sum_{j=i+1}^{n}a_{ij}x^{(k)}_j\right),\, i=1,2,\ldots,n.

Uwaga!

  • W powyższych wzorach zakłada się, że w razie potrzeby kolejność równań została zmieniona tak, by dominujące (tj. największe co do modułu w danym równaniu) współczynniki równania znajdowały się na głównej przekątnej macierzy \ A.
  • Jeżeli \ A jest macierzą nieosobliwą, to zawsze można tak przestawić jej wiersze i kolumny, by macierz \ D też była nieosobliwa.
  • Metodę Gaussa-Seidla stosuje się niemal wyłącznie do układów z macierzą przekątniowo dominującą, gdyż w wielu praktycznych zastosowaniach (np. przy rozwiązywaniu eliptycznych równań różniczkowych cząstkowych) jest to łatwy do spełnienia warunek gwarantujący zbieżność metody.
  • Metodę Gaussa-Seidla można stosować także do układów równań liniowych, w których macierz układu nie jest przekątniowo dominująca, ale poza nielicznymi wyjątkami zwykle nie ma gwarancji, że w tym przypadku metoda będzie zbieżna[1].

Warunki zbieżności[edytuj | edytuj kod]

Warunki wystarczające[edytuj | edytuj kod]

1a. Kryterium silnej dominacji w rzędach[edytuj | edytuj kod]

Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy \ A spełniającej warunek ścisłej dominacji przekątniowej w rzędach:

  |a_{ii}| > \sum_{j=1,j\neq i}^{n} |a_{ij}| \qquad \mathrm{dla~wszystkich~}i=1,2,\ldots,n

(źródło: Stoer i Bulirsch 1993).

1b. Kryterium silnej dominacji w kolumnach[edytuj | edytuj kod]

Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy \ A spełniającej warunek ścisłej dominacji przekątniowej w kolumnach:

 |a_{ii}| > \sum_{j=1,j\neq i}^{n} |a_{ji}| \qquad \mathrm{dla~wszystkich~}i=1,2,\ldots,n

(źródło: Stoer i Bulirsch 1993).

2a. Kryterium słabej dominacji w rzędach[edytuj | edytuj kod]

Kolejne kryterium dotyczy nieredukowalnych układów równań liniowych, tj. takich układów, których nie można uporządkować przez permutacje wierszy i kolumn w ten sposób, by niektóre niewiadome można było wyznaczyć poprzez rozwiązanie mniejszej liczby równań niż \ n.

Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej \ A dominują rzędami w sensie słabym

 |a_{ii}| \geqslant \sum_{j=1,j\neq i}^{n} |a_{ij}| \qquad \mathrm{dla~wszystkich~}i=1,2,\ldots,n

oraz jeżeli dla co najmniej jednego wiersza i \in\{1,2,\ldots,n\} zachodzi dominacja silna:

 |a_{ii}| > \sum_{j=1,j\neq i}^{n} |a_{ij}|,

to ciąg iteracji Gaussa-Seidla jest zbieżny (źródło: Stoer i Bulirsch 1993, Tannehill i in., 1997).

Uwagi

  • Macierz \ A jest nieredukowalna, jeżeli poprzez przestawienie wierszy i kolumn nie można jej sprowadzić do postaci blokowej górnej trójkątnej.
  • Nieredukowalność macierzy kwadratowej \ A o rozmiarze \ n\times n można sprawdzić za pomocą grafu skierowanego \ G(A) mającego \ n węzłów \ P_1,P_2,\ldots,P_n, w którym para \ (P_i, P_j) jest połączona (skierowanym) łukiem biegnącym od \ P_i do \ P_j wtedy i tylko wtedy, gdy \ A_{ij} \neq 0. Macierz \ A jest nieredukowalna wtedy i tylko wtedy, gdy między dowolnymi dwoma różnymi węzłami \ P_i, P_j w grafie \ G(A) istnieje połączenie łukami skierowanymi (połączenie to nie musi być bezpośrednie).

2b. Kryterium słabej dominacji w kolumnach[edytuj | edytuj kod]

Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej \ A dominują kolumnami w sensie słabym

 |a_{ii}| \geqslant \sum_{j=1,j\neq i}^{n} |a_{ji}| \qquad \mathrm{dla~wszystkich~}i=1,2,\ldots,n

oraz jeżeli dla co najmniej jednej kolumny i \in\{1,2,\ldots,n\} zachodzi dominacja silna:

 |a_{ii}| > \sum_{j=1,j\neq i}^{n} |a_{ji}|,

to ciąg iteracji Gaussa-Seidla jest zbieżny (źródło: Stoer i Bulirsch 1993).

3. Kryterium dodatniej określoności[edytuj | edytuj kod]

Jeżeli macierz \ A jest dodatnio określona, to metoda Gaussa-Seidla jest zbieżna dla dowolnego wektora początkowego (źródło: Ralston 1983, Stoer i Bulirsch 1993).

Warunek konieczny[edytuj | edytuj kod]

Niech

\ B = -(D+L)^{-1}U

Metoda Gaussa-Seidla jest zbieżna wtedy i tylko wtedy, gdy moduły wszystkich wartości własnych \ B są mniejsze od 1 (źródło: Ralston 1983).

Uwaga: powyższe kryterium jest niepraktyczne i nie jest wykorzystywane w obliczeniach numerycznych.

Warunek zakończenia iteracji[edytuj | edytuj kod]

W praktyce iteracje Gaussa-Seidla kończy się wtedy, gdy dla iteracji o numerze \ k > 1 maksymalna względna zmiana składowej przybliżonego rozwiązania nie przekracza pewnego z góry zadanego małego parametru \varepsilon (np. 10^{-10}):


  \max_i \{|x^{k}_i - x^{k-1}_i|\} < \varepsilon \max_{i}\{|x^{k}_i|\}.

Alternatywny sposób polega na śledzeniu wektora reszt:


  \ r = b-Ax.

Obliczenia przerywa się, gdy \max_i\{|r_i|\} osiągnie wartość mniejszą od pewnego z góry ustalonego małego parametru \ \epsilon.

Uwagi:

  • W metodzie Gaussa-Seidla w każdym kroku modyfikuje się pewną składową rozwiązania (\ x_j), tak by wyzerować odpowiadającą mu składową wektora reszt (\ r_j).
  • Sukcesywne zerowanie jednej lub kilku składowych wektora reszt stanowi istotę wszystkich metod relaksacyjnych.
  • Aktualizacja wektora reszt w kolejnych krokach może być przeprowadzona stosunkowo niewielkim nakładem obliczeń.

Przykłady[edytuj | edytuj kod]

Układ trzech równań liniowych[edytuj | edytuj kod]

Rozważmy następujący układ równań liniowych:

\ x_1 + 6x_2 + x_3 = 9
\ 4x_1 - x_2 + x_3 = 4
\ -x_1 +2x_2 + 5x_3 = 2

W pierwszym i drugim równaniu wyrazy dominujące (\ 6x_2 i \ 4x_1) leżą poza główną przekątną. Po zamianie kolejności tych równań otrzymujemy układ, w którym wartości dominujące leżą na głównej przekątnej:

\ 4x_1 - x_2 + x_3 = 4
\ x_1 + 6x_2 + x_3 = 9
\ -x_1 +2x_2 + 5x_3 = 2

Układ ten spełnia warunek zbieżności metody (macierz układu jest dominująca przekątniowo). Układ zapisujemy w postaci równań na wyrazy dominujące:

\ x_1 = \frac{1}{4}(4 + x_2 - x_3)
\ x_2 = \frac{1}{6}(9 -x_1- x_3)
\ x_3 = \frac{1}{5}(2+x_1 -2x_2)

Dokonujemy wyboru ("zgadujemy") wartości x_2 i x_3, np. x_2 = 0 i x_3 = 0. Następnie podstawiamy te wartości do równania na x_1, uzyskując początkową wartość x_1. Tak uzyskaną wartość podstawiamy do równania na x_2, uzyskując nowe przybliżenie tej niewiadomej. Iteracje kontynuujemy do osiągnięcia określonej dokładności względnej.

Dla x_2 = 0 i x_3 = 0 powyższa procedura daje następujące wyniki (dwie pełne iteracje):

\ x_1 = \frac{1}{4}(4 + 0 - 0) = 1
\ x_2 = \frac{1}{6}(9 - 1 - 0) = 4/3 \approx 1.333
\ x_3 = \frac{1}{5}(2 + 1 - 2\times 4/3) = 1/15 \approx 0.0667
\ x_1 = \frac{1}{4}(4 + 4/3 - 1/15) = 79/60 \approx 1.317
\ x_2 = \frac{1}{6}(9 - 79/60- 1/15)= 457/360 \approx 1.269
\ x_3 = \frac{1}{5}(2 + 79/60 -2\times 457/360) =7/45 \approx 0.1556

Dokładne rozwiązanie: x_1 = 23/18 \approx 1.278, x_2 = 53/42 \approx 1.262, x_3= 19/126 \approx 0.1508.

Jak łatwo sprawdzić, gdyby na początku nie zmieniono kolejności równań, iteracje Gaussa-Seidla byłyby rozbieżne.

Jednowymiarowe równanie Laplace'a[edytuj | edytuj kod]

Jednowymiarowe równanie Laplace'a ma postać \ Ax = b, gdzie \ A jest macierzą trójprzekątniową:


\displaystyle  A =
\begin{pmatrix}
   2  & -1     &    &        &  \\
   -1 & 2      &  -1& \ddots & \\
      & \ddots &  \ddots   &\ddots & -1 \\
      &        &    & -1    & 2
\end{pmatrix}

Macierz \ A , jako pełna macierz trójprzekątniowa, jest nieredukowalna. Wszystkie elementy dominujące znajdują się na głównej przekątnej. Wartość bezwzględna każdego elementu dominującego jest co najmniej równa sumie wartości bezwzględnych pozostałych elementów w danym wierszu. Istnieją dwa elementy dominujące (w pierwszym i ostatnim wierszu, czyli na brzegach układu), których wartość bezwzględna jest większa od sumy wartości bezwzględnych pozostałych elementów wiersza. Dlatego na mocy kryterium dominacji przekątniowej metoda Gaussa-Seidla jest w przypadku tego równania zbieżna.

Ten sam wniosek można wyciągnąć z kryterium dodatniej określoności macierzy \ A, ale wymaga to bardziej zaawansowanych rachunków.

Równanie niezbieżne[edytuj | edytuj kod]

Rozpatrzmy układ równań Ax=0, gdzie


\displaystyle  A =
\begin{pmatrix}
   2  & -1     &  -1   \\
   -1 & 2      &  -1 \\
    -1  &  -1  & 2
\end{pmatrix}

Układ ten ma nieskończenie wiele rozwiązań postaci \ x=c[1,1,1], gdzie \ c jest dowolną liczbą rzeczywistą. Macierz \ A nie spełnia żadnego z opisanych powyżej warunków dostatecznych zbieżności metody Gaussa-Seidla. Mimo tego, jak łatwo sprawdzić, metoda Gaussa-Seidla jest w tym przypadku zbieżna dla każdego wektora początkowego \ x_0; problem w tym, że wartość graniczna zależy od wyboru rozwiązania próbnego \ x_0.

Algorytm[edytuj | edytuj kod]

Wybierz początkowe przybliżenie x^0
for k := 1 step 1 until oczekiwane przybliżenie do
for i := 1 step 1 until n do
 \sigma = 0
for j := 1 step 1 until i-1 do
 \sigma  = \sigma  + a_{ij} x_j^{(k)}
end (j-for)
for j := i+1 step 1 until n do
 \sigma  = \sigma  + a_{ij} x_j^{(k-1)}
end (j-for)
 x_i^{(k)}  = {{\left( {b_i  - \sigma } \right)} \over {a_{ii} }}
end (i-for)
sprawdź czy osiągnięty zostało oczekiwane przybliżenie
end (k-for)
x\approx x^{(k)}

Przypisy

  1. metoda iteracyjna wyrażona równaniem {Q}x^{(k)}  = \left( {Q - A} \right) {x^{(k-1)}  + b} jest zbieżna, gdy ciąg \{ x^{(k)} \} jest zbieżny do x dla dowolnego wektora początkowego x^{(0)}

Bibliografia[edytuj | edytuj kod]

  • David Kincaid, Ward Cheney Analiza Numeryczna ISBN 83-204-3078-X
  • Anthony Ralston, Wstęp do analizy numerycznej, PWN, Warszawa 1983, ISBN 83-01-01626-4.
  • J. Stoer R. Bulirsch, Introduction to Numerical Analysis, Second Edition, Springer-Verlag, New York 1993, ISBN 0-387-97878-X
  • John C. Tannehill, Dale A. Anderson i Richard H. Pletcher, Computational Fluid Mechanics and HeatTransfer (Second Edition), Francis & Taylor, Philadelphia, 1997, ISBN 1-56032-046-X

Linki zewnętrzne[edytuj | edytuj kod]