CALCUL MATRICIAL. Metoda de eliminare Gauss, Factorizarea LU, Metoda lui Jacobi, Metoda...

8/7/2019 CALCUL MATRICIAL. Metoda de eliminare Gauss, Factorizarea LU, Metoda lui Jacobi, Metoda Gauss-Seidel, Metoda d

1/8

Carmen-Sanda Georgescu, Tudor Petrovici, Radu PopaMetode numerice. Fisa de laborator nr. 8:

CALCUL MATRICIAL(I). Metoda de eliminare Gauss, FactorizareaLU, Metoda lui Jacobi, Metoda Gauss-Seidel, Metoda de rezolvare

a sistemelor liniare in sensul CMMP. Conditionarea sistemelorliniare

5.1. METODE DIRECTE DE REZOLVARE A SISTEMELOR DE ECUATII LINIAREUn sistem de n ecuatii liniare cu n necunoscute, scris sub forma matriceala A*X = Bunde matricea coeficientilor necunoscutelor A are proprietatea ca det(A)0, este un sistem liniar compatibilsi determinat. Acest sistem are o solutie unica X=A-1*B. Vom remarca ca formula care da solutia folosindA-1 nu este recomandata; la fel nu sunt recomandate in calculele numerice solutiile care pot fi obtinute curegula lui Cramer.5 .1 .1 .ME T ODA DE E L I MI NA R E GA U SSFie un sistem de ecuatii liniare A*X=B, unde ARn x n si BRn, unde toate submatricile lider principale alematricii A, adica A[k]=[ai j]1 i, j k sunt nesingulare . Se demonstreaza ca in acest caz matricea A poate fiadusa prin transformari elementare la o forma superior triunghiulara (elementele de sub diagonala suntegale cu zero). Algoritmul care aduce matricea A la forma superior diagonala se bazeaza pe aplicareasuccesiva a 3 operatii elementare: normalizarea= inmultirea unei ecuatii cu o constanta; reducerea=scaderea unei ecuatii din alta si inlocuirea celei de-a doua ecuatii cu rezultatul scaderii; rearanjarea=

schimbarea ordinii ecuatiilor. Se exemplifica acest algoritm pentru un sistem de 3 ecuatii cu 3 necunoscutex1, x2, x3: (sistemul 8.1)

=++

=++

=++

3333232131

2323222121

1313212111

bxaxaxa

bxaxaxa

bxaxaxa

(8.1)

Algoritmul implica urmatorii pasi succesivi:

Pasul 1 (normalizare): In ipoteza ca a110, se inmulteste prima ecuatie din (8.1) cu 11a1 (unde a11 senumeste pivot) si coeficientii noii ecuatii se noteaza i1a (i= 2, 3), respectiv 1b :

=++

=++

=++

3333232131

2323222121

13132121

bxaxaxa

bxaxaxa

bxaxax

(8.2)

Pasul 2 (reducere): Se inmulteste prima ecuatie din (8.2) cu 21a si se aduna la a doua ecuatie, apoirezultatul obtinut inlocuieste cea de-a doua ecuatie a sistemului. Similar, se inmulteste prima ecuatie din

(8.2) cu si se aduna la a treia ecuatie, apoi rezultatul obtinut inlocuieste cea de-a 3-a ecuatie a

sistemului:31a

=+

=+

=++

3333232

2323222

13132121

bxaxa

bxaxa

bxaxax

(8.3)

Pasul 3 (normalizare): Se inmulteste a 2-a ecuatie din (8.3) cu 22a1 si coeficientii noii ecuatii se noteaza

, respectiv :23a 2b

=+

=+

=++

3333232

23232

13132121

bxaxa

bxax

bxaxax

(8.4)

Pasul 4 (reducere): Se inmulteste a doua ecuatie din (8.4) cu 32a si se aduna la a 3-a ecuatie, apoirezultatul obtinut inlocuieste cea de-a 3-a ecuatie a sistemului - vezi sistemul (8.5):


2/8

=

=+

=++

3333

23232

13132121

bxa

bxax

bxaxax

(8.5)

==

=

31321211

32322

3333

xaxabx

xabx

abx

(8.6)

Pornind de la ultima ecuatie catre primain (8.5), se obtine solutia sistemului sub forma (8.6).Nota: Daca la o operatie de normalizare, pivotul corespunzator este nul (sau apropiat de zero), se

vor rearanja ecuatiile urmatoare din sistem, aducandu-se pe pozitia respectiva una din ecuatiile cu pivotnenul.REMARCI

1) Algoritmul Gauss poate fi aplicat pentru determinarea rangului unei matrici oarecare ARn x m2) Pentru rezolvarea sistemelor compatibile si determinate A*X = B dar si pentru determinarea

inversei unei matrici nesingulare ARn x n pot fi aplicate si alte variante ale metodei lui Gauss sianume metodele de tip Gauss-Jordan dupa schema

A B InTransformarielementare

. . . . . .

In X A-1In schema prezentata mai sus: A este matricea coeficientilor necunoscutellor, B este coloanatermenilor liberi, In este matricea unitate, X este solutia sistemului A*X=B si A-1 este inversa matricii A(in ipoteza ca matricea A este inversabila)3) Atat pentru metoda lui Gauss cat si pentru variantele Gauss-Jordan calculele por fi optimizate prinREGULA PIVOTULUI

Exemplul 1: Sa se rezolve sistemul liniar si inversa matricii coeficientilor sistemului prin metoda Gauss-Jordan si regula pivotului4*x1-x2+x3=74*x1-8*x2+x3=-21-2*x1+x2+5*x3=15

1 0 1/40 1 00 0 11/2

11/44

33/2

2/7 -1/28 01/7 -1/7 06/77 1/77 2/11

1 0 0

0 1 00 0 1

2

43

-47/308 3/77 -1/22

1/7 -1/7 06/77 1/77 2/11

4 -1 14 -8 1-2 1 5

7-2115

1 0 00 1 00 0 1

Regula pivotului

1 -1/4 1/44 -8 1-2 1 5

7/4-2115

1/4 0 00 1 00 0 1

-8*1-(1/4*4)=-71*1-(4*1/4)=01*1-(-2)*(1/4)=1/2

1 -1/4 1/40 -7 00 1/2 11/2

7/4-28

37/2

1/4 0 0-1 1 01/2 0 1

(-21*1)-(4*7/4)=-2815*1-(-2)*(7/4)=37/20*1-4*(1/4)=-1

1 -1/4 1/40 1 00 1/2 11/2

7/44

37/2

1/4 0 01/7 -1/7 01/2 0 1

1*1-4*0=10*1-(4*0)=00*1-(1/4)*(-2)=1/2

1 0 1/40 1 00 0 11/2

11/443

2/7 -1/28 01/7 -1/7 03/7 1/14 1

0*1-(-2)*0=01*1-(-2)*0=1

Regula pivotului: Se alege linia 1 (normalizata) drept linie pivot cu primul element desemnat drept pivot.Se inlocuiesc elementele de pe coloana 1 (sub pivot) cu 0. Celelalte elemente ale liniilor 2 si 3 secalculeaza cu regula pivotului; astfel pentru inlocuirea elementului 8 se formeaza dreptunghiul cuvarfurile opuse in 8 si pivot si se calculeaza dupa regula determinantului -8*1-(-1/4*4)=-7;pentrutoate celelalte elemente de pe linii 2 si 3 se calculeaza in acelasi mod.

5 .1 .2 . ME T ODA GA U SS C U FA C T OR I ZA R E T R I U NGH I U L A R


3/8

Metoda cu factorizare triunghiular, bazat pe procedeul de eliminare Gauss, reprezint o altmetod direct de rezolvare a sistemelor de ecuaii algebrice liniare de forma:

=

n

2

1

n

2

1

nn2n1n

n22221

n11211

b

b

b

x

x

x

aaa

aaa

aaa

KK

K

KKKK

K

K

, adic BXA = , (8.7)

unde Aeste matricea ptrat a coeficienilor sistemului, Xeste vectorul coloan al necunoscutelor

i Beste vectorul coloan al membrilor drepi (ambii vectori fiind n-dimensionali).

( nn )

Se tie c orice matrice ptrat A, de dimensiuni ( )nn , se poate exprima ca produsul a doifactori matriceali Li U, unde L=(lij) este o matrice inferior triunghiular, iar U=(uij) este o matrice superiortriunghiular, cu elementele diagonalei principale egale cu unitatea, adic:

ULA =

=+

+

+

100000

uu1000

uuuu10

uuuuu1

llll

0lll

00ll

000l

in1ii

n21i2i223

n11i1i11312

nnni2n1n

ii2i1i

2221

11

KK

KKKKKKKK

KK

KKKKKKKK

KK

KK

KK

KKKKKK

KK

KKKKKK

KK

KK

,

,

,

(8.8)Descompunerea matricei ptrate A, definit n (8.8), se mai numete descompunerea LU1.

nlocuind matricea Adin sistemul iniial (8.7) cu produsul LU, rezult:

= XA BXUL = (8.9)i, dac se noteaz prin Yun vector intermediar, definit de ecuaia matriceal:

XUY = , (8.10)atunci ecuaia matriceal (8.9) devine:

BYL = . (8.11)Algoritm: Sistemul (8.11) se rezolv n raport cu Y, prin substituie progresiv, rezultnd:

11

11

l

by = i apoi, n2iylb

l

1y

1i

1j

jijiii

i ,pentru, =

=

=

. (8.12)

n continuare, vectorul Y fiind cunoscut, se rezolv sistemul (8.10) prin substituie regresiv, rezultndastfel componentele soluiei X:

nn yx = i apoi, 11nixuyxn

1ij

jijii ),(pentru, == +=

. (8.13)

Pentru determinarea elementelor lij i uij ale factorilor matriceali, ordinea calculelor esteurmtoarea: Prima coloan din L: din produsele liniilor lui Lcu prima coloan din Ui echivalarea cu elementele lui Arezult:

n1ial1i1i ,,==

Prima linie a lui U: din produsele primei linii a lui Lcu coloanele de la 2 la nale lui Ui echivalarea cuelementele lui Arezult:

n2jlau 11j1j1 ,, == A doua coloan a lui L: din produsele liniilor de la 2la nale lui Lcu a doua coloan din Ui echivalarea cuelementele lui Arezult:

1 Aceste simboluri deriv din cuvintele scrise n englez: L = lower (inferior; de jos), U = upper (superior;de sus).


4/8

n2iulal 121i2i2i ,, == A doua linie a lui U: .a.m.d.Relaiile generale de calcul sunt urmtoarele:

n1kkiulal1k

1j

jkijikik ,,pentru, ==

=

n1kkjulal

1u

1k

1i

ijkikjkk

kj ,,pentru, =>

=

=

(8.14)

iar operaiile se execut n ordinea menionat mai sus.

octave#1> function X = lufactori(A,B)[N,N]=size(A);

X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));

C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;

if A(p,p)==0'A este singulara. Nu exista solutie unica'break

endiffor k=p+1:N

mult=A(k,p)/A(p,p);A(k,p) = mult;

A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endfor

endforY(1) = B(R(1));for k=2:N

Y(k)= B(R(k))-A(k,1:k-1)*Y(1:k-1);endforX(N)=Y(N)/A(N,N);for k=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endforendfunction

Exemplul 2: Sa se rezolve sistemul liniar prin metoda factorizarii LUx1+x2+5*x3=-72*x1+x2+x3=2x1+3*x2+x3=5Program Octave: octave#1>a=[1 1 5;2 1 1;1 3 1], b=[-7 2 5], X=lufactori(a,b)

Se obtine solutia: X=[ 1 2 -2]

5.2. METODE ITERATIVE DE REZOLVARE A SISTEMELOR DE ECUATII LINIARE

5.2 .1 . ME T ODA I T E R A T I VA JA C OB I ( ME T ODA I T E R A T I I L OR SI MU L T A NE )

Atentie! Condiia de convergen implic dominan diagonal n cadrul matricei Aadica(max|ai j/ai i|


5/8

= =

+n

ij1j

kjiji

ii

1ki xab

a

1x

)()(, K,,,;, 210kn1i == , (8.15)

corespunztoare metodei iterative Jacobi. Aceasta formula permite calcularea componentelor vectorului

soluie X(k+1)

, folosind componentele soluiei la iteraia precedent X(k)

. Vectorul iniial X(0)

este ales dectre utilizator.Calculul iterativ poate nceta la ndeplinirea urmtorului criteriu de convergen:

K,,,,,,max)()(

210kn1ixx xk

i1k

ii

==+ , (8.16)

unde eroarea xeste impus de ctre utilizator.

octave#1> function X=jacobi(A,B,P,delta, max)N = length(B);for k=1:max

for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);

endforerr=abs(norm(X'-P));

relerr=err/(norm(X)+eps);P=X';if (err


6/8

se obin componentele vectorului soluie la iteraia (k+1), acestea sunt folosite la calculul celorlaltecomponente.

Algoritm: Sistemul de ecuatii liniare (8.7) se rezolva iterativ cu formula:

= +=

=

++n

1ij

kjij

1i

1j

1kjiji

ii

1ki xaxab

a

1x

)()()(, K,,,;, 210kn1i ==

(8.17)Pornind cu aproximaia iniial2: n1ix

0i ,,

)( = , formula iterativ (8.17), corespunztoare metodeiGauss-Seidel, permite calcularea componentelor vectorului soluie X(k+1), folosind o parte din componentelesoluiei de la iteraia curent (k+1), respectiv cealalt parte a componentelor soluiei de la iteraiaprecedent (k). Calculul iterativ se efectueaz pn la ndeplinirea criteriului de convergen (8.16).

octave#1> function X=gaussseidel(A,B,P,delta, max)N = length(B);for k=1:max

for j=1:Nif j==1

X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==N

X(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);

else X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)'-A(j,j+1:N)*P(j+1:N))/A(j,j);endif

endforerr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';

if (err


7/8

o pseudosolutie X* in sensul celor mai mici patrate (CMMP) consideram reziduul r=B-A*X. Se cauta solutiaX* astfel incat ||B-A*X*||2=min||B-A*X||2 pentru xRnPseudosolutia X* poate fi scrisa sub forma X*=A+*B , unde matricea A+Rn x m , definita prin A+=(AT*A)-1*ATsi se numeste pseudoinversa lui A(Vezi Exemplul 3 din Fisa_4) SI ST E ME SU B DE T E R MI NA TEConsideram un sistem de n ecuatii liniare cu m necunoscute, A*X=B, unde ARn x m si BRn .Presupunem can < m si ca toate liniile sunt independente. Un astfel de sistem este compatibil nedeterminat (in sensulteoremei de compatibilitate a lui Kroneker-Capelli). Problema este abordata in sensul unei solutii normale(de norma minima). Se cauta o solutie X* in sensul celor mai mici patrate (CMMP) cu proprietatea||X*||2=min||X||2 pentru xRmSolutia X* poate fi scrisa sub forma X*=A+*B , unde matricea A+Rn x m , definita prin A+= AT *(A* AT)-1 si senumeste pseudoinversa normala a lui AExemplul 5. Sa se rezolve sistemul liniar subdeterminat4*x1+x2-x3=2x1+5*x2+4*x3=1Solutie teoretica: X1=(9/19)*(1+X3), X2=(2-17*X3)/19, X3Solutie numerica in sensul CMMP:octave#1>A=[4 1 1;1 5 4];B=[2 1];x=A*inv(A*A)*Boctave#2>X1=A\BReturneaza X=X1 =0.21751 0.25581 0.12585

5.4. CONDITIONAREA SISTEMELOR LINIARE

Algoritmii de calcul folositi pentru rezolvarea sistemelor liniare de tip A*X = B, trebuie sa fie numericstabili ,in sensul ca produc o solutie calculata X* care coincide cu solutia exacta a unui sistem de aceiasiforma dar cu datele foarte putin perturbate (A+E)*X* = B unde E este perturbatia.O masura a preciziei solutiei obtinute este data de conditionarea numericaa sistemului evaluata prinnumarul de conditionare la inversare, cond(A)=||A||*||A-1||1, al matricii A. Se apreciaza ca la numerede conditionare mici (apropriate de 1) sistemul este bine conditionat; pentru numere de conditionare marisistemul este prost conditionat (Vezi Exemplul 2 din Fisa_1)Exemplul 5. Fie sistemul liniarx1+1.001*x2=2.001x1+x2=2Solutia exacta a acestui sistem este: x1=1, x2=1. Se observa ca solutia x1=2, x2=0 verifica sistemul cu oaproximatie de ordinul 10-3 . Cauza aparitiei acestei anomalii este numarul de conditionare al matriciiatasata sistemului si anume cond(A)= 4002.0

Functii Octave/Matlab pentru rezolvarea sistemelor liniareOperatorul \ (backslesh) inv(A), inverse(A), ginv(A), det(A) pentru matrici patratice nesingulare rank(A), cond(A), pinv(A) calculul pseudoinversei lu(A),qr(A), chol(A), hess(A), schur(A), svd(A), qzhess(A), qz(A)- factorizari syl(A,B,C) - rezolva ecuatia matriciala a lui Sylvester : A X + X B + C = 0 are, dare rezolva ecuatii algebrice de tip Ricatti

A P L I C AII

Problema 1. Rezolvati sistemul de ecuatii liniare de mai jos, prin metoda de eliminare Gauss.

=

=+

=++

8xx8

5x4x2x

2xx5x

21

321

321

Problema 2. Gasiti solutiile sistemului de la Problema 1, cu metoda Jacobi, dupa 3 iteratii, cu aproximatia

initiala . Comparati cu rezultatele obtinute cu metoda de eliminare Gauss. Sa se

afle de asemenea solutia sistemului folosind functia Octave jacobi cu delta = 1e-6 si maximum 50 deiteratii

( ) ( ) ( )0xxx

03

02

01 ===

Problema 3. Gasiti solutiile sistemului de la Problema 1, cu metoda Gauss-Seidel, dupa 3 iteratii, cu

aproximatia initiala . Comparati cu rezultatele obtinute cu metoda Jacobi si cu( ) ( ) ( )

0xxx03

02

01 ===


8/8

metoda de eliminare Gauss. Sa se afle de asemenea solutia sistemului folosind functia Octave gaussseidelcu delta = 1e-6 si maximum 50 de iteratii

Problema 4. Gasiti descompunerea LU a matricei A, unde si apoi solutia sistemului

=

131

112

511

A

BXA = , folosind factorizarea triunghiulara, daca

=

5

2

7

B

Rezolvare: Pentru A se cauta o descompunere in doua matrici L si U astfel incat:L=[l11 0 0; l21 l22 0; l31 l32 l33] , U=[1 u12 u13; 0 1 u23; 0 0 1]. Din identitatea L*U=A se obtine unsisteme de 9 ecuatii cu 9 necunoscute.Odata determinate matricile L si U se rezolva sistemul astfel: A*X=B devine L*U*X=B. Rezolvam pe randsistemele liniare U*X=Y si U*Y=BProblema 5. Pentru sistemul de la Problema 1 sa se calculeze numarul de conditionare al matriciicoeficientilor necunoscutelor folosind functia Octave cond. Comentarii.Aceiasi problema pentru matricea A1=[ 1 2 3;4 5 6;7 8 9]. Indicatie: Pentru matricea sistemului de laproblema 1 se obtine cond(A)= 2.1389 iar pentru matricea A1, cond(A1)=3.8131e+16

Problema 6. Folosind procedura Gauss si regula pivotului sa se determine rangul matriciiA=[1 6 11;2 7 12;3 8 13;4 9 14;5 10 15]

1 6 112 7 123 8 134 9 145 10 15

1 6 110 -5 -100 -10 -200 -13 -300 -20 -40

1 6 110 -5 -100 0 00 0 00 0 0

CALCUL MATRICIAL. Metoda de eliminare Gauss, Factorizarea LU, Metoda lui Jacobi, Metoda...

Documents

Transcript of CALCUL MATRICIAL. Metoda de eliminare Gauss, Factorizarea LU, Metoda lui Jacobi, Metoda...