INTEGRAREA NUMERIC Ă - Math · Acest polinom are gradul m-1, deoarece prezint ă o sum ă de...

CALCUL NUMERIC. Integrarea numerică 1

INTEGRAREA NUMERICĂ

1. APROXIMAREA FUNCłIILOR1 Deseori în cadrul experienŃelor apar serii de rezultate obŃinute pentru anumite valori iniŃiale fixe. Apare

problema prognozării rezultatelor pentru careva condiŃii ini Ńiale, realizarea cărora este dificilă sau imposibilă. Problema se rezolvă prin modelarea experienŃei în baza rezultatelor obŃinute anterior.

Exemplul 1

Fie că a fost măsurată salinitatea apei oceanice într-un punct cu coordonate date, la adâncimile 0, 100, 200, ... , 500 m. Cercetătorii însă trebuie să cunoască salinitatea apei la orice adâncime între 0 şi 500 m. Problema poate fi rezolvată prin modelarea graficului salinităŃii în baza datelor obŃinute experimental. Evident, graficul modelat poate diferi de cel real, dar la adâncimile, în care salinitatea este cunoscută exact, graficele vor coincide.

Exemplul 2

Temperatura unui pacient este măsurată pe parcursul a 24 ore cu intervalul de 4 ore. Se cere să se construiască graficul de temperaturi al pacientului. Graficul de temperaturi se modelează în baza datelor din fişa bolnavului. Graficul modelat al temperaturii poate diferi de cel real, dar în momentele măsurărilor de temperatură ele vor coincide.

În termini matematici problema se formulează în felul următor:

Fie dată o funcŃie f(x), descrisă prin valorile sale în punctele x1, x2, ..., xn (ordonate). Se cere să se estimeze valoarea funcŃiei în oricare din punctele interioare ale segmentului [x1, xn]

Problema aproximării funcŃiilor apare destul de des atât ca problemă aparte, cât şi în calitate de subproblemă în cadrul rezolvării altor probleme.

Exemplul 3

Un careva proces este descris de funcŃia f cu structură complicată, dificilă pentru cercetare. Pentru a simplifica studierea procesului, se aproximează funcŃia iniŃială f cu o altă funcŃie g, de exemplu - un polinom, care poate fi cercetat relativ simplu.

Există diferite metode de aproximare a funcŃiilor. Toate au drept scop o apropiere cât mai mare a graficului funcŃiei care aproximează de graficul funcŃiei aproximate. În cele ce urmează va fi expusă doar una dintre numeroasele metode de aproximare a funcŃiilor - interpolarea funcŃiei prin un polinom.

Fie că se cercetează o funcŃie liniară ( y ax b= + ). Pentru a-i construi graficul, e necesar să se

cunoască două puncte, care îi aparŃin. Într-adevăr, dacă sunt cunoscute 2 puncte ale graficului funcŃiei de gradul 1 (x1, y1) şi (x2, y2) se construieşte sistemul de ecuaŃii:

din care în mod univoc se determină coeficienŃii a şi b. Prin urmare poate fi determinată valoarea funcŃiei în orice alt punct.

Pentru funcŃia polinomială de gradul 2 ( 2y ax bx c= + + ), metoda este aceeaşi: se creează un sistem

din trei ecuaŃii, din care se determină univoc coeficienŃii a, b, c ai funcŃiei iniŃiale. Valoarea funcŃiei în orice alt punct poate fi calculată.

1 OpŃional

=+=+

22

11

ybax

ybax


În caz general, dacă sunt cunoscute m puncte distincte care aparŃin graficului unei funcŃii polinomiale de grad m-1, valoarea funcŃiei poate fi calculată în orice alt punct al domeniului de definiŃie. Se caută un polinom P de forma

1 20 1 2 1( ) ...m m

m mP x a x a x a x a− −− −= + + +

Conform condiŃiei puse

( ) 1,..i iP x y i m= ∀ =

Se construieşte sistemul de m ecuaŃii:

1 20 1 1 1 2 1 1 1

1 20 2 1 2 2 2 1 2

1 20 1 2 1

...

...

...

...

m mm m

m mm m

m mm m m m m m

a x a x a x a y

a x a x a x a y

a x a x a x a y

− −− −

− −− −

− −− −

+ + + + =

+ + + + = + + + + =

Necunoscutele sistemului sunt valorile a0, a1, a2, ... am-1.CoeficienŃii sistemului din determinantul principal formează o serie de puteri, deci determinantul principal, format de ei este diferit de 0. Rangul lui este m (după numărul de ecuaŃii în sistem). În acest caz soluŃia sistemului există şi e unică. Prin urmare, în baza soluŃiilor sistemului dat putem reconstrui polinomul căutat P. Deoarece soluŃia sistemului e unică, rezultă că e unic sistemul de coeficienŃi ai polinomului, prin urmare şi polinomul e unic.

Polinomul astfel determinat este numit polinom de interpolare al punctelor (xi, yi), i = 1, 2,…, m. În cazul, când perechile (xi, yi), i=1, 2,…, m aparŃin graficului unei funcŃii, spunem că polinomul interpolează funcŃia dată.

Dacă funcŃia, graficului căreia îi aparŃin punctele, nu este polinomială, atunci polinomul de interpolare este o aproximare a ei, iar interpolarea poate fi realizată doar pe segmentul determinat de nodurile de interpolare.

Determinarea polinomului de interpolare a unui sistem de puncte prin metoda directă este destul de dificil ă, în special în cazul insuficienŃei de noduri. În cele ce urmează se va discuta o metodă de construcŃie a polinoamelor de interpolare, care poate fi aplicată pentru orice număr de noduri.

Polinoamele simetrice (Polinoame de interpolare Lagrange)

Fie dat un sistem de noduri (xi, yi), i = 1, 2,…,m care aparŃin graficului unei funcŃii. Se va încerca generarea unei serie de polinoame Fi , i=1,...,m, care interpolează punctele

(x1, 0), (x2, 0) … (xi-1, 0), (xi, 1),(xi+1, 0), … (xm, 0) .

Cu alte cuvinte Fi(xi) = 1 , Fj(xi) = 0 ∀ j≠ i.

Pentru ca polinomul să devină 0 în punctele xj diferite de xi vom introduce în el factorii

(x – x1), (x – x2), … (x – xi-1),(x – xi+1), … (x - xm), care îl transformă în 0 pentru x = x1, x2, …xi-1, xi+1, …xm. Deci polinomul va conŃine fragmentul

(x – x1) × (x – x2) ×… × (x – xi-1) × (x – xi+1) × … × (x - xm)

Pentru a transforma acest produs în 1 pentru x = xi, x se înlocuieşte prin xi , iar expresia obŃinută este plasată la numitorul fracŃiei:

(xi – x1) × (xi – x2) ×… (xi – xi-1) × (xi – xi+1) × … × (xi - xm)

în rezultat se obŃine

Cu ajutorul polinoamelor construite mai sus polinomul de interpolare P(x) se defineşte prin formula

1

( ) ( )m

i ii

P x F x y=

=∑

)())(())((

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

1121

1121

miiiiiii

miii xxxxxxxxxx

xxxxxxxxxxxF

−−−−−−−−−−

=+−

+−

KK

KK


Acest polinom are gradul m-1, deoarece prezintă o sumă de termeni, fiecare având gradul m-1, iar pentru valorile xi ale lui x primeşte valorile yi ale funcŃiei interpolate. Prin urmare P(x) este polinom de interpolare pentru funcŃia F(x).

În cazul, când se aproximează o funcŃie, care trece prin punctele (x1, f(x1)), (x2, f(x2)), (x3, f(x3)) …(xm, f(xm)), după excluderea termenilor identic egali cu 0 formula capătă forma

Această formulă poartă denumirea de polinom de interpolare Lagrange.

În cazul când funcŃia f(x) este de m ori derivabilă eroarea de interpolare a funcŃiei f(x) prin polinomul de interpolare Lagrange P(x) este determinată de formula

1

1

( )

[ , ]

( ) ( ) ( )!

unde

sup ( ( ) )m

m

ii

m

x x x

Mf x P x x x

m

M f x

=

∈

− ≤ −

=

∏

m – numărul de noduri de interpolare

Exemplul 4

Să se determine polinomul de interpolare Lagrange pentru funcŃia dată prin perechile (1, 1), (2, 5), (3, 7).

x1 = 1, y1 = 1; x2 = 2, y2 = 5; x3 = 3, y3 = 7;

2 3 1 3 1 21 2 3

1 2 1 3 2 1 2 3 3 1 3 2

2

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

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

( 2)( 3) ( 1)( 2)5( 1)( 3) 7 7 5.

2 2

x x x x x x x x x x x xP x y y y

x x x x x x x x x x x x

x x x xx x x x

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

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

Exemplul 5

Temperatura unui pacient a fost măsurată pe parcursul a n ore cu intervalul de 1 oră. ScrieŃi un program, care să permită determinarea temperaturii în orice moment de timp în intervalul în care au fost făcute măsurările

t 0 1 2 3 4 5 6 f(t) 38 37,1 37.8 39.1 38,4 37,4 36,9

Intrare: prima linie a fişierului in.txt conŃine un număr natural (0< n < 10 – numărul de măsurări. Fiecare din următoarele n linii conŃine câte două numere reale – timpul măsurării temperaturii şi valoarea măsurată (x, y).

Ieşire: temperatura determinată se afişează la ecran.

Rezolvare: pentru comoditate vom elabora un program care va permite introducerea repetată a datelor de la tastatură.

Program cn013; const nmax=10; {num ărul maxim de noduri de interpolare} uses crt; type nod= record {structura nodului} nx,ny :real; end; sistem=array[1..nmax] of nod; var u:sistem; {tabloul nodurilor de interpolare}

i

m

i miiiiiii

mii yxxxxxxxxxx

xxxxxxxxxxxP ×

−−−−−−−−−−

=∑= +−

+−

1 1121

1121

)())(())((

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

KK

KK


xmax,xmin,x:real; {extremit ăŃile segmentului pe care poate fi}

{realizat ă interpolarea şi un punct interior} n:integer; g:text; procedure citirenoduri; var I : integer; begin {deschiderea fi şierului de intrare pentru citire} assign(g,'in1.txt'); reset(g); {citirea datelor} readln(g,n); for i:=1 to n do readln(g, u[i].nx, u[i].ny); close(g); {determinarea extremit ăŃilor segmentului de interpolare} xmax:=u[1].nx; xmin:=u[1].nx; for i:=2 to n do begin if xmax<u[i].nx then xmax:= u.[i].nx; if xmin>u[i].nx then xmin:= u.[i].nx; end; end;

function lag(q:real):real; var i,j:integer; s,p1,p2:real; begin s:=0;{ini Ńializarea valorii totale} for i:=1 to n do begin {calculul urm ătorului termen} p1:=1; p2:=1; {ini Ńializarea num ăr ătorului şi numitorului } for j:=1 to n do {recalculul numitorului şi a num ăr ătorului} if i<>j then begin p1:=p1*(q-u[j].nx); p2:=p2*(u[i].nx-u[j]. nx); end; {ad ăugarea termenului curent la suma total ă} s:=s+(p1/p2)*u[i].ny; end; lag:=s; end; begin clrscr; citirenoduri; repeat writeln ('introdu un punct de pe segmental [',xmin:0:2,';'xmax:0:2,']'); readln(x); writeln (x:0:2,' 'lag(x):0:4); writeln('Alt nod (D)a/(N)u'); until upcase(readkey)='N'; end.

Rezultate: Pentru fluxul de date de intrare


0.5

1.5

2.5

3.5

4.5

5.5

au fost obŃinute rezultatele

37.7390

37.1300

38.6233

38.9878

37.7296

37.4425


Desenul 5. Pe segmentul elementar [xi, xi+1] funcŃia f(x) este aproximată de funcŃia constantă g, iar integrala - prin aria dreptunghiului cu baza l şi înălŃimea g.

2. METODA DREPTUNGHIURILOR PENTRU CALCULUL APROXIMATIV AL INTEGRALEI

Una dintre cele mai des aplicate implementări ale calculului numeric este calcularea integralei definite prin metode aproximative. Metodele directe nu întotdeauna permit calculul analitic al integralei, de multe ori nici nu e cunoscută formula care defineşte funcŃia ce trebuie integrată. De obicei sunt date doar o serie de puncte în care este cunoscută valoarea funcŃiei. În aceste cazuri integrala poate fi calculată doar prin metode aproximative (în presupunerea că funcŃia de sub integrală este continuă pe segmentul pe care se face integrarea).

Din cursul de analiză matematică se ştie, că sensul geometric al integralei definite este aria trapezului curbiliniu, determinat de axa 0X, dreptele x = a şi x = b, şi graficul funcŃiei f(x) pe segmentul [a,b] (desenul 3).

Fie f o funcŃie derivabilă pe [a,b] şi

],[

)('supbax

xfM∈

=

Considerăm o diviziune a segmentului [a,b] în forma a =x0 < x1 < x2 < …< xn-1 < xn = b

Pe fiecare din segmentele [xi, xi+1] vom aproxima supra-faŃa trapezului curbiliniu prin aria dreptunghiului cu baza l = |xi, - xi+1| iar înălŃimea egală cu valoarea funcŃiei în punctul

1

2i i

i

x xz ++=

Pentru simplitate, segmentele [xi, xi+1] se vor considera egale. Atunci lungimea segmentului elementar va fi de

n

abl

|| −=

iar punctele xi vor putea fi determinate de formula xi =a + i * l

FuncŃia constantă g = (xi,+ xi+1)/2 pe fiecare segment elementar [xi, xi+1] este de fapt un polinom de interpolare de grad 0, care aproximează funcŃia f(x) într-un singur punct al segmentului - mijlocul lui.

În acest caz valoarea integralei definite pe segmentul [a, b] pentru funcŃia f se va aproxima conform formulei:

Desenul 4 Divizarea segmentului ini-Ńial în segmente elementare şi aproxi-marea integralei pe fiecare din ele

Desenul 3 Sensul geometric al integ-ralei definite este aria trapezului curbiliniu mărginit de liniile OX, x=a, x=b, y=f(x)


∑∑∫ ∑−

=

+−

=

+−

=+

++=

+=

+−≈

1

0

11

0

11

01 2

)12(22

)()(n

i

iin

i

iib

a

n

iii

liafl

xxfl

xxfxxdxxf

∑−

=

++=1

0 2)12(

n

i

liaflS

n – numărul de divizări a segmentului iniŃial

l – lungimea segmentului elementar

S – valoarea calculată a integralei.

Astfel, integrala se transformă într-o expresie aritmetică care depinde doar de numărul de diviziuni ale segmentului şi de valoarea funcŃiei în punctele de mijloc ale segmentelor elementare. Metoda care reduce calculul integralei la calculul unei sume de arii a dreptunghiurilor este numită metoda dreptunghiurilor.

Determinarea erorii de calcul

Eroarea de interpolare pe un interval elementar se estimează ca modulul diferenŃei între valoarea exactă a funcŃiei pe interval şi funcŃia de interpolare. Conform paragrafului 5.1 această eroare este:

1( ) ( )2

i ix xf x g x M x + − − ≤ −

Corespunzător, eroarea de la integrare va fi integrala erorii de interpolare. Prin urmare, pe un segment elementar eroarea de integrare va fi:

2

)(|)

2(|)()(

211

11 1

ii

x

x

ii

x

x

x

x

xxMdx

xxxMdxxgdxxf

i

i

i

i

i

i

−≤+−≤− ++∫∫ ∫++ +

M – supremul )(xf ′ pe segmentul elementar.

O formulă mai exactă de estimare a erorii se obŃine prin înlocuirea variabilei x:

22)( 11 iiii xx

txx

txx−

++

== ++

( )

2

1111

1

11

22222

)(1 1

−=−+−−++

≤−

+++

−

++∫

∫ ∫+ +

iiiiiiiiii

x

x

x

x

xxMdt

xxxxt

xxxxM

dxxgdxxfi

i

i

i

Prin urmare, eroarea de calcul pa segmentul de integrare [a,b]va fi

unde M – supremul )(xf ′ pe [a,b].

Din formula precedentă rezultă un fapt important: eroarea de calcul este dependentă de numărul de divizări a segmentului. Astfel, pentru a obŃine o eroare de calcul, ce nu depăşeşte o valoare prestabilită ε , este suficient să se realizeze divizarea segmentului de integrare în n segmente elementare, valoarea lui n fiind determinată prin secvenŃa de transformări

4)(

42

22

1 lMab

lM

l

abxxnM ii −=−=

−+


⇒≤−⇒≤− εε

n

MablMab

4

)(

4)(

2

−≥ε4

)( 2 Mabn

Algoritmul

Deoarece în cazul calculelor cu o eroare ce nu depăşeşte valoarea prestabilită numărul necesar de divizări poate fi dedus apriori, este suficient să se realizeze un singur algoritm – pentru un număr fixat de divizări n:

Pas 1 Se introduc datele iniŃiale.

Pentru valoarea dată e a erorii de calcul se determină numărul de divizări din

formula

−≥ε4

)( 2 Mabn

Pas 2 Se calculează pasul de deplasare l. S ⇐ 0.

Pas 3 Pentru toŃi i de la 0 la n-1 :

Se calculează valorile

liax

ilax

i

i

)1(1 ++=+=

+

Se calculează aria dreptunghiului elementar: lxx

fS iii ×

+= +

21

Aria calculată se sumează cu ariile precedente: S ⇐⇐⇐⇐ S+ Si

Pas 4 Se afişează aria totală calculată S. SFÂRŞIT.

Exemplul 1 Fie dată integrala definită dxx

x∫ +

2

1

4

1. Să se scrie un program care

calculează integrala dată prin divizarea segmentului de integrare în 20, 40, 80, 160 divizări. Pentru fiecare număr de divizări să se indice valoarea calculată cu şase semne după virgulă.

Rezolvare: deoarece numărul de divizări este indicat în condiŃiile problemei, preprocesarea matematică nu este necesară.

Program:

program cn014; const r=4; var S, a, b, x, l : real; j,i,n:integer; function f(x:real): real; begin f:=x*x*x*x/sqrt(1+x); end; begin a:=1; b:=2; n:=10; for j:=1 to r do begin S:=0; n:=n*2; l:=(b-a)/n; for i:=0 to n-1 do S:=S+ l*f(a + i*l + l/2); writeln (' n=',n,' I=',S:0:6); end;


readln; end.

Rezultate

n=20 I=3.788513

n=40 I=3.789629

n=80 I=3.789908

n=160 I=3.789977

Exemplul 2 Să se scrie un program, care calculează valoarea aproximativă a integralei

∫−

2

2

2

2

)sin(cosπ

πdx

x prin metoda dreptunghiurilor de mijloc. Calculul va fi oprit, dacă se

îndeplineşte condiŃia: n

abll

−=≤ ;

8επ

, n – numărul curent de divizări

Valorile a, b, ε se definesc direct în program.

Rezolvare: Din condiŃiile problemei se deduce numărul necesar de divizări:

−=⇒≤

−ε

πε

π8

;8

abn

n

ab

Program2:

program cn015; var S, a, b, x, e, l : real; j,i,n:integer; function f(x:real): real; begin f:=sin((cos (x*x))/2); end; begin a:=-pi/2; b:=pi/2; e:=0.0001; n:=round(pi*(abs(b-a))/(8*e))+1; S:=0; l:=(b-a)/n; for i:=0 to n-1 do S:=S+ l*f(a + i*l + l/2); writeln ('n=',n,' I=',S:0:6); readln; end.

Rezultate

n = 12338 I=0.818240

2 Aici şi în alte programe din prezenta ediŃie se foloseşte constanta predefinită pi, utilizată în mediul Turbo Pascal. Utilizatorii altor limbaje de programare sau a altor compilatoare Pascal vor defini constanta respectivă în cadrul programului. Ex: const pi=3.141


VariaŃii ale metodei dreptunghiurilor

În unele cazuri punctele determină înălŃimea dreptunghiului care aproximează integrala pe un segment elementar se iau nu la mijlocul segmentelor elementare, ci la extremităŃile lor. În acest caz, dacă funcŃia se aproximează prin extremităŃile stângi ale segmentelor elementare, integrala definită se aproximează prin suma:

)()2/)()()()(1

0

1

0

1

01 ilaflxflxfxxdxxf

n

ii

n

ii

b

a

n

iii +==−≈ ∑∑∫ ∑

−

=

−

=

−

=+

∑−

=

+=1

0

)(n

i

ilaflS

iar metoda este numită metoda dreptunghiurilor de stânga.

Dacă funcŃia se aproximează prin extremităŃile drepte ale segmentelor elementare, integrala definită se aproximează prin suma:

)()2/)()()()(11

1

1

01 ilaflxflxfxxdxxf

n

ii

n

ii

b

a

n

iii +==−≈ ∑∑∫ ∑

==+

−

=+

∑=

+=n

i

ilaflS1

)(

iar metoda este numită metoda dreptunghiurilor de dreapta.


În cazul variaŃiilor metodei dreptunghiurilor eroarea de interpolare pe un segment elementar este:

))(()()( 1 ii xxxMxgxf −−≤− +

Corespunzător, eroarea de la integrare va fi:

2

)()()(

21

1 1

ii

x

x

x

x

xxMdxxgdxxf

i

i

i

i

−≤− +∫ ∫+ +

M – supremul )(xf ′ pe segmentul elementar

Desenul 6. Aproximarea integralei definite prin dreptunghiuri de stânga

Desenul 7. Aproximarea integralei definite prin dreptunghiuri de dreapta


Prin urmare eroarea de calcul pa segmentul de integrare ],[ ba va fi

unde M – supremul )(xf ′ pe ],[ ba .

La fel ca în cazul dreptunghiurilor de mijloc numărul de divizări necesare pentru a obŃine o soluŃie calculată cu o eroare ce nu depăşeşte o valoare prestabilită ε poate fi dedus direct din formula erorii:

Algoritmul (dreptunghiuri de stânga)3

Deoarece în cazul calculelor cu o eroare ce nu depăşeşte valoarea prestabilită ε numărul necesar de divizări poate fi dedus apriori, este suficient să se realizeze un singur algoritm – pentru un număr fixat de divizări n :

Pas 1 Se introduc datele iniŃiale.

Pentru valoarea dată e a erorii de calcul se determină numărul de divizări din

formula

−≥ε2

)( 2 Mabn

Pas 2 Se calculează pasul de deplasare l. S ⇐ 0.

Pas 3 Pentru toŃi i de la 0 la n-1 :

Se calculează valorile ilaxi +=

Se calculează aria dreptunghiului elementar: lxfS ii ×= )(

Aria calculată se sumează cu ariile precedente: S ⇐⇐⇐⇐ S+ Si

Pas 4 Se afişează aria totală calculată S. SFÂRŞIT.

Exemplul 3. Să se calculeze integrala definită ( )∫π

0

2)sin( dxxx utilizând metoda

dreptunghiurilor, variaŃia dreptunghiurilor de stânga, pentru 10, 100, 1000 divizări. Atribuirea valorilor iniŃiale se face direct în program. Pentru fiecare număr de divizări la ecran se va afişa valoarea calculată şi numărul de divizări, separate prin spaŃiu.

program cn016; const r=3; var S, a, b, x, l : real; j,i,n:integer; function f(x:real): real; begin

3 Pentru dreptunghiurile de dreapta algoritmul se deosebeşte doar în pas 3 prin diapazonul valorilor pe care le primeşte indicele i: de la 1 la n

( )2

)(22

221 l

Mabl

Ml

abxxnM ii −=−=

−+

;;2

)(n

abl

lMab

−=≤− ε ε≤−−n

abMab

2)(

( )

−=ε2

2abMn


f:=sqr((x* sin(x))); end; begin a:=0; b:=pi; n:=1; for j:=1 to r do begin S:=0; n:=n*10; l:=(b-a)/n; for i:=0 to n-1 do S:=S+ l*f(a + i*l); writeln ('n=',n,' I=',S:0:6); end; end.

Rezultate:

n=10 I=4.381796

n=100 I=4.382315

n=1000 I=4.382315

Exemplul 4. Să se calculeze integrala definită ∫1

0

3sin2sinsin xdxxx utilizând metoda

dreptunghiurilor, variaŃia dreptunghiurilor de dreapta, cu o eroare de calcul, ce nu depăşeşte valoarea ε=0.001. Atribuirea valorilor iniŃiale se face direct în program. La ecran se va afişa valoarea calculată a integralei şi numărul de divizări necesare pentru a obŃine precizia cerută, separate prin spaŃiu.

Rezolvare: numărul necesar de divizări poate fi determinat analitic:

Prin calcule elementare se determină valoarea M = 6.

Program:

program cn017; var S, a, b, x, e, l, M : real; j,i,n:integer; function f(x:real): real; begin f:=sin(x)* sin(2*x)* sin(3*x); end; begin a:=0; b:=1; e:=0.001; M:=6; n:=trunc(M*(b-a)*(b-a)/(2*e))+ 1; S:=0; l:=(b-a)/n; for i:=1 to n do S:=S+ l*f(a + i*l); writeln ('n=',n,' I=',S:0:6); end.

Rezultat:

n=3001 I=0.278729

( )

−=ε2

2abMn


Întreb ări şi exerciŃii

1. MotivaŃi necesitatea calculului numeric al integralei definite.

2. Ce procedeu stă la baza calculului numeric al integralei?

3. DemonstraŃi corectitudinea formulelor de integrare.

4. Care este eroarea la integrarea numerică? Cum este ea legată de eroarea de aproximare a funcŃiei?

5. ScrieŃi un program care calculează integrala definită

prin metoda dreptunghiurilor de stânga, de mijloc, de dreapta. UtilizaŃi un număr variabil de divizări (de exemplu 10, 100, 1000). ComparaŃi rezultatele obŃinute prin diferite metode pentru acelaşi număr de divizări.

6. Care dintre integralele definite de mai jos este mai mare?

( ) dxxx∫3

1

2ln dxxex

∫π

0

2cos dxxx∫ −6

1

1

dxxx∫π

0

2)sin( dxxx∫ +1

0

22 31 ∫ ++

π2

0 )sin3)(cos2( xx

dx

dxesaudxe xx

∫∫−−

1

0

1

0

2

dxxesaudxxe xx

∫∫−−

π

π

π 22

0

2 coscos22


Desenul 8. Pe fiecare din segmentele elementare în care este divizat [a,b] valoarea integralei definite este aproximată de aria unui trapez. Latura superioară a trapezului, g(x) este în acelaşi timp aproximarea funcŃiei f(x) pe segmentul elementar[xi, xi+1]

Desenul 9. Aproximarea integralei pe un seg-ment elementar prin aria trapezului determinat de axa 0X, dreptele x=xi, x=xi+1 şi funcŃia liniară g(x)

3. FORMULA TRAPEZELOR

În paragraful precedent pentru aproximarea valorii integralei au fost folosite dreptunghiuri – figuri geometrice, suprafaŃa cărora se calculează prin formule elementare. Un neajuns al acestei metode este numărul mare de divizări, necesare pentru obŃinerea unor rezultate cu erori de calcul suficient de mici. Este logică aproximarea mai exactă a integralei prin alte figuri geometrice, aria cărora poate fi calculată prin formule. Una din figurile, care permit acest lucru este trapezul. Aproximarea integralei definite prin un set de trapeze corespunde aproximării funcŃiei f(x) cu un polinom Lagrange de gradul 1 (o funcŃie liniară care trece prin nodurile de interpolare). EcuaŃia acestei funcŃii este dată de formula:

unde xi şi xi+1 sunt punctele în care valorile funcŃiilor f(x) şi g(x) coincid (desenul 8).

Pentru aplicarea metoda trapezelor, este necesar ca funcŃia integrală să fie de două ori derivabilă. Atunci:

Fie f o funcŃie de două ori derivabilă pe ],[ ba şi

],[

)(''supbax

xfM∈

=

Se consideră o diviziune a intervalului ],[ ba în forma

a = x0 < x1 < x2 < … < xn-1 < xn = b. Pe fiecare din segmentele [xi, xi+1] suprafaŃa trapezului curbiliniu se aproximează cu aria trapezului cu baza l = |xi, - xi+1| , laturile laterale cu lungimile f(xi,) şi f(xi+1) şi latura superioară determinată de segmentul cu extremităŃile (xi, f(xi)) (xi+1, f(xi+1)) (desenul 9).

La fel în cazul formulei dreptunghiurilor se va considera că segmentele [xi, xi+1] sunt egale. Lungi-mea segmentului elemen-tar va fi de:

n

abl

−=

iar punctele xi vor putea fi determinate de formula

xi =a + i * l

Atunci valoarea integralei definite pe un segment elementar [xi, xi+1] pentru funcŃia f(x) se va aproxima conform formulei

2

)()(

2

)()()()()( 11

1

11

+++

+=

+−=≈ ∫∫

++iiii

ii

x

x

x

x

xfxfl

xfxfxxdxxgdxxf

i

i

i

i

iar pe întreg segmentul [a,b]:

)()()(1

11

1i

ii

ii

ii

i xfxx

xxxf

xx

xxxg

+

++

+ −−

+−

−=


∑∑∫∫−

=

+−

=

++

+=

+−=≈

1

0

11

0

11 2

)()(

2

)()()()()(

n

i

iin

i

iiii

b

a

b

a

xfxfl

xfxfxxdxxgdxxf

ii

sau

( )∑−

=++=

1

01)()(

2

n

iii xfxf

lS


La fel ca în cazul formulei dreptunghiurilor eroarea de calcul la integrare va fi cercetată ca integrala erorii de aproximare a funcŃiei f(x) prin g(x) pe fiecare segment elementar:

!2

))(()()( 1 ii xxxxM

xgxf−−

≤− +

unde xi , xi+1 – punctele de interpolare, M – supremul )(xf ′′ pe [xi , xi+1].

Atunci eroarea de integrare4 pe segmentul elementar e de

( ) 31 )(

12)(

1 1

ii

x

x

x

x

xxM

dxxgdxxfi

i

i

i

−≤− +∫ ∫+ +

Eroarea de integrare pe segmentul [a, b] este suma erorilor de integrare pe segmentele elementare, prin urmare poate fi estimată prin formula:

unde M – supremul )(xf ′′ pe [a, b], l – lungimea unui segment elementar

În cazul când în enunŃul problemei se cere calculul integralei prin metoda trapezelor cu o eroare de calcul care nu depăşeşte o valoare dată ε , numărul necesar de divizări a segmentului de integrare poate fi determinat apriori din relaŃia:

ε≤− 2

12)( l

Mab

echivalentă cu ε12

)( 3 Mabn

−≥ ,

sau

−=ε12

)( 3 Mabn

Algoritmul

Pas 1. Se introduc datele iniŃiale: limitele de integrare ,,ba numărul de divizări n . În cazul când în enunŃ se cere obŃinerea unei erori de calcul care nu depăşeşte valoarea predefinită ε, numărul necesar de divizări se calculează cu ajutorul formulei

−=ε12

)( 3 Mabn

4 Demonstrarea corectitudinii acestei formule poate fi propusă în calitate de exerciŃiu realizat în afara orei de curs

231 12

)()(12

)()( lM

abxxM

ndxxgdxxf ii

b

a

b

a

−=−≤−− +∫ ∫


Pas 2. Se calculează pasul de deplasare l după formula n

abl

−= . Se iniŃializează valoarea

integralei S ⇐ 0, a contorului i ⇐ 0.

Pas 3. Pentru toŃi i de la 0 la n -1

Se calculează valorile ilaxi +=

liaxi )1(1 ++=+

Se calculează aria trapezului elementar:

2

)()( 1++= ii

i

xfxflS

Aria calculată se sumează cu ariile precedente:

S ⇐⇐⇐⇐ S + Si

Pas 4 Se afişează aria totală calculată S. SFÂRŞIT

Exemplul 1.

Să se calculeze integrala definită ∫ +

π2

044 cossin

1dx

xx, utilizând metoda trapezelor,

pentru 20, 40, 80 divizări a segmentului de integrare. Atribuirea valorilor iniŃiale se face direct în program. Pentru fiecare număr de divizări la ecran se va afişa valoarea calculată şi numărul de divizări, separate prin spaŃiu.

Rezolvare: numărul de divizări ale segmentului de integrare se modifică de la caz la caz prin formula n ⇐ 2n. Celelalte calcule rezultă direct din algoritm.

Program

program cn018; const r=3; var S, a, b, x, l : real; j,i,n:integer; function f(x:real): real; begin f:=1/(sqr(sqr(sin(x)))+sqr(sqr(cos(x)))); end; begin a:=0; b:=2*pi; n:=10; for j:= 1 to r do begin S:=0; n:=n*2; l:=(b-a)/n; for i:=0 to n-1 do S:=S+ l*(f(a + i*l)+f(a+(i+1)*l) )/2; writeln ('n=',n,' I=',S:0:11); end; end.


Rezultate:

n=20 I=8.88312405500

n=40 I=8.88576626920

n=80 I=8.88576587630

Exemplul 2.

Să se calculeze integrala definită ∫ +−3

1

23 )273( dxxx , utilizând metoda trapezelor, cu

o eroare de calcul, ce nu depăşeşte valoarea ε=0.001. Atribuirea valorilor iniŃiale se realizează direct în program. La ecran se va afişa valoarea calculată a integralei şi numărul de divizări necesar pentru a obŃine exactitatea cerută.

Rezolvare: numărul de divizări ale segmentului de integrare se calculează direct cu ajutorul

formulei

−=ε12

)( 3 Mabn . Prin calcule elementare se stabileşte M=42.

Program program cn019; var S, a, b, x, e, l, M : real; j,i:integer; n:longint; function f(x:real): real; begin f:=3*x*x*x-7*x*x+2; end; begin a:=1; b:=3; e:=0.001; M:=42; n:=trunc(sqrt(M*(b-a)*(b-a)*(b-a)/(12*e)))+ 1; S:=0; l:=(b-a)/n; for i:=0 to n-1 do S:=S+ l*(f(a + i*l)+ f(a + (i+1)*l))/2; writeln ('n=',n,' I=',S:0:9); end.

Rezultate:

n=168 I=3.598890846

Întreb ări şi exerciŃii

1. EnumeraŃi condiŃiile în care poate fi aplicată metoda trapezelor.

2. Care este deosebirea între metoda dreptunghiurilor şi cea a trapezelor?

3. Care dintre metodele studiate de calculul numeric al integralei definite este mai exactă? MotivaŃi.

4. CalculaŃi integralele definite, utilizând metoda trapezelor pentru 10, 20, 30 divizări ale segmentului de integrare:


∫+

4

1

4dxe x dxxx∫ +−2

1

2 453 ∫

++

4

22

3

1

2dx

x

xx

∫π

0

2 cos)sin( xdxxx ∫ +1

0

cossin dxxx ∫ +π

π3

2

23 )sincos( dxxxtgx

INTEGRAREA NUMERIC Ă - Math · Acest polinom are gradul m-1, deoarece prezint ă o sum ă de...

Documents

Transcript of INTEGRAREA NUMERIC Ă - Math · Acest polinom are gradul m-1, deoarece prezint ă o sum ă de...