Tuesday, November 16, 2010

Fortran program to calculate an integral using "crude" Monte Carlo

The following program calculates the integral I

using "crude" Monte Carlo. In this case, n=3, xia=0, xib=1 and f(x)=exp(a·x).


Program intmc

double precision fsum,fsum2,a,xa,xb,V,x,f,I,IA,varf
integer M,n,j,k
dimension a(3),xa(3),xb(3),x(3)
c n: dimension of the array x
c M: number of random numbers we use
c V: volum of the hipervolume
n=3
M=1000000
V=1
do j=1,n
a(j)=rand()
xa(j)=0
xb(j)=1
V=V*(xb(j)-xa(j))
enddo
c let's find the value of the integral that we get using monte-carlo
fsum=0
c fsum2 will help us to find the error
fsum2=0
do k=1,M
do j=1,n
x(j)=rand()
enddo
call func(a,x,f,n)
fsum=fsum+f
fsum2=fsum2+f**2
enddo
c I is the value of the integral we wanted to find
c IA is the value of the integral calculated analiticaly
I=V*fsum/M
varf=(fsum2/M-((fsum/M)**2))
sigma=sqrt(varf/M)
IA=1
do j=1,n
IA=IA*((exp(a(j))-1)/a(j))
enddo
write(*,*) I,sigma*3,IA
stop
end


SUBROUTINE func(a,x,f,n)
double precision a,x,f,ax
integer n
dimension a(n),x(n)
ax=0
do j=1,n
ax=ax+a(j)*x(j)
enddo
f=exp(ax)
end

1 comment: