Search code examples
pythonjupyter-notebookmontecarlo

How to use Monte Carlo method to find the volume of a sphere in python?


I need to use the Monte Carlo method to find the volume of a sphere in a python3 environment in jupyter notebook. But the examples of the method in my lecturer's code doesn't seem to make sense.

I've been instructed to modify my lecturer's code which uses the monte carlo method to find pi using a circle in a unit square; however when I run their code I get the value of pi to be 0.04 which is far too small. I know how to adapt it to find the volume of a sphere, but again the answer I get is far too small. Below is the lecturer's code.

import numpy as np

n = 100
x = np.random.random(n)
y = np.random.random(n)

z = np.sqrt(x**2 + y**2)

ninside = len( np.where(z < 1.) )
pi = 1.*ninside/n * 4.

print(pi)

When you run the code and print pi, you get 0.04 which is obviously not correct, it should be roughly 3.14. How is this code supposed to work and how am I supposed to use the Monte Carlo method to find the volume of a sphere?


Solution

  • Replace

    ninside = len( np.where(z < 1.) )
    

    with

    ninside = len( np.where(z < 1.)[0] )
    

    np.where() returns a tuple, so you need to take the array out of the tuple before getting the length.

    When debugging something like this, I recommend running it line-by-line in an interpreter and checking the value of every new variable.

    That's how I found that len(np.where(z < 1.)) was always 1, and then I looked at the value of np.where(z < 1.), and found it was a tuple with an array inside it.