Integrali generalizzati o impropri

Integrali generalizzati o impropri
Finora si `
e considerato il caso di integrali di una funzione sufficientemente regolare su un
intervallo chiuso e limitato [a, b].
Estendiamo l’analisi al caso di integrali generalizzati o impropri per i quali succede che
la funzione integranda presenta discontinuit`
a di prima specie;
la funzione integranda ha discontinuit`
a di seconda specie (per esempio `
e illimitata in
uno degli estremi o in entrambi);
il dominio `
e illimitato.
Non esiste un modo generale per trattare queste situazioni. Piuttosto esistono alcune idee da
adottare caso per caso.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
1 / 53
Integrali di funzioni con discontinuit`a di prima specie
Sia c un punto interno ad [a, b] e sia f (x) una funzione continua e limitata in [a, c) e
(c, b] con un salto finito in c, ossia f (c + ) − f (c − ) `e finito. Si ha allora che
Z b
Z c
Z b
f (x)dx =
f (x)dx +
f (x)dx
a
a
c
Si pu`
o procedere applicando una formula di quadratura a ciascuno dei due integrali.
E’ immediato generalizzare al caso in cui ci sono pi`
u punti di discontinuit`
a di
prima specie. Anche in questo caso si decompone il dominio di integrazione in modo
che i punti di discontinuit`
a risultino estremi degli intervalli di integrazione. E’ come
applicare una formula composta.
Ci sono vari casi:
si conosce la collocazione dei punti di discontinuit`
a; basta risolvere gli integrali
sui sottointervalli;
se la posizione dei punti di discontinuit`
a non `e nota a priori, ma la funzione `e
nota analiticamente, occorre procedere a un esame preliminare della funzione
(uno studio grafico) per localizzare almeno approssimativamente i punti di
discontinuit`
a;
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
2 / 53
Integrali di funzioni con discontinuit`a di prima specie (cont.)
se la funzione non `e nota in forma esplicita o `e troppo complessa da studiare, si
pu`
o affidare a un metodo di integrazione il compito di scoprire punti di
discontinuit`
a; quando si applica un metodo adattivo, la presenza di tali punti `e
segnalata da un eccessivo ridursi del passo di integrazione.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
3 / 53
Esempio
x2 − 1
e x sin x
Z
0.5
f (x) =
1
Z
f (x)dx =
0
0
0 ≤ x < 0.5
0.5 ≤ x ≤ 1
(x 2 − 1)dx +
Z
1
e x sin xdx
0.5
Usando la formula di Simpson composta relativa a 10 sottointervalli per il secondo
integrale (il primo si risolve analiticamente), si ottiene
I ' 0.53 /3 − 0.5 + 0.737555547941281 = 0.279222214607947
.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
4 / 53
Esempio (cont.)
clear all;close all
%I CASO
f2=@(x)exp(x).*sin(x);
a=0.5;b=1;
[i,itf]=simpson(f2,a,b,10);
i=i+ 0.5^3/3-0.5;
fprintf(’I metodo con Simpson \n i=%20.13e\n’,i);
pause
%II CASO
fprintf(’Metodo adattivo\n’);
tol=1e-5;lev=0;maxlev=10;
[s1,itc,vet_val]=quad_ada(’dis1’,0,1,0,0,0,0,tol,lev,maxlev);
fprintf(’i=%g\n’,s1);
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
function y = dis1(x)
y=zeros(size(x));
ij=find(x>=0&x<0.5);
y(ij)=x(ij).^2-1;
ij=find(x>=0.5 & x<=1);
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
5 / 53
Esempio (cont.)
y(ij)=exp(x(ij)).*sin(x(ij));
end
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
6 / 53
Integrali di funzioni con discontinuit`a di seconda specie
Consideriamo il caso in cui limx→a+ f (x) = ∞; l’analisi `e simile se
limx→b− f (x) = ∞; se poi il punto di discontinuit`
a `e un punto c interno ad [a, b], ci
si pu`
o ricondurre a uno dei due casi precedenti dopo aver decomposto l’integrale
come somma di due integrali sui due sottodomini [a, c) e (c, b].
Per fissare le idee supponiamo che l’integrale sia del tipo:
Z b
g (x)
I =
dx 0 < p < 1
(x
− a)p
a
con g (x) funzione regolare nell’intervallo [a, b], limitata in valore assoluto da una
costante M.
La restrizione 0 < p < 1 assicura l’esistenza dell’integrale come limite di integrali di
Riemann su sottointervalli che tendono all’intervallo [a, b]. Infatti si ha
b
Z
|I | ≤ M lim
t→a+
t
(b − a)−p+1
1
dx
=
M
(x − a)p
−p + 1
Si vuole approssimare I a meno di una tolleranza prefissata.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
7 / 53
Un esempio di integrale di questo tipo `e dato da:
1
Z
0
x
Poich`e e `e limitata in [0, 1] e
procedere con vari metodi.
√1
x
ex
√ dx
x
√
ammette primitiva (2 x), l’integrale esiste. Si pu`
o
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
8 / 53
Metodo 1
Si decompone I come somma di due integrali I = I1 + I2 , con
a+
Z
I1 =
a
g (x)
dx
(x − a)p
Z
b
I2 =
a+
g (x)
dx
(x − a)p
con 0 < < b − a.
L’integrale I2 si integra con un metodo di quadratura usuale anche se occorre
tener conto che per → 0 il comportamento della funzione integranda tende a
peggiorare.
L’integrale I1 pu`
o essere approssimato nel seguente modo: si utilizza lo sviluppo
in serie di Taylor di punto iniziale a di g (x):
g (x) = g (a) + (x − a)g 0 (a) + ... +
(x − a)r r
(x − a)r +1 r +1
g (a) +
g (ξ(x))
r!
(r + 1)!
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
9 / 53
Metodo 1 (cont.)
Dunque
a+
Z
g (a) + (x − a)g 0 (a) + ... +
I1 =
(x−a)r
r!
(x −
a
g r (a) +
(x−a)r +1
(r +1)!
g r +1 (ξ(x))
dx
a)p
Integrando si ha
I1
=
+
g (a)
g 0 (a)
2 g 00 (a)
r g r (a)
+
+
+ ... +
1−p
2−p
2!(3 − p)
r !(r + 1 − p)
Z a+
1
r +1−p r +1
(x − a)
g (ξ(x))dx
(r + 1)! a
1−p
!
+
Si pu`
o approssimare I1 con il primo termine con un errore maggiorato da
|E1 | ≤
r +2−p
maxx∈[a,a+] |g r +1 (x)|
(r + 1)!(r + 2 − p)
Supponendo di fissare < 1, se le derivate di g sono limitate e non aumentano
troppo con r , E1 `e decrescente al crescere di r . Allora data una tolleranza tol, si
devono trovare r ed e una formula di quadratura per il calcolo di I2 tale che
|E1 | + |E2 | ≤ tol
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
10 / 53
Metodo 1 (cont.)
Se si risolve il secondo integrale con una formula composta, occorre trovare r , ed h (passo per selezionae i nodi) in modo conveniente.
Si ha un problema numerico–pratico di bilanciamento: diminuire pu`
o
migliorare l’approssimazione di I1 , ma pu`
o richiedere un h pi`
u piccolo per I2 , con
un alto numero di suddivisioni. D’altra parte aumentare r implica il calcolo di un
numero maggiore di derivate.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
11 / 53
Esempio
1
Z
0
x
e =1+x +
x2
2!
+ ... +
Z
|E1 | = |
0
xr
r!
+
x r +1 ξ
e ,
(r +1)!
ex
√ dx
x
con ξ ∈ (0, 1), si ha che
x r +1−1/2 e ξ(x)
|
(r + 1)!dx
≤
r +3/2
maxx∈[0,] |e x |
(r + 1)!(r + 3/2)
|E2 |
≤
(1 − )
h4
maxx∈[,1] |f 4 (x)|
180
ex
f (x) = √
x
ove si assume di risolvere I2 con la regola di Simpson composta e dunque h = 1−
.
2m
Si vuol risolvere l’integrale con una precisione dell’ordine di 10−5 , equidistribuendo
l’errore sul calcolo dei due integrali.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
12 / 53
Occorre determinare , r e m in modo che entrambe gli errori non superino tale
tolleranza.
Posto = 0.1 (se `e troppo piccolo, il valore maggiorante del valore assoluto della
ex
derivata quarta di √
diventa molto grande e produce un passo h piccolo e molte
x
valutazioni di funzione), si determinano m e r di conseguenza. Per r = 3 e n = 108,
si ha
|E1 | ≤
0.19/2 e 0.1
−7
= 3.2360
4!9/2
|E2 | ≤
0.95
0.1
4
3
2
9/2
−6
(1/16)e (160.1 − 320.1 + 720.1 − 1200.1 + 105)/0.1
= 4.9345
1084 · 180
(la derivata quarta della funzione `e 1/16e x (16x 4 − 32x 3 + 72x 2 − 120x + 105)/x 9/2
ed `e decrescente.)
Pertanto si calcola
1
0.1
0.12
0.13
+
+
+
1 − 1/2
2 − 1/2
2!(3 − 1/2)
3!(4 − 1/2)
I1
=
√
0.1
I2
=
2.271118
I
=
2.925303
!
= 0.654185
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
13 / 53
clear all
close all
epsilon=0.1;
tol=1e-5; %% tolleranza prefissata
%%%%%%%%%%%%%%%%%%
syms t
f=exp(t)/sqrt(t);
f4=diff(f,4);
x=linspace(epsilon,1,200);
M2=max(abs(subs(f4,x)));
%%%%%%%%%%%%%%%%%%%%%
m= ((1-epsilon)^5*M2/(180*16*tol/2))^(1/4);
n=2*(fix(m)+1);
fprintf(’numero di suddivisioni per calcolare I2=%g\n’,n);
f=inline(’exp(x)./sqrt(x)’);
[i2,itf]=simpson(f,epsilon,1,n);
fprintf(’i2=%g\n’,i2);
M1=exp(epsilon);
r=0;
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
14 / 53
(cont.)
pp=1-1/2;
i1=exp(0)/(pp);
while epsilon^(r+3/2)/(factorial(r+1)*(r+3/2))*M1>tol/2
r=r+1;
i1=i1+epsilon^r*exp(0)/(factorial(r)*(pp+r));
end;
i1=epsilon^pp*i1;
fprintf(’numero di elementi dello sviluppo di Tayor r=%g\n’,r);
fprintf(’i1=%g\n’,i1);
i=i1+i2;
fprintf(’I=%g\n’,i);
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
15 / 53
Metodo 2
b
Z
I =
a
g (x)
dx
(x − a)p
0<p<1
Si pu`
o scrivere g (x) come
g (x) = pr (x) +
(x − a)r +1 r +1
g (ξ(x))
(r + 1)!
ove pr (x) `e il polinomio di Taylor di grado r . Allora
Z b
Z b
pr (x)
g (x) − pr (x)
dx
+
dx
I = I1 + I2 =
p
(x − a)p
a
a (x − a)
Z b
Z b
pr (x)
(x − a)r +1−p r +1
=
g (ξ(x)) +
dx
(r
+
1)!
(x
− a)p
a
a
Ora I2 pu`
o essere calcolato in modo esplicito:
I2 =
r
X
k=0
g k (a)
(b − a)k+1−p
k!(k + 1 − p)
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
16 / 53
Metodo 2 (cont.)
Rb
r +1−p
L’integrale I1 = a (x−a)
g r +1 (ξ(x))dx ha una funzione integranda con una
(r +1)!
singolarit`
a su derivate di ordine pi`
u elevato e pu`
o essere approssimata con formule di
quadratura usuali. La scelta di r dipender`
a dalla formula usata. Se si usa una
formula con grado di precisione p, `e opportuno che r ≥ p.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
17 / 53
Esempio
1
Z
0
x
2
3
ex
√ dx
x
4
e ξ x4!
e = 1 + x + x2! + x3! +
Allora
Z 1
1
x
x2
x3
I2 =
( √ + √ + √ + √ )dx
x
x
2! x
3! x
0
1
1
1
1
+
+
+
= 2.914285714285714
=
1/2
2 − 1/2
2!(3 − 1/2)
3!(4 − 1/2)
R 1 ξ(x) 4−1/2
Inoltre I1 = 0 e x4!
dx.
Si pu`
o usare una formula composta del punto medio (grado di precisione 1), che non
richiede il calcoo in 0. Occorre che
(b − a)
h2
max|f 00 (x)| ≤ 10−5
6
(x)
√ 3
(calcolabile con le istruzioni simboliche di Matlab) `e
La derivata seconda di g (x)−p
x
crescente in [0, 1] e assume valore massimo in 1, ove vale 0.53871137.
Allora il
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
18 / 53
Esempio (cont.)
numero di suddivisioni dell’intervallo deve essere pari a 96. In tal modo si trova
I1 = 0.011014296962322 e infine I = 2.9252999 con 48 valutazioni di funzione.
syms t
f=exp(t);
r=3;
tol=1e-5;
ff=taylor(f,’Order’,r+1);
fun=(f-ff)/sqrt(t);
f2=diff(fun,2);
x=linspace(0,1,200);
M=max(abs(subs(f2,x)));
m=fix(sqrt(M/(6*4*tol)))+1;
n=2*m;
fnum=inline(’(exp(x)-(1+x+x.^2/2+x.^3/6))./sqrt(x)’);
[i1,itf]=puntomedio(fnum,0,1,n);
i2=2+1/(3/2)+1/(2*(3-1/2))+1/(6*(4-1/2));
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
19 / 53
Metodo 3
b
Z
I =
a
g (x)
dx
(x − a)p
0<p<1
Si opera il cambiamento di variabile (x − a) = z q , dx = qz q−1 dz.
Dunque z ha esponente z (−p)q+q−1 = z (1−p)q−1 ; se (1 − p)q = k, con k intero
positivo opportuno, allora:
(b−a)1/q
Z
g (a + z q )z k−1 dz
I =q
0
Poich`e k ≥ 1, la funzione integranda `e continua per z = 0 e si pu`
o applicare
direttamente una formula di quadratura.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
20 / 53
Esempio
Z
I =
0
1
ex
√ dx
x
Deve essere (1 − p)q = (1 − 1/2)q = k ≥ 1. Se k = 1, si ha (1 − 1/2)2 = 1, ossia
q = 2. Allora il cambio di variabile x = z 2 implica dx = 2zdz e
Z
I =2
0
1
2
ez
zdz = 2
z
1
Z
2
e z dz
0
Integrando con la regola di Cavalieri Simpson relativa a 20 sottointervalli si ottiene
I ' 2.925307.. con un errore dell’ordine di 10−5 .
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
21 / 53
Il tipo di procedimento 3 `e specifico del particolare tipo di singolarit`
a, mentre i
Metodi 1 e 2 possono essere applicati a integrali del tipo
Z
I =
b
g (x)w (x)dx
a
ove w (x) `e singolare in a; `e sufficiente che le funzioni (x − a)k w (x) siano integrabili
esplicitamente per k = 0, 1, 2, ....
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
22 / 53
Metodo 4
Un caso interessante `e quando la funzione w (x) che ha una singolarit`
a in uno degli
estremi, non cambia segno in [a, b]. Allora si possono utilizzare formule di Gauss
relative alla funzione peso w (x), facendo finire la singolarit`
a nei pesi della formula di
quadratura. Un esempio `e quello delle formule di Gauss-Chebyshev; ma il principio `e
generale. In questo caso occorre considerare polinomi ortogonali rispetto alla
funzione peso w (x):
Z 1
g (x)
√
dx
1 − x2
−1
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
23 / 53
Integrali su intervalli illimitati
∞
Z
I =
f (x)dx
a
Si assume f (x) continua e tale che limx→∞ f (x)x 1+ = 0 con > 0.
Metodo 1. Si scrive I = I1 + I2 , con
Z b
I1 =
f (x)dx
∞
Z
I2 =
a
f (x)dx
b
Si cerca b in modo che esista una funzione che maggiora la f (x) su [b, ∞) e per
la quale sia possibile calcolare l’integrale in forma esplicita, in modo da dare una
maggiorazione dell’errore. Se l’errore complessivo ammissibile vale tol, occorre
determinare b in modo che |I2 | ≤ tol/2. Poi si approssima I con I1 , che si pu`
o
calcolare con formule di quadratura usuali, in modo che l’errore commesso sia
non superiore aR tol/2.
∞
Esempio. I = 0 (cos x)2 e −x dx = 3/5
b
Z
I1 =
0
(cos x)2 e −x dx
∞
Z
I2 =
b
(cos x)2 e −x dx ≤
Z
∞
e −x dx = e −b
b
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
24 / 53
Integrali su intervalli illimitati (cont.)
Si cerca b tale che e −b = tol/2 = 10−5 /2. Da cui b ' 12.206072.
Applicando la regola di Cavalieri Simpson, poich`e max[0,b] |f 4 (x)| ' 7, per avere
una precisione pari a 10−5 /2 sono necessarie 130 suddivisioni dell’intervallo. Si
ottiene I1 = 0.59999279. Pertanto l’approssimazione ottenuta vale
I1 + I2 = 0.599994296.
Metodo 2. Come prima, ma per il secondo integrale si esegue il cambiamento di
variabile x = t −1 che trasforma l’integrale in modo che sia definito e sia
calcolabile almeno una stima dell’errore:
Z 1/a
I =
t −2 f (1/t)dt
0
Esempio.
R∞
1
sin(x)
dx
x2
=
Rc
sin(x)
dx
1
x2
+
R∞
c
sin(x)
2 dx
Rx∞ sin(x)
dx
c
x2
R 1/c
Con il cambiamento di variabile si ottiene
= 0 sin 1/tdt. Se
cR = 1000, con la formula di Simpson composta si ottiene che
c sin(x)
dx =' 0.504179. per cui si ha un errore dell’ordine di 10−3 .
1
x2
Metodo 3. Si usano formule di Gauss-Laguerre o Gauss-Hermite.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
25 / 53
Integrali multipli
Sia Ω un dominio limitato di R2 con frontiera ∂Ω sufficientemente regolare. Si
considera il caso di calcolo approssimato dell’integrale doppio
Z
f (x, y )dxdy
Ω
con f (x, y ) funzione continua in Ω ∪ ∂Ω. Il problema viene considerato in alcuni casi
particolari:
dominio rettangolare
dominio normale
dominio poligonale decomponibile in triangoli
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
26 / 53
Integrale doppio su un dominio rettangolare
Sia Ω = [a, b] × [c, d] = {(x, y ), a ≤ x ≤ b, c ≤ y ≤ d} e si vuole calcolare
Z bZ d
f (x, y )dydx
I =
a
c
L’osservazione cruciale per risolvere il problema `e la formula di riduzione degli
integrali doppi, che riporta l’integrale al caso monodimensionale:
Z b
Z bZ d
f (x, y )dydx =
F (x)dx
a
c
a
ove la funzione F (x) `e a sua volta calcolabile con una formula di quadratura in una
dimensione
Z d
F (x) =
f (x, y )dy
c
Rb
Si approssima l’integrale a F (x)dx con una formula composta che usa nx
suddivisioni di [a, b] di ampiezza hx = b−a
:
nx
I ' hx
nx
X
i=0
V.Ruggiero (UNIFE)
ai F (xi )
Analisi Numerica II - a.a. 2012-’13
March 28, 2015
27 / 53
Integrale doppio su un dominio rettangolare (cont.)
La valutazione di F (x) in un nodo xi , i = 0, .., nx , richiede l’uso di una formula
composta basata su una suddivisione di [c, d] in ny sottointervalli di ampiezza
hy = d−c
:
ny
F (xi ) ' hy
ny
X
bj f (xi , yj )
j=0
Pertanto la formula di quadratura risultante `e:
I ' hx hy
ny
nx X
X
ai bj f (xi , yj )
i=0 j=0
che comporta la valutazione della funzione f (x, y ) negli (nx + 1)(ny + 1) punti di una
griglia sovrapposta al dominio Ω.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
28 / 53
Formula di Simpson composta bidimensionale
Consideriamo il caso speciale in cui si usa una formula di Simpson composta per
entrambe le dimensioni (si pu`
o usare qualsiasi altra formula composta). Si assume
che nx = 2mx , ny = 2my . Allora
Z d
F (x) =
f (x, y )dy =
c
hy
3

f (x, y0 ) + 2
my −1
X
m
f (x, y2j ) + 4
j=1
y
X

f (x, y2j−1 ) + f (x, yny ) −
(d − c)hy4 fy4 (x, µ)
180
j=1
con µ ∈ [c, d]. Pertanto si ha
Z bZ d
I
f (x, y )dydx =
=
a
c
hy
Z b
3
a
f (x, y0 )dx +
m
y Z
4hy X b
+
3
j=1
a
f (x, y2j−1 )dx +
hy
Z b
3
a
2hy
3
my −1 Z
X
b
j=1
a
f (x, y2j )dx +
f (x, yny )dx +
(d − c)hy4 Z b 4
fy (x, µ)dx
180
a
−
Adesso si approssima ogni integrale con la regola di Simpson composta; dunque per
ogni j = 0, 1, .., ny si ha:
Z b
a
f (x, yj )dx =
hx
3

f (x0 , yj ) + 2
mX
x −1
i=1
f (x2i , yj ) + 4
mx
X
i=1

f (x2i−1 , yj ) + f (xnx , yj ) −
(b − a)hx4 fx4 (ξj , yj )
180
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
29 / 53
Formula di Simpson composta bidimensionale (cont.)
per ξj ∈ [a, b].
La formula risultante `e:
I

hx hy
=
f (x0 , y0 ) + 2
9
2
X
f (x0 , y2j ) + 4
j=1
m
4
y
X
f (x2i−1 , y0 ) + f (xnx , y0 )+
f (x2i , y2j ) + 8
i=1
my −1 mx
X X
j=1
m
f (x0 , y2j−1 ) + 8
j=1
+
mx
X
i=1
my −1 mx −1
X X
j=1
+
f (x2i , y0 ) + 4
i=1
my −1
+
mX
x −1
y mX
x −1
X
j=1
f (x0 , yny ) + 2
mX
x −1
f (x2i−1 , y2j ) + 2
my −1
X
i=1
m
f (x2i , y2j−1 ) + 16
i=1
f (x2i , yny ) + 4
i=1
mx
y X
X
m
f (x2i−1 , y2j−1 ) + 4
j=1 i=1
mx
X
f (xnx , y2j ) +
j=1
y
X
f (xnx , y2j−1 ) +
j=1

f (x2i−1 , yny ) + f (xnx , yny ) + E [f ]
i=1
Il termine errore `e dato da
E [f ]
=
−
−
hy (b − a)hx4
(d −
180
f 4 (ξ0 , y0 ) + 2
my −1
X
x
540
c)hy4

j=1
4
fx (ξ2j , y2j ) + 4

my
X
4
4
fx (ξ2j−1 , y2j−1 ) + fx (ξny , yny ) +
j=1
Z b
4
fy (x, µ)dx
a
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
30 / 53
Applicando ripetutamente il teorema del valor medio discreto e il teorema integrale
della media, si ha
E [f ]
=
=
(d − c)(b − a)hy4 4
4
6my fx (η, ζ) −
fy (ψ, µ)
540
180
(b − a)(d − c)
4 4
4 4
−
hx fx (η, ζ) + hy fy (ψ, µ)
180
−
hy (b − a)hx4 (η, ζ), (ψ, µ) ∈ Ω.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
31 / 53
Esempio
1
Z
0
Z
1
cos(10(x 2 + y )) exp(y 2 + x)dxdy
−1
In questo caso a = 0, b = 1, c = −1, d = 1. Si usa una regola di Simpson composta
relativa a 50 suddivisioni nell’intervallo [0, 1] e 50 suddivisioni nell’intervallo [−1, 1].
Segue il codice Matlab.
i=simpson1(’cal_int’, 0,1,50);
fprintf(’i=%g \n’,i);
%
function y=cal_int(x);
c=-1; d=1;
y=simpson2(’f11’,c,d,50,x);
%
function z=f11(x,y);
z=cos(10*(x.^2+y)).*exp(y.^2+x);
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
32 / 53
Esempio (cont.)
function [i,itf]=simpson1(f,a,b,n);
if mod(n,2)==0
h=(b-a)/n;
x=linspace(a,b,n+1);
for i=1:n+1
f1(i)=feval(f,x(i));
end;
i=(f1(1)+f1(n+1)+4*sum(f1(2:2:n))+...
2*sum(f1(3:2:n-1)))*h/3;
itf=n+1;
else
error(’n. di sottointervalli dispari’);
end;
%
function [i,itf]=simpson2(f,a,b,n,t);
if mod(n,2)==0
h=(b-a)/n;
x=linspace(a,b,n+1);
for i=1:n+1
f1(i)=feval(f,t,x(i));
V.Ruggiero (UNIFE)
Analisi Numerica II - a.a. 2012-’13
March 28, 2015
33 / 53
Esempio (cont.)
end;
i=(f1(1)+f1(n+1)+4*sum(f1(2:2:n))+...
2*sum(f1(3:2:n-1)))*h/3;
itf=n+1;
else
error(’n. di sottointervalli dispari’);
end;
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
34 / 53
Per approssimare un integrale doppio su un dominio rettangolare si possono usare
anche formule di Romberg, formule di quadratura adattiva, formule di Gauss.
Esempio. Calcolo del valore approssimato di
Z 2 Z 1.5
log(x + 2y )dxdy
1.4
1
usando una formula di Gauss-Legendre relativa a tre punti
u0 = v0 = −0.7745966692, u1 = v1 = 0, u2 = v2 = 0.7745966692 e pesi
w0 = w2 = 0.5555555556, w1 = 0.88888888889.
Occorre prima di tutto trasformare gli intervalli [1.4, 2] e [1, 1.5] a [−1, 1]:
2x − 1.4 − 2.0
2.0 − 1.4
2y − 1.0 − 1.5
v=
1.5 − 1.0
u=
x = 0.3u + 1.7
y = 0.25v + 1.25
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
35 / 53
(cont.)
Z 1 Z 1
Z 2 Z 1.5
log(0.3u + 0.5v + 4.2)dvdu
log(x + 2y )dxdy = 0.075
1.4
−1
1
−1
Z 1
(w0 log(0.3u + 0.5(−0.7745966692) + 4.2) + w1 log(0.3u + 4.2) +
' 0.075(
1
+w2 log(0.3u + 0.5(0.7745966692) + 4.2))du
2
' 0.075(w0 log(0.3u0 + 0.5(−0.7745966692) + 4.2) + w0 w1 log(0.3u1 + 0.5(−0.7745966692) + 4.2) +
+w0 w2 log(0.3u2 + 0.5(−0.7745966692) + 4.2) + w1 w0 log(0.3u0 + 4.2) +
2
+w1 log(0.3u1 + 4.2) + w1 w2 log(0.3u2 + 4.2) +
+w2 w0 log(0.3u0 + 0.5(0.7745966692) + 4.2) + w2 w1 log(0.3u1 + 0.5(0.7745966692) + 4.2) +
2
+w2 log(0.3u2 + 0.5(0.7745966692) + 4.2)) =
= 0.4295545313
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
36 / 53
(cont.)
close all
clear all
F=@(u,v)log(0.3*u+0.5*v+4.2);
nodi=[-0.7745966692; 0; 0.7745966692];
w=[0.55555555556; 0.88888888889;0.55555555556 ];
[X,Y]=meshgrid(nodi,nodi); % griglia di nodi
Z=F(X,Y);
W=w*w’;
%griglia di pesi
I=0.075*sum(sum(Z.*W));
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
37 / 53
Integrale doppio su un dominio normale
La tecnica discussa pu`
o essere facilmente estesa all’approssimazione di integrali doppi
della forma
Z b Z d(x)
f (x, y )dxdy
a
c(x)
con Ω = {(x, y ), a ≤ x ≤ b, c(x) ≤ y ≤ d(x)} (dominio normale all’asse x). In
questo caso la formula di riduzione `e
Z
b
Z
d(x)
b
Z
f (x, y )dxdy =
a
con F (x) =
R d(x)
c(x)
c(x)
f (x, y )dy . In questo caso hx =
F (x)dx
a
b−a
nx
mentre l’intervallo di
suddivisione lungo l’altra dimensione dipende da x ed `e hy (x) =
d(x)−c(x)
.
ny
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
38 / 53
Per esempio nel caso di applicazione di una formula di Simpson semplice
(hx = b−a
, hy (x) = d(x)−c(x)
) si ha:
2
2
Z b Z d(x)
I
=
f (x, y )dydx =
a
'
+
+
hx
c(x)
hy (a)
3
3
4hy (a + hx )
3
hy (b)
3
Z b
hy (x)
3
a
f (x, c(x)) + 4f (x, c(x) + hy (x)) + f (x, d(x)) dx =
f (a, c(a)) + 4f (a, c(a) + hy (a)) + f (a, d(a)) +
f (a + hx , c(a + hx )) + 4f (a + hx , c(a + hx ) + hy (a + hx )) + f (a + hx , d(a + hx )) +
f (b, c(b)) + 4f (b, c(b) + hy (b)) + f (b, d(b))
Lo stesso approccio si pu`
o usare per domini normali rispetto alla variabile y :
Z
d
Z
b(y )
f (x, y )dxdy
c
a(y )
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
39 / 53
Esempio
1
Z
Z √1−(x−1)2
√
0
1−
cos(10(x 2 + y )) exp(y 2 + x)dydx
1−x 2
dominio normale
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
V.Ruggiero (UNIFE)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Analisi Numerica II - a.a. 2012-’13
1
March 28, 2015
40 / 53
Esempio (cont.)
Si fornisce il codice Matlab relativo. In questo caso si applica una formula di Simpson
composta relativa a 100 suddivisioni per quanto riguarda l’intervallo [0, 1] e una
formula di Simpson composta relativa a 50 suddivisioni lungo la dimensione y .
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
41 / 53
s1=simpson1(’int1’,0,1,100);
fprintf(’s1=%g \n’,s1);
%
function y=int1(x);
y=simpson2(’f11’,lim1(x),lim2(x),50,x);
%
function y=lim1(t);
y=1-sqrt(1-t.^2);
%
function y=lim2(t);
y=sqrt(1-(t-1).^2);
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
42 / 53
Estensione a integrali tripli
La stessa tecnica vista per integrali doppi su domini rettangolari o normali pu`
o essere
usata per integrali tripli:
Z
b
Z
d(x)
Z
v (x,y )
b
Z
f (x, y , z)dzdydx =
a
c(x)
Z
F (x)dx
u(x,y )
d(x)
Z
a
v (x,y )
Z
d(x)
f (x, y , z)dzdy =
F (x) =
c(x)
u(x,y )
F (x, y )dy
c(x)
Z
v (x,y )
F (x, y ) =
f (x, y , z)dz
u(x,y )
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
43 / 53
Integrale doppio su un poligono convesso
Se Ω `e un poligono convesso, un altro modo per approssimare un integrale doppio
consiste nel costruire una triangolazione del dominio e nel calcolare una
approssimazione dell’integrale su ciascun triangolo.
Se si considera come approssimazione dell’integrale sul triangolo la formula
!
3
X
(xi , yi )
IT = AT f
3
i=1
ove AT `e l’area di un triangolo
e l’indice i corrisponde all’indice di uno degli
P
estremi del triangolo e 3i=1 (xi 3,yi ) `e il baricentro del triangolo T , si ottiene la
formula del punto medio composta:


3
X
X
(x
,
y
)
ih
ih

I '
ATh f 
3
i =1
Th
h
(Area del prisma di base triangolare e altezza f
P
(xi ,yi )
3
h
h
ih =1
3
)
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
44 / 53
Integrale doppio su un poligono convesso (cont.)
Se si considera come approssimazione dell’integrale sul triangolo la formula
IT =
3
AT X
f (xi , yi )
3 i=1
si ottiene la formula del trapezio composta:
I '
3
X
1X
A Th
f (xih , yih )
3
i =1
Th
h
(Volume sotteso dal piano passante per i tre punti (xi , yi , f (xi , yi )))
Nel seguito si riporta una function che nel caso dell’esempio precedente esegue una
triangolazione del dominio (triangolazione di Delaunay) e implementa la formula del
trapezio composta.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
45 / 53
Integrale doppio su un poligono convesso (cont.)
function z=trireg(dx,dy);
X=[];
Y=[];
for x=0:dx:1
c=lim1(x);
d=lim2(x);
Ynew=linspace(c,d,floor((d-c)/dy+1));
Xnew=x*ones(size(Ynew));
X=[X,Xnew];
Y=[Y,Ynew];
end;
tri=delaunay(X,Y);
plot(x,lim1(x),’r’,x,lim2(x),’r’);
hold on
Z=zeros(size(X));
trimesh(tri,X,Y,Z);
nt=length(tri);
s=0;
At=abs(((X(tri(:,2))-X(tri(:,1))).*(...
Y(tri(:,3))-Y(tri(:,1)))-...
V.Ruggiero (UNIFE)
Analisi Numerica II - a.a. 2012-’13
March 28, 2015
46 / 53
Integrale doppio su un poligono convesso (cont.)
(X(tri(:,3))-X(tri(:,1))).*...
(Y(tri(:,2))-Y(tri(:,1)))))/2;
s=sum(At.*(f11(X(tri(:,1)),Y(tri(:,1)))+ ...
f11(X(tri(:,2)),Y(tri(:,2)))+f11(X(tri(:,3)),Y(tri(:,3)))));
s=s/3;
figure
Z=f11(X,Y);
trimesh(tri,X,Y,Z);
z=s;
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
47 / 53
Metodi di integrazione Montecarlo
I metodi Montecarlo sono un efficace strumento di calcolo approssimato di integrali
multidimensionali quando la dimensione spaziale del dominio diventa grande. La
scelta dei nodi in questo caso usa un approccio statistico.
Richiami. Si dice variabile casuale X ∈ R n il processo di assegnare una n-upla di reali
X1 (ζ), ..., Xn (ζ) ad ogni risultato ζ di un esperimento. Fissato un vettore x ∈ R n , la
probabilit`
a P(X ≤ x) dell’evento casuale {X 1 ≤ x1 , ..., Xn ≤ xn } `e data da
Z x1
Z xn
P(X ≤ x) =
...
g (X1 , .., Xn )dX1 ....dXn
−∞
−∞
ove g (X1 , .., Xn ) `e la densit`
a di probabilit`
a della variabile casuale X ed `e tale che
Z
g (X1 , .., Xn ) ≥ 0
g (X1 , .., Xn )dX = 1
Rn
Data una funzione φ(X ) della variabile casuale X si dice media (valore atteso) di
φ(X )
Z
g (X1 , .., Xn )φ(X1 , ..., Xn )dX
Rn
V.Ruggiero (UNIFE)
Analisi Numerica II - a.a. 2012-’13
March 28, 2015
48 / 53
Metodi di integrazione Montecarlo (cont.)
L’idea `e quella di interpretare l’integrale come valore medio statistico:
Z
Z
|Ω|−1 χΩ (x)f (x)dx = |Ω|µ(f )
f (x)dx = |Ω|
Ω
Rn
ove |Ω| indica il volume n-dimensionale di Ω, χΩ (x) `e la funzione caratteristica di Ω e
dunque
vale 1 per gli x ∈ Ω e 0 altrove.
R
−1
|Ω|
χΩ (x)f (x)dx si interpreta come valor medio µ(f ) della funzione f (X ) ove
n
R
X `e una variabile casuale con densit`
a di probabilit`
a uniforme |Ω|−1 χΩ su R n .
Il calcolo numerico del valor medio µ(f ) viene effettuato prendendo N campioni
indipendenti x1 , ..., xN ∈ R n con densit`
a di probabilit`
a |Ω|−1 χΩ e valutando la media
aritmetica:
N
1 X
f (xi )χΩ (xi ) = IN (f )
f¯N =
N i=1
I campioni x1 , ..., xN si possono pensare come le realizzazioni di N variabili casuali
X1 , ..., XN indipendenti tra loro con la medesima densit`
a di probabilit`
a |Ω|−1 χΩ . Per
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
49 / 53
Metodi di integrazione Montecarlo (cont.)
tale successione di valori la legge dei grandi numeri assicura con probabilit`
a 1 la
convergenza della media aritmetica f¯N al valor medio µ(f ) per N → ∞.
Nella pratica la successione dei campioni x1 , ..., xN viene costruita in modo
deterministico da un generatore di numeri casuali, dando luogo alle formule di
integrazione numerica pseudo–casuali.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
50 / 53
L’errore di quadratura EN (f ) = µ(f ) − f¯N si pu`
o caratterizzare mediante la varianza
p
σ(IN (f )) = µ(IN (f ) − µ(f ))2
Di nuovo interpretando f come una funzione della variabile casuale X distribuita con
densit`
a di probabilit`
a uniforme e varianza σ(f ) si ha
σ(f )
σ(IN (f )) = √
N
da cui risulta per N → ∞ un ordine di convergenza O(N −1/2 ) della stima statistica
dell’errore. Tale ordine di convergenza non dipende dalla dimensione n del
dominio di integrazione e prescinde dalla regolarit`
a della f e costituisce la
propriet`
a pi`
u significativa del metodo Montecarlo.
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
51 / 53
Il seguente codice utilizza il metodo di Montecarlo per il calcolo dell’integrale
dell’esempio precedente.
N=500000;
rand(’state’,994243);
x=rand(N,1);y=rand(N,1);
ind=find(1-sqrt(1-x.^2)<=y & y<=sqrt(1-(x-1).^2));
x=x(ind);
y=y(ind);
I=sum(f11(x,y))/N;
fprintf(’I=%g\n’,I);
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
52 / 53
Calcolo dell’area del cerchio con il metodo di Montecarlo
N=input(’n. punti=’);
x=rand(N,1);
y=rand(N,1);
it=0;
for i=1:N
if y(i)<=sqrt(1-x(i)^2)
it=it+1
end;
end;
a=it/N*4
Analisi Numerica II - a.a. 2012-’13
V.Ruggiero (UNIFE)
March 28, 2015
53 / 53

Download Report