MN Indr Final

Mitic TEMNEANU

METODE NUMERICE - Indrumar de laborator -

Iai, 2009

Metode Numerice Indrumar de laborator

PREFA Prezenta lucrare se adreseaz studenilor 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 cnd o anumit tem nu se poate finaliza prin utilizarea metodelor clasice, este necesar gsirea unei modaliti de calcul care s permit totui obinerea rezultatelor problemei n cauz. n aceast situaie intervine calculul numeric, adic implementarea unei metode numerice pe un sistem de calcul, utiliznd un anumit limbaj de programare. Pentru fixarea mai clar a metodei de calcul abordate, au fost readuse n atenia studenilor elementele de baz specifice fiecrei metode de calcul, cu evidenierea pailor teoretici ce trebuie parcuri pentru rezolvarea unei probleme printr-o anumit metod. Programele au fost elaborate n MATLAB, avnd n vedere faptul c acest mediu de programare este cunoscut studenilor, el fiind utilizat i n cadrul lucrrilor practice aferente altor discipline. Fiind destinate nelegerii metodelor de calcul i a modului efectiv de realizare a modulelor soft, fiecare dintre aceste programe este conceput a se desfura pas cu pas, studentul avnd 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 urmtor. Acest lucru este deosebit de util, deoarece impune studentului aprofundarea teoretic a metodei respective de calcul. Avnd 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 studenilor 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 ecuaiilor prin metoda biseciei (njumtirii intervalului) 3 L 2. Rezolvarea ecuaiilor prin metoda lui Newton (tangentei) 7 L 3. Rezolvarea ecuatiilor prin metoda aproximaiilor succesive 11 L 4. Rezolvarea ecuatiilor prin metoda secantei 15 L 5. Rezolvarea sistemelor de ecuaii liniare prin metoda lui Gauss 19 L 6. Rezolvarea sistemelor de ecuaii liniare prin metoda Gauss-Jordan 27 L 7. Descompunerea unei matrice n produs TS 35 L 8. Rezolvarea sistemelor de ecuaii liniare prin metoda iterativ Jacobi 40

Indrumar de laborator 2 L 9. Rezolvarea sistemelor de ecuaii liniare prin metoda iterativ Gauss-Seidel 44 L 10. Polinomul de interpolare Lagrange 48 L 11. Polinomul de interpolare Newton cu diferene divizate 51 L 12. Polinoame de interpolare Newton cu pas constant 55 L 13. Aproximarea funciilor prin metoda celor mai mici ptrate. 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 ecuaiilor prin metoda biseciei (njumtirii intervalului)

Considerm c ecuaia 0)x(f = are o singur rdcin n intervalul i c funcia f este continu pe acest interval. Aceast presupunere este valabil, n condiiile parcurgerii primei etape, aceea de separare a unei singure rdcini ntr-un anumit interval.

)b,a( 00

Fie eroarea admis pentru soluia ecuaiei. Din punct de vedere grafic, rezolvarea ecuaiei prin aceast metod, este ilustrat n Fig.1.1. Intervalul iniial se mparte n dou pri egale prin punctul: )b,a( 00

2

bac 000

+= (1.1) f(x)

x

A0

B0

a0a1

xR

c0

b1

b0

Fig.1.1.

Se efectueaz apoi produsul . Intervalul care conine n continuare soluia se noteaz . Situaiile care pot aprea sunt urmtoarele:

)c(f)a(f 00)b,a( 11

(1.2)

==>==

=== 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 ecuaiilor prin metoda lui Newton (tangentei)

Considerm c ecuaia 0)x(f = conine n intervalul o singur soluie . Grafic, rezolvarea ecuaiei prin aceast metod este ilustrat n Fig.2.1. De asemenea, considerm c pe acest interval derivatele f' i f" pstreaz semn constant, deci f este strict monoton i

)b,a( 00Rx

f(x)

x

A0

B0

x0=a0 x1

A1

xR b0

x2

Fig.2.1.

nu are punte de inflexiune. Ea presupune aproximarea soluiei exacte printr-un ir de valori , ... obinute prin intersecia tangentelor duse la graficul funciei f n punctele A

Rx

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

Punctul iniial se alege ca fiind una din extremitile intervalului i anume aceea care ndeplinete condiia:

0x)b,a( 00

(2.1) 0)x("f)x(f 00 > Aceast condiie ne asigur c intersecia tangentei la grafic cu axa Ox se va afla n interiorul intervalului iniial. Considerm punctul generic ))x(f,x(A kkk situat pe graficul funciei f. Ecuaia tangentei n acest punct la graficul lui f este:

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

Indrumar de laborator 8 Intersecia cu axa Ox este punctul 1kx + obinut 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)

Funcia se numete funcia de iterare a lui Newton. Cu noua notaie, relaia de recuren dat de (2.3) se poate pune sub forma:

,...2,1,0k),x(x k1k ==+ (2.5) Notnd Rkk xx = , eroarea de aproximare la pasul k, se obine:

.!2

)x(" R2k1k

= + (2.6) Aceast relaie arat c eroarea de aproximare evolueaz dup o lege ptratic. Metoda lui Newton are deci o vitez mare de convergen. Dezavantajul metodei const n faptul c n cadrul fiecrui pas este necesar calculul derivatei funciei 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 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 aproximaiilor succesive

Pentru rezolvarea ecuaiei 0)x(f = prin aceast metod, printr-un artificiu de calcul se aduce aceast ecuaie la forma echivalent:

)x(x = . (3.1) Dac este soluie pentru ecuaia Rx 0)x(f = , ea este soluie i pentru ecuaia (3.1), deci vom avea:

)x(x RR = . (3.2) Fie o aproximaie iniial pentru soluia ecuaiei. )b,a(x 000 Metoda const n aproximarea soluiei exacte printr-un ir de valori , ..., construit pe relaia iterativ

Rx

1x 2x

,...2,1,0k),x(x k1k ==+ (3.3) Considerm c funcia , derivabil pe intervalul satisface condiia:

)b,a( 00

)b,a(x)(1)x(' 00

Indrumar de laborator 12 Din punct de vedere grafic, rezolvarea ecuaiei (3.1) prin aceast metod se desfoar aa 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

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

Considerm ecuaia 0)x(f = pentru care am separat n intervalul o singur soluie. )b,a( 00Din punct de vedere grafic, rezolvarea ecuaei prin aceast metod este ilustrat n Fig.4.1.

A0

a0 x0

f(x)

xRx1

b0

B0 Fig.4.1.

Se unesc punctele i printr-o dreapt de ecuaie: 0A 0B

)ax(ab

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

00

000

= . (4.1) Intersecia 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 conine n continuare soluia .

0x )b,a( 00Rx

Pentru a determina care subinterval conine soluia se efectueaz produsul:

Rx

Indrumar de laborator 16

. (4.3)

==>== max_er & 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 ecuaii liniare prin metoda lui Gauss

Se consider un sistem de n ecuaii cu n necunoscute de forma:

. (5.1.)

=+++

=+++=+++

nnnn22n11n

2nn2222121

1nn1212111

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

bxa...xaxabxa...xaxa

Dac se fac notaiile:

(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 numr de pai, a matricei ptratice A a sistemului ntr-o matrice superior triunghiular. Se consider 0a11 . n caz contrar, se reordoneaz i renumeroteaz ecuaiile sistemului astfel nct aceast condiie s fie ndeplinit. Numrul de pai pentru un sistem de n ecuaii cu n necunoscute este (n-1). Pentru a putea urmri modul n care se efectueaz triunghiularizarea sistemului este necesar descrierea operaiilor care se efectueaz n cadrul fiecrui pas n parte. Pasul 1: Elementul 0a11 se numete pivot n cadrul acestui pas, iar ecuaia "1" care rmne nemodificat se numete ecuaie de pivotare. Necunoscuta este eliminat din ultimele (n-1) ecuaii prin 1x


nmulirea ecuaiei 1 a sistemului cu rapoartele 11

1iaa i scderea acesteia

din ecuaia "i", cu n,2i = . Pasul k, 1n,2k = : Primele k ecuaii ale sistemului obinut la pasul anterior rmn nemodificate. Ecuaia "k" este ecuaia de pivotare iar elementul situat pe diagonala principal n aceast ecuaie este pivotul. Dac acest element este nul, este necesar reordonarea i renumerotarea ultimelor n-(k-1) ecuaii ale sistemului de la pasul precedent astfel nct s avem un pivot nenul, dup care se procedeaz la transcrierea primelor k ecuaii. n cadrul acestui pas, necunoscuta kx este eliminat din ultimele (n-k) ecuaii ntr-o manier asemntoare cu cea de la pasul 1. De aceast dat, ecuaia de pivotare se nmulete cu un raport care are la numrtor coeficientul lui kx din ecuaia curent iar la numitor pivotul din cadrul acestui pas i se scade din ecuaia curent. n acest fel necunoscuta kx este eliminat din ultimele (n k) ecuaii. La ultimul pas, (n-1), necunoscuta este eliminat numai din ultima ecuaie a sistemului obinut la pasul anterior, primele (n-1) ecuaii rmnnd nemodificate.

1nx

Pentru a putea scrie relaiile de recuren care leag matricele sistemului la trecerea de la un pas la cellalt, convenim s atam matricei iniiale 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. Notm:

(5.5) )k()k( bxA =forma sistemului nainte de eliminarea necunoscutei kx din ultimele (n-k) ecuaii. 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 funcie 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)

Relaiile (5.7) i (5.8) reprezint de fapt rezultatul nmulirii n cadrul pasului (k-1) a ecuaiei (k-1) a sistemului:

(5.9) )1k()1k( bxA =cu raportul )1k( 1k,1k

)1k(1k,i aa

, i scznd rezultatul din ecuaia i, pentru

i k. Este eliminat astfel necunoscuta xk-1 din ultimele n-(k-1) ecuaii. La sfritul 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 ecuaie, obinndu-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

Gsirea efectiv a soluiei 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 obine soluia 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 creterea preciziei calculelor numerice este necesar reordonarea ecuaiilor n sistem n cadrul fiecrui pas, astfel nct, de fiecare dat, pivotul s aib o valoare maxim, n modul. Reordonarea ecuaiilor se face prin nlocuirea de tip pereche a ecuaiei de pivotare ce rezult n mod normal la finele unui anumit pas cu o ecuaie situat sub aceasta i care are, n modul, cea mai mare valoare a primului coeficient nenul din ecuaie. Acest coeficient devine noul pivot iar ntreaga ecuaie va nlocui ecuaia de pivotare obinut 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 ecuaii liniare prin metoda Gauss-Jordan

Aceast metod permite rezolvarea sistemului (5.1) prin transformarea sa ntr-un sistem echivalent a crui matrice este diagonal. Rezolvarea presupune parcurgerea a n pai. n cadrul pasului generic kP se transcrie numai ecuaia avnd indicele k din sistemul de la pasul anterior. Necunoscuta kx este eliminat din toate celelalte (n-1) ecuaii, printr-un procedeu asemntor 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) rm

ne

nem

odifi

cat

Coloana (k-1) va avea zerouri n

afara elementului )1k( 1k,1ka

Coloanemodificate

Primele (k-2)coloane rmnnemodificate

(6.1)

Matricea sistemului va purta la finele fiecrui pas un indice superior egal cu numrul pasului care urmeaz. n cadrul fiecrui pas toate elementele matricei sistemului i indexeaz indicele superior cu o unitate.

Indrumar de laborator 28 Relaiile de recuren ce permit trecerea de la matricea A(k-1) la matricea A(k) pot fi scrise mai uor dac se evideniaz, pe elementele matricei A(k-1) operaiile care se execut, aa 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 obine un sistem diagonal de forma:

. (6.4)

=

==

++

++++

)1n(nn

)1n(nn

)1n(22

)1n(22

)1n(11

)1n(11

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

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

Soluia sistemului este evident dat de relaiile:

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 urmtor: - se formeaz o matrice de dimensiuni n(2n+1) compus din:

pe primele n coloane matricea A;

Metode Numerice 29

pe urmtoarele n coloane matricea ; nI pe ultima coloan vectorul b;

- se aplic procedeul de rezolvare dat de metoda Gauss-Jordan numai c operaiile descrise pentru matricea sistemului se extind asupra ntregii matrice de dimensiuni n(2n+1). Se efectueaz aceste operaii pn cnd n locul matricei A apare matricea diagonal. Ultima operaie care se efectueaz asupra acestei matrice extinse este aceea de a mpri 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 obine matricea unitate . Pe urmtoarele n coloane vom obine matricea invers A

nI-1 iar pe ultima coloan, a (2n+1) - a, vom obine chiar

soluia sistemului. Sintetic, metoda poate fi descris astfel:

[ ] [ ]x|A|Ib|I|A 1nn (6.6) x fiind soluia 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 creterea preciziei calculelor numerice este necesar reordonarea ecuaiilor n sistem n cadrul fiecrui pas, astfel nct, 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 acelai ordin i prin aplicarea algoritmului de rezolvare a sistemului prin metoda Gauss-Jordan, fr schimbarea pivotului. Odat cu soluia 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, introducnd n acelai timp un algoritm de calcul pentru elementele celor dou matrice n funcie 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

>==

==

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

. (7.5) =

=k

1ppkipik sta

Separnd elementul ikt din (7.5) se obine:

)ki(,stas1t

1k

1ppkipik

kkik >

=

=. (7.6)

Relaiile (7.4) i (7.6) permit determinarea tuturor elementelor matricelor T i S presupunnd 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 urmtoarea: o linie din S, o coloan din T, o linie din S, o coloan din T,..., pn la epuizarea elementelor celor dou matrice. Avnd determinate elementele matricelor produs se poate proceda la rezolvarea sistemului:

bTSx,bAx == . (7.9) Se noteaz:

ySx = , (7.10) astfel nct rezolvarea sistemului iniial (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 relaiile:

=

=

=

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

t

ytb

y

tby

kk

1k

1ppkpk

k

11

11

. (7.13)

Avnd 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

Soluia se afl prin eliminarea invers:

=

=

=

+= 1,...,2nk;

s

xsy

x

syx

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, corespondena ntre diversele notaii 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(kj) 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

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 ecuaii liniare prin metoda iterativ Jacobi

Se consider sistemul de ecuaii (5.1), pentru care avem ndeplinit condiia:

n,1i,0a ii = . (8.1) Rezolvarea sistemului prin metoda Jacobi ncepe prin a pune n eviden n partea stng a semnului egal a necunoscutei din ecuaia i: ix

( )( )( )

=

==

nn1n1n,n22n11nnn

22nn232312122

11nn121211

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

axa...xaxabxaxa...xabx

. (8.2)

Pentru uurina scrierii se trece la forma matriceal, notnd:

=

=

=

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 aproximaie iniial x(0) a soluiei sistemului.

Metode Numerice 41 n practic, este mai util scrierea relaiei de recuren referitoare la componente:

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

n

ij1j

)k(jiji

ii

)1k(i ==

=

=+ (8.6)

Metoda este convergent cu condiia ca toate valorile proprii ale matricei A s fie mai mici dect unitatea n modul. Aceast condiie cere ca:

n,1i,1aan

ij1j ii

ij = 0 n obinerea soluiei. Oprirea calculelor se poate face atunci cnd valorile pe componente n cadrul a dou iteraii succesive sunt suficient de apropiate, printr-un test de forma:

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 ij 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 ecuaii liniare prin metoda iterativ Gauss-Seidel

Aceast metod de rezolvare este asemntoare cu metoda iterativ Jacobi, ea putnd fi aplicat doar sistemelor care ndeplinesc condiia de convergen dat de relaia (8.7). Etapele de calcul presupun evidenierea n partea stng a semnului egal a necunoscutei din ecuaia i: ix

n,1i,0a;xaba1x ii

n

ij1j

jijiii

i =

=

= (9.1)

n contrast cu metoda Jacobi, care nu folosete n calculul lui dect componente ale vectorului , relaia (8.6), metoda

Gauss-Seidel folosete pentru determinarea diferitelor componente ale vectorului att componente ale vectorului de la pasul anterior ct i componente ale vectorului din pasul prezent , pe msur ce acestea au fost determinate, accelernd astfel convergena 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, ct 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 condiii, relaia de recuren pentru metoda Gauss-Seidel este de forma:

=

+=

=++ n

1ij

)k(jij

1i

1j

)1k(jiji

ii

)1k(i xaxaba

1x . (9.2)

Se apreciaz c, la acelai nivel de precizie impus n calculul soluiei sistemului, metoda Gauss-Seidel este de dou ori mai rapid dect metoda Jacobi.

Metode Numerice 45 Pentru metodele de rezolvare prin tehnici iterative este necesar s se impun un nivel de precizie > 0 n obinerea soluiei. Oprirea calculelor se poate face atunci cnd valorile pe componente n cadrul a dou iteraii succesive sunt suficient de apropiate, prin testul:

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 ij 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 Considerm c funcia f este cunoscut pe baza unui ansamblu de valori ntr-o reea de noduri nu neaprat echidistante: n10 x,...,x,x

n,0i),x(fy ii == . (10.1) Se pune problema determinrii unui polinom de grad n, de forma:

)x(y...)x(y)x(y)x(L nn1100n +++= (10.2) care s coincid cu funcia f n nodurile de interpolare:

n,0i,y)x(f)x(L iiin === . (10.3) Funciile n10 ,...,, sunt polinoame de grad n ce urmeaz a fi determinate, fiind valorile funciei n punctele

n10 y,...,y,y

n10 x,...,x,x .Relaiile (10.3) pot fi puse sub forma:

n,0i,y)x(yn

0kiikk ==

=. (10.4)

Ansamblul de relaii (10.4) poate fi satisfcut 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 condiia

Metode Numerice 49 1)x( kk = (10.7) provenit din (10.5), rezultnd pentru k o expresie de forma:

=

==

n

ki0i ik

ik n,0k,xx

xx)x( . (10.8)

n aceste condiii, polinomul Lagrange capt 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 diferene divizate

Se consider o funcie f ale crei valori se cunosc ntr-un ansamblu de noduri n,0i,x i = , nu neaprat 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 diferene divizate de ordinul nti ale lui f. Aceste diferene se noteaz prin:

. (11.2) )x,x(f),...,x,x(f),x,x(f n1n2110 Diferenele divizate de ordinul I se constituie ca fiind operatori liniari fa de f, fiind simetrice n raport cu punctele . Se pot defini recursiv diferenele divizate de ordinul II prin:

ji x,x

ik

jikjkji xx

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

= . (11.3)

Diferenele divizate de ordinul k se scriu pornind de la diferenele divizate de ordin (k-1) astfel:

1iki

1ki1ikiikii1i xx

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

+++ = . (11.4)

Pentru aproximarea funciei f printr-un polinom Newton cu diferene divizate se pornete de la:

)xx)...(xx(c

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

1n0n

102010n

+++++=

(11.5)

Se impune coincidena polinomului cu funcia n nodurile de interpolare:

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

Indrumar de laborator 52 relaii din care se pot determina coeficienii n,0k,ck = i anume:

. (11.7)

===

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

n10n

10100

Expresia final a polinomului Newton cu diferene divizate este:

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

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

1n0n10

0100n

++++=

(11.8)

Diferenele divizate de diverse ordine pot fi uor calculate dac se construiete 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)

relaie ce se determin urmrind liniile ngroate din (11.9). Avantajul utilizrii acestui tip de polinom este dat de faptul c nu se pornete de la nceput cu o anumit valoare a lui n ci se crete treptat gradul polinomului pn cnd se atinge precizia dorit. Atunci cnd se dorete aproximarea valorii funciei f ntr-un punct oarecare t, prin creterea gradului polinomului Newton nu ntotdeauna se obine o cretere semnificativ a preciziei de aproximare. Ar fi de dorit, pentru a obine 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 msur ce folosim un numr mai mare de astfel de noduri (atunci cnd gradul polinomului crete) ne deprtm de punctul t, fapt care conduce la "pierderea" din precizie a aproximrii.

,...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 spea I i a II-a, folosind diferenele finite.

310 x,...,x,x

Pasul h al unei astfel de reele de noduri reprezint diferena dintre dou puncte consecutive:

1n,0i,xxh i1i == + . (12.1) Se consider funcia y = f(x) ale crei valori sunt cunoscute ntr-un ansamblu de noduri echidistante, de pas h. Se numete diferen finit la dreapta pentru funcia f(x) expresia:

)x(f)hx(f)x(f += . (12.2) Diferenele finite de ordin superior se definesc recursiv astfel:

. (12.3) ))x(f())x(f()x(f 1n1nn == Aceste relaii se simplific dac se calculeaz diferenele finite n nodurile reelei, 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 diferenele finite la dreapta dac se construiete tabelul diagonal urmtor:


kx ky ky k2y k3y ..... k4y 0x 0y

0y 1x 1y 0

2y 1y 03y

2x 2y 12y 04y ..... (12.5)

2y 13y 3x 3y 2

2y ........ ........ 3y ........

4x 4y ..... n mod asemntor se pot defini diferenele finite la stnga printr-o relaie de forma:

)hx(f)x(f)x(f = . (12.6) i aceste diferene se pot calcula recursiv prin relaia:

. (12.7) ))x(f())x(f()x(f 1n1nn ==Dac se efectueaz calculul diferenelor finite la stnga n nodurile reelei se obine:

=

=

+===

n

0jjk

jn

jk

n

2k1kkkk2

1kkk

yC)1(y

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

yyy

. (12.8)

Tabelul diagonal, care este mai util n determinarea diferenelor la stnga de diverse ordine, se construiete pornind de la baz, el avnd aspectul:

Metode Numerice 57

kx ky ky k2y k3y ..... k4y ...... ......

4nx 4ny ......... 3ny ...........

3nx 3ny 2n2y ........... 2ny 03y ...........

2nx 2ny 1n2y n4y ..... (12.9) 1ny 13y

1nx 1ny n2y ny

nx ny n practic, pentru un ansamblu de noduri dat din punctul de vedere al calculelor efective aceste dou tabele sunt identice. Diferena dintre cele dou tabele este la nivel de semnificaii i anume un anumit element calculat al tabelului poate s aib o anumit notaie prin prisma primului tabel i o alt notaie prin prisma celui de al doilea.

n10 x,...,x,x ,

De exemplu,

3232 yyyy == , (12.10) i deci diferena (care este unic) este notat prin intermediul primului tabel i

23 yy 2y3y prin cel de-al doilea.

a) Polinomul de interpolare Newton de spea I Acest tip de polinom folosete diferenele finite la dreapta. Se pune problema determinrii unui polinom care s coincid cu funcia f n nodurile reelei 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 coeficienilor se face impunnd condiia: n10 c,...,c,c

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

n relaia (12.11) se dau pe rnd 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 condiii, 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)

Coeficienii 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 aplicaii care s nu depind n mod explicit de nodurile reelei 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 dorete 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 funciei f ntr-un punct x situat n vecintatea punctului

)thx(N 0dn +

0x . b) Polinomul de interpolare Newton de spea a II-a Acest tip de polinom utilizeaz diferenele finite la stnga. Modul de construcie al polinomului este oarecum asemntor. Se pornete de la:

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

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

11nnn

1nn2n10sn

+++++=

Se impune condiia de coinciden a polinomului cu funcia n nodurile de interpolare:

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

Se dau pe rnd lui x valorile pentru determinarea coeficienilor

01nn x,...,x,x n10 c,...,c,c :

nn

n

n

n1nn1

n0

h!nyc

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

yhyyc

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)

Coeficienii 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 reelei de interpolare, dac se noteaz:

hxxt n= . (12.22)

Se obine:

)).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 funciei f ntr-un punct x situat n vecintatea 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 funciilor prin metoda celor mai mici ptrate. Dreapta de regresie

Se consider funcia f cunoscut ntr-un ansamblu de noduri

avnd valorile: n10 x,...,x,x

n,0i),x(fy ii == . (13.1) Se pune problema aproximrii acestei funcii printr-un polinom de grad m, m < n, de forma:

(13.2) mm10m xc...xcc)x(F +++=astfel nct abaterea medie ptratic i deci suma:

(13.3) [=

=n

0k

2kmk )x(F)x(fS ]

s aib o valoare minim. Cazul cel mai simplu l constituie aproximarea funciei f printr-un polinom de gradul I reprezentnd geometric o dreapt de ecuaie

bax)x(F1 += . (13.4) Se impune condiia ca suma:

(13.5) [=

+=n

0k

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

s aib o valoare minim. Condiiile 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 relaii reprezint un sistem de 2 ecuaii cu 2 necunoscute care mai poate fi scris sub forma:

=+=++

kkk2k

kk

yxxbxa

y)1n(bxa. (13.7)

Pentru uurina scrierii s-a renunat la limitele de sumare. Rezolvarea sistemului permite obinerea 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 uor, folosind derivatele de ordinul al II-lea pentru S(a,b) c valoarea obinut este o valoare minim. Procedura descris poate fi generalizat cutnd n locul funciei liniare o funcie polinomial de grad m, cu m

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

Soluia acestui sistem reprezint coeficienii ai polinomului de aproximare cutat.

m0 c,...,c

Aceast metod este utilizat frecvent n prelucrarea datelor experimentale. Atunci cnd se dispune de un numr mare de date (funcia f cunoscut ntr-un numr mare de puncte) iar la o simpl reprezentare grafic aceste date se grupeaz sub forma unui nor avnd o anumit dezvoltare, se poate face aproximarea funciei prin aceast metod. Nu se va putea gsi un polinom de aproximare care s treac prin toate punctele norului (deci nu se poate folosi interpolarea) dar se va putea gsi o funcie liniar (sau polinom de grad superior) care s urmreasc direcia norului de puncte ndeplinind condiia dat de relaia (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 obin 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

2bax;

2abh 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)

Relaia 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 , extremitile intervalului [a,b].

1x 2x

Coeficienii 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 numete "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 creterea preciziei se mbuntete prima formula a trapezelor prin mprirea 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~ i1i1iii =+= . (14.10)

Dac pasul reelei de noduri este:

ihxx,n

abh 0i +== (14.11)

Metode Numerice 69 atunci:

+++==

==

1n

1i

n

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

abI~I~ . (14.12)

Pentru evaluarea restului, presupunnd c exist M>0 astfel nct )b,a(x)(,M)x(f '' atunci se poate scrie c:

23

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 obine particulariznd formulele Newton-Cotes pentru n = 3 i , de tip nchis. Coeficienii 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 mbuntit prin mprirea 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 cte 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

( ))i21i22i22i2i2i x(f)x(f4)x(f6xxI~ ++= . (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 asemntor la formula trapezelor, determinnd mai nti 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 nct M)x(f )4( , )b,a(x)( atunci se poate scrie restul pentru ntreaga formul:

Mn2880)ab()f(R 45

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 crei 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 att nodurile kx , ct i coeficienii kA , n,1k = , astfel nct s se obin un grad de exactitate ct mai mare pe

clasa polinoamelor algebrice. Dac aproximarea funciei 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)

Relaia (16.3) reprezint un sistem de (m+1) ecuaii cu 2n necunoscute (nodurile kx i coeficienii kA , n,1k = ). Pentru a avea o soluie unic trebuie s avem n21m =+ , adic 1n2m = . Rezult c gradul de exactitate al formulelor Gauss este 1n2m = . Determinarea nodurilor kx i a coeficienilor 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 privina nodurilor i se formeaz polinomul:

)xx)...(xx)(xx()x(W n21n = . (16.5) Se face mprirea 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 oricrui 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] oricrui polinom de grad cel mult (n-1) sunt polinoamele Legendre de forma:

)x(Wnnx

( ) ( )[ ]nnnnn bxaxdxd)!n2( !n)x(W = . (16.8) n plus, aceste polinoame au proprietatea c toate rdcinile lor sunt reale i cuprinse n intervalul (a,b). Polinoamele fiind determinate, sunt determinate i nodurile

)x(Wnn,1k,xk = , ca fiind rdcinile acestui polinom.

Determinarea coeficienilor n,1k,Ak = se poate face acum rezolvnd 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). Ptratul su, )x(Q2k , de grad (2n 2), se anuleaz n punctele

, deoarece din expresia lui ki,xx i = )x(Qk lipsete paranteza )xx( k .

Avnd 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 ==

=Integrnd prin pri, obinem:

+

+==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 folosete formula de derivare Leibniz:

)x(Wn

( )[ ] ( )[ ]nkn knnkkn0k

knn bxdx

daxdxdC

)!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 innd cont de expresia lui )x(Qk din (16.9) se poate scrie:

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

Se obine:

[ ] ( )[ ] [ ]2k'nkkk21n242

k'nk )x(WA2

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

=+

(16.16)

De aici coeficienii lui Gauss:

( )[ ] [ ]2k'nkk21n24

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

)ab()!n(A=

+. (16.17)

De obicei, polinoamele Legendre i rdcinile 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)

Coeficienii kA se pot determina independent de funcie i de limitele de integrare i pot fi tabelai. 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 numr mic de noduri. De multe ori, atunci cnd valorile lui f sunt dificil de calculat, formulele lui Gauss pentru un numr mic de noduri sunt de preferat altor formule de cuadratur. Pentru calculul rdcinilor polinomului Legendre de diverse grade se poate folosi urmtorul 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 utiliznd urmtorul 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

MN Indr Final

Documents

Transcript of MN Indr Final