METODE NUMERICE - Indrumar de laborator

Mitică TEMNEANU

METODE NUMERICE - Indrumar de laborator -

Iaşi, 2009

Metode Numerice – Indrumar de laborator

PREFAŢĂ Prezenta lucrare se adresează studenţilor şi îşi propune să exemplifice modul în care pot fi puse în practică diversele metode de calcul numeric prezentate în cadrul cursului de Metode Numerice. Este evident că, atunci când o anumită temă nu se poate finaliza prin utilizarea metodelor clasice, este necesară găsirea unei modalităţi de calcul care să permită totuşi obţinerea rezultatelor problemei în cauză. În această situaţie intervine calculul numeric, adică implementarea unei metode numerice pe un sistem de calcul, utilizând un anumit limbaj de programare. Pentru fixarea mai clară a metodei de calcul abordate, au fost readuse în atenţia studenţilor elementele de bază specifice fiecărei metode de calcul, cu evidenţierea paşilor teoretici ce trebuie parcurşi pentru rezolvarea unei probleme printr-o anumită metodă. Programele au fost elaborate în MATLAB, având în vedere faptul că acest mediu de programare este cunoscut studenţilor, el fiind utilizat şi în cadrul lucrărilor practice aferente altor discipline. Fiind destinate înţelegerii metodelor de calcul şi a modului efectiv de realizare a modulelor soft, fiecare dintre aceste programe este conceput a se desfăşura pas cu pas, studentul având posibilitatea de a analiza modul în care a fost efectuat calculul în cadrul unui anumit pas şi de a analiza şi chiar de a anticipa ce presupune parcurgerea pasului următor. Acest lucru este deosebit de util, deoarece impune studentului aprofundarea teoretică a metodei respective de calcul. Având în vedere modul în care au fost concepute aceste programe, studentul are posibilitatea de a aplica metodele de calcul prezentate fie pe unele exemple predefinite, fie pe probleme concrete pe care le are de rezolvat. Programele prezentate sunt utile nu numai studenţilor ci şi tuturor celor care au de rezolvat o anumită problemă iar metodele clasice de lucru nu le sunt de folos.

Autorul

Metode Numerice 1 CUPRINS L 1. Rezolvarea ecuaţiilor prin metoda bisecţiei (înjumătăţirii intervalului) 3 L 2. Rezolvarea ecuaţiilor prin metoda lui Newton (tangentei) 7 L 3. Rezolvarea ecuatiilor prin metoda aproximaţiilor succesive 11 L 4. Rezolvarea ecuatiilor prin metoda secantei 15 L 5. Rezolvarea sistemelor de ecuaţii liniare prin metoda lui Gauss 19 L 6. Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss-Jordan 27 L 7. Descompunerea unei matrice în produs TS 35 L 8. Rezolvarea sistemelor de ecuaţii liniare prin metoda iterativă Jacobi 40

Indrumar de laborator 2 L 9. Rezolvarea sistemelor de ecuaţii liniare prin metoda iterativă Gauss-Seidel 44 L 10. Polinomul de interpolare Lagrange 48 L 11. Polinomul de interpolare Newton cu diferenţe divizate 51 L 12. Polinoame de interpolare Newton cu pas constant 55 L 13. Aproximarea funcţiilor prin metoda celor mai mici pătrate. Dreapta de regresie 63 L 14. Integrare numerică. Formulele trapezelor 67 L 15. Integrare numerică. Formula lui Simpson 71 L 16. Integrare numerică. Formulele lui Gauss 74

Metode Numerice 3 LUCRAREA 1

Rezolvarea ecuaţiilor prin metoda bisecţiei (înjumătăţirii intervalului)

Considerăm că ecuaţia 0)x(f = are o singură rădăcină în intervalul şi că funcţia f este continuă pe acest interval. Această presupunere este valabilă, în condiţiile parcurgerii primei etape, aceea de separare a unei singure rădăcini într-un anumit interval.

)b,a( 00

Fie ε eroarea admisă pentru soluţia ecuaţiei. Din punct de vedere grafic, rezolvarea ecuaţiei prin această metodă, este ilustrată în Fig.1.1. Intervalul iniţial se împarte în două părţi egale prin punctul: )b,a( 00

2

bac 00

0+

= (1.1)

f(x)

x

A0

B0

a0a1

xR

c0

b1

b0

Fig.1.1.

Se efectuează apoi produsul . Intervalul care conţine în continuare soluţia se notează . Situaţiile care pot apărea sunt următoarele:

)c(f)a(f 00

)b,a( 11

(1.2) ⎪⎩

⎪⎨

⎧

==∈>==

==∈<

.bb,canoteazase),b,c(x,0cx,0

cb,aanoteazase),c,a(x,0:)c(f)a(f

010100R

0R

010100R

00

Indrumar de laborator 4 În situaţia prezentată în grafic avem < 0. Cu intervalul

se procedează în mod asemănător. )c(f)a(f 00

)b,a( 11 Rezultă două şiruri, { } Nnna ∈ crescător (sau constant pe porţiuni) şi descrescător (sau constant pe porţiuni) soluţia aflându-se în permanenţă în intervalul .

{ }nb Rx)b,a( nn

Lungimea acestui interval este:

00nnn ab21ab −=− . (1.3)

Numărul minim de iteraţii nmin se determină funcţie de precizia ε impusă calculelor:

ε<− nn ab . (1.4)

Rezultă pentru nmin o expresie de forma:

,1ab

logn 002min +⎥

⎦

⎤⎢⎣

⎡

ε

−= (1.5)

unde [⋅] reprezintă funcţia parte întreagă. Valoarea nmin nu depinde de complexitatea ecuaţiei care se rezolvă ci numai de lungimea intervalului iniţial şi de precizia ε impusă. Orice valoare cuprinsă în intervalul final poate fi considerată ca fiind soluţie aproximativă pentru ecuaţia dată. De obicei se consideră

)b,a( nn

),ba(21x nnR +≅ (1.6)

mijlocul ultimului interval determinat.

Metode Numerice 5 Programul în MATLAB: L1_bisectie.m function rad_ec = bisectie_ec(a0,b0,max_er,max_it,index_f) % Rezolvarea ecuatiei f(x)=0 prin metoda bisectiei % ( metoda injumatatirii intervalului ); % a0, b0 = intervalul initial ce contine o radacina a ecuatiei; % max_er = valoarea maxima admisa a erorii de aproximare; % max_it = numarul maxim de iteratii admis; % index_f = selectorul de functii; ecuatia este f(x)=0. if a0 >= b0 disp('a0 si b0 sunt incorect date; STOP!') return end format long a = a0; b = b0; fa = f(a,index_f); fb = f(b,index_f); if sign(fa*fb) > 0 disp('f(a0) si f(b0) nu au semne contrare; STOP!') return end c = (a+b)/2; contor = 0; while b-c > max_er & contor < max_it Pasul = contor + 1; contor = contor + 1; fc = f(c,index_f); Pasul disp('Cele trei puncte a, c, b sunt:'), [a c b] disp('Valorile functiei in aceste puncte sunt:'), [fa fc fb] disp('Apreciati pozitia radacinii ecuatiei si eroarea ') pause if fc == 0 disp(' ') disp('Radacina este exact in mijlocul intervalului, STOP! ')

Indrumar de laborator 6 format short Pasul_final = Pasul format long Radacina_aprox = c format short e Eroarea = 0 return end if sign(fa*fc) < 0 b = c; fb = fc; else a = c; fa = fc; end c = (a+b)/2; eroarea = b-c pause end format short Pasul_final = Pasul format long Radacina_aprox = c format short e Eroarea = b-c %%%%%%%%%%%%%%%%%%%%%%%%%%%% function fct = f(x,index) % functia f, ce defineste ecuatia f(x)=0 : switch index case 1 fct = x.^6 - x - 1; case 2 fct = x - exp(-x); case 3 fct = x.^3-x-2; case 4 fct = x.^3-2; end


Rezolvarea ecuaţiilor prin metoda lui Newton (tangentei)

Considerăm că ecuaţia 0)x(f = conţine în intervalul o singură soluţie . Grafic, rezolvarea ecuaţiei prin această metodă este ilustrată în Fig.2.1. De asemenea, considerăm că pe acest interval derivatele f' şi f" păstrează semn constant, deci f este strict monotonă şi

)b,a( 00

Rx

f(x)

x

A0

B0

x0=a0 x1

A1

xR b0

x2

Fig.2.1.

nu are punte de inflexiune. Ea presupune aproximarea soluţiei exacte printr-un şir de valori , ... obţinute prin intersecţia tangentelor duse la graficul funcţiei f în punctele A

Rx

1x 2x0, A1,... cu axa absciselor.

Punctul iniţial se alege ca fiind una din extremităţile intervalului şi anume aceea care îndeplineşte condiţia:

0x)b,a( 00

(2.1) 0)x("f)x(f 00 >

Această condiţie ne asigură că intersecţia tangentei la grafic cu axa Ox se va afla în interiorul intervalului iniţial. Considerăm punctul generic ))x(f,x(A kkk situat pe graficul funcţiei f. Ecuaţia tangentei în acest punct la graficul lui f este:

)xx)(x('f)x(fy kkk −=− . (2.2)

Indrumar de laborator 8 Intersecţia cu axa Ox este punctul 1kx + obţinut pentru y = 0:

,....2,1,0k;)x('f)x(f

xxk

kk1k =−=+ (2.3)

Această expresie este formula iterativă a lui Newton. Se notează:

)x('f)x(fx)x( −=ϕ . (2.4)

Funcţia ϕ se numeşte funcţia de iterare a lui Newton. Cu noua notaţie, relaţia de recurenţă dată de (2.3) se poate pune sub forma:

,...2,1,0k),x(x k1k =ϕ=+ (2.5)

Notând Rkk xx −=ε , eroarea de aproximare la pasul k, se obţine:

.!2

)x(" R2k1k

ϕ⋅ε=ε + (2.6)

Această relaţie arată că eroarea de aproximare evoluează după o lege pătratică. Metoda lui Newton are deci o viteză mare de convergenţă. Dezavantajul metodei constă în faptul că în cadrul fiecărui pas este necesar calculul derivatei funcţiei în punctul respectiv. Programul în MATLAB: L2_newton.m function rad_ec = newton_ec(x0,max_er,max_it,index_f) % Rezolvarea ecuatiei f(x)=0 prin metoda lui Newton (tangentei) % x0 = aproximarea initiala pentru radacina ecuatiei; % max_er = valoarea maxima admisa a erorii de aproximare; % max_it = numarul maxim de iteratii admis; % index_f = selectorul de functii; ecuatia este f(x)=0. eroare = 1; contor = 1;

Metode Numerice 9 while abs(eroare) > max_er & contor <= max_it fx = f(x0, index_f); dfx = deriv_f(x0, index_f); if dfx == 0 disp('Valoarea derivatei este zero. STOP !') return end format short Pasul = contor disp('Valoarea aproximativa a radacinii este :') format long x1 = x0 - fx/dfx; x=x1 format short e eroare = x1 - x0 pause x0 = x1; contor = contor + 1; end if contor > max_it disp('Numarul de iteratii stabilit initial a fost atins.') disp('Acesta este prea mic pentru a se gasi o radacina') disp('cu precizia impusa pentru ecuatia data !') disp(' ') else disp('Nivelul de precizie impus a fost atins :') format short Pasul format long Radacina_aprox = x1 format short e Eroarea = eroare end %%%%%%%%%%%%%%%%%%%%%%%%%%%% function fct = f(x, index) % Definirea functiei f ce a generat ecuatia : switch index

Indrumar de laborator 10 case 1 fct = x.^6 - x - 1; case 2 fct = x - exp(-x); case 3 fct = x.^3-x.^2-x-1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%% function d_fct = deriv_f(x, index) % Derivata functiei f este: switch index case 1 d_fct = 6*x.^5 - 1; case 2 d_fct = 1 + exp(-x); case 3 d_fct = 3*x.^2-2*x-1; end


Rezolvarea ecuatiilor prin metoda aproximaţiilor succesive

Pentru rezolvarea ecuaţiei 0)x(f = prin această metodă, printr-un artificiu de calcul se aduce această ecuaţie la forma echivalentă:

)x(x ϕ= . (3.1)

Dacă este soluţie pentru ecuaţia Rx 0)x(f = , ea este soluţie şi pentru ecuaţia (3.1), deci vom avea:

)x(x RR ϕ= . (3.2)

Fie o aproximaţie iniţială pentru soluţia ecuaţiei. )b,a(x 000 ∈ Metoda constă în aproximarea soluţiei exacte printr-un şir de valori , ..., construit pe relaţia iterativă

Rx

1x 2x

,...2,1,0k),x(x k1k =ϕ=+ (3.3)

Considerăm că funcţia ϕ, derivabilă pe intervalul satisface condiţia:

)b,a( 00

)b,a(x)(1)x(' 00∈∀<λ≤ϕ . (3.4)

În aceste condiţii, şirul de iterare definit de relaţia (3.3) este convergent către soluţia a ecuaţiei. Rx Funcţia ϕ este continuă iar prin trecere la limită în relaţia de recurenţă (3.3) se obţine:

)xlim(xlim kk

1kk ∞→

+∞→

ϕ= (3.5)

adică:

).x(x RR ϕ= (3.6)

Rezultă că limita şirului de puncte , ,... reprezintă soluţia ecuaţiei.

0x , 1x 2x

Indrumar de laborator 12 Din punct de vedere grafic, rezolvarea ecuaţiei (3.1) prin această metodă se desfăşoară aşa cum este prezentat în Fig.3.1 şi în Fig.3.2.

xR x2 x1 x0 x

yy=x

y=φ(x)

xR x2x1 x0 x

yy=xy=φ(x)

(a) (b) 0<λ<1, (Convergenţă)

Fig.3.1.

xR x2x1 x0 x

y y=φ(x)y=x

xRx1 x0

yy=ϕ(x)

y=x

(a) (b) λ>1 (Divergenţă)

Fig.3.2.

Notând eroarea Rkk xx −=ε obţinem:

).x(' Rk1k ϕ⋅ε=ε + (3.7)

Se obţine un grad de convergenţă liniar; eroarea de la pasul (k+1) este proporţională cu eroarea de la pasul k.

Metode Numerice 13 Programul în MATLAB: L3_aprox_succ.m function rad_ec = aprox_succ_ec(x0,max_er,max_it,index_fi) % Rezolvarea ecuatiei f(x)=0 prin metoda aproximatiilor succesive % x0 = aproximarea initiala pentru radacina ecuatiei; % max_er = valoarea maxima admisa a erorii de aproximare; % max_it = numarul maxim de iteratii admis; % index_fi = selectorul de functii; ecuatia este x=fi(x). neconv = 0; eroare = 1; eroare_0 = 10; contor = 0; while abs(eroare) > max_er & contor < max_it contor = contor + 1; Pasul = contor format long x1 = fi(x0,index_fi); disp('Valoarea aproximativa a radacinii este :') x = x1 format short e eroare = x1-x0 if abs(eroare)>abs(eroare_0) neconv = neconv + 1; disp('Metoda este neconvergenta') neconv if neconv == 3 disp(' ') disp(' Va rugam sa abandonati ! ') disp(' Nu luati in considerare rezultatele ! ') disp(' ') return end end eroare_0=eroare; x0=x1; pause end if contor >= max_it disp('Numarul de iteratii stabilit initial, atins.') disp('Acesta este prea mic pentru a se gasi o radacina')

Indrumar de laborator 14 disp('cu precizia impusa pentru ecuatia data !') disp(' ') else disp('Nivelul de precizie impus a fost atins :') format short Pasul format long Radacina_aprox = x1 format short e Eroarea = eroare end %%%%%%%%%%%%%%%%%%%%%%%%%%%% function fct = fi(x,index) % Definirea functiei fi: switch index case 1 fct = x.^3-2; case 2 fct = (x+2)^(1/3); end


Rezolvarea ecuatiilor prin metoda secantei

Considerăm ecuaţia 0)x(f = pentru care am separat în intervalul o singură soluţie. )b,a( 00Din punct de vedere grafic, rezolvarea ecuaţei prin această metodă este ilustrată în Fig.4.1.

A0

a0 x0

f(x)

xR

x1

b0

B0 Fig.4.1.

Se unesc punctele şi printr-o dreaptă de ecuaţie: 0A 0B

)ax(ab

)a(f)b(f)a(fy 0

00

000 −⋅

−−

=− . (4.1)

Intersecţia dreptei cu axa Ox (y = 0) ne dă punctul 0x :

)a(f)b(f

)a(fb)b(fax

00

00000 −

−= (4.2)

Punctul împarte intervalul în două subintervale din care numai unul va conţine în continuare soluţia .

0x )b,a( 00

Rx Pentru a determina care subinterval conţine soluţia se efectuează produsul:

Rx

Indrumar de laborator 16

. (4.3) ⎩⎨⎧

==∈>==∈<

010100R

010100R00 bb,xa.notse),b,x(x,0

xb,aa.notse),x,a(x,0:)x(f)a(f

Cu intervalul se reia procedeul descris pentru . Rezultă un şir de puncte , ,…, , unde:

)b,a( 11 )b,a( 00

1x 2x nx

)a(f)b(f

)a(fb)b(fax

nn

nnnnn −

−= . (4.4)

Aceste două şiruri, { } Nnna ∈ şi { } Nnnb ∈ au o limită comună, soluţia a ecuaţiei. Rx La toate metodele iterative prezentate procedeul de calcul se opreşte atunci când există un indice n pentru care:

ε<− −1nn xx , (4.5)

ε fiind eroarea admisă. Acest test nu ne asigură că am determinat o valoare aproximativă pentru soluţia situată faţă de valoarea exactă la o distanţă mai mică decât ε. Pentru a fi mai siguri, în această situaţie se consideră o condiţie de oprire a calculelor de forma:

.10

xx 1nnε

<− − (4.6)

Programul în MATLAB: L4_secanta.m function rad_ec = secanta_ec(a0,b0,max_er,max_it,index_f) % Rezolvarea ecuatiilor prin metoda secantei % a0, b0 = intervalul initial ce contine o radacina a ecuatiei; % max_er = valoarea maxima admisa a erorii de aproximare; % max_it = numarul maxim de iteratii admis; % index_f = selectorul de functii; ecuatia este f(x)=0. eroare = 1; fa0 = f(a0,index_f); fb0 = f(b0,index_f);

Metode Numerice 17 contor = 0; while abs(eroare) > max_er & contor <= max_it contor = contor + 1; Pasul = contor disp(' Punctul de intersectie cu abscisa este :') format long x0=(a0*fb0-b0*fa0)/(fb0-fa0); x = x0 fx0 = f(x0,index_f); if fa0*fx0 < 0 b0 = x0; else a0 = x0; end fa0 = f(a0,index_f); fb0 = f(b0,index_f); % Eroarea se calculeaza incepand de la pasul 2. if contor > 1 format short e eroare = x0 - x1 end pause x1 = x0; end if contor > max_it disp('Numarul de iteratii stabilit initial, atins.') disp('Acesta este prea mic pentru a se gasi o radacina') disp('cu precizia impusa pentru ecuatia data !') disp(' ') else disp('Nivelul de precizie impus a fost atins :') format short Pasul format long Radacina_aprox = x0 format short e Eroarea = eroare end

Indrumar de laborator 18 %%%%%%%%%%%%%%%%%%%%%%%%%%%% function fct = f(x,index) % Definirea functiei f : switch index case 1 fct = x.^6 - x - 1; case 2 fct = x - exp(-x); end


Rezolvarea sistemelor de ecuaţii liniare prin metoda lui Gauss

Se consideră un sistem de n ecuaţii cu n necunoscute de forma:

. (5.1.)

⎪⎪⎩

⎪⎪⎨

⎧

=+++

=+++=+++

nnnn22n11n

2nn2222121

1nn1212111

bxa...xaxa..................................................

bxa...xaxabxa...xaxa

Dacă se fac notaţiile:

(5.2)

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

n

2

1

n

2

1

nn2n1n

n22221

n11211

b...bb

b;

x...xx

x;

a...aa.............a...aaa...aa

A

sistemul poate fi scris în formă matriceală:

.bAx = (5.3) Această metodă permite rezolvarea sistemului (5.1) prin transformarea, în cadrul unui anumit număr de paşi, a matricei pătratice A a sistemului într-o matrice superior triunghiulară. Se consideră 0a11 ≠ . În caz contrar, se reordonează şi renumerotează ecuaţiile sistemului astfel încât această condiţie să fie îndeplinită. Numărul de paşi pentru un sistem de n ecuaţii cu n necunoscute este (n-1). Pentru a putea urmări modul în care se efectuează triunghiularizarea sistemului este necesară descrierea operaţiilor care se efectuează în cadrul fiecărui pas în parte. Pasul 1: Elementul 0a11 ≠ se numeşte pivot în cadrul acestui pas, iar ecuaţia "1" care rămâne nemodificată se numeşte ecuaţie de pivotare. Necunoscuta este eliminată din ultimele (n-1) ecuaţii prin 1x


înmulţirea ecuaţiei 1 a sistemului cu rapoartele 11

1iaa şi scăderea acesteia

din ecuaţia "i", cu n,2i = . Pasul k, 1n,2k −= : Primele k ecuaţii ale sistemului obţinut la pasul anterior rămân nemodificate. Ecuaţia "k" este ecuaţia de pivotare iar elementul situat pe diagonala principală în această ecuaţie este pivotul. Dacă acest element este nul, este necesară reordonarea şi renumerotarea ultimelor n-(k-1) ecuaţii ale sistemului de la pasul precedent astfel încât să avem un pivot nenul, după care se procedează la transcrierea primelor k ecuaţii. În cadrul acestui pas, necunoscuta kx este eliminată din ultimele (n-k) ecuaţii într-o manieră asemănătoare cu cea de la pasul 1. De această dată, ecuaţia de pivotare se înmulţeşte cu un raport care are la numărător coeficientul lui kx din ecuaţia curentă iar la numitor pivotul din cadrul acestui pas şi se scade din ecuaţia curentă. În acest fel necunoscuta kx este eliminată din ultimele (n – k) ecuaţii. La ultimul pas, (n-1), necunoscuta este eliminată numai din ultima ecuaţie a sistemului obţinut la pasul anterior, primele (n-1) ecuaţii rămânând nemodificate.

1nx −

Pentru a putea scrie relaţiile de recurenţă care leagă matricele sistemului la trecerea de la un pas la celălalt, convenim să ataşăm matricei iniţiale A indicele superior "1", ea devenind:

( )n,1j,i

)1(ij

)1( aA=

= . (5.4)

Numai elementele care se modifică în cadrul unui pas îşi indexează indicele superior cu o unitate. Notăm:

(5.5) )k()k( bxA =⋅

forma sistemului înainte de eliminarea necunoscutei kx din ultimele (n-k) ecuaţii. Astfel, pentru k=1 avem:

AA )1( = ; . (5.6) bb )1( =

Metode Numerice 21 Pentru n,2k = , elementele matricei A(k) şi a vectorului b(k) se vor calcula în funcţie de elementele matricelor de la pasul precedent astfel:

⎪⎪⎪⎪

⎩

⎪⎪⎪⎪

⎨

⎧

≥≥−

−≤≥

−≤

=

−−−

−−

−−−

−

kj,kipentru,a

aaa

1kj,kipentru,0

1kipentru,a

a

)1k(1k,1k

)1k(j,1k

)1k(1k,i)1k(

ij

)1k(ij

)k(ij (5.7)

⎪⎪⎩

⎪⎪⎨

⎧

≥⋅−

−≤

=−−−

−−

−−−

−

kipentru,ba

ab

1kipentru,b

b )1k(1k)1k(

1k,1k

)1k(1k,i)1k(

i

)1k(i

)k(i . (5.8)

Relaţiile (5.7) şi (5.8) reprezintă de fapt rezultatul înmulţirii în cadrul pasului (k-1) a ecuaţiei (k-1) a sistemului:

(5.9) )1k()1k( bxA −− =⋅

cu raportul )1k(1k,1k

)1k(1k,i aa −

−−−− , şi scăzând rezultatul din ecuaţia i, pentru

i ≥ k. Este eliminată astfel necunoscuta xk-1 din ultimele n-(k-1) ecuaţii. La sfârşitul acestui pas, matricea A(k) şi vectorul termenilor liberi b(k) arată astfel:


. (5.10) ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

=

−−

−−

−−

−−−

−

−

)k(n

)k(k

)1k(1k

)2(2

)1(1

)k(

)k(nn

)k(nk

)k(kn

)k(kk

)1k(n,1k

)1k(k,1k

)1k(1k,1k

)2(n2

)2(k2

)2(1k,2

)2(22

)1(n1

)1(k1

)1(1k,1

)1(12

)1(11

)k(

b...

bb

...bb

b

;

a...a0...00.....................

a...a0...00

a...aa...00.....................

a...aa...a0

a...aa...aa

A

În tot acest proces am presupus n,2k,0a )1k(1k,1k =≠−−− .

În cadrul ultimului pas avem k = n-1, necunoscuta este eliminată din ultima ecuaţie, obţinându-se sistemul echivalent sub formă triunghiulară:

1nx −

(5.11)

⎪⎪⎪⎪

⎩

⎪⎪⎪⎪

⎨

⎧

=

=++

=++++

=+++++

)n(nn

)n(nn

)k(kn

)k(knk

)k(kk

)2(2n

)2(n2k

)2(k22

)2(22

)1(1n

)1(n1k

)1(k12

)1(121

)1(11

bxa

..........................................bxa...xa

..............................................................bxa...xa...xa

bxa...xa...xaxa

Găsirea efectivă a soluţiei sistemului se face printr-un procedeu de eliminare inversă aplicat sistemului triunghiular (5.11):

Metode Numerice 23

⎪⎪⎪⎪⎪⎪

⎩

⎪⎪⎪⎪⎪⎪

⎨

⎧

−−=

−

=

−=

=

∑+=

−−−

−−

−−

−

1,2,...3n,2nk,a

xab

x

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

a

xabx

abx

)k(kk

j

n

1kj

)k(kj

)k(k

k

)1n(1n,1n

n)1n(n,1n

)1n(1n

1n

)n(nn

)n(n

n

. (5.12)

În acest fel se obţine soluţia sistemului. Programul în MATLAB: L5_gauss_1.m function sol_sist = gauss_1(A,b) % Rezolvarea sistemelor de ecuatii liniare prin metoda % Gauss de eliminare, fara schimbarea pivotului [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n=length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end format disp('Matricea sistemului este :') A disp('Vectorul termenilor liberi este :') b pause

Indrumar de laborator 24 for k=1:n-1, for i=k+1:n, if A(k,k)==0 disp('Pivotul este 0: EROARE') return else P(i)=A(i,k)/A(k,k); A(i,k)=0; end for j=k+1:n, A(i,j)=A(i,j)-P(i)*A(k,j); end b(i)=b(i)-P(i)*b(k); end Pasul=k; Pasul disp('Matricea sistemului este acum :') A disp('Vectorul termenilor liberi este :') b pause end x(n)=b(n)/A(n,n); for i=n-1:-1:1, T=0; for j=i+1:n, T=T+A(i,j)*x(j); end x(i)=(b(i)-T)/A(i,i); end disp('Solutia sistemului este :') sol_sist=x; Pentru creşterea preciziei calculelor numerice este necesară reordonarea ecuaţiilor în sistem în cadrul fiecărui pas, astfel încât, de fiecare dată, pivotul să aibă o valoare maximă, în modul. Reordonarea ecuaţiilor se face prin înlocuirea de tip „pereche” a ecuaţiei de pivotare ce rezultă în mod normal la finele unui anumit pas cu o ecuaţie situată sub aceasta şi care are, în modul, cea mai mare valoare a primului coeficient nenul din ecuaţie. Acest coeficient devine noul pivot iar întreaga ecuaţie va înlocui ecuaţia de pivotare obţinută în mod natural.

Metode Numerice 25 Programul în MATLAB:L5_gauss_2.m function sol_sist = gauss_2(A,b) % Rezolvarea sistemelor de ecuatii liniare prin metoda % Gauss de eliminare, cu schimbarea pivotului [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n=length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end format piv = (1:n)'; disp('Matricea sistemului este ;') A disp('Vectorul termenilor liberi este ;') b pause % Valoarea maxima a elementului de pe coloana si % linia careia ii corespunde este: for k=1:n-1, [col_max index] = max(abs(A(k:n,k))); index = index + k-1; if index ~= k % Schimba liniile k si index intre ele temp_A = A(k,k:n); A(k,k:n) = A(index,k:n); A(index,k:n) = temp_A; temp_A pause temp_b = b(k);

Indrumar de laborator 26 b(k) = b(index); b(index) = temp_b; temp = piv(k); piv(k) = piv(index); piv(index) = temp; end disp('Matricea sistemului cu liniile invesate este :') A pause for i=k+1:n, if A(k,k)==0 disp('Pivotul este 0: EROARE') return else P(i)=A(i,k)/A(k,k); A(i,k)=0; end for j=k+1:n, A(i,j)=A(i,j)-P(i)*A(k,j); end b(i)=b(i)-P(i)*b(k); end Pasul=k; Pasul disp('Matricea sistemului este acum :') A disp('Vectorul termenilor liberi este :') b pause end x(n)=b(n)/A(n,n); for i=n-1:-1:1, T=0; for j=i+1:n, T=T+A(i,j)*x(j); end x(i)=(b(i)-T)/A(i,i); end disp('Solutia sistemului este :') sol_sist = x(1:n)


Rezolvarea sistemelor de ecuaţii liniare prin metoda Gauss-Jordan

Această metodă permite rezolvarea sistemului (5.1) prin transformarea sa într-un sistem echivalent a cărui matrice este diagonală. Rezolvarea presupune parcurgerea a n paşi. În cadrul pasului generic kP se transcrie numai ecuaţia având indicele k din sistemul de la pasul anterior. Necunoscuta kx este eliminată din toate celelalte (n-1) ecuaţii, printr-un procedeu asemănător celui de la metoda lui Gauss.

k-1

i

)1k(1k,1ka −−−

)1k(1k,ia −−

)1k(j,1ka −

−

)1k(j,ia −

k-1 j

Liniile

:A )1k( − Coloanele

Lini

a (k

-1) răm

âne

nem

odifi

cată

Coloana (k-1) va avea zerouri în

afara elementului )1k(1k,1ka −−−

Coloanemodificate

Primele (k-2)coloane rămânnemodificate

(6.1)

Matricea sistemului va purta la finele fiecărui pas un indice superior egal cu numărul pasului care urmează. În cadrul fiecărui pas toate elementele matricei sistemului îşi indexează indicele superior cu o unitate.

Indrumar de laborator 28 Relaţiile de recurenţă ce permit trecerea de la matricea A(k-1) la matricea A(k) pot fi scrise mai uşor dacă se evidenţiază, pe elementele matricei A(k-1) operaţiile care se execută, aşa cum este prezentat în (6.1): Ele au forma:

⎪⎪⎪⎪

⎩

⎪⎪⎪⎪

⎨

⎧

−≤=

−≥=

−≥−≠

−=

−

−−−

−−−

−−

−−−

2kjpentruaa

1kjpentru,aa

1kj1ki

pentru,a

aaaa

)1k(ij

)k(ij

)1k(j,1k

)k(j,1k

)1k(1k,1k

)1k(j,1k

)1k(1k,i)1k(

ij)k(

ij

(6.2)

şi respectiv:

⎪⎪⎩

⎪⎪⎨

⎧

=

−≠⋅−=

−−−

−−−

−−

−−−

)1k(1k

)k(1k

)1k(1k)1k(

1k,1k

)1k(1k,i)1k(

i)k(

i

bb

1kipentru,ba

abb

. (6.3)

În final se obţine un sistem diagonal de forma:

. (6.4)

⎪⎪⎩

⎪⎪

⎨

⎧

=

==

++

++

++

)1n(nn

)1n(nn

)1n(22

)1n(22

)1n(11

)1n(11

bxa............

b......xab.........xa

Soluţia sistemului este evident dată de relaţiile:

n,1ipentru,a

bx

)1n(ii

)1n(i

i ==+

+

. (6.5)

Metoda Gauss-Jordan permite şi calculul matricei inverse A-1 dacă se procedează în felul următor: - se formează o matrice de dimensiuni n×(2n+1) compusă din:

• pe primele n coloane matricea A;

Metode Numerice 29

• pe următoarele n coloane matricea ; nI• pe ultima coloană vectorul b;

- se aplică procedeul de rezolvare dat de metoda Gauss-Jordan numai că operaţiile descrise pentru matricea sistemului se extind asupra întregii matrice de dimensiuni n×(2n+1). Se efectuează aceste operaţii până când în locul matricei A apare matricea diagonală. Ultima operaţie care se efectuează asupra acestei matrice extinse este aceea de a împărţi fiecare linie la elementul situat pe diagonala principală a matricei de pe primele n coloane. Ca urmare, în locul matricei diagonale cu elemente oarecare ce se află pe primele n coloane, vom obţine matricea unitate . Pe următoarele n coloane vom obţine matricea inversă A

nI-1 iar pe ultima coloană, a (2n+1) - a, vom obţine chiar

soluţia sistemului. Sintetic, metoda poate fi descrisă astfel:

[ ] [ ]x|A|Ib|I|A 1nn

−→ (6.6)

x fiind soluţia sistemului cu componentele n,1i,x i = . Programul in MATLAB: L6_G_J_1.m function sol_sist = G_J_1(A,b) % Rezolvarea sistemelor de ecuatii liniare prin % metoda Gauss-Jordan de eliminare, % fara schimbarea pivotului [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n=length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end

Indrumar de laborator 30 format disp('Matricea sistemului este :') A disp('Vectorul termenilor liberi este :') b pause for k=1:n, for i=1:n, if i~=k if A(k,k)==0 disp('Pivotul este 0: EROARE') return else P(i)=A(i,k)/A(k,k); A(i,k)=0; end for j=k+1:n, A(i,j)=A(i,j)-P(i)*A(k,j); end b(i)=b(i)-P(i)*b(k); end end Pasul=k; Pasul disp('Matricea sistemului este acum :') A disp('Vectorul termenilor liberi este :') b pause end x(n)=b(n)/A(n,n); for i=n-1:-1:1, T=0; for j=i+1:n, T=T+A(i,j)*x(j); end x(i)=(b(i)-T)/A(i,i); end disp('Solutia sistemului este :') sol_sist=x(1:n)

Metode Numerice 31 Ca şi la metoda Gauss, pentru creşterea preciziei calculelor numerice este necesară reordonarea ecuaţiilor în sistem în cadrul fiecărui pas, astfel încât, de fiecare dată, pivotul să aibă o valoare maximă, în modul. Programul in MATLAB: L6_G_J_2.m function sol_sist = G_J_2(A,b) % Rezolvarea sistemelor de ecuatii liniare prin % metoda Gauss-Jordan de eliminare, % cu schimbarea pivotului [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n=length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end format piv = (1:n)'; disp('Matricea sistemului este ;') A disp('Vectorul termenilor liberi este ;') b pause for k=1:n, [col_max index] = max(abs(A(k:n,k))); index = index + k-1; if index ~= k % Schimba liniile k si index intre ele temp_A = A(k,k:n); A(k,k:n) = A(index,k:n);

Indrumar de laborator 32 A(index,k:n) = temp_A; temp_A pause temp_b = b(k); b(k) = b(index); b(index) = temp_b; temp = piv(k); piv(k) = piv(index); piv(index) = temp; end disp('Matricea sistemului cu liniile inversate :') A pause for i=1:n, if i~=k if A(k,k)==0 disp('Pivotul este 0: EROARE') return else P(i)=A(i,k)/A(k,k); A(i,k)=0; end for j=k+1:n, A(i,j)=A(i,j)-P(i)*A(k,j); end b(i)=b(i)-P(i)*b(k); end end Pasul=k; Pasul disp('Matricea sistemului este acum :') A disp('Vectorul termenilor liberi este :') b pause end x(n)=b(n)/A(n,n); for i=n-1:-1:1, T=0; for j=i+1:n, T=T+A(i,j)*x(j); end

Metode Numerice 33 x(i)=(b(i)-T)/A(i,i); end disp('Solutia sistemului este :') sol_sist=x(1:n) Se poate determina matricea inversă matricii A a sistemului, prin alipirea la aceasta, în partea dreaptă, a unei matrice unitate de acelaşi ordin şi prin aplicarea algoritmului de rezolvare a sistemului prin metoda Gauss-Jordan, fără schimbarea pivotului. Odată cu soluţia sistemului se obtine şi matricea inversă. Programul in MATLAB: L6_G_J_inv.m function sol_sist = g_j_inv_1(A,b) % Matricea A se extinde, prin alipirea matricii unitate % de acelasi ordin cu A format disp('Matricea initiala a sistemului este ;') A disp('Vectorul termenilor liberi este ;') b pause n=length(b); % Se formeaza matricea extinsa : for i=1:n for j=n+1:2*n if j==i+n A(i,j)=1; else A(i,j)=0; end end end disp('Matricea extinsa asociata sistemului este :') A_extinsa = A disp('Vectorul termenilor liberi este :') b pause

Indrumar de laborator 34 for k=1:n, for i=1:n, if i==k else if A(k,k)==0 disp('Pivotul este 0: EROARE') else P(i)=A(i,k)/A(k,k); A(i,k)=0; end for j=k+1:2*n, A(i,j)=A(i,j)-P(i)*A(k,j); end b(i)=b(i)-P(i)*b(k); end end Pasul=k; Pasul disp('Matricea sistemului este acum :') A disp('Vectorul termenilor liberi este :') b pause end for i=1:n, for j=1:2*n inv_A(i,j)=A(i,j)/A(i,i); end x(i)=b(i)/A(i,i); end disp('Matricea finala extinsa este :') A_final=A pause disp('Matricea inversa este urmatoarea:') inv_A = A(:,n+1:2*n) disp('Solutia sistemului este :') sol_sist = x(1:n)


Descompunerea unei matrice în produs TS

Această metodă se bazează pe faptul că orice matrice nesingulară A poate fi descompusă într-un produs de două matrice, una inferior triunghiulară, cealaltă superior triunghiulară, introducând în acelaşi timp un algoritm de calcul pentru elementele celor două matrice în funcţie de elementele matricei A. Se consideră sistemul (5.1), matricea A fiind matricea sistemului, nesingulară. Descompunerea va fi de forma:

STA ⋅= (7.1)

unde:

jipentru0s,)s(Sn,1ipentru1t

jipentru0t,)t(T

ijn,1j,iij

ii

ijn,1j,iij

>==

⎪⎩

⎪⎨⎧

==

<==

=

= . (7.2)

Se pot determina relaţiile pentru calculul elementelor matricelor T şi S dacă se impune coincidenţa între elementele matricei A şi cele ale matricei produs TS. Considerăm un element oarecare al matricei A, kja având k ≤ j. Acest element se obţine prin înmulţirea liniei k din matricea T cu coloana j din matricea S. Deoarece k ≤ j se poate scrie:

. (7.3) ∑=

=k

1ppjkpkj sta

Separând elementul din (7.3) vom avea: kjs

. (7.4) )jk(,stas1k

1ppjkpkjkj ≤−= ∑

−

=

Indrumar de laborator 36 Considerăm acum un element ki,aik > . Se poate scrie:

. (7.5) ∑=

=k

1ppkipik sta

Separând elementul ikt din (7.5) se obţine:

)ki(,stas1t

1k

1ppkipik

kkik >

⎥⎥⎦

⎤

⎢⎢⎣

⎡−= ∑

−

=. (7.6)

Relaţiile (7.4) şi (7.6) permit determinarea tuturor elementelor matricelor T şi S presupunând că se cunosc elementele primei linii din S şi ale primei coloane din T. Este necesar să se definească deci elementele primei linii din S:

n,1j,as j1j1 == (7.7)

şi ale primei coloane din T:

)0a(n,1i,aat 11

11

1i1i ≠== . (7.8)

Ordinea în care se pot determina elementele matricelor T şi S este următoarea: o linie din S, o coloană din T, o linie din S, o coloană din T,..., până la epuizarea elementelor celor două matrice. Având determinate elementele matricelor produs se poate proceda la rezolvarea sistemului:

bTSx,bAx =⇒= . (7.9)

Se notează:

ySx = , (7.10)

astfel încât rezolvarea sistemului iniţial (5.1) presupune rezolvarea a două sisteme triunghiulare de forma:

. (7.11) ⎩⎨⎧

==

bTyySx

Al doilea sistem din (7.11) poate fi scris detaliat astfel:

Metode Numerice 37

. (7.12)

⎪⎪⎩

⎪⎪⎨

⎧

=+++

=+=

nnnn22n11n

2222121

1111

byt...ytyt..........................................

b..................ytytb..............................yt

Aflarea vectorului necunoscut y se face printr-un procedeu de eliminare directă, cu relaţiile:

⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

=

−

=

=

∑−

= n,...,3,2k,t

ytb

y

tby

kk

1k

1ppkpk

k

11

11

. (7.13)

Având cunoscut vectorul y se poate trece la rezolvarea primului sistem din (7.11) care se poate scrie:

. (7.14)

⎪⎪⎩

⎪⎪⎨

⎧

=

=++=+++

nnnn

2nn2222

1nn1212111

yxs...........................

yxs...xsyxs...xsxs

Soluţia se află prin eliminarea inversă:

⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

−=

−

=

=

∑+= 1,...,2nk;

s

xsy

x

sy

x

kk

n

1kppkpk

k

nn

nn

. (7.15)

Această metodă, care presupune descompunerea matricei A în produs TS, mai apare în literatura de specialitate ca fiind descompunerea LR, corespondenţa între diversele notaţii fiind evidentă. Programul in MATLAB: L7_desc_TS.m

Indrumar de laborator 38 function [T,S] = desc_ts(A); % Matricea sistemului se descompune in produs TS, % matricea T inferior triunghiulara, % matricea S superior triunghiulara n=length(A); disp(' Matricea initiala a sistemului este :') A pause for i=1:n S(1,i)=A(1,i); end for j=1:n T(j,1)=A(j,1)/A(1,1); end for k=2:n for j=2:n if j>=2 if(k<j) T(k,j)=0; elseif k==j T(k,k)=1; end if(k>j) Q=0; for p=1:j-1 Q=Q+T(k,p)*S(p,j); T(k,j)=(A(k,j)-Q)/S(j,j); end end end end for j=1:n if k>=2 if k>j S(k,j)=0; elseif k<=j R=0; for p=1:k-1 R=R+T(k,p)*S(p,j); end S(k,j)=A(k,j)-R;

Metode Numerice 39 end end end end disp(' Matricea inferior triunghiulara, T este :') T disp(' Matricea superior triunghiulara, S este :') S disp(' Verificare T*S=A :' ) T*S

Indrumar de laborator 40 LUCRAREA 8

Rezolvarea sistemelor de ecuaţii liniare prin metoda iterativă Jacobi

Se consideră sistemul de ecuaţii (5.1), pentru care avem îndeplinită condiţia:

n,1i,0a ii =≠ . (8.1)

Rezolvarea sistemului prin metoda Jacobi începe prin a pune în evidenţă în partea stângă a semnului egal a necunoscutei din ecuaţia i: ix

( )( )

( )⎪⎪⎩

⎪⎪⎨

⎧

−−−−=

−−−−=−−−=

−− nn1n1n,n22n11nnn

22nn232312122

11nn121211

axa...xaxabx................................

axa...xaxabxaxa...xabx

. (8.2)

Pentru uşurinţa scrierii se trece la forma matriceală, notând:

⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

⎟⎟⎟⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜⎜⎜⎜

⎝

⎛

−=

nn

n

22

211

1

n

2

1

nn

2n

nn

1n

22

n2

22

2111

n1

11

12

ab...

abab

u;

x...xx

x;

0...aa

aa

............aa...0

aa

aa...

aa0

A . (8.3)

Sistemul se scrie sub forma:

uAxx += . (8.4)

Metoda Jacobi constă în faptul că pune sistemul (8.4) într-o formă ce permite calculul iterativ:

(8.5) uAxx )k()1k( +=+

pornind de la o aproximaţie iniţială x(0) a soluţiei sistemului.

Metode Numerice 41 În practică, este mai utilă scrierea relaţiei de recurenţă referitoare la componente:

,....1,0k,n,1i,xaba1x

n

ij1j

)k(jiji

ii

)1k(i ==

⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜

⎝

⎛

−= ∑≠=

+ (8.6)

Metoda este convergentă cu condiţia ca toate valorile proprii ale matricei A să fie mai mici decât unitatea în modul. Această condiţie cere ca:

n,1i,1aan

ij1j ii

ij =<∑≠=

. (8.7)

Relaţia (8.7) reprezintă pentru metoda iterativă Jacobi o condiţie necesară de convergenţă. În matricea sistemului iniţial, trebuie ca fiecare element de pe diagonala principală, în modul, să fie mai mare decât suma modulelor celorlalte elemente din linia respectivă. Este posibil ca sistemul, în forma sa iniţială să nu respecte această condiţie necesară de convergenţă. Înainte de a începe aplicarea algoritmului metodei, este necesară o reordonare şi renumerotare a ecuaţiilor în sistem astfel încât condiţia de convergenţă să fie îndeplinită. Sistemele pentru care se pretează metoda Jacobi au matricea diagonal dominantă. Pentru metodele de rezolvare prin tehnici iterative este necesar să se impună un nivel de precizie ε > 0 în obţinerea soluţiei. Oprirea calculelor se poate face atunci când valorile pe componente în cadrul a două iteraţii succesive sunt “suficient de apropiate”, printr-un test de forma:

ε<−+ )k(i

)1k(ii

xxmax (8.8)

Pentru aceste metode, convergenţa nu depinde de alegerea soluţiei iniţiale . )0(x

Indrumar de laborator 42 Programul in MATLAB: L8_Jacobi_1.m function sol_sist = Jacobi_1(A,b,x0,max_er,max_it) % Sistemul este de forma A*x=b % x0 este aproximatia initiala a solutiei sistemului % max_er este valoarea maxima admisa pentru componentele solutiei % max_it este numarul maxim de iteratii impus [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n = length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end % Se verifica mai intai conditia de convergenta; cc=1; for i=1:n, S=0; for j=1:n, if i==j else S=S+abs(A(i,j)); end if S>abs(A(i,i)) cc=0; end end end if cc==0 disp('Sistemul nu indeplineste conditia de converg.!') end % Se construieste procedeul iterativ pentru calculul necunoscutelor; k=1;

Metode Numerice 43 for i=1:n, x(i,k)=x0(i); end while k < max_it+1 k = k+1; Pasul=k-1; Pasul format long max=0; for i=1:n, S=0; for j=1:n, if i<j S=S+A(i,j)*x(j,k-1); end if i>j S=S+A(i,j)*x(j,k-1); end end x(i,k)=(b(i)-S)/A(i,i); err(i)=abs(x(i,k)-x(i,k-1)); if err(i)>max max=err(i); end end eroarea=max; format short e eroarea if max < max_er disp('Limita erorii a fost atinsa !') disp(' ') disp('Solutia sistemului este :') format long x return end if k == max_it+1 disp('Numarul maxim de iteratii, atins') disp(' ') disp('Solutia partiala a sistemului este :') format long x return end end


Rezolvarea sistemelor de ecuaţii liniare prin metoda iterativă Gauss-Seidel

Această metodă de rezolvare este asemănătoare cu metoda iterativă Jacobi, ea putând fi aplicată doar sistemelor care îndeplinesc condiţia de convergenţă dată de relaţia (8.7). Etapele de calcul presupun evidenţierea în partea stângă a semnului egal a necunoscutei din ecuaţia i: ix

n,1i,0a;xaba1x ii

n

ij1j

jijiii

i =≠⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜

⎝

⎛

−= ∑≠=

(9.1)

În contrast cu metoda Jacobi, care nu foloseşte în calculul lui decât componente ale vectorului , relaţia (8.6), metoda

Gauss-Seidel foloseşte pentru determinarea diferitelor componente ale vectorului atât componente ale vectorului de la pasul anterior cât şi componente ale vectorului din pasul prezent , pe măsură ce acestea au fost determinate, accelerând astfel convergenţa procesului.

)1k(x + )k(x

)1k(x + )k(x)1k(x +

Astfel, pentru calculul lui , se folosesc componentele noi

, deja calculate în cadrul acestui pas, cât şi

componentele de la pasul anterior, necalculate încă în pasul prezent.

)1k(ix +

)1k(1i

)1k(2

)1k(1 x,...,x,x +

−++

)k(n

)k(2i

)k(1i x...,,x,x ++

În aceste condiţii, relaţia de recurenţă pentru metoda Gauss-Seidel este de forma:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛−−= ∑∑

+=

−

=

++n

1ij

)k(jij

1i

1j

)1k(jiji

ii

)1k(i xaxab

a1x . (9.2)

Se apreciază că, la acelaşi nivel de precizie impus în calculul soluţiei sistemului, metoda Gauss-Seidel este de două ori mai rapidă decât metoda Jacobi.

Metode Numerice 45 Pentru metodele de rezolvare prin tehnici iterative este necesar să se impună un nivel de precizie ε > 0 în obţinerea soluţiei. Oprirea calculelor se poate face atunci când valorile pe componente în cadrul a două iteraţii succesive sunt “suficient de apropiate”, prin testul:

ε<−+ )k(i

)1k(ii

xxmax (9.3)

Pentru aceste metode, convergenţa nu depinde de alegerea soluţiei iniţiale . Este totuşi importantă alegerea acestui vector în ceea ce priveşte viteza de convergenţă. Cu cât o mai mare parte dintre componentele vectorului iniţial se află mai aproape de soluţia exactă a sistemului, cu atât viteza de convergenţă este mai mare.

)0(x

)0(x

În practică, este bine să ţinem seama de datele problemei practice care a generat un asemenea sistem astfel încât să nu alegem chiar la întâmplare vectorul iniţial . )0(x Programul in MATLAB: L9_G_Seidel_1.m function sol_sist = G_Seidel_1(A,b,x0,max_er,max_it) % Sistemul este de forma A*x=b % x0 este aproximatia initiala a solutiei sistemului % max_er este valoarea maxima admisa pentru componentele vectorului solutie % max_it este numarul maxim de iteratii impus [m,p]=size(A); if m~=p disp('Matricea sistemului nu este patratica !'); return end n = length(b); if m~=n disp('Vectorul b are dimensiune gresita !'); return end % Se verifica mai intai conditia de convergenta;

Indrumar de laborator 46 cc=1; for i=1:n, S=0; for j=1:n, if i==j else S=S+abs(A(i,j)); end if S>abs(A(i,i)) cc=0; end end end if cc==0 disp('Sistemul nu indeplineste conditia de converg.!') end % Se construieste procedeul iterativ pentru calculul necunoscutelor; k=1; for i=1:n, x(i,k)=x0(i); end while k < max_it+1 k = k+1; Pasul=k-1; Pasul format long max=0; for i=1:n, S=0; for j=1:n, if i<j S=S+A(i,j)*x(j,k-1); end if i>j S=S+A(i,j)*x(j,k); end end x(i,k)=(b(i)-S)/A(i,i); err(i)=abs(x(i,k)-x(i,k-1)); if err(i)>max max=err(i);

Metode Numerice 47 end end eroarea=max; format short e eroarea if max < max_er disp('Limita erorii a fost atinsa !') disp(' ') disp('Solutia sistemului, pas cu pas, este :') format long x return end if k == max_it+1 disp('Numarul maxim de iteratii, atins') disp(' ') disp('Solutia partiala a sistemului este :') format long x return end end


Polinomul de interpolare Lagrange Considerăm că funcţia f este cunoscută pe baza unui ansamblu de valori într-o reţea de noduri nu neapărat echidistante: n10 x,...,x,x

n,0i),x(fy ii == . (10.1)

Se pune problema determinării unui polinom de grad n, de forma:

)x(y...)x(y)x(y)x(L nn1100n ϕ++ϕ+ϕ= (10.2)

care să coincidă cu funcţia f în nodurile de interpolare:

n,0i,y)x(f)x(L iiin === . (10.3)

Funcţiile n10 ,...,, ϕϕϕ sunt polinoame de grad n ce urmează a fi determinate, fiind valorile funcţiei în punctele

n10 y,...,y,y

n10 x,...,x,x .Relaţiile (10.3) pot fi puse sub forma:

n,0i,y)x(yn

0kiikk ==ϕ∑

=. (10.4)

Ansamblul de relaţii (10.4) poate fi satisfăcut dacă:

. (10.5) ⎩⎨⎧

=≠

=ϕki,1ki,0

)x( ik

Polinomul kϕ se anulează în toate punctele ki,xi ≠ , deci este de forma:

∏≠=

=−=ϕn

ki0i

ikk n,0k),xx(a)x( . (10.6)

Determinarea coeficientului ka se poate face din condiţia

Metode Numerice 49 1)x( kk =ϕ (10.7)

provenită din (10.5), rezultând pentru kϕ o expresie de forma:

∏≠=

=−−

=ϕn

ki0i ik

ik n,0k,

xxxx)x( . (10.8)

În aceste condiţii, polinomul Lagrange capătă forma:

∏∑≠== −

−=

n

ki0i ik

in

0kkn xx

xxy)x(L . (10.9)

Programul în MATLABL10_Lagrange.m function val_functie = Lagrange(xk,yk,x) % Aproximarea functiilor prin polinom Lagrange % xk este vectorul nodurilor de interpolare % yk este vectorul valorilor functiei in nodurile xk % Vectorii xk si yk trebuie sa aiba aceeasi lungime % x este punctul in care se face aproximarea functiei n=length(xk); n disp('Functia este cunoscuta prin valorile :') xk yk pause L=0; disp('Valoarea polinomului, adunand termen cu termen :') for k=1:n F(k)=1; for i=1:n if i==k else F(k)=F(k)*(x-xk(i))/(xk(k)-xk(i)); end end

Indrumar de laborator 50 L=L+F(k)*yk(k) pause end val_functie=L; disp('Valoarea functiei in punctul dat este :') disp(' ') val_functie


Polinomul de interpolare Newton cu diferenţe divizate

Se consideră o funcţie f ale cărei valori se cunosc într-un ansamblu de noduri n,0i,x i = , nu neapărat echidistante. Expresiile:

1nn

1nn

12

12

01

01xx

)x(f)x(f,...,xx

)x(f)x(f,xx

)x(f)x(f

−

−

−−

−−

−−

(11.1)

poartă numele de diferenţe divizate de ordinul întâi ale lui f. Aceste diferenţe se notează prin:

. (11.2) )x,x(f),...,x,x(f),x,x(f n1n2110 −

Diferenţele divizate de ordinul I se constituie ca fiind operatori liniari faţă de f, fiind simetrice în raport cu punctele . Se pot defini recursiv diferenţele divizate de ordinul II prin:

ji x,x

ik

jikjkji xx

)x,x(f)x,x(f)x,x,x(f

−

−= . (11.3)

Diferenţele divizate de ordinul k se scriu pornind de la diferenţele divizate de ordin (k-1) astfel:

1iki

1ki1ikiikii1i xx

)x,...,x(f)x,...,x(f)x,...,x,x(f−+

−+−++− −

−= . (11.4)

Pentru aproximarea funcţiei f printr-un polinom Newton cu diferenţe divizate se porneşte de la:

)xx)...(xx(c

...)xx)(xx(c)xx(cc)x(N

1n0n

102010n

−−−+++−−+−+=

(11.5)

Se impune coincidenţa polinomului cu funcţia în nodurile de interpolare:

n,0i),x(f)x(N iin == (11.6)

Indrumar de laborator 52 relaţii din care se pot determina coeficienţii n,0k,ck = şi anume:

. (11.7) ⎩⎨⎧

===

)x,...,x,x(fc;.........)x,x(fc;)x(fc

n10n

10100

Expresia finală a polinomului Newton cu diferenţe divizate este:

).xx)...(xx)(x,...,x,x(f

...)xx)(x,x(f)x(f)x(N

1n0n10

0100n

−−−+++−+=

(11.8)

Diferenţele divizate de diverse ordine pot fi uşor calculate dacă se construieşte un tabel de forma:

kx )x(f k )x,x(f ji )x,x,x(f pji

0x )x(f 0 )x,x(f 10

1x )x(f 1 )x,x,x(f 210 )x,x(f 21 )x,x,x,x(f 3210 (11.9)

2x )x(f 2 )x,x,x(f 321 )x,x(f 32

3x )x(f 3 De exemplu,

13

2132321 xx

)x,x(f)x,x(f)x,x,x(f

−−

= , (11.10)

relaţie ce se determină urmărind liniile îngroşate din (11.9). Avantajul utilizării acestui tip de polinom este dat de faptul că nu se porneşte de la început cu o anumită valoare a lui n ci se creşte treptat gradul polinomului până când se atinge precizia dorită. Atunci când se doreşte aproximarea valorii funcţiei f într-un punct oarecare t, prin creşterea gradului polinomului Newton nu întotdeauna se obţine o creştere semnificativă a preciziei de aproximare. Ar fi de dorit, pentru a obţine o precizie foarte bună, ca punctele

să se situeze destul de aproape de punctul t în care se face ,...x,x,x 210

Metode Numerice 53 aproximarea. Numai că, în general, nodurile sunt fixate dinainte, iar pe măsură ce folosim un număr mai mare de astfel de noduri (atunci când gradul polinomului creşte) ne depărtăm de punctul t, fapt care conduce la "pierderea" din precizie a aproximării.

,...x,x,x 210

Programul în MATLAB L11_Newton_DD.m function val_functie = Newton_DD(xk,yk,x) % Aproximare prin Polinoame Newton cu diferente divizate % xk este vectorul nodurilor de interpolare % yk este vectorul valorilor functiei in nodurile xk % Vectorii xk si yk trebuie sa aiba aceeasi lungime % x este punctul in care se doreste aproximarea n=length(xk); n1=length(yk); if n==n1 else disp('Functia este incorect data !') disp('Cei doi vectori nu au aceeasi lungime !') return end n disp('Functia este cunoscuta prin valorile :') xk yk pause % Calcul diferente divizate for i=0:n-2, DD(1,i+1)=(yk(i+2)-yk(i+1))/(xk(i+2)-xk(i+1)); end for k=2:n-1, m=n-k; for i=1:m, DD(k,i)=(DD(k-1,i+1)-DD(k-1,i))/(xk(i+k)-xk(i)); end end disp('Tabloul cu diferente divizate de diverse ordine:') disp(' ') DD' pause

Indrumar de laborator 54 % Calculul polinomului NW_DD disp('Polinoamele NW_DD, calculate pas cu pas, sint :') P=yk(1); NW_DD=P T=1; for k=1:n-1, pause T=T*(x-xk(k)); P=P+DD(k,1)*T; NW_DD=P end pause disp('Valoarea functiei in punctul indicat') disp('aproximata prin Polinom Newton cu') disp('diferente divizate este :') NW_DD


Polinoame de interpolare Newton cu pas constant

Dacă nodurile de interpolare sunt echidistante, se pot construi polinoame de interpolare Newton, de speţa I şi a II-a, folosind diferenţele finite.

310 x,...,x,x

Pasul h al unei astfel de reţele de noduri reprezintă diferenţa dintre două puncte consecutive:

1n,0i,xxh i1i −=−= + . (12.1)

Se consideră funcţia y = f(x) ale cărei valori sunt cunoscute într-un ansamblu de noduri echidistante, de pas h. Se numeşte diferenţă finită la dreapta pentru funcţia f(x) expresia:

)x(f)hx(f)x(f −+=∆ . (12.2)

Diferenţele finite de ordin superior se definesc recursiv astfel:

. (12.3) ))x(f())x(f()x(f 1n1nn ∆∆=∆∆=∆ −−

Aceste relaţii se simplifică dacă se calculează diferenţele finite în nodurile reţelei, deci pentru valorile kxx = :

∑=

−+

++

+

−=∆

+−=∆∆=∆

−=∆

n

0jjkn

jn

jk

n

k1k2kkk2

k1kk

yC)1(y

................yy2y)y(y

yyy

. (12.4)

În practică, este mai util a se calcula diferenţele finite la dreapta dacă se construieşte tabelul diagonal următor:


kx ky ky∆ k

2y∆ k

3y∆ ..... k

4y∆

0x 0y 0y∆

1x 1y 0

2y∆

1y∆ 0

3y∆

2x 2y 1

2y∆ 0

4y∆ ..... (12.5)

2y∆ 1

3y∆

3x 3y 2

2y∆ ........ ........

3y∆ ........

4x 4y ..... În mod asemănător se pot defini diferenţele finite la stânga printr-o relaţie de forma:

)hx(f)x(f)x(f −−=∇ . (12.6)

Şi aceste diferenţe se pot calcula recursiv prin relaţia:

. (12.7) ))x(f())x(f()x(f 1n1nn −− ∇∇=∇∇=∇

Dacă se efectuează calculul diferenţelor finite la stânga în nodurile reţelei se obţine:

∑=

−

−−

−

−=∇

+−=∇∇=∇

−=∇

n

0jjk

jn

jk

n

2k1kkkk2

1kkk

yC)1(y

............yy2y)y(y

yyy

. (12.8)

Tabelul diagonal, care este mai util în determinarea diferenţelor la stânga de diverse ordine, se construieşte pornind de la bază, el având aspectul:

Metode Numerice 57

kx ky ky∇ k

2y∇ k

3y∇ .....k

4y∇

...... ...... 4nx − 4ny − .........

3ny −∇ ...........

3nx − 3ny − 2n

2y −∇ ...........

2ny −∇ 0

3y∆ ...........

2nx − 2ny − 1n

2y −∇ n

4y∇ ..... (12.9)

1ny −∇ 1

3y∆

1nx − 1ny − n

2y∇

ny∇

nx ny În practică, pentru un ansamblu de noduri dat din punctul de vedere al calculelor efective aceste două tabele sunt identice. Diferenţa dintre cele două tabele este la nivel de semnificaţii şi anume un anumit element calculat al tabelului poate să aibă o anumită notaţie prin prisma primului tabel şi o altă notaţie prin prisma celui de al doilea.

n10 x,...,x,x ,

De exemplu,

3232 yyyy ∇=−=∆ , (12.10)

şi deci diferenţa (care este unică) este notată prin intermediul primului tabel şi

23 yy − 2y∆

3y∇ prin cel de-al doilea. a) Polinomul de interpolare Newton de speţa I Acest tip de polinom foloseşte diferenţele finite la dreapta. Se pune problema determinării unui polinom care să coincidă cu funcţia f în nodurile reţelei echidistante. Acest polinom se consideră de forma:

(12.11) )xx)...(xx)(xx(c...

)xx)(xx(c)xx(cc)x(N

1n10n

102010dn

−−−−+++−−+−+=

Indrumar de laborator 58 Determinarea coeficienţilor se face impunând condiţia: n10 c,...,c,c

n,0i,y)x(f)x(N iiidn === . (12.12)

În relaţia (12.11) se dau pe rând lui x valorile n10 x,...,x,x .Rezultă:

n0

n

n

20

2012

2

0011

00

h!ny

c

.........h!2y

hyy2y

c

h!1y

hyy

c

yc

∆=

∆=

+−=

∆=

−=

=

. (12.13)

În aceste condiţii, expresia polinomului devine:

).xx)...(xx)(xx(h!ny

...

)xx)(xx(h!2y

)xx(h!1

yy)x(N

1n10n0

n

1020

2

00

0dn

−−−−∆

++

+−−∆

+−∆

+= (12.14)

Coeficienţii se iau din tabelul (12.5) de pe latura de sus. 0n

00 y...y,y ∆∆ Expresia acestui polinom poate fi pusă sub formă mai comodă pentru aplicaţii care să nu depindă în mod explicit de nodurile reţelei de interpolare, dacă se notează:

hxx

t 0−= . (12.15)

Atunci:

n,1k)),1k(t)...(1t(th

)xx)...(xx)(xx(k

1k10 =−−−=−−− − . (12.16)

Rezultă pentru polinom expresia:

Metode Numerice 59

.y

!n))1n(t)...(1t(t

...y!2

)1t(tyty)thx(N

0n

02

000dn

∆−−−

+

++∆−

+∆+=+ (12.17)

Dacă se doreşte interpolarea într-un punct oarecare x diferit de noduri, se calculează t şi apoi se evaluează cu (12.17). Acest tip de polinom de interpolare se utilizează de obicei pentru aproximarea funcţiei f într-un punct x situat în vecinătatea punctului

)thx(N 0dn +

0x . b) Polinomul de interpolare Newton de speţa a II-a Acest tip de polinom utilizează diferenţele finite la stânga. Modul de construcţie al polinomului este oarecum asemănător. Se porneşte de la:

(12.18) ).xx)...(xx)(xx(c...

)xx)(xx(c)xx(cc)x(N

11nnn

1nn2n10sn

−−−+++−−+−+=

−

−

Se impune condiţia de coincidenţă a polinomului cu funcţia în nodurile de interpolare:

n,0i,y)x(f)x(N iiisn === . (12.19)

Se dau pe rând lui x valorile pentru determinarea coeficienţilor

01nn x,...,x,x −

n10 c,...,c,c :

nn

n

n

n1nn1

n0

h!ny

c

.........h!1

yhyy

c

yc

∇=

∇=

−=

=

−

. (12.20)

Polinomul are deci forma:


).xx)...(xx)(xx(h!ny

...

)xx)(xx(h!2y

)xx(h!1y

y)x(N

11nnnn

n

1nn2n

2

nn

nsn

−−−∇

++

+−−∇

+−∇

+=

−

−

(12.21)

Coeficienţii se iau din tabelul (12.9) de pe latura de jos.

nn

nn y,...,y,y ∇∇

Şi această expresie a polinomului poate fi pusă într-o formă care să nu depindă în mod explicit de nodurile reţelei de interpolare, dacă se notează:

hxxt n−

= . (12.22)

Se obţine:

)).1n(t)...(1t(t!ny

...)1t(t!2yt

!1yy)thx(N

nn

n2

nnn

sn

−++∇

+

+++∇

+∇

+=+ (12.23)

Acest polinom se utilizează de obicei pentru aproximarea funcţiei f într-un punct x situat în vecinătatea punctului . nx Programul în MATLAB L12_Newton_DF.m function val_functie = Newton_DF(xk,yk,x) % Aproximare prin Polinoame Newton cu diferente finite % xk este vectorul nodurilor de interpolare, echidistante % yk este vectorul valorilor functiei in nodurile xk % Vectorii xk si yk trebuie sa aiba aceeasi lungime % x este punctul in care se doreste aproximarea n=length(xk); n1=length(yk); if n==n1 else disp('Functia este incorect data !')

Metode Numerice 61 disp('Cei doi vectori nu au aceeasi lungime !') return end n disp('Functia este cunoscuta prin valorile :') xk yk pause % Se verifica mai intai echidistanta nodurilor a=1; for i=1:n-2, if xk(i+1)-xk(i)==xk(i+2)-xk(i+1) else a=0; end end if a==1 h=xk(2)-xk(1); else disp('Punctele date nu sunt echidistante.') disp('Corectati reteaua de noduri xk !') disp(' ') return end % Calcul diferente finite for i=0:n-2, Dyk(1,i+1)=yk(i+2)-yk(i+1); end for k=2:n-1, m=n-k; for i=1:m, Dyk(k,i)=Dyk(k-1,i+1)-Dyk(k-1,i); end end disp('Tabloul cu diferente finite de diverse ordine :') disp(' ') Dyk' pause % Calculul polinomului NW1 disp('Polinoamele NW1, calculate pas cu pas, sint:') P=yk(1);

Indrumar de laborator 62 NW1=P T=1; for k=1:n-1, pause T=T*(x-xk(k))/(k*h); P=P+Dyk(k,1)*T; NW1=P end disp('Valoarea functiei in punctul indicat') disp('aproximata prin Polinom Newton de speta I este :') NW1 pause % Calculul polinomului NW2 disp('Polinoamele NW2, calculate pas cu pas, sint:') R=yk(n); NW2=R S=1; for k=1:n-1, pause S=S*(x-xk(n-k+1))/(k*h); R=R+Dyk(k,n-k)*S; NW2=R end disp('Valoarea functiei in punctul indicat') disp('aproximata prin polinom Newton') disp('de speta a II-a este :') NW2


Aproximarea funcţiilor prin metoda celor mai mici pătrate. Dreapta de regresie

Se consideră funcţia f cunoscută într-un ansamblu de noduri

având valorile: n10 x,...,x,x

n,0i),x(fy ii == . (13.1)

Se pune problema aproximării acestei funcţii printr-un polinom de grad m, m < n, de forma:

(13.2) mm10m xc...xcc)x(F +++=

astfel încât abaterea medie pătratică şi deci suma:

(13.3) [∑=

−=n

0k

2kmk )x(F)x(fS ]

să aibă o valoare minimă. Cazul cel mai simplu îl constituie aproximarea funcţiei f printr-un polinom de gradul I reprezentând geometric o dreaptă de ecuaţie

bax)x(F1 += . (13.4)

Se impune condiţia ca suma:

(13.5) [∑=

+−=n

0k

2kk )bax(y)b,a(S ]

să aibă o valoare minimă. Condiţiile necesare de extrem se scriu:

[ ]

[ ]⎪⎪

⎩

⎪⎪

⎨

⎧

=+−−=∂

∂

=⋅+−−=∂

∂

∑

∑

=

=n

0kkk

n

0kkkk

0)bax(y2b

)b,a(S

0x)bax(y2a

)b,a(S

. (13.6)

Indrumar de laborator 64 Aceste relaţii reprezintă un sistem de 2 ecuaţii cu 2 necunoscute care mai poate fi scris sub forma:

⎪⎩

⎪⎨⎧

=+

=++

∑∑∑∑∑

kkk2k

kk

yxxbxa

y)1n(bxa. (13.7)

Pentru uşurinţa scrierii s-a renunţat la limitele de sumare. Rezolvarea sistemului permite obţinerea necunoscutelor a şi b de forma:

( )( )( )

( )( ) ( )(( )

)⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

−+

−=

+−

+−=

∑∑∑∑∑∑

∑∑∑∑∑

2k

2k

kkkk2k

2k

2k

kkkk

xx)1n(

yxxyxb

x)1n(x

yx)1n(yxa

. (13.8)

Se poate verifica uşor, folosind derivatele de ordinul al II-lea pentru S(a,b) că valoarea obţinută este o valoare minimă. Procedura descrisă poate fi generalizată căutând în locul funcţiei liniare o funcţie polinomială de grad m, cu m<n. În această situaţie, trebuie ca suma:

(13.9) [∑=

+++−=n

0k

2mkmk10k )xc...xcc(yS ]

să aibă o valoare minimă. Condiţiile necesare de extrem se scriu sub forma:

[ ] 0x)xc...xcc(y2cS n

0k

ik

mkmk10k

i=⋅+++−−=

∂∂ ∑

=. (13.10)

Dacă se notează:

Metode Numerice 65

(13.11)

⎪⎪

⎩

⎪⎪

⎨

⎧

=

=

∑

∑

=

=n

0kk

iki

n

0k

iki

yxv

xu

atunci sistemul (13.10) poate fi scris sub forma:

. (13.12)

⎪⎪⎩

⎪⎪⎨

⎧

=+++

=+++=+++

+

+

mm2m1m1m0

11mm2110

0mm1100

vuc...ucuc................................................

vuc...ucucvuc...ucuc

Soluţia acestui sistem reprezintă coeficienţii ai polinomului de aproximare căutat.

m0 c,...,c

Această metodă este utilizată frecvent în prelucrarea datelor experimentale. Atunci când se dispune de un număr mare de date (funcţia f cunoscută într-un număr mare de puncte) iar la o simplă reprezentare grafică aceste date se grupează sub forma unui nor având o anumită dezvoltare, se poate face aproximarea funcţiei prin această metodă. Nu se va putea găsi un polinom de aproximare care să treacă prin toate punctele norului (deci nu se poate folosi interpolarea) dar se va putea găsi o funcţie liniară (sau polinom de grad superior) care să urmărească direcţia norului de puncte îndeplinind condiţia dată de relaţia (13.3). Programul în MATLAB L13_regresie.m function drp_regresie = regresie(xk,yk) % Aproximarea functiilor prin Metoda celor mai mici patrate % xk este vectorul nodurilor in care se cunosc valorile functiei % yk este vectorul valorilor functiei in nodurile xk % Vectorii xk si yk trebuie sa aiba aceeasi lungime n=length(xk); n

Indrumar de laborator 66 disp('Functia este cunoscuta prin valorile :') xk yk pause A=sum(xk); B=sum(xk.^2); C=sum(yk); D=sum(xk.*yk); disp('Ecuatia dreptei construita prin') disp('aceasta metoda este y=a*x+b, unde:') a=(n*D-A*C)/(n*B-A*A) b=(C*B-D*A)/(n*B-A*A) % Reprezentarea grafica a dreptei de regresie si a punctelor disp('Reprezentarea grafica arata astfel : (vezi Fig.1)') clf i=1:n; x=min(xk)-10:0.1:max(xk)+10; plot(xk(i),yk(i),'b*',x,a*x+b,'r.') grid


Integrare numerică. Formulele trapezelor Aceste formule de cuadratură numerică se obţin prin particularizarea lui n în formulele lui Newton-Cotes şi folosind formule de tip închis şi deschis.

• n = 1, , formulă de tip deschis: 1=l

2

bax;2

abh 1+

=−

= . (14.1)

Singurul nod de calcul este punctul . Acest nod este un punct interior intervalului [a,b].

1x

Se calculează coeficientul Newton-Cotes:

(14.2) ).ab(I)ab(A;1I )1(111

)1(11 −=⋅−==

Rezultă:

⎟⎠⎞

⎜⎝⎛ +

⋅−==2

baf)ab()x(fAI~ 11 . (14.3)

Relaţia dată de (14.3) poartă numele de "a doua formulă a trapezelor". Restul acestei formule are expresia:

12

)ab(dx)xx(x!2

)(f)f(R3b

a1

''

1−

=−ξ

= ∫ (14.4)

• n=2, , formulă de tip închis: 0=l

bx;ax;)ab(h 21 ==−= . (14.5)

În acest caz avem două noduri, punctele şi , extremităţile intervalului [a,b].

1x 2x

Coeficienţii Newton-Cotes vor fi:


2abA;

21dt)1t(I

2abA;

21dt)2t(I

2

2

1

)2(20

1

2

1

)2(10

−==−=

−==−−=

∫

∫ (14.6)

Formula de cuadratură este:

))b(f)a(f(2

ab)x(fAI~2

1kkk +⋅

−== ∑

=. (14.7)

Această formulă, dată de (14.7), se numeşte "prima formulă a trapezelor". Restul formulei este, în acest caz, :

)(f12

)ab()bx)(ax(!2

)(f)f(R ''b

a

3''

2 ξ−

−=−−ξ

= ∫ . (14.8)

Aceste formule de cuadratură, aplicate o singură dată pentru întregul interval, conduc la erori mari de aproximare. Pentru creşterea preciziei se îmbunătăţeşte prima formula a trapezelor prin împărţirea intervalului [a,b] în n subintervale egale, prin punctele echidistante .bx,...,x,x,xa n210 == Se poate scrie:

. (14.9) ∑∑ ∫∫==

===−

n

1ii

n

1i

x

x

b

a

Idx)x(fdx)x(fIi

1i

Pe fiecare interval, integrala de tipul dată de (14.9) se aproximează prin:

iI

( ) n,1i,)x(f)x(f2xxI~ i1i

1iii =+⋅

−= −

− . (14.10)

Dacă pasul reţelei de noduri este:

ihxx,n

abh 0i +=−

= (14.11)

Metode Numerice 69 atunci:

⎟⎟⎠

⎞⎜⎜⎝

⎛+++⋅

−== ∑∑

−

==

1n

1i

n

1ii )b(f)iha(f2)a(f

n2abI~I~ . (14.12)

Pentru evaluarea restului, presupunând că există M>0 astfel încât )b,a(x)(,M)x(f '' ∈∀≤ atunci se poate scrie că:

2

3

2 n12)ab(M)f(R −

≤ . (14.13)

Programul în MATLAB : L14_trapez_int.m function integr = trapez_int(a,b,n,index_f) % Calculul integralei prin metoda trapezelor % Intervalul [a,b] se imparte in n subintervale format h=(b-a)/n int=0; n format long for i=0:n-1, int=int+h*(f(a+i*h,index_f)+f(a+h+i*h,index_f))/2; end disp('Valoarea integralei este :') int % Definirea functiei de integrat. function fct = f(x,index) switch index case 1 fct = exp(-x.^2); case 2 fct = 1 ./(1+x.^2); case 3 fct = 1 ./(2+sin(x)); case 4

Indrumar de laborator 70 fct = exp(cos(x)); case 5 fct = exp(x); end


Integrare numerică. Formula lui Simpson Acest tip de formulă se obţine particularizând formulele Newton-Cotes pentru n = 3 şi , de tip închis. Coeficienţii Newton-Cotes vor fi: 0=l

3)ab(2A;

32dt)3t)(1t(

21I

6abA;

61dt)3t)(2t(

41I

2

3

1

)3(20

1

3

1

)3(10

−==−−−=

−==−−=

∫

∫ (15.1)

6

abA;61dt)2t)(1t(

41I 3

3

1

)3(30

−==−−= ∫

Cele trei noduri sunt:

bx;2

bax;ax 321 =+

== (15.2)

deoarece

2

abh −= . (15.3)

Rezultă o valoare aproximativă a integralei de forma:

⎥⎦

⎤⎢⎣

⎡+⎟

⎠⎞

⎜⎝⎛ +

⋅+⋅−

== ∑=

)b(f2

baf4)a(f6

ab)x(fAI~3

1kkk . (15.4)

Şi această formulă de cuadratură introduce erori mari. Ea poate fi îmbunătăţită prin împărţirea intervalului [a,b] în 2n subintervale egale prin punctele bx,...,x,xa n10 == , de pas n2/)ab(h −= şi prin aplicarea formulei lui Simpson pe câte două subintervale consecutive, scriind integrala:

. (15.5) ∑∑ ∫∫==

===−

n

1ii

n

1i

x

x

b

a

Idx)x(fdx)x(fIi2

2i2

Indrumar de laborator 72 Folosind valoarea integralei dată de (15.4), se poate determina valoarea aproximativă pentru : iI

( ))i21i22i22i2i2

i x(f)x(f4)x(f6xx

I~ ++⋅−

= −−− . (15.6)

Rezultă:

⎟⎟⎠

⎞⎜⎜⎝

⎛+++⋅== ∑∑∑

−

==−

=)b(f)x(f2)x(f4)a(f

3hI~I~

1n

1ii2

n

1i1i2

n

1ii . (15.7)

unde n2,1i,ihxx 0i =+= . Pentru evaluarea restului se recurge în mod asemănător la formula trapezelor, determinând mai întâi restul pentru un interval de forma [ ] : i22i2 x,x −

∫−

−−−ξ

= −−

i2

2i2

x

xi21i22i2

i2)4(

i3 dx)xx)(xx)(xx(x!4

)(f)f(R . (15.8)

Rezultă:

)(f2880

)xx()f(R i2)4(

52i2i2

i3 ξ−

−= − . (15.9)

unde . )x,x( i22i2i2 −∈ξ

Dacă există M>0 astfel încât M)x(f )4( ≤ , )b,a(x)( ∈∀ atunci se poate

scrie restul pentru întreaga formulă:

Mn2880)ab()f(R 4

5

3−

≤ . (15.10)

Programul în MATLAB : L15_simpson_int.m function integr = simpson_int(a,b,n,index_f) % Calculul integralei prin metoda lui Simpson % Intervalul [a,b] se imparte in n subintervale

Metode Numerice 73 h=(b-a)/(2*n) int=0; n format long for i=0:n-1, int=int+h*(f(a+2*i*h,index_f)+4*f(a+h+2*i*h,index_f)+ +f(a+2*h+2*i*h,index_f))/3; end disp('Valoarea integralei este :') int % Definirea functiei de integrat. function fct = f(x,index) switch index case 1 fct = exp(-x.^2); case 2 fct = 1 ./(1+x.^2); case 3 fct = 1 ./(2+sin(x)); case 4 fct = exp(cos(x)); case 5 fct = exp(x); end


Integrare numerică. Formulele lui Gauss Se consideră integrala:

(16.1) ∫=b

a

dx)x(f)x(pI

a cărei valoare exactă nu poate fi determinată analitic. Formulele lui Gauss de cuadratură numerică sunt tot de forma:

, (16.2) ∑=

=≅m

1kkk )x(fAI~I

dar se caută să se determine atât nodurile kx , cât şi coeficienţii kA , n,1k = , astfel încât să se obţină un grad de exactitate cât mai mare pe

clasa polinoamelor algebrice. Dacă aproximarea funcţiei f se consideră un polinom construit pe baza şirului fundamental , vom avea: m2 x,...,x,x,1

m,0i,xAdx)x(f)x(pn

1k

ikk

b

a

==∑∫=

. (16.3)

Relaţia (16.3) reprezintă un sistem de (m+1) ecuaţii cu 2n necunoscute (nodurile kx şi coeficienţii kA , n,1k = ). Pentru a avea o soluţie unică trebuie să avem n21m =+ , adică 1n2m −= . Rezultă că gradul de exactitate al formulelor Gauss este 1n2m −= . Determinarea nodurilor kx şi a coeficienţilor kA , n,1k = prin rezolvarea sistemului (16.3) este dificilă, sistemul fiind neliniar. Se caută altă modalitate de determinare a necunoscutelor şi anume:

• Se consideră un polinom oarecare de grad , notat pentru care formula de cuadratură este exactă:

1n2 −)x(P 1n2 −

. (16.4) ∑∫=

−− =n

1kk1n2k

b

a1n2 )x(PAdxP)x(p

Metode Numerice 75

• Se presupune problema rezolvată în privinţa nodurilor şi se formează polinomul:

)xx)...(xx)(xx()x(W n21n −−−= . (16.5)

• Se face împărţirea cu rest a lui la : )x(P 1n2 − )x(Wn

)x()x(q)x(W)x(P n1n2 γ+=− (16.6)

cu grad q ≤ n - 1, grad γ < n. • pentru )x()x(Wn γ , de grad 1n2 − , formula de cuadratură este

exactă, adică:

. (16.7) 0)x()x(WAdx)x()x(W)x(pn

1kkknk

b

an =γ=γ ∑∫

=

deoarece n,1k,0)x(W kn == . Rezultă polinomul este ortogonal oricărui polinom γ(x) de grad (n-1), cu ponderea p(x).

)x(Wn

Se poate demonstra că pentru p(x) ≡ 1, polinoamele care

au coeficientul lui egal cu 1 şi care sunt ortogonale pe [a,b] oricărui polinom de grad cel mult (n-1) sunt polinoamele Legendre de forma:

)x(Wnnx

( ) ( )[ ]nnn

n

n bxaxdxd

)!n2(!n)x(W −−= . (16.8)

În plus, aceste polinoame au proprietatea că toate rădăcinile lor sunt reale şi cuprinse în intervalul (a,b). Polinoamele fiind determinate, sunt determinate şi nodurile

)x(Wn

n,1k,xk = , ca fiind rădăcinile acestui polinom. • Determinarea coeficienţilor n,1k,Ak = se poate face acum

rezolvând sistemul (16.3) care devine liniar. Mai comod, se poate proceda astfel: Se consideră polinomul:

k

nk xx

)x(W)x(Q−

= (16.9)

Indrumar de laborator 76 de grad (n – 1). Pătratul său, )x(Q2

k , de grad (2n – 2), se anulează în punctele , deoarece din expresia lui ki,xx i ≠= )x(Qk lipseşte paranteza

)xx( k− . Având gradul (2n-2), pentru acest polinom formulele de cuadratură sunt exacte, adică:

. (16.10) )x(QA)x(QAdx)x(Q k2kki

2k

n

1ii

b

a

2k ==∑∫

=

Integrând prin părţi, obţinem:

∫

∫∫

−+

+−

−=−

=

b

a k

nn

b

ak

2n

b

a2

k

2n

b

a

2k

dxxx

)x('W)x(W2

xx)x(W

dx)xx()x(W

dx)x(Q

. (16.11)

Primul termen poate fi calculat dacă în scrierea lui se foloseşte formula de derivare Leibniz:

)x(Wn

( )[ ] ( )[ ]nkn

knn

k

kn

0k

knn bx

dxdax

dxdC

)!n2(!n)x(W −−= −

−

=∑ (16.12)

Rezultă:

n

2

n

n2

nn

)ab()!n2()!n()b(W

)ab()!n2()!n()1()a(W

−=

−−= (16.13)

Ultima integrală din (16.11) este exactă deoarece polinomul de sub integrală este de grad (2n-2). Vom avea:

. (16.14) )x(W)x(QA2dx)x(W)x(Q2 i'nik

n

1ii

b

a

'nk ∑∫

==

Metode Numerice 77 Ţinând cont de expresia lui )x(Qk din (16.9) se poate scrie:

)x(W)x(Q k'nkk = . (16.15)

Se obţine:

[ ]( )[ ]

[ ]2k'nk

kk2

1n242k

'nk )x(WA2

)xb)(xa()!n2)ab()!n()x(WA +

−−

−=

+ (16.16)

De aici coeficienţii lui Gauss:

( )[ ] [ ]2k

'nkk

2

1n24

k)x(W)xb)(ax(!n2

)ab()!n(A−−

−=

+. (16.17)

De obicei, polinoamele Legendre şi rădăcinile lor se calculează pentru intervalul [-1,1]. Trecerea de la acest interval la intervalul [a,b] se poate face prin:

]1,1[t,2

)ab(t)ab(x −∈++−

= . (16.18)

Coeficienţii kA se pot determina independent de funcţie şi de limitele de integrare şi pot fi tabelaţi. Se poate demonstra că restul formulei de cuadratură Gauss este de forma:

( )[ ]

),(f)!1n2(!n2

)ab()!n()f(R )n2(2

1n24

n2 ξ+

−=

+ (16.19)

unde ξ ∈ (a,b). Aceste formule de cuadratură conduc la cele mai exacte rezultate, chiar şi pentru un număr mic de noduri. De multe ori, atunci când valorile lui f sunt dificil de calculat, formulele lui Gauss pentru un număr mic de noduri sunt de preferat altor formule de cuadratură. Pentru calculul rădăcinilor polinomului Legendre de diverse grade se poate folosi următorul program MATLAB:

Indrumar de laborator 78 Programul în MATLAB : L16_rad_legendre.m function set_pct = rad_legendre(n); % Functia permite determinarea radacinilor % polinomului Legendre de grad n folosit in % Formula de integrare Gauss a=[1 0 -1]; b=[1]; for i=1:n, b=conv(a,b); end % Aranj(n,k)=Comb(n,k)*Perm(k) % Noul coef b este vechiul coef b * Aranj(n,k) for i=1:n+1, b(i)=nchoosek((2*n-i+1),n)*factorial(n)*b(i); b=b(1:n+1); end format long rad_leg=roots(b); set_pct=sort(rad_leg); disp(' ') disp('Setul de puncte pe [-1, 1] este :') set_pct

Calculul integralei prin formule Gauss se face utilizând următorul program MATLAB:

Programul în MATLAB : L16_gauss_int.m function integr = gauss_int(a,b,n,index_f); % [int,baza,pond] = gauss_int(fun,a,b,n) integreaza % functia f intre a si b folosind metoda lui % Gauss in n puncte, care este o metoda exacta pentru % un polinom de grad 2*n-1. Conceput dupa modelul % din 'Methods of Numerical Integration', autori % Philip Davis si Philip Rabinowitz, dupa care au fost

Metode Numerice 79 % alese punctele de integrare si factorii % de ponderare. % a,b - limitele de intergrare % n - ordinul formulei lui Gauss (implicit n = 20) % int - valoarea integralei % baza, pond - sistemul de puncte si factorii de ponderare pe [-1,1] % index_f - indexul functiei de integrat. format long if nargin < 4, n=20; end u = (1:n-1)./sqrt((2*(1:n-1)).^2-1); [vc,baza] = eig(diag(u,-1)+diag(u,1)); [baza,k] = sort(diag(baza)); pond = 2*vc(1,k)'.^2; x = (a+b)/2+((b-a)/2)*baza; % Punctele in care se face calculul valorilor functiei: disp(' ') disp('Baza de noduri din intervalul [-1, 1] ') disp('pe care se bazeaza formula de calcul Gauss :') baza pause disp('Punctele reale in care se face calculul valorilor functiei :') x pause disp('Ponderile valorilor functiei f in aceste puncte :') pond pause f = fcn(x,index_f)*(b-a)/2; int = pond(:)'*f(:); disp('Valoarea integralei este: ') int % Definirea functiei de integrat function fct = fcn(x,index)

Indrumar de laborator 80 switch index case 1 fct = exp(-x.^2); case 2 fct = 1 ./(1+x.^2); case 3 fct = 1 ./(2+sin(x)); case 4 fct = exp(cos(x)); case 5 fct = exp(x); end

METODE NUMERICE - Indrumar de laborator

Documents

Transcript of METODE NUMERICE - Indrumar de laborator