Computational Physics - Department of Physics

(Axel Boer) #1

374 11 Outline of the Monte Carlo Strategy


}
fx=gaussian_MC(x);
int_mc += fx;
sum_sigma += fxfx;
}
int_mc = int_mc/((double) n );
sum_sigma = sum_sigma/((double) n );
variance=sum_sigma-int_mc
int_mc;
// final output
cout << setiosflags(ios::showpoint | ios::uppercase);
cout <<" Monte carlo result= "<< setw(10) << setprecision(8) << jacobidetint_mc;
cout <<" Sigma= "<< setw(10) << setprecision(8) << volume
sqrt(variance/((double) n
)) << endl;
return0;
}// end of main program


// this function defines the integrand to integrate


doublegaussian_MC(double*x)
{
// evaluate the different terms of the exponential
doublexy=pow((x[0]-x[3]),2)+pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2);
returnxy;
}// end function for the integrand


// random numbers with gaussian distribution
doublegaussian_deviate(longidum)
{
static intiset = 0;
static doublegset;
doublefac, rsq, v1, v2;
if( idum < 0) iset =0;
if(iset == 0){
do {
v1 = 2.
ran0(idum) -1.0;
v2 = 2.ran0(idum) -1.0;
rsq = v1
v1+v2v2;
}while(rsq >= 1.0 || rsq == 0.);
fac = sqrt(-2.
log(rsq)/rsq);
gset = v1fac;
iset = 1;
returnv2
fac;
}else{
iset =0;
returngset;
}
}// end function for gaussian deviates


The following table lists the results from the above two programs as function of the number of
Monte Carlo samples. The suffixcrstands for the brute force approach whilegdstands for the
use of a Gaussian distribution function. One sees clearly that the approachwith a Gaussian
distribution function yields a much improved numerical result, with fewer samples.

Free download pdf