1 / 31

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.

kimberly-weaver
Download Presentation

Integrazione numerica

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Integrazione numerica Gabriella Puppo

  2. Integrazione numerica • Formula dei trapezi • Formula di Simpson • Stima numerica dell’accuratezza • Funzioni di quadratura di Matlab • Esempi

  3. 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

  4. 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;

  5. 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

  6. 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.

  7. 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;

  8. 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

  9. 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

  10. 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

  11. 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 ...

  12. 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 ...

  13. 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

  14. Ottengo i seguenti risultati: >> 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 4.071645e+00 3 1.978968e+00 4.017791e+00 4 1.994853e+00 4.004440e+00 5 1.998720e+00 4.001111e+00 6 1.999680e+00 4.000253e+00 7 1.999920e+00 4.000329e+00 8 1.999980e+00 4.005884e+00 Accuratezza delle formule di quadratura per la funzione: f(x) = exp(x/10) * sin(x)

  15. Se invece integro una funzione non regolare, l’accuratezza si deteriora: >> 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 Accuratezza delle formule di quadratura per la funzione: f(x) = abs( cos(x) ) sull’intervallo [0,3].

  16. Funzioni di quadratura di Matlab • Uso di base di quad e quadl • Impostare le tolleranze • Visualizzare la costruzione dell’integrale

  17. 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 .

  18. 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

  19. 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.

  20. 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.

  21. Tolleranza e numero operazioni Abbassando la tolleranza, aumenta il numero di operazioni: >> [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 >> 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

  22. 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) int = 4.20867153987319 c = 29 >> [int,c]=quad(f,1,3,1e-5) int = 4.20867006658840 c = 45 Nella stima del primo integrale ci sono circa 5 cifre significative, perché le prime 5 cifre restano invariate abbassando la tolleranza. Confrontare con il risultato esatto: 4.20867001966441

  23. 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

  24. >> int=quad(f,1,3,1e-6,1) 5 1.0000000000 2.00000000e+00 4.4267775593 7 1.0000000000 1.00000000e+00 1.7184979290 9 1.0000000000 5.00000000e-01 0.6124073532 11 1.0000000000 2.50000000e-01 0.1380928507 13 1.0000000000 1.25000000e-01 0.0313242305 15 1.1250000000 1.25000000e-01 0.1067683274 17 1.2500000000 2.50000000e-01 0.4743125490 19 1.2500000000 1.25000000e-01 0.1979664452 21 1.3750000000 1.25000000e-01 0.2763463708 23 1.5000000000 5.00000000e-01 1.1121040946 25 1.5000000000 2.50000000e-01 0.6317457320 27 1.5000000000 1.25000000e-01 0.3181821137 29 1.6250000000 1.25000000e-01 0.3135640796 31 1.7500000000 2.50000000e-01 0.4803951449 33 2.0000000000 1.00000000e+00 2.4845509621 35 2.0000000000 5.00000000e-01 0.7576837835 37 2.0000000000 2.50000000e-01 0.3130982172 39 2.0000000000 1.25000000e-01 0.1624283564 41 2.1250000000 1.25000000e-01 0.1506694146 ……………………………………..

  25. Uso di quad e quadl Se la funzione da integrare è regolare, di solito quadl dà risultati più accurati di quad, con meno operazioni: La funzione da integrare è: f(x) = x + sin(5x) e l’integrale esatto è: 4.20867001966441 >> [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 Notare che la tolleranza è la stessa, ma i due risultati hanno accuratezza diversa.

  26. 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

  27. 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.

  28. 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à

  29. Esercizio 3

  30. 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?

  31. 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!

More Related