Rezolvarea Sistemelor de Ecuatii I

7/25/2019 Rezolvarea Sistemelor de Ecuatii I

http://slidepdf.com/reader/full/rezolvarea-sistemelor-de-ecuatii-i 1/17

METODE MUMERICE – LUCRĂRI DE LABORATOR #06

77

LUCRAREA #06

REZOLVAREA SISTEMELOR DE ECUAŢI I I

REZUMATUL ŞI SCOPUL LUCRĂRII

Rezolvarea sistemelor de ecuaţii liniare se poate face prin două tipuri de metode: metode directe: constau în transformarea sistemului într–un sistem

particular (diagonal, triunghiular) echivalent, care se rezolvă uşor prinmetode elementare. Cele mai cunoscute metode directe sunt: metoda Gauss,

metoda Cholesky (utilizată pentru sistemele în care matricea A estesimetrică şi pozitiv definită) şi metoda Householder. Pentru sisteme cu unnumăr de ecuaţii mai mare de 100, metodele directe devin inutilizabiledatorită acumulării erorilor de rotunjire care alterează soluţia.

metode indirecte sau iterative: pornind de la o aproximaţie iniţială x0 şi dela o relaţie de iteraţie, se construieşte un şir de aproximaţii xk, care converge

în anumite condiţii către soluţia exactă a sistemului. Procesul iterativ seopreşte atunci când aproximaţia de ordinul k se încadrează între limiteleunei precizii stabilite iniţial. Cele mai cunoscute sunt metoda Jacobi,metoda Gauss-Seidel, metoda relaxării.

Lucrarea are drept scop prezentarea metodelor directe de rezolvare a sistemelorde ecuaţii algebrice liniare: metoda eliminării Gauss şi metoda pentru rezolvareasistemelor tridiagonale. Se vor elabora funcţiile specifice în MatLab, iar î nsuşireametodelor se va face pe exemple concrete de sisteme de ecuaţii.

CUPRINSUL LUCRĂRII

1. Rezolvarea sistemelor de ecuaţii liniare

2. Condiţionarea sistemelor. Norma unei matrice 3. Metoda eliminării lui Gauss 4. Rezolvarea sistemelor tridiagonale5. Aplicaţii de laborator



#06 Rezolvarea sistemelor de ecuaţii I

78

1. Rezolvarea sistemelor de ecuaţii liniare

Un sistem de ecuaţii liniare are următoarea formă:

11 1 12 2 1 1

21 1 22 2 2 2

1 1 2 2

1 1 2 2

n n

n n

k k kn n k

n n nn n n

a x a x a x b

a x a x a x b

a x a x a x b

a x a x a x b

(6.1)

Acest sistem se mai poate scrie mai simplu astfel: A x b , unde A este matricea

coeficienţilor, x vectorul necunoscutelor şi b vectorul termenilor liberi sau:

11 12 1 1 1

21 22 2 2 2

` 2

1 2

n

n

k k kn k k

n n nn n n

a a a x b

a a a x b

a a a x b

a a a x b

(6.2)

Dacă matricea A este nesingulară, adică determinantul său este diferit de 0, atunciexistă o matrice A-1, astfel încât 1

x A b .

Dacă însă det( A) = 0 atunci sistemul de ecuaţii este degenerat şi poate avea oinfinitate de soluţii sau nici una.

Deci, cea mai simplă metodă de rezolvare a sistemelor de ecuaţii liniare esterezolvarea prin împărţire la stânga sau folosirea matricei inverse, condiţia fiind impusădeterminantului matricei.

Exemplul 6.1. de rezolvare în MatLab pentru sistemul:2 2 1

2 2 2 210 5 5

x y z

x y y x y z

(6.3)

A=[ 1 2 3; 2 2 2; 10 5 1]b=[ 1 2 5]d=det ( A)i f d~=0

x1=A\ b' %i mpart i r e l a st angax2=i nv(A) *b' %i nver sa mat r i cei

el sedi sp( ' Mat r i ce si ngul ar a. . . Si st emul nu ar e sol ut i e' )

end

d =- 2

x1 =




79

- 4. 000010. 0000- 5. 0000

x2 =- 4. 000010. 0000

- 5. 0000

Relaţia de exprimare a soluţiei prezentată mai sus nu este o metodă practică derezolvare a sistemelor de ecuaţii liniare, deoarece inversarea unei matrice de dimensiunimari este dificilă.

De aceea, se utilizează diferite metode numerice de rezolvare a acestor sisteme,metode ce se clasifică în două categorii:

metode directe: constau în transformarea sistemului într–un sistemparticular (diagonal, triunghiular) echivalent, care se rezolvă uşor prinmetode elementare. Cele mai cunoscute metode directe sunt: metoda Gauss,

metoda Cholesky (utilizată pentru sistemele în care matricea A estesimetrică şi pozitiv definită) şi metoda Householder. Pentru sisteme cu unnumăr de ecuaţii mai mare de 100, metodele directe devin inutilizabiledatorită acumulării erorilor de rotunjire care alterează soluţia.

metode indirecte sau iterative: pornind de la o aproximaţie iniţială x0 şi dela o relaţie de iteraţie, se construieşte un şir de aproximaţii xk, care converge

în anumite condiţii către soluţia exactă a sistemului. Procesul iterativ seopreşte atunci când aproximaţia de ordinul k se încadrează între limiteleunei precizii stabilite iniţial. Cele mai cunoscute sunt metoda Jacobi,metoda Gauss-Seidel, metoda relaxării.

În cele ce urmează se vor aprofunda şi se vor scrie algoritmi în MatLab pentrumetoda eliminării Gauss şi pentru metoda eliminării pentru sistemele tridiagonale.

2. Condiţionarea sistemelor. Norma unei matrice

Înainte de a trece efectiv la metodele de rezolvare a sistemelor de ecuaţii liniare,este necesar a se trece în revistă câteva noţiuni despre norme şi condiţionarea matricelor.

Importanţa cunoaşterii condiţionării unei matrice decurge din faptul că rezultateleobţinute prin calcul numeric corespund întotdeauna unei probleme “perturbate”, caurmare a efectelor erorilor de rotunjire.

Numărul de condiţionare al unei matrice indică sensibilitatea soluţiei unui sistemde ecuaţii liniare şi dă indicaţii despre precizia rezultatelor la inversarea matricelor şi larezolvarea sistemelor de ecuaţii liniare. O matrice bine condiţionată este o matrice relativinsensibilă la mici perturbaţii.

Numărul de condiţionare al unei matrice se defineşte astfel:

1( )cond A A A (6.4)

unde prin s-a notat norma matricei.




80

Pentru matrice sunt definite următoarele tipuri de norme: norma 2, returnează cea mai mare valoare singulară a lui A şi este definită de

relaţia:

2

* max , A A A A (6.5)

norma 1, returnează cea mai mare sumă a elementelor de pe coloană şi estedefinită de relaţia:

1 11

maxn

iji n

j

A x

(6.6)

norma infinită, returnează cea mai mare sumă a elementelor de pe linie, curelaţia:

11

maxn

jii n

j

A x

(6.7)

norma Frobenius, calculată cu relaţia: 2

, 1

n

ijF i j

A x

(6.8)

În MatLab, aceste norme se pot calcula direct prin apelarea următoarelor funcţii, înordinea enumerată mai sus:

norm( x, 2) – nor ma 2norm( x, 1) – nor ma 1nor m( x, i nf ) – norma infinită nor m( x, ’ f r o’ ) – nor ma f r obeni us

Numărul de condiţionare al matricelor are următoarele proprietăţi: ncond(I ) 1

-1cond(A)=cond A

cond( A)=cond A , pentru orice 0

1

2

n

cond A

, 1 2 ... 0

n sunt valorile proprii la matricei

T B A A

Dacă A este simetrică, atunci 2

max

min

i

i

cond A

, unde 1,..., n

sunt valorile

proprii ale matricei A.

MatLab are definite funcţii pentru calcularea numărului de condiţionare al uneimatrice cu următoarele sintaxe:

c = cond( A)

– numărul de condiţionare calculat ca raportul dintre cea mai mare şi cea mai micăvaloare singulara




81

c = r cond( A)

- dacă o matrice este bine condiţionată, valoarea funcţiei este aproximativ 1.0, dacănu valoarea este aproximativ 0.0.

Exemplul 6.2. Se studiază condiţionarea sistemelor: [B.1]1 2

1 2

6 6.917 6.543

1.152 1.095

x x

x x

şi sistemul perturbat 1 2

1 2

6 6.911 6.543

1.152 1.095

x x

x x

Soluţiile sistemelor sunt (rulare program ex6_2):

x1 =7. 3158

- 5. 4000x2 =

- 30. 009027. 0000

Se observă slaba condiţionare a sistemului, prin diferenţa mare a soluţiilor.

3. Metoda eliminării lui Gauss

În cazul metodei eliminării lui Gauss, necunoscutele se elimină pe rând, sistemultransformându-se într-un sistem particular echivalent superior triunghiular şi soluţiile secalculează prin retrosubstituţie sau substituţie regresivă; din ultima ecuaţie se determinăxn, care ulterior se foloseşte la găsirea necunoscutei xn-1 şi aşa mai departe până la x1.

La pasul k, k=1..n-1, se împarte ecuaţia numărul k cu factorul akk (pivot) şi se înmulţeşte cu aki, pentru i=2..n şi prin adunare la linia i se elimină necunoscuta x i din toateecuaţiile de sub cea cu numărul k. Se obţine un sistem de forma:

' ' ' ' '

11 1 12 2 1 1 1

2 2 2 2

22 2 2 2 2

... ...

... ...

..........................................................

...

.....................................

k k n n

k k n n

k k k

kk k kn n k

k

ik

a x a x a x a x b

a x a x a x b

a x a x b

a

...

......................................

...

k k

k in n i

k k k

nk k nn n n

x a x b

a x a x b

(6.9)

Coeficienţii necunoscutelor se modifică astfel: Pentru linia pivot k avem:

1 1

1 1

1

/ , 1../

k

kk

k k k

kj kj kk k k k

k k kk

a

a a a j k nb b a

(6.10)




82

Pentru liniile de sub linia pivot avem:

1 1

1 1

0

, 1.. , 1..

k

ik

k k k k

ij ij ik kj

k k k k

i i ik k

a

a a a a j k n i k n

b b a b

(6.11)

Acestea sunt de fapt formulele de iteraţie pentru eliminarea gaussiană.

Repetând paşii pentru toate ecuaţiile, se obţine în final sistemul superiortriunghiular, a cărui rezolvare este foarte facilă şi care are următoarea formă:

' ' ' ' '

11 1 12 2 1 1 1

2 2 2 2

22 2 2 2 2

... ...

... ...

.............................................................

......................................

k k n n

k k n n

k k k

kk k kn n k

nn

a x a x a x a x b

a x a x a x b

a x a x b

a

n n

n n x b

(6.12)

Din ultima ecuaţie rezultă:

n

n x n

nn

b x

a (6.13)

Cu această valoare se poate stabili o relaţie de iteraţie pentru celelaltenecunoscute:

1

1nk k

k k kj j k j k kk

x b a xa

, pentru k=n-1..1 (6.14)

Această metodă, fără a se efectua alte operaţii asupra matricei sistemului, senumeşte metoda de eliminare a lui Gauss fără pivotare.

Există cazuri pentru care pe diagonala principală a matricei există un element egalcu zero. În acest caz, împărţirea ecuaţiei la pivotul akk nu este posibilă şi ca atare metodanu se poate aplica, chiar dacă sistemul admite soluţii. În acest caz se recurge la aşa zisapivotare, care poate fi parţială sau totală.

Pivotarea parţială constă în permutarea liniilor matricei între ele astfel încât pediagonala principală să avem elementele cele mai mari, în valoare absolută, acest lucru şipentru a minimiza erorile ce pot apărea în calcule.

Acest lucru presupune ca la pasul k, să se găsească elementul cel mai mare înmodul, care se află pe coloana k, sub elementul akk. Fie acest element aik, pentru k<i<n. Înacest caz se permută element cu element linia k şi i a matricei. Sistemul nu se schimbă, se

schimbă doar ordinea ecuaţiilor în sistem, iar soluţia nu este afectată. Se face observaţia, că şi în acest caz, este posibil ca pe o linie a matricei pivotată săexiste un element akk egal cu zero. În acest caz metoda lui Gauss nu se poate aplica.




83

Pentru a elimina şi acest inconvenient, se aplică pivotarea totală, care constă încăutarea pivotului de modul maxim printre toţi coeficienţii aflaţi sub linia k şi la dreaptacoloanei k şi aducerea acestui element pe poziţia (k;k) prin permutări de linii şi de coloane(se schimbă atât ordinea liniilor cât şi cea a coloanelor). La metoda cu pivotare totală,datorită permutării coloanelor, se modifică şi ordinea necunoscutelor în sistem. De aceea,se impune memorarea permutărilor de coloane făcute pe parcursul algoritmului, pentru ase putea reconstitui la sfârşit ordinea iniţială a necunoscutelor.

Pentru implementarea în MatLab a metodei lui Gauss, o soluţie ar fi următoarea: După introducerea matricei sistemului şi a vectorului termenilor liberi, se poate

face un test simplu pentru a se verifica dacă vreun element al matricei de pe poziţia (k,k)este egal cu zero. Dacă acest lucru nu se întâmplă se poate aplica eliminarea lui Gauss fărăpivotare, dar în cazul unui sistem rău condiţionat este posibil ca erorile să fie destul demari.

În caz că se găseşte un element egal cu zero, se poate aplica pivotarea parţială sirezolvarea sistemului echivalent.

Dacă şi după pivotare un element de pe poziţia (k,k) rămâne 0, atunci ultimasoluţie este aplicarea pivotării totale.

Următoarea soluţie de implementare în MatLab, consta în scrierea a 3 funcţii, carevor face pivotarea parţială, eliminarea lui Gauss şi calcularea soluţiei prin retrosubstituţie.S-a ales acest mod, pentru ca funcţiile să poată fi sau nu apelate funcţie de cazurileprezentate mai sus. Mai mult, dacă de la început sistemul este superior triunghiular,atunci nu mai are rost a se face pivotare sau eliminare, ci se poate apela direct funcţia decalcul a soluţiei sistemului.

Să le prezentăm pe rând: Funcţia de pivotare parţială:

Date de intrare:

A – matricea sistemului

b – vectorul termenilor liberi

n – numărul de ecuaţii Date de ieşire:

A – matricea sistemului schimbată prin permutări de linii b – vectorul termenilor liberi schimbat prin permutări de linii

pentru k=1 până la n % numărul de iteraţii

max= kk

a %abs(akk)

i=k %linia pe care se presupune că se găseşte elem. maxim în modul sw=0 %var i abi l a swi t ch cu doua val ori : sw=0 nu s- a gasi t al t el ement maimare decat el ement ul de%pe di agonal a pr i nci pal a; sw=1 s- a gasi t al t el ement mai mare si se vaf ace schi mbar ea de l i ni i

pentru j=k+1:n %parcurge col oana de sub el ement ul A( k, k} - el em. de pedi agonal a pr i nci pal a

dacă max< kk

a

atunci max= kk

a

i=j; sw=1 %s- a gasi t un el ement mai mar e dacă sw=1 %schi mbar e l i ni i i nt r e el e i n mat r i cea A si r espect i v vect or ul b

aux=b( k) ;b( k) =b( i ) ;

b( i ) =aux;pent r u j =k: naux=A( k, j ) ;




84

A(k, j )=A( i , j ) ;A( i , j ) =aux;

Afişează A şi b

Pentru permutarea a două linii între ele se foloseşte un ciclu repetitiv, iar pentru

schimbarea efectivă a elementelor corespondente pe de liniile în cauză se foloseşte metodacelor 3 pahare. Această metodă constă în utilizarea unei variabile auxiliare aux pentru a nuse pierde valorile iniţiale ale elementelor: la primul pas se “goleşte” conţinutul primeivariabile în variabila aux, se atribuie primei variabile valoarea celei de-a doua şi la altreilea pas, a doua variabilă ia valoarea aux, care de fapt este prima variabilă.

Ex: fie x şi y cele două variabile:1: aux = x2: x = y3: y = aux.Se observă simetria paşilor metodei.

Codul sursa pentru funcţia MatLab este următorul:

f unct i on [ A, b] =pi vot ar e_par t i al a( A, b)%f uncti e car e r eal i zeaza pi vot ar ea par t i al a l a un si st em l i ni ar de ecuat i i%i n ur ma pi vot ar i i par t i al e, pe di agonal a pr i nci pal a se vor gasi el ement el e cel emai mar e i n modul%ast f el l a i mpar t i r ea pe l i ni a pi vot , er or i l e sunt mi ni me, deci sol ut i asi st emul ui est e mai exact a% Dat e de i nt r ar e:% A - mat r i cea si st emul ui% b - vector ul t er meni l or l i ber i%Date de i esi r e:

% A - mat r i cea si st emul ui schi mbat a pr i n per mut ar i de l i ni i% b - vectorul t er meni l or l i ber i schi mbat pr i n per mut ar i de l i ni i

n=l engt h( b) ;for k=1: n

max=abs( A( k, k) ) ;i =k; %l i ni a pe car e se gasest e el ement ul maxi m i n modulsw=0; %var i abi l a swi t ch cu doua val or i : sw=0 nu s- a gasi t al t el ement mai

mare decat el ement ul de%pe di agonal a pr i nci pal a; sw=1 s- a gasi t al t el ement mai mare si se va

f ace schi mbar ea de l i ni ifor j =k+1: n %par curge col oana de sub el ement ul A( k, k} - el em. de pe di agonal a

pri nci pal a

i f max<abs( A( j , k) )max=abs( A( j , k) ) ;i =j ;sw=1; %s- a gasi t un el ement mai mar e

endendi f sw==1 %schi mbar e l i ni i i nt r e el e i n mat r i cea A si r espect i v vect or ul b

aux=b( k) ;b( k) =b( i ) ;b( i ) =aux;for j =k: n

aux=A( k, j ) ;A( k, j )=A( i , j ) ;

A( i , j ) =aux;endend




85

%af i sare matr i ce si si st eme i nt ermedi areAbpause

end

Funcţia de eliminare gaussiană s-a scris folosind următorul algoritm: elementele de pe linia pivotului se împart la pivot:

1 , 1..k

k kiki

kk

aa i k n

a

; 1 1/k k k

k k kk b b a (6.15)

elementele de pe poziţie (k,k) se vor înlocui cu 1 (prin împărţire cu

elementul de pe poziţia (k,k) devin 1) elementele de pe coloana pivotului, cu excepţia pivotului, se înlocuiesc cu 0.

ai k=0, i= k+1…n elementele din submatricea delimitată de ultimele n-k linii se transformă curegula dreptunghiului:

1 1

1 1

, 1.. , 1..

, 1..

k k k k

ij ij ik kj

k k k k

i i ik k

a a a a j k n i k n

b b a b i k n

(6.16)

În urma acestor operaţii, sistemul se reduce la sistemul superior triunghiular. Algoritmul funcţiei este următorul: funcţia se bazează pe relaţiile iterative

prezentate şi operarea doar a elementelor matricei aflate deasupra diagonalei principale(pentru a micşora numărul de operaţii) şi înlocuirea elementelor de sub diagonalaprincipală, exclusiv, cu 0.

Date de intrare:

A – matricea sistemului

b – vectorul termenilor liberi

n – numărul de ecuaţii Date de ieşire:

A – matricea sistemului schimbată prin eliminare gaussiana b - vectorul termenilor liberi schimbat prin eliminare gaussiana

pentru k=1 până la n-1 % etapele eliminariib( k) =b( k)/ A( k, k) pentru j=k+1 pană la n+1 % elementele liniei pivot

A( k, j ) =A( k, j ) / A( k, k) ; A( k, k)=1;pent r u i =k+1 pana l a n %el ement el e l i ni i l or de sub l i ni a pi vot

pentru j=k+1 pana la n+1

A( i , j ) =A( i , j ) - A( k, j ) * A( i , k) ;b( i ) =b( i ) - b( k) *A( i , k)

pentru i=k+1 pana la n

A( i , k) =0; Afişează A

Codul sursă al funcţiei MatLab este următorul:




86

f unct i on [ A, b] =el i mi nar e_gauss( A, b)%f unct i e pent r u el i mi nar ea gaussi ana pent r u un si st em de ecuat i i l i ni ar%Date de i nt r are:% A - mat r i cea si st emul ui% b - vectorul t ermeni l or l i beri%Date de i esi r e:

% A - mat r i cea si st emul ui schi mbat ? pr i n el i mi nar e gaussi ana% b - vect or ul t er meni l or l i ber i schi mbat pr i n el i mi nar e gaussi ana

n=l engt h( b) ;for k=1: n- 1 %et apel e el i mi nar i i

i f A( k, k) ==0di sp( ' Nu se poat e apl i ca metoda Gauss. . . ' )return;

end%el ement el e l i ni ei pi votb( k) =b( k) / A( k, k) ;for i =k+1: n

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

endA( k, k)=1;%el ement el e l i ni i l or de sub l i ni a pi votfor i =k+1: n

for j =k+1: nA( i , j ) =A( i , j ) - A( k, j ) * A( i , k) ;

endb( i )=b( i ) - b(k)*A( i , k) ;

end%coef i ci ent i i necunoscut el or el i mi nat efor i =k+1: n

A( i , k) =0;endAbpause

end

Funcţia pentru calculul soluţiilor prin retrosubstituţie utilizează de asemeneacicluri repetitive, folosind relaţiile prezentate la începutul subcapitolului.

Se prezintă direct codul sursă al funcţiei MatLab:

f unct i on x=sol ut i e( A, b)%cal cul ar e sol ut i i pr i n retr osubsti t ut i e pent r u un si stemde ecuat i i l i ni arsuper i or t r i unghi ul ar

% Dat e de i nt r ar e% A - mat r i cea si st emul ui super i or t r i unghi ul ar% b - vectorul t ermel i l or l i ber i% n - numar ul de ecuat i i% Dat e de i esi r e% x - sol ut i a si stemul ui super i or t r i unghi ul ar

n=l engt h( b) ; i f A( n, n) ~=0

x( n) =b( n) / A( n, n) ;for i =n- 1: - 1: 1

i f A( i , i ) ~=0s=b( i ) ;

for j =n: - 1: i +1s=s- A( i , j ) * x( j ) ;

end




87

x( i )=s/ A( i , i ) ;el se

di sp( ' i mpar t i r e l a zer o' )r et ur n

endend

el sedi sp( ' i mpar t i r e l a zer o' )return

end

Se observă că la fiecare iteraţie se testează dacă elementul akk este egal cu zero,pentru a evita o împărţire imposibilă.

Practic, pentru a calcula soluţiile unui sistem de ecuaţii folosind această metodă,respectiv aceste funcţii definite mai sus, se parcurg următoarele etape:

1. se definesc elementele sistemului, respectiv matricea A, vectorul termenilorliberi şi numărul de ecuaţii;

2. se calculează determinantul matricei A; dacă acesta este egal cu 0, sistemul nuare soluţie;

3. se stabileşte dacă este sau nu nevoie de pivotare (adică, dacă avem pe diagonalaprincipală elemente egale cu 0);

4. se apelează pe rând funcţiile pivotare_partiala (dacă este cazul) şi eliminare_gauss,transformându-se sistemul într-un sistem superior triunghiular;

5. se apelează funcţia soluţie pentru a calcula soluţia sistemului.

Exemplul 6.3. Să se rezolve sistemul următor folosind metoda eliminării lui Gauss,cu şi fără pivotare parţială:

1 2 3

1 2 3

1 2 3

2 29

8 2 3 1

2 5 1

x x x

x x x

x x x

Se definesc elementele sistemului:

» A=[ 1 0. 222 1; 8 2 - 3; 1 2 - 5]» b=[ 2 - 1 1]» n=3

Determinantul matricei A este diferit de 0:

» det ( A)ans =

18. 2140

Nu este nevoie de pivotare. Pentru exemplificare se va efectua această operaţie:

» [ A, b] =pi vot ar e_par t i al a( A, b)A =

8. 0000 2. 0000 - 3. 00001. 0000 0. 2220 1. 00001. 0000 2. 0000 - 5. 0000

b =




88

- 1 2 1A =

8. 0000 2. 0000 - 3. 00001. 0000 2. 0000 - 5. 00001. 0000 0. 2220 1. 0000

b =

- 1 1 2A =8. 0000 2. 0000 - 3. 00001. 0000 2. 0000 - 5. 00001. 0000 0. 2220 1. 0000

b =- 1 1 2

Eliminare Gauss:

» [ A, b] =el i mi nar e_gauss( A, b)A =

1. 0000 0. 2500 - 0. 37500 1. 7500 - 4. 62500 - 0. 0280 1. 3750

b =- 0. 1250 1. 1250 2. 1250

A =1. 0000 0. 2500 - 0. 3750

0 1. 0000 - 2. 64290 0 1. 3010

b =- 0. 1250 0. 6429 2. 1430

A =1. 0000 0. 2500 - 0. 3750

0 1. 0000 - 2. 64290 0 1. 3010

b =- 0. 1250 0. 6429 2. 1430

Soluţia sistemului:

» x=sol ut i e( A, b)x =

- 0. 7563 4. 9962 1. 6472

4. Rezolvarea sistemelor tridiagonale

Un astfel de sistem poate fi scris astfel:

1 1k k k k k k ia x b x c x d (6.17)

unde a1=0 şi cn=0.

În formă matriceală, sistemul poate fi scris astfel:




89

1 1 1 1

2 2 2 2 2

3 3 3

1 1 1 1 1

. .

. . . . . . .

n n n n n

n n n n

b c x d

a b c x d

a b c

a b c x d

a b x d

(6.18)

Se observă acum că matricea sistemului are toate elementele egale cu zero, maipuţin diagonala principală şi cele două diagonale de deasupra şi dedesubtul celeiprincipale, de unde şi numele de sistem tridiagonal.

Un astfel de sistem este un caz particular al unui sistem de ecuaţii liniar. De aceearezolvarea acestuia se face printr-o eliminare gaussiană simplificată, metodă ce se mainumeşte şi metoda Thomas.

În primă fază se elimină necunoscutele care au coeficienţii pe ai, elementele bi şi dimodificându-se corespunzător, iar ci rămân neschimbaţi. Formulele sunt următoarele:

'

'

1 1

'

1 1

'

0, 1: 1

, 1: 1

, 1: 1

, 1: 1

k

k

k k k

k

k

k k k

k

k k

a k n

ab b c k n

b

ad d d k n

b

c c k n

(6.19)

Sistemul va avea o formă bidiagonală, din care din ultima ecuaţie se poate deducexn. Apoi prin substituţie inversă se pot calcula şi celelalte necunoscute, folosindu -senecunoascuta determinată la pasul anterior:

'

'

'

1

', 1:1

nn

n

k k k k

k

d x

b

d c x x k n

b

(6.20)

Se poate scrie astfel algoritmul

Date de intrare:

a - subdiagonala cu n-1 componente

b – diagonala principala cu n componente

c – supra diagonala cu n-1 componente

d – termenii liberi cu n componante

Date de ieşire:

x – soluţia sistemului tridiagonal

1: pentru k=1 până la n-1

Calculează'

1 1k

k k k

k

ab b c

b




90

Calculeaza'

1 1k

k k k

k

ad d d

b

' 0k

a

'

k k c c

2: Calculează

'

'

nn

n

d x

b

Pentru k=n-1 pana la 1 cu pasul -1

Calculeaza

'

1

'

k k k k

k

d c x x

b

3: afişează valorile obţinute

Cu cele expune mai sus se poate scrie o funcţie care să calculeze soluţia unuisistem tridiagonal.

Algoritmul de mai sus este un algoritm fără pivotare, care se poate aplica dacătoate elementele bk sunt diferite de zero, astfel vom avea împărţire la zero. În interiorul funcţiei vom pune toate condiţiile necesare pentru a evita astfel de

neplăceri. Vom lucra cu vectori care conţin elementele de pe subdiagonala, diagonala,

supradiagonala şi termenii liberi pentru lejeritate şi vom forma pe parcurs matriceasistemului.

f unct i on x=t r i sys(a, b, c, d)% el i mi nar e gaussi ana f ar a pi vot , si st em t r i di agonal% Dat e de i nt ar e

% a - subdi agonal a% b - di agonal a% c - supr adi agonal a% d - t ermeni i l i ber i% Dat e de i esi r e% x - sol ut i a si stemul ui t r i di agonal

% t est si stem t r i di agonali f ( l engt h( d) ~=l engt h( b) ) | ( l engt h( a) ~=l engt h( c) )

di sp( ' Si st emul nu est e t r i di agonal ' )return

endn=l engt h( d) ; x=zer os( n, 1) ;

% f or mare mat r i cea si st emul ui A=di ag( b) +di ag( a, - 1) +di ag( c, 1)% t est condi t i e pent r u el emi nar e i f ~al l ( b)

di sp( ' I mpar t i r e l a zer o. Cel put i n un el ement pe pe di agonal a pr i nci pal a est eegal cu zer o. . . ' )

returnendi f det( A) =0 er r or ( ' Si st emul nu ar e sol ut i e' ) ; return; end%el i mi nar efor k=1: n- 1

m=a( k) / b( k) ; %mul t i pl i cat or ulb( k+1) =b(k+1)- m*c( k) ;d( k+1) =d(k+1)- m*d( k) ;

end% f or mare mat r i ce a si st emul ui dupa el i mi nar e A1=di ag( b) +di ag( c, 1)




91

% t est pent r u cal cul ar e sol ut i i i f ~al l ( b)

di sp( ' I mpar t i r e l a zer o. Cel put i n un el ement pe pe di agonal a pr i nci pal a est eegal cu zero. . . ' )

returnend

%subst i t ut i e i nver sax( n) =d( n) / b( n) ;for k=n- 1: - 1: 1

x( k) =( d( k) - c(k) *x( k+1) ) / b( k) ;end

Evident că, dacă şi după eliminare avem elemente egale cu zero pe diagonalaprincipală, atunci se poate aplica pivotarea.

Exemplul 6.4 Să se rezolve sistemul:

1 2

1 2 3

2 3 4

3 4 5

4 5

2 7

3 3 3 2

2 1

4 6 4 9

5 5 4

x x

x x x

x x x

x x x

x x

Matriceal sistemul se scrie astfel:

1

2

3

4

5

1 2 7

3 3 3 2

2 1 1 1

4 6 4 9

5 5 4

x

x

x

x

x

Nu avem elemente pe diagonala principală egale cu zero, deci putem aplicaeliminarea.

Primul pas definim vectorii care conţin elementele sistemului, respectiv:

» a=[ 3 2 4 5]» b=[ 1 3 1 6 5]» c=[ 2 3 1 4]» d=[ 7 2 1 9 4]

Apelăm funcţia trisys:

» x=t r i sys( a, b, c, d)

A =1 2 0 0 03 3 3 0 00 2 1 1 00 0 4 6 40 0 0 5 5

A1 =1. 0000 2. 0000 0 0 0

0 - 3. 0000 3. 0000 0 0




92

0 0 3. 0000 1. 0000 00 0 0 4. 6667 4. 00000 0 0 0 0. 7143

x =23. 4667

- 8. 2333- 14. 566732. 0333

- 31. 2333

5. Aplicaţii de laborator

Aplicaţia 6.1. Să se studieze sistemul:

1 2 3

1 2 3

1 2 3

8 2 15

10 4 21

50 25 8 124

x x x

x x x

x x x

Cu metoda eliminării lui Gauss se obţine soluţia:

» [ A, b] =pi vot ar e_par t i al a( A, b)A =

50 25 810 4 18 2 1

b =124 21 15

A = 50 25 810 4 18 2 1

b =124 21 15

A =50 25 810 4 18 2 1

b =124 21 15

» [ A, b] =el i mi nare_gauss( A, b)A =1. 0000 0. 5000 0. 1600

0 - 1. 0000 - 0. 60000 - 2. 0000 - 0. 2800

b =2. 4800 - 3. 8000 - 4. 8400

A =1. 0000 0. 5000 0. 1600

0 1. 0000 0. 60000 0 0. 9200

b =2. 4800 3. 8000 2. 7600

A = 1. 0000 0. 5000 0. 16000 1. 0000 0. 60000 0 0. 9200




93

b =2. 4800 3. 8000 2. 7600

» [ x] =sol ut i e( A, b, n)x =

1. 0000 2. 0000 3. 0000

Rezolvarea Sistemelor de Ecuatii I

Documents

Transcript of Rezolvarea Sistemelor de Ecuatii I