...

Integrazione numerica

by user

on
Category: Documents
12

views

Report

Comments

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!
Fly UP