130 likes | 440 Views
f(x) dx = F(b) – F(a). b. . . a. 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:.
E N D
f(x) dx = F(b) – F(a) b a 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: 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.
b b b b n n n a a a a k =1 k =1 k =1 • 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 g(x) = akgk(x) cu k = 1,2,…,n 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: I = g(x) dx = ak gk(x) = akIk in care Ik = gk(x) dx Eroarea de integrare este: = R(x) dx si se cauta o minimizare a erorii.
f(x) g(x) R(x) x1 x2 x3 x4 x5 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. • Doua din cele mai importante clase de metode folosite pentru integrarea numerica sunt: • 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. • 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.
h = xk+1 – xk = cu proprietatile: Hk = 1 si Hk = Hn-k+1 b b b n n n n n a a a k =1 k =1 k =1 k =1 k =1 (-1)n-k·h dq Ak = (k-1)!(n-k)! b-a n-1 n-1 n-1 q(n) q(n) q – (k-1) q – (k-1) 0 0 a) Metodele Newton-Côtes Se considera ca aproximanta g(x) este un polinom de interpolare Lagrange, scris sub forma: yk·Lk(x) g(x) = iar intervalul [a,b] este impartit in (n-1) subintervale egale prin punctele x1 = a, x2, …, xn = b. Pasul diviziunii va fi: I = g(x)dx = [ yk Lk(x) dx ] = yk·Ak de unde rezulta ca: Ak = (b-a)·Hk in care Hk se numesc coeficientii lui Cotes. 1 (-1)n-k dq Hk = n-1 (k-1)!(n-k)! astfel integrala aproximantei g(x) se poate scrie in forma generala: I = g(x) dx = (b-a) yk·Hk(x) , k = 1, 2,…, n
I = g(x) dx = (b-a) ykHk h = xk+1 – xk = b x2 b b 1 n n 0 a a a x1 k =1 k =1 b-a n-1 h 2 I = g(x) dx = h (0.5 y1 + y2 + y3 + … + yn-1 + 0.5 yn) = 0.5 h (y1 + yn) + h yk a1) Metoda trapezelor Se particularizeaza formula de integrare numerica 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 H2 = q dq = 0.5 g(x) y1 Se obtine formula trapezului: (y1 + y2) g(x) dx = x1 x2 h Pentru a calcula g(x) dx se divide domeniul de integrare [a, b] in (n-1) domenii egale [x1, x2], [x2, x3],…, [xn-1, xn], se aplica formula trapezelor pentru fiecare 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. in care,
q =g(x) dx /4 /4 tgx tgx 0 0 b 1 1 a 0 0 UTILIZAREA COMENZILOR MATLAB PENTRU INTEGRAREA NUMERICA In Matlab pot fi calculate integrale definite de forma urmatoare: Sunt folosite doua functii:trapz(x, y), care calculeaza integrala prin metoda trapezelor, si quad(). Sa se calculeze integralele : a) q1 = ln(x+1)dx si b) q2 = 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: >>rez = quad(‘q1',0,1) APLICATIA 1 dx Rezultatul obtinut va fi: 0.3863 1 1 1 In mod analitic, ln(x+1) dx = xln(x+1) - x + ln (x+1) = 2ln2 – 1 = 0.3863 0 0 0 b) In mod similar, function y = q2(x) y = tan(x) Rezultatul obtinut va fi: 0.3466 /4 In fereastra principala se lanseaza comanda: >>rez = quad(‘q2',0,pi/4) = - ln cosx = - ln (1/ 2 ) = 0.3465 dx 0
1.3 1 b n a k =1 I = g(x) dx = h (0.5 y1 + y2 + y3 + … + yn-1 + 0.5 yn) = 0.5 h (y1 + yn) + h yk APLICATIA 2 Valorile functiei f(x) = x sunt date tabelar, pe intervalul [1, 1.3], cu pasul h = 0.05. Sa se calculeze f(x) dx si sa se estimeze eroarea. a) Cu ajutorul procedurii MATLAB >> x = [1:.05:1.3]; >> y = sqrt(x); >> z = trapz (x,y) obtinandu-se rezultatul: 0.3215 (eroare prin adaos = 1.463·10-5) b) Cu ajutorul formulei trapezelor generalizata 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) = = 0.32147 (eroare prin lipsa = 1.536·10-5) 1 Valoarea exacta a integralei este: 0.3214853
I = g(x) dx = (b-a) ykHk , particularizata pt. n = 3, 2 1 b b x3 x3 H1 = · (q-1)(q-2) dq = n 2 2 0 x1 a x1 a k =1 h 1 1 1 1 h 1 2 1 2 6 3 3 2 2 6 3 1 n-1 I = g(x) dx q(n) q – (k-1) 0 a2) Metoda Simpson Se aplica formula de integare numerica 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). 1 (-1)n-k Aplicand formula: dq Hk = n-1 (k-1)!(n-k)! I = g(x) dx = (x3 – x1) (H1y1 + H2y2 + H3y3), in care: Rezulta: H2 = - · (q-1)(q-2) dq = 2 care se inlocuiesc in integrala de mai sus. 0 2 1 H3 = · q (q-1) dq = 0 Valorile lui H1, H2, H3 se introduc in calculul integralei I, obtinandu-se formula lui Simpson: g(x) dx = (y1 + 4y2 + y3) , unde 2h = x3 – x1 Generalizand, = [(y1 + y2n+1) + 4 (y2 + y4 +…+y2n) + 2 (y3 + y5 +…+ y2n-1)]
0.125 3 1 1 1 - 0 9 - 1 0 0 APLICATIA 3 Sa se calculeze numeric integrala (4x3 + 3x2 + 2x +5) dx folosind: • Calculul analitic; • Formula trapezelor cu 9 puncte; • Formula Simpson cu 9 puncte; • Sa se arate care dintre cele trei variante este cea mai apropiata de valoarea exacta. 1 REZOLVARE a) (4x3 + 3x2 + 2x +5) dx = (x4 + x3 +x2 + 5x) = 8 0 b,c) Pasul de discretizare este: h = = 0.125. Punctele xk si valorile functiei f(xk) sunt prezentate in tabelul de mai jos: Itrapez = 0.125 (2.5 + 5.3046 + 5.75 + 6.3828 + 7.25 + 8.3984 + 9,875 + 11.7265 + 7) = 8.0234 ISimpson = [(5 + 14) + 4·(5.3046 + 6.3828 + 8.3984 + 11.7625) + 2·(5.75 + 7.25 + 9.8750] = = 7.9999 Eroarea de calcul cea mai mare ( = 0.2875%) este in cazul folosirii metodei trapezelor.
3x2 - 1 2 1 1 dn 8 [(x2 – 1)n] dxn 2n·n! 1 -1 5x3 – 3x 2 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: Pn(x) = n =0,1,2,… • Aceste polinoame au urmatoarele proprietati: • Pn(1) = 1 si Pn(-1) = (-1)n. • Satisfac conditiile de ortogonalitate: Pn(x)·Qm(x) = 0 unde Qm(x) este orice polinom de grad m <n • Polinomul Legendre Pn(x) are n radacini reale si distincte in intervalul (-1,1). • Primele cinci polinoame Legendre au expresiile: • P0(x) = 1; P1(x) = x; P2(x) = , P3(x) = , P4(x) = (35x4 – 30x2 +3)
I = f(x) dx = Ak·F(zk) + in care: Ak = 1 1 1 (b-a)2n+1 (n!)4 f(2n)() [(2n)!]3 (2n+1) 0 -1 -1 b n a k =1 d3 1 [(z2 – 1)3] 23·3! dx3 1 + 2x 1 2 Formula lui Gauss-Legendre este de forma: b-a Lk(z) dz se numesc ponderi, 2 b+a b-a iar F(zk) = f ( zk+ ) 2 2 In care zk sunt radacinile polinoamelor Legendre de grad n iar eroarea este data de: (a,b) = APLICATIA 4 Sa se calculeze integrala dx folosind formula de cuadratura Gauss-Legendre cu n = 3. REZOLVARE Se porneste de la expresia generala a polinomului Legendre, particularizandu-se pentru n = 3. P3(z) = = (5z3 – 3z) Anuland P3(z) se obtin solutiile: z1 = - 0.774597; z2 = 0; z3 = 0.774597 Pentru determinarea ponderilor Ak se aplica relatia: Ak = Lk(z) dz, particularizata pt. n = 1, 2, 3.
x1 = (- 0.774597) + = 0.11270 3 1 1+0 1+0 1+0 1 0 2 2 2 b 1-0 b+a b-a 1-0 1-0 3 2 2 2 2 2 a (z-z1)(z-z3) (z-z2)(z-z3) (z-z1)(z-z2) k =1 L3(z) = L2(z) = L1(z) = (z3-z1)(z3-z2) (z2-z1)(z2-z3) (z1-z2)(z1-z3) 1 + 2x 5 1 8 1 9 3 2 9 3 in care: Rezulta asadar: A1 = A3 = ; A2 = Se face substitutia de variabila x = z + x2 = (0) + = 0.50000 x3 = (0.774597) + = 0.88730 b-a I = f(x) dx = Ak·f(xk) =[A1f(x1) +A2f(x2) + A3f(x3)] = 1.39870 2 Calculand integrala prin metoda analitica, efectuand substitutia 2x+1 = t rezulta dx = tdt, iar integrala devine: dx = t2dt = t3 = 1.39872 1 obtinandu-se o eroare de 2·10-5.