Integrarea Numerica a Functiilor

14. INTEGRAREA NUMERIC A FUNCIILOR

Integrarea numeric este utilizat foarte des n calculele inginereti, fie pentru calculul unor integrale propriu-zise, fie pentru integrarea numeric a ecuaiilor difereniale.

Calculul numeric al unei integrale, denumit i cuadratur, se poate face prin aproximarea funciei de integrat printr-o alt funcie, pentru care integrala se cunoate sau se poate calcula uor, sau prin aproximarea funciei printr-un set de funcii liniare sau parabolice, pe poriuni.

14.1. Calculul cu aproximaie al integralelor definite

Metodele aproximative de calcul ale integralei definite au ca principiu nlocuirea curbei )(xfy = , n intervalul de integrare [a, b], cu o curb mai simpl

)(xgy = , deci :

)d )(d )(b

a

b

a

xxgxxf

unde curba )(xgy = , este o funcie n scar, o linie poligonal (metoda dreptunghiului, metoda trapezelor, metoda tangentelor), un lan de parabole, (metoda lui Simpson) sau polinomul de interpolare al lui Lagrange. Ideea fundamental, de obinere a unei soluii aproximative, este de a scrie integrala sub forma,

j

jj hwxfxxf )(d )(b

a

unde jw , sunt ponderi asociate cu procesul de integrare. Relaia este cunoscut ca formula prin cuadraturi i un procedeu posibil de generare este cu polinoame Lagrange.

14.1.1. Metoda dreptunghiurilor Pentru funcia de integrat, )(xfy = , pe intervalul de integrat [a, b], se pot lua

subintervale egale, n

abxx kk

=+1 , astfel nct

bxxxxxxxannkk

Integrarea numeric a funciilor 293

[ ])(...)()()( 1211 ++++ nxfxfxfafn

abI ,

[ ])()(...)(( 1212 bfxfxfxfn

abIn

++++

,

astfel c, dac funcia este cresctoare pe interval, atunci I1 aproximeaz prin lips, iar I2 prin exces.

n metoda dreptunghiurilor, se nlocuiete arcul )(xfy = , cu o funcie n trepte, iar eroarea de calcul ( )

n

abA2

, cu )('sup xfAbxa

= .

14.1.2. Metoda trapezelor

Metoda trapezelor const n a aproxima integrala definit =b

a

d )( xxfI , prin

semisuma valorilor I1 i I2, calculate anterior,

[ ])()(2...)(2)(2)(22 121

21 bfxfxfxfafn

abIIIntr +++++

=

+

care, dup prelucrri, se poate scrie :

++

=

1

1)(2)()(

2n

iitr xfbfaf

n

abI

n metoda trapezelor, se nlocuiete arcul de curb )(xfy = , cu o coard, iar eroarea de calcul ( )2

3

12 nabB

, cu )(sup '' xfBbxa

= .

14.1.3. Metoda tangentelor n aceast metod, se aproximeaz pe fiecare subinterval [ ]1, +kk xx , arcul de

curb )(xfy = , cu tangenta la curb, lundu-se numrul diviziunilor pare (n=2m) i se duce tangenta la curb n punctele de abscis x2k+1, tangent pe care o mrginim la dreptele paralele cu axa Oy, duse prin punctele (x2k , 0) i (x2k+1 , 0), astfel nct, aria trapezului este dat de ( ) )()()(

2 12222 ++

=+

kkk xfm

abxfxf

m

ab.

nsumnd pentru toate ariile elementare, obinem, ( ) ( ) ( )[ ]1231

b

a

tan ...d )( +++

= mxfxfxfm

abxxfI

14.1.4. Metoda lui Simpson n metoda lui Simpson se ia o medie ponderat a valorilor aproximative,

obinute prin metoda trapezelor i a tangentelor, anume 3

2 tanIII trSims+

= .

SISTEME DE PROGRAMARE PENTRU MODELARE I SIMULARE 294

Formula lui Simpson, lund numrul de diviziuni par, n=2m, este : ( ) ( ) ( ) ( ) ( ) ( ) ( )[ ]bfxfxfxfxfxfaf

m

abImmSims +++++++

1222321 42...4246

n metoda lui Simpson se aproximeaz curba n intervalul [ ]222 , +kk xx , printr-un arc de parabol CxBxAy ++= 2 , care trece prin punctele curbei de abscise

x2k , x2k+1 , x2k+2 , iar eroarea de calcul ( )

4

5

2880 mabC

, cu )(sup xfC IVbxa

= .

14.1.5. Aproximarea prin interpolare Considerm c intervalul de interpolare [a, b], se poate mpri n n

subintervale oarecare, bxxxxxa

nn=


( ) ( )( ) ( ) ( ) ( ) ( )

nkkkkkkk

kk tttttttttt

xfabY

=

+ KK 1110

,

( ) ( ) ( ) ( ) ( ) ( )nkkk tttttttttttT = + KK 1110 ,

n care ( )tTk , sunt polinoame de gradul n , n t. Pentru diverse valori ale lui n, s-au calculat integralele ( )=

1

0

' d ttTT kk i

coeficienii ( ) ( ) ( ) ( ) ( )nkkkkkkk

kk tttttttttt

TA

=

+ KK 1110

'

, astfel nct, prin

intermediul acestora, ce se gsesc calculai, aproximarea integralei se reduce la calculul sumei ( ) ( )

=

n

kkk xfAab

0.

14.1.6. Calculul numeric al integralelor pe intervale infinite

n definiia anterioar a integralei definite b

a

d )( xxf , am presupus c limitele a i b sunt finite, iar funcia f(x) este mrginit pe [a, b].

Dac unul sau amndou numerele a i b sunt infinite, integralele sunt de forma,

a

d )( xxf ,

b

-

d )( xxf , +

-

d )( xxf ,

care se pot reduce cu uurin la cazul

=

a

d )( xxfI . Pentru calculul integralelor de acest tip, presupunem c integrala are sens sau

c este convergent, dac

A

a

d )(lim xxfA

exist i este finit, pentru funcia f(x), definit pe intervalul [ ]+,a i integrabil pe intervalul [a, A], pentru orice A>0.

O integral care nu este convergent, se spune c este divergent sau c nu are sens.

Calculul integralelor convergente de tipul

=

a

d )( xxfI , se poate face : prin schimbarea variabilei, astfel nct intervalul infinit [ ]+,a , s se

transforme ntr-un interval finit, de exemplu, [ ]ba, , i apoi, se evalueaz prin una dintre metodele expuse anterior;

prin separarea integralei n dou integrale,

+=baa

d )(d )(d )( xxfxxfxxfb

, iar

apoi calculul se poate face: neglijnd termenul integral infinit, dac b este suficient de mare,

astfel nct >>bb

xxfxxf2

bad )(d )( ;


aproximnd forma asimptotic a lui f(x) , printr-o funcie )()( xfxg , care, pentru x suficient de mare, permite calculul analitic al integralei funciei g(x).

14.1.7. Calculul numeric al integralelor funciilor cu singulariti Singularitile sunt puncte ale unui interval mrginit, pentru care funcia este

infinit, ca de exemplu integrala

1

0 1d

x

x. Pentru aceste cazuri se recomand :

eliminarea singularitilor, atunci cnd este posibil, prin schimbarea de variabile, integrarea prin pri etc;

folosirea metodei Simpson, cu excluderea singularitilor, prin nlocuirea limitei 1 cu 1 , cu suficient de mic;

folosirea cuadraturii, bazate pe aproximarea funciei cu polinoame.

14.1.8. Integrarea dubl numeric n multe situaii, se ntlnesc integrale duble definite, avnd aplicabilitate la

calculul ariei unui domeniu plan, ariei unei suprafee din spaiu, volumului corpurilor, centre de greutate ale plcilor, momente de inerie.

O integral dubl este

===b a

DD

xyxfySyxfyxyxfIIa 0

d ),(d d),(dd),( unde f, este funcia de integrat, x i y sunt variabilele de integrare, D este domeniul de integrare, considerat dreptunghiular, [ ] [ ]bax 0,y , ,0 .

Vom considera pasul de integrare egal cu h, pe ambele subdomenii ale domeniului D. Considerm c I este numrul de subintervale spaiale pe intervalul [0, a] i J numrul de subintervale spaiale pe intervalul [0, b], care, evident, vor avea lrgimea, h

Jb

Iahh yx ==== .

Dac integrala dubl se calculeaz prin metoda trapezului, atunci :

( ) = =

=

=

++++ =+++I

i

J

jjiji

J

j

I

ijijijiji hwfffffhII

0 0

2,,

1

0

1

0,11,11,,

2

1 4

n care ponderile jiw , , asociate procesului de integrare, sunt 1/4, 2/4 i 4/4, dup cum nodul respectiv este punct de col, punct de margine sau punct interior.

Aplicnd formula lui Simpson, integrala dubl este :

( )( )

= =

=

=

+++++++

=

++++++++

I

i

J

jjiji

J

j

I

ijijijijijijijijiji

hwf

fffffffffhII

0 0

2,,

1

1

1

1,11,1,11,,111,11,11,11,1

2

2 1649

unde ponderile jiw , , iau acum valori diferite de cele anterioare.


14.2. Funcii Matlab pentru integrarea numeric

Funcia trapz, realizeaz calculul numeric al integralei prin metoda trapezelor. Se apeleaz cu una dintre sintaxele:

Z=trapz(Y) estimeaz aproximativ integrala lui Y pe intervalul unitate, iar pentru a calcula integrala pe alte intervale, se nmuleste intervalul cu Z; pentru vectori, Z este un vector ce conine integrala lui Y, iar dac Y este o matrice, Z este un vector coloan, ce integreaz fiecare coloan;

Z=trapz(X,Y) calculeaz integrala lui Y, innd cont de X; dac X este un vector coloan i Y un ir a crui prim dimensiune este length(X), trapz(X,Y) opereaz pe dimensiune;

Z=trapz(,dim) integreaz pe dimensiunea lui Y, specificat de scalarul dim; lungimea lui X trebuie s fie la fel cu size(Y,dim).

Functia trapz presupune c integrantul este dat prin valori numerice, n noduri echidistante ale intervalului de integrare. Dac funcia de integrat este descris sub forma analitic )(xfy = , n primul rnd, se creeaz cei doi vectori, care conin valorile perechilor abscis-ordonat (x-y).Acest lucru se realizeaz prin evaluarea funciei f(x), cu un pas adecvat pentru x i apoi se aplic funcia trapz.

Pentru exemplificare, se calculeaz numeric integrala x

dxx0

)sin( , care are valoarea exact 2. Secvena Matlab este :

X = 0:pi/100:pi; Y = sin(X); Z = pi/100*trapz(Y), rezultnd 1.9998. Pentru un interval neuniform avem: X = sort(rand(1,101)*pi); Y=sin(X); Z = trapz(X,Y), rezultnd 1.9981.

Pentru integrarea numeric cumulativ se utilizeaz funcia cumtrapz, care utilizeaz tot metoda trapezelor i se apeleaz cu una dintre sintaxele :

Z=cumtrapz(Y) estimeaz aproximativ integrala cumulativ a lui Y, prin metoda trapezelor, pe intervalul unitate i pentru a calcula integrala pe alte intervale, se nmulete intervalul cu Z; pentru vectori, Z este un vector ce conine integrale cumulative ale lui Y, iar pentru matrice, Z este o matrice de aceeai dimensiune cu Y, cu integralele cumulative calculate pe fiecare coloan;

Z=cumtrapz(X,Y) calculeaz integrale cumulative ale lui Y, respectiv ale lui X, folosind metoda trapezelor; X si Y trebuie s fie vectori de aceeai lungime sau X trebuie s fie un vector coloan i Y un ir de dimensiune length(X); dac X este un vector coloan i Y este un ir, a crui prim dimensiune este length(X), cumtrapz(X,Y) opereaz n aceast dimensiune;

Z= cumtrapz(X,Y,dim) sau Z=cumtrapz(Y,dim) integreaz n dimensiunea lui Y, specificat de dimensiunea scalar dim; lungimea lui X trebuie s fie aceeai cu size(Y,dim).


Pentru exemplificare, se consider Y = [0 1 2; 3 4 5], aplicnd Z=cumtrapz(Y,1) rezult,

Z=[ 0 0 0 ; 1.5000 2.5000 3.5000], iar aplicnd Z=cumtrapz(Y,2) rezult,

Z=[ 0 0.5000 2.0000 ; 0 3.5000 8.0000 ].

Dac integrantul este exprimat sub forma unei funcii analitice, se utilizeaz pentru integrare funciile Matlab quad i quadl sau mai vechi, quad8. Dac funcia de integrat este descris prin valori numerice, se aproximeaz mai nti acestea, printr-o funcie, utiliznd o metod de interpolare i apoi aceasta va fi scris n sintaxa MATLAB, ca fiier de tip M.

Funcia quad utilizeaz metoda Simpson aplicat recursiv adaptiv, avnd o eroare de 10-6, n versiunile mai noi i de 10-3 n versiunile mai vechi, unde, pentru mbuntirea toleranei, se utiliza quad8, care se bazeaz pe o cuadratur prin metoda Newton-Cotes, n 8 puncte. Acum, aceasta cheam practic funcia la apelare, quadl.

Mai nou, pentru o evaluare mai bun, cu eroare relativ de 10-6, se utilizeaz funcia quadl, ce are la baz metoda de cuadratur recursiv adaptiv Lobatto.

Funciile quad i quadl pot trata unele singularitai, care se gsesc la unul dintre capetele intervalului de integrare, ns nu rezolv singularitaile care se gsesc n interiorul intervalului. Nici una dintre funciile quad, quad8, quadl nu poate rezolva integrale cu singulariti, de genul :

1

0

1 dxx

. Acestea trebuie aduse la

forma unor funcii simple, fr singulariti.

Sintaxele de apelare a funciilor quad, nelegnd prin aceasta quad, quadl sau eventual quad8, sunt :

q = quad...('fun',a,b), returneaz rezultatul integrrii numerice a funciei fun, ntre limitele a i b; funcia fun poate fi dat sub forma unei expreii, funcie obiect, funcie de tip M sau vector i, n acest caz, returneaz un vector de valori;

q = quad...('fun',a,b,tol), itereaz pn cnd eroarea relativ este mai mic dect eroarea relativ tol, care n variantele noi este implicit 1e-6, iar n cele vechi 1e-3, iar dac tol este sub forma unui vector de dou dimensiuni, tol=[rel_tol abs_tol], atunci afieaz tolerana relativ, rel_tol i tolerana absolut, abs_tol;

q = quad...('fun',a,b,tol,trace), n plus, controleaz, prin trace, afiarea pe ecran a valorilor intermediare;

q = quad...('fun',a,b,tol,trace,p1,p2,...) adaug argumentele p1, p2 funciei fun(x,p1,p2), iar pentru a utiliza valorile implicite pentru tol i trace se folosete matrice vid, ca de exemplu, quad ...(('fun',a,b,[ ],[ ],p1,p2).

Pentru a preveni un ciclu infinit, att quad ct i quadl au limita de


recursivitate 10. Atingerea acestei limite duce la afiarea mesajului "Recursion level limit reached in quad. Singularity likely " i seteaz q=inf.

Funcia fun se poate da : direct ca o expresie, Q = quad('1./(x.^3+x.^2-x-2)',-1,2), rezultnd -0.2632; sub forma unei funcii obiect, F=inline('1./(x.^3-2*x-5)'); Q= quadl(F,0,2),

rezultnd Q= - 0.4605; sub forma unei funcii de tip M, definind anterior funcia function y = myf(x); y = 1./(x.^3-2*x-5); i apoi Q = quadl('myf',0,2);

Pentru rezolvarea numeric a integralelor duble se folosete funcia dblquad, care se poate apela cu una dintre sintaxele :

result=dblquad('fun',xmin,xmax,ymin,ymax), returneaz rezultatul integrrii duble a funciei fun(x,y), unde x, aparine intervalului [xmin xmax] i y, aparine intervalului [ymin ymax], funcia fun putnd fi precizat prin una dintre metodele de la funcia quad;

result=dblquad('fun',xmin,xmax,ymin,ymax,tol) specific tolerana tol, care are valoarea implicit 1e-6;

result=dblquad('fun',xmin,xmax,ymin,ymax,tol,metod) specific metoda de integrare, care n mod normal este quadl sau se poate face o funcie de integrare de ctre utilizator, avnd grij ca argumentele s fie transmise n aceeai ordine ca i la quadl;

result=dblquad('fun',xmin,xmax,ymin,ymax,tol,metod,p1,p2,...) adaug argumentele p1, p2 funciei fun(x,p1,p2), iar pentru a utiliza valorile implicite pentru tol i metod se folosete matrice vid, ca, de exemplu, dblquad ...(('fun',a,b,[ ],[ ],p1,p2..) sau

dblquad(fun,xmin,xmax,ymin,ymax,1.e-6,@quad,p1,p2,...).

Integrarea funciei [ ] [ ]pipipi ,0 ,2, ,cossin + yxyxxy se face utiliznd urmtoarea secven :

Q = dblquad(inline('y*sin(x)+x*cos(y)'), pi, 2*pi, 0, pi), rezultnd Q = -9.8696, x fiind un vector, iar y un scalar.

Integrarea Numerica a Functiilor

Documents

Transcript of Integrarea Numerica a Functiilor