Lucrarea de laborator nr. 7 I. Scopul lucrăriiMetode Numerice - Lucrarea de laborator 7 1 Lucrarea...
Transcript of Lucrarea de laborator nr. 7 I. Scopul lucrăriiMetode Numerice - Lucrarea de laborator 7 1 Lucrarea...
Metode Numerice - Lucrarea de laborator 7
1
Lucrarea de laborator nr. 7
I. Scopul lucrării
Derivarea aproximativă
Integrarea numerică
II. Conținutul lucrării
1. Formule de derivare aproximativă folosind dezvoltări în serie Taylor
2. Proceduri MAPLE pentru implementarea formulelor de derivare aproximativă.
Exemple
3. Integrarea numerică
- Formula dreptunghiurilor.
- Formula trapezelor
III. Prezentarea lucrării
III.1. Formule de derivare aproximativă folosind dezvoltări în serie
Taylor
Fie f : [a, b] R, derivabilă. Se recurge la aproximarea derivatei unei funcții f
sau a derivatelor ei de ordin superior atunci când expresia lui f este prea complicată
sau când nu se cunoaște expresia analitică a funcției f (f este dată prin intermediul unui
tabel de valori).
Fie f : [a, b] R, derivabilă. Presupunem că se dau n+1 puncte distincte în
intervalul [a, b], x0, x1, …., xn, în care se cunosc valorile funcției f.
Prezentăm o tehnică de găsire a unor formule de aproximare pentru valorile
derivatelor de orice ordin ale funcției f. Această tehnică are la bază formula lui Taylor.
Mădălina Roxana Buneci Metode Numerice –Laborator
2
Reamintim această formulă.
Fie I un interval de numere reale, a I și f : I R o funcție de n ori derivabilă
în a. Polinomul Taylor de ordin n atașat lui f în punctul a este funcția polinomială Ta,n
: I R, definită prin:
Ta,n(x) = f(a) + !1
a'f(x-a) +
!2
a"f(x-a)2 +
!n
af n
(x-a)n.
Restul formulei Taylor de ordin n atașat funcției f în punctul a este funcția Ra,n
: I R definită prin Ra,n(x) = f(x) - Ta,n(x). Egalitatea
f(x) = Ta,n(x) + Ra,n(x)
valabilă pentru orice x I se numește formulă Taylor de ordin n atașată funcției f în
punctul a. Se demonstrează că n
x,a
ax ax
Rlim
= 0. Dacă f este de clasă Cn+1, atunci există
u strict cuprins între a și x (sau echivalent există (0, 1) astfel încât u = a+ (x – a))
cu proprietatea că
Ra,n(x) = !1n
uf 1n
(x - a)n+1.
Ne propunem să aproximăm valorile derivatelor de ordinul întâi și doi ale
funcției f în punctul xi [a, b] când sunt cunoscute valorile funcției în punctele
xi-1, xi, xi+1
care nu sunt neapărat echidistante, adică:
xi-1 = xi –h, xi+1 = xi + h.
Presupunem că f este de 3 ori derivabilă și scriem formula lui Taylor de ordinul 2
atașată lui f în xi:
2 3i i xi i i i
f (x ) f (x ) f (w )f (x) f (x ) (x x ) (x x ) (x x )
1! 2! 3!
.
cu wx stric cuprins între x și xi. Înlocuim în această formulă x = xi-1, respectiv x = xi+1
și obținem relațiile:
23i i
i 1 if (x ) f (x )h f (u)
f (x ) f (x ) ( h) ( h)1! 2! 3!
(1)
Metode Numerice - Lucrarea de laborator 7
3
2 2 3 3i ii 1 i
f (x ) f (x ) f (v)f (x ) f (x ) h h h
1! 2! 3!
(2)
cu u strict cuprins între xi-1 și xi și v strict cuprins între xi și xi+1. Scăzând din a doua
relație prima relație înmulțită cu 2 obținem:
2 2 3 3 3 2i 1 i i 1
i
f (v) f (u)f (x ) f x 1 f (x ) h h
3! 3!f (x )h 1
Putem aproxima valoarea derivatei funcției f în xi prin:
2 2i 1 i i 1
i
f (x ) f x 1 f (x )f (x )
h 1
cu eroarea:
22f (v) f (u)
h1 3! 1 3!
care tinde la zero odată cu h2.
Adunând la a doua relație prima înmulțită cu obținem:
3 3 3i 1 i i 1
i 2
f (v) f (u)2 f (x ) 1 f (x ) f (x ) h h
3! 3!f (x )
h 1
.
Putem astfel aproxima valoarea derivatei de ordinul doi al lui f în xi prin
i 1 i i 1i 2
2 f (x ) 1 f (x ) f (x )f (x )
h 1
cu eroarea
2h f (v) f (u)
1 3 3
care tinde la zero odată cu h.
Dacă nodurile sunt echidistante (adică dacă = 1) se obțin următoarele formule
de aproximare:
i 1 i 1i
f (x ) f (x )f (x )
2h
Mădălina Roxana Buneci Metode Numerice –Laborator
4
i 1 i i 1i 2
f (x ) 2f (x ) f (x )f (x )
h
.
III.2. Proceduri MAPLE pentru implementarea formulelor de derivare
aproximativă. Exemple
Procedură MAPLE pentru determinarea valorilor aproximative ale
derivatei de ordinul 1
Procedura d1 admite drept parametri lista x ce conține punctele x1, x2, …., xn
și lista y ce conține valorile y1 = f(x1), y2 = f(x2), …., yn= f(xn). Procedura întoarce lista
aproximațiilor derivatei de ordinul I a lui f în punctele x2, x3, ..., xn-1.
> d1:=proc(x,y)
local df,h,alpha,i,n;
n:=nops(x);df:=[seq(1,i=1..n-2)];
for i from 2 to n-1 do
h:=x[i]-x[i-1];
alpha:=(x[i+1]-x[i])/h;
df[i-1]:=(y[i+1]+y[i]*(alpha^2-1)-alpha^2*y[i-1])/
(h*alpha*(alpha+1))
end do;
return df;
end proc;
Procedură MAPLE pentru determinarea valorilor aproximative ale
derivatei de ordinul 2
Procedura d2 admite drept parametri lista x ce conține punctele x1, x2, …., xn
și lista y ce conține valorile y1 = f(x1), y2 = f(x2), …., yn= f(xn). Procedura întoarce lista
Metode Numerice - Lucrarea de laborator 7
5
aproximațiilor derivatei de ordinul al II-lea a lui f în punctele x2, x3, ..., xn-1.
> d2:=proc(x,y)
local d2f,h,alpha,i,n;
n:=nops(x);d2f:=[seq(1,i=1..n-2)];
for i from 2 to n-1 do h:=x[i]-x[i-1];
alpha:=(x[i+1]-x[i])/h;
d2f[i-1]:=2*(y[i+1]-(alpha+1)*y[i]+alpha*y[i-
1])/(h^2*alpha*(alpha+1))
end do;
return d2f;
end proc;
Exemple
> x1:=[seq(-1+2*i/5,i=0..5)];
> f1:=t->ln(1+t^2);
> y1:=[seq(evalf(f1(-1+2*i/5)),i=0..5)];
> d1(x1,y1);
> d2(x1,y1);
> z1:=[seq(D(f1)(-1+2*i/5),i=1..4)];
> z2:=[seq((D@@2)(f1)(-1+2*i/5),i=1..4)];
:= x1
, , , , ,-1
-3
5
-1
5
1
5
3
51
:= f1 t ( )ln 1 t2
y1 0.6931471806 0.3074846997 0.03922071315 0.03922071315 0.3074846997, , , , ,[ :=
0.6931471806]
[ ], , ,-0.8174080840 -0.3353299832 0.3353299832 0.8174080840
[ ], , ,0.7337405900 1.676649916 1.676649916 0.7337405900
:= z1
, , ,
-15
17
-5
13
5
13
15
17
Mădălina Roxana Buneci Metode Numerice –Laborator
6
> map(evalf,d1(x1,y1));
> map(evalf,z1);-.8174080840
> map(evalf,d2(x1,y1));
> map(evalf,z2);
III.3. Integrarea numerică
Fie f : [a, b] R o funcție continuă. Ne punem problema să calculăm valoarea
aproximativă a integralei
b
a
f (x) (x)dx , unde : [a, b]R este o funcție continuă strict
pozitivă numită pondere. Considerăm x0, x1, …, xn n+1 puncte distincte din intervalul
[a, b], și notăm yi = f(xi) pentru orice i = 0,1,…n. Fie Ln polinomul Lagrange asociat
lui f și punctelor considerate:
Ln(x) =
n0 1 i 1 i 1 n
ii 0 i 1 i i 1 i i 1 i ni 0
x x x x .... x x x x ... x xy
x x x x .... x x x x ... x x
Înlocuind f prin Ln, obținem formula de aproximare
b b
n
a a
f (x) (x)dx L (x) (x)dx
Reprezentăm punctele sub forma
xi = ia b b a
t2 2
, ti [-1, 1], i 0, 1, ..., n
și obținem
:= z2
, , ,
200
289
300
169
300
169
200
289
[ ], , ,-0.8174080840 -0.3353299832 0.3353299832 0.8174080840
[ ], , ,-0.8823529412 -0.3846153846 0.3846153846 0.8823529412
[ ], , ,0.7337405900 1.676649916 1.676649916 0.7337405900
[ ], , ,0.6920415225 1.775147929 1.775147929 0.6920415225
Metode Numerice - Lucrarea de laborator 7
7
b
a
f (x) (x)dx
1n0 i 1 i 1 n
ii 0 i i 1 i i 1 i ni 0 1
(t t ) (t t )(t t ) (t t )b a a b b af (x ) ( t) dt.
2 2 2 (t t ) (t t )(t t ) (t t )
( formula generală de cuadratură)
Restul formulei generale de cuadratură (eroarea absolută) este:
Rn(f)
M
n 1 !
1n 2
0 1 n
1
b a a b b at (t t )(t t ) (t t )dt
2 2 2
,
unde M = x [a,b]
sup
f(n+1)(x).
III.3.1. Formula dreptunghiurilor
Fie f : [a, b] R o funcție de clasă C1. Aplicând formula generală de
cuadratură pentru 1, n=0, x0 = a b
2
(deci t0 = 0) obținem
b
a
f (x)dx (b-a)a b
f2
cu o eroare
2b a
4
x a,b
sup f ' x
.
Considerăm o diviziune (x0, x1, …., xn) a intervalului [a, b] cu puncte
echidistante (xi = a + n
abi, i = 0, 1, 2, …, n.) și aplicăm pe fiecare subinterval [xi,xi+1]
formula de aproximare
x
i 1x
i
f x dx (xi+1-xi) i i 1x xf
2
Se obține astfel formula dreptunghiurilor:
Mădălina Roxana Buneci Metode Numerice –Laborator
8
b
a
f (x)dx b a
n
n 1i i 1
i 0
x xf
2
.
Restul (eroarea) este dat de:
n 1
b i i 1
ai 0
x xb af x dx f
n 2
2
b a
4n
x a,b
sup f ' x
Interpretarea geometrică a formulei dreptunghiurilor
Fie Di dreptunghiul cu o dimensiune dată intervalul [xi, xi+]] și cu cealaltă
dimensiune dată de i i 1x xf
2
Atunci aria dreptunghiului Di este
(xi+1-xi) i i 1x xf
2
=b a
n
i i 1x xf
2
,
și deci formula dreptunghiurilor presupune aproximarea
b
a
f (x)dx prin suma ariilor
xi xi+1
Metode Numerice - Lucrarea de laborator 7
9
dreptunghiurilor Di, i = 0, 1, …n-1.
Proceduri MAPLE pentru calculul valorii aproximative a unei integrale
definite folosind formula dreptunghiurilor
Procedura dreptunghiuri1 are drept parametri funcția care se integrează,
limitele de integrare, și numărul de subintervale din diviziune. Procedura întoarce
valoarea aproximativă a integralei obținută aplicând formula dreptunghiurilor.
Procedura dreptunghiuri2 este similară, cu singura deosebire că în locul numărului de
subintervale se introduce un număr pozitiv eps ce reprezintă eroarea maximă.
> dreptunghiuri1:=proc(f,a,b,n)
local i,iab,h,h0;
iab:=0;h:=(b-a)/n;h0:=a+1/2*h;
for i from 0 to n-1 do
iab:=iab+evalf(f(h0+i*h))
end do;
iab:=iab*h;
return evalf(iab)
end proc;
> dreptunghiuri2:=proc(f,a,b,eps)
local i,iab,h,h0,n;n:=floor(1/4*(b-a)^2/eps)+1;
print(`Numar de pasi`,n); h:=(b-a)/n;iab:=0;h0:=a+1/2*h;
for i from 0 to n-1 do
iab:=iab+evalf(f(h0+i*h))
end do;
iab:=iab*h;
return evalf(iab)
Mădălina Roxana Buneci Metode Numerice –Laborator
10
end proc;
Exemple
> f:=(x->x^7*ln(x)+x*cos(x));
> evalf(int(f(x),x=2..3));
> dreptunghiuri1(f,2,3,5);
> dreptunghiuri1(f,2,3,50);
> dreptunghiuri1(f,2,3,500);
> dreptunghiuri2(f,2,3,0.01);
> dreptunghiuri2(f,2,3,0.001);
> dreptunghiuri2(f,2,3,0.0001);
> g:=(x->exp(-x^2));
> evalf(int(g(x),x=-1..1));
> dreptunghiuri1(g,-1,1,5);
:= f x x7 ( )ln x x ( )cos x
778.3339881
768.4434052
778.2346340
778.3329936
,Numar de pasi 26
777.9666004
,Numar de pasi 251
778.3300446
,Numar de pasi 2501
778.3339476
:= g x e( )x
2
1.493648266
Metode Numerice - Lucrarea de laborator 7
11
> dreptunghiuri1(g,-1,1,10);
> dreptunghiuri1(g,-1,1,500);
> dreptunghiuri2(g,-1,1,0.01);
> dreptunghiuri2(g,-1,1,0.001);
> dreptunghiuri2(g,-1,1,0.0001);
> with(student):
> middlesum(g(x),x=-1..1,5);
> evalf(middlesum(g(x),x=-1..1,5));
> middlebox(g(x),x=-1..1,5);
1.503548970
1.496106505
1.493649246
,Numar de pasi 101
1.493672307
,Numar de pasi 1001
1.493648512
,Numar de pasi 10001
1.493648267
2
5
i 0
4
e
/4 5
2 i
5
2
1.503548970
Mădălina Roxana Buneci Metode Numerice –Laborator
12
> middlesum(g(x),x=-1..1,10);
> evalf(middlesum(g(x),x=-1..1,10));
> middlebox(g(x),x=-1..1,10);
Comanda middlesum(g(x), x=a..b,n) din pachetul student întoarce aproximația
integralei
b
a
g(x)dx obținută prin aplicarea formulei dreptunghiurilor utilizând n
subintervale. Comanda middlebox(g(x), x=a..b,n) reprezintă grafic dreptunghiurile
utilizate în formulei dreptunghiurilor.
1
5
i 0
9
e
/9 10
i
5
2
1.496106505
Metode Numerice - Lucrarea de laborator 7
13
III.3.2. Formula trapezelor
Fie f : [a, b] R o funcție de clasă C2. Aplicăm formula generală de cuadratură
pentru 1, n=1, x0 = a, x1 =b (deci t0 = -1 și t1 =1). Obținem
b
a
f (x)dx b a
2
(f(a) +f(b)).
cu eroarea
b
a
f (x)dx - b a
2
(f(a) +f(b))
3
b a
12
x a,b
sup f " x
.
Considerăm o diviziune (x0, x1, …., xn) a intervalului [a, b] cu puncte
echidistante (xi = a + b a
n
i, i = 0, 1, 2, …, n.) și aplicăm pe fiecare subinterval [xi,xi+1]
formula de aproximare
x
i 1x
i
f x dx i 1 ix x
2
( f(xi) + f(xi+1))
Obținem următoarea formulă de cuadratură numită formula trapezelor:
b
a
f (x)dx b a
n
n 1
ii 1
f a f bf x
2
.
Restul (eroarea) este dat de:
n 1
b
iai 1
f a f bb af x dx f x
n 2
3
2
b a
12n
x a,b
sup f " x
Interpretarea geometrică a formulei trapezelor
Fie Ti trapezul dreptunghic cu înălțimea egală cu lungimea intervalului [xi, xi+]]
și cu bazele f(xi) și f(xi+1).
Mădălina Roxana Buneci Metode Numerice –Laborator
14
Atunci aria trapezului Ti este
i 1 ix x
2
( f(xi) + f(xi+1)) =
b a
2n
( f(xi) + f(xi+1)),
și deci formula trapezelor arată că b
af x dx se poate aproxima prin suma ariilor
trapezelor Ti, i = 0, 1, …n-1.
Proceduri MAPLE pentru calculul valorii aproximative a unei integrale
definite folosind formula trapezelor
Procedura trapeze1 are drept parametri funcția care se integrează, limitele
de integrare, și numărul de subintervale din diviziune. Procedura întoarce valoarea
aproximativă a integralei obținută aplicând formula trapezelor. Procedura trapeze2
este similară, cu singura deosebire că în locul numărului de subintervale se introduce
un număr pozitiv eps ce reprezintă eroarea maximă.
> trapeze1:=proc(f,a,b,n)
local i,iab,h;
iab:=0;h:=(b-a)/n;;
for i from 1 to n-1 do
iab:=iab+evalf(f(a+i*h))
end do;
f(xi+1)
f(xi)
xi xi+1
Metode Numerice - Lucrarea de laborator 7
15
iab:=iab+(f(a)+f(b))/2;iab:=iab*h;
return evalf(iab)
end proc;
> trapeze2:=proc(f,a,b,eps)
local i,iab,h,n;n:=floor(abs((b-
a)^3/(12*eps))^(1/2))+1;
print(`Numar de pasi`,n); h:=(b-a)/n;iab:=0;
for i from 1 to n-1 do
iab:=iab+evalf(f(a+i*h))
end do;
iab:=iab+(f(a)+f(b))/2;iab:=iab*h;
return evalf(iab)
end proc;
Exemple
> f:=(x->x^7*ln(x)+x*cos(x));
> evalf(int(f(x),x=2..3));
> trapeze1(f,2,3,5);
> trapeze1(f,2,3,50);
> trapeze1(f,2,3,500);
> trapeze2(f,2,3,0.01);
> trapeze2(f,2,3,0.001);
:= f x x7 ( )ln x x ( )cos x
778.3339881
798.1539466
778.5326999
778.3359754
,Numar de pasi 3
833.1348363
Mădălina Roxana Buneci Metode Numerice –Laborator
16
> trapeze2(f,2,3,0.0001);
> g:=(x->exp(-x^2));
> evalf(int(g(x),x=0..1));
> trapeze1(g,0,1,5);
> trapeze1(g,0,1,50);
> trapeze1(g,0,1,500);
> trapeze2(g,0,1,0.01);
> trapeze2(g,0,1,0.001);
> trapeze2(g,0,1,0.0001);
> trapeze2(g,0,1,10^(-8));
,Numar de pasi 10
783.2986759
,Numar de pasi 29
778.9246586
:= g x e( )x
2
0.7468241330
0.7443683397
0.7467996064
0.7468238866
,Numar de pasi 3
0.7399864752
,Numar de pasi 10
0.7462107961
,Numar de pasi 29
0.7467512252
,Numar de pasi 2887
0.7468241295
Metode Numerice - Lucrarea de laborator 7
17
> with(student):
> trapezoid(g(x),x=-1..1,5);
> evalf(trapezoid(g(x),x=-1..1,5));
> trapezoid(g(x),x=-1..1,10);
> evalf(trapezoid(g(x),x=-1..1,10));
Comanda trapezoid(g(x), x=a..b,n) din pachetul student întoarce aproximația integraleib
a
g(x)dx obținută prin aplicarea formulei trapezelor utilizând n subintervale.
2
5e
( )-1 2
5
i 1
4
e
1
2 i
5
2
1.473924388
1
5e
( )-1 1
5
i 1
9
e
1
i
5
2
1.488736679