Metoda Gaussa-Seidla – iteracyjna metoda numerycznego rozwiązywania układów równań liniowych. Nazwa upamiętnia niemieckich matematyków: Carla Friedricha Gaussa i 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ącą. Równania tego typu występują powszechnie podczas rozwiązywania równań różniczkowych cząstkowych, np. 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).
Metoda Gaussa-Seidla jest metodą relaksacyjną, w której poszukiwanie rozwiązania rozpoczyna się od dowolnie wybranego rozwiązania próbnego
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
równań liniowych z
niewiadomymi:
![{\displaystyle \mathbf {Ax} =\mathbf {b} .}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7a38890b1321031020678764f373239161242ad2)
Pojedynczą iterację metody Gaussa-Seidla można zapisać algebraicznie jako
![{\displaystyle \mathbf {x} ^{(k+1)}=(\mathbf {D} +\mathbf {L} )^{-1}\left(-\mathbf {Ux} ^{(k)}+\mathbf {b} \right),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a3fe0c8f3139460490c286587f13595b9de6115f)
gdzie:
– nieosobliwa macierz diagonalna,
i
– odpowiednio macierz dolnotrójkątna i górnotrójkątna macierzy
(tzn.
oraz
mają zera na głównej przekątnej oraz
),- indeks
– numer porządkowy iteracji.
Po rozpisaniu na składowe wzór ten przyjmuje postać używaną w implementacjach numerycznych:
![{\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad (i=1,2,\dots ,n).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/901e77a32f2aa62284c80841c43f8e41312f154d)
- 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
![{\displaystyle \mathbf {A} .}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ae3d0814d4c81a814b240ffd86eff80b36945131)
- Jeżeli
jest macierzą nieosobliwą, to zawsze można tak przestawić jej wiersze i kolumny, by macierz
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[a].
Warunki wystarczająceedytuj kod
Kryterium silnej dominacji w rzędachedytuj kod
Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy
spełniającej warunek ścisłej dominacji przekątniowej w rzędach[1]:
![{\displaystyle |a_{ii}|>\sum _{j=1,j\neq i}^{n}|a_{ij}|\qquad \mathrm {dla~wszystkich~} i=1,2,\dots ,n.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f0f3cc997d95dda6b4ce74a20f9a58ab0ded6d53)
Kryterium silnej dominacji w kolumnachedytuj kod
Metoda Gaussa-Seidla jest zbieżna dla każdej macierzy
spełniającej warunek ścisłej dominacji przekątniowej w kolumnach[1]:
![{\displaystyle |a_{jj}|>\sum _{i=1,i\neq j}^{n}|a_{ij}|\qquad \mathrm {dla~wszystkich~} j=1,2,\dots ,n.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/007d057ff6fb9cc9fe15c4714a347909bb0a63b0)
Kryterium słabej dominacji w rzędachedytuj kod
Kolejne kryterium dotyczy nieredukowalnych układów równań liniowych[b]. Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej[c][d]
dominują rzędami w sensie słabym
![{\displaystyle |a_{ii}|\geqslant \sum _{j=1,j\neq i}^{n}|a_{ij}|\qquad \mathrm {dla~wszystkich~} i=1,2,\dots ,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/674b3b1c63f23b7b4e221b824dea15b550135496)
oraz jeżeli dla co najmniej jednego wiersza
zachodzi dominacja silna:
![{\displaystyle |a_{ii}|>\sum _{j=1,j\neq i}^{n}|a_{ij}|,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/204cf2e0d3ac551014b79e3a3e16e3cf8e50f734)
to ciąg iteracji Gaussa-Seidla jest zbieżny[1][2].
Kryterium słabej dominacji w kolumnachedytuj kod
Jeżeli wszystkie wyrazy diagonalne macierzy nieredukowalnej[c]
dominują kolumnami w sensie słabym
![{\displaystyle |a_{jj}|\geqslant \sum _{i=1,i\neq j}^{n}|a_{ij}|\qquad \mathrm {dla~wszystkich~} i=1,2,\dots ,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8d02c6f6f14cce166b9d21395e2e1b001912dd4a)
oraz jeżeli dla co najmniej jednej kolumny
zachodzi dominacja silna:
![{\displaystyle |a_{jj}|>\sum _{i=1,i\neq j}^{n}|a_{ij}|,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d23960c4edae5b5c8c8cb402051b72ca99caf92f)
to ciąg iteracji Gaussa-Seidla jest zbieżny[1].
Kryterium dodatniej określonościedytuj kod
Jeżeli macierz
jest dodatnio określona, to metoda Gaussa-Seidla jest zbieżna dla dowolnego wektora początkowego[1][3].
Warunek konieczny i wystarczającyedytuj kod
Niech
![{\displaystyle \mathbf {B} =-(\mathbf {D} +\mathbf {L} )^{-1}\mathbf {U} }](https://wikimedia.org/api/rest_v1/media/math/render/svg/b54c5c0d9253cb710f1ee59ebe478acb65d52fcc)
Metoda Gaussa-Seidla jest zbieżna wtedy i tylko wtedy, gdy moduły wszystkich wartości własnych
są mniejsze od 1[3].
Uwaga: powyższe kryterium jest niepraktyczne i nie jest wykorzystywane w obliczeniach numerycznych.
Warunek zakończenia iteracjiedytuj kod
W praktyce iteracje Gaussa-Seidla kończy się wtedy, gdy dla iteracji o numerze
maksymalna względna zmiana składowej przybliżonego rozwiązania nie przekracza pewnego z góry zadanego małego parametru
(np.
):
![{\displaystyle \max _{i}|x_{i}^{k}-x_{i}^{k-1}|<\varepsilon \max _{i}|x_{i}^{k}|.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/88f4b60be8f8a149014a964f76b93ce9153c3d7f)
Alternatywny sposób polega na śledzeniu wektora reszt:
![{\displaystyle \mathbf {r} =\mathbf {b} -\mathbf {Ax} .}](https://wikimedia.org/api/rest_v1/media/math/render/svg/eda4d1733c4f7bf88e9dac73f0ebb0230815159c)
Obliczenia przerywa się, gdy
osiągnie wartość mniejszą od pewnego z góry ustalonego małego parametru ![{\displaystyle \epsilon .}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fd314041977406051e589626aa6b0a1d74301fa3)
Uwagi:
- W metodzie Gaussa-Seidla w każdym kroku modyfikuje się pewną składową rozwiązania (
), tak by wyzerować odpowiadającą mu składową wektora reszt (
). - 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ń.
Układ trzech równań liniowychedytuj kod
Rozważmy następujący układ równań liniowych:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a6918c810a2d8a6504b28807528ac31c277996e9)
W pierwszym i drugim równaniu wyrazy dominujące (
i
) 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:
![{\displaystyle 4x_{1}-x_{2}+x_{3}=4}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a6d218d635257235fc865b550b1b6136c5a3f98e)
![{\displaystyle x_{1}+6x_{2}+x_{3}=9}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aa173e7df4bf7590e6be1c2f6c39dd08286b9052)
![{\displaystyle -x_{1}+2x_{2}+5x_{3}=2}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f694620ffbd5566cea7c871d509482fa11b9545e)
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:
![{\displaystyle x_{1}={\tfrac {1}{4}}(4+x_{2}-x_{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9e37f34b3ec797f30eba762445d854b52a59022b)
![{\displaystyle x_{2}={\tfrac {1}{6}}(9-x_{1}-x_{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/97c2d4e4e2e1327b8b679651bedf2b4ceebc2129)
![{\displaystyle x_{3}={\tfrac {1}{5}}(2+x_{1}-2x_{2})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8bd873431cf3c2c7367400faad48d773de26c8d5)
Dokonujemy wyboru („zgadujemy”) wartości
i
np.
i
Następnie podstawiamy te wartości do równania na
uzyskując początkową wartość
Tak uzyskaną wartość podstawiamy do równania na
uzyskując nowe przybliżenie tej niewiadomej. Iteracje kontynuujemy do osiągnięcia określonej dokładności względnej.
Dla
i
powyższa procedura daje następujące wyniki (dwie pełne iteracje):
![{\displaystyle x_{1}={\tfrac {1}{4}}(4+0-0)=1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/56a2a025b6561a481faa2e693c88ec6449796f53)
![{\displaystyle x_{2}={\tfrac {1}{6}}(9-1-0)={\tfrac {4}{3}}\approx 1{,}333}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fee367b617e93dcb0c82557962bf0addfd950a37)
![{\displaystyle x_{3}={\tfrac {1}{5}}(2+1-2\cdot {\tfrac {4}{3}})={\tfrac {1}{15}}\approx 0{,}0667}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ba02687a2020aaf5624031fe9312dba9e228f013)
![{\displaystyle x_{1}={\tfrac {1}{4}}(4+{\tfrac {4}{3}}-{\tfrac {1}{15}})={\tfrac {79}{60}}\approx 1{,}317}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fe1cd5e6b8a24fbd82c1f160250a6d1eb52ac6c8)
![{\displaystyle x_{2}={\tfrac {1}{6}}(9-{\tfrac {79}{60}}-{\tfrac {1}{15}})={\tfrac {457}{360}}\approx 1{,}269}](https://wikimedia.org/api/rest_v1/media/math/render/svg/444f7e3cc110d967950ffe320c363f62c13aef3b)
![{\displaystyle x_{3}={\tfrac {1}{5}}(2+{\tfrac {79}{60}}-2\cdot {\tfrac {457}{360}})={\tfrac {7}{45}}\approx 0{,}1556}](https://wikimedia.org/api/rest_v1/media/math/render/svg/35281473536f920b55ab91a8038344016bd8909d)
Dokładne rozwiązanie:
![{\displaystyle x_{3}={\tfrac {19}{126}}\approx 0{,}1508.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c25dae523760d39660d4650fecea43246b4bac22)
Jak łatwo sprawdzić, gdyby na początku nie zmieniono kolejności równań, iteracje Gaussa-Seidla byłyby rozbieżne.
Jednowymiarowe równanie Laplace’aedytuj kod
Jednowymiarowe równanie Laplace’a ma postać
gdzie
jest macierzą trójprzekątniową:
![{\displaystyle \mathbf {A} ={\begin{bmatrix}2&-1&&&\\-1&2&-1&\ddots &\\&\ddots &\ddots &\ddots &-1\\&&&-1&2\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6180719032cc9a2500057eacc35680228a329211)
Macierz
jako pełna macierz trójprzekątniowa, jest nieredukowalna[c]. 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
ale wymaga to bardziej zaawansowanych rachunków.
Rozpatrzmy układ równań
gdzie
![{\displaystyle \mathbf {A} ={\begin{bmatrix}2&-1&-1\\-1&2&-1\\-1&-1&2\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/81142455fa70322c4b585cbb77edaa63ece3d27d)
Układ ten ma nieskończenie wiele rozwiązań postaci
gdzie
jest dowolną liczbą rzeczywistą. Macierz
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
problem w tym, że wartość graniczna zależy od wyboru rozwiązania próbnego ![{\displaystyle \mathbf {x} ^{(0)}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dd2c065bf74a50c04ecef017b25d7bfa8834df98)
- Wybierz początkowe przybliżenie
![{\displaystyle \mathbf {x} ^{(0)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/420304397f5d97b8e25a6954a2016c515fb4222b)
- for k := 1 step 1 until oczekiwane przybliżenie do
- for i := 1 step 1 until n do
![{\displaystyle \sigma =0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1eb4831f1e0ca1ba7d007dc6b973e54787e1a4b4)
- for j := 1 step 1 until i-1 do
![{\displaystyle \sigma =\sigma +a_{ij}x_{j}^{(k)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ee0c47c2b53e85c5de51c3c7040fb0fa970bc6e7)
- end (j-for)
- for j := i+1 step 1 until n do
![{\displaystyle \sigma =\sigma +a_{ij}x_{j}^{(k-1)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7de122eccb4f79177aff10e0a9ee8f3f6d0477c8)
- end (j-for)
![{\displaystyle x_{i}^{(k)}={\frac {\left({b_{i}-\sigma }\right)}{a_{ii}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/af9877f3923a77eca13544c2b1c45d2fcda9bd7b)
- end (i-for)
- sprawdź, czy osiągnięto oczekiwane przybliżenie
- end (k-for)
![{\displaystyle \mathbf {x} \approx \mathbf {x} ^{(k)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/74533ba5ef8cd919e6b65f45f25d0144ccb26133)
Przykład w Python 3 i pakiecie NumPyedytuj kod
import numpy as npITERATION_LIMIT = 1000# initialize the matrixA = np.array([[10., -1., 2., 0.], [-1., 11., -1., 3.], [2., -1., 10., -1.], [0.0, 3., -1., 8.]])# initialize the RHS vectorb = np.array([6., 25., -11., 15.])# prints the systemprint("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_newprint("Solution:")print(x)error = np.dot(A, x) - bprint("Error:")print(error)
Powyższy przykład wyświetla wynik:
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.02.0*x1 + -1.0*x2 + 10.0*x3 + -1.0*x4 = -11.00.0*x1 + 3.0*x2 + -1.0*x3 + 8.0*x4 = 15.0Current 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]
- AnthonyA. Ralston AnthonyA., Wstęp do analizy numerycznej, RomanR. Zuber (tłum.), StefanS. Paszkowski (oprac.), Warszawa: PWN, 1983, ISBN 83-01-01626-4, OCLC 749823556 . Brak numerów stron w książce
- JosefJ. Stoer JosefJ., Introduction to Numerical Analysis, Second Edition, RolandR. Bulirsch, wyd. 2nd ed, New York: Springer-Verlag, 1993, ISBN 0-387-97878-X, OCLC 26097604 . Brak numerów stron w książce
- 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.