Error in the value Montecarlo cru subrutine

Help
1

Hello! I amb taking a begginer’s class into Fortran. One of the exercises they asked was solving an integral using the “Montecarlo-cru” method. I believe the subrutine I wrote is correct and so is the rest of the code, and i haven’t got any errors when executing, but the value it gives me isn’t the correct one. Could there be an error when calling the function? I have tried with some of my collegues subrutines and it still give mes the wrong value.

PROGRAM PREPA6
IMPLICIT NONE
DOUBLE PRECISION:: pi, a, b, integral1, integral2, error1, error2, func1, func2, x
INTEGER:: N, iseed
EXTERNAL:: func1, func2

pi = 4d0*atan(1d0)
iseed=20506813
call srand(iseed)

OPEN(10, file=‘prepa6-res.dat’)

a = -pi
b=pi
CALL montecarlo(a, b, 2000, func1, integral1, error1)
CALL montecarlo(2d0a, 2d0b, 2000, func2, integral2, error2)
print*, integral1, integral2

CLOSE(10)
END PROGRAM PREPA6

SUBROUTINE montecarlo(a, b, N, func, integral, error)
IMPLICIT NONE
DOUBLE PRECISION:: a, b, func, t, integral, error, suma, sumaquad, h
INTEGER:: i, N, iseed
iseed=20506813
call srand(iseed)

!h=(b-a)*func((b-a)*t+a)

!Inicialitzem
suma = 0d0
sumaquad = 0d0


do i = 1, N
	t = rand()
	suma = suma + (b-a)*func(((b-a)*t)+a)
	sumaquad = sumaquad + ((b-a)*func(((b-a)*t)+a))**2
enddo

integral = suma / dble(N)
error = (1/sqrt(dble(N)))*sqrt(abs(sumaquad/dble(N) -(suma/dble(N))**2))

END SUBROUTINE montecarlo

DOUBLE PRECISION FUNCTION func1(x)
IMPLICIT NONE
DOUBLE PRECISION:: pi, x
func1=sqrt(abs((pi2)-(x2)))

END FUNCTION func1

DOUBLE PRECISION FUNCTION func2(x)
IMPLICIT NONE
DOUBLE PRECISION:: x
func2=(x2 * sin(x)-x3)(cos(x))**2sin(x)

END FUNCTION func2

I program in a MAC OS. Please help would really be appreciated!!!

1 Like
2

In function func1, you should set pi=4*atan(1.d0), and for each monte carlo method, step = 2000 is too small

1 Like
3

Thank you!!!

4

For code readability I would also suggest that you specify input and output arguments using the intent in and intent out qualifications. Also, better to put all subroutines into a module that you will USE in the main program. Something like this:

module mymod
implicit none

contains

subroutine mysub(output1, input1)
implicit none ! Redundant

real(8), intent(in) :: input1
real(8), intent(out) :: output1

output1 = 2d0*input1

end subroutine mysub

end module mymod

program main
use mymod, only: mysub
implicit none

real(8) :: x, y

x = 1d0

call mysub(y,x)

write(* , * ) "x = ", x
write(* , * ) "y = ", y

end program main

1 Like
5

Maybe I overlooked it, but i do not see rand() declared anywhere. I have no idea what your compiler is doing with this function reference. Is it returning a REAL and converting to DOUBLE PRECISION, or is it returning a DOUBLE PRECISION?

Fortran has its own intrinsic pseudorandom number generator, random_number(), that works for all real kinds. You might consider switching to that to avoid any uncertainty or ambiguity in how the nonstandard rand() function is handled. In actual application codes, programmers often use their own PRNG codes, just so they control the algorithm and know its characteristics; for your coursework, you probably don’t need to go that far.

In this function the variable pi does not have a value. This is a local variable, so it does not assume the value of the identically named variable in the main program.

1 Like
6

Hello! I’m taking a begginer’s course on Fortran and I have been asked to program using the Montecarlo’s method. Specifically I have been asked to use Montercarlo with importance sampling.
I had to generate random numbers with distribuition p(x)= (10/3)* e^-x * sin(x)**2/(1+e^-pi) and a gaussian distribution with mu=0 and sigma=1/sqrt(2)).
I have generated the numbers using an acceptance/discarrd method for the first one and a box-muller method. Using these numbers I have been asked to calculate some integrals
Captura de Pantalla 2023-11-19 a les 21.42.27

Even though when executing it takes a bit of time, both first integrals give me the correct value (but I don’t know why it takes 2 minutes to execute). The third one I don’t know if i’m defining correctly the infiniy index aswell as I believe the gaussian function is incorrect 'cause it doesn’t give me any values, only NaNs or 0.00.
I don’t exactly know how to fix this, help would be very much appreciated. Thank you!!

PROGRAM PREPA6
IMPLICIT NONE
DOUBLE PRECISION:: pi, a, b, integral1,  integral2, error1, error2, func1, func2, xlow, xhigh, cotasup, distribution, sigma, mu
DOUBLE PRECISION:: func3, func4, func5, integral3, error3, integral4, integral5, error4, error5, gaussiana
DOUBLE PRECISION, DIMENSION(1100000):: nums, nums_gauss
INTEGER:: N, iseed, ndades
EXTERNAL:: func1, func2, distribution, func3, func4, func5, gaussiana

pi = 4d0*atan(1d0)
iseed=20506813
call srand(iseed)

OPEN(10, file='prepa6-res.dat')
!.... code

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!APARTAT 1b
ndades=1100000
xlow=0d0
xhigh=pi
cotasup=0.8d0
CALL accepta(ndades,nums,xlow,xhigh,cotasup,distribution)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!APARTAT 1c
ndades=1100000
mu=0
sigma=1d0/sqrt(2d0)
CALL BOX_MULLER (ndades, nums_gauss, mu, sigma)
!print*, gaussiana(100d0, mu, sigma)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!APARTAT 1d
WRITE(10, *) "#INTEGRALS I ERRORS AMB LLISTA DE NUMS. ALEATÒRIS"
do N=5000, 1100000, 5000
mu=0
sigma=1d0/sqrt(2d0)
CALL montecarlo_sam(0d0, pi, N, func3, distribution, integral3, error3, nums)
CALL montecarlo_sam(0d0, pi, N, func4, distribution, integral4, error4, nums)
CALL montecarlo_sam(-1.0D12, 1.0D12, N, func5, gaussiana, integral5, error5, nums_gauss)
WRITE(10, '(I8, 4(2X, F20.12, E20.12))') N, integral3, error3, integral4, error4, integral5, error5
ENDDO

CLOSE(10)
END PROGRAM PREPA6


!Subrutina que genra una distribució gausianna aleatòria utilitzant el mèode de BOX-MULLER
SUBROUTINE BOX_MULLER (ndades, numeros, mu, sigma)
    IMPLICIT NONE
    DOUBLE PRECISION:: y1, y2, mu, sigma, pi, numeros(ndades), r, phi, x1, x2
    INTEGER:: ndades, i

    pi = 4.d0 * atan(1.d0)

     do i = 1, ndades, 2

        r = sqrt(-2.d0*log(rand()))
        phi = 2.d0 * pi * rand()

        x1 = r*cos(phi)
        x2 = r*sin(phi)

        y1 = mu + x1*sigma
        y2 = mu + x2*sigma

        numeros(i) = y1
        numeros(i+1) = y2

    enddo

END SUBROUTINE

SUBROUTINE  accepta(ndades,numeros,xlow,xhigh,cotasup,funcio)
    IMPLICIT NONE
    !Declarem les variables
    DOUBLE PRECISION:: xhigh, xlow, cotasup, funcio, x, p, numeros(ndades),dens
    INTEGER:: i, ndades, iseed
    
       

    iseed=20506813
    call srand(iseed)
       i=1
do while (i.LT.ndades)
    !fer x i p. x[xlow,xhigh] i p[0,cotasup] 
!a partir d'una distribució uniforme x[0,1] passam a una y[a,b] segons y=(b-a)x+a 
    x = (xhigh - xlow)* rand() + xlow
    p = cotasup * rand()
    dens=funcio(x)

    if (dens.GE.p) then
        numeros(i) = x
        i=i+1
    endif

enddo
return

END SUBROUTINE accepta


SUBROUTINE montecarlo_sam(a, b, N, func, densitat, integral, error, numeros)
	IMPLICIT NONE
	DOUBLE PRECISION:: a, b, func, t, integral, error, suma, sumaquad, numeros(N), densitat
	INTEGER:: i, N, iseed
	iseed=20506813
	call srand(iseed)

	!h=(b-a)*func((b-a)*t+a)
    
    !Inicialitzem
suma = 0d0
sumaquad = 0d0

	
do i = 1, N
	t = numeros(i)
	suma = suma + func(t)/densitat(t)
	sumaquad = sumaquad + (func(t)/densitat(t))**2
enddo

integral = suma / dble(N)
error = (1/sqrt(dble(N)))*sqrt(abs(sumaquad/N -(suma/N)**2))
return

END SUBROUTINE montecarlo_sam


!SUBROUTINE montecarlo_multidim(a, b, N, func, integral, error, numeros)
!DOUBLE PRECISION:: x1, x2, x3, x4, x5

!!Funcions densitat de probabilitat
DOUBLE PRECISION FUNCTION distribution(x)
DOUBLE PRECISION::x, pi
pi=4d0*atan(1d0)
distribution= ((10d0/3d0)*exp(-x)*(sin(x)**3))/(1d0+exp(-pi))
END FUNCTION distribution

DOUBLE PRECISION FUNCTION gaussiana(x, mu, sigma)
DOUBLE PRECISION:: mu, sigma, x, pi
pi = 4d0*atan(1d0)
gaussiana = exp(-(x-mu)**2/(2d0*sigma**2))/(sigma*sqrt(2d0*pi))
END FUNCTION gaussiana


!Funció corresponent a l'integral 3 [0, pi]
DOUBLE PRECISION FUNCTION func3(x)
DOUBLE PRECISION:: x
func3=exp(-abs(x))*(x**2)*(sin(x)**2)
END FUNCTION func3

!Funció corresponent a l'integral 4 [0, pi]
DOUBLE PRECISION FUNCTION func4(x)
DOUBLE PRECISION:: x, pi
pi=4d0*atan(1d0)
func4=exp(-(x**2)/2)*(pi+(4d0*(x**2)))*(cos(x)**2)
END FUNCTION func4

!Funció corresponent a l'integral 5 [-infinit, infinit]
DOUBLE PRECISION FUNCTION func5(x)
DOUBLE PRECISION:: x
func5=exp(-(x**2))*(x**2)*(sin(x)**2)
END FUNCTION func5
2 Likes
7

You can put your Fortran codes between two pairs of triple quotes like this:

```fortran

```

2 Likes