Metoda lui Gauss

Metoda lui Gauss

Breviar teoretic

Metoda lui Gauss presupune transformarea sistemului Ax = b intr-un sistem superior triunghiular si apoi rezolvarea acestuia prin substitutie inversa.

Constructia sistemului superior triunghiular se face astfel: la pasul k se elimina x_k din ecuatiile k + 1, , n, prin inmultirea ecuatiei k cu (elementul a_kk se numeste pivot) si adunarea acestora la ecuatia i (i > k).

In functie de alegerea pivotului, exista urmatoarele variante ale metodei lui Gauss:

metoda lui Gauss clasica - in care la fiecare pas, pivotul este elementul a_kk, ;

metoda lui Gauss cu semipivot - in care la fiecare pas, se alege ca pivot elementul a_ik maxim in valoare absoluta pe coloana, pentru i > k, permutandu-se linia k cu linia i;

metoda lui Gauss cu pivot total - in care la fiecare pas, se alege ca pivot elementul maxim atat pe linie, cat si pe coloana, pentru i > k, j > k, permutandu-se linia k cu linia i si coloana k cu coloana j;

In acest fel, sistemul (1) se reduce la forma superior triunghiulara

(2)

iar rezolvarea sistemului (2) se face prin substitutie inversa:

(3)

, k = n-1, n-2, ., 1

Observatia 1.1. Cu ajutorul eliminarii gaussiene se poate determina si inversa unei matrice. Redam in continuare algoritmul de aflare a inversei unei matrice A.

generarea matricei B prin concatenarea matricelor A (de dimensiune n) cu matricea I_n.

pentru

pentru

pentru ,

pentru

prin stergerea primelor n coloane ale matricei B astfel transformate, se obtine inversa matricei A.

Probleme rezolvate

Exercitiul 1.1. Sa se rezolve urmatorul sistem folosind cele trei variante ale eliminarii Gauss:

Matricea sistemului este

iar este matricea sa extinsa:

Deoarece numarul ecuatiilor este egal cu cel al necunoscutelor si det A = −3 ≠ 0, sistemul este compatibil determinat (de tip Cramer), si deci metoda eliminarii a lui Gauss este aplicabila.

In continuare, pentru a efectua operatiile asupra matricei extinse a sistemului vom nota linia i cu L_i , iar coloana j cu C_j .

Rezolvare utilizand metoda lui Gauss clasica

A. Constructia sistemului superior triunghiular

Pasul 1

pivot:

Pasul 2

pivot:

In acest moment am ajuns la un sistem de forma , echivalent cu sistemul initial, in care matricea este superior triunghiulara, unde:

, .

B. Rezolvarea sistemului superior triunghiular

Prin metoda substitutiei inverse, avem:

de unde obtinem solutia sistemului: x = 1, y = 2, z = 3.

Rezolvare cu metoda lui Gauss cu semipivot

A. Constructia sistemului superior triunghiular

Pasul 1

• Ca pivot se ia elementul de modul maxim de pe coloana 1. In cazul nostru, pivotul este a₁₂, deci se permuta linia 1 cu linia 2, si se fac zerouri pe coloana 1 pentru i > 1:

Pasul 2

• Ca pivot se ia elementul a_i2 de modul maxim de pe coloana 2, pentru i ≥ In cazul nostru, pivotul este a₃₂, deci se permuta linia 2 cu linia 3 si se fac zerouri pe coloana 2, pentru i > 2:

In acest moment am ajuns la un sistem de forma , echivalent cu sistemul initial, unde matricea este superior triunghiulara, iar:

,

B. Rezolvarea sistemului superior triunghiular se face ca si in cazul metodei lui Gauss clasice, si conduce la solutia x = 1, y = 2, z = 3.

Rezolvare cu metoda lui Gauss cu pivot total

A. Constructia sistemului superior triunghiular

Pasul 1

• ca pivot se alege elementul a_ij de modul maxim pentru i, j ≥ 1. In cazul nostru pivotul este a₃₂, deci se permuta linia 3 cu linia 1, si coloana 2 cu coloana 1:

• pentru corectitudinea rezultatului final este necesar ca, ori de cate ori se permuta coloanele matricei extinse, sa se permute si elementele corespunzatoare ale vectorului x. Astfel, avem:

• in final, obtinem:

Pasul 2

• ca pivot se alege elementul a_ij de modul maxim pentru i, j ≥ Deoarece pivotul este a₂₃, se permuta coloana 3 cu coloana 2:

In acest moment am ajuns la un sistem de forma , echivalent cu sistemul initial, unde matricea este superior triunghiulara, iar:

, .

B. Rezolvarea sistemului superior triunghiular se face ca si in cazul metodei lui Gauss clasice, si conduce la solutia x = 1, z = 3, y =

Exercitiul 1. Sa se gaseasca inversa matricei

Rezolvare

Consideram matricea B obtinuta prin concatenarea matricei A cu matricea unitate I₃:

Folosind metoda eliminarii a lui Gauss, transformam matricea B dupa cum urmeaza:

Inversa matricei A va fi matricea C, obtinuta prin stergerea primelor 3 coloane ale matricei B:

Intr-adevar, se verifica usor ca

A C = C A = I₃.

Implementare

A. Algoritm

Algoritmii pentru cele 3 metode sunt asemanatori, diferenta dintre ei aparand (asa cum se poate vedea si din exemplul rezolvat) in modul de rezolvare a eliminarii Gauss.

Date de intrare: un sistem de ecuatii (scris ca multime de ecuatii)

Date de iesire: solutia sistemului

Algoritmul consta din urmatoarele etape:

1. generarea matricei extinse a sistemului, A =
• n = numarul de ecuatii (numarul de linii ale matricei A);
2. a) eliminarea Gauss pentru metoda lui Gauss clasica

pentru

daca , atunci se cauta r pentru care , si se schimba linia k cu linia r;

daca toti , atunci se returneaza eroare;

pentru , , unde ;

pentru ,

b) eliminarea Gauss pentru metoda lui Gauss cu semipivot

pentru se cauta elementul de modul maxim pe linie, i.e. daca |a_kr | > |a_kk |, , se schimba linia k cu linia r

pentru , , unde ;

pentru ,    ;

c) eliminarea Gauss pentru metoda lui Gauss cu pivot total

pentru se cauta elementul de modul maxim pe linie, i.e. daca |a_pr | > |a_kk |, p, , se schimba coloana p cu coloana r si linia r cu linia k

pentru , , unde ;

pentru ,    ;

3. rezolvarea sistemului superior triunghiular prin substitutie inversa

- pentru

B. Programe MAPLE si rezultate

Deoarece diferitele variante ale metodei lui Gauss se deosebesc doar prin modul in care se realizeaza eliminarea Gauss, in cele ce urmeaza am implementat separat cele trei variante de eliminare, folosind procedurile cgauss, spgauss tpgauss. Aceste proceduri vor fi folosite apoi ca optiuni in procedura finala gauss

restart: with(linalg):

cgauss:=proc(A::matrix)

local A1, A2, n, k, r, i, m, j;

n:=rowdim(A);

A1:=A;

A2:=delcols(A1,n+1..n+1);

if(det(A2)=0) then ERROR('sistemul nu are solutie unica!') fi;

for k from 1 to n-1 do

if A1[k,k]=0 then

for r from k+1 to n

while A1[k,r]=0 do r=r+1 od;

if r>n then ERROR('sistemul nu are solutie unica!')

else A1:=swaprow(A1,k,r);

fi;

fi;

for i from k+1 to n do

m:=A1[i,k]/A1[k,k];

for j from k to n+1 do

A1[i,j]:=A1[i,j]-m*A1[k,j];

od;

od;

od;

RETURN(evalm(A1));

end:

spgauss:=proc(A::matrix)

local A1, A2, n, k, r, i, m, j, mx;

n:=rowdim(A);

A1:=A;

A2:=delcols(A1,n+1..n+1);

if(det(A2)=0) then ERROR('sistemul nu are solutie unica!') fi;

for k from 1 to n-1 do

mx:=k;

for r from k to n do

if (abs(A1[r,k])>abs(A1[k,k])) then mx := r

fi;

od;

if mx<>k then A1:=swaprow(A1,k,mx); fi;

for i from k+1 to n do

m:=A1[i,k]/A1[k,k];

for j from k to n+1 do

A1[i,j]:=A1[i,j]-m*A1[k,j];

od;

od;

od;

RETURN(evalm(A1));

end:

tpgauss:=proc(A::matrix)

local A1, A2, n, k, r, i, m, j, mx, my, l;

n:=rowdim(A);

A1:=A;

A2:=delcols(A1,n+1..n+1);

l:=[seq(i), i=1..n];

if(det(A2)=0) then ERROR('sistemul nu are solutie unica!') fi;

for k from 1 to n-1 do

mx:=k; my:=k;

for i from k to n do

for j from k to n do

if (abs(A1[i,j])>abs(A1[k,k])) then mx:=i; my:=j;

fi;

od;

od;

if mx<>k then A1:=swaprow(A1,k,mx); fi;

if my<>k then

A1:=swapcol(A1,k,my);

l:=subsop(k=l[my],my=l[k],l);

fi;

for i from k+1 to n do

m:=A1[i,k]/A1[k,k];

for j from k to n+1 do

A1[i,j]:=A1[i,j]-m*A1[k,j];

od;

od;

od;

RETURN(evalm(A1),l);

end:

gauss:=proc(eqn::set(equation), opt::symbol)

local A,A1,l,n,r,k,i,m,j,s,x,rez;

l:=[op(indets(eqn))];

n:=nops(l);

A:=genmatrix(eqn, l, flag);

if opt=clasic then A1:=cgauss(A);

elif opt=semipivot then A1:=spgauss(A);

elif opt=totalpivot then

rez:=tpgauss(A);

A1:=rez[1];

l:=[seq(l[rez[2][i]],i=1..n)];

else ERROR('optiunile sunt: clasic, semipivot sau totalpivot');

fi;

x[n]:=A1[n,n+1]/A1[n,n];

for i from n-1 by -1 to 1 do

s:=0;

for j from i+1 to n do

s:=s+A1[i,j]*x[j];

od;

x[i]:=1/A1[i,i]*(A1[i,n+1]-s);

od;

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

end:

Observatia 1. Instructiunea indets(set_eq) returneaza multimea nedeterminatelor sistemului set_eq. Deoarece ordinea elementelor acestei multimi nu este neaparat aceeasi cu ordinea nedeterminatelor din prima ecuatie a sistemului, pot aparea diferente intre rezultatele furnizate cu ajutorul codului MAPLE si rezultatele calculate pe hartie. Desi matricea sistemului generata cu ajutorul instructiunii indets nu este intotdeauna aceeasi cu matricea sistemului scrisa pe hartie, rezultatele furnizate de program vor fi aceleasi (eventual ordinea solutiilor va fi schimbata).

Observatia 1.3. Pentru a urmari executia unei proceduri, se foloseste instructiunea debug. In cazul programelor din exemplele de mai sus, se poate folosi urmatorul set de instructiuni:

debug(cgauss):

debug(spgauss):

debug(tpgauss):

debug(gauss):

Redam mai jos, cu titlu de exemplu, rezultatul urmaririi procedurilor gauss cgauss spgauss si tpgauss pentru acelasi sistem de ecuatii.

> gauss(,clasic);

, clasic

l [x, y, z]

n

x₃ := 3 s := 0 s := 3 x₂ := 2

s := 0 s := 2 s := 5 x₁ := 1

<-- exit gauss (now at top level) = x = 1, y = 2, z = 3}

x = 1, y = 2, z = 3

> gauss(, semipivot);

, semipivot

l [x, y, z]

n

x₃ := 3     s := 0 s := x₂ := 2

s := 0 s := −2 s := 7 x₁ := 1

<-- exit gauss (now at top level) = x = 1, y = 2, z = 3 }

x = 1, y = 2, z = 3

> gauss ( , totalpivot);

, totalpivot

l := [x, y, z]

n

l := [y, z, x]

x

s

s

x

s

s

s

x

<-- exit gauss (now at top level) y = 2, z = 3, x = 1}

y = 2, z = 3, x = 1

Observatia 1.4. Pachetul linalg furnizeaza procedurile gausselim si backsub. Astfel, procedura gausselim efectueaza eliminarea gaussiana cu pivot partial asupra unei matrice n m. Procedura backsub ia ca argument rezultatul procedurii gausselim si furnizeaza solutia sistemului. Astfel, pentru matricea din exemplul precedent, avem:

> A := matrix ([[1, 1, 1, 6], [2, -1, 3, 9], [1, 4, 1, 12]]);

> gausselim(A);

> backsub(%);

Redam programul Maple pentru determinarea inversei unei matrice, precum si un exemplu:

invmatrix := proc(A :: matrix)

local n, A1, i, B, j, m, m1, k, C;

n := rowdim(A);

if coldim(A) <> n then

ERROR('Matricea trebuie sa fie patratica!')

fi;

if det(A)=0 then

ERROR('Matricea trebuie sa fie nesingulara!')

fi;

A1: = matrix (n, n, 0);

for i from 1 to n do

A1 [i, i] :=1

od;

B := concat (A, A1);

for i from 1 to n do

m := B[i, i];

for j from 1 to 2_*n do B[i, j] := B[i, j] / m od;

for j from 1 to n do

if i<>j then

m1 := B[j, i];

for k from 1 to 2_*n do

B[j, k] = B[j, k] - m1_*B[i, k];

od;

fi;

od;

od;

evalm(B); # rezultatul eliminarii gaussiene

C := delcols(B,1..n); # inversa matricei A

end:

> A := matrix([[2, 1, 1], [2, -1, 3], [1, 4, 1]]):

> debug (invmatrix):

> C = invmatrix(A);

> multiply(A,C);

> multiply(C,A);

Programul pentru metoda eliminarii a lui Gauss cu limbajul C++

Programul prezentat rezolva un sistem de n ecuatii liniare cu n necunoscute compatibil determinate. Pentru reducerea erorilor de rotunjire pe care le face calculatorul, inainte de prelucrarea sistemului si a subsistemelor se inlocuieste ecuatia care trebuie sa contina pivotul cu o ecuatie din (sub)sistem astfel incat pivotul obtinut sa aiba valoare absoluta maxima.

Datele de intrare sunt: dimensiunea sistemului, coeficientii si termenii liberi.

# include<iostream.h>

# include<math.h>

# include<conio.h>

int n, i, j, k, p, t;

double a[10][11], x[10], max, s;

void main (void)

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

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

k = 1;

t = 0;

do

if (max = = 0)

t = 1;

else

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

for (j = k+1; j <= n+1; j++)

a[i][j] - = a[i][k]*a[k][j]/a[k][k];

k++;}}

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

if ((t = =1) | | (a[n][n] = = 0.0))

cout<<"Sistemul nu este compatibil determinat!"<<endl'

else

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

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

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

}

getch ( );

}