Curs_04

Capitolul 2 SISTEME DETERMINATE DE ECUAŢII

ALGEBRICE LINIARE

2.1 Formularea problemei

Se consideră sistemul de n ecuaţii algebrice liniare cu n necunoscute:

1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅ . (2.1)

Problema de calcul este determinarea unei soluţii 1nx ×ℜ∈ a sistemului de ecuaţii (2.1) sau, altfel spus, determinarea unui vector care să satisfacă ecuaţia dată.

Definiţie:

Oricare ar fi matricea nnA ×ℜ∈ , aceasta se numeşte inversabilă (nesingulară) dacă există o matrice nnX ×ℜ∈ astfel încât să fie îndeplinită relaţia: nIAXXA =⋅=⋅ .

În relaţia anterioară, nI este notaţia pentru matricea unitate de ordinul n. În continuare, se notează cu 1A − matricea X care satisface relaţia din definiţia anterioară.

Teoremă de existenţă şi unicitate:

Dacă matricea nnA ×ℜ∈ este inversabilă, atunci oricare ar fi vectorul 1nb ×ℜ∈ există şi este unică soluţia x a ecuaţiei (2.1).

Soluţia se poate scrie sub forma:

bAx 1 ⋅= − . (2.2)

Observaţii:

1. 1A − este notaţia pentru inversa matricei A. În practică nu se recomandă calculul matricei inverse şi apoi aplicarea relaţiei (2.2). Un exemplu care să ilustreze aceasta este următorul: se consideră ecuaţia bxa =⋅ , în care 7a = şi 21b = . Se mai consideră o aritmetică a virgulei mobile cu 10=β , 6t = şi

36 2. Sisteme determinate de ecuaţii algebrice liniare

rotunjire prin tăiere. Atunci soluţia ecuaţiei este:

99997.221142857.021)7(bax 11 =⋅=⋅=⋅= −− .

Acelaşi sistem se poate rezolva direct şi anume: 37/21x == , aceasta reprezentând soluţia exactă.

2. Nu se recomandă rezolvarea ecuaţiei prin regula Cramer: n,...,1i,/x ii =ΔΔ= ,

unde Δ reprezintă determinantul matricei A, iar iΔ reprezintă determinanţii matricelor obţinute prin înlocuirea coloanei numărul i a matricei A cu termenul liber asociat ecuaţiei (2.1). Un exemplu pentru această situaţie este acela în care considerând 20n = , rezultă că trebuie estimaţi 21 de determinanţi care, dacă sunt calculaţi după definiţie, necesită calculul a !20 termeni care implică 19 înmulţiri/termen, deci !2019× operaţii în virgulă mobilă. O operaţie în virgulă mobilă înseamnă o înmulţire şi o adunare: zwvq ⊗⊕= . Pentru un calculator cu 100000 înmulţiri/secundă, rezultă că numai pentru efectuarea înmulţirilor sunt necesari 8103 ⋅ ani, la care se adaugă erorile de rotunjire!

Exemplul 2.1:

Problemele de tipul (2.1) sunt foarte des întâlnite în practică. Ca exemplu generic, se consideră un proces dinamic cu n mărimi de intrare ( 1u , ..., nu ) şi o mărime de ieşire (y). Corespunzător funcţionării acestuia în regim staţionar sau în regim dinamic, pentru care se realizează o liniarizare a modelului neliniar după un punct de funcţionare, modelul intrare-ieşire al procesului poate fi descris de o ecuaţie de forma:

nn11n1 ucuc)u,,u(y ⋅++⋅= …… . (2.3)

Se consideră un număr de n experimente prin care se impun anumite valori mărimilor de intrare şi se măsoară valorile pe care le ia mărimea de ieşire a procesului. Rezultatele obţinute pot fi tabelate astfel:

nr. experiment u1 un y

1 a11 a1n b1

i ai1 ain bi

2.1 Formularea problemei 37

n an1 ann bn

Notând cu A matricea formată din valorile pe care le iau mărimile de intrare ale procesului,

nj1ni1ij ]a[A

≤≤≤≤= şi cu ni1i ]b[b ≤≤= vectorul valorilor mărimii de

ieşire, relaţia (2.3) se poate rescrie sub forma 1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅ . (2.1(2.1), unde T

n21 ]ccc[x …= , T fiind notaţia pentru operaţia de transpunere vectorială. Rezolvarea ecuaţiei rezultate constituie, în acest caz, o identificare a modelului pentru regimul static de funcţionare a procesului. Problema aceasta va avea o soluţie, dacă şi numai dacă liniile sau coloanele matricei experimentelor [ ]b|A sunt vectori liniari independenţi, după cum se va specifica ulterior.

Metodele numerice de rezolvare a unui sistem de ecuaţii algebrice liniare se împart în următoarele două categorii: metode directe şi metode iterative (indirecte).

Metodele directe se bazează pe reducerea sistemului (2.1) la un sistem echivalent, direct rezolvabil prin mijloace elementare. Principial, aceasta foloseşte eliminarea progresivă a necunoscutelor, numită şi eliminare gaussiană.

Practic, prin transformări elementare de echivalenţă, se aduce matricea A a sistemului la anumite forme tipice: a) forma superior triunghiulară:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

nn

n111

u00

0uu

U

…

……

;

b) forma inferior triunghiulară:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

nn1n

11

ll0

00l

L

……

…

.

Procedura de transformare se numeşte triangularizare. Astfel, se furnizează soluţia exactă a sistemului de ecuaţii (2.1), în cazurile (ideale) în care erorile de rotunjire sunt absente. Numărul de operaţii în virgulă mobilă, necesare triangularizării unei matrice pătratice de ordin n, este de ordinul lui 3n . De aceea, aceste metode se recomandă pentru rezolvarea unor sisteme de ecuaţii algebrice liniare de ordin mai mic decât 210 .


Metodele iterative au la bază construirea unui şir de aproximaţii pentru soluţia sistemului (2.1), …,1,0k,x ]k[ = care să fie convergent pentru ∞→k la soluţia adevărată, x :

0||xx||lim ]k[

k=− α∞→

.

Practic, calculele se opresc la un index de iterare [s], atunci când este îndeplinită o condiţie de forma:

impus]1s[]s[ ||xx|| ε≤− α

− ,

sau, altfel spus, ]s[x constituie o aproximare satisfăcătoare a soluţiei calculate. Având în vedere faptul că pentru o singură iteraţie numărul de operaţii în

virgulă mobilă este de ordinul lui 2n , asemenea metode se folosesc pentru sisteme de ordin mare şi anume 42 1010 ÷ .

2.2 Rezolvarea sistemelor prin triangularizare directă

2.2.1 Principiul metodei

Se consideră problema (2.1) şi matricea sistemului: nj1ni1ij ]a[A

≤≤≤≤= .

Definiţie:

Matricele ]a[ 11 , …, 1nj11ni1ij ]a[−≤≤−≤≤ se numesc submatrice principale ale lui A

sau minori principali directori.

Teoremă:

Dacă matricea nnA ×ℜ∈ are toate submatricele principale inversabile (nesingulare), atunci există matricele nnU,D,L ×ℜ∈ astfel încât:

UDLA ⋅⋅= (2.4)

unde L este o matrice inferior triunghiulară, D este o matrice diagonală şi U este o matrice superior triunghiulară.

Se pot face, în cele ce urmează, următoarele observaţii: 1. relaţia (2.4) se numeşte factorizare L-D-U a matricei A; 2. uzuale sunt factorizările: ULA ⋅= ; 3. demonstraţia teoremei enunţate este constructivă, constituind însuşi

algoritmul de descompunere L-U a matricei A;

2.2 Rezolvarea sistemelor prin triangularizare directă 39

4. algoritmul de descompunere este în fond procedeul de eliminare gaussiană, prin care matricea A este adusă la forma superior triunghiulară în urma unui şir de transformări de asemănare. Transformările efectuate asupra matricei A se “acumulează” într-o matrice inferior triunghiulară, cu elementele de pe diagonala principală egale cu 1. Acest tip de descompunere se numeşte descompunere Doolittle.

5. considerând descompunerea L-U a matricei sistemului A, ULA ⋅= , atunci rezolvarea sistemului (2.1) implică două subetape: a. rezolvarea sistemului byL =⋅ , etapă numită şi substituţie înainte,

obţinând soluţia intermediară y . Determinarea componentelor vectorului ni1i ]y[y ≤≤= are loc din aproape în aproape: se începe cu 1y

(prima ecuaţie), se înlocuieşte în a doua ecuaţie determinând pe 2y şi aşa mai departe.

b. rezolvarea sistemului yxU =⋅ , în care necunoscuta este x , etapă

numită şi substituţie inversă. În acest caz, determinarea componentelor vectorului x are loc pornind de la ultima ecuaţie.

Această manieră de descompunere şi de rezolvare se încadrează în aşa-numita rezolvare a sistemelor determinate de ecuaţii algebrice liniare prin triangularizare simplă (directă).

În continuare se prezintă câteva rezultate referitoare la algebra matricelor triunghiulare.

Definiţie:

Se numeşte matrice triunghiulară unitate o matrice inferior sau superior triunghiulară care are elementele de pe diagonala principală egale cu 1.

În ceea ce priveşte matricele triunghiulare, sunt valabile următoarele rezultate. R1: Inversa unei matrice superior (inferior) triunghiulară este o matrice superior

(inferior) triunghiulară. R2: Produsul a două matrice superior (inferior) triunghiulare este o matrice

superior (inferior) triunghiulară. R3: Inversa unei matrice superior (inferior) triunghiulară unitate este o matrice

superior (inferior) triunghiulară unitate. R4: Produsul a două matrice superior (inferior) triunghiulare unitate este o

matrice superior (inferior) triunghiulară unitate. În cele ce urmează, se enunţă şi demonstrează următorul rezultat esenţial.

Propoziţie:


Dacă matricea A admite o descompunere L-U, atunci această descompunere este unică.

Demonstraţia se realizează prin reducere la absurd, presupunând că matricea A admite două descompuneri L-U şi anume: 11 ULA ⋅= , 22 ULA ⋅= . De aici rezultă că 2211 ULUL ⋅=⋅ , ceea ce conduce la:

1121

12 UULL −− ⋅=⋅ ,

adică o matrice inferior triunghiulară unitate ( 11

2 LL ⋅− ) este identică cu o matrice superior triunghiulară ( 1

12 UU −⋅ ). Acest lucru este posibil numai dacă ambele matrice sunt diagonale şi au diagonala principală unitară, adică:

n11

2 ILL ≡⋅− şi n1

12 IUU ≡⋅ − , ceea ce implică 21 LL ≡ , 21 UU ≡ .

Procedura de triangularizare directă necesită un număr de operaţii în virgulă mobilă de ordinul lui 3/n 3 . Numărul total de operaţii în virgulă mobilă pentru rezolvarea unui sistem determinat de ecuaţii algebrice liniare, folosind triangularizarea simplă, este de ordinul lui 23 n)3/n( + , 2n operaţii fiind necesare pentru parcurgerea celor două etape din rezolvarea propriu-zisă a sistemului, anume substituţia înainte şi substituţia înapoi.

Dacă matricea A este simetrică ( TAA = ) şi pozitiv definită ( 0xAx,0x,x T

n1n >⋅⋅≠ℜ∈∀ × şi n

T 0x0xAx ≡⇔=⋅⋅ ), atunci A se

descompune sub forma TLLA ⋅= . Aceasta se numeşte descompunerea Cholesky. În acest caz, algoritmul necesită mai puţine operaţii în virgulă mobilă şi anume 6/n 3 , exploatând faptul că matricea A este simetrică.

2.2.2 Procedura de triangularizare directă a unei matrice

Principiul triangularizării simple poate fi prezentat prin următorul algoritm, descris principial în limbajul pseudocod: atribuie AA1 ← pentru 1n,1k −= execută ⎡ * determinare matrice kM astfel încât matricea kk1k AMA ⋅=+ să ⎢ aibă elementele: ⎢ n,...,1ki,0a ]1k[

k,i +==+ şi 1k,...,1j;n,...,1i,aa ]k[j,i

]1k[j,i −===+

⎢ atribuie kk1k AMA ⋅←+ ⎣


În final se obţine matricea UAn = . Acest algoritm parcurge )1n( − etape, la fiecare etapă zerorizându-se

elementele de sub diagonala principală şi păstrând nealterate transformările care s-au efectuat în coloanele anterioare ale matricei A.

Notând cu ξ vectorul conţinând elementele coloanei k a matricei kA ,

anume: T]k[

k,n]k[

k,1k]k[k,k

]k[k,1

Tn1kk1 ]aaaa[][ ………… ++ =ξξξξ=ξ ,

atunci matricea de transformare kM se construieşte astfel încât vectorul ξ⋅kM

să aibă elementele: T

k1k ]00[M …… ξξ=ξ⋅ .

Se consideră vectorul km de forma: T

k,nk,1kk ]00[m μμ= + …… , elementele n,...,1ki,k,i +=μ numindu-se multiplicatori.

Vectorul km se numeşte vector Gauss sau vector de multiplicatori. Subvectorul care conţine strict numai multiplicatorii se numeşte subvector Gauss:

Tk,nk,1kk ][t μμ= + … .

Definiţie:

Matricea kM se numeşte matrice de transformare elementară de ordin n şi indice k sau matrice Gauss şi este definită prin:

Tkknk emIM ⋅−= ,

în care TTk ]00100[e ……= , elementul egal cu 1 fiind în poziţia k.

Definită astfel, matricea kM este o matrice inferior triunghiulară unitate, este nesingulară şi deci admite inversă. Inversa acesteia este de forma:

Tkkn

1k emIM ⋅+=− . (2.5)

Efectul aplicării matricei kM asupra vectorului definit anterior, ξ , este:

Tkk,nnkk,1k1kk1

kkTkk

Tkknk

][

mem)emI(M

ξ⋅μ−ξξ⋅μ−ξξξ=

ξ⋅−ξ=ξ⋅⋅−ξ=ξ⋅⋅−=ξ⋅

++ …… ultimele kn − elemente trebuind a fi zerorizate. Presupunând că 0k ≠ξ şi alegând n,...,1ki,/ kik,i +=ξξ=μ , va rezulta:

Tk1k ]00[M …… ξξ=ξ⋅ .


Dacă la etapa k a triangularizării, elementul 0a ]k[k,k ≠ şi

n,...,1ki,a/a ]k[k,k

]k[k,ik,i +==μ , atunci se obţine:

T]k[k,k

]k[k,1k

]k[k,1kkk ]00aaa[)A(cM …… −=⋅ ,

în care )A(c kk reprezintă notaţia pentru coloana k a matricei kA . Acest rezultat evidenţiază faptul că primele k elemente din coloana k a matricei kA rămân neschimbate, iar ultimele kn − elemente devin zero. Elementul

]k[k,kk a=ξ se numeşte pivot.

Observaţii:

1. În practică, pe calculator, etapa k descrisă mai sus se poate realiza testând condiţia:

ε>|a| ]k[k,k ,

în loc de a verifica 0a ]k[k,k ≠ , unde ε este o constantă impusă, de valoare

mică sau foarte mică. De exemplu, constanta ε poate fi egală cu epsilonul-maşină. Aceasta se realizează datorită faptului că, dacă în aritmetica reală (exactă) pivotul este nul, în aritmetica virgulei mobile, datorită erorilor de calcul, această situaţie este echivalentă cu:

ε≤|a| ]k[k,k .

2. Când pivotul este în modul mai mic sau egal cu ε , eliminarea gaussiană eşuează. Aceasta corespunde situaţiei când matricea iniţială A are submatricea principală de ordin k singulară, deci conform teormei enunţate anterior, descompunerea L-U a matricei A nu există.

Efectul aplicării transformării kM asupra celorlaltor coloane ale matricei

kA este următorul. Se consideră un vector ξ≠η de forma Tn1 ][ ηη=η … .

Aplicând transformarea kM vectorului η , se obţine: T

kk,nnkk,1k1kk1k ][M η⋅μ−ηη⋅μ−ηηη=η⋅ ++ …… .

Concluzii:

a) Matricea kM lasă nemodificate primele 1k − coloane ale matricei kA . Considerând vectorul η ca fiind coloana numărul j a matricei kA :

1k,...,1j,kj,]000**[)A(c T

nk1jj1kj −=<==η+

……… ,

atunci se obţine:


)A(c

]000000**[)A(cM

kj

T

nk,n

1kk,1kk1jj1kjk

=

⋅μ−⋅μ−=⋅++

+………

b) Matricea kM transformă coloana k a matricei kA , zerorizând liniile n,...,1k + .

c) Matricea kM transformă coloanele n,...,1k + ale lui kA în liniile n,...,1k + . Considerând vectorul η ca fiind coloana j a matricei kA :

n,...,1kj),A(c kj +==η ,

atunci se obţine: T]k[

j,kk,n]k[j,n

]k[j,kk,1k

]k[j,1kkjk ]aaaa**[)A(cM ⋅μ−⋅μ−=⋅ ++ …… ,

notaţia * semnificând faptul că elementele implicate rămân nemodificate.

Sumarizând, asupra matricei iniţiale A se aplică transformările 1M , ...,

1nM − , obţinându-se în final forma superior triunghiulară U:

UAMMM 121n =⋅⋅⋅⋅− … . (2.6)

Înmulţind la stânga relaţia (2.6) cu 11n

12

11 MMM −

−−− ⋅⋅⋅ … , se obţine:

UMMMA 11n

12

11 ⋅⋅⋅⋅= −

−−− … . (2.7)

Notând produsul 11n

12

11 MMM −

−−− ⋅⋅⋅ … cu L, atunci relaţia (2.7) devine:

ULA ⋅= . Ţinând cont de forma matricelor 1n,...,1k,M 1

k −=− dată de relaţia (2.5), atunci matricea L este egală cu:

∑−

=

−−

−− ⋅+=⋅⋅⋅=1n

1k

Tkkn

11n

12

11 emIMMML … .

Matricea L este inferior triunghiulară unitate şi conţine în fiecare coloană, sub elementul unitar de pe diagonala principală, subvectorii Gauss.

Prima sub-etapă de rezolvare a sistemului (2.1) este substituţia înainte aplicată sistemului de ecuaţii byL =⋅ . Vectorul y rezultat este, de fapt,

vectorul care se obţine aplicând la stânga, în aceeaşi ordine, transformările elementare care s-au aplicat matricei A:

bMMMbLy 121n1 ⋅⋅⋅⋅=⋅= −− … .

Exemplul 2.2:

Se consideră o aritmetică a virgulei mobile cu 10=β , 5t = şi rotunjire prin tăiere.


?x,bxA;6901.37

b;5156099.230710

A ==⋅⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

−=

;1015005.0**

6101*0710

A;55.2*6101*0710

A;AA5

33

321

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

⋅⋅−−

=⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡⋅−−

== −−

a5

x11

0

99993.05.135.0

x;1015004.0

001.67

y =⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡−≠

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−

=⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

⋅= ,

unde ax reprezintă soluţia adevărată a sistemului. Rezultă, aşadar, erori mari în soluţia calculată, x . Cauza care a determinat apariţia acestor erori este aceea că la pasul al doilea al triangularizării s-a lucrat cu un pivot foarte mic în modul ( 3101 −⋅− ), pentru aritmetica virgulei mobile folosite. Multiplicatorul corespunzător este 33

2,3 105.2)10/(5.2 ⋅=−=μ − , deci are o valoare foarte mare în modul pentru aceeaşi aritmetică. Aceasta a condus, mai departe, la apariţia fenomenului de omitere catastrofală în calculele care s-au efectuat pentru obţinerea vectorilor y şi x .

Concluzie:

La triangularizarea simplă, unde elementele matricei A se modifică corespunzător relaţiei:

n,...,kj;n,...,1ki,aaa ]k[kjik

]k[ij

]1k[ij =+=⋅μ−=+ ,

multiplicatorii ikμ pot avea, în principiu, orice valoare. Dacă aceste valori sunt mari sau foarte mari, atunci pot apare fenomenele de omitere catastrofală şi/sau de neutralizare a termenilor. Mai mult, dacă aceşti multiplicatori au valori supraunitare în modul, atunci ei amplifică erorile prezente în termenii ]k[

kja , în

felul acesta triangularizarea simplă fiind instabilă numeric, în general. Altfel spus, nu există nici un control asupra stabilităţii numerice a algoritmului triangularizării simple.

2.3 Rezolvarea sistemelor prin triangularizare cu pivotare parţială

În cazul triangularizării cu pivotare parţială (Figura 2.1), la pasul k se caută pivotul kξ printre elementele din coloana k, pornind de la elementul de pe

2.3 Rezolvarea sistemelor prin triangularizare cu pivotare parţială 45

diagonala principală în jos, alegându-se elementul care are cea mai mare valoare în modul:

k]k[

k,inik

]k[k,i |}a{|max|a|

kξ==

≤≤.

Ak =

|k

⎯ k

⎯ ik

⎯ n

0

Fig. 2.1 Principiul triangularizării cu pivotare parţială a unei matrice: pivotul se

găseşte în coloana k, liniile k÷ n; k=1,...,n

Dacă ki k ≠ , elementul maxim în modul nu se găseşte pe diagonala principală, atunci se interschimbă (permută) liniile k şi ki . Lucrul acesta se realizează automat cu ajutorul unei matrice de permutare de linii kP care multiplică la stânga matricea kA : kk AP ⋅ . Pentru matricea care rezultă astfel, se determină apoi matricea de transformare kM ca şi în cazul traingularizării simple, obţinând matricea:

)AP(MA kkk1k ⋅⋅=+ . În felul acesta, multiplicatorii calculaţi sunt subunitari în modul

n,...,1ki,1|| k,i +=≤μ , iar algoritmul triangularizării devine stabil numeric. Matricea kk PM ⋅ se numeşte matrice de transformare elementară stabilizată.

Matricea de permutare kP se obţine din matricea unitate de ordinul n, nI , schimbând adecvat între ele liniile k şi ki . Această matrice are următoarele proprietăţi:

1kkk PP;1)Pdet( −=−= .

Teoremă:

Dacă matricea nnA ×ℜ∈ este nesingulară, atunci există o matrice nnP ×ℜ∈ , numită matrice generală de permutare de linii, astfel încât:

ULAP ' ⋅=⋅ ,


în care U este o matrice superior triunghiulară şi 'L este o matrice inferior triunghiulară unitate. cu elementele ji,1|l| j,i >≤ .

Demonstraţia acestei teoreme este constructivă, constituind însuşi algoritmul de triangularizare cu pivotare parţială a unei matrice. Se poate face observaţia că matricea 'L conţine în coloana k, sub elementul de pe diagonala principală, subvectorul Gauss (de la triangulaizarea simplă) având liniile permutate.

Algoritmul de triangularizare cu pivotare parţială este următorul: atribuie AA1 ← pentru 1n,1k −= execută ⎡ * determinare pivot ]k[

k,ik ka←ξ astfel încât |}a{|max|a| ]k[

k,inik

]k[k,ik ≤≤=

⎢ dacă )ki( k ≠ atunci ⎢ ⎡ * determină kP ⎢ ⎢altfel ⎢ ⎢ atribuie nk IP ← ⎢ ⎣ ⎢ * calcul kk AP ⋅ ⎢ * determinare matrice kM astfel încât matricea ⎢ )AP(MA kkk1k ⋅⋅=+ să îndeplinească condiţiile de la ⎢ triangularizarea simplă ⎢ atribuie )AP(MA kkk1k ⋅⋅←+ ⎣ atribuie nAU ←

În ansamblu, asupra matricei A sunt aplicate următoarele transformări:

UAPMPMPM 11221n1n =⋅⋅⋅⋅⋅⋅⋅ −− … . (2.8)

Folosind faptul că 1kk PP −= , sau altfel spus nkk IPP =⋅ , relaţia (2.8) poate fi

scrisă sub forma:

UA)PPP()PPPPMPM( 121n1n21111n1n =⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ −−−− ……… (2.9)

Matricea produs 121n PPP ⋅⋅⋅− … se notează cu P. Ea este numită matrice generală de permutare de linii. Produsul din prima paranteză din relaţia (2.9) se notează cu 1'L− . Relaţia (2.9) se scrie sub forma: ULAP ' ⋅=⋅ , în care matricea

11n1n

122

1121n

' MPMPMPPL −−−

−−− ⋅⋅⋅⋅⋅⋅⋅⋅= …… este o matrice inferior

triunghiulară unitate având în fiecare coloană, sub diagonala principală, subvectori Gauss cu liniile permutate.

2.3 Rezolvarea sistemelor prin triangularizare cu pivotare parţială 47

Definiţie:

Matricea nnA ×ℜ∈ se numeşte diagonal dominantă pe coloane dacă în fiecare coloană a sa elementul de pe diagonala principală este, în modul, mai mare decât suma modulelor celorlaltor elemente:

n,...,1j,|a||a|n

ji1i

j,ij,j =∀≥∑≠=

.

În acest caz, se poate enunţa următorul rezultat.

Propoziţie:

Dacă matricea A este diagonal dominantă pe coloane, atunci ea admite factorizarea:

ULAP ' ⋅=⋅ ,

în care nIP = şi elementele matricei 'L sunt în modul subunitare.

Altfel spus, nu este implicată permutarea de linii în cazul traingularizării unei matrice diagonal dominante pe coloane.

Rezolvarea sistemului 1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅ se realizează în două etape: a) descompunerea 'L -U a matricei A:

ULAP ' ⋅=⋅ b) rezolvarea propriu-zisă a sistemului care comportă trei subetape, şi anume:

b.1.) calculul vectorului bPc ⋅= ; b.2.) rezolvarea sistemului cyL' =⋅ prin substituţie înainte; b.3.) rezolvarea sistemului yxU =⋅ prin substituţie înapoi.

Această modalitate de rezolvare se bazează pe următoarele relaţii: bPxAP ⋅=⋅⋅ , ULAP ' ⋅=⋅ , bPxUL' ⋅=⋅⋅ , xUy ⋅= , bPc ⋅= .

Exemplul 2.3:

Se consideră problema de la Exemplul 2.2. Se aplică, de această dată, triangularizarea cu pivotare parţială. În acest caz, permutarea de linii intervine la pasul al doilea al triangularizării, permutându-se liniile 2 şi 3: 2PP = . Se obţine, în final:

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡−=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡ −=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

11

0x;

002.60055.200710

U;002.6

5.27

y .


Dacă algoritmul de triangularizare cu pivotare parţială eşuează, în sensul că pivotul găsit la o anumită etapă [k] este nul sau foarte mic în modul, aceasta corespunde situaţiei când în aritmetica reală primele k coloane ale matricei A sunt liniar dependente. Dacă însă pivotul găsit este în modul foarte mic în sensul preciziei dorite (mai mic decât un anumit parametru ε impus), atunci se recurge la traingularizarea matricei sistemului prin pivotare totală, ceea ce implică atât permutări de linii cât şi de coloane.

2.4 Rezolvarea sistemelor prin triangularizare cu pivotare totală

În cazul triangularizării cu pivotare totală (Figura 2.2), la pasul [k] al triangularizării se alege drept pivot elementul maxim în modul din submatricea formată din liniile de la k la n, coloanele de la k la n:

|}a{|max|a| ]k[j,i

njknik

]k[j,i kk

≤≤≤≤

= .

Dacă acest element nu se află în linia şi/sau coloana k, atunci are loc permutarea adecvată de linii şi/sau coloane în scopul aducerii acelui element pe diagonala principală, anume în linia k şi coloana k.

Ak =

|k

⎯ k

⎯ ik

⎯ n

0

|jk

Fig. 2.2 Principiul triangularizării cu pivotare totală a unei matrice: pivotul se găseşte în submatricea determinată de coloanele k ÷ n şi liniile k÷ n; k=1,...,n

În continuare, se enunţă şi demonstrează următorul rezultat.

Teoremă:

Pentru orice matrice nnA ×ℜ∈ nesingulară, există două matrice generale de permutare, P – matrice generală de permutare de linii şi S – matrice generală de permutare de coloane, astfel încât:

ULSAP ' ⋅=⋅⋅ ,

2.4 Rezolvarea sistemelor prin triangularizare cu pivotare totală 49

unde U este o matrice superior triunghiulară, iar 'L este o matrice inferior triunghiulară unitate având elementele ji,1|l| j,i ≥≤ , în fiecare coloană a

matricei L’, sub elementul de pe diagonala principală, găsindu-se subvectori Gauss având liniile permutate între ele. Matricele generale de permutare P şi S sunt:

1n21121n SSSS;PPPP −− ⋅⋅⋅=⋅⋅⋅= …… , unde matricele 1n,...,1k,S,P kk −= sunt matrice de permutare de linii şi, respectiv, de coloane.

Demonstraţia teoremei este constructivă, reprezentând însuşi algoritmul triangularizării cu pivotare totală a matricei A. Acesta este descris în limbajul pseudocod, după cum urmează:

atribuie AA1 ← pentru 1n,1k −= execută ⎡ * determinare pivot ]k[

j,ik kka←ξ care satisface: |}a{|max|a| ]k[

j,injknik

]k[j,i kk

≤≤≤≤

=

⎢ dacă )ki( k ≠ atunci ⎢ ⎡ * determinare kP (permutarea liniilor ki şi k) ⎢ ⎢altfel ⎢ ⎢ atribuie nk IP ← ⎢ ⎣ ⎢ dacă )kj( k ≠ atunci ⎢ ⎡ * determinare kS (permutarea coloanelor kj şi k) ⎢ ⎢altfel ⎢ ⎢ atribuie nk IS ← ⎢ ⎣

⎢ atribuie kkk'

1k SAPA ⋅⋅←+ ⎢ * traingularizare matrice '

1kA + : ⎢ )SAP(MAMA kkkk

'1kk1k ⋅⋅⋅=⋅← ++

⎣ atribuie nAU ←

Tabloul general al transformărilor este:

USSSAPMPMPM 1n2111221n1n =⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ −−− …… . (2.10)

În relaţia (2.10) se notează cu S produsul 1n21 SSS −⋅⋅⋅ … şi se ţine cont de faptul că nkk IPP =⋅ . Atunci relaţia (2.10) devine:


USA)PP(

)PPPPMPMPM(

11n

1n2111221n1n

=⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

−

−−−

………

(2.11)

Matricea 121n PPP ⋅⋅⋅− … se notează cu P, iar cu 1' )L( − se notează produsul

1n2111221n1n PPPPMPMPM −−− ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ …… . În felul acesta, relaţia (2.11) devine: ULSAP ' ⋅=⋅⋅ .

Observaţii:

1. Matricea de permutare de linii, kP , se obţine din matricea unitate nI schimbând adecvat între ele liniile k şi ki . Deoarece se schimbă linii între ele, matricea kP se aplică la stânga matricei kA . Matricea kP are proprietăţile:

1kkk PP;1)Pdet( −=−= .

2. Matricea de permutare de coloane, kS , se obţine din matricea unitate nI schimbând adecvat între ele coloanele k şi kj . Deoarece se schimbă coloane între ele, matricea kS se aplică la dreapta matricei kA . Matricea kS are proprietăţile:

1kkk SS;1)Sdet( −=−= .

Pentru înţelegerea etapelor rezolvării sistemului (2.1), se înmulţesc ambii membri ai ecuaţiei, la stânga, cu matricea P şi se inserează între A şi x produsul n

1 ISS =⋅ − , obţinându-se:

bPxSSAP 1 ⋅=⋅⋅⋅⋅ − . (2.12)

Aşadar, rezolvarea sistemului comportă următoarele etape: a) triangularizarea cu pivotare totală a matricei siatemului A:

ULSAP ' ⋅=⋅⋅ ; b) rezolvarea propriu-zisă a sistemului (2.12), cu următoarele subetape:

b.1.) calculul vectorului bPc ⋅= ; b.2.) rezolvarea sistemului cyL' =⋅ prin substituţie înainte;

b.3.) rezolvarea sistemului yzU =⋅ prin substituţie înapoi; b.4.) determinarea soluţiei: zSx ⋅= .

Observaţii:

1) Permutările de linii efectuate asupra matricei A implică permutări de linii asupra termenului liber b .

2.5 Aplicaţii ale descompunerilor L-U 51

2) Permutările de coloane efectuate asupra matricei A implică permutări de linii în soluţia calculată a sistemului z .

3) Triangularizarea cu pivotare totală asigură, la fiecare iteraţie a sa, pivoţii cei mai mari în valoare absolută. Astfel, multiplicatorii vor fi subunitari în modul, de valoarea cea mai mică posibil, 1|| k,i ≤μ , iar elementele care se transformă devin:

n,...,1ki;n,...,kj,aaa ]k[j,kk,i

]k[j,i

]1k[j,i +==⋅μ−←+ .

Ca urmare, triangularizarea cu pivotare totală reprezintă procedura de triangularizare cea mai precisă şi stabilă numeric. Dezavantajul ei este acela că necesită un timp de calcul mai mare. De regulă, se foloseşte triangularizarea cu pivotare parţială, recurgându-se la triangularizarea cu pivotare totală numai când cea parţială eşuează.

4) Dacă A este o matrice singulară, atunci pivotarea totală va eşua. În aritmetica reală exactă aceasta corespunde situaţiei când pivotul este nul, matricea A având rangul egal cu k-1, daca algoritmul eşuează la iteraţia k. În aritmetica în virgulă mobilă, datorită erorilor de rotunjire, un pivot nul înseamnă îndeplinirea condiţiei:

ε≤|a| ]k[j,i kk

,

şi se spune că matricele A, ca şi kA , sunt algoritmic singulare. 5) Dacă matricea A este diagonal dominantă pe linii şi pe coloane şi în plus

elementele de pe diagonală satisfac relaţiile: |a||a||a| nn2211 ≥≥≥ … ,

atunci descompunerea 'L -U cu pivotare totală este:

nn' IS;IP:unde;ULSAP ==⋅=⋅⋅ .

2.5 Aplicaţii ale descompunerilor L-U

2.5.1 Calculul determinantului

Considerând o matrice nnA ×ℜ∈ nesingulară, pentru care s-a calculat una din descompunerile de tip L-U, calculul determinantului acesteia poate fi făcut după cum urmează.

Descompunerea L-U bazată pe triangularizarea simplă

În acest caz, descompunerea L-U a matricei A este:

ULA ⋅= , (2.13)


unde L este o matrice inferior triunghiulară unitate, iar matricea nj,i1j,i ]u[U ≤≤=

este o matrice superior triunghiulară. Aplicând funcţia determinant det(.) relaţiei (2.13) se obţine:

)Udet()Ldet()ULdet()Adet( ⋅=⋅= . (2.14)

Matricea L fiind o matrice inferior triunghiulară unitate, determinantul său este egal cu 1, iar matricea U fiind o matrice superior triunghiulare, determinantul său este egal cu produsul elementelor de pe diagonala principală. Ţinând cont de acestea, relaţia (2.14) devine:

∏=

=n

1ii,iu)Adet( .

Descompunerea L-U bazată pe triangularizarea cu pivotare parţială


ULAP ' ⋅=⋅ , (2.15)

în care P este matricea generală de permutare de linii (matrice nesingulară), 'L este o matrice inferior triunghiulară unitate, iar nj,i1j,i ]u[U ≤≤= este o matrice

superior triunghiulară. Ţinând cont de faptul că matricea P este inversabilă, relaţia (2.15) poate fi scrisă sub forma:

ULPA '1 ⋅⋅= − . (2.16)

Aplicând funcţia determinant det(.) relaţiei (2.16), se obţine: )Udet()Ldet()Pdet()ULPdet()Adet( '1'1 ⋅⋅=⋅⋅= −− .

Ţinând cont de următoarele: • 1n21

1 PPPP −− ⋅⋅⋅= … , )Pdet()Pdet()Pdet()Pdet( 1n21

1−

− ⋅⋅⋅= … ; • }1n,...,1{i,1)Pdet( i −∈−= , deoarece dacă ni IP ≠ este realizată o

permutare de linii la iteraţia [i] a algoritmului de triangularizare cu pivotare parţială;

• 1)Ldet( ' = ;

• ∏=

=n

1ii,iu)Udet( ,

atunci relaţia (2.16) devine:

∏=

⋅−=n

1ii,i

Npl u)1()Adet( ,

unde Npl reprezintă numărul de permutări de linii efectiv realizate în procesul de triangularizare a matricei A.


Descompunerea L-U bazată pe triangularizarea cu pivotare totală


ULSAP ' ⋅=⋅⋅ , (2.17)

în care P este matricea generală de permutare de linii (matrice nesingulară), S este matricea generală de permutare de coloane (matrice nesingulară), 'L este o matrice inferior triunghiulară unitate, iar nj,i1j,i ]u[U ≤≤= este o matrice superior

triunghiulară. Ţinând cont de faptul că matricele P şi S sunt inversabile, relaţia (2.17) devine:

1'1 SULPA −− ⋅⋅⋅= . (2.18)

Dar 1n211 PPPP −− ⋅⋅⋅= … şi 121n

1 SSSS ⋅⋅⋅= −− … . Ca urmare, relaţia (2.18)

devine:

∏∏==

⋅⋅⋅=n

1ii

'n

1ii )Sdet()Udet()Ldet()Pdet()Adet( . (2.19)

Ţinând cont de următoarele: • }1n,...,1{i,1)Pdet( i −∈−= , deoarece dacă ni IP ≠ este realizată o

permutare de linii la iteraţia [i] a algoritmului de triangularizare cu pivotare totală;

• }1n,...,1{i,1)Sdet( i −∈−= , deoarece dacă ni IS ≠ este realizată permutare de coloane la iteraţia [i] a algoritmului de triangularizare cu pivotare totală;

• 1)Ldet( ' = ;

• ∏=

=n

1ii,iu)Udet( ,

atunci relaţia (2.19) devine:

∏=

+ ⋅−=n

1ii,i

NpcNpl u)1()Adet( ,

unde Npl este numărul permutărilor de linii, iar Npc numărul permutărilor de coloane ce s-au realizat efectiv pe parcursul procedurii de triangularizare cu pivotare totală a matricei A.

2.5.2 Rezolvarea ecuaţiilor matriciale

Se consideră ecuaţia matricială de forma:

?X;X;B;A,BXA pnpnnn =ℜ∈ℜ∈ℜ∈=⋅ ××× . (2.20)


Se scriu matricele B şi X pe coloane, sub forma:

p,...,1k,x];xxx[X

p,...,1k,b];bbb[B1n

kpk1

1nkpk1

=ℜ∈=

=ℜ∈=×

×

……

……. (2.21)

Apelând la descompunerea L-U cu traingularizare simplă (relaţia (2.13)) a matricei A a sistemului (2.20), se obţine relaţia:

BXUL =⋅⋅ . (2.22)

Notând produsul XU ⋅ cu Y, relaţia (2.22) devine: BYL =⋅ sau, folosind relaţiile (2.21), se poate scrie:

p,...,1k,yxU

p,...,1k,byL

kk

kk

==⋅

==⋅. (2.23)

Aşadar, rezolvarea ecuaţiei matriciale (2.20) implică descompunerea L-U a matricei sistemului şi apoi determinarea succesivă a coloanelor matricei necunoscutelor, X cu ajutorul relaţiilor (2.23). Numărul de operaţii în virgulă mobilă necesar rezolvării unui astfel de sistem este de ordinul lui

23 np)3/n( ⋅+ , unde 3/n3 operaţii în virgulă mobilă sunt necesare descompunerii L-U, iar 2np ⋅ operaţii sunt necesare rezolvării sistemului în maniera descrisă.

În cazul în care se utilizează descompunerea L-U cu pivotare parţială (relaţia (2.15)) a matricei A a sistemului, atunci rezolvarea sistemului de ecuaţii matriciale (2.20) parcurge următoarele etape: a) descompunerea 'L -U a matricei A (relaţia (2.15)); b) calculul coloanelor matricei X; pentru p,...,1k = se execută:

b.1.) calculul vectorului kbPc ⋅= ;

b.2.) rezolvarea sistemului cyLk

' =⋅ prin substituţie înainte;

b.3.) rezolvarea sistemului kk yxU =⋅ prin substituţie inversă.

În cazul în care se apelează la descompunerea L-U cu pivotare totală (relaţia (2.17)) a matricei A a sistemului, atunci rezolvarea sistemului de ecuaţii matriciale (2.20) parcurge următoarele etape: a) descompunerea 'L -U a matricei A (relaţia (2.17)); b) calculul coloanelor matricei X; pentru p,...,1k = se execută:

b.1.) calculul vectorului kbPc ⋅= ;

b.2.) rezolvarea sistemului cyLk

' =⋅ prin substituţie înainte;


b.3.) rezolvarea sistemului kk yzU =⋅ prin substituţie inversă;

b.4.) ordonare soluţie calculată: kk zSx ⋅= .

2.5.3 Calculul inversei unei matrice

Fie o matrice nnA ×ℜ∈ nesingulară. Se doreşte aflarea inversei acesteia, notată cu 1A − . Acest tip de problemă se încadrează în problematica rezolvării ecuaţiilor matriciale de tipul (2.20), considerând nIB = .

Astfel, în prima fază se utilizează una din descompunerile L-U ale matricei A, anume traingularizare simplă, triangularizare cu pivotare parţială sau triangularizare cu pivotare totală, urmată de o a doua fază de rezolvare propriu-zisă a unui sistem de tipul (2.20). În final, matricea inversă 1A− este egală cu matricea ale cărei coloane sunt vectorii rezultaţi la faza a doua menţionată şi anume: [ ]nk1

1 xxxXA ……==− . Detaliile acestei proceduri sunt următoarele, în funcţie de tipul de triangularizare a matricei de inversat A caee este folosit:

triangularizare simplă

1. descompunere L-U (relaţia (2.13)); 2. kk

eyL =⋅ , Tk ]00100[e ……= , unde 1 apare în poziţia k;

kk yxU =⋅ , n,...,1k =

triangularizare cu pivotare parţială

1. descompunere L-U (relaţia (2.15)); 2. kePc ⋅= , T

k ]00100[e ……= , unde 1 apare în poziţia k;

cyLk

' =⋅

kk yxU =⋅ , n,...,1k =

triangularizare cu pivotare totală

1. descompunere L-U (relaţia (2.17)); 2. kePc ⋅= , T

k ]00100[e ……= , unde 1 apare în poziţia k;

cyLk

' =⋅

kk yzU =⋅

kk zSx ⋅= , n,...,1k = .


Această procedură de determinare a inversei unei matrice necesită un număr de operaţii în virgulă mobilă de ordinul lui 23 nn)3/n( ⋅+ , în cazul utilizării triangularizării simple. Concluzia care se poate desprinde este aceea că inversarea unei matrice necesită un număr mare de operaţii în virgulă mobilă. Ca urmare, în practică, nu se recomandă rezolvarea sistemelor prin metoda bazată pe calculul explicit al inversei matricei sistemului:

bAxbxA 1 ⋅=⇒=⋅ − , deoarece există posibilitatea afectării rezultatului obţinut de către erorile de rotunjire acumulate.

2.5.4 Rezolvarea sistemelor în corpul numerelor complexe

Fie sistemul de ecuaţii:

1nnn b,A,bzA ×× ∈∈=⋅ CC . (2.24)

În principiu, se poate aplica metodologia de rezolvare a sistemelor de ecuaţii având matricea sistemului şi termenul liber cu elemente numere reale, dacă operaţiile cu numere complexe sunt definite (implementate) în limbajul de programare folosit. Altfel, trebuie scrise funcţii sau rutine, stabile din punct de vedere numeric, care să implementeze operaţiile cu numere complexe.

De regulă, problema rezolvării unui sistem complex de ordinul n se transformă în problema rezolvării unui sistem real de dimensiune n2 ⋅ . Astfel, se pot rescrie matricele implicate în (2.24) sub forma următoare, unde 1i2 −= :

nn2121 A,A,AiAA ×ℜ∈⋅+= ,

1n2121 b,b,bibb ×ℜ∈⋅+= , (2.25)

1n2121 z,z,zizz ×ℜ∈⋅+= .

Înlocuind relaţiile (2.25) în (2.24) se obţine:

212121 bib)ziz()AiA( ⋅+=⋅+⋅⋅+ .

Efectuând calculele se obţine:

2121122211 bib)zAzA(i)zAzA( ⋅+=⋅+⋅⋅+⋅−⋅ . (2.26)

Din relaţia (2.26), identificând partea reală şi partea imaginară pentru cei doi membri ai egalităţii, se obţine:

12211 bzAzA =⋅−⋅ ,

22112 bzAzA =⋅+⋅ ,

ceea ce se poate scrie sub formă matricială astfel:

2.6 Metode iterative 57

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡⋅

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡ −

2

1

2

1

12

21

b

b

z

z

AA

AA. (2.27)

Notând cu C matricea de blocuri ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡ −

12

21

AA

AA, cu x vectorul

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

2

1

z

z şi cu

d vectorul ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

2

1

b

b, relaţia (2.27) se reduce la:

dxC =⋅ . (2.28)

Rezolvând sistemul (2.28), se obţine o soluţie care se poate rescrie sub forma: T

n21nn1 ]xxxx[x ⋅+= …… ,

iar soluţia sistemului complex este: T

n21nT

n1 ]xx[i]xx[z ⋅+⋅+= …… .

2.6 Metode iterative

2.6.1 Principiul şi convergenţa metodelor iterative

Fie sistemul de ecuaţii algebrice liniare:

1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅ , (2.29)

în care A este o matrice nesingulară. Metodele iterative se bazează pe construcţia unui şir de aproximaţii ale

soluţiei, ,...1,0k,x ]k[ = , convergent la soluţia adevărată:

]0[]k[

kx,xxlim ∀=

∞→.

Pentru construcţia acestui şir, se consideră rescrierea sau descompunerea matricei A a sistemului sub forma: PNA −= , în care N este o matrice nesingulară. De regulă, se alege matricea N cu o formă simplă. Ecuaţia (2.29) devine:

bxPxNbx)PN( +⋅=⋅⇔=⋅− . (2.30)


Şirul de aproximaţii se construieşte cu ajutorul relaţiei:

,...1,0k,bxPxN ]k[]1k[ =+⋅=⋅ + (2.31)

în care estimaţia iniţială ]0[x este dată (cunoscută). În particular, n]0[ 0x = .

Din relaţia (2.31) rezultă:

,...1,0k,bNxPNx 1]k[1]1k[ =⋅+⋅⋅= −−+ . (2.32)

Se face următoarea notaţie: nn1 G,PNG ×− ℜ∈⋅= .

Matricea G are valori proprii în general complexe, C∈λ )G(i , care formează mulţimea numită spectrul matricei G.

Definiţie:

Se numeşte rază spectrală a matricei G mărimea:

|})G({|max)G( ini1λ=ρ

≤≤.

În cele ce urmează, se enunţă şi demonstrează următorul rezultat.

Propoziţie:

Condiţia necesară şi suficientă ca şirul de soluţii aproximative, defint prin relaţia (2.32), să fie convergent către soluţia adevărată a sistemului de ecuaţii (2.29) este ca matricea PNG 1 ⋅= − să aibă toate valorile proprii în modul subunitare sau, altfel spus, raza spectrală a matricei G să fie subunitară.

Demonstraţia porneşte de la expresia erorii la iteraţia [k], care este: ]k[]k[ xxe −= , ,...1,0k = . Eroarea la pasul [ 1k + ] este:

bNxGxxxe 1]k[]1k[]1k[ ⋅−⋅−=−= −++ . (2.33)

Înlocuind în relaţia (2.33) expresia vectorului b din relaţia (2.30), se obţine:

]k[]k[1]k[]1k[ eG)xx(G)xPxN(NxGxe ⋅=−⋅=⋅−⋅⋅−⋅−= −+ . (2.34)

Exprimând eroarea la pasul [k] în funcţie de eroarea la pasul [k-1], prin folosirea repetată a relaţiei (2.34), se obţine următorul rezultat:

]0[1k]1k[2]k[]1k[ eGeGeGe ⋅==⋅=⋅= +−+ … .

Cum n]0[ 0e ≠ în general, atunci condiţia referitoare la limita şirului de

aproximaţii, prezentată la început, este îndeplinită dacă şi numai dacă


nnk

kGlim ×∞→

= 0 . Această condiţie este satisfăcută dacă matricea G are valorile

proprii subunitare în modul, altfel spus, dacă rază spectrală a matricei G este subunitară.

Observaţii:

1. Cu cât raza spectrală subunitară a matricei G este mai mică, cu atât viteza de convergenţă a şirului de soluţii aproximative (2.32) va fi mai mare.

2. În practică, de multe ori, condiţia necesară şi suficientă prezentată anterior se verifică (înlocuieşte) printr-o condiţie suficientă, dacă este posibil, şi anume:

dacă 1||G|| <α atunci 1)G( <ρ .

De regulă se foloseşte norma matricială infinit, rezultând condiţia:

1}|g|{max||G||n

1jj,ini1

<= ∑=≤≤

∞ .

Dacă această ultimă condiţie este îndeplinită, atunci metoda iterativă este sigur convergentă şi nu mai este necesar să se calculeze valorile proprii ale matricei G de caracterizare a convergenţei. Dacă, însă, condiţia suficientă nu este satisfăcută, atunci nu se poate afirma nimic în ceea ce priveşte convergenţa metodei şi se recurge la calculul valorilor proprii ale matricei G, în scopul verificării condiţiei necesare şi suficiente.

În continuare se consideră următoarea descompunere a matrcii A a sistemului (2.29):

UDLA ++= , în care L este o matrice inferior triunghiulară cu diagonala principală nulă având elementele de sub diagonala principală egale cu elementele matricei A de acelaşi rang, D este o matrice diagonală având elementele de pe diagonala principală egale cu elementele de pe diaonala principală a matricei A, iar U este o matrice superior triunghiulară având elementele de deasupra diagonalei principale egale cu elementele matricei A de acelaşi rang.

Se mai face presupunerea că elementele diagonalei principale a matricei A sunt nenule. În caz contrar, se pot face permutări de linii şi/sau coloane astfel încât să fie îndeplinită şi această condiţie.

2.6.2 Metoda Jacobi şi metoda Gauss-Seidel

În cazul metodei Jacobi, matricele N şi P sunt:

)UL(P,DN +−== . (2.35)


Din relaţia (2.31), folosind relaţia (2.35), rezultă:

,...1,0k,bx)UL(xD ]k[]1k[ =+⋅+−=⋅ + . (2.36)

Relaţia (2.36) se poate scrie pe linii astfel:

n,...,1i,bxaxa i

n

ij1j

]k[jj,i

]1k[ii,i =+⋅−=⋅ ∑

≠=

+ . (2.37)

Dacă elementul 0a i,i ≠ , atunci relaţia (2.37) se poate rescrie sub forma:

n,...,1i,x)a/a()a/b(xn

ij1j

]k[ji,ij,ii,ii

]1k[i =⋅−= ∑

≠=

+ . (2.38)

Matricea G corespunzătoare metodei Jacobi este:

nj,i1j,i11

Jacobi ]g[)UL(DPNG ≤≤−− =+⋅−=⋅= ,

în care ⎩⎨⎧

≠−=

=ji,a/aji,0

gi,ij,i

j,i .

Condiţia suficientă care se impune pentru ca metoda Jacobi să fie convergentă este:

1}|a/a|{max}|g|{maxn

ij1j

i,ij,ini1

n

1jj,ini1

<= ∑∑≠=≤≤=≤≤

. (2.39)

În relaţia (2.39) fiind implicate doar numere pozitive, dacă maximul lor este mai mic decât 1, atunci toate sunt subunitare:

∑∑≠=

≠=

>⇔<n

ij1j

j,ii,i

n

ij1j

i,ij,i |a||a|1|a/a| . (2.40)

O matrice care satisface relaţia (2.40) se numeşte matrice diagonal dominantă pe linii.

Propoziţie:

Dacă matricea A este diagonal dominantă pe linii, atunci metoda Jacobi este convergentă, oricare ar fi estimaţia iniţială a soluţiei sistemului de ecuaţii (2.29).

Observaţie:


Condiţia (2.40) înseamnă că ji,n,...,1j,n,...,1i ≠==∀ , 1|a/a| i,ij,i < . Revenind la relaţia (2.38), se observă că i,ij,i a/a sunt coeficienţii care

multiplică estimaţiile anterioare, deci erorile ce afectează aceste estimaţii sunt micşorate pe măsură ce procesul iterativ avansează. Ca urmare, dacă matricea A este diagonal dominantă pe linii, atunci procedura este sigur stabilă numeric.

Metoda Gauss-Seidel porneşte de la descompunerea PNA −= în care:

UP,DLN −=+= . Relaţia (2.31) se poate rescrie, în acest caz, sub forma:

,...1,0k,bxUx)DL( ]k[]1k[ =+⋅−=⋅+ + . (2.41)

Relaţia (2.41) se poate scrie, pe linii, sub forma:

n,...,1i,bxaxa i

n

1ij

]k[jj,i

i

1j

]1k[jj,i =+⋅−=⋅ ∑∑

+==

+ ,

din care rezultă, dacă 0a i,i ≠ :

n,...,1i,x)a/a(x)a/a()a/b(xn

1ij

]k[ji,ij,i

1i

1j

]1k[ji,ij,ii,ii

]1k[i =⋅+⋅−= ∑∑

+=

−

=

++ .

Se poate demonstra că, şi în cazul metodei Gauss-Seidel, dacă matricea A este diagonal dominantă pe linii, atunci metoda este convergentă.

În general, se demonstrează că între raza spectrală subunitară a matricei )UL(DG 1

Jacobi +⋅−= − şi raza spectrală subunitară a matricei U)DL(G 1

SeidelGauss ⋅+−= −− există relaţia:

1)G()G( SeidelGaussJacobi2 <ρ≅ρ − .

Rezultă, în general, că dacă metoda Gauss-Seidel este convergentă, atunci viteza sa de convergenţă este mai mare decât cea corespunzătoare metodei Jacobi.

Există situaţii în care, pentru ambele metode, razele spectrale sunt subunitare, dar apropiate de valoarea 1. În acest caz, convergenţa ambelor metode este extrem de lentă, recomandându-se ca în descompunerea matricei A să se utilizeze unul sau mai mulţi parametri de accelerare a convergenţei.

2.6.3 Metoda relaxărilor succesive

Se porneşte de la descompunerea matricei A corespunzătoare metodei Gauss-Seidel. Relaţia (2.41) se poate scrie sub forma:


,...1,0k,b)DL(xU)DL(x 1]k[1]1k[ =⋅++⋅⋅+−= −−+ . În acestă relaţie, în membrul drept se adună şi se scade cantitatea ]k[x .

Ţinând cont de faptul că )DL()DL(I 1n +⋅+= − rezultă:

]bxA[)DL(x

]bx)DL(xU[)DL(xx]k[1]k[

]k[]k[1]k[]1k[

−⋅⋅+−=

−⋅++⋅⋅+−=−

−+

Se face următoarea notaţie: ]k[]k[ xAbr ⋅−= ,

în care ]k[r se numeşte reziduu corespunzător iteraţiei [k]. În acest caz, se poate scrie:

]k[1]k[]1k[ r)DL(xx ⋅++= −+ . (2.42)

Se transformă relaţia (2.42), înmulţind termenul care conţine reziduul cu un parametru ω de accelerare a convergenţei:

]k[1]k[]1k[ r)DL(xx ⋅+⋅ω+= −+ . (2.43)

Al doilea termen din suma exprimată în relaţia (2.43) poate fi interpretat ca un factor de corecţie pentru estimaţia de la iteraţia [k]. Astfel, soluţia de la iteraţia [ 1k + ] se obţine prin corectarea soluţiei de la iteraţia anterioară [k], corecţia realizându-se în funcţie de reziduul ]k[r .

Relaţia (2.43) se poate rescrie sub forma:

]b)DL(xU)DL([

x)DL()DL(xx1]k[1

]k[1]k[]1k[

⋅++⋅⋅+−⋅ω+

⋅+⋅+⋅ω−=−−

−+

(2.44)

Paranteza dreptunghiulară pune în evidenţă aproximaţia la iteraţia [ 1k + ], obţinută prin metoda Gauss-Seidel. Relaţia (2.44) devine:

]1k[SeidelGauss

]k[]1k[ xx)1(x +−

+ ⋅ω+⋅ω−= . (2.45)

Se demonstrează că şirul de aproximaţii obţinut cu ajutorul relaţiei de recurenţă (2.45) este convergent pentru )2,0(∈ω . Dacă )1,0(∈ω metoda se numeşte a subrelaxărilor succesive, iar dacă )2,1(∈ω ea se numeşte a suprarelaxărilor succesive. Dacă 1=ω , se obţine metoda Gauss-Seidel.

În general, dacă metoda este convergentă, atunci numărul de iteraţii ca funcţie de parametrul ω atinge un punct de minim corespunzător unei valori

optimω (a se vedea Figura 2.3). Se poate arăta că optimω şi optimρ au aproximativ

valorile:


411;

11

2 2optim

SeidelGauss2Jacobi

Jacobioptim2

Jacobi

optim

ω⋅ρ≅

ρ−+

ρ≅ρ

ρ−+≅ω − .

În practică, la nivelul fiecărei componente, se pot folosi relaţiile de calcul (2.45):

]a/)xaxab[(x)1(x i,i

n

1ij

]k[jj,i

1i

1j

]1k[jj,ii

]k[i

]1k[i ∑∑

+=

−

=

++ ⋅−⋅−⋅ω+⋅ω−= .

Relaţiile prezentate corespund unei descompuneri a matricei A care depinde de parametrul ω , şi anume:

)(P)(NAbxA;bxA ω−ω=⋅ω⇒⋅ω=⋅⋅ω=⋅ în care DL)(N +⋅ω=ω şi UD)1()(P ⋅ω−⋅ω−=ω .

numar iteratii

ω

ωoptim1 2

numar minimiteratii

Fig. 2.3 Metoda suprarelaxării succesive

***

Algoritmii prezentaţi lucrează până când este îndeplinită o condiţie de tipul:

x]1s[]s[ ||xx|| ε≤− α

− .

Se foloseşte, de regulă, norma infinit şi atunci condiţia de stop este:

x]1s[

i]s[

ini1|}xx{|max ε≤− −

≤≤.

Notând ]1k[n xx += şi ]k[

v xx = , atunci următorul algoritm descrie procedura generală aplicată în cazul unei metode iterative:

citeşte ]0[x x,ε

atribuie ]0[n xx ←

atribuie 1vn ←∞


atribuie 0iter ← cât timp )vn( xε>∞ execută ⎡ atribuie nv xx ← ⎢ atribuie 1iteriter +← ⎢ * calculează nx ⎢ * calculează |})i(x)i(x{|maxvn vnni1

−=≤≤

∞

⎣

În final, se obţin următoarele: numărul de iteraţii (iter) şi vectorul nx care aproximează soluţia sistemului cu precizia impusă xε . În descrierea anterioară,

]0[x reprezintă estimaţia iniţială a soluţiei, iar ∞vn este notaţia pentru valoarea normei vectoriale infinit.

2.7 Caracterizarea problemei rezolvării unui sistem de ecuaţii algebrice liniare. Precizia soluţiei calculate

Precizia soluţiei calculate a problemei: 1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅

depinde critic de buna sa condiţionare numerică. Această proprietate este caracterizată de numărul de condiţie al matricei A a sistemului. Aceste concepte fac obiectul acestui subcapitol.

2.7.1 Numărul de condiţie al unei matrice

Definiţie: Oricare ar fi matricea nnA ×ℜ∈ , se defineşte numărul său de condiţie în

raport cu norma vectorială α , notat )A(kα , ca fiind: m/M)A(k =α ,

}||x||/||xA{||minm},||x||/||xA{||maxMnn 0x0x αα≠αα

≠⋅=⋅= ,

unde α||.|| este una din normele vectoriale uzuale (a se vedea Anexa A).

Cum în general mM ≥ , rezultă că şi 1)A(k ≥α . Se poate demonstra că

α−= ||A||m/1 1 . Astfel, se poate defini numărul de condiţie în raport cu

operaţia de inversare matricială ca fiind:

α−

αα ⋅= ||A||||A||)A(k 1 .

2.7 Caracterizarea problemei rezolvării unui sistem de ecuaţii algebrice liniare. Precizia soluţiei calculate 65

Numărul de condiţie depinde în general de norma matricială folosită. Între diferitele numere de condiţie corespunzătoare aceleiaşi matrice există relaţii care reflectă relaţiile dintre normele matriciale echivalente utilizate. Astfel, dacă α şi β sunt două norme matriciale (vectoriale), există constantele ℜ∈21 c,c astfel încât să aibă loc relaţia:

)A(kc)A(k)A(kc 21 αβα ⋅≤≤⋅ .

În multe abordări, este folosită norma infinit pentru a calcula numărul de condiţie al unei matrice:

∞−

∞∞ ⋅= ||A||||A||)A(k 1 . Cum calculul normei infinit este imediat (a se vedea Anexa A), problema

este de a determina norma infinit pentru inversă. O soluţie o reprezintă calculul inversei matricei A prin triangularizare şi apoi calculul normei infinit. În acest fel, se obţine o aproximare a numărului de condiţie al matricei. Dezavantajul acestei modalităţi este că este necesar un volum mare de calcule.

Discuţia asupra numărului de condiţie al unei matrice este reluată în capitolul 5 destinat studiului valorilor singulare ale unei matrice. Se prezintă acolo modalitatea uzual folosită pentru calculul eficient al numărului de condiţie.

Exemplul 2.4:

• dacă nIA = (matrice unitate de ordinul n), atunci 1)A(k =α ; • dacă kPA = (matrice elementară de permutare de linii), atunci 1)A(k =α ; • dacă kSA = (matrice elementară de permutare de coloane), atunci

1)A(k =α ; • dacă }d,...,d,...,d{diagA nnii11= (matrice diagonală), atunci

|}d{|min/|}d{|max)A(k iini1iini1 ≤≤≤≤α = .

În funcţie de valorile numărului de condiţie, matricea A poate fi caracterizată după cum urmează: a) dacă numărul de condiţie este egal cu 1, matricea A se spune că este foarte

bine condiţionată; b) dacă numărul de condiţie este mic, apropiat de valoarea 1, atunci matricea

A este bine condiţionată; c) dacă numărul de condiţie este mare, atunci matricea A este prost

condiţionată; d) dacă numărul de condiţie este foarte mare, atunci matricea A este aproape

singulară; e) dacă numărul de condiţie este infinit, atunci matricea A este singulară.


Se mai utilizează şi inversul numărului de condiţie: )A(k/1)A(rk αα = .

În funcţie de valorile pe care le poate lua acesta, matricea A poate fi caracterizată după cum urmează, unde mε este notaţia pentru epsilonul-maşină: a) matrice foarte bine condiţionată, dacă 1)A(rk =α ; b) matrice bine condiţionată, dacă 1)A(rk <α , dar apropiat de valoarea 1:

1)A(rkm <<<ε α ; c) matrice prost condiţionată (aproape singulară) dacă m)A(rk ε≤α ; d) matrice singulară (foarte prost condiţionată) dacă 0)A(rk =α . Conform acestor caracterizări se poate afirma, în general, că numărul de condiţie al unei matrice reprezintă inversul distanţei dintre matricea A şi mulţimea matricelor singulare.

Exemplul 2.5:

De regulă se preferă caracterizarea singularităţii pe baza numărului de condiţie sau a inversului său, decât folosirea valorii determinantului. Fie matricea:

1001001.00

01.0D

×⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

…

….

Atunci 010)Ddet( 100 ≅= − în simplă precizie, deşi 1)D(k =∞ , ceea ce arată că matricea D este o matrice foarte bine condiţionată.

Exemplul 2.6:

În general, nu există o corelaţie între numărul de condiţie şi valoarea determinantului unei matrice. Fie matricea:

nj,i1j,i

nn

n ]b[

1001

0111

B ≤≤

×

=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−

−−

=

…

…

cu ⎪⎩

⎪⎨

⎧=

>=<−

= n,...,1j,i,ji,0ji,1

ji,1b j,i .

Atunci 1)Bdet( n = şi nn 2n)B(k ⋅=∞ . Pentru n având valori mari, matricea

nB este prost condiţionată. Proasta condiţionare a unei matrice se poate reflecta în rezultatele triangularizării sale. Astfel, dacă în urma triangularizării unei matrice, matricea U superior triunghiulară are pe diagonala principală elemente nule sau foarte mici, atunci sigur matricea iniţială este prost condiţionată. Reciproca, în general, nu este valabilă: matricea poate fi prost condiţionată, însă


datorită erorilor de calcul, elementele de pe diagonala matricei U pot să nu fie neapărat foarte mici.

2.7.2 Precizia soluţiei calculate

Fie sistemul de ecuaţii: 1nnn b,A,bxA ×× ℜ∈ℜ∈=⋅ . (2.46)

Caracterizarea (condiţionarea) acestei probleme este dată de numărul de condiţie al matricei A, aşa cum a fost el definit la începutul acestui capitol. Pot exista două situaţii generale prezentate în continuare: I. Soluţia este calculată fără erori (exact), dar pot exista erori în matricea

A şi termenul liber b . Se disting următoarele două subcazuri, I.1 şi I.2, prezentate în continuare.

I.1 Se presupune că termenul liber al sistemului (2.46) este uşor perturbat, lucrându-se cu vectorul b~ în loc de vectorul b . De fapt, se rezolvă (exact) problema:

b~x~A =⋅ . Dacă x~ este soluţia exactă a acestei noi probleme, atunci se demonstrează că eroarea relativă satisface la relaţia:

bx )A(k||b~||/||b~b||)A(k||x||/||x~x|| ε⋅=−⋅≤−=ε αααααα .

Rezultă, aşadar, că eroarea relativă în termenul liber este amplificată în soluţia calculată de )A(kα ori, în absenţa erorilor de rotunjire.

Exemplul 2.7:

Fie sistemul (2.46) cu următoarele matrice:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=

70.911.4

b~,7.91.4

b,6.67.98.21.4

A .

Soluţia exactă a sistemului bxA =⋅ este ⎥⎦

⎤⎢⎣

⎡=

01

x . Soluţia exactă a sistemului

b~x~A =⋅ este x97.034.0

x~ ≠⎥⎦

⎤⎢⎣

⎡= .

Realizeazând următoarele notaţii: b~bb −=Δ , x~xx −=Δ , rezultă pentru acest exemplu:

4.2294)||b||/||b/(||)||x||/||x(||)A(k 11111 =ΔΔ≥ ,


deci matricea A este prost condiţionată.

I.2 De această dată, se consideră matricea sistemului ca fiind afectată de perturbaţii:

αα <<+= ||A||||E||,EAA~ . De fapt, se rezolvă exact problema:

bx~A~ =⋅ . (2.47)

Dacă x~ este soluţia exactă a ecuaţiei (2.47), atunci se demonstrează că eroarea relativă satisface la relaţia:

Ax )A(k)||A||/||E(||)A(k ε⋅=⋅≤ε αααα .

Rezultă faptul că perturbaţiile în matricea sistemului, în ipoteza absenţei erorilor de rotunjire la calculul soluţiei, se regăsesc în soluţia calculată amplificate de

)A(kα ori.

II. Datele iniţiale ale problemei (2.46), A şi b , sunt exacte, dar procedura de calcul a soluţiei este afectată de eroare (erori de rotunjire). În acest caz se disting două situaţii, II.1 şi II.2, prezentate în cele ce urmează.

II.1 Apar erori datorate decompunerii gaussiene. Se obţin următoarele rezultate: • ULA ⋅= , în cazul triangularizării simple; • ULAP ' ⋅=⋅ , în cazul triangularizării cu pivotare parţială; • ULSAP ' ⋅=⋅⋅ , în cazul triangularizării cu pivotare totală.

Plecând de la matricea A, căreia i se aplică o procedură de triangularizare, se obţin matricele L ( 'L ), U, P, (S). Dacă se efectuează operaţia inversă, se obţine:

• A~UL =⋅ , corespunzător triangularizării simple; • A~ULP '1 =⋅⋅− , corespunzător triangularizării cu pivotare parţială; • A~SULP 1'1 =⋅⋅⋅ −− , corespunzător triangularizării cu pivotare totală,

unde A~ este matricea iniţială plus o matrice de eroare: cEAA~ += . În continuare, discuţia se referă la descompunerea cu pivotare parţială şi la

cea cu pivotare totală, deoarece procedurile corespunzătoare lor se bazează pe matrice stabilizate, iar algoritmii sunt stabili din punct de vedere numeric.

Având in vedere volumul de calcule implicat (numărul de operaţii în virgulă mobilă), se poate afirma că dacă ordinul n al sistemului este mic, atunci:

αα << ||A||||E|| c , iar dacă ordinul n este mare, atunci: αα < ||A||||E|| c . Aşadar, norma matricei de eroare cE se poate apropia de cea a matricei A a sistemului.


Ca urmare, o posibilitate de a caracteriza precizia descompunerii este de a calcula raportul:

c,Ac ||A||/||AA~||||A||/||E|| ε=−= αααα . De regulă se foloseşte norma 1. În general, se demonstrează că:

c,Ax )A(k ε⋅≤ε α .

Dacă 0d,10 dc,A >≅ε − , atunci se poate afirma că triangularizarea matricei

A se face cu d cifre zecimale exacte. Rezultă că eroarea în soluţia calculată satisface la relaţia:

)pd(dpc,Ax 1010)A(k −−−

α =≤ε⋅≤ε .

O matrice foarte bine condiţionată este caracterizată de 0p = . Dacă matricea este bine condiţionată, atunci 21p ÷= , iar dacă ea este prost condiţionată, atunci 2p > . Se desprinde următoarea concluzie: dacă matricea A este prost condiţionată, atunci se pierd p cifre ca precizie în soluţia calculată.

Dacă, însă, matricea A este diagonal dominantă pe linii sau coloane, atunci ea este bine condiţionată. Matricea A poate fi adusă la o astfel de formă înainte de a aplica procedura de descompunere, dacă este posibil, prin permutări adecvate de linii şi/sau coloane.

II.2 Erorile afectează procedura de rezolvare prpriu-zisă a sistemului (fazele de substituţie înainte şi înapoi).

În acest caz, reziduul asociat soluţiei calculate x~ este: x~Abr ⋅−= . Se demonstrează că, în această situaţie, eroarea relativă în soluţia calculată satisface la relaţia:

)||b||/||r(||)A(kx ααα ⋅≤ε .

În concluzie, un reziduu mic în normă faţă de norma termenului liber nu garantează o precizie bună a soluţiei calculate, dacă matricea sistemului este prost condiţionată.

Exemplul 2.8:

Fie sistemul bxA =⋅ cu:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=

2001.2

b,11001.11

A .

Soluţia exactă a acestui sistem este T]11[x = . În urma calculelor rezultă: x]02[x~ T ≠= , T]0001.0[xAbr =⋅−= . Aceasta se datorează faptului că:


311111 10)||b||/||r/(||)||x||/||x~x(||)A(k ≈−≥ ,

deci matricea A este prost condiţionată.

O concluzie generală care se poate desprinde este următoarea: buna condiţionare a problemei de calcul ( 100101)A(k ÷÷=α ), împreună cu stabilitatea numerică a algoritmului de triangularizare ( αα << ||A||||E|| c ), garantează faptul că eroarea relativă a soluţiei calculate faţă de soluţia exactă este de acelaşi ordin de mărime cu erorile relative în datele iniţiale ale problemei (matricea A şi termenul liber b ).

O procedură folosită pentru estimarea numărului de cifre exacte ale soluţiei calculate este următoarea:

• triangularizarea matricei A; • determinarea soluţiei x~ folosind rezultatele triangularizării matricei A: • calculul reziduului asociat x~Abr ⋅−= ; • rezolvarea sistemului reA =⋅ , folosind rezultatele procedurii de

triangularizare; • dacă 0||e|| 1≠ , atunci d

11 10||x~||/||e|| −≅ . Numărul de cifre exacte ale soluţiei calculate este:

)]||x~||/||elg(||[d 11−≅ ,

unde lg(.) semnifică funcţia logaritm zecimal, iar [.] semnifică funcţia parte întreagă.

O procedură de îmbunătăţire (rafinare) a soluţiei calculate este următoarea: 1. rezolvare sistem bxA =⋅ folosind descompunerea L-U a matricei A şi

determinarea soluţiei x~ ; 2. calcul reziduu asociat: x~Abr ⋅−= ; 3. rezolvare sistem reA =⋅ folosind rezultatele descompunerii L-U a

matricei A şi determinarea lui e ;

4. rafinarea soluţiei: ex~x~~ += ;

5. dacă 11 ||x~~||||x~|| ≠ atunci x

~~x~ = şi reluare de la pasul 2, altfel x~x = . 6. stop

Îmbunătăţirea soluţiei se face într-un număr maxim de iteraţii notat ITMAX. Dacă după ITMAX iteraţii, testul de la pasul 5 continuă să fie satisfăcut (şi se reia de la pasul 2), acesta este un semn că matricea A este prost condiţionată deoarece:

2.8 Exerciţii propuse 71

)||b||/||r(||)A(k||x~||/||e|| 11111 ⋅≤ .

Astfel, dacă numărul de condiţie )A(k1 este mare, atunci şi raportul

11 ||x~||/||e|| este mare şi, ca urmare, valoarea 1||e|| nu se va micşora dacă matricea A este prost condiţionată.

2.8 Exerciţii propuse

E2.1 Să se realizeze, în mediul de programare MATLAB, un program pentru rezolvarea unui sistem de n ecuaţii liniare cu n necunoscute, prin metoda directă folosind triangularizarea matricei sistemului cu pivotare totală. Programul va compara soluţia găsită cu cea furnizată de mediul MATLAB.

În realizarea programului, se vor parcurge următoarele etape descrise parţial în limbajul pseudocod.

1. Introducere parametru de control al execuţiei procedurii de triangularizare, în variabila EPS.

2. Introducere ordin sistem, în variabila n, 2n ≥ .

3. Alocare spaţiu de memorie pentru matricea sistemului, în variabila a, prin iniţializarea acesteia cu o matrice nulă de dimensiune n. Introducerea elementelor matricei sistemului se face pe linii (a se vedea observaţia 2 care urmează). Salvarea valorilor introduse pentru matricea a, în matricea de lucru aa:

atribuie aa ← a

Reluare introducere, în caz de eroare.

4. Alocare spaţiu de memorie pentru termenul liber al sistemului, în variabila b, prin iniţializarea acesteia cu un vector coloană nul, cu n componente. Introducerea elementelor termenului liber se realizează pe linii. Salvarea valorilor din vectorul b, în vectorul de lucru bb:

atribuie bb ← b


5. Faza I-a: eliminare înainte cu pivotare totală:

pentru i = 1,n execută


⎡ atribuie iord(i) ← i ⎣

* pivotare totală:

pentru k = 1,n-1 execută ⎡ atribuie [mc,imc] ← max ( abs ( a(k:n,k:n) ) ) ⎢ atribuie [ml,jk] ← max ( mc ) ⎢ atribuie ik ← imc ( jk ) ⎢ atribuie ik ← ik + (k-1) ⎢ atribuie jk ← jk + (k-1) ⎢ dacă ( ik ~= k ) atunci ⎢ ⎡ atribuie wmx ← zeros (1,n) ⎢ ⎢ atribuie wmx ← a(k,:) ⎢ ⎢ atribuie a(k,:) ← a(ik,:) ⎢ ⎢ atribuie a(ik,:) ← wmx ⎢ ⎢ atribuie t ← b(k) ⎢ ⎢ atribuie b(k) ← b(ik) ⎢ ⎢ atribuie b(ik) ← t ⎢ ⎣ ⎢ dacă ( jk ~= k ) atunci ⎢ ⎡ atribuie zmx ← zeros (n,1) ⎢ ⎢ atribuie zmx ← a(:,k) ⎢ ⎢ atribuie a(:,k) ← a(:,jk) ⎢ ⎢ atribuie a(:,jk) ← zmx ⎢ ⎢ atribuie it ← iord(k) ⎢ ⎢ atribuie iord(k) ← iord(jk) ⎢ ⎢ atribuie iord(jk) ← it ⎢ ⎣ ⎢ * triangularizare: ⎢ dacă ( abs ( a(k,k) ) > EPS ) atunci ⎢ ⎡ atribuie a(k+1:n,k) ← a(k+1:n,k) / a(k,k) ⎢ ⎢ atribuie a(k+1:n,k+1:n) ← a(k+1:n,k+1:n) – ⎢ ⎢ a(k+1:n,k) * a(k,k+1:n) ⎢ ⎢ atribuie b(k+1:n) ← b(k+1:n) – a(k+1:n,k) * b(k) ⎢ ⎢ atribuie a(k+1:n,k) ← zeros (n-k,1) ⎢ ⎢altfel ⎢ ⎢ scrie ‘pivot nul sau foarte mic’ ⎢ ⎢ scrie ‘STOP: algoritm’ ⎢ ⎢ * STOP program (funcţia MATLAB return) ⎢ ⎣


⎢ scrie ‘k = ’, k ⎢ scrie ‘a = ’, a ⎢ scrie ‘b = ’, b ⎣

6. Faza a-II-a: calcul soluţie sistem:

* substituţie înapoi:

pentru k = 1,n execută ⎡ dacă ( abs ( a(i,i) ) < EPS ) atunci ⎢ ⎡ scrie ‘matrice prost condiţionată’ ⎢ ⎣ ⎣

atribuie b(n) ← b(n) / a(n,n) pentru i = n-1,(pas = -1),1 execută ⎡ atribuie sum ← a(i,i+1:n) * b(i+1:n) ⎢ atribuie b(i) ← ( b(i) – sum ) / a(i,i) ⎣

* ordonare soluţie:

atribuie x ← zeros (n,1) pentru i = 1,n execută ⎡ atribuie x ( iord(i) ) ← b(i) ⎣

7. Calcul soluţie folosind funcţiile MATLAB, pentru comparaţie:

atribuie x_bs ← aa \ bb atribuie x_inv ← inv (aa) * bb

8. Afişare soluţii x, x_bs, x_inv ( mod de afişare MATLAB dat de funcţia: format long e; ).

9. Calcul reziduuri asociate soluţiilor calculate, precum şi norme euclideene:

atribuie r ← bb – aa * x atribuie nr ← norm (r,2)

atribuie r_bs ← bb – aa * x_bs atribuie nr_bs ← norm (r_bs,2)

atribuie r_inv ← bb – aa * x_inv atribuie nr_inv ← norm (r_inv,2)


Afişare rezultate n, nr, r_bs, nr_bs, r_inv, nr_inv ( mod de afişare MATLAB dat de funcţia: format long e; ).

10. Reluare program, dacă se doreşte, cu alt termen liber (punctul 4). Dacă are loc reluarea, înainte de aceasta se reface matricea a, utilizând valorile salvate în matricea de lucru aa:

atribuie a ← aa

11. Reluare program, dacă se doreşte, cu alt sistem de acelaşi ordin n (punctul 3).

12. Reluare program, dacă se doreşte, cu un sistem de alt ordin n (punctul 2).

13. Reluare program, dacă se doreşte, cu alt parametru EPS (punctul 1).

Observaţii:

1.) Descrierea anterioară a programului foloseşte următoarele funcţii MATLAB: max (determinare valoare maximă dintr-un tablou de numere), zeros (iniţializare tablou cu elemente nule), abs (calcul valoare absolută), return (părăsire comandă structurată, de exemplu ciclu cu contor), \ (împărţire matricială la stânga, rezolvată prin triangularizarea cu pivotare parţială a matricei deîmpărţit), inv (calcul inversă matrice pătratică), norm (calcul normă vectorială/matricială).

2.) Se recomandă următoarea modalitate generală de introducere a unei matrice, pe linii (n linii, m coloane), în variabila a:

* secvenţă de comenzi MATLAB: a = zeros (n,m); for i = 1:n, fprintf ( ‘Linia %.0f : ’, i ); a(i,:) = input ( ‘ ’ ); end;

* exemplu de introducere (n=m=3):

Linia 1: [1 2 3] <RETURN>




În descrierea anterioară, s-au subliniat mesajele scrise de către program iar <RETURN> semnifică apăsarea tastei respective.

3.) Programul se va rula pentru următoarele date de intrare, rezulatele analizându-se şi comparându-se:

EPS = 1.e-10

I. n = 3; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

−=

6901.37

b;5156099.230710

a ;

II. n = 3; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−=

2102105

b;121

1003110023

a ;

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−=

00.102.105.1

b;50.000.150.0

00.103.001.000.102.003.0

a

(acelaşi exemplu, cu scalare pe linii)

III. n = 3; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−

−=

275

b;11141

142321

a ;

IV. n = 4; ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−

=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−

=

672

7

b;

5616103423221020

a .

E2.2 Să se realizeze, în mediul de programare MATLAB, un program pentru rezolvarea unui sistem de n ecuaţii liniare cu n necunoscute, prin metoda iterativă Jacobi. Programul va compara soluţia găsită cu cea furnizată de mediul MATLAB, aplicând o metodă directă de rezolvare.

În realizarea programului, se vor parcurge următoarele etape descrise parţial în limbajul pseudocod.

1. Introducere parametru de control al execuţiei procedurii iterative de rezolvare, în variabila EPS. Introducere număr maxim de iteraţii, în variabila max_iter.

2. Introducere ordin sistem, în variabila n, 2n ≥ .


3. Alocare spaţiu de memorie pentru matricea sistemului, în variabila a, prin iniţializarea acesteia cu o matrice nulă de dimensiune n. Introducerea elementelor matricei sistemului se face pe linii (a se vedea observaţia 2 de la exerciţiul E2.1). Testarea elementelor de pe diagonala matricei a:

pentru k = 1,n execută ⎡ dacă ( abs ( a(i,i) ) < EPS ) atunci ⎢ ⎡ scrie ‘Element pe diagonală nul sau foarte mic’ ⎢ ⎢ scrie ‘Rearanjaţi şi reintroduceţi matricea a’ ⎢ ⎣ ⎣


4. Alocare spaţiu de memorie pentru termenul liber al sistemului, în variabila b, prin iniţializarea acesteia cu un vector coloană nul, cu n componente. Introducerea elementelor termenului liber se realizează pe linii. Reluare introducere, în caz de eroare.

5. Alocare spaţiu de memorie pentru aproximaţia iniţială a soluţiei, în variabila xn, prin iniţializarea acesteia cu un vector coloană nul, cu n componente. Introducerea elementelor vectorului xn se realizează pe linii. Reluare introducere, în caz de eroare.

6. Studiu de convergenţă:

atribuie nn ← diag ( diag (a) ) atribuie p ← nn - a atribuie g ← inv (nn) * p atribuie valp ← eig (g) atribuie ro ← max ( abs (valp) ) scrie ‘=> rază spectrală matrice convergenţă, ro = ’, ro

dacă ( ro < 1 ) atunci ⎡ scrie ‘metoda converge’ ⎢ altfel ⎢ scrie ‘metoda NU converge!’ ⎣

7. Calcul soluţie iterativă:

atribuie vninf ← 1 atribuie iter ← 0

* comutare mod de afişare (limbaj MATLAB): format long e;


cât timp ( (vninf > EPS) şi (iter < max_iter) ) execută ⎡ atribuie iter ← iter + 1 ⎢ atribuie xv ← zeros (n,1) ⎢ atribuie xv ← xn ⎢ * adaptare soluţie metoda Jacobi: ⎢ pentru i = 1,n execută ⎢ ⎡ atribuie sum ← 0 ⎢ ⎢ pentru j = 1,i-1 execută ⎢ ⎢ ⎡ atribuie sum ← sum + a(i,j) * xv(j) ⎢ ⎢ ⎣ ⎢ ⎢ pentru j = i+1,n execută ⎢ ⎢ ⎡ atribuie sum ← sum + a(i,j) * xv(j) ⎢ ⎢ ⎣ ⎢ ⎢ atribuie xn(i) ← ( b(i) - sum ) / a(i,i) ⎢ ⎣ ⎢ * sfârşit adaptare soluţie metoda Jacobi ⎢ atribuie vninf ← max ( abs (xn - xv) ) ⎢ scrie ‘* iter = ’, iter ⎢ scrie ‘xn = ’, xn ⎢ scrie ‘vninf = ’, vninf ⎣

8. Calcul soluţie folosind o metodă directă, în mediul MATLAB:

atribuie x ← a \ b

9. Afişare rezultate: iter, xn, x. Comutare mod de afişare (limbaj MATLAB): format short;

10. Reluare program, dacă se doreşte, cu altă estimaţie iniţială a soluţiei (punctul 5).

11. Reluare program, dacă se doreşte, cu alt termen liber al sistemului (punctul 4).

12. Reluare program, dacă se doreşte, cu alt sistem de acelaşi ordin (punctul 3).

13. Reluare program, dacă se doreşte, cu alt sistem de alt ordin n (punctul 2).

14. Reluare program, dacă se doreşte, cu alţi parametri EPS şi max_iter (punctul 1).

Observaţii:


1.) Descrierea anterioară a programului foloseşte următoarele funcţii MATLAB: abs (calcul valoare absolută), diag (producerea unei matrice diagonale; extragerea diagonalei principale a unei matrice), inv (calcul inversă matrice pătratică), eig (calcul valori proprii ale unei matrice), max (determinare valoare maximă dintr-un tablou de numere), zeros (iniţializare tablou cu elemente nule), \ (împărţire matricială la stânga, rezolvată prin triangularizarea cu pivotare parţială a matricei deîmpărţit).

2.) Operatorul logic şi este implementat în limbajul MATLAB prin simbolul: &

3.) Programul se va rula pentru următoarele date de intrare, rezulatele analizându-se şi comparându-se:

EPS = 1.e-10; max_iter = 30

I. n = 3; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−

−=

649

b;510151

015a ;

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

000

xn ;

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

122

xn ;

II. n = 3; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−

−=

4128

b;271912118

a ; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

000

xn ;

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−

−=

4128

b;721

192118

a ; ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

000

xn .

E2.3 Să se realizeze, în mediul de programare MATLAB, un program pentru rezolvarea unui sistem de n ecuaţii liniare cu n necunoscute, prin metoda iterativă Gauss-Seidel. Programul va compara soluţia găsită cu cea furnizată de mediul MATLAB, aplicând o metodă directă de rezolvare.

În realizarea programului, se vor parcurge exact aceleaşi etape descrise în cadrul exerciţiului E2.2. Deosebirile constau în modul de calcul al matricelor nn şi p, care definesc matricea de conveegenţă g (punctul 6), precum şi în modalitatea de adaptare a soluţiei iterative (punctul 7) de la exerciţiul E2.2.

2.9 Probleme propuse 79

Matricele nn şi p se calculează cu comenzile:

atribuie nn ← tril (a) atribuie p ← nn - a

Porţiunea din descrierea în limbaj pseudocod (punctul 7) de la exrciţiul E2.2: ⎢ * adaptare soluţie metoda Jacobi: ⎢ .................................................................... ⎢ * sfârşit adaptare soluţie metoda Jacobi

se înlocuieşte cu următoarea secvenţă de comenzi:

⎢ * adaptare soluţie metoda Gauss-Seidel: ⎢ pentru i = 1,n execută ⎢ ⎡ atribuie sum ← 0 ⎢ ⎢ pentru j = 1,i-1 execută ⎢ ⎢ ⎡ atribuie sum ← sum + a(i,j) * xn(j) ⎢ ⎢ ⎣ ⎢ ⎢ pentru j = i+1,n execută ⎢ ⎢ ⎡ atribuie sum ← sum + a(i,j) * xv(j) ⎢ ⎢ ⎣ ⎢ ⎢ atribuie xn(i) ← ( b(i) - sum ) / a(i,i) ⎢ ⎣ ⎢ * sfârşit adaptare soluţie metoda Gauss-Seidel

Observaţii:

1.) Descrierea anterioară a programului foloseşte aceleaşi funcţii MATLAB ca la exerciţiul E2.2. Singura deosebire este reprezentată de înlocuirea funcţiei diag cu funcţia tril (extragerea parţii inferior triunghiulare a unei matrice, inclusiv diagonala, restul de elemente fiind completate cu zerouri).

2.) Programul se va rula pentru datele de intrare menţionate la exerciţiul E2.2, rezulatele analizându-se şi comparându-se cele două metode iterative, Jacobi şi Gauss-Seidel.

2.9 Probleme propuse

P2.1. Se consideră sistemul de ecuaţii algebrice liniare:


ℜ∈α⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

α=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡⋅

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−

−,7

5

xxx

11141142

321

3

2

1

.

Folosind transformări elementare directe (fără pivotare), să se determine pentru ce valoare a parametrului α sistemul (a) nu are soluţii? (b) are o infinitate de soluţii?

P2.2. Să se realizeze programul pentru rezolvarea unui sistem de ecuaţii algebrice liniare:

pnnn B,A,BXA ×× ℜ∈ℜ∈=⋅ ,

prin metoda directă bazată pe triangularizarea cu pivotare parţială a matricei sistemului.

Exemple numerice:

I. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−

−=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−

=

242424101

700

B,

50112230

02123104

A ;

II. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−−

=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

44332211

B,

6824572346223521

A .

P2.3. Să se realizeze programul pentru inversarea unei matrice pătratice reală, de ordinul n, prin metoda directă bazată pe triangularizarea cu pivotare parţială a matricei. Programul va verifica relaţia: n

1 IAA =⋅ − , unde 1A − este inversa calculată, iar nI este matricea unitate de ordinul n.

Exemple numerice:

I. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−

=

50112230

02123104

A ; II. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

6824572346223521

A .


P2.4. Să se realizeze programul pentru inversarea unei matrice pătratice reală, de ordinul n, inferior triunghiulară, exploatând structura particulară a acesteia. Se va verifica rezultatul obţinut efectând produsele:

1AA −⋅ , AA 1 ⋅− .

Exemple numerice:

I. ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

171012006

A ; II. ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−=

812001006

A .

P2.5. Să se realizeze programul pentru calculul determinantului unei matrice pătratice reală, de ordin n, prin metoda directă bazată pe triangularizarea cu pivotare totală a matricei.

Exemple numerice:

I. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

6543876522224321

A ; II. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−

−

=

3221210340223241

A .

P2.6. Să se realizeze programul pentru rezolvarea, în corpul numerelor complexe, a sistemului de ecuaţii:

1nnn b,A,bzA ×× ∈∈=⋅ CC ,

prin transformarea acestuia într-un sistem de ecuaţii echivalent:

n2m,d,C,dxC 1mmm ⋅=ℜ∈ℜ∈=⋅ ×× .

Sistemul transformat se rezolvă utilizând programul realizat pentru rezolvarea problemei 8P2.2. Exemple numerice:

I. 1i,i24

i2b,

i12i2i1

A 2 −=⎥⎦

⎤⎢⎣

⎡−

=⎥⎦

⎤⎢⎣

⎡−

+= ;

II. 1i,i42

i64b,

i21i2i32i23

A 2 −=⎥⎦

⎤⎢⎣

⎡+−+

=⎥⎦

⎤⎢⎣

⎡−+−−−

= . 8


P2.7. Să se realizeze programul pentru aducerea unei matrice nm,A nm >ℜ∈ × , la forma superior triunghiulară prin eliminare gaussiană cu pivotare parţială (de linii). Exemple numerice:

I.

⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

=

10594837261

A ; II.

⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

=

151051494138312721161

A .

P2.8. Să se realizeze programul pentru estimarea numărului de condiţie al unei matrice pătratice reală, de ordin n, folosind relaţia de calcul:

11

11 ||A||||A||)A(k −⋅= .

Inversa matricei A se calculează folosind programul realizat pentru rezolvarea problemei P2.3. Norma 1 a matricelor implicate se va calcula conform definiţiei (Anexa A), prin program.

Exemple numerice:

I. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−−

=

3243210220241321

A ; II. ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

7/16/15/14/16/15/14/13/15/14/13/12/14/13/12/11

A .

P2.9. Se consideră sistemul de ecuaţii algebrice liniare:

ℜ∈⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡==⋅ a,

a42432234

A,bxA .

Care este cea mai mică valoare a parametrului a pentru care metoda iterativă Jacobi va fi convergentă?

P2.10. Să se realizeze programul pentru rezolvarea unui sistem determinat de ecuaţii algebrice liniare, de ordinul n, prin metoda iterativă a suprarelaxării succesive de tip Gauss-Seidel. Considerând ω drept parametru de accelerare a convergenţei metodei iterative, ]2,1[∈ω , programul va determina, crescând succesiv parametrul ω cu pasul 1.0 , valoarea optimă optimω pentru care numărul de iteraţii este minim.


Exemplu numeric:

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡=

695

b,311122112

A .

Curs_04

Documents

Transcript of Curs_04