Transcript MNMSD04
CURS 4
INTEGRAREA NUMERICA
a) Notiuni introductive
Daca o functie f este continua pe [a,b] si se cunoaste o primitiva a sa F(x), atunci integrala
definita a functiei f(x) intre limitele [a,b] se poate calcula cu formula lui Leibniz-Newton:
b
f(x) dx = F(b) – F(a)
a
unde F’(x) = f(x).
Pentru o serie de functii a caror expresie analitica este cunoscuta, calculul integralei poate fi
efectuat prin aplicarea unor relatii explicite. Exista insa numeroase situatii practice cand
evaluarea unei integrale definite cu metode analitice cunoscute este extrem de complicata sau
chiar imposibil de efectuat. Astfel de situatii apar atunci cand:
• functia de integrat are o expresie foarte complicata, pentru care nu se cunoaste nici o
primitiva;
• functia de integrat se prezinta sub forma unei serii de valori numerice tabelate, expresia sa
analitica fiin necunoscuta;
• functia de integrat reprezinta un produs matriceal sau tensorial greu de evaluat analitic, asa
cum se intampla in modelarea cu elemente finite.
In aceste situatii se apeleaza la metodele numerice. Schema generala pentru a obtine o formula
de integrare numerica este urmatoarea:
se introduce o diviziune a intervalului de calcul [a,b] prin punctele xi (i = 1, 2,…, n);
se scrie functia de integrat f(x) punandu-se in evidenta aproximanta g(x) si restul R(x), adica:
f(x) = g(x) + R(x)
se integreaza relatia de mai sus termen cu termen.
Daca aproximarea lui g(x) se face prin interpolare, in mod uzual aproximanta se alege de forma
n
g(x) = akgk(x) cu k = 1,2,…,n
k =1
Este de preferat utilizarea unor polinoame de grad mic (unu, doi sau trei) pe subintervale,
ceea ce revine de fapt la integrarea unor functii spline, racordate in noduri numai prin valorile
functiei, nu si prin derivate.
Integrala I a aproximantei se evalueaza numeric pe fiecare subinterval al diviziunii:
b
b
I = g(x) dx = ak gk(x) = akIk
a
k =1
k =1
n
n
a
b
in care Ik = gk(x) dx
a
b
Eroarea de integrare este: = R(x) dx si se cauta o minimizare a erorii.
a
In figura de mai jos se ilustreaza modul de aproximare al unei functii f(x) de catre polinomul g(x)
care reproduce exact valoarea functiei date in punctele de interpolare. Valoarea adevarata a
integralei este aria de sub curba f(x), delimitata de verticalele duse prin punctele extreme x1 si
x5 ale intervalului de integrare.
f(x)
g(x)
R(x)
x1 x2
x3
x4
x5
Doua din cele mai importante clase de metode folosite pentru integrarea numerica sunt:
a) Metodele Newton-Côtes, dezvoltate pentru o baza de puncte echidistante. Aceste metode se
folosesc pentru intervale [a,b] impartite in subintervale egale, punctele de integrare fiind
echidistante (pas de integrare constant). Formulele de integrare numerica obtinute cu
ajutorul acestor metode pot fi inchise sau deschise. In cazul formulelor inchise, punctele
extreme ale intervalului de interpolare coincid cu punctele extreme ale intervalului de
integrare [a,b]. In cazul formulelor deschise, valoarea functiilor la extremitatile intervalului de
integrare [a,b] nu este necesara si deci extremitatile intervalului de interpolare [x1, xn] nu
coincid cu cele ale intervalului de integrare numerica. Dintre cele mai folosite metode din
aceasta categorie fac parte metoda trapezului, metoda Simpson, metoda Romberg.
b) Metodele Gauss-Legendre, dezvoltate pentru o baza de puncte neechidistante. Pozitia
punctelor de evaluare a functiilor de integrat decurge dintr-un proces de optimizare, in sensul
realizarii integrarii numerice cu un numar minim de puncte pentru aceeasi eroare data de
calcul. Din aceasta categorie fac parte metoda Gauss cu polinoame Legendre, met. Cebisev.
a) Metodele Newton-Côtes
Se considera ca aproximanta g(x) este un polinom de interpolare Lagrange, scris sub forma:
n
g(x) = yk·Lk(x)
k =1
iar intervalul [a,b] este impartit in (n-1) subintervale egale prin punctele x1 = a, x2, …, xn = b.
Pasul diviziunii va fi:
b
b-a
h = xk+1 – xk = n-1
b
n
n
I = g(x)dx = [ yk Lk(x) dx ] = yk·Ak
a
a
k =1
k =1
n-1
q(n)
(-1)n-k·h
dq
Ak =
(k-1)!(n-k)! 0 q – (k-1)
n-1
q(n)
1
(-1)n-k
dq
Hk =
n-1 (k-1)!(n-k)! 0 q – (k-1)
de unde rezulta ca: Ak = (b-a)·Hk
in care Hk se numesc coeficientii lui Cotes.
astfel integrala aproximantei g(x) se poate scrie in forma generala:
b
n
I = g(x) dx = (b-a) yk·Hk(x) , k = 1, 2,…, n
a
k =1
n
cu proprietatile: Hk = 1 si Hk = Hn-k+1
k =1
a1) Metoda trapezelor
b
Se particularizeaza formula de integrare numerica
n
I = g(x) dx = (b-a) ykHk
a
k =1
pentru n = 2., deci pentru un interval de integrare cu doua puncte de interpolare (k = 1 si k = 2).
Aproximanta g(x) devine o dreapta care trece prin punctele M1 (x1, y1) si M2 (x2, y2).
f(x)
y2
g(x)
y1
1
H2 = q dq = 0.5
0
Se obtine formula trapezului:
x1
x2
x2
h
g(x)
dx = (y1 + y2)
2
x1
h
b
Pentru a calcula g(x) dx se divide domeniul de integrare [a, b] in (n-1) domenii egale
a [x , x ], [x , x ],…, [x x ], se aplica formula trapezelor pentru fiecare
1
2
2
3
n-1, n
subdomeniu in parte si apoi se insumeaza rezultatele partiale.
Tinand cont de faptul ca extremele domeniului de integrare y1 si yn se iau o singura data in
calcul, in timp ce valorile functiilor pentru celelalte pozitii se iau in calcul de doua ori se obtine
formula trapezelor generalizata.
b
n
I = g(x) dx = h (0.5 y1 + y2 + y3 + … + yn-1 + 0.5 yn) = 0.5 h (y1 + yn) + h yk
a
k =1
b-a
in care, h = xk+1 – xk = n-1
UTILIZAREA COMENZILOR MATLAB PENTRU
INTEGRAREA NUMERICA
b
In Matlab pot fi calculate integrale definite de forma urmatoare: q = g(x) dx
a
Sunt folosite doua functii: trapz(x, y), care calculeaza integrala prin metoda trapezelor, si
quad().
APLICATIA 1
1
/4
Sa se calculeze integralele : a) q1= ln(x+1) dx si b) q2 = tgx dx
0
0
REZOLVARE
a) Se creeaza un fisier (spre exemplu q1.m) cu ajutorul comenzii Edit:
function y = q1(x)
y = log(x+1)
In fereastra principala se lanseaza comanda:
Rezultatul obtinut va fi: 0.3863
>> rez = quad(‘q1',0,1)
1
1
1
1
In mod analitic, ln(x+1) dx = xln(x+1) - x + ln (x+1) = 2ln2 – 1 = 0.3863
0
0
0
0
b) In mod similar,
function y = q2(x)
y = tan(x)
In fereastra principala se lanseaza comanda:
>> rez = quad(‘q2',0,pi/4)
Rezultatul obtinut va fi: 0.3466
/4
/4
tgx dx = - ln cosx = - ln (1/ 2 ) = 0.3465
0
0
APLICATIA 2
Valorile functiei f(x) = x sunt date tabelar, pe intervalul [1, 1.3], cu pasul h = 0.05.
x
1.00
f(x)
1.00
1.05
1.10
1.15
1.20
1.25
1.30
1.02470 1.04881 1.07238 1.09544 1.11803 1.14017
1.3
Sa se calculeze f(x) dx si sa se estimeze eroarea.
1
a) Cu ajutorul procedurii MATLAB
>> x = [1:.05:1.3];
>> y = sqrt(x);
obtinandu-se rezultatul: 0.3215 (eroare prin adaos = 1.463·10-5)
>> z = trapz (x,y)
b) Cu ajutorul formulei trapezelor generalizata
b
n
I = g(x) dx = h (0.5 y1 + y2 + y3 + … + yn-1 + 0.5 yn) = 0.5 h (y1 + yn) + h yk
a
k =1
1.3
I = g(x) dx = 0.05 (0.5·1 + 1.02470 + 1.04881 + 1.07238 + 1.09544 + 1.11803 + 0.5· 1.14017) =
1
= 0.32147 (eroare prin lipsa = 1.536·10-5)
Valoarea exacta a integralei este: 0.3214853
a2) Metoda Simpson
b
n
Se aplica formula de integare numerica I = g(x) dx = (b-a) ykHk , particularizata pt. n = 3,
a
k =1
deci pentru un interval de interpolare cu trei puncte (k =1, 2, 3). Aproximanta g(x) reprezinta o
parabola care trece prin punctele A(x1, y1), B(x2, y2), C(x3, y3).
n-1
q(n)
1
(-1)n-k
Aplicand formula: Hk =
dq
n-1 (k-1)!(n-k)! 0 q – (k-1)
x
3
Rezulta: I = g(x) dx = (x3 – x1) (H1y1 + H2y2 + H3y3), in care:
x1
2
1 1
1
H1 = · (q-1)(q-2) dq =
2 2
6
0
2
1 1 2
H2 = - ·
(q-1)(q-2) dq =
3
2 1
0
care se inlocuiesc in integrala de mai sus.
2
1
1
H3 = · 1 q (q-1) dq =
6
2 2
0
Valorile lui H1, H2, H3 se introduc in calculul integralei I, obtinandu-se formula lui Simpson:
x
3
g(x)
dx = h (y1 + 4y2 + y3) , unde 2h = x3 – x1
x1
3
b
Generalizand, I = g(x) dx = h [(y1 + y2n+1) + 4 (y2 + y4 +…+y2n) + 2 (y3 + y5 +…+ y2n-1)]
3
a
APLICATIA 3
1
Sa se calculeze numeric integrala (4x3 + 3x2 + 2x +5) dx folosind:
0
a) Calculul analitic;
b) Formula trapezelor cu 9 puncte;
c) Formula Simpson cu 9 puncte;
Sa se arate care dintre cele trei variante este cea mai apropiata de valoarea exacta.
1
a) (4x3 + 3x2 + 2x +5) dx = (x4 + x3 +x2 + 5x)
0
REZOLVARE
1
=8
0
1-0
b,c) Pasul de discretizare este: h =
= 0.125. Punctele xk si valorile functiei f(xk) sunt
9-1
prezentate in tabelul de mai jos:
k
1
2
3
4
5
6
7
8
9
xk
0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1
yk
5
5.3046
5.75
6.3828
7.25
8.3984
9.875
11.7265
14
Itrapez = 0.125 (2.5 + 5.3046 + 5.75 + 6.3828 + 7.25 + 8.3984 + 9,875 + 11.7265 + 7) = 8.0234
ISimpson =
0.125
[(5 + 14) + 4·(5.3046 + 6.3828 + 8.3984 + 11.7625) + 2·(5.75 + 7.25 + 9.8750] =
3
= 7.9999
Eroarea de calcul cea mai mare ( = 0.2875%) este in cazul folosirii metodei trapezelor.
b) Metode Gauss-Legendre
Aceste metode reduc eroarea de integrare printr-o alegere optimizata a punctelor de diviziune.
In aceste conditii, pasul nu mai poate fi considerat constant, iar valorile functiei la capetele
intervalului nu mai sunt folosite. Se pune deci problema alegerii unor noduri astfel incat eroarea
de integrare sa fie cat mai mica. Acest lucru se poate realiza printr-o interpolare cu un polinom
de grad mai mare, fara a mari insa numarul de noduri. Astfel, se aleg ca puncte de diviziune
radacinile unor polinoame ortogonale: Legendre, Cebisev, Hermite.
b1) Metoda Gauss cu polinoame Legendre
Se numesc polinoame Legendre expresii de forma:
1
dn [(x2 – 1)n]
Pn(x) = n
n =0,1,2,…
2 ·n! dxn
Aceste polinoame au urmatoarele proprietati:
1) Pn(1) = 1 si Pn(-1) = (-1)n.
1
2) Satisfac conditiile de ortogonalitate: Pn(x)·Qm(x) = 0 unde Qm(x) este orice polinom de
-1
grad m <n
3) Polinomul Legendre Pn(x) are n radacini reale si distincte in intervalul (-1,1).
Primele cinci polinoame Legendre au expresiile:
5x3 – 3x
1
3x2 - 1
P0(x) = 1; P1(x) = x; P2(x) =
, P3(x) =
, P4(x) = (35x4 – 30x2 +3)
2
8
2
Formula lui Gauss-Legendre este de forma:
b
b-a
I = f(x) dx = 2
a
1
n
Ak·F(zk) + in care: Ak = Lk(z) dz se numesc ponderi,
-1
k =1
iar F(zk) = f ( b-a zk + b+a )
2
2
In care zk sunt radacinile polinoamelor Legendre de grad n iar eroarea este data de:
(b-a)2n+1 (n!)4 f(2n)()
=
[(2n)!]3 (2n+1)
(a,b)
APLICATIA 4
1
Sa se calculeze integrala
dx
0 1 + 2x
folosind formula de cuadratura Gauss-Legendre cu n = 3.
REZOLVARE
Se porneste de la expresia generala a polinomului Legendre, particularizandu-se pentru n = 3.
1
d3 [(z2 – 1)3] 1
P3(z) = 3
= (5z3 – 3z)
3
2 ·3! dx
2
Anuland P3(z) se obtin solutiile: z1 = - 0.774597; z2 = 0; z3 = 0.774597
1
Pentru determinarea ponderilor Ak se aplica relatia: Ak = Lk(z) dz, particularizata pt. n = 1, 2, 3.
-1
(z-z1)(z-z3)
(z-z2)(z-z3)
L2(z) =
in care: L1(z) =
(z2-z1)(z2-z3)
(z1-z2)(z1-z3)
5
8
Rezulta asadar: A1 = A3 = ; A2 =
9
9
Se face substitutia de variabila x =
(z-z1)(z-z2)
L3(z) =
(z3-z1)(z3-z2)
b-a
b+a
z+
2
2
1-0
1+0
(- 0.774597) +
= 0.11270
2
2
1-0
1+0
x2 =
(0) +
= 0.50000
2
2
1-0
1+0
x3 =
(0.774597) +
= 0.88730
2
2
x1 =
b
I = f(x) dx = b-a
2
a
3
Ak·f(xk) =
k =1
1
[A f(x ) +A2f(x2) + A3f(x3)] = 1.39870
2 1 1
Calculand integrala prin metoda analitica, efectuand substitutia 2x+1 = t rezulta dx = tdt,
iar integrala devine:
1
3
1
1 + 2x dx = t2dt = t3
0
3
1
3
= 1.39872
1
obtinandu-se o eroare de 2·10-5.