Metode numerice

Determinarea matricei inverse prin metoda Gauss-Jordan

LABORATOR 6

DETERMINAREA MATRICEI INVERSE PRIN METODA GAUSS-JORDAN

1. Consideraţii teoretice

Rezolvarea sistemelor de ecuaţii liniare constituie unul dintre subiectele larg abordate de calculul numeric. În parte, acest fapt se datorează numeroaselor clase de fenomene care fie au ca model matematic sisteme de ecuaţii algebrice liniare, fie generează modele de acest gen.

Metodele de rezolvare a sistemelor de ecuaţii liniare se pot grupa, în principal, în trei mari clase:

- metode ce utilizează calculul de determinanţi- metode directe- metode iterative

Din punct de vedere numeric, metodele care utilizează calcul de determinanţi sunt improprii procesării automate, fapt care face ca algoritmii numerici să ocolească acest tip de abordare. Celelalte două clase de metode, însă, sunt utilizate în ponderi similare.

Înainte de a trece la o analiză detaliată a metodelor amintite mai sus, să precizăm câteva specificităţi relative la ecuaţiile algebrice liniare. În scopul facilitării analizei, vom lua în considerare un sistem cu dimensiunea 2, fără a reduce, în felul acesta, generalitatea rezultatele obţinute. Fie astfel:

un sistem de două ecuaţii liniare cu două necunoscute. Fiecare din cele două ecuaţii poartă numele de formă liniară. O formă liniară se defineşte prin intermediul a două proprietăţi:

cu a constantă arbitrară.

2. Metode directe

Metodele directe permit obţinerea soluţiei exacte a sistemului considerat, făcând abstracţie de erorile de rotunjire şi trunchiere, folosind un număr finit de operaţii elementare.

În cazul metodelor directe, sistemul iniţial este adus la unul echivalent, cu matricea caracteristică mult mai simplă, de regulă de formă triunghiulară. Într-o astfel de situaţie, rezolvarea sistemului este imediată.2.1 Metoda Gauss

1

Metode numerice implementate C

Este cunoscută şi sub numele de metoda eliminării, fiind una dintre cele mai vechi şi mai utilizate metode în rezolvarea sistemelor de ecuaţii liniare şi în inversarea matricelor.

Metoda Gauss reduce sistemul de ordin n la forma superior triunghiulară în n-1

paşi. Dacă sistemul Ax=b se renotează A0x=b0, cu şi ,

atunci la pasul k se obţine sistemul echivalent Akx=bk, cu , .

Pasul k constă în eliminarea lui xk din ecuaţiile de rang k+i, i=1(n-k). Înaintea pasului k matricea Ak-1 are primele k-1 coloane aduse la forma superior triunghiulară, adică pentru i>j şi j=1,2,...,k-1. La pasul k, pentru ik+1 se scade din ecuaţia

de rang i pe cea de rang k înmulţită cu factorul . Primele k ecuaţii rămân

neschimbate. De remarcat faptul că:

(5.1)

Aceste elemente prezintă un rol deosebit de important, fiind numite pivoţi. Dacă nici unul dintre pivoţii obţinuţi în cei n-1 paşi nu este nul sau mai mic decât cea mai mică valoare reprezentabilă pe sistemul de calcul automat, se obţine în final sistemul superior triunghiular An-1x=bn-1. Dacă un pivot este nenul dar mult mai mic în

comparaţie cu numărătorul, factorul va fi afectat de erori mari de rotunjire,

fapt care va duce la pierderea preciziei calculelor din acest moment. Pentru a împiedica un astfel de comportament, se procedează la aşa numita operaţie de pivotare. În acest sens, la fiecare pas k, înaintea eliminării lui xk se alege coeficientul cel mai mare în modul al acestei variabile, relativ la toate ecuaţiile de rang ik în care ea intervine, adică:

(5.2)

Apoi, ecuaţia de rang k se interverteşte cu ecuaţia de rang s găsită mai sus. În felul acesta se asigură imposibilitatea utilizării unui pivot nul (dacă toţi consideraţi ar

fi nuli atunci rezultă ). Un pivot nul, fără posibilitate de înlocuire a lui cu un alt parametru nenul, indică faptul că matricea sistemului este singulară şi sistemul nu are soluţie unică sau este incompatibil.

Algoritmul de pivotare poate fi implementat uşor pe un sistem programabil, elementele matriceale şi termenii liberi succesivi fiind stocaţi în memorie în aceleaşi spaţii cu cei iniţiali (prin reînlocuire).

Pasul k al algoritmului Gauss este de forma:

1. Se determină pivotul după relaţia (5.2). Dacă ( fiind precizia calculelor – spre exemplu 10-p dacă se lucrează cu p zecimale), atunci matricea sistemului este singulară şi algoritmul este oprit.

2


2. Dacă , se intervertesc elementele şi , şi respectiv bs şi bk.

3. Pentru i=k+1,...,n se pune , . Pentru se pune

.

Uneori, dacă matricea sistemului este slab condiţionată, se apelează la fiecare pas k la aşa numita pivotare totală. În cadrul acestei proceduri se calculează:

(5.3)

În acest caz, pentru a aduce pivotul maxim în poziţia (k,k) trebuie intervertite atât liniile s şi k, cât şi coloanele p şi k ale matricei sistemului. În acest fel, însă, ordinea necunoscutelor se schimbă. Dacă se doreşte ca în final ordinea necunoscutelor să fie păstrată, este necesar ca după rezolvarea sistemului să se efectueze o procedură de reobţinere a ordinii soluţiilor conform cu ordinea iniţială a necunoscutelor sistemului. Pentru aceasta, va fi necesar ca la fiecare interschimbare de coloane să se memoreze indicii p corespunzători fiecărui pas k.

Să pătrundem mai adânc în mecanismul algoritmului Gauss. Pentru aceasta să considerăm, fără a reduce generalitatea analizei, un sistem particular, de ordin 3, de forma:

Presupunem că matricea este nesingulară, deci . Atunci, sistemul are soluţie unică ce poate fi obţinută prin algoritmul Gauss.

Algoritmul Gauss implică o serie de paşi succesivi, fiecare dintre aceştia fiind compus din două operaţii: una de normalizare şi alta de reducere.

Pasul 1

a) Normalizare

Se înmulţeşte prima ecuaţie cu . Celelalte ecuaţii rămân neschimbate.

Sistemul devine:

unde

b) Reducere

Se înmulţeşte prima ecuaţie cu şi se adună cu a doua, rezultatul înlocuind cea de-a doua ecuaţie. Se înmulţeşte apoi prima ecuaţie cu şi se adună cu a treia, ultima fiind înlocuită cu rezultatul astfel obţinut. Sistemul devine:

3


unde

Pasul 2

c) Normalizare

Se înmulţeşte cea de-a doua ecuaţie cu . Sistemul devine:

unde

d) Reducere

Se înmulţeşte cea de-a două ecuaţie cu şi se adună cu a treia. Rezultatul înlocuieşte cea de-a treia ecuaţie. Sistemul devine:

unde

Observăm că, aşa cum precizează teoria, avem de procesat un număr de doi paşi (n-1 cu n=3). După încheierea acestei secvenţe de paşi, sistemul este adus la forma superior triunghiulară. Acum, determinarea celor trei necunoscute nu mai ridică nici o problemă:

2.2 Metoda Gauss-Jordan

Metoda Gauss poate fi adaptată pentru a opera asupra matricei extinse a sistemului:

Acest fapt se poate dovedi util în special când se implementează algoritmul într-un program rulabil pe un sistem de calcul, limbajele actuale de programare deţinând facilităţi de procesare a matricelor.

Operaţiile de eliminare gaussiană se realizează acum pe liniile matricei S. Aceste operaţii pot fi:

- înmulţirea unei linii cu o constantă- scăderea unei linii din alta şi înlocuirea celei de-a doua cu rezultatul scăderii- rearanjarea liniilor

4


În urma celor doi paşi, matricea extinsă devine:

Dacă se acţionează în continuare asupra acestei matrice pentru a forma în primele trei coloane matricea unitate, se obţine aşa numita metodă Gauss-Jordan. În acest caz coloana a patra va reprezenta tocmai vectorul soluţie:

Metoda Gauss-Jordan poate fi utilizată şi pentru calculul inversei unei matrice date. Astfel, dacă se cere inversa matricei:

atunci se formează matricea:

prin alăturarea matricei unitate la matricea A. Algoritmul Gauss-Jordan aplicat matricei C va transforma primele trei coloane ale acesteia în coloanele matricei unitate şi va avea ca efect apariţia în următoarele trei coloane a matricei inverse A-1.

Exemplu de aplicaţieSe dă matricea nesingulară:

Se cere să se determine matricea inversă acesteia, folosind algoritmul Gauss-Jordan.

RezolvareSe formează matricea:

Aplicând algoritmul Gauss-Jordan se ajunge în final la matricea:

Rezultă:

5


3. Implementarea algoritmului

Algoritmul descris se poate implementa sub limbajul de programare C conform cu listingul 1 prezentat în continuare.

Listingul 1

#include<stdio.h>#include<conio.h>#include<math.h>

void main(void){ int o,i,j,p=1,r,c,k,flag=0,grad,linii,coloane,t; double v[10][10],a[10][10],u[10][10],s[10][20],cod[10],linie[10][20],\ pivot[10],lin[10],col[10],sum[10][20],piv;

clrscr();

//Citirea matricei directe printf("Dati ordinul matricei de inversat:"); scanf("%d",&o); linii=o;coloane=2*o; for(i=0;i<o;i++) for(j=0;j<o;j++) { printf("v[%d,%d]=",i,j); scanf("%lf",&v[i][j]); }

//Afisarea matricei directe printf("\nMatricea directa este:"); for(i=0;i<o;i++) { printf("\n"); for(j=0;j<o;j++) { printf("\t%2.2f",v[i][j]); } }

//Construirea matricei unitare for(i=0;i<o;i++) for(j=0;j<o;j++) { if(i==j) u[i][j]=1; else u[i][j]=0; }

//Matricea extinsa printf("\n\nMatricea extinsa este:\n"); for(i=0;i<o;i++) { printf("\n"); for(j=0;j<2*o;j++) { if(j<o) { s[i][j]=v[i][j];

linie[i][j]=v[i][j]; if(i==j) pivot[i]=v[i][j];

6


printf("%2.2f\t",linie[i][j]); }

else { s[i][j]=u[i][j-o]; linie[i][j]=u[i][j-o]; printf("%2.2f\t",linie[i][j]); }

} }

for(i=0;i<o;i++) for(j=0;j<o;j++) { if(i==j) if(s[i][j]==0) p=10*p+i; }

c=p; i=0; while(p>=10) { r=p%10;cod[i]=r; p=p/10;i++; } k=i; printf("\n\nAfiseaza numarul de linie pe care pivotul este nul:");

if(i==0) printf("\nNU EXISTA PIVOT NUL!"); else for(i=k-1;i>=0;i--) { printf("\n\tlinia %1.0f",cod[i]); flag=1; }

printf("\nCodul de pivot:\t%d\n\n",c);

if(flag) goto stop; else { for(i=0;i<o;i++) { piv=linie[i][i]; for(j=0;j<2*o;j++) { linie[i][j]=(1/piv)*linie[i][j]; } for(k=i+1;k<o;k++) { piv=linie[k][k-1]; for(j=0;j<2*o;j++)

sum[i][j]=(-1)*piv*linie[i][j]; for(j=0;j<2*o;j++)

linie[k][j]=linie[k][j]+sum[i][j]; } }

//Matricea rezultanta printf("\nMatricea rezultanta pas 1 este:\n"); for(i=0;i<o;i++) {

7


printf("\n"); for(j=0;j<2*o;j++) printf("%2.2f\t",linie[i][j]); }

for(i=1;i<o;i++) { for(k=i-1;k>=0;k--) { piv=linie[k][i]; for(j=0;j<2*o;j++)

sum[i][j]=(-1)*piv*linie[i][j]; for(j=0;j<2*o;j++)

linie[k][j]=linie[k][j]+sum[i][j]; } }

//Matricea rezultanta printf("\n\nMatricea rezultanta pas 2 este:\n"); for(i=0;i<o;i++) { printf("\n"); for(j=0;j<2*o;j++) printf("%2.2f\t",linie[i][j]); } } goto end;stop: printf(”\nA fost intalnit un pivot nul!!!”);end: getch();}

Analiza programului conduce la câteva concluzii interesante. Astfel, o primă constatare ar fi aceea că în cazul furnizării în intrarea aplicaţiei a unei matrice singulare, programul generează o eroare de împărţire cu zero. Datele de ieşire pentru un astfel de caz, sunt prezentate mai jos.

Dati ordinul matricei de inversat:2v[0,0]=1v[0,1]=1v[1,0]=1v[1,1]=1

Matricea directa este: 1.00 1.00 1.00 1.00

Matricea extinsa este:

1.00 1.00 1.00 0.001.00 1.00 0.00 1.00

Afiseaza numarul de linie pe care pivotul este nul:NU EXISTA PIVOT NUL!Codul de pivot: 1

8


Floating point error: Divide by 0.Abnormal program termination

Într-adevăr, matricea este singulară, având determinantul nul.

În cazul în care cel puţin unul dintre pivoţi este nul, programul furnizează o ieşire de felul următor:

Dati ordinul matricei de inversat:2v[0,0]=1v[0,1]=0v[1,0]=-1v[1,1]=0

Matricea directa este: 1.00 0.00 -1.00 0.00


1.00 0.00 1.00 0.00-1.00 0.00 0.00 1.00

Afiseaza numarul de linie pe care pivotul este nul: linia 1Codul de pivot: 11

A fost intalnit un pivot nul!!!

Este de asemenea posibil ca nici unul dintre pivoţi să nu fie nul însă cel puţin unul să fie de valoare mică. În acest caz programul se comportă astfel:

Dati ordinul matricei de inversat:2v[0,0]=0.0000001v[0,1]=2v[1,0]=1v[1,1]=1



0.00 2.00 1.00 0.001.00 1.00 0.00 1.00

9



Matricea rezultanta pas 1 este:

1.00 20000000.00 10000000.00 0.00-0.00 1.00 0.50 -0.00


1.00 0.00 -0.50 1.00-0.00 1.00 0.50 -0.00

Matricea directă este afişată ca fiind deoarece în formatul de

afişare a datelor a fost făcută specificaţia 2.2%f, dorindu-se ca elementele matriceale să apară scrise doar cu primele două cifre zecimale. Cu toate acestea, în memorie este stocată, pentru elementul de pe prima linie şi prima coloană, valoarea 0.0000001. Dacă nu ar fi aşa, s-ar afişa mesajul erorii de împărţire prin zero.

Semnificativ aici este că putem sesiza faptul că operaţiile de calcul matematic s-au desfăşurat la limită, observând apariţia unui semn minus în faţa unui element nul al blocului matriceal unitate din matricea finală. Un astfel de semn va trebui să creeze suspiciune asupra corectitudinii rezultatului furnizat de program. Într-adevăr, calculând “la mână” produsul matricei directe cu matricea rezultat, declarată ca matrice inversă, constatăm că nu mai regăsim matricea unitate, aşa cum ar fi trebuit. Şi mai simplu, putem verifica eroarea produsă lansând încă o dată programul, având în intrare acum matricea de ieşire din cazul de mai sus. Evident, nu vom putea utiliza matricea de ieşire aşa cum este ea afişată, deoarece am obţine mesajul că avem un pivot nul. De aceea, mai întâi va trebui să rulăm programul afiţând un număr mai mare de cifre zecimale, de exemplu 8.

Dati ordinul matricei de inversat:2v[0,0]=0.0000001v[0,1]=2v[1,0]=1v[1,1]=1



0.00 2.00 1.00 0.001.00 1.00 0.00 1.00


10



1.00 20000000.00 10000000.00 0.00-0.00 1.00 0.50 -0.00


1.00000000 0.00000000 -0.50000003 1.00000005-0.00000000 1.00000000 0.50000003 -0.00000005

Acum putem introduce ca intrare matricea de ieşire de mai sus, găsind:

Dati ordinul matricei de inversat:2v[0,0]=-0.50000003v[0,1]=1.00000005v[1,0]=0.50000003v[1,1]=-0.00000005

Matricea directa este: -0.50000003 1.00000005 0.50000003 -0.00000005


-0.50 1.00 1.00 0.000.50 -0.00 0.00 1.00



1.00 -2.00 -2.00 -0.000.00 1.00 1.00 1.00


1.00000000 0.00000000 0.00000010 1.999999980.00000000 1.00000000 1.00000000 1.00000000

Concluzia este că va trebui să fim extrem de atenţi la orice informaţie ne poate furniza programul.

4. Teme de studiu

11


1. Să se scrie un program C care să ceară introducerea de la tastatură a dimensiunii unei matrice pătratice şi să calculeze matricea inversă prin metoda Gauss-Jordan. Matricele alese ca exemplu de rulare vor fi de dimensiune rezonabilă astfel încât matricea inversă corespunzătoare să poată fi determinată fără efort, pe cale manuală.

2. Să se facă observaţii asupra problemelor care pot să intervină în cazul acestei metode de inversare matriceală.

3. Să se explice necesitatea introducerii în program a variabilei piv.4. Să se explice mecanismul de determinare a codului de pivot şi semnificaţia

acestuia.5. Să se experimenteze funcţionarea programului cu următoarele intrări:

; ; ; ;

; ;

6. Să se imagineze o metodă de determinare a singularităţii matriceale prin calcul automat.

12

Metode numerice

Documents

Transcript of Metode numerice