# Error in the value Montecarlo cru subrutine

**URL:** <https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825>\
**Category:** Help\
**Created:** [November 19, 2023, 9:28am UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825 "2023-11-19T09:28:38Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [November 19, 2023, 9:28am UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/1 "2023-11-19T09:28:38Z")

</div>

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(2d0_a, 2d0_b, 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((pi**2)-(x**2)))

END FUNCTION func1

DOUBLE PRECISION FUNCTION func2(x)  
IMPLICIT NONE  
DOUBLE PRECISION:: x  
func2=(x**2 \* sin(x)-x**3)_(cos(x))\*\*2_sin(x)

END FUNCTION func2

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

---

<div class="post-metadata">

**Author:** ![Euler-37](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@Euler-37](https://fortran-lang.discourse.group/u/Euler-37)\
**Post date:** [November 19, 2023, 10:39am UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/2 "2023-11-19T10:39:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [November 19, 2023, 10:48am UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/3 "2023-11-19T10:48:46Z")

</div>

Thank you!!!

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [November 19, 2023, 12:15pm UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/4 "2023-11-19T12:15:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [November 19, 2023, 3:42pm UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/5 "2023-11-19T15:42:55Z")

</div>

> [@4ulalia](#):
>
> `t = rand()`

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.

> [@4ulalia](#):
>
> `func1=sqrt(abs((pi**2)-(x**2)))`

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.

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [November 19, 2023, 8:47pm UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/6 "2023-11-19T20:47:03Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/c/c47c4908ab8ac50d45aa2267909a7ce1210c6e31.png)

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!!

```auto
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

```

---

<div class="post-metadata">

**Author:** ![vmagnin](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/vmagnin/32/28_2.png) [@vmagnin](https://fortran-lang.discourse.group/u/vmagnin)\
**Post date:** [November 19, 2023, 8:53pm UTC](https://fortran-lang.discourse.group/t/error-in-the-value-montecarlo-cru-subrutine/6825/7 "2023-11-19T20:53:28Z")

</div>

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

````fortran`

`````