1

I've been trying to compute the integral of (10*6)*sin(x0+x1+x2+x3+x4+x5+x6+x7)dx0dx1dx2dx3dx4dx5dx6dx7 using Monte-Carlo methods as part of my computing coursework at uni. I've managed to get the code to work after some trying, but it doesn't seem to give the correct answer. I've been given a reference analytical result of 537.187.... but my code converges to about 360

I'm not sure what else to try to resolve it. Is it something to do with my random numbers, could they be correlated by defining them seperately? Should I use a single 8 dimensional generator? Or do I need to use a more stratified sampler?

Apologies if it's something small, I'm quite new to Python!

UPDATE: I'm sorry I forgot to say what limits I was using. I'm integrating from 0 to pi/8 in each dimension

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=20000    #Number of random numbers generated in each dimension

x0rand, x1rand, x2rand, x3rand, x4rand, x5rand, x6rand, x7rand = [[0]*N]*8
for i in range(N):
   x0rand[i]= np.random.uniform(a,b) #Generates N random numbers 8D between 0 and pi/8
   x1rand[i]= np.random.uniform(a,b)
   x2rand[i]= np.random.uniform(a,b)
   x3rand[i]= np.random.uniform(a,b)
   x4rand[i]= np.random.uniform(a,b)
   x5rand[i]= np.random.uniform(a,b)
   x6rand[i]= np.random.uniform(a,b)
   x7rand[i]= np.random.uniform(a,b)

def func(x0,x1,x2,x3,x4,x5,x6,x7):  #Defining the integrand
return (10**6)*np.sin(x0+x1+x2+x3+x4+x5+x6+x7)

integral = 0.0
for i in range(N):
   integral += func(x0rand[i],x1rand[i],x2rand[i],x3rand[i],x4rand[i],x5rand[i],x6rand[i],x7rand[i])

answer=(((b-a)**8)/N)*integral
print ('The Integral is:', answer)
Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
Jake A
  • 25
  • 4
  • np.random.uniform takes a third argument which is the number of samples. You could just do: `x0rand = np.random.uniform(a,b,N)` instead of the for loop. It's more pythonic and orders of magnitude faster. – GuillemB Mar 19 '21 at 13:05
  • There are a few things wrong with the initial code but the most critical one and why you converge to the wrong value is that all your `xNrand` lists contain all the same values! Which I don't quite understand why. – GuillemB Mar 19 '21 at 13:31

1 Answers1

3

This codes converges to your analytical answer. I've made it more compact by generating all the samples in one go. Basically, you generate a matrix of random numbers and each row is one sample.

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=200000    #Number of random numbers generated in each dimension

samples = np.random.uniform(a,b,(N,8))

def func(x):  #Defining the integrand
    return (10**6)*np.sin(np.sum(x))

integral = 0.0
for i in range(N):
    integral += func(samples[i,:])

answer=(((b-a)**8)/N)*integral
print ('The Integral is:', answer)

EDIT: Faster more efficient version

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=2000000    #Number of random numbers generated in each dimension

samples = np.random.uniform(a,b,(N,8))

integrand_all_samples = (10**6)*np.sin(np.sum(samples, axis=1))
sum_all_samples = np.sum(integrand_all_samples)


answer=(((b-a)**8)/N)*sum_all_samples
print ('The Integral is:', answer)

EDIT 2: What is going wrong

The problematic code can be reduced to:

N=4

a,b = [[0]*N]*2
# [[0]*N]*2: [[0,0,0,0],[0,0,0,0]]
# a: [0,0,0,0]
# b: [0,0,0,0]
for i in range(N):
    a[i]= 1
    b[i] = 2
# a: [2,2,2,2]
# b: [2,2,2,2]

a and b are pointing to the same list in memory despite having the different names. You can check this by looking at the output of id(a) and id(b). You truly only have one list. You can find the details of why here: List of lists changes reflected across sublists unexpectedly

GuillemB
  • 540
  • 4
  • 13