CODEWITHSUNDEEP
pythonpython-3.xmontecarlonumerical-integration
Vishnu
Vishnu

Reputation: 551

Monte-Carlo Integration with dependent parameters

I would like to integrate the function given below using the monte-carlo method:

f(x,y,z) = (x**2 * y * z) + (y**2 * z**2 * x) + (x**3 * y**3 *z)

with the limits:

0 <= x < 4,

0 <= y < 2*(x**2),

0 <= z < 1.

This is what I tried so far.

from skmonaco import mcquad

def g(t):
   """
   The integrand.
   """
   x = t[0]
   y = t[1]
   z = t[2]

   f = (x**2 * y * z) + (y**2 * z**2 * x) + (x**3 * y**3 *z)

   return f

mcquad(g,npoints=100000,xl=[0., 0., 0.],xu=[4., 10., 1.],nprocs=4)

I am getting the correct answer if I set the limits of y between two constants. For eg: 0 <= y < 10. Can someone please tell me how do I set the limits of y depending on x?

Upvotes: 1

Views: 485

Answers (1)

Siong Thye Goh
Siong Thye Goh

Reputation: 3586

Let h(x,y,z) be equal to f(x,y,z) when (x,y,z) is inside the domain and 0 otherwise.

Integrating h over a bigger region 

0 <= x < 4,
0 <= y < 32,
0 <= z < 1

gives us the same result.

That is we can run the following code:

from skmonaco import mcquad

def g(t):
   """
   The integrand.
   """
   x, y, z = t
   if y < 2 * (x**2):
       return (x**2 * y * z) + (y**2 * z**2 * x) + (x**3 * y**3 *z)
   return 0

print(mcquad(g,npoints=100000,xl=[0., 0., 0.],xu=[4., 32, 1.],nprocs=4))

Upvotes: 1

Related Questions