Matematica e statistica
Versione didascalica: parte 3
•
•
•
Sito web del corso
http://www.labmat.it
Docente: Prof. Sergio Invernizzi, Università di Trieste
e-mail: [email protected]
2.8. Il problema dell’area
Supponiamo che una funzione sia positiva (o zero) f(x)  0 per tutti
i valori x di un intervallo [a, b], con a < b.
L’integrale della funzione f da a a b è l’area della parte di piano
compresa fra l’asse X ed il grafico di f, entro le ascisse a e b
Simbologia
area
Esempi:
(
)= 
b
f ( x)dx
a
1
 (2  x)dx  2.5
(2b )(2 a )
 (2  x)dx  2  (b  a)
0
b
a
2.9. Calcolo numerico degli integrali
I metodi per il calcolo di integrali che qui trattiamo sono:
• Metodi di interpolazione: il metodo dei rettangoli, ed il metodo
dei trapezi, basati sulla interpolazione di Lagrange;
• Metodi probabilistici: il metodo di Monte-Carlo, un metodo
molto generale basato sulla simulazione di variabili aleatorie;
• Metodi esatti (o metodi formali), basati sul Teorema
Fondamentale del Calcolo.
2.9.1. Metodi di interpolazione
Si supponga di conoscere una tabulazione della funzione f(x)  0 a
passo costante h sull’intervallo [a, b], con a < b.
a  x1  x2  ...  xn 1  b
h  (b  a ) / n
xk  x1  (k  1)h
k  1, 2,..., n  4
yk  f ( xk )
2.9.2/3. Metodi dei rettangoli e dei trapezi
La funzione può essere interpolata con la interpolazione costante o
con la interpolazione lineare

b
a
f ( x)dx può essere approssimato con le aree verdi
• Sommando aree di rettangoli (metodo dei rettangoli), oppure
• Sommando aree di trapezi (metodo dei trapezi)
Regole
a  x1  x2  ...  xn 1  b
h  (b  a ) / n
xk  x1  (k  1)h
k  1, 2,..., n
yk  f ( xk )
Regola dei rettangoli:

b
a
f ( x)dx  h { y1  y2  ...  yn  yn1}
qua no
 h   k 1 yk
n
Regola dei trapezi:
Qui salta
l’ultimo punto

b
a
f ( x)dx  h { 12 y1  [ y2  ...  yn ]  12 yn1}
 h { ( y1  yn 1 )   k 2 yk }
1
2
n
La regola dei rettangoli qui definita può
essere detta dei rettangoli destri, in quanto
tali rettangoli stanno a destra della ascissa in
cui sono calcolate le loro altezze.
Chi studia può ricavare le regola dei rettangoli sinistri, che stanno
a sinistra della ascissa in cui sono calcolate le loro altezze: viene
conseguentemente saltato il primo punto:

b
a
f ( x)dx  h { y1  y2  ...  yn  yn1}
 h   k 2 yk
n 1
In questo corso utilizziamo di default i rettangoli destri.
Cenno storico
Etimologia di “trapezio”
trapezio = banco, tavolo
τραπεζα = banca (Τραπεζα τησ Ελλαδοσ, Banco di Napoli,
Banco Monte dei Paschi di Siena)
Regola dei rettangoli su R

5
1
2(1  e
 15 ( x  4)2
f ( x)  2(1  e
>
>
>
>
>
>
>
>
a 1
)dx
 15 ( x  4)
2
)
b5
n  100
f <- function(x) 2*(1+exp((-1/5)*(x-4)^2))
a <- 1
b <- 5
a, x2 , x3 ,..., xn , b
n <- 100
h <- (b-a)/n
xtab <- a + h*c(0:n-1)
h*(sum(f(xtab))) -> integrale
integrale
Regola dei trapezi su R

5
1
2(1  e
 15 ( x  4)2
f ( x)  2(1  e
>
>
>
>
>
>
>
>
a 1
)dx
 15 ( x  4)2
)
b5
n  100
f <- function(x) 2*(1+exp((-1/5)*(x-4)^2))
a <- 1
a, x2 , x3 ,..., xn , b
b <- 5
n <- 100
h <- (b-a)/n
xtab <- a + h*c(0:n)
h*(sum(f(xtab))-(f(a)+f(b))/2) -> integrale
integrale
Regola dei rettangoli su TI-82

5
1
2(1  e
 15 ( x  4)2
f ( x)  2(1  e
a 1
b5
n  100
)dx
 15 ( x  4)2
)
ClrHome
Input "A= ",A
Input "B= ",B
Input "N= ",N
Input "F(X)= ",Y1
(B-A)/ N -> H
Y1(A)-> S
For(K,2,N,1)
S+Y1(A+(K-1)*H) -> S
End
S*H -> S
Disp "Integrale= ",S
Regola dei trapezi su TI-82

5
1
2(1  e
 15 ( x  4)2
f ( x)  2(1  e
a 1
b5
n  100
)dx
 15 ( x  4)2
)
ClrHome
Input "A= ",A
Input "B= ",B
Input "N= ",N
Input "F(X)= ",Y1
(B-A)/N -> H
(Y1(A)+ Y1(B))/2 -> S
For(K,2,N,1)
S+Y1(A+(K-1)*H) -> S
End
S*H -> S
Disp "Integrale= ",S
Stima dell’errore
• Metodo dei rettangoli
| f '( x) | M 1
( b  a )2
n
| I  Rn | M 1
1
2
• Metodo dei trapezi
| f "( x) | M 2
| I  Tn | M 2
1
8
Ad esempio per il calcolo di


0
(b a )3
n2
sin( x)dx si può assumere M1 = M2 = 1
per cui con “sole” n = 250 suddivisioni si ha

 sin( x)dx - Rn |  0.0197392

• |  sin( x)dx - Tn |  0.00024805
•|
(NB: si divide per n )
0
0
(NB: si divide per n² )
Esercizio
• Calcolare

5
1
2(1  e
 15 ( x 4)2
)dx con il metodo dei trapezi
con n = 250 e fornire una stima dell’errore esaminando graficamente
la derivata seconda con R (o con la funzione TRACE della TI-82)
• Suggerimento:
f ( x)  2(1  e
f "( x) 
4
25
 15 ( x  4)2
)
(27  16 x  2 x )e
2
 15 ( x  4)2
I matematici hanno dimostrato che per questo integrale (di una
funzione importantissima: la gaussiana) non esistono metodi esatti.
| f "( x) | 
4
25
| ( 27  16 x  2 x ) | e
2
 15 ( x  4)2
> f <- function(x)
2*(1+exp(-(1/5)*(x-4)^2))
> a <- 1
> b <- 5
> n <- 250
> h <-(b-a)/n> x<-a+c(0:n)*h
> y <-f(x)
> plot(x,y)
> (sum(y)-(y[1]+y[n+1])/2)*h
[1] 13.60861
> f2 <- function(x)
(4/25)*(27-16*x+2*x^2)*exp(-(1/5)*(x-4)^2)
> curve(abs(f2(x)),1,5)
> M2 <- 0.8
> (1/8)*M2*(b-a)^3/n^2
[1] 0.0001024
>
2.9.4. Metodo Monte-Carlo
Calcolare l’area della parte di piano definita dalla disuguaglianza
2 x 2 ( x  1)4 
y2
2
1
 20
Fissiamo un rettangolo [a, b] × [c, d ] che contenga tutto il “pesce” e
spariamo n = 50000 (cinquantamila) punti a caso nel rettangolo.
Contiamo il numero k dei punti che colpiscono il pesce
(bordo del pesce compreso). Allora sarà:
area cercata = (k / n ) × area totale del rettangolo
Esempio
{{a, b}, {c, d}} = {{-0.15, 1.37}, {-0.35, 0.35}};
area totale = (b - a) (d - c) = 1.064;
n = 50000;
k = 34280;
Integrale (o area) = 0.729478
Formalmente sarebbe
(ma gli estremi sono comunque calcolati numericamente!)
Metodo Monte-Carlo su R
ycaso
(insuccesso)
ycaso
(successo)
xcaso
xcaso
Metodo Monte-Carlo su R
Il comando which
> x <- c(1,3,4,1,3,5,6,3)
> z <- which(x < 4)
> z
[1] 1 2 4 5 8
x[1]=1
x[2]=3
x[5]=3
x[8]=3
x[4]=1
z ha 5 elementi, che non sono i cinque elementi di x minori di 4,
bensì i cinque indici di tali elementi, ordinati come lo sono in x.
La lista c(1,3,1,3,3)degli elementi di x minori di 4 è x[z].
Metodo Monte-Carlo su R
>
>
>
>
>
>
>
>
>
>
>
>
>
>
f <- function(x) formula di f(x)
a <- valore di a
b <- valore di b
c <- valore di c
d <- valore di d
prove <- numero delle prove
xcaso <- runif(prove,min=a,max=b)
ycaso <- runif(prove,min=c,max=d)
z <- which(ycaso < f(xcaso))
successi <- length(z)
p <- successi/prove
integrale <- p*(b-a)*(d-c)
plot(xcaso[z],ycaso[z], col="red")
plot(f,0,b, add=TRUE, col="blue")
Metodo Monte-Carlo su TI-82, I
(in blu i comandi essenziali)
ClrHome
Esempio:
ClrDraw
PlotsOff
FnOff
Disp "----------------"
Disp "Integrazione"
Disp "di f(x)  0 "
Disp "Met. Monte Carlo"
Input "A= ",A
Input "B= ",B
Input "N= ",N
Input "f(x)= ", Y1

5
0
xe x dx
Metodo Monte-Carlo su TI-82, II
0 -> C
(B-A)/100 -> W
max(seq(Y1(X),X,A,B,W) -> D
A -> Xmin
B -> Xmax
C -> Ymin
D -> Ymax
AxesOff
DrawF Y1(X)
0 -> S
Text(55,1,"Successi =")
Text(47,1,"Prove =")
Metodo Monte-Carlo su TI-82, III
For(K,1,N)
Text(47,28,K)
A+rand*(B-A) -> X
C+rand*(D-C) -> Y
If Y  Y1(X)
Then
Pt-On(X,Y)
S+1 -> S
Text(55,40,S)
End
End
Text(5,60,"[ ENTER ]")
Pause
ClrHome
(B-A)*(D-C)*S/N -> I
Disp "Integrale =",I