CALCUL MATRICEALINVERSARE. FACTORIZARE

Tema la disciplina Matematici asistate de calculator

1. Partea teoretica

Chestiuni generale

Definirea si proprietatile inversei:

Matricea inversa A‾¹ a matricei A(A e o matrice patratica reala nesingulara de ordinul n), se defineste ca fiind acea matrice care satisface:

A* A‾¹ = A‾¹*A=I (I este matricea unitate de ordinul n)

A‾¹=(1/detA)*A˚

unde: - detA = numarul real nenul, numit determinantul matricei A

A˚ = matricea adjuncta. A˚ este transpusa matricei obtinute prin inlocuirea elementelor lui A cu complementii lor algebrici(cofatorii lor); cofactirul lui a(i,j) este minorul inmultit cu (-1)^(i+j).

Modalitatea de calcul a inversei este greu algoritmizabila si in plus necesita un volum mare de calcule.

O matrice A este o matrice nesingulara in cazul in care detA≠ 0.

Proprietati ale operatiei de inversare:

det A‾¹=1/detA;

(A*B) ‾¹=B‾¹*A‾¹;

(A‾¹)=(a)‾¹; (a = transpusa matricei A)

In Matlab, inversa unei matrici se calculeaza cu functia "inv", cu apelul : Y=inv(X).

Functia se utilizeaza des pt rezolvarea sistemelor de ecuatii liniare Ax=b. Solutia sistemului este x=inv(A)*b.

Exemplu:

Calculati inversa matricei: A=[ 2 1 ]

[ -1 1 ]

In Matlab: A=[2 1;-1 1]

Se apeleaza B=inv(A) => B=[0.3333 -0.3333]

[0.3333 0.6667];

Inversarea prin metoda completarii matricei unitate

Avand o matrice initiala, A, de dimensiune n x n se cauta gasirea inversei acesteia.

Pasul 1: Se completeaza la A matricea unitate. Va rezulta o matrice de dimensiune n x 2n

a₁₁ a₁₂ a₁₃ ....a_1n a₁₁ a₁₂ a₁₃ ....a_1n 1 0 0 ...0

a₂₁ a₂₂ a₂₃ ....a_2n a₂₁ a₂₂ a₂₃ ....a_2n 0 1 0....0

A= ........ B= .........0 0 1...0

a_n1 a_n2 a_n3 ...a_nn a_n1 a_n2 a_n3 ...a_nn 0 0 0...1

Pasul 2: Se cauta un element a_ii ≠ 0. Daca nu se gaseste, se cauta un element

a_{t i} ≠ 0 i=1,n ; t=i+1,n

Pasul 3: Daca nu se gaseste a_{t i} ≠ 0, matricea A este neinversabila, deci algoritmul se opreste. Daca se gaseste a_ii ≠ 0 sau a_{t i} ≠ 0, se schimba intre ele liniile t si i di matricea extinsa, B.

Pasul 4: Se imparte linia i din matricea E la elementul a_ii

Pasul 5: Pentru j ≠ i calculaeaza coeficientii noi:

a_{j t}=a_{j t}^vechi - a_{i t}*a_{j i} i=1,n ; j=1,n ; t=1,2n

Pasul 6: Coeficientii a_{i j} cu i=1,n si j=n+1,2n reprezinta coeficientii matricei unitate.

Inversarea unei matrici folosind algoritmul GAUSS- JORDAN

Algoritmul Gauss-Jordan se bazeaza pe o teorema din algebra matriceala: Daca o matrice nesingulara A poate fi redusa la matricea unitate I prin inmultirea la stanga cu un sir de matrice, atunci inversa A^-1 se poate calcula prin inmultirea lui I la stanga cu acelasi sir de matrice in ordine inversa.

Algoritmul de inversare cuprinde 2 procese de calcul derulate in paralel si are n pasi:

I. Se fac initializarile:

A =A

D =I

II. La pasul k, k=1, 2, ., n, se calculeaza elementele matricelor A^k si D^k, utilizand formulele:

a_kj^k = a_kj^k-1 / a_kk^k-1 j = k+1,n

d_kj^k = d_kj^k-1 / a_kk^k-1 j = 1,k

a_ii^k = 1 i = 1,k

d_ii^k =1    i= k+1,n

a_ij^k = a_ij^k-1 - a_ik^k-1 a_kj^k i=1,n i ≠ k j= k+1,n

d_ij^k = d_ij^k-1 - a_ik^k-1 d_kj^k i=1,n i ≠ k j= 1,k

a_ij^k = 0 i=1,n i ≠ j j= 1,k

d_ij^k = 0 i=1,n i ≠ j j= k+1,n

III. In final se obtin matricea unitate, inversa si valoarea determinantului:

I = Aⁿ

A = Dⁿ

Observatii:

1. Indicele superior corespunde pasului de calcul.

2. Daca elementul a_kk^k-1, numit pivot, are valoare nula sau modulul sau este sub un prag de zero

prestabilit, rezulta probleme la efectuarea impartirilor. Pivot nul nu implica neaparat ca matricea este singulara; poate fi adus in pozitia de pivot orice element a_ij^k-1, cu schimbarea intre elea liniilor k si i si a coloanelor k si j. Daca toate aceste posibilitati duc la esec rezulta ca matricea A este singulara.

3. Relatiile de mai sus sunt echivalente cu urmatoarele relatii, care vor fi implementate in matlab:

a_kj^k = a_kj^k-1 / a_kk^k-1 j = 1,n

d_kj^k = d_kj^k-1 / a_kk^k-1 j = 1,n

a_ij^k = a_ij^k-1 - a_ik^k-1 a_kj^k i=1,n j=1,n i ≠ k

d_ij^k = d_ij^k-1 - a_ik^k-1 d_kj^k i=1,n j=1,n i ≠ k

cu urmatoarele interpretari:

- liniile corespunzatoare pivotului (liniile k) ale matricelor A^k si D^k se obtin prin impartirea liniilor matricelor A^k-1, respectiv D^k-1, corespunzatoare pivotului, la pivot (a_kk^k-1);

- liniile i, cu i = 1,n , i ≠ k, ale matricelor A^k si D^k se obtin adunand la liniile i ale matricelor A^k-1, respectiv D^k-1, liniile k ale matricelor A^k si D^k inmultite cu (-a_ik^k-1) .

Factorizarea LU

Prin factorizarea LU, o matrice patratica este descompusa sub forma produsului a doua matrice triunghiulare: una inferior triunghiulara, cu elemente 1 pe diagonala principala(matricea L) si cealalta superior triunghiulara (matricea U). In Matlab, factorizarea LU este utilizata pentru obtinerea inversei unei matrice cu functia 'inv' si pentru calculul determinantului cu functia 'det'. Este totodata baza rezolvarii ecuatiilor liniare prin "impartirea matricelor" obtinuta cu operatorii '' si '/'.

Factorizarea LU a unei matrice se face cu functia lu ; se apeleaza cu una dintre sintaxele:

* [L,U]=lu(X) - returneaza o matrice superior triunghiulara U si o matrice inferior triunghiulara permutata L, astfel incat X=L*U;

* [L,U,P]=lu(X) - returneaza o matrice superior triunghiulara in U, inferior triunghiulara in L, si permutarea matriceala in P: L*U=P*X.

Exemplu:

Sa se factorizeze LU matricea A=[ 1 2 3]

[ 4 5 6]

[ 7 8 0]

Executand in Matlab:

A=[1 2 3;4 5 6;7 8 0];

[L,U]=lu(A);

se obtin rezultatele:

L= 0.1429 1.0000 0 U= 7.0000 8.0000 0

0.5714 0.5000 1.0000 0 0.8571 3.0000

1.0000 0 0 0 0 4.5000

L este o permutare a matricei inferior triunghiulara, iar U este o matrice superior triunghiulara.

Executand in Matlab:

[L,U,P]=lu(A)

vor rezulta urmatoarele:

L=    U= P=

Inversa unei matrice, X=inv(A) poate fi calculata si cu expresia:

X=inv(U)*inv(L)

iar determinantul:

D=det(A);

poata fi calculat si cu expresia: D=det(L)*det(U)

Rezolvarea sistemelor de ecuatii Ax=b prin factorizarea LU presupune solutionarea succesiva a sistemelor:

y=Lb x=Uy

Factorizare QR

Factorizarea QR este o descompunere a unei matrice ca produs a unei matrice ortonormale Q, cu o matrice superior triunghiulara R.

Factorizarea "QR" se utilizeaza pt rezolvarea sistemelor de ecuatii liniare cu mai multe ecuatii decat necunoscute. Cea mai buna solutie a unor astfel de sisteme, in sensul celor mai mici patrate, este calculata cu instructiunea: x=Ab, care utilizeaza aceasta factorizare. Prin utilizarea factorizarii QR, solutia x este calculata in doi pasi:

y=Q'*b; x=Ry

Descompunerea ortogonal-triunghiulara(factorizarea QR) se realizeaza cu functia "qr"; se apeleaza cu una dintre sintaxele:

[Q,R]=qr(X) - returneaza matricea triunghiular superioara R de aceeasi dimensiune cu X si matricea Q, astfel incat X=Q*R;

[Q,R,E]=qr(X) - returneaza matricea permutata E a matricei superior triunghiulara R, cu elementele diagonalei descrescatoare si matricea Q, astfel incat X*E=Q*R

Exemplu:

Sa se determine descompunerea QR a matricei: A=[1 2 3]

[4 5 6]

[7 8 0]

Executand in Matlab:

A=[1 2 3;4 5 6;7 8 0];

[Q,R]=qr(A)

se obtin rezultatele:

Q=    R=

Factorizarea Cholesky

Factorizarea Choleskey este o metoda directa de rezolvare a unui sistem de ecuatii liniare de forma:

Ax=b; cu A matrice pozitiv definita.

Urmatoarele afirmatii sunt echivalente:

matricea A este pozitiv definita

A= a si det(A)>0 (a este notatia matricei transpuse);

exista o matrice unica R, superior triunghiulara, nesingulara, care satisface relatia: A=R*R

Factorizarea Choleskey consta in descompunerea matricei coeficientilor sistemului conform pct. 3 unde R este numit factorul Cholesky. Sistemul initial devine: R*R*x=b, iar solutia acestuia presupune rezolvarea succesiva a sistemelor: Ry=b si apoi Rx=y.

Factorizarea Cholesky se calculeaza cu functia Matlab 'chol', se apeleaza cu una dintre sintaxele:

R=chol(A) [R,p]=chol(A);

unde: A - este matrice pozitiv definita;

R - este o matrice superior triunghiulara, astfel incat R*R=A;

p - este un scalar de test, egal cu zero daca matricea A este pozitiv definita si un intreg pozitiv in caz contrar.

Rezolvarea in Matlab a unui sistem de ecuatii liniare presupune etapele:

se calculeaza factorul Cholesky, cu relatia : R=chol(A);

se rezolva sistemul Ry=b cu relatia: y=Rb;

se rezolva sistemul Rx=y cu relatia: x=Ry;

Aceste etape se pot scrie compact sub forma:

R=chol(A); x=RRb

Exemplu:

Calculati factorul Cholesky al matricei:

A=[ 9 1 -1 0]

[ 1 9 0 -1]

[-1 0 9 1]

[ 0 -1 1 9]

Executand in Matlab secventele:

A=[9 1 -1 0;1 9 0 -1;-1 0 9 1;0 -1 1 9];

[R,p]=chol(A);

se obtine rezultatul:

R=

0 2.9814 0.0373 -0.3354

0 0 2.9812 0.3396

0 0 0 2.9618

p= 0 (matrice pozitiv definita)

2. Exemple

2.1. Exemplu in cazul implementarii inversei matriciale, pentru metoda clasica:

consideram matricea patratica A(de dimensiuni 3/3):

aplicam algoritmul implementat si rezulta iteratiile:

la inceput se calculeaza transpusa care va fi

At =

apoi utilizand o matrice temporara Actemp, vom determina minorii matricei:

Actemp = 1 1

1 6

Actemp = 1 1

1 3

(doar cativa)

iar apoi folosindu-se de Actemp va calcula conjugata matricei A dupa formula Ac(i,j)=(-1)^(i+j)*det(Actemp) i,j=1:n

Ac =

in final utilizanduse formula inversei clasice:

Ai=(1/det(A))*Ac va rezulta inversa:

Ai =

2.2. Exemplu pentru inversa matrice cu metoda gauss-Jordan:

Se initializeaza o matrice oarecare A:

>> A=rand(5)

A =

>> inversare_gauss(A)

deter = 0.0066 (se calculeaza determinantul matricii; daca determinantul e nul programul se opreste si se afiseaza un mesaj de eroare)

m = 5 (se determina nr de linii)

n = 5 (se determina nr de coloane)

L = (se initializeaza L cu matricea unitate)

I = (se initializeaza I cu matricea unitate)

k = 1 (se initializeaza k cu 1)

( Daca pivotul e diferit de zero (a_kk) se calculeaza elementele matricilor L,A,I pentru fiecare valoare a lui k pana la n )

( Calculul matricilor L, A e necesar pentru obtinerea matricii intermediare I, care ca fi in final inversa lui A )

L =    (matricile L,A,I pentru k=1, respectiv pivotul este a₁₁)

1.0000 0 0 0 0 Se folosesc formulele:

0.8984 1 0 0 0 L(i,k)=A(i,k)/A(k,k)

0 0 1 0 0 A(i,:)=A(i,:)-L(i,k)*A(k,:)

0 0 0 1 0 I(i,:)=I(i,:)-L(i,k)*I(k,:)

0 0 0 0 1.0000

A =

0 0.2343 -0.2445 0.0230 0.5624

I =

0 0 1.0000 0 0

0 0 0 1.0000 0

0 0 0 0 1.0000

U = (se initializeaza U cu elementele finale ale matricii A)

0 0.2343 -0.2445 0.0230 0.5624

0 0 -2.0543 0.4560 2.9872

0 0 0 -0.0340 0.0178

0 0 0 0 0.8684

Z = (se initializeaza Z cu matricea unitate)

k = 1

( Daca pivotul e diferit de zero (U_kk) se calculeaza elementele matricilor Z,U,I pentru fiecare valoare a lui k pana la n )

( Calculul matricilor Z, U e necesar pentru obtinerea matricii intermediare I, care ca fi in final inversa lui A )

Z = (matricile U,Z,I pentru k=1, respectiv pivotul este U₁₁)

1.0000 2.8691 0 0 0 Se folosesc formulele:

0 1.0000 0 0 0 Z(k,i)=U(k,i)/U(i,i)

0 0 1.0000 0 0 U(k,:)=U(k,:)-Z(k,i)*U(i,:)

0 0 0 1.0000 0 I(k,:)=I(k,:)-Z(k,i)*I(i,:)

0 0 0 0 1.0000

U =

0 0.2343 -0.2445 0.0230 0.5624

0 0 -2.0543 0.4560 2.9872

0 0 0 -0.0340 0.0178

0 0 0 0 0.8684

I =

I =    (se obtine matricea I dupa n iteratii)

d = (se copiaza in vectorul d elementele de pe diagonala principala a lui U)

I = (dupa fiecare incrementare j se calculeaza matricea intermediara

I,cu formula: I(j,:)=I(j,:)/d(j) )

elapsed_time =

0.0630

INVERSA MATRICEI:

I =

2.3. Exemplu pentru inversarea prin completarea cu matrice unitate

luam o matrice oarecare, A:

>> A=rand(5)

A =

Pasul 1: se completeaza matricea unitate la matricea A:

B= 0.3200 0.2679 0.8392 0.6299 0.0272 0 0 1 0 0

Pasul 2: se ia a₁₁=0.21 (≠ 0) si se sare la pasul 3, astfel:

Pasul 3: se sare peste

Pasul 4: se imparte prima linie a lui B la a₁₁=0.21

1 1.9252 3.1929 0.9677 2.1093 4.6729 0 0 0 0

B= 0.3200 0.2679 0.8392 0.6299 0.0272 0 0 1 0 0

Pasul 5: se calculeaza noii coeficienti ai matricei B :

a₂₁ = 0.6435- 1*0.6435=0. a₂₂=0.7446-1*0.6235=0.1211

asadar se parcurge algoritmul in continuare, dar datorita numarului mare de iteratii facute, ne vom opri aici.

In final va rezulta matricea B:

B= 0 0 1 0 0 0.1269 -5.3521 1.6627 0.2752 3.6063

iar coeficientii de pe linia 1->5 si coloana 6 -> 10 reprezinta coeficientii matricei inverse a lui A.

deci A^-1 = 0.1269 -5.3521 1.6627 0.2752 3.6063

1.8173 8.2791 -2.4390 -0.1018 -6.6582

in final programul va afisa:

elapsed_time =

iteratii =

305

inversa lui a este :

ans =

2.4. Exemplu pentru factorizarea prin metoda Doolitle

consideram matricea patratica A(de dimensiuni 3/3):

aplicam algoritmul implementat si rezulta iteratiile:

U = 1 0 0

U = 1 1 0

U = 1 1 1

L = 1 0 0

L = 1 0 0

U = 1 1 1

U = 1 1 1

L = 1 0 0

U = 1 1 1

L = 1 0 0

ultimele doua iteratii fiind si matricile rezultante

2.5. Exemplu pentru factorizarea prin metoda Crout

consideram matricea patratica A(de dimensiuni 3/3):

aplicam algoritmul implementat si rezulta iteratiile:

U = 1 0 0

U = 1 1 0

U = 1 1 1

L = 1 0 0

L = 1 0 0

U = 1 1 1

U = 1 1 1

L = 1 0 0

U = 1 1 1

L = 1 0 0

2.6. Exemplu pt factorizare prin metoda Choleskey :

Se genereaza o matrice simetrica si pozitiv definita A:

>> A=pascal(9)

A =

careia i se aplica algoritmul Cholesky.

>> fact_cholesky(A)

n = 9 (se determina marimea matricei A)

L =    (se initializeaza L cu matricea nula)

( Incepe bucla principala cu implementarea algoritmului Cholesky )

v = 1 (v,d sunt vectori intermediari necesari pentru calculul elementelor matricilor L si D)

d = 1 ( v(1) = A(1,1), d(1) = v(1) )

L = ( se calculeaza termenii de pe prima coloana a lui L cu formula:

L(2:n,1)=A(2:n,1)/v(1) )

se termina bucla corespunzatoare functiei else , activata pentru j=1

( In continuare la fiecare incrementare a lui j, se calculeaza vectorii intermediari v si d, necesari pentru calculul elemetelor matricilor L si D )

v = 1 (vectorii v si d pentru j=2)

Vectorii se calculeaza aplicand formulele:

v(1:j-1)=L(j,1:j-1).*d(1:j-1);

v(j)=A(j,j)-L(j,1:j-1)*v(1:j-1)';

d(j)=v(j);

v = 1 1

d = 1 1

L =    (se completeaza coloana 2 a matricii L)

Termenii lui L se calculeaza cu formula:

L(j+1:n,j)=(A(j+1:n,j)-L(j+1:n,1:j-1)*v(1:j-1)')/v(j);

v = 1 2 (vectorii v si d pentru j=3)

v = 1 2 1

d = 1 1 1

L =    (se completeaza coloana 3 a matricii L)

v = 1 3 3 (vectorii v si d pentru j=4)

v = 1 3 3 1

d = 1 1 1 1

L = (se completeaza coloana 4 a matricii L)

v = 1 4 6 4 (vectorii v si d pentru j=5)

v = 1 4 6 4 1

d = 1 1 1 1 1

L = (se completeaza coloana 5 a matricii L)

v = 1 5 10 10 5 (vectorii v si d pentru j=6)

v = 1 5 10 10 5 1

d = 1 1 1 1 1 1

L = (se completeaza coloana 6 a matricii L)

v = 1 6 15 20 15 6 (vectorii v si d pentru j=7)

v = 1 6 15 20 15 6 1

d = 1 1 1 1 1 1 1

L = (se completeaza coloana 7 a matricii L)

v = 1 7 21 35 35 21 7 ( vectorii v si d pentru j=8)

v = 1 7 21 35 35 21 7 1

d = 1 1 1 1 1 1 1 1

L = (se completeaza coloana 8 a matricii L)

v = 1 8 28 56 70 56 28 8 ( vectorii v si d pentru j=9)

v = 1 8 28 56 70 56 28 8 1

d = 1 1 1 1 1 1 1 1 1

elapsed_time = 0.0160 (se afiseaza timpul necesar aplicarii algoritmului)

( Se afiseaza matricea D implementata cu formula:D=diag(d) )

D =

( In final elementele de pe diagonala lui L sunt 1 )

L =

3. Prezentarea implementarii

Functia principala, de unde se acceseaza meniul:

function []=proiect()

optiune=' ';

while optiune~='x'

clc;

disp('c/C = Citire')

disp('i/I = Inversare')

disp('f/F = Factorizare')

disp('x/X = eXit')

optiune=input('Dati optiunea: ')

switch lower(optiune)

case 'c'

disp('Citirea matricei')

A=input('de exemplu [1 2;3 4] ' )

disp('press any key when ready')

pause

case 'i'

disp('invesare')

inversare(A)

disp('press any key when ready')

pause

case 'f'

disp('Factorizare')

factorizare(A)

disp('press any key when ready')

pause

case 'x'

disp('Iesirea')

otherwise

disp('COMANDA NECUNOSCUTA repetati comanda ')

disp('press any key when ready')

pause

end

end

Functia pt determinarea inversei:

function []=inversare(A)

opt1=' ';

while opt1~='x'

clc;

disp('a/A = Inversarea clasica')

disp('b/B = Inversarea prin metoda completarii matricii unitate')

disp('c/C = Inversarea prin metoda Gauss Jordan')

disp('x/X = eXit la meniul principal')

opt1=input('Dati optiunea: ')

[m,n] = size(A);

switch lower(opt1)

case 'a'

%!!!!!!Metoda intai ! ! ! ! ! ! ! ! ! ! ! !

B=A;

tic;

iteratii=0;

%calculul transpusei

if det(A)~=0

for i=1:m

for j=1:n

At(i,j)=A(j,i);

iteratii=iteratii+1;

end

iteratii=iteratii+1;

end

disp('Transpusa')

At

%calculul conjugatei

for i=1:m

for j=1:n

for i1=1:(m-1)

for j1=1:(m-1)

if (i1<i)&(j1<j) Actemp(i1,j1)=At(i1,j1);

end

if (j1>=j) Actemp(i1,j1)=At(i1,j1+1);

end

if (i1>=i) Actemp(i1,j1)=At(i1+1,j1);

end

if (i1>=i)&(j1>=j) Actemp(i1,j1)=At(i1+1,j1+1);

end

iteratii=iteratii+1;

end

iteratii=iteratii+1;

end

Ac(i,j)=(-1)^(i+j)*det(Actemp);

iteratii=iteratii+1;

end

iteratii=iteratii+1;

end

disp('Conjugata')

Ac

%calculul inversei

Ai=(1/det(A))*Ac;

timp1=toc;

disp('Inversa matricii este')

Ai

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp1

disp(' secunde ')

disp('press any key when ready')

pause

else

disp('Determinantul matricii este 0, nu se poate calcula inversa dupa metoda clasica. press amy key when ready');

pause

end

A=B;

case 'b'

%! !!!Metoda a doua ! ! ! ! ! ! ! ! ! ! !!!!!!

B=A;

temp2=1;

tic;

iteratii=0;

if det(A)~=0

for x=1:n

for j=1:n

iteratii=iteratii+1;

if x==j A(x,j+n)=1;

else A(x,j+n)=0;

end

end

end

for i=1:n

iteratii=iteratii+1;

if A(i,i)==0

for t=i+1:n

if A(t,i)~=0

auxiliar=A(i,:); A(i,:)=A(t,:); A(t,:)=auxiliar;

iteratii=iteratii+1;

end

end

if t==n

disp('matricea initiala nu se poate inversa')

temp2=0;

break

end

end

aux=A(i,i);

for t=1:2*n

A(i,t)=A(i,t)/aux;

iteratii=iteratii+1;

end

for j=1:n

if (j~=i)

aux=A(j,i);

for t=1:2*n

A(j,t)=A(j,t)-A(i,t)*aux;

iteratii=iteratii+1;

end

end

iteratii=iteratii+1;

end

end

timp2=toc;

if temp2==1

disp('inversa lui A este : ')

Ai=A(1:n,n+1:2*n)

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp2

disp(' secunde ')

disp('press any key when ready')

pause

end

else %if det(A)~=0

disp('Determinantul matricii este 0, nu se poate calcula inversa dupa metoda completarii matricii unitate');

pause

end

A=B;

case 'c'

%!!!!!!Metoda a treia ! ! ! ! ! ! ! ! ! ! !

B=A;

tic;

iteratii=0;

if det(A)~=0

tic;

U=A; I=A; B=A;

L=eye(m);

I=eye(m);

k=1;

while k<n;

if A(k,k)==0; % Pivot

disp('!!! EROARE - PIVOTUL NU POATE FI ZERO !!!')

break

else

for i=k+1:m;

L(i,k)=A(i,k)/A(k,k);

A(i,:)=A(i,:)-L(i,k)*A(k,:);

I(i,:)=I(i,:)-L(i,k)*I(k,:);

iteratii=iteratii+3;

end

k=k+1;

iteratii=iteratii+1;

end

end

U=A(:,1:n);

Z=eye(m);

k=1;

while k<n;

if U(k,k)==0;

disp('!!! EROARE - PIVOTUL NU POATE FI ZERO !!!')

break

else

for i=k+1:m;

Z(k,i)=U(k,i)/U(i,i);

U(k,:)=U(k,:)-Z(k,i)*U(i,:);

I(k,:)=I(k,:)-Z(k,i)*I(i,:);

iteratii=iteratii+3;

end

k=k+1;

iteratii=iteratii+1;

end

end

d=(diag(U))';

for j=1:n;

if d(j)==1;

j=j+1;

else

I(j,:)=I(j,:)/d(j);

j=j+1;

end

iteratii=iteratii+1;

end

timp3=toc;

disp('inversa lui A este : ')

I

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp3

disp(' secunde ')

disp('press any key when ready')

pause

else

disp('Determinantul matricii este 0, nu se poate calcula inversa dupa metoda Gauss Jordan');

pause

end

A=B;

case 'x'

disp('Am Iesit')

otherwise

disp('press any key when ready')

disp('COMANDA NECUNOSCUTA repetati comanda ')

pause

end

end

Functia pt determinarea factorizarilor:

function []=factorizare(A)

opt2=' ';

while opt2~='x'

clc;

disp(' Factorizari L U (lower triangle upper triangle)')

disp(' a Metoda Cholesky')

disp(' b Metoda Doolitlel ')

disp(' c Metoda Crout')

disp(' x eXit in meniul principal')

opt2=input('Dati optiunea: ')

[m,n] = size(A);

switch lower(opt2)

%! !!!metoda Cholesky! ! !

case 'a'

B=A;

temp=1;

for i=1:n %test de simetrie

for j=1:n

if A(i,j)~=A(j,i)

temp=0;

end

end

end

if temp==1

tic;

iteratii=0;

L=zeros(n,n);

for j=1:n

if (j > 1)

v(1:j-1)=L(j,1:j-1).*d(1:j-1);

v(j)=A(j,j)-L(j,1:j-1)*v(1:j-1)';

d(j)=v(j);

iteratii=iteratii+3;

if (j < n)

L(j+1:n,j)=(A(j+1:n,j)-L(j+1:n,1:j-1)*v(1:j-1)')/v(j);

iteratii=iteratii+1;

end

else

v(1)=A(1,1);

d(1)=v(1);

L(2:n,1)=A(2:n,1)/v(1);

iteratii=iteratii+3;

end;

end;

D=diag(d);

L=L+eye(n);

timp1=toc;

disp('matricea originala')

A

disp('matricea recompusa din L si L-transpus')

A1=L*(L')

disp('matricea triunghiulara')

L

disp('matricea triunghiulara transpusa')

L'

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp1

disp(' secunde ')

disp('press any key when ready')

pause

end

if temp==0

disp('!!!!!Imi pare rau nu pot efectua metoda de factorizare Cholesky')

disp('!!!!! deoarece matricea nu este simetrica')

disp('press any key when ready')

pause

end

A=B;

%! !!!metoda Doolitle! ! !

case 'b'

B=A;

tic;

iteratii=0;

L=zeros(n); U=L;

for k=1:n

L(k,k)=1;

for j=k:n

sum1=0;

for m=1:(k-1) sum1=L(k,m)*U(m,j)+sum1;

iteratii=iteratii+1;

end

U(k,j)=A(k,j)-sum1;

iteratii=iteratii+1;

end

for i=(k+1):n

sum2=0;

for m=1:(k-1) sum2=L(i,m)*U(m,k)+sum2;

iteratii=iteratii+1;

end

L(i,k)=(A(i,k)-sum2)/U(k,k);

iteratii=iteratii+1;

end

iteratii=iteratii+1;

end

timp2=toc;

disp('matricea originala')

A

disp('matricea recompusa din L si U')

A1=L*U

disp('matricea inferior triunghiulara')

L

disp('matricea superios triunghiulara')

U

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp2

disp(' secunde ')

disp('press any key when ready')

pause

A=B;

%! !!!metoda Crout! ! !

case 'c'

B=A;

tic;

iteratii=0;

L=zeros(n,n); U=L;

for k=1:n

sum1=0

for m=1:(k-1) sum1=L(k,m)*U(m,k)+sum1;

iteratii=iteratii+1;

end

L(k,k)=A(k,k)-sum1;

for j=k:n

sum2=0;

for m=1:(k-1) sum2=L(k,m)*U(m,j)+sum2;

iteratii=iteratii+1;

end

U(k,j)=(A(k,j)-sum2)/L(k,k);

iteratii=iteratii+1;

end

for i=(k+1):n

sum3=0;

for m=1:(k-1) sum3=sum3+L(i,m)*U(m,k);

iteratii=iteratii+1;

end

L(i,k)=(A(i,k)-sum3)/U(k,k);

iteratii=iteratii+1;

end

iteratii=iteratii+1;

end

timp3=toc;

disp('matricea originala')

A

disp('matricea recompusa din L si U')

A1=L*U

disp('matricea inferior triunghiulara')

L

disp('matricea superios triunghiulara')

U

disp('numarul de iteratii a fost')

iteratii

disp('timpu executiei a fost de ')

timp3

disp(' secunde ')

disp('press any key when ready')

pause

A=B;

case 'x'

disp('Am Iesit')

otherwise

disp('COMANDA NECUNOSCUTA repetati comanda ')

disp('press any key when ready')

pause

end

end

4. Conluzii

4.1. Pentru algoritmul inversei in prima varianta (metoda clasica) numarul de iteratii este mare. Timpul de executie este redus crescand odata cu dimensiunea matricii (de exemplu pentru o matrice random de 20 pe 20, timpul de executie a fost de 1.8280 secunde pe un sistem AthlonXP 2000+ Palomino cu 512 Mbram, acelasi sistem fiind folosit si in urmatoarele exemple)

Din punct de vedere al complexitatii algoritmului acesta este unul usor, el fiind bazat pe algoritmul clasic de calcul al inversei a unei matrici.

4.2. Pentru algoritmul prin metoda lui Gauss-Jordan observam ca numarul de pasi este cel mai mic, fata de celelalte metode, astfel constatind ca el e cel mai eficient, insa nu si din punct de vedere al timpului.

4.3 Algoritmul de completare cu matricea inversa e un algoritm care necesita foarte putin timp, cel mai putin dintre metodele de fata, el fiind insa impementat pt un numar foarte mare de pasii.

4.4.Pentru algoritmul Doolitle(bazat pe metoda de factorizare LU) timpul de executie este de 0.0310 secunde(acest timp este calculat in cazul unei matrice patratice 20 pe 20), numarul de iteratii fiind foarte redus.

4.5. Pentru algoritmul Crout(bazat pe metoda de factorizare LU) timpul de executie este de 0.0310 secunde(acest timp este calculat in cazul unei matrice patratice 20 pe 20), numarul de iteratii fiind foarte redus.

Se observa ca are acelasi timp de executie, datorita faptului ca au formula generala asemanatoare, avand astfel si numar de iteratii foarte reduse.

La metoda Doolitle matricea U are 1 pe diagonala principala, iar la metoda Crout matricea L are 1 pe diagonala principala.

4.6. Metoda Cholesky

O forma particulara a factorizarii Cholesky este descompunerea matricei initiale A dupa formula:

A = L D L^T

unde L este o matrice inferior triunghiulara, iar D este o matrice dreptunghiulara.

5. Functii utilizate

In proiectul de fata am avut de implementat metode de rezolvare a inversei unei matrice, precum si de factorizare a unei matrice:

In Matlab, inversa unei matrice se calculeaza cu functia: inv(A), A^(-1)

Conditia ca inversa sa poata fi calculata este ca matricea sa fie patratica si ca determinantul ei sa fie diferit de 0.

Pentru facorizare in Matlab se aplica functiile:

[L,U]=lu(A) (L=matrice inferior triunghiulara, U=matrice superior triunghiulara)

[Q,R]=qr(A) (Q=matrice ortogonala, R= matrice superior triunghiulara)

Conditia ca factorizarea sa se poata aplica este ca matricea sa fie patratica.

6. Bibliografie

Ghinea M., Fireteanu V., "Matlab. Calcul Numeric - Grafica - Aplicatii", Editura Teora,1995

Curteanu S., "Calcul numeric si simbolic in Mathcad", Editura MatrixRom, Bucuresti 2001

Precup R., Dragomir L., Blavitchi I., "Matematici asistate pe calculator. Aplicatii", Editura Politehnica, Timisoara, 2002

Kilyeni St., 'Metode numerice. Algoritme, programe Turbo Pascal, aplicatii in energetica', Editura Orizonturi Universitare, Timisoara, 2001

Internet(adresele sunt mentionate in directorul Bibliografie)