Questo programma illustra il funzionamento dei principali metodi numerici di integrazione : metodo dei rettangoli, dei trapezi, simpson e montecarlo. La funzione “funzione” viene integrata nell’intervallo [A,B] (modificabile tramite #define) con precisione scelta (sempre modificabile nel #define).
/********************
*
* Il programma restituisce l’integrale di una funzione in un
* intervallo scelto calcolato con diversi metodi di integrazione
*
* Autore : snoomerz
* (Web: http://snoomerz.netsons.org)
*
* Rilasciato con licenza Creative Commons Attribuzione 2.5 Italia
* (http://creativecommons.org/licenses/by/2.5/it/)
*
*********************/
#include <math.h>
#define A 1 /* Estremi dell’intervallo */
#define B 2 /* di integrazione */
#define N_INT 10000 /* Numero di suddivisioni intervallo per metodi dei rettangoli, dei trapezi, simpson */
#define N_MNT 100000 /* Numero suddivisioni intervallo per metodo Montecarlo */
#define SEED 987654321
double funzione(double x);
double rettangoli(double(*fp)(double), double a, double b, int n_int);
double trapezi(double(*fp)(double), double a, double b, int n_int);
double simpson(double(*fp)(double), double a, double b, int n_int);
double montecarlo(double(*fp)(double), double a, double b, int n_mnt);
int main(void){
printf(”\nL’ integrale della funzione vale :\n\n”);
printf(”Metodo dei rettangoli: %f\n”, rettangoli(funzione, A, B, N_INT));
printf(”Metodo dei trapezi : %f\n”, trapezi(funzione, A, B, N_INT));
printf(”Metodo simpson : %f\n”, simpson(funzione, A, B, N_INT));
printf(”Metodo montecarlo : %f\n\n”, montecarlo(funzione, A, B, N_MNT));
return 0;
}
double funzione(double x){
return exp(x)+log(x); /* Funzione da integrare */
}
double rettangoli(double(*fp)(double), double a, double b, int n_int){
double h=0, value=0, i=0;
h=(b-a)/n_int;
for(i=a;i<b;i+=h) value+=h*(*fp)(i+h/2);
return value;
}
double trapezi(double(*fp)(double), double a, double b, int n_int){
double h=0, value=0, i=0;
h=(b-a)/n_int;
for(i=a;i<b;i+=h) value+=(((*fp)(i)+(*fp)(i+h))*h/2);
return value ;
}
double montecarlo(double(*fp)(double), double a, double b, int n_mnt){
double value=0, i=0;
srand(SEED);
for(i=0;i<n_mnt;i++) value+=((b-a)/n_mnt)*(*fp)((b-a)*(double)rand()/RAND_MAX + a); /* (b-a)*(double)rand()/RAND_MAX + a restituisce un valore casuale compreso fra a e b*/
return value;
}
double simpson(double(*fp)(double), double a, double b, int n_int){
double h=0, i=0, par=0, dis=0;
int j=0;
if ( n_int % 2 ) ++n_int;
h=(b-a)/n_int;
for(i=a+h;i<b;i+=h){
if(j%2==0) par+=(*fp)(i);
else dis+=(*fp)(i);
j++;
}
return h* ((*fp)(a)/3 + 4*par/3 + 2*dis/3 + (*fp)(b)/3);
}
! Anche se a dire la verità mi manca un po’ sudo …
