Comments
Description
Transcript
Integrazione numerica
Integrazione numerica Gabriella Puppo Integrazione numerica Formula dei trapezi Formula di Simpson Stima numerica dell’accuratezza Funzioni di quadratura di Matlab Esempi Formula dei trapezi Per costruire una function che applichi il metodo dei trapezi, devo costruire una function con le caratteristiche seguenti: La function deve avere in input: il nome della funzione integranda, l’intervallo [a,b] di integrazione, il numero di intervalli in cui suddividere [a,b] La function deve dare in output il risultato del calcolo dell’integrale All’interno, ci deve essere un ciclo, nel quale si applica la formula composita dei trapezi function s=trapezi(fun,a,b,n) % %TRAPEZI Calcola integrali usando il metodo dei trapezi % TRAPEZI(FUN,A,B,N): Calcola l'integrale di FUN fra A e B % usando la formula dei trapezi composita e dividendo % [A,B] in N intervalli uguali % h=(b-a)/n; x=a:h:b; f=feval(fun,x); s=f(1)/2; for i=2:n s=s+f(i); end s=s+f(n+1)/2; s=s*h; Esempio Calcolo l’integrale fra 0 e di f(x) = sin(x) (Risultato: 2) >> f=inline('sin(x)'); >> trapezi(f,0,pi,10) ans = 1.9835 Se uso più punti, ottengo un’approssimazione migliore: >> trapezi(f,0,pi,20) ans = 1.9959 >> trapezi(f,0,pi,40) ans = 1.9990 Formula di Simpson Per costruire una function che applichi il metodo di Simpson devo soddisfare i requisiti seguenti: La function deve avere in input: il nome della funzione integranda, l’intervallo [a,b] di integrazione, il numero di intervalli in cui suddividere [a,b] La function deve dare in output il risultato del calcolo dell’integrale All’interno, ci deve essere un ciclo, nel quale si applica la formula composita di Simpson. function s=simpson(fun,a,b,n) % %SIMPSON Calcola integrali definiti usando il metodo di Simpson % SIMPSON(FUN,A,B,N): Calcola l'integrale di FUN fra A e B % usando la formula di Simpson composita e dividendo % [A,B] in N intervalli uguali % h=(b-a)/n; x=a:h:b; f=feval(fun,x); for i=1:n xmezzi(i)=0.5*(x(i)+x(i+1)); %Calcola i punti medi di [xj,xjp1] end fmez=feval(fun,xmezzi); s=0; for i=1:n s=s+f(i)+4*fmez(i)+f(i+1); end s=s*h/6; Esempio Calcolo l’integrale fra 0 e di f(x) = sin(x) >> f=inline('sin(x)'); Questa volta i risultati sono molto più precisi >> format long >> simpson(f,0,pi,10) ans = 2.00000678444180 >> simpson(f,0,pi,20) ans = 2.00000042309318 >> simpson(f,0,pi,40) ans = 2.00000002642876 Valutare l’accuratezza di una formula di quadratura Considero un problema di cui conosco il risultato esatto Calcolo l’integrale con una griglia di ampiezza h. Ottengo l’errore e(h) Calcolo l’integrale con una griglia di ampiezza h/2. Ottengo l’errore e(h/2) Uso la stima e(h) ~ C h^p per ricavare p Stima numerica dell’accuratezza Stimo l’accuratezza in due casi: Funzione regolare: f(x)=exp(x/10)*sin(x) su [5,5] Funzione non regolare: f(x)= abs( cos(x) ) su [0,3] Nel primo caso, otterrò la stima teorica (per h abbastanza piccolo) Nel secondo caso, otterrò una velocità di convergenza più lenta Script acc_smooth.m Imposta la funzione ed il calcolo dell’integrale esatto: % % calcola l'accuratezza del metodo dei trapezi e del metodo di Simpson % usando la funzione exp(x/10)*sin(x) su (-5,5) % f = inline('exp(x/10).*sin(x)'); exa=inline('10/101*exp(x/10)*(sin(x)-10*cos(x))'); exa_int=exa(5)-exa(-5); n=10; nmax=8; % Fissa i parametri iniziali fprintf(' n trapezi Simpson \n') continua ... Calcola l’errore su diverse griglie per il metodo dei trapezi ed il metodo di Simpson for k=1:nmax trap=trapezi(f,-5,5,n); simp=simpson(f,-5,5,n); err_trap(k)=abs(trap - exa_int); err_simp(k)=abs(simp - exa_int); n=n*2; % Raddoppia il numero di intervalli fprintf('%4.0f %12.6e %12.6e \n',n,err_trap(k),err_simp(k)) end continua ... Usa gli errori calcolati sulle diverse griglie per stimare l’accuratezza per il metodo dei trapezi ed il metodo di Simpson %stampa l'accuratezza fprintf('\n accuratezza \n') for k=2:nmax acc_trap=log(err_trap(k-1)/err_trap(k))/log(2); acc_simp=log(err_simp(k-1)/err_simp(k))/log(2); fprintf('%4.0f %12.6e %12.6e \n',k,acc_trap,acc_simp) end Ottengo i seguenti risultati: Accuratezza delle formule di quadratura per la funzione: f(x) = exp(x/10) * sin(x) >> acc_smooth n trapezi Simpson 20 6.086964e-03 1.334895e-04 40 1.621858e-03 7.938893e-06 80 4.114187e-04 4.900996e-07 160 1.032223e-04 3.053709e-08 320 2.582847e-05 1.907099e-09 640 6.458547e-06 1.191728e-10 1280 1.614726e-06 7.446599e-12 2560 4.036871e-07 4.635181e-13 accuratezza 2 1.908075e+00 3 1.978968e+00 4 1.994853e+00 5 1.998720e+00 6 1.999680e+00 7 1.999920e+00 8 1.999980e+00 4.071645e+00 4.017791e+00 4.004440e+00 4.001111e+00 4.000253e+00 4.000329e+00 4.005884e+00 Se invece integro una funzione non regolare, l’accuratezza si deteriora: Accuratezza delle formule di quadratura per la funzione: f(x) = abs( cos(x) ) sull’intervallo [0,3]. >> acc_sing n trapezi Simpson 20 2.286208e-03 2.068921e-03 40 2.123243e-03 1.472841e-03 80 5.738197e-04 8.744110e-05 160 7.787410e-05 3.487619e-05 320 6.688622e-06 8.602118e-06 640 8.123744e-06 4.534397e-06 1280 1.369862e-06 5.218913e-07 2560 4.895294e-08 1.439013e-07 accuratezza 2 1.066872e-01 4.902775e-01 3 1.887600e+00 4.074146e+00 4 2.881382e+00 1.326069e+00 5 3.541363e+00 2.019479e+00 6 -2.804357e-01 9.237810e-01 7 2.568114e+00 3.119090e+00 8 4.806491e+00 1.858669e+00 Funzioni di quadratura di Matlab Uso di base di quad e quadl Impostare le tolleranze Visualizzare la costruzione dell’integrale Funzioni quad e quadl Le functions quad e quad8 (Matlab 5) o quad e quadl (Matlab 6 e 7) sono le functions per l’integrazione numerica di Matlab. Il passo di integrazione viene determinato in maniera adattiva, in modo da soddisfare una tolleranza, che può essere fornita in input. La function quad è basata sulla formula di Simpson, mentre quad8 usa una formula di tipo Newton-Cotes di accuratezza 8. In Matlab 6 e 7, quad8 è stata sostituita dalla function quadl, che usa una formula gaussiana di accuratezza elevata con nodi di Gauss Lobatto . Sintassi La chiamata più semplice è int = quad (fun,a,b) fun è una stringa, che contiene il nome della funzione integranda. La funzione fun deve essere scritta in modo da accettare un vettore in input, e fornire un vettore in output. a e b sono gli estremi di integrazione Esempio >> f=inline('log(x).^2'); >> int = quad(f,1,3) int = 1.02917317609112 Per avere una stima del numero di operazioni che sono state effettuate, uso un secondo argomento (opzionale) in output: >> f=inline('log(x).^2'); >> [int, op] = quad(f,1,3) int = 1.02917317609112 op = 25 Quindi, per ottenere il risultato, sono state necessarie 25 valutazioni funzionali, cioè 25 chiamate della funzione f. Impostare la tolleranza Per impostare la tolleranza sulla quale è basato l’algoritmo adattivo, devo fornire un ulteriore argomento in input: >>int = quad(f,a,b, tol) Il valore di default è 1e-3*C in Matlab5 e 1e-6*C in Matlab6 e 7, dove C è una stima numerica del valore dell’integrale. Se tol viene impostato dall’esterno, quad usa tol come una tolleranza assoluta. Sta quindi all’utente scegliere un valore appropriato per il parametro tol. Tolleranza e numero operazioni Abbassando la tolleranza, aumenta il numero di operazioni: >> f=inline('x+sin(5*x)') >> [int,c]=quad(f,1,3,1e-3) int = 4.20906240989543 c= 13 La risposta esatta è: >> 4-0.2*(cos(15)-cos(5)) ans = 4.20867001966441 >> [int,c]=quad(f,1,3,1e-4) int = 4.20867153987319 c= 29 >> [int,c]=quad(f,1,3,1e-5) int = 4.20867006658840 c= 45 Cifre significative Calcolando l’integrale due volte, con due valori diversi della tolleranza, posso avere una stima del numero di cifre significative. Per esempio >> [int,c]=quad(f,1,3,1e-4) Nella stima del primo int = integrale ci sono circa 5 cifre 4.20867153987319 significative, c= perché le prime 5 cifre 29 restano invariate abbassando >> [int,c]=quad(f,1,3,1e-5) la tolleranza. Confrontare int = con il risultato esatto: 4.20867006658840 4.20867001966441 c= 45 Visualizzare la costruzione dell’integrale La chiamata quad(f, a, b, tol, trace) dove trace è un qualunque numero diverso da zero, produce un output complesso, nel quale vengono evidenziate le vicissitudini dell’algoritmo adattivo. Si ottengono 4 colonne: 1) Numero di valutazioni funzionali calcolate finora; 2) Estremo sinistro di integrazione corrente a k; 3) Passo di integrazione corrente h k; 4) Valore stimato dell’integrale fra a k e a k + h k >> int=quad(f,1,3,1e-6,1) 5 1.0000000000 2.00000000e+00 7 1.0000000000 1.00000000e+00 9 1.0000000000 5.00000000e-01 11 1.0000000000 2.50000000e-01 13 1.0000000000 1.25000000e-01 15 1.1250000000 1.25000000e-01 17 1.2500000000 2.50000000e-01 19 1.2500000000 1.25000000e-01 21 1.3750000000 1.25000000e-01 23 1.5000000000 5.00000000e-01 25 1.5000000000 2.50000000e-01 27 1.5000000000 1.25000000e-01 29 1.6250000000 1.25000000e-01 31 1.7500000000 2.50000000e-01 33 2.0000000000 1.00000000e+00 35 2.0000000000 5.00000000e-01 37 2.0000000000 2.50000000e-01 39 2.0000000000 1.25000000e-01 41 2.1250000000 1.25000000e-01 …………………………………….. 4.4267775593 1.7184979290 0.6124073532 0.1380928507 0.0313242305 0.1067683274 0.4743125490 0.1979664452 0.2763463708 1.1121040946 0.6317457320 0.3181821137 0.3135640796 0.4803951449 2.4845509621 0.7576837835 0.3130982172 0.1624283564 0.1506694146 Uso di quad e quadl Se la funzione da integrare è regolare, di solito quadl dà risultati più accurati di quad, con meno operazioni: >> [int,fcn]=quadl(f,1,3,1e-6) int = 4.20867001970245 fcn = 48 EDU>> [int,fcn]=quad(f,1,3,1e-6) int = 4.20866998933645 fcn = 61 La funzione da integrare è: f(x) = x + sin(5x) e l’integrale esatto è: 4.20867001966441 Notare che la tolleranza è la stessa, ma i due risultati hanno accuratezza diversa. Se invece la funzione integranda non è regolare... >> f=inline('abs(cos(x))'); >> [int,fcn]=quad(f,0,2,1e-6) int = 1.09070291894070 fcn = 41 >> [int,fcn]=quadl(f,0,2,1e-6) int = 1.09070190779511 fcn = 108 In questo caso, quadl richiede più operazioni a parità di tolleranza, e inoltre il risultato è meno accurato, infatti il risultato esatto è: >> exa=2-sin(2) exa = 1.09070257317432 In ogni caso, qui è meglio spezzare l’integrale su due intervalli Esercizio 1 Calcolare l’area delle regioni di piano delimitate dalle curve f(x) = exp( -(x-1)^2 ) e: g(x) = -1/8 x + 1/2 Suggerimento: Qui devo prima disegnare un grafico delle due curve, calcolare i punti di intersezione, stabilire come è formata la regione di integrazione, e infine stimare l’area calcolando gli integrali richiesti con il segno corrretto. Esercizio 2 Calcolare l’integrale di f(x) = |cos(x)| sull’intervallo [0, 6] con 5 cifre decimali significative Suggerimento: Qui la funzione integranda non è regolare, è bene spezzare il dominio di integrazione in modo da integrare su intervalli che non contengano punti di discontinuità Esercizio 3 Esercizio 4 Calcolare l’integrale della funzione f(x) = sin(1/x) su [0.1, 1] 1) Usare la function che utilizza il metodo di Simpson su una griglia uniforme. Quanti intervalli è necessario usare per ottenere un risultato con 5 cifre significative? 2) Usare la function quad di Matlab. Che tolleranza è necessario impostare per avere un risultato con 5 cifre significative? Qual è il numero di valutazioni funzionali richiesto in quel caso? Esercizio 5 Suggerimento: Suddividere l’intervallo [0,4] (per esempio) in 400 punti, e formare un vettore t = 0:4/400:4. In corrispondenza dei valori t(i), calcolare x(i) e y(i), osservando che: x(i) = integrale (da t(i) a t(i+1)) cos(s^2) ds + x(i-1) Quanto è scomodo Word per scrivere formule di matematica!