Integrazione numerica con il metodo Monte Carlo

Il metodo Monte Carlo può anche essere utilizzato per effettuare il calcolo approssimato dell'integrale definito. Supponiamo di individuare un rettangolo che contiene la funzione da integrare nell'intervallo acelto facendo attenzione a selezionare con cura i massimi e i minimi della funzione nell'intervallo.

Adesso supponiamo di generare dentro il quadrato tanti punti casuali; Contiamo quanti punti abbiamo generato nel rettangolo e quanti punti sono caduti sotto la curva all'interno del segmento (b-a); per motivi di probabilità non c'è nessuna ragione che la frequenza dei punti generati sotto la curva abbiano una frequenza diversa dai punti generati dentro il quadrato; quindi in formula avremo la seguente proporzione;

Area sotto la curva : Area del rettangolo = (Integrale tra i punti a e b ) : (b-a)*f(b)

dove (b-a) è la base del rettangolo ed f(b) è l'altezza del rettangolo che coincide con il valore di f(x) calcolato per x=b.

Per calcolare il numero dei punti appartenenti all'area sottesa dalla curva è sufficiente confrontare il valore di Yn con il valore della funzione f(Xn)  e se Yn<= f(Xn) allora il punto appartiene all'area sottesa dalla curva.

Esempio

Esercizio

Adesso scriviamo un programma in dev C che attraverso la funzione rand() generi dei numeri casuali e conti quanti punti sono generati sotto la curva rispetto a quelli generati nel quadrato. In questo modo risalire al valore dell'integrale tra a e b.

(NOTA: SE vuoi il programma comleto scrivi nel box sottostante!!!)

Il programma in dev C è il seguente

#include<...........>   //includi le librerie necessarie

creiamo la funzione

float f(float x){

float y;

y=3*x*x+2*x+1;              //funzione di cui vogliamo calcolare l'integrale

return y;

}

scriviamo il main

main(){

.......................                                            

-------------------                                        

                                                                      //calcoliamo l'integrale con gli estremi fissati dall'esercizio

I=4*86*double(dentro)/double(tanti);             //questo perché  nel nostro esempio

                                                                          //(b-a)=4  e MAX=f(5)=86

                                                                          //nota: se cambiamo gli estremi di integrazione

                                                                          //cambierà sia (b-a) che il valore MAX

//adesso abbiamo calcolato l'integrale e dobbiamo solo stamparlo

printf("\nIl valore dell'integrale calcolato con il metodo Monte Carlo\n");

printf("              I=%8.16f\n",I);

printf("  Numeri dentro=%d\n",dentro);

printf("Numeri generati=%d\n",tanti);

getch();

}

 In output avremo:

I=151,9897308349609400

Con 10000000 numeri generati che si avvicina all’integrale cercato.

 

Scrivere con il dev C un programma che calcoli l' integrale con il metodo Monte Carlo

ESERCIZIO 1:

Calcolare l'integrale  della funzione y=3x*x+2*x+1 con il metodo Monte Carlo inserendo gli estremi a e b dall'esterno.

ESERCIZIO 2:

Calcolare l'integrale della funzione y=x*x*x+3 con gli estremi inseriti dall'esterno.

ESERCIZIO 3:
Calcolare l'integrale della funzione  y=x*x-sin(x) con il metodo monteCarlo con a=2 e b=10.

Svolgimento:

Per svolgere i precedenti esercizi, consideriamo quindi in generale una funzione y=f(x) e il relativo integrale.

In questo caso, con il nostro generatore di numeri pseudo casuali uniformemente distribuito, dovremo generare le coordinate dei punti Pi che hanno coordinate xi, e yi cioè Pi =[xi;yi] tali che:

- l’ascissa xi sia  compresa nell’intervallo a e b


- l’ordinata yi sia  compresa tra f(a) e f(b) (nel nostro caso tra 0 e f(b), se consideriamo l'integrale sotteso alla curva nel primo quadrante, dove è necessario valutare, 

- caso per caso, la dimensione verticale del rettangolo, cioè il suo valore massimo.

Tutti questi punti dovranno essere divisi in due gruppi:

- quelli appartenenti a tutto il rettangolo
- quelli appartenenti all’area sottesa dalla curva.

Essendo il generatore equiprobabile, è ragionevole ipotizzare che il numero di punti presenti nelle due aree risulti proporzionale al valore delle aree stesse e cioè:

(Numero dei punti sottesi)/(numeri dei punti totali) = (area sottesa alla curva)/(area del rettangolo)   = (integrale di f(x))/(B-A)*f(B)

Per confrontare il numero dei punti appartenenti all’area sottesa dalla curva è sufficiente confrontare il valore della funzione nel punto xi con il valore yi e se yi <= f(xi) allora il punto appartiene all’area sottesa dalla curva.