Tekst źródłowy strony Metoda Gaussa-Seidla

'''Metoda Gaussa-Seidla''' – [[algorytm iteracyjny|iteracyjna metoda]] [[Metoda numeryczna|numerycznego]] rozwiązywania [[Układ równań liniowych|układów równań liniowych]]. Nazwa upamiętnia [[Niemcy|niemieckich]] [[matematyk]]ów: [[Carl Friedrich Gauss|Carla Friedricha Gaussa]] i [[Philipp Ludwig von Seidel|Philippa Ludwiga von Seidla]].

Metoda stosowana jest głównie do rozwiązywania układów o dużej liczbie równań i niewiadomych (nawet rzędu milionów), których macierz główna jest [[macierz przekątniowo dominująca|macierzą przekątniowo dominującą]]. Równania tego typu występują powszechnie podczas rozwiązywania [[równanie różniczkowe cząstkowe|równań różniczkowych cząstkowych]], np. [[równanie różniczkowe Laplace’a|równania Laplace’a]]. 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 ==
{{zobacz też|układ równań liniowych}}
Metoda Gaussa-Seidla jest metodą relaksacyjną, w której poszukiwanie rozwiązania rozpoczyna się od dowolnie wybranego rozwiązania próbnego <math>\mathbf x^{(0)},</math> 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 [[metoda Jacobiego|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 <math>n</math> równań liniowych z <math>n</math> niewiadomymi:
: <math>\mathbf{Ax} = \mathbf b.</math>

Pojedynczą iterację metody Gaussa-Seidla można zapisać algebraicznie jako
: <math>\mathbf x^{(k+1)} = (\mathbf D + \mathbf L)^{-1} \left(-\mathbf{Ux}^{(k)} + \mathbf b\right),</math>

gdzie:
: <math>\mathbf D</math> – [[Macierz odwrotna#Odwracalność a nieosobliwość|nieosobliwa]] [[macierz diagonalna]],
: <math>\mathbf L</math> i <math>\mathbf U</math> – odpowiednio macierz [[macierz trójkątna|dolnotrójkątna]] i [[macierz trójkątna|górnotrójkątna]] macierzy <math>\mathbf A</math> (tzn. <math>\mathbf L</math> oraz <math>\mathbf U</math> mają zera na głównej przekątnej oraz <math>\mathbf A \equiv \mathbf D + \mathbf L + \mathbf U</math>),
: indeks <math>k = 0, 1, 2, \dots</math> – numer porządkowy iteracji.

Po rozpisaniu na składowe wzór ten przyjmuje postać używaną w implementacjach numerycznych:
: <math>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), \quad (i = 1, 2, \dots, n).</math>

; 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 <math>\mathbf A.</math>
* Jeżeli <math>\mathbf A</math> jest [[Macierz odwrotna#Odwracalność a nieosobliwość|macierzą nieosobliwą]], to zawsze można tak przestawić jej wiersze i kolumny, by macierz <math>\mathbf D</math> też była nieosobliwa.
* Metodę Gaussa-Seidla stosuje się niemal wyłącznie do układów z macierzą [[macierz przekątniowo dominująca|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<ref group="uwaga">Metoda iteracyjna wyrażona równaniem <math>\mathbf{Qx}^{(k)} = (\mathbf Q - \mathbf A) {\mathbf x^{(k-1)} + \mathbf b}</math> jest zbieżna, gdy ciąg <math>(\mathbf x^{(k)})</math> jest zbieżny do <math>\mathbf x</math> dla dowolnego wektora początkowego <math>\mathbf x^{(0)}.</math></ref>.

== Warunki zbieżności ==
=== Warunki wystarczające ===
==== Kryterium silnej dominacji w rzędach ====
Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy <math>\mathbf A</math> spełniającej warunek ścisłej dominacji przekątniowej w rzędach<ref name="St93">Stoer i Bulirsch, 1993.</ref>:
: <math>|a_{ii}| > \sum_{j = 1, j \neq i}^n |a_{ij}| \qquad \mathrm{dla~wszystkich~}i = 1, 2, \dots, n.</math>

==== Kryterium silnej dominacji w kolumnach ====
Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy <math>\mathbf A</math> spełniającej warunek ścisłej dominacji przekątniowej w kolumnach<ref name="St93" />:
: <math>|a_{jj}| > \sum_{i = 1, i \neq j}^n |a_{ij}| \qquad \mathrm{dla~wszystkich~}j = 1, 2, \dots, n.</math>

==== Kryterium słabej dominacji w rzędach ====
Kolejne kryterium dotyczy nieredukowalnych układów równań liniowych<ref group="uwaga">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ż <math>n.</math></ref>. Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej<ref group="uwaga" name="red">Macierz <math>\mathbf A</math> jest nieredukowalna, jeżeli poprzez przestawienie wierszy i kolumn nie można jej sprowadzić do postaci blokowej górnej trójkątnej.</ref><ref group="uwaga">Nieredukowalność macierzy kwadratowej <math>\mathbf A</math> stopnia <math>n</math> można sprawdzić za pomocą [[Graf (matematyka)|grafu skierowanego]] <math>G(\mathbf A)</math> mającego <math>n</math> węzłów <math>P_1, P_2, \dots, P_n,</math> w którym para <math>(P_i, P_j)</math> jest połączona ścieżką od <math>P_i</math> do <math>P_j</math> wtedy i tylko wtedy, gdy <math>a_{ij} \neq 0.</math> Macierz <math>\mathbf A</math> jest nieredukowalna wtedy i tylko wtedy, gdy między dowolnymi dwoma różnymi węzłami <math>P_i, P_j</math> w grafie <math>G(\mathbf A)</math> istnieje połączenie ścieżkami.</ref> <math>\mathbf A</math> dominują rzędami w sensie słabym
: <math>|a_{ii}| \geqslant \sum_{j = 1, j \neq i}^n |a_{ij}| \qquad \mathrm{dla~wszystkich~}i = 1, 2, \dots, n</math>

oraz jeżeli dla co najmniej jednego wiersza <math>i \in \{1, 2, \dots, n\}</math> zachodzi dominacja silna:
: <math>|a_{ii}| > \sum_{j = 1, j \neq i}^n |a_{ij}|,</math>

to ciąg iteracji Gaussa-Seidla jest zbieżny<ref name="St93" /><ref name="Ta97">Tannehill i in., 1997.</ref>.

==== Kryterium słabej dominacji w kolumnach ====
Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej<ref group="uwaga" name="red" /> <math>\mathbf A</math> dominują kolumnami w sensie słabym
: <math>|a_{jj}| \geqslant \sum_{i = 1, i \neq j}^n |a_{ij}| \qquad \mathrm{dla~wszystkich~}i = 1, 2, \dots, n</math>

oraz jeżeli dla co najmniej jednej kolumny <math>j \in \{1, 2, \dots, n\}</math> zachodzi dominacja silna:
: <math>|a_{jj}| > \sum_{i = 1, i \neq j}^n |a_{ij}|,</math>

to ciąg iteracji Gaussa-Seidla jest zbieżny<ref name="St93" />.

==== Kryterium dodatniej określoności ====
Jeżeli macierz <math>\mathbf A</math> jest [[Określoność formy|dodatnio określona]], to metoda Gaussa-Seidla jest zbieżna dla dowolnego wektora początkowego<ref name="St93" /><ref name="Ra83">Ralston, 1983.</ref>.

=== Warunek konieczny i wystarczający ===
Niech
: <math>\mathbf B = -(\mathbf D + \mathbf L)^{-1} \mathbf U</math>

Metoda Gaussa-Seidla jest zbieżna wtedy i tylko wtedy, gdy moduły wszystkich wartości własnych <math>\mathbf B</math> są mniejsze od 1<ref name="Ra83" />.

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

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

Alternatywny sposób polega na śledzeniu wektora reszt:
: <math>\mathbf r = \mathbf b - \mathbf{Ax}.</math>

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

Uwagi:
* W metodzie Gaussa-Seidla w każdym kroku modyfikuje się pewną składową rozwiązania (<math>x_j</math>), tak by wyzerować odpowiadającą mu składową wektora reszt (<math>r_j</math>).
* 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 ==
=== Układ trzech równań liniowych ===
Rozważmy następujący [[układ równań]] liniowych:
: <math>\begin{cases} x_1 + 6x_2 + x_3 = 9 \\ 4x_1 - x_2 + x_3 = 4 \\ -x_1 +2x_2 + 5x_3 = 2 \end{cases}</math>

W pierwszym i drugim równaniu wyrazy dominujące (<math>6x_2</math> i <math>4x_1</math>) 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:
: <math>4x_1 - x_2 + x_3 = 4</math>
: <math>x_1 + 6x_2 + x_3 = 9</math>
: <math>-x_1 +2x_2 + 5x_3 = 2</math>

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:
: <math>x_1 = \tfrac{1}{4}(4 + x_2 - x_3)</math>
: <math>x_2 = \tfrac{1}{6}(9 -x_1- x_3)</math>
: <math>x_3 = \tfrac{1}{5}(2+x_1 -2x_2)</math>

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

Dla <math>x_2 = 0</math> i <math>x_3 = 0</math> powyższa procedura daje następujące wyniki (dwie pełne iteracje):
: <math>x_1 = \tfrac{1}{4}(4 + 0 - 0) = 1</math>
: <math>x_2 = \tfrac{1}{6}(9 - 1 - 0) = \tfrac{4}{3} \approx 1{,}333</math>
: <math>x_3 = \tfrac{1}{5}(2 + 1 - 2 \cdot \tfrac{4}{3}) = \tfrac{1}{15} \approx 0{,}0667</math>
: <math>x_1 = \tfrac{1}{4}(4 + \tfrac{4}{3} - \tfrac{1}{15}) = \tfrac{79}{60} \approx 1{,}317</math>
: <math>x_2 = \tfrac{1}{6}(9 - \tfrac{79}{60}- \tfrac{1}{15})= \tfrac{457}{360} \approx 1{,}269</math>
: <math>x_3 = \tfrac{1}{5}(2 + \tfrac{79}{60} -2 \cdot \tfrac{457}{360}) = \tfrac{7}{45} \approx 0{,}1556</math>

Dokładne rozwiązanie: <math>x_1 = \tfrac{23}{18} \approx 1{,}278,</math> <math>x_2 = \tfrac{53}{42} \approx 1{,}262,</math> <math>x_3= \tfrac{19}{126} \approx 0{,}1508.</math>

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 ===
Jednowymiarowe równanie Laplace’a ma postać <math>\mathbf{Ax} = \mathbf b,</math> gdzie <math>\mathbf A</math> jest macierzą trójprzekątniową:
: <math>\mathbf A = \begin{bmatrix}
   2 & -1 & & & \\
  -1 & 2 & -1 & \ddots & \\
     & \ddots & \ddots & \ddots & -1 \\
     & & & -1 & 2
\end{bmatrix}</math>

Macierz <math>\mathbf A,</math> jako pełna macierz trójprzekątniowa, jest nieredukowalna<ref group="uwaga" name="red" />. 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 <math>\mathbf A,</math> ale wymaga to bardziej zaawansowanych rachunków.

=== Równanie niezbieżne ===
Rozpatrzmy układ równań <math>\mathbf{Ax} = \boldsymbol\Theta,</math> gdzie
: <math>\mathbf A =
\begin{bmatrix}
   2 & -1 & -1 \\
  -1 & 2 & -1 \\
  -1 & -1 & 2
\end{bmatrix}</math>

Układ ten ma nieskończenie wiele rozwiązań postaci <math>\mathbf x = c[1, 1, 1],</math> gdzie <math>c</math> jest dowolną [[Liczby rzeczywiste|liczbą rzeczywistą]]. Macierz <math>\mathbf A</math> 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 <math>x_0;</math> problem w tym, że wartość graniczna zależy od wyboru rozwiązania próbnego <math>\mathbf x^{(0)}.</math>

== Algorytm ==
: ''Wybierz początkowe przybliżenie'' <math>\mathbf x^{(0)}</math>
: '''for''' k := 1 '''step''' 1 '''until''' ''oczekiwane przybliżenie'' '''do'''
:: '''for''' i := 1 '''step''' 1 '''until''' n '''do'''
::: <math>\sigma = 0</math>
::: '''for''' j := 1 '''step''' 1 '''until''' i-1 '''do'''
:::: <math>\sigma = \sigma + a_{ij} x_j^{(k)}</math>
::: '''end''' (j-for)
::: '''for''' j := i+1 '''step''' 1 '''until''' n '''do'''
:::: <math>\sigma = \sigma + a_{ij} x_j^{(k-1)}</math>
::: '''end''' (j-for)
::: <math>x_i^{(k)} = \frac{\left( {b_i - \sigma } \right)}{a_{ii}}</math>
:: '''end''' (i-for)
:: ''sprawdź, czy osiągnięto oczekiwane przybliżenie''
: '''end''' (k-for)
: <math>\mathbf x \approx \mathbf x^{(k)}</math>

=== Przykład w Python 3 i pakiecie NumPy ===
<syntaxhighlight lang="numpy">
import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
       [-1., 11., -1., 3.],
       [2., -1., 10., -1.],
       [0.0, 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

# prints the system
print("System:")
for i in range(A.shape[0]):
  row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
  print(" + ".join(row), "=", b[i])
print()

x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
  print("Current solution:", x)
  x_new = np.zeros_like(x)

  for i in range(A.shape[0]):
    s1 = np.dot(A[i, :i], x_new[:i])
    s2 = np.dot(A[i, i + 1:], x[i + 1:])
    x_new[i] = (b[i] - s1 - s2) / A[i, i]

  if np.allclose(x, x_new, rtol=1e-8):
    break

  x = x_new

print("Solution:")
print(x)
error = np.dot(A, x) - b
print("Error:")
print(error)
</syntaxhighlight>

Powyższy przykład wyświetla wynik:

<syntaxhighlight lang="python">
System:
10.0*x1 + -1.0*x2 + 2.0*x3 + 0.0*x4 = 6.0
-1.0*x1 + 11.0*x2 + -1.0*x3 + 3.0*x4 = 25.0
2.0*x1 + -1.0*x2 + 10.0*x3 + -1.0*x4 = -11.0
0.0*x1 + 3.0*x2 + -1.0*x3 + 8.0*x4 = 15.0

Current solution: [ 0. 0. 0. 0.]
Current solution: [ 0.6     2.32727273 -0.98727273 0.87886364]
Current solution: [ 1.03018182 2.03693802 -1.0144562  0.98434122]
Current solution: [ 1.00658504 2.00355502 -1.00252738 0.99835095]
Current solution: [ 1.00086098 2.00029825 -1.00030728 0.99984975]
Current solution: [ 1.00009128 2.00002134 -1.00003115 0.9999881 ]
Current solution: [ 1.00000836 2.00000117 -1.00000275 0.99999922]
Current solution: [ 1.00000067 2.00000002 -1.00000021 0.99999996]
Current solution: [ 1.00000004 1.99999999 -1.00000001 1.    ]
Current solution: [ 1. 2. -1. 1.]
Solution:
[ 1. 2. -1. 1.]
Error:
[ 2.06480930e-08 -1.25551054e-08  3.61417563e-11  0.00000000e+00]
</syntaxhighlight>

== Uwagi ==
{{Uwagi}}

== Przypisy ==
{{Przypisy}}

== Bibliografia ==
* {{Cytuj|autor = Anthony Ralston |tytuł = Wstęp do analizy numerycznej |data = 1983 |isbn = 83-01-01626-4 |inni = Roman Zuber (tłum.); Stefan Paszkowski (oprac.) |miejsce = Warszawa |wydawca = PWN |oclc = 749823556 }}
* {{Cytuj|autor = Josef Stoer |tytuł = Introduction to Numerical Analysis, Second Edition |data = 1993 |isbn = 0-387-97878-X |wydanie = 2nd ed |inni = Roland Bulirsch |miejsce = New York |wydawca = Springer-Verlag |oclc = 26097604 }}
* John C. Tannehill, Dale A. Anderson i Richard H. Pletcher, ''Computational Fluid Mechanics and Heat Transfer (Second Edition)'', Francis & Taylor, Philadelphia, 1997, {{ISBN|1-56032-046-X}}.

== Literatura dodatkowa ==
* David Kincaid, Ward Cheney, ''Analiza Numeryczna'' {{ISBN|83-204-3078-X}}.

== Linki zewnętrzne ==
* [http://wazniak.mimuw.edu.pl/index.php?title=MN08#Metoda_Gaussa-Seidla Metoda Gaussa-Seidla]

[[Kategoria:Metody iteracyjne]]
[[Kategoria:Układy równań liniowych]]