MetodeNumerice_cap2

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:

. (2.)

Problema de calcul este determinarea unei soluţii 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 , aceasta se numeşte inversabilă (nesingulară) dacă există o matrice astfel încât să fie îndeplinită

relaţia: .

În relaţia anterioară, 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 este inversabilă, atunci oricare ar fi vectorul

există şi este unică soluţia x a ecuaţiei (2.).

Soluţia se poate scrie sub forma:

. (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.). Un exemplu care să ilustreze aceasta este următorul: se consideră ecuaţia , în care şi . Se mai consideră o aritmetică a virgulei mobile cu , şi rotunjire prin tăiere. Atunci soluţia ecuaţiei este:

36 2. Sisteme determinate de ecuaţii algebrice liniare

.

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

2. Nu se recomandă rezolvarea ecuaţiei prin regula Cramer:

,

unde reprezintă determinantul matricei A, iar reprezintă determinanţii

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

Exemplul 2.1:

Problemele de tipul (2.) sunt foarte des întâlnite în practică. Ca exemplu generic, se consideră un proces dinamic cu n mărimi de intrare ( , ..., ) ş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:

. (2.)

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

n an1 ann bn

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

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

valorilor mărimii de ieşire, relaţia (2.) se poate rescrie sub forma . (2.(2.), unde , 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 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ă:

;

b) forma inferior triunghiulară:

.

Procedura de transformare se numeşte triangularizare. Astfel, se furnizează soluţia exactă a sistemului de ecuaţii (2.), î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 . De

aceea, aceste metode se recomandă pentru rezolvarea unor sisteme de ecuaţii algebrice liniare de ordin mai mic decât .

Metodele iterative au la bază construirea unui şir de aproximaţii pentru

soluţia sistemului (2.), care să fie convergent pentru la

soluţia adevărată, :

.


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

,

sau, altfel spus, 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 , asemenea metode se folosesc pentru sisteme de ordin mare şi anume .

2.2 Rezolvarea sistemelor prin triangularizare directă

2.2.1 Principiul metodei

Se consideră problema (2.) şi matricea sistemului: .

Definiţie:

Matricele , …, se numesc submatrice

principale ale lui A sau minori principali directori.

Teoremă:

Dacă matricea are toate submatricele principale inversabile (nesingulare), atunci există matricele astfel încât:

(2.)

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.) se numeşte factorizare L-D-U a matricei A;2. uzuale sunt factorizările: ;3. demonstraţia teoremei enunţate este constructivă, constituind însuşi

algoritmul de descompunere L-U a matricei A;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, , atunci rezolvarea sistemului (2.) implică două subetape:a. rezolvarea sistemului , etapă numită şi substituţie înainte,

obţinând soluţia intermediară . Determinarea componentelor vectorului are loc din aproape în aproape: se începe cu

(prima ecuaţie), se înlocuieşte în a doua ecuaţie determinând pe

şi aşa mai departe.b. rezolvarea sistemului , în care necunoscuta este , etapă

numită şi substituţie inversă. În acest caz, determinarea componentelor vectorului 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: , . De aici

rezultă că , ceea ce conduce la:

,


adică o matrice inferior triunghiulară unitate ( ) este identică cu o matrice

superior triunghiulară ( ). Acest lucru este posibil numai dacă ambele

matrice sunt diagonale şi au diagonala principală unitară, adică: şi ,

ceea ce implică , .

Procedura de triangularizare directă necesită un număr de operaţii în virgulă mobilă de ordinul lui . 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 , 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ă ( ) şi pozitiv definită (

şi ), atunci A se

descompune sub forma . Aceasta se numeşte descompunerea

Cholesky. În acest caz, algoritmul necesită mai puţine operaţii în virgulă mobilă şi anume , 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

pentru execută * determinare matrice kM astfel încât matricea să

aibă elementele:

şi

atribuie

În final se obţine matricea .

Acest algoritm parcurge 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 ,

anume:


,atunci matricea de transformare se construieşte astfel încât vectorul

să aibă elementele:

.

Se consideră vectorul de forma:

,

elementele numindu-se multiplicatori.

Vectorul se numeşte vector Gauss sau vector de multiplicatori.

Subvectorul care conţine strict numai multiplicatorii se numeşte subvector Gauss:

.

Definiţie:

Matricea se numeşte matrice de transformare elementară de ordin n şi

indice k sau matrice Gauss şi este definită prin:,

în care , elementul egal cu 1 fiind în poziţia

k.Definită astfel, matricea este o matrice inferior triunghiulară unitate,

este nesingulară şi deci admite inversă. Inversa acesteia este de forma:

. (2.)

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

ultimele elemente trebuind a fi zerorizate. Presupunând că şi

alegând , va rezulta:

.

Dacă la etapa k a triangularizării, elementul şi

, atunci se obţine:

,

în care reprezintă notaţia pentru coloana k a matricei . Acest rezultat

evidenţiază faptul că primele k elemente din coloana k a matricei rămân


neschimbate, iar ultimele elemente devin zero. Elementul se

numeşte pivot.

Observaţii:

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

,

în loc de a verifica , 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:

.

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 asupra celorlaltor coloane ale matricei

este următorul. Se consideră un vector de forma .

Aplicând transformarea vectorului , se obţine:

.

Concluzii:

a) Matricea lasă nemodificate primele coloane ale matricei .

Considerând vectorul ca fiind coloana numărul j a matricei :

,

atunci se obţine:

b) Matricea transformă coloana k a matricei , zerorizând liniile .

c) Matricea transformă coloanele ale lui în liniile

. Considerând vectorul ca fiind coloana j a matricei :

,


atunci se obţine:,

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

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

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

. (2.)

Înmulţind la stânga relaţia (2.) cu , se obţine:

. (2.)

Notând produsul cu L, atunci relaţia (2.) devine:

.

Ţinând cont de forma matricelor dată de relaţia (2.),

atunci matricea L este egală cu:

.

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.) este substituţia înainte aplicată sistemului de ecuaţii . Vectorul 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:

.

Exemplul 2.2:

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

,


unde reprezintă soluţia adevărată a sistemului. Rezultă, aşadar, erori mari în

soluţia calculată, . 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 (), pentru aritmetica virgulei mobile folosite. Multiplicatorul

corespunzător este , 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 şi .

Concluzie:

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

,

multiplicatorii 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 , î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 printre elementele din coloana k, pornind de la elementul de pe

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

.


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ă , elementul maxim în modul nu se găseşte pe diagonala

principală, atunci se interschimbă (permută) liniile k şi . Lucrul acesta se

realizează automat cu ajutorul unei matrice de permutare de linii care

multiplică la stânga matricea : . Pentru matricea care rezultă astfel, se

determină apoi matricea de transformare ca şi în cazul traingularizării

simple, obţinând matricea:.

În felul acesta, multiplicatorii calculaţi sunt subunitari în modul , iar algoritmul triangularizării devine stabil numeric.

Matricea se numeşte matrice de transformare elementară stabilizată.

Matricea de permutare se obţine din matricea unitate de ordinul n, ,

schimbând adecvat între ele liniile k şi . Această matrice are următoarele

proprietăţi:.

Teoremă:

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

,

în care U este o matrice superior triunghiulară şi este o matrice inferior

triunghiulară unitate. cu elementele .

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 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 pentru execută

* determinare pivot astfel încât

dacă atunci

* determină

altfel atribuie

* calcul

* determinare matrice kM astfel încât matricea

să îndeplinească condiţiile de la

triangularizarea simplă atribuie

atribuie

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

. (2.)

Folosind faptul că , sau altfel spus , relaţia (2.) poate fi

scrisă sub forma:

(2.)

Matricea produs se notează cu P. Ea este numită matrice

generală de permutare de linii. Produsul din prima paranteză din relaţia (2.) se notează cu . Relaţia (2.) se scrie sub forma: , în care matricea

este o matrice inferior

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

Definiţie:

Matricea 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 acest caz, se poate enunţa următorul rezultat.


Propoziţie:

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

,

în care şi elementele matricei sunt în modul subunitare.

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

Rezolvarea sistemului se realizează în două

etape:a) descompunerea -U a matricei A:

b) rezolvarea propriu-zisă a sistemului care comportă trei subetape, şi anume:b.1.) calculul vectorului ;b.2.) rezolvarea sistemului prin substituţie înainte;b.3.) rezolvarea sistemului prin substituţie înapoi.

Această modalitate de rezolvare se bazează pe următoarele relaţii:, , , , .

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: . Se

obţine, în final:

.

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:

.

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.

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 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:

,

unde U este o matrice superior triunghiulară, iar este o matrice inferior

triunghiulară unitate având elementele , î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:

,

unde matricele 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

pentru execută

* determinare pivot care satisface:

dacă atunci

* determinare (permutarea liniilor şi k)

altfel atribuie

dacă atunci

* determinare (permutarea coloanelor şi k)

altfel atribuie

atribuie

* traingularizare matrice :

atribuie

Tabloul general al transformărilor este:

. (2.)

În relaţia (2.) se notează cu S produsul şi se ţine cont de faptul

că . Atunci relaţia (2.) devine:

(2.)

Matricea se notează cu P, iar cu se notează produsul

. În felul acesta, relaţia (2.)

devine: .

Observaţii:


1. Matricea de permutare de linii, , se obţine din matricea unitate

schimbând adecvat între ele liniile k şi . Deoarece se schimbă linii între

ele, matricea se aplică la stânga matricei . Matricea are

proprietăţile:.

2. Matricea de permutare de coloane, , se obţine din matricea unitate

schimbând adecvat între ele coloanele k şi . Deoarece se schimbă coloane

între ele, matricea se aplică la dreapta matricei . Matricea are

proprietăţile:.

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

, obţinându-se:

. (2.)

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

;b) rezolvarea propriu-zisă a sistemului (2.), cu următoarele subetape:

b.1.) calculul vectorului ;b.2.) rezolvarea sistemului prin substituţie înainte;b.3.) rezolvarea sistemului prin substituţie înapoi;b.4.) determinarea soluţiei: .

Observaţii:

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

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

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, , iar elementele care se

transformă devin:.

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:

,

şi se spune că matricele A, ca şi , 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:

,

atunci descompunerea -U cu pivotare totală este:

.

2.5 Aplicaţii ale descompunerilor L-U

2.5.1 Calculul determinantului

Considerând o matrice 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:

, (2.)

unde L este o matrice inferior triunghiulară unitate, iar matricea este o matrice superior triunghiulară. Aplicând funcţia

determinant det(.) relaţiei (2.) se obţine:

. (2.)

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.) devine:


.

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


, (2.)

în care P este matricea generală de permutare de linii (matrice nesingulară),

este o matrice inferior triunghiulară unitate, iar este o matrice

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

. (2.)

Aplicând funcţia determinant det(.) relaţiei (2.), se obţine:.

Ţinând cont de următoarele: , ;

, deoarece dacă este realizată o

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

;

,

atunci relaţia (2.) devine:

,

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ă


, (2.)

î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 este o matrice

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


. (2.)

Dar şi . Ca urmare, relaţia (2.)

devine:

. (2.)

Ţinând cont de următoarele: , deoarece dacă este realizată o

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

, deoarece dacă este realizată

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

;

,

atunci relaţia (2.) devine:

,

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:

. (2.)

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

. (2.)

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

. (2.)


Notând produsul cu Y, relaţia (2.) devine: sau, folosind relaţiile (2.), se poate scrie:

. (2.)

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

, unde operaţii în virgulă mobilă sunt necesare descompunerii L-U, iar 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.)) a matricei A a sistemului, atunci rezolvarea sistemului de ecuaţii matriciale (2.) parcurge următoarele etape:a) descompunerea -U a matricei A (relaţia (2.));b) calculul coloanelor matricei X; pentru se execută:

b.1.) calculul vectorului ;

b.2.) rezolvarea sistemului prin substituţie înainte;

b.3.) rezolvarea sistemului prin substituţie inversă.

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

b.1.) calculul vectorului ;

b.2.) rezolvarea sistemului prin substituţie înainte;

b.3.) rezolvarea sistemului prin substituţie inversă;

b.4.) ordonare soluţie calculată: .

2.5.3 Calculul inversei unei matrice

Fie o matrice nesingulară. Se doreşte aflarea inversei acesteia, notată cu . Acest tip de problemă se încadrează în problematica rezolvării

ecuaţiilor matriciale de tipul (2.), considerând .

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.). În final, matricea inversă este egală cu

matricea ale cărei coloane sunt vectorii rezultaţi la faza a doua menţionată şi

anume: . 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.));2. , , unde 1 apare în poziţia k;

,

triangularizare cu pivotare parţială


,

triangularizare cu pivotare totală


, .

Această procedură de determinare a inversei unei matrice necesită un număr de operaţii în virgulă mobilă de ordinul lui , î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:

, 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:

. (2.)


Î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 . Astfel, se pot rescrie matricele implicate în (2.) sub forma următoare, unde :

,

, (2.)

.

Înlocuind relaţiile (2.) în (2.) se obţine:

.

Efectuând calculele se obţine:

. (2.)

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

,

,

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

. (2.)

Notând cu C matricea de blocuri , cu vectorul şi cu

vectorul , relaţia (2.) se reduce la:

. (2.)

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


,

iar soluţia sistemului complex este:

.

2.6 Metode iterative

2.6.1 Principiul şi convergenţa metodelor iterative

Fie sistemul de ecuaţii algebrice liniare:

, (2.)

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

soluţiei, , convergent la soluţia adevărată:

.

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

. (2.)

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

(2.)

în care estimaţia iniţială este dată (cunoscută). În particular, .

Din relaţia (2.) rezultă:

. (2.)

Se face următoarea notaţie:

.

Matricea G are valori proprii în general complexe, , care

formează mulţimea numită spectrul matricei G.

Definiţie:

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


.

Î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.), să fie convergent către soluţia adevărată a sistemului de ecuaţii (2.) este ca matricea 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: , . Eroarea la pasul [ ] este:

. (2.)

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

. (2.)

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

.

Cum în general, atunci condiţia referitoare la limita şirului de

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

. 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.) 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ă atunci .

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

.

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.):

,

î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:

. (2.)

Din relaţia (2.), folosind relaţia (2.), rezultă:

. (2.)

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

. (2.)

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

. (2.)

Matricea G corespunzătoare metodei Jacobi este:,

în care .

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


. (2.)

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

. (2.)

O matrice care satisface relaţia (2.) 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.).

Observaţie:

Condiţia (2.) înseamnă că , .

Revenind la relaţia (2.), se observă că 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 în care:

.Relaţia (2.) se poate rescrie, în acest caz, sub forma:

. (2.)

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

,

din care rezultă, dacă :

.

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 şi raza spectrală subunitară a matricei

există relaţia:


.

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.) se poate scrie sub forma:

.

În acestă relaţie, în membrul drept se adună şi se scade cantitatea .

Ţinând cont de faptul că rezultă:

Se face următoarea notaţie:,

în care se numeşte reziduu corespunzător iteraţiei [k].

În acest caz, se poate scrie:

. (2.)

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

. (2.)

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

Relaţia (2.) se poate rescrie sub forma:

(2.)


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

. (2.)

Se demonstrează că şirul de aproximaţii obţinut cu ajutorul relaţiei de recurenţă (2.) este convergent pentru . Dacă metoda se numeşte a subrelaxărilor succesive, iar dacă ea se numeşte a suprarelaxărilor succesive. Dacă , 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

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

valorile:

.

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

.

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

în care şi .

Fig. 2.3 Metoda suprarelaxării succesive


***

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

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

.

Notând şi , atunci următorul algoritm descrie procedura

generală aplicată în cazul unei metode iterative:

citeşte

atribuie

atribuie atribuie cât timp execută

atribuie

atribuie * calculează

* calculează

În final, se obţin următoarele: numărul de iteraţii (iter) şi vectorul care

aproximează soluţia sistemului cu precizia impusă . În descrierea anterioară,

reprezintă estimaţia iniţială a soluţiei, iar 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:

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 , se defineşte numărul său de condiţie în

raport cu norma vectorială , notat , ca fiind:

,

,

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

Cum în general , rezultă că şi . Se poate demonstra că

. Astfel, se poate defini numărul de condiţie în raport cu

operaţia de inversare matricială ca fiind:.

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

astfel încât să aibă loc relaţia:.

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

.

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ă (matrice unitate de ordinul n), atunci ;

dacă (matrice elementară de permutare de linii), atunci ;

dacă (matrice elementară de permutare de coloane), atunci

;

dacă (matrice diagonală), atunci

.


Î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:.

În funcţie de valorile pe care le poate lua acesta, matricea A poate fi caracterizată după cum urmează, unde este notaţia pentru epsilonul-maşină:

a) matrice foarte bine condiţionată, dacă ;

b) matrice bine condiţionată, dacă , dar apropiat de valoarea 1:

;

c) matrice prost condiţionată (aproape singulară) dacă ;

d) matrice singulară (foarte prost condiţionată) dacă .

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:

.

Atunci în simplă precizie, deşi , 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:


cu

.

Atunci şi . Pentru n având valori mari, matricea

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:. (2.)

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 . Se disting următoarele două subcazuri, I.1 şi I.2, prezentate în continuare.

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

.Dacă este soluţia exactă a acestei noi probleme, atunci se demonstrează că eroarea relativă satisface la relaţia:

.

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

Exemplul 2.7:


Fie sistemul (2.) cu următoarele matrice:

.

Soluţia exactă a sistemului este . Soluţia exactă a sistemului

este .

Realizeazând următoarele notaţii: , , rezultă pentru acest exemplu:

,

deci matricea A este prost condiţionată.

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

.

De fapt, se rezolvă exact problema:

. (2.)

Dacă x~ este soluţia exactă a ecuaţiei (2.), atunci se demonstrează că eroarea

relativă satisface la relaţia:

.

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

ori.

II. Datele iniţiale ale problemei (2.), A şi , 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: , în cazul triangularizării simple; , în cazul triangularizării cu pivotare parţială; , î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 ( ), U, P, (S). Dacă se efectuează operaţia inversă, se obţine:

, corespunzător triangularizării simple; , corespunzător triangularizării cu pivotare parţială; , corespunzător triangularizării cu pivotare totală,


unde este matricea iniţială plus o matrice de eroare: .

Î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:

, iar dacă ordinul n este mare, atunci: . Aşadar,

norma matricei de eroare se poate apropia de cea a matricei A a sistemului.

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

.

De regulă se foloseşte norma 1. În general, se demonstrează că:.

Dacă , 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:

.

O matrice foarte bine condiţionată este caracterizată de . Dacă matricea este bine condiţionată, atunci , iar dacă ea este prost condiţionată, atunci . 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 este: .Se demonstrează că, în această situaţie, eroarea relativă în soluţia calculată satisface la relaţia:

.

Î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 cu:

.

Soluţia exactă a acestui sistem este .

În urma calculelor rezultă: , .

Aceasta se datorează faptului că:

,

deci matricea A este prost condiţionată.

O concluzie generală care se poate desprinde este următoarea: buna condiţionare a problemei de calcul ( ), împreună cu

stabilitatea numerică a algoritmului de triangularizare ( ),

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 ).

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

triangularizarea matricei A; determinarea soluţiei folosind rezultatele triangularizării matricei A: calculul reziduului asociat ; rezolvarea sistemului , folosind rezultatele procedurii de

triangularizare; dacă , atunci .

Numărul de cifre exacte ale soluţiei calculate este:

,

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 folosind descompunerea L-U a matricei A şi

determinarea soluţiei ;2. calcul reziduu asociat: ;3. rezolvare sistem folosind rezultatele descompunerii L-U a

matricei A şi determinarea lui ;4. rafinarea soluţiei: ;

5. dacă atunci şi reluare de la pasul 2, altfel .


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:

.

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

este mare şi, ca urmare, valoarea 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, .

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 \ bbatribuie 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_bsatribuie nr_bs norm (r_bs,2)

atribuie r_inv bb – aa * x_invatribuie 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>

Linia 2:[4 5 6] <RETURN>

Linia 3:[7 8 9] <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; ;

II. n = 3; ;

(acelaşi exemplu, cu scalare pe linii)

III. n = 3; ;

IV. n = 4; .

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, .

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 - aatribuie g inv (nn) * patribuie 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 1atribuie 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; ; ;

;

II. n = 3; ; ;

; .

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.

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:

.

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:

,

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

Exemple numerice:

I. ;

II. .

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: , unde

este inversa calculată, iar este matricea unitate de ordinul n.

Exemple numerice:


I. ; II. .

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:

, .

Exemple numerice:

I. ; II. .

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. ; II. .

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

,

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

.

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

I. ;


II. .

P2.7. Să se realizeze programul pentru aducerea unei matrice , la forma superior triunghiulară prin eliminare

gaussiană cu pivotare parţială (de linii).Exemple numerice:

I. ; II. .

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:

.

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. ; II. .

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

.

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, ,


programul va determina, crescând succesiv parametrul cu pasul , valoarea optimă pentru care numărul de iteraţii este minim.

Exemplu numeric:

.

MetodeNumerice_cap2

Documents

Transcript of MetodeNumerice_cap2