Integrazione numerica

Gabriella Puppo

Integrazione numerica

Formula dei trapeziFormula di SimpsonStima numerica dell’accuratezzaFunzioni di quadratura di MatlabEsempi

Formula dei trapezi

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

Per costruire una function che applichi il metodo dei trapezi, devo costruire una function con le caratteristiche seguenti:

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);ends=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

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.

Per costruire una function che applichi il metodo di Simpson devo soddisfare i requisiti seguenti:

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]endfmez=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

Funzione regolare: f(x)=exp(x/10)*sin(x) su [-5,5]

Funzione non regolare: f(x)= abs( cos(x) ) su [0,3]

Stimo l’accuratezza in due casi:

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

%% 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 inizialifprintf(' n trapezi Simpson \n')

Imposta la funzione ed il calcolo dell’integrale esatto:

continua ...

Calcola l’errore su diverse griglie per il metodo dei trapezi ed ilmetodo 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’accuratezzaper il metodo dei trapezi ed il metodo di Simpson

%stampa l'accuratezzafprintf('\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: >> 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-101280 1.614726e-06 7.446599e-122560 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)

Se invece integro una funzione non regolare, l’accuratezza sideteriora:

>> 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-061280 1.369862e-06 5.218913e-072560 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].

Funzioni di quadratura di Matlab

Uso di base di quad e quadlImpostare le tolleranzeVisualizzare 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 integrazioneEsempio

>> 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.02917317609112op = 25

Quindi, per ottenere il risultato, sono state necessarie 25 valutazionifunzionali, 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.

Abbassando la tolleranza, aumenta il numero di operazioni:

>> f=inline('x+sin(5*x)')>> [int,c]=quad(f,1,3,1e-3)int = 4.20906240989543c = 13

>> [int,c]=quad(f,1,3,1e-4)int = 4.20867153987319c = 29>> [int,c]=quad(f,1,3,1e-5)int = 4.20867006658840c = 45

Tolleranza e numero operazioni

>> 4-0.2*(cos(15)-cos(5))ans = 4.20867001966441

La risposta esatta è:

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.20867153987319c = 29>> [int,c]=quad(f,1,3,1e-5)int = 4.20867006658840c = 45

Nella stima del primo integrale ci sono circa 5 cifre significative, perché le prime 5 cifre restano invariate abbassandola tolleranza. Confrontare con il risultato esatto: 4.20867001966441

Visualizzare la costruzione dell’integrale

La chiamata

quad(f, a, b, tol, trace)

dove trace è un qualunque numero diverso da zero, produce unoutput 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 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……………………………………..

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.20867001970245fcn = 48EDU>> [int,fcn]=quad(f,1,3,1e-6)int = 4.20866998933645fcn = 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.09070291894070fcn = 41>> [int,fcn]=quadl(f,0,2,1e-6)int = 1.09070190779511fcn = 108

>> exa=2-sin(2)exa = 1.09070257317432

In questo caso, quadl richiede più operazioni a parità di tolleranza, e inoltre il risultato è meno accurato, infatti il risultato esatto è:

In ogni caso, qui è meglio spezzare l’integrale su due intervalli

Esercizio 1

Calcolare l’area delle regioni di piano delimitate dalle curvef(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 400punti, 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!