Multidimensional Monte Carlo Integral

[4ulalia]

Hello, I’m a beginner Fortran student and i have been asked to compute this integral using numbers I generated with the box muller method. I don’t understand how the 5 dimensional integral works since the numbers I generated are only 1 dimensional, how can i have 5 dimensions?

I generated a code but the value comes to infinity so I believe it’s incorrect. Can someone explain and fix the code?

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, dis2, func6, &
integral6, error6
DOUBLE PRECISION, DIMENSION(1100000):: nums, nums_gauss
INTEGER:: N, iseed, ndades
EXTERNAL:: func1, func2, distribution, func3, func4, func5, gaussiana, dis2, func6

!... other code

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

CALL montecarlo_5dim(-10d0, 10d0, 210000, func6, dis2, integral6, error6, nums_gauss)
print*, integral6
!DO N=1500, 210000 ,1500

CLOSE(10)
END PROGRAM PREPA6
UBROUTINE 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
DOUBLE PRECISION FUNCTION dis2(x1, x2, x3, x4, x5)
DOUBLE PRECISION:: x1, x2, x3, x4, x5, pi
pi = 4d0*atan(1d0)
dis2=exp(x2*cos(x2+x3-x5))*((pi*(x3*x4*x5)**2)+((cos(x3+x4-(2d0*x1)))**2)*x3*sin(x5))
END FUNCTION dis2



SUBROUTINE montecarlo_5dim(a, b, N, func, densitat, integral, error, numeros)
  IMPLICIT NONE
  DOUBLE PRECISION, DIMENSION(N, 5) :: numeros
  DOUBLE PRECISION:: x1, x2, x3, x4, x5
  DOUBLE PRECISION :: func, integral, error, suma, sumaquad,f, a , b, g, densitat
  INTEGER :: N, i

  ! Initialize variables
  suma = 0.0d0
  sumaquad = 0.0d0

  ! Perform Monte Carlo integration
  do i = 1, N
    !Genera les llistes de valors x_i a partir dels nums aleatoris generats
    x1= numeros(i, 1)
    x2= numeros(i, 2)
    x3= numeros(i, 3)
    x4= numeros(i, 4)
    x5= numeros(i, 5)

    f=func(x1,x2,x3,x4,x5)
    g=densitat(x1,x2,x3,x4,x5)

    suma = suma+ f/g
    sumaquad= sumaquad +(f/g)**2
    
  end do

  ! Calculate integral and error
  integral = suma / dble(N)
  error = sqrt((sumaquad / dble(N) - integral**2) / dble(N))

END SUBROUTINE montecarlo_5dim
!Funció per la integral multidimensional
DOUBLE PRECISION FUNCTION func6(x1,x2,x3,x4,x5)
DOUBLE PRECISION:: x1, x2, x3, x4, x5
func6=exp(-(x1**2)+(x2**2)+(2*(x3**2))+(x4**2)+(2*(x5**2)))
END FUNCTION func6```

[Arjen]

What is the outcome if you use a simple function g(x1,x2,x3,x4,x5), such as g(…) = 1? You can evaluate that analytically. And therefore you get a chance to see where things go wrong.

(By the looks of it, the original integral should indeed be finite )

[conradoat]

Is this a homework?

[Arjen]

Yes, it is - see previous messages.

[conradoat]

What is the outcome if you use a simple function g(x1,x2,x3,x4,x5), such as g(…) = 1? You can evaluate that analytically. And therefore you get a chance to see where things go wrong.

This is an excellent example of how to start debugging the code of your homework yourself, @4ulalia.

[feenberg]

The densitat function is called 5 times for each evaluation of g. That gives you your 5 dimensions. The suggestion to use g()=1 is problematic - over the interval +/- infinity that will be infinite. I suggest setting g()=1 only for the range +/-1 while you experiment with the code.

I sense you are uncomfortable with Monte-Carlo integration. It might help to start with an integral over a finite range. In that case you can see that the integral is equal to the average value of the function times the range. The range is known, the average value of the function is approximated by sample mean. Another way to calculate the same thing is to sum the product of the function and the probability of those x values generating that result over a sample of x values. This extends to the case of infinite ranges. All you need is a source of random numbers over that range, with the probability for each.

[Arjen]

Well, the weight is Gaussian, so there should be a quick decay, but the infinite range is not easy to deal with. Indeed, limit the range to [-1,1] and proceed from there.