Matematica
Metoda lui GaussMetoda lui Gauss 1 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 xk din ecuatiile k + 1, , n, prin inmultirea ecuatiei k cu (elementul akk se numeste pivot) si adunarea acestora la ecuatia i (i > k). In functie de alegerea pivotului, exista urmatoarele variante ale metodei lui Gauss: 1. metoda lui Gauss clasica - in care la fiecare pas, pivotul este elementul akk, ; 2. metoda lui Gauss cu semipivot - in care la fiecare pas, se alege ca pivot elementul aik maxim in valoare absoluta pe coloana, pentru i > k, permutandu-se linia k cu linia i; 3. 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 (2.1) se reduce la forma superior triunghiulara (2.2) iar rezolvarea sistemului (2.2) se face prin substitutie inversa:
(2.3) , k = n-1, n-2, ., 1 Observatia 1. Cu ajutorul eliminarii gaussiene se poate determina si inversa unei matrice. Redam in continuare algoritmul de aflare a inversei unei matrice A. 1. generarea matricei B prin concatenarea matricelor A (de dimensiune n) cu matricea In. 2. pentru
pentru
pentru ,
pentru
3. prin stergerea primelor n coloane ale matricei B astfel transformate, se obtine inversa matricei A. 2 Probleme rezolvate Exercitiul 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 Li , iar coloana j cu Cj . 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
Pasul 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
Pasul 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 = 2. Exercitiul 2. Sa se gaseasca inversa matricei . Rezolvare Consideram matricea B obtinuta prin concatenarea matricei A cu matricea unitate I3:
Folosind metoda eliminarii a lui Gauss, transformam matricea B dupa cum urmeaza:
InversamatriceiAvafimatriceaC,obtinutaprin stergereaprimelor3 coloane ale matricei B:
Intr-adevar, se verifica usor ca A · C = C · A = I3. 3 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:
pentru daca,atuncisecautarpentru 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 |akr | > |akk |, , 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|apr | > |akk |, 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 imple-mentat separat cele trei variante de eliminare, folosind procedurile cgauss, spgauss tpgauss. Aceste proceduri vor fi folosite apoi ca optiuni in proce-dura 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 2. Instructiunea indets(set_eq) returneaza multimea nede-terminatelor sistemului set_eq. Deoarece ordinea elementelor acestei multimi nu este neaparat aceeasi cu ordinea nedeterminatelor din prima ecuatie a sis-temului, 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 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 3
x3 := 3 s := 0 s := 3 x2 := 2 s := 0 s := 2 s := 5 x1 := 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 = 3
x3 := 3 s := 0 s := x2 := 2 s := 0 s := −2 s := 7 x1 := 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 := 3
l := [y, z, x] x3 := 1 s := 0 s x2 := 3 s := 0 s := 3 s := 4 x1 := 2 <-- exit gauss (now at top level) = y = 2, z = 3, x = 1} y = 2, z = 3, x = 1 Observatia 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(%); [1, 2, 3] 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);
|