Metoda Gauss-Seidel

Metoda Gauss-Seidel

Breviar teoretic

Metoda Gauss-Seidel este o metoda de rezolvare numerica a sistemelor de tip Cramer, prin aproximatii succesive. Matricea A a sistemului se descompune in suma L + D + U , unde L este matrice triunghiulara subdiagonala, D matrice diagonala si U matrice triunghiulara supradiagonala.

Pentru un vector x⁽⁰⁾ Rⁿ, sirul de vectori x^(k) definit prin:

x^(k+1) = (L + D)⁻¹(b − Ux^(k))

se numeste traiectoria Gauss-Seidel a vectorului x⁽⁰⁾.

Traiectoria Gauss-Seidel a vectorului x⁽⁰⁾ converge daca si numai daca raza spectrala ρ a matricei

−(L + D)⁻¹ U (2.35)

este strict subunitara.

In caz de    convergenta, componentele ale vectorului, situat pe traiectoria Gauss-Seidel a vectorului x⁽⁰⁾, sunt date de relatiile:

(2.36)

(2.37)

Problema rezolvata

Sa se determine primele 3 puncte de pe traiectoria Gauss-Seidel a vectorului (0, 0)^T pentru sistemul urmator:

Rezolvare

Sistemul se poate scrie sub forma Ax = b, unde

iar matricea A se descompune in suma L + D + U, dupa cum urmeaza:

.

Algoritmul converge daca raza spectrala a matricei

M = −(L + D)^-1U =

este strict subunitara. Efectuand calculele, obtinem

(M) =

si deci algoritmul este convergent.

In continuare, aplicam formulele (2.36)-(2.37) si, plecand de la punctele x⁽⁰⁾ = 0, y⁽⁰⁾ = 0, obtinem:

x y⁽¹⁾ = −0.3333333333

x y⁽²⁾ = −0.4444444443

x y⁽³⁾ = −0.4814814813

Pentru comparatie, determinam solutia exacta a sistemului considerat:

x = , y =

adica x = −0.125 , y = −0.5.

Implementare

A. Algoritm

Observatia. Deoarece metoda Gauss-Seidel este o metoda iterativa, trebuie specificata o conditie de oprire a algoritmului. In continuare vom folosi aceeasi conditie de oprire a algoritmului ca si cea prezentata in paragraful 2.6.3:

(2.38)

Date de intrare: un sistem de ecuatii (o multime de ecuatii), un punct initial, x⁽⁰⁾, o eroare ε.

Date de iesire: solutia aproximativa a sistemului, obtinuta in urma aplicarii traiectoriei Gauss-Seidel vectorului x⁽⁰⁾ pana cand este indeplinita conditia (2.38).

Algoritmul consta in urmatoarele etape:

1. generarea matricei A a sistemului si a vectorului b

. n - numarul de necunoscute (numarul de linii ale matricei A)

2. generarea matricelor L, D, U si verificarea convergentei metodei:

(−(L + D)⁻¹U) < 1

3. constructia traiectoriei Gauss-Seidel

repeta

pana cand

B. Programe MAPLE si rezultate

gseidel := proc (eq::set(equation), init::vector, eps::float)

local var, n, AA, A, b, l, d, u, i, j, m, lst, xo, test, k, x;

var := [op(indets(eq))];

n := nops(var);

if vectdim(init) <> n then

ERROR('numarul de necunoscute nu este egal cu dimensiunea vectorului initial')

fi;

AA := genmatrix(eq, var, flag);

A := delcols(AA, n+1 .. n+1);

b := col(AA, n+1);

l := matrix(n, n, 0):

u := matrix(n, n, 0):

d := matrix(n, n, 0):

for i from 1 to n do

for j from 1 to i-1 do

l[i, j] := A[i, j];

od;

d[i, i] := A[i, i];

for j from i+1 to n do

u[i, j] := A[i, j];

od;

od;

# conditia de convergenta

m := multiply(inverse(matadd(l, d)), -u);

lst := [eigenvals(m)];

if evalf(max(seq(abs(lst[k]), k=1..nops(lst))))>=1 then

ERROR('Algoritmul nu converge');

fi;

# algoritmul propriu-zis

for i from 1 to n do

xo[i] := init[i]

od;

test := 1;

while test>=evalf(eps*sqrt(n)) do

x[1] := evalf(1/A[1, 1]*( b[1]-sum(A[1, k]*xo[k], k=2..n)));

for i from 2 to n do

x[i] := evalf(1/A[i, i]*( b[i]-sum(A[i, k]*x[k], k=1..i-1)-sum(A[i, k]*xo[k], k=i+1..n)));

od;

test := evalf(sqrt(sum((x[k]-xo[k])^2, k=1..n )));

for i from 1 to n do

xo[i] := x[i];

od;

od;

RETURN(seq(var[i]=x[i], i=1..n));

end:

debug(gseidel):

gseidel(, vector(2, [0, 0]), 0.01);

, array(1 ..2, [(1)=0, (2)=0]), .1e-1

var := [x, y]

n

b

d u_1,2 := 1 l_2,1 := 1 d_2,2 := 2

xo xo₂ := 0

test

x x₂ := 1.666666666 test := 2.357022604

xo xo₂ := 1.666666666

x x₂ := 1.944444444

test

xo xo₂ := 1.944444444

x x₂ := 1.990740740

test

xo xo₂ := 1.990740740

x x₂ := 1.998456790

test

xo xo₂ := 1.998456790

x x₂ := 1.999742798

test

xo xo₂ := 1.999742798

<-- exit gseidel (now at top level) = x = 1.000514403, y = 1.999742798}

x = 1.000514403, y = 1.999742798

Comparativ, prezentam rezultatele obtinute cu ajutorul metodei Gauss-Seidel pentru acelasi sistem si vector initial al traiectoriei, dar pentru diferite valori ale erorii:

> gseidel( , vector(2, [0,0]), 0.01);

x = 1.000514403, y = 1.999742798

> gseidel( , vector(2, [0,0]), 0.001);

x = 1.000085734, y = 1.999957133

> gseidel( , vector(2, [0,0]), 0.00001);

x = 1.000002381, y = 1.999998810

Programul pentru metoda Gauss-Seidel prin limbajul C++

# include<iostream.h>

# include<math.h>

# include<conio.h>

int n, i, j, t;

double a[10][10], b[10], x[10], xa, s, eps=0.00001;

void main (void)

cout<<"Introduceti termenii liberi;"<<endl;

for (i = 1; i <= n; i++)

i = 1;

t = 0;

do

while ((i <= n)&&(t = = 0));

if (t = = 1)

cout<<"Metoda nu converge!"<<endl;

else

do

s = sqrt (s);

}

while (s>=eps);

cout<<"Solutia sistemului:"<<end;

for (i = 1; i <= n; i++)

cout<<"x("<<i<<")="<<x[i]<<endl;

}

getch ( );

}