INTEGRARE ŞI DERIVARE NUMERICĂ - andrei.clubcisco.roandrei.clubcisco.ro/cursuri/1mn/curs/c09mn...
Transcript of INTEGRARE ŞI DERIVARE NUMERICĂ - andrei.clubcisco.roandrei.clubcisco.ro/cursuri/1mn/curs/c09mn...
1
INTEGRARE ŞI DERIVARE NUMERICĂ
Ne propunem în acest capitol să calculăm în mod aproximativ valorile
dxxffIb
a ,
0pxffD .
în condiţiile în care
funcţia f este continuă pe b,aCf:b,a şi derivabilă în 0x
primitiva F nu este cunoscută
funcţia f este cunoscută numai prin valorile f(xi) pe care le ia într-un număr restrîns de
puncte xi, i=0 : N
Definim o metodă aproximativă de integrare ca
N
1i
iNiNN xfAfI ,
Metoda aproximativă de integrare este slab convergentă dacă
0fIfIlim NN
.
In mod similar se defineşte o metodă aproximativă de derivare.
Teorema 7.1. Condiţia necesară şi suficientă ca metoda de integrare IN[f] să conveargă slab către
I[f] se exprimă prin relaţiile
a) există M>0 astfel încât Ma
N
1i
iN
, pentru toţi N=1,2,...
b)
b
a
kk
Nn
,dxxxIlim pentru toţi k=0,1,...
1. Metode de tip Newton-Cotes
In general, pentru o formulă de integrare aproximativă putem scrie
N
N
1i
iNiN
b
aRxfAdxxwxf
.
funcţia pondere w:[a,b]R+, nu modifică problema (1), întrucât putem lua g(x)=f(x).w(x), iar
Rn este eroarea (sau restul) formulei aproximative de integrare.
Metodele de tip Newton-Cotes se bazează pe integrarea polinomului de interpolare, utilizând ca suport al
interpolării nodurile xiN echidistante în intervalul [a,b], adică
N:0i,N
abiaxiN
.
Metodele de integrare de tip Fejer integrează polinomul de interpolare folosind ca noduri xiN - rădăcinile
polinomului ortogonal Pn(x), definit relativ la ponderea w(x).
Coeficienţii aiN se determină impunând ca formula aproximativă să fie exactă (R=0), dacă f aparţine
unei anumite clase de funcţii (de exemplu polinoame de grad ≤N, f N,).
Cum funcţia este cunoscută numai în nodurile xi, i=1:N, o vom aproxima prin polinomul ei de
interpolare Lagrange
iN
N
1i
i1N xfxlxPxf
,
cu care putem scrie
2
N
1i
iNiN
b
a1N xfAdxxwxP ,
sau
dxxwxfxldxxwxP iN
b
a
N
1i
iN
b
a1N
iN
N
1i
iN
b
aiN
N
0i
iN xfAdxxwxlxf
,
de unde
b
aiNiN dxxwxlA .
Printr-o schimbare liniară de variabilă, coeficienţii aiN pot fi făcuţi independenţi de intervalul de
integrare; ei sunt totuşi inutilizabili, fiind de valori mari şi de semne contrarii, ceea ce conduce la
instabilitate numerică.
Expresia erorii în metodele de tip Newton-Cotes se deduce integrând expresia erorii din polinomul de
interpolare.
xExPxf 1N1N ,
obţinându-se
NNR
b
a1N
fI
b
a1N
fI
b
adxxwxEdxxwxPdxxwxf ,
deci
b,aξ,dxxwxxxx
!N
ξfR
b
aN1
N
1N ,
cu majorarea
dxxwxxxx
!N
ξfR
b
aN1
N
1N .
Datorită instabilităţii interpolării polinomiale se folosesc polinoame de interpolare cu grad mic.
Astfel pentru N=1 se obţine formula trapezelor
12
fhbfaf
2
hdxxf
3b
a
.
în care abN
abh
şi b,a
Pentru N=2 se obţine formula lui Simpson
90
ξfhbf
2
baf4af
3
hdxxf
iv5b
a
.
Aceste formule folosesc puţine puncte ceeace ne determină să aproximăm funcţia f, local, pe intervale
dxxfdxxfdxxfn
1n
1
0
n
0
x
x
x
x
x
x
,
Se obţin în acest fel
formula compusă a trapezelor
3
1N
1i
ihaf2bfaf2
hT .
cu N
abh
.
function I = Trapez(a, b, n, f)
% Intrări:
% a, b = capetele intervalului de integrare
% n = ordinul metodei
% f = funcţia de integrat
% Ieşiri:valoare integrala definita
h = (b-a) / n;
s = 0;
for i = 1 : n-1
s = s + f(a+i*h);
end
I = h*(f(a) + f(b) + 2*s) / 2;
formula compusă Simpson
N
1i
1N
1i
i21i2 xf2xf4bfaf3
hS .
cu N2:0i,hiax,N2
abh i
function I= Simpson(n, a, b, f)
/* Intrări:
/* a,b = capete interval de integrare
/* n = ordinul metodei
/* f = funcţia de integrat
/* Ieşiri: valoarea integralei definite
h =(b-a) / (2*n);
s1 = 0;
s2 = 0;
for i = 1 : n
s1 = s1 + f(a+(2*i-1)*h);
end
for i = 1 : n-1
s2 = s2 + f(a+2*i*h);
end
I = h*(f(a) + f(b) + 4*s1 + 2*s2)/3;
7.2. Formule de integrare bazate pe integrarea prin părţi
Metoda lui Euler
Metodele din această categorie sunt aplicabile, dacă se cunosc informaţii privind derivatele funcţiei de
integrat f(x) la capetele intervalului de integrare.
Considerăm funcţia f(x) continuă, împreună cu derivatele până la ordinul r inclusiv
( b,aCf,Rb,a:fr ) şi polinomul monic de grad ≤r, rr xP .
Întrucât Pr(x) este monic, avem evident !rxPr
r şi putem scrie
b
a
r
r
b
a
dxxfxP!r
1dxxffI
Integrăm prin părţi de r ori
4
dxxfxP
!r
1xfxP1
!r
1fI
b
a
r
r
r
a
br
1k
1kkr
r
1k
Considerăm dezvoltarea în serie
0n
nn
t
tx
t!n
xB
1e
et
relaţie scrisă şi sub forma
tx
0n
nntett
!n
xB1e
, în care exponenţialele se dezvoltă de asemenea în
serie Taylor
0n
nn
0n 1n
n1nn
t!n
xB
!n
t
!n
tx
Identificând termenii din cei doi membri obţinem
n
1k
1n
kn0 xnxBk
n,1xB ,
relaţie care ne arată că Bn(x) este un polinom de grad n numit polinomul lui Bernoulli.
Folosind relaţia cu n=1,2,… obţinem
2
xx
2
3xxB,
6
1xxxB,
2
1xxB
23
3
2
21 ,…
Se obţin relativ uşor relaţiile
1n
nn xnxB1xB
xBnxB 1nn , n=1,2,…
xB1x1B n
n
n
1
0
n 0dxxB
Numerele lui Bernoulli se definesc ca
0BBnn
Dacă se notează
!n
BC n
n
particularizăm relaţia de recurenţă pentru x=0
01Cn
1nC
2
1nC
1
1n
11nn
,
din care se determină, luând pe rând n=1,2,… valorile
,6
1C,
2
1C 21 şi
0B 1n2
Particularizăm relaţia pentru
!n2
xBxP,n2r
n2
r şi facem schimbarea de variabilă
x=a+t.h cu h=b-a
obţinem
5
b,aξ,ξfBhfR
afbfBh2
bfafhfI
n2
n2
1n2
1n
1k
1k21k2
k2
k2
n
Pentru obţinerea unei formule compuse, împărţim intervalul (a,b) cu N-1 puncte echidistante
N
abiaxi
şi aplicăm formula de mai sus în fiecare subinterval. Obţinem în final
o relaţie cunoscută sub numele de formula lui Euler
1n
1k
1k21k2k2
k21N
1i
n afbf!k2
Bhihaf
2
bfafhfI
3. Integrare gaussiană
Vom indexa nodurile de la 1, xin, i=1:n este justificată, deoarece prezenţa a n noduri impune un grad
al polinomului de interpolare egal cu n-1. În cazul integrării gaussiene, acelaşi suport reprezintă cele n
rădăcini ale unui polinom ortogonal de grad n. Aşadar, pentru calculul unei integrale de forma
dxxwxffI
b
a
cu Rb,a:w , o funcţie continuă pe intervalul finit sau infinit (a,b), vom folosi o metodă de
integrare numerică de forma
N
1i
iNiN1N xfafI
În plus se impun condiţiile
,1,0k,dxxwx
b
a
k să fie absolut convergente
b
a
2dxxwxf
Metodele de tip Newton-Cotes (pentru suportul NNN1 x,,x ) au gradul de valabilitate N-1 (sunt
exacte pentru polinoame până la gradul N-1 inclusiv).
Ne punem problema determinării unor metode aproximative de integrare cu grad de valabilitate mai
ridicat, pe seama alegerii corespunzătoare a nodurilor xiN.
Determinarea celor 2N necunoscute aiN şi xiN cu i=1:N, necesită 2N ecuaţii. Formula este, prin urmare
exactă pentru polinoame de la grad 0 până la 1N2 .
1N2N P,PIPI
Sistemul format, considerând polinoamele 1N2x,,x,1P
este
dxxwAb
a
N
1i
iN
,
b
a
N
1i
iNiN dxxwxxA ,
dxxwxxAb
a
1N2N
1i
1N2
iNiN
.
fiind neliniar în raport cu xiN, şi nu poate fi rezolvat uşor în mod direct.
Pentru obţinerea nodurilor xiN, i=1:N, Gauss a folosit metoda prezentată în cele ce urmează.
Seconsideră polinomul
6
NNN1N xxxxxπ ,
Orice polinom xP 1qN de grad N21qN se poate exprima sub forma
xRxQxπxP 1N1qN1qN ,
pe baza teoremei împărţirii cu rest.
Dorim ca formula aproximativă de integrare să aibă grad de valabilitate > N-1, adică să fie exactă
pentru polinoamele xP 1qN , ceeace implică
iN
N
1i
1qNiN
b
a1qN xPAdxxwxP
,
dar:
iN1NiN1qN xRxP ,
sau
iN
N
1i
1NiN
b
a1qN xRAdxxwxP
.
Pe de altă parte
b
a1N1q
b
aN
b
a1qN dxxwxRdxxwxQxπdxxwxP ,
dar formula are cel puţin gradul de valabilitate N-1, adică
N
1i
iN1NiN
b
a1N xRAdxxwxR ,
de unde
0dxxwxQxπ 1q
b
aN .
Relaţia arată că xπN este ortogonal cu orice polinom de grad < N. Se ştie că, în raport cu o funcţie
pondere w(x) definită pe un interval (a,b), există un polinom ortogonal unic; aşadar xπN este acest
polinom ortogonal.
In concluzie, nodurile xiN (abscisele Gauss) sunt rădăcinile polinomului ortogonal, definit în mod unic în
raport cu funcţia pondere w(x).
Coeficienţii AiN se vor determina apoi, rezolvând N ecuaţii, din cele 2N ale sistemului liniar.
Coeficienţii AiN se pot exprima şi prin intermediul polinomului de interpolare Lagrange. Formula
gaussiană, având gradul de valabilitate 2N-1, este exactă şi pentru funcţia
xlxf2
N,k ,
pentru care
kNkN
2
N,kkNiN
N
1i
2
N,kiN
b
a
2
N,k AxlAxlAdxxwxl
,
deci
dxxwxlAb
a
2
N,kkN .
Eroarea integrării în metoda gaussiană este
dxxwxπ
!N2
ξfR
b
a
2
N
N2
N .
Particularizăm relaţia pentru polinoamele ortogonale uzuale, obţinând formule de integrare utile
-Cebâşev ordin 1
7
N
1i
1
1 2 N2
π1i2cosf
N
πdx
x1
xf.
O funcţie MATLAB care calculează o integrală, folosind relaţia precedentă este:
function y=GCeb(f,n)
%integrare Gauss-Cebasev
k=1:n;
x=cos(2*k-1)*pi/(2*n);
y=pi/n*sum(feval(f,x));
În scriptul de test, funcţia integrand f se defineşte cu inline. De exemplu:
F=inline(x*sqrt(x))
Z=Gceb(F,5)
-Cebâşev ordin 2
1N
iπcosf
1N
iπsin
1N
πdxx1xf
N
1i
22
1
1
;
-Legendre
N
1i
iNiN
1
1xfAdxxf ,
cu
iN
2
1N
2
2
iN
iNxLN
x12A
.
Funcţia MATLAB care calculează o integrală definită Gauss-Legendre este
function y=GLeg(f,a,b,n)
%integrare Gauss-Legendre
z =Leg(n+1); %coeficienti Ln(x)
x=roots(z); %radacini Ln(x)
xx=(b+a)/2+(b-a)/2*x;
yy=feval(f,xx);
z1=Leg(n-1);
A=2*(1-x.^2)/n^2;
A=A./polyval(z1,x).^2;
y=(b-a)/2*A'*yy;
Dacă limitele de integrare sunt altele decât -1 şi 1, se face o schimbare de variabilă:
-Laguerre
N
1i
iNiN0
xxfAdxxfe ,
cu
iN
2
1N
2
iN
iNxG1N
xA
.
1
1
b
a
dtt2
ab
2
abf
2
abdxxfI
n
1k
kk t2
ab
2
abfA
2
abI
8
-Hermite
N
1i
iNiN
xxfAdxxfe
2
,
cu
iN
2
1N
1N
iNxH
π!N2A
.
Nodurile xiN reprezintă rădăcinile polinomului ortogonal corespunzător, iar coeficienţii AiN se obţin fie
prin rezolvarea sistemelor de ecuaţii liniare deduse, impunând ca formula de integrare să fie exactă
pentru polinoame până la gradul N, fie direct cu formulele de mai sus.
Metodele gaussiene au o precizie mai mare decât metodele de tip Newton-Cotes (având gradul de
valabilitate mai ridicat) şi pot fi utilizate şi pentru calculul unor integrale improprii de tipul de mai sus.
4. Integrare Romberg
Fie IN şi EN valoarea aproximativă a integralei şi estimarea erorii în metoda compusă a trapezelor.
Valoarea exactă a integralei este
NN EII ,
b,aξ,ξf12
NhE
3
N .
Pentru două valori diferite ale lui N avem
2211
NNNN EIEII ,
de unde
1
2
12
1
1
12
1
21
N
N
NN
N
N
NN
N
NN
E
E1
IIE
E
II
E
EE
,
2
2
2
1
3
1
3
1
3
2
3
2
1
3
11
2
3
22
N
N
N
N
N
abN
N
abN
ξf12
hN
ξf12
hN
E
E
1
2
,
(dacă se aproximează 21 ξfξf ). Prin urmare
2
2
2
1
NN
N
2
2
2
1
NN
N
N
N1
IIII
N
N1
IIE 12
1
12
1
,
Dacă se alege N1=2N1 atunci
3
II4I
4
11
IIII
NN2
N
NN2
N
,
de unde, în final
3
II4I
NN2
Dacă în această formulă înlocuim pe IN şi I2N cu estimările lor din formula compusă a trapezelor, atunci
se obţine pentru I estimarea din formula compusă a lui Simpson.
Estimând acum integrala cu formula compusă a lui Simpson, în care eroarea are expresia
ξf90
NhE
iv5
N ,
9
se obţine
4
2
4
1
NN
N
N
N1
IIII 12
1
,
Dacă se urmează tactica dublării numărului de puncte, atunci
15
II16
14
II4
14
II4II NN2
2
NN2
2
2
NN2
2
N
.
Notăm cu 1N1101 I,,I,I estimările integralelor calculate cu formula compusă a trapezelor
considerând N102,,2,2 intervale.
Se demonstrează uşor că 1NI poate fi obţinut prin recurenţă din 1,1NI cu formula
12
2i,1iN1N1,1N1,N
N
i2
abaf
2
abI
2
1I .
Integralele Simpson pentru N10
2,,2,2 intervale se obţin cu formulele
14
II4I
1,k1,1k
2,k
.
Analog
14
II4I
2
2,k2,1k
2
3,k
, şi
14
II4I
1j
1j,k1j,1k
1j
j,k
.
Se formează matricea inferior triunghiulară
NN2N1N
333231
2221
11
III
III
II
I
Fiecare coloană converge către I , cu atât mai rapid cu cât este situată mai la dreapta.
Pentru o coloană j , calculul iterativ este oprit în momentul în care
j,kj,1kj,k IεII
unde ε este toleranţa impusă. Calculul integralei pe această bază face obiectul algoritmului 7.4.
function R = Romberg(a, b, n, f)
% Intrări:
%a, b = intervalul de integrare
% f = funcţia de integrat
% 2n = numărul intervalelor de diviziune
% Ieşiri:
% R = matricea cu valorile integralelor
h = b-a;
R(1,1) = h(f(a)+f(b)) / 2;
l = 1
for i = 2 : n
s = R(i-1, 1);
for k = 1 : l
10
s = s + h.f(a+(k-0.5)h);
end
R(i,1) = s / 2;
l = 2*l;
p = 1;
for j = 2 : i
p = 4*p;
R(i,j)=(p*R(i,j-1)-R(i-1,j-1))/(p-1);
end
h = h / 2;
end
Metoda Romberg converge pentru orice funcţie integrabilă în sens Riemann.
5. Metoda seriei generatoare
Se utilizează cea de-a treia formulă de interpolare Newton-Gregory, considerând punctele echidistante
k10 x,,x,x .
uEupuhxfxf k0 ,
unde
0
k
00k fk
1kuf
1
ufup
,
cu eroarea interpolării
ξfku1uu!1k
huE
1k1k
.
Prin integrarea formulei se obţine
,Rduk
1kufdu
1
ufdufhdxxf
1
0
x
x
1
0
1
0
1
00
k
00
în care expresia restului este
du1k
kuξfkR
1
0
1k2k
.
Dacă se introduce notaţia
du!m
1mu1uuc
1
0m
,
formula de integrare devine
ξfhcfcfcfhdxxf1k2k
1k0
k
k010
x
x
1
0
.
Pentru calculul integralei folosim seria
m
0m
mtctC
,
care este absolut convergentă pentru 1t , întrucât 1cm .avem
du!m
1mu1uuttC
1
00m
m
1
0
u1
00m
mdut1dut
!m
1mu1uu ,
11
t1lnt1
ttC
,
care se rescrie
t1
1t1ln
t
tC
,
Dezvoltăm în serie
2
10
2
tt1tcc3
t
2
t1 ,
şi identificăm coeficienţii, rezultând relaţia de recurenţă
11m
c
2
cc
01m
m
.
Valorile coeficienţilor astfel calculaţi sunt 1, ,188
95,
720
251,
8
3,
12
5,
2
1 etc.
Metoda este utilizată şi pentru integrarea ecuaţiilor diferenţiale prin metode multipas.
6. Derivare numerică
Pentru o funcţie fC([a,b]) cunoscută numai prin valorile n10 xf,,xf,xf se cere
aproximarea derivatei αfp într-un punct b,aα .
Aproximarea derivatei se exprimă printr-o funcţională liniară discretă
αRxfαAαf i
n
0i
i
p
.
Coeficienţii αAi se determină impunând ca formula să fie exactă 0αR pentru o anumită clasă
de funcţii.
De exemplu pentru baza xu,,xu,xu n10 se obţine
n
0i
iki
p
k xuαAαu .
sistem de ecuaţii care ne permite să determinăm necunoscutele αAi .
Utilizarea bazei polinomiale n2
x,,x,x,1 conduce la sistemul liniar
n
0i
pkk
ii
n
0i
k
ii
.n:pk,α1pk1kkxαA
,1p:1k,0xαA
Funcţia f poate fi înlocuită cu polinomul de interpolare Lagrange, derivata funcţiei fiind estimată prin
derivata polinomului de interpolare.
αEαlxfαf
pn
0i
p
ii
p
,
de unde prin identificare se obţine
n:0i,αlαAp
ii .
Eroarea formulei de derivare este
12
pn02nn0
px,x,xFxxxxαEαR ,
)!1kpn(
)ξ(f)α(π
!k
!p)α(R
kp
)1kpn(
)k(p
0k
.
Formulele de aproximare obţinute pentru prima derivată, considerând punctele ix echidistante sunt de
forma
formule în 3 puncte
ξf
3
h
h2
xfxf4xf3xf
32
101
1
,
ξf
6
h
h2
xfxfxf
32
11
0
,
ξf
3
h
h2
xf3xf4xfxf
32
101
1
,
ξf
5
h
h12
xf3xf16xf36xf48xf25xf
54
21012
2
.
formule în 5 puncte
ξf
5
h
h12
xf3xf16xf36xf48xf25xf
54
21012
2
,
ξf
20
h
h12
xfxf6xf18xf10xf3xf
54
21012
1
,
ξf
30
h
h12
xfxf8xf8xfxf
54
2112
0
,
ξf
20
h
h12
xf3xf10xf18xf6xfxf
54
21012
1
,
ξf
5
h
h12
xf25xf48xf36xf16xf3xf
54
21012
2
.
Eroarea este minimă dacă derivata este calculată într-un punct central.