Konrad Krowicki

Piotr Krupa

L07

Metoda Gaussa-Seidla

Metoda Gaussa-Seidla - iteracyjna metoda numeryczna rozwiązywania układów równań liniowych. Stosowana jest głównie do rozwiązywania ogromnych układów równań postaci
, w których
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:

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:

gdzie

• 
  - macierz układu (
  ,

• 
  - zadany wektor (
  składowych),

• 
  - poszukiwane rozwiązanie (wektor o
  składowych).

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

gdzie
jest nieosobliwą macierzą diagonalną, a
i
są odpowiednio macierzą dolnotrójkątną i górnotrójkątną macierzy
(tzn.
oraz
mają zera na głównej przekątnej oraz
), natomiast indeks
oznacza numer porządkowy iteracji.

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

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
  .

• 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.

Opis programu :

GaussSeidel (n,a,b,mit,eps,x,it,st)

n - liczba równań (równań liczbie niewiadomych)

a - tablica zawierająca wartość elementów macierzy A (element a[i,k] powinien zawierać wartość a_ij, gdzie i,j = 1,2,…,n),

b- tablica zawierająca wartości składowych wektora b(element b[i] powinien zawierać wartość b_i, gdzie i = 1,2,…,n);

mit - maksymalna liczba iteracji w procesie(3),

eps - względna dokładność rozwiązania,

x- tablica zawierająca początkowe przybliżenia wartości x_i (i=1,2,…,n).

Inne parametry

ST- zmienna w procedurze GaussSeidel przypisuje się jedną z następujących wartości :

1, jeżeli n<1,

2, gdy macierz A jest osobliwa,

3, jeśli wymagana dokładność rozwiązania nie jest osiągnięta po mit interajach,

0, w przeciwnym wypadku.

Typy parametrów

Integer : it,mit,n,ST

Extended : eps

matrix: a

vector : b,x

Schemat rozwiązywania:

Na początku macierz współczynników A rozłożymy na sumę trzech macierzy A = L + D + U, gdzie L jest macierzą w której znajdują się elementy których numer wiersza jest większy od numeru kolumny, D to macierz diagonalna z elementami tylko na głównej przekątnej, a U to macierz, w której znajdują się elementy których numery wiersza są mniejsze od numerów kolumny. Można to zapisać następująco:

Następnie obliczymy macierz odwrotną do macierzy D, czyli D^-1. Otrzymamy ją po podniesieniu do potęgi -1 wszystkich wartości na głównej przekątnej macierzy. Po tych operacjach możemy przystąpić już do iteracyjnego obliczania kolejnych przybliżeń rozwiązania według następującego wzoru:

xⁿ⁺¹ = D^-1b - D^-1Lxⁿ⁺¹ - D^-1Uxⁿ
(indeksy n oznaczają tutaj numer iteracji).

Przykład:

Obliczymy następujący układ równań:

4x₁ - x₂ - 0.2x₃ + 2x₄ = 30
-1x₁ + 5x₂ - 2x₄ = 0
0.2x₁ + x₂ + 10x₃ - x₄ = -10
- 2x₂ - x₃ + 4x₄ = 5

Zapiszmy go teraz w postaci Ax = b



Podzielmy teraz macierz A na sumę macierzy L + D + U



Obliczmy teraz macierz D^-1



A teraz kolejno D^-1b, D^-1L, D^-1U



Rozpoczynamy od zerowego przybliżenia czyli x₁⁰ = 0, x₂⁰ = 0, x₃⁰ = 0, x₄⁰ = 0

Obliczmy pierwszą iterację metody, według przytoczonego na początku wzoru:

x₁¹ = 7,5 + 0,25x₂⁰ + 0,05x₃⁰ - 0,5x₄⁰
x₁¹ = 7,5
x₂¹ = 0 + 0,2x₁¹ + 0,4x₄⁰
x₂¹ = 1,5
x₃¹ = -1 - 0,02x₁¹ - 0,1x₂¹ + 0,1x₄⁰
x₃¹ = -1,30
x₄¹ = 1,25 + 0,5x₂¹ + 0,25x₃¹
x₄¹ = 1,675

Kolejna iteracja

x₁² = 7,5 + 0,25x₂¹ + 0,05x₃¹ - 0,5x₄¹
x₁² = 6,9725
x₂² = 0 + 0,2x₁² + 0,4x₄¹
x₂² = 2,0645
x₃² = -1 - 0,02x₁² - 0,1x₂² + 0,1x₄¹
x₃² = -1,1784
x₄² = 1,25 + 0,5x₂² + 0,25x₃²
x₄² = 1,98765

Można teraz obliczyć kolejną iterację. Każda z nich przybliża nas do dokładnego wyniku.

Wyniki otrzymane przy użyciu algorytmu:

a) 1 iteracja:

b) 2 iteracje:

c) 100 iteracji:

Przykład 2 (układ równań z 4.3.2:

a)

Dla danych:

Wyniki:

x[1] = -3.87215061601966E-0021

x[2] = 2.00000000000000E-0001

x[3] = 2.00000000000000E-0001

x[4] = 4.00000000000000E-0001

it = 46, st = 0

Komentarz: W tym układzie równań do rozwiązania go, wystarczyło tylko 46 iteracji.

b)

Dla danych:

Wyniki:

x[1] = -8.53655928498797E-0001

x[2] = -7.75175766736041E+0000

x[3] = 6.86615393897085E-0002

x[4] = 1.07951328548122E+0000

it = 10, st = 3

Komentarz: Przy rozwiązywaniu tego układu równań, okazało się, że 10 iteracji to za mało i dokładność rozwiązania nie została osiągnięta. Powyższe wyniki to ostatnie obliczone przybliżenie rozwiązania.

c)

Dla danych:

Wyniki:

x[1] = -8.53655931618867E-0001

x[2] = -7.75175767876022E+0000

x[3] = 6.86615398239442E-0002

x[4] = 1.07951328720871E+0000

it = 18, st = 0

Komentarz: W tym podpunkcie zadania rozwiązywaliśmy ponownie przykład z podpunktu b, tym razem zwiększyliśmy ilość iteracji z 10 na 50, program rozwiązał układ równań z wymaganą dokładnością po 18 iteracjach.

Wnioski:

Metoda Gaussa-Seidla pozwala uzyskać tym dokładniejsze wyniki im więcej wykonamy iteracji. Jest to o tyle dobre, że dzięki temu sami możemy ustalić z jaką dokładnością otrzymamy końcowy wynik.

Wynik

x[1] = 7,5

x[2] = 1,5

x[3] = -1,3

x[4] = 1,675

Wynik

x[1] = 6,9725

x[2] = 2,0645

x[3] = -1,1784

x[4] = 1,98765

Wynik

x[1] = 6,96156056445949

x[2] = 2,2211225336824

x[3] = -1,1541408594598

x[4] = 2,07202605197625

---