# Use Monte Carlo integration to compute the volume and centre of mass in Python

In particular, I want to focus on finding the volume $$V$$ because I will need it to start working on solving the centre of mass

This $$3D$$ homogenous body (Torus section) is defined by

$$x^2 + \left(\sqrt{y^2+z^2}-5\right)^2 < 4$$ where,

$$y<1, z>2$$

I confined the $$3D$$ shape to a rectangular region where

$$-2

The volume is defined as follows $$V = \iiint dxdydz$$

Below is my attempt at coding the integral to find the volume.

I am new to Monte Carlo coding so any feedback would be appreciated.

import numpy as np
from scipy import random

#Limits of Integration
x_min = -2
x_max = 2

y_min = -7
y_max = 1

z_min = 2
z_max = 7

#Number of Iterations
N = 1000

Xrand = random.uniform(x_min,x_max,N)
Yrand = random.uniform(y_min,y_max,N)
Zrand = random.uniform(z_min,z_max,N)

for i in range(len(Xrand)):
Xrand[i] = random.uniform(x_min,x_max)
for i in range(len(Yrand)):
Yrand[i] = random.uniform(y_min,y_max)
for i in range(len(Zrand)):
Zrand[i] = random.uniform(z_min,z_max)

def func(x,y,z):
return x**2 + (np.sqrt(y**2 + z**2)-5)**2

integral = 0.0

for i in range(N):
integral += func(Xrand[i],Yrand[i],Zrand[i])

Volume =(1/N)*float(N)*integral
print(Volume)

• To calculate the volume with MC, you should test if a random point is inside the volume, if it is then add it to the sum. Oct 14 '20 at 18:00
• Can you add a more explicit question? What is the problem? Oct 15 '20 at 14:18
• @MPIchael The problem is just to use the Monte Carlo integration to compute the volume and centre of mass of the 3-dimensional homogeneous body (torus section) defined by the equation given above Oct 15 '20 at 15:02

Here is a corrected and slightly improved code, set up here for calculation of the full torus volume to verify the result.

import numpy as np
from scipy import random

def func(x,y,z):
return x**2 + (np.sqrt(y**2 + z**2)-5)**2

def MCvolume(N=1000):
#Limits of Integration
x_min = -2
x_max = 2

y_min = -7
y_max = 7 #1

z_min = -7 #2
z_max = 7

Xrand = random.uniform(x_min,x_max,N)
Yrand = random.uniform(y_min,y_max,N)
Zrand = random.uniform(z_min,z_max,N)

for i in range(len(Xrand)):
Xrand[i] = random.uniform(x_min,x_max)
for i in range(len(Yrand)):
Yrand[i] = random.uniform(y_min,y_max)
for i in range(len(Zrand)):
Zrand[i] = random.uniform(z_min,z_max)

integral = 0.0

for i in range(N):
if (func(Xrand[i],Yrand[i],Zrand[i]) < 4.0):
integral += 1.0

VolumeBox=(x_max-x_min)*(y_max-y_min)*(z_max-z_min)
Volume = VolumeBox*integral/float(N)

return Volume

#-stats from M trials, in each one using N random points
Mtry=10
Nrand=1000
vols=np.zeros(Mtry)

for k in range(0,Mtry):
vols[k]=MCvolume(N=Nrand)

print("For Mtry, Nrand =", Mtry, ", ", Nrand, "; mean=", np.mean(vols), ", standard deviation=", np.std(vols))
print("Exact torus volume = ", 2*np.pi*5.0*np.pi*2.0**2)


You are missing an integral part of what makes up a Monte Carlo integration.

The steps to integrate the area (in your case volume) are, to generate random samples, across the whole domain which includes the wanted volume, (which you did). Now, you want to sum up all the dots that are within the region you want to integrate , and discard those that do not lie in it. As I see it, you currently do not check whether your random samples are within the torus! Two check whether they lie within it, you just test if the torus definition you stated is true for the sample or not.

Secondly, you have to think about the weights you are going to assign to those hits that make up your integral. You are distributing a total of 1000 random points through your domain, which you define by x_max,x_min etc. The weights should then be: 1/Samples * Volume_sampled. When you fix these two things in your code, you should get to the correct volume integral.

Also, I suggest you increase the sample size dramatically, MC methods converge slowly.