How to find the intersecting volume of four spheres using Monte-Carlo method

Computer Science Level 2

Use the Monte-Carlo simulation and random-numbers to find the volume of intersection of the following four spheres. Given a = 2 − 3 , b = 2 − 2 3 3 , c = 8 3 + 2 a=2-\sqrt3, b=2- \frac23 \sqrt3, c=\sqrt{\frac83} +2

with centers ( 1 , a , 2 ) , ( 2 , 2 , 2 ) , ( 3 , a , 2 ) , ( 2 , b , c ) (1,a,2) , (2,2,2) ,(3,a,2), (2,b,c) and radii of 2.


The answer is 3.37.

Frank Petiprin by Frank Petiprin , 81, USA

This section requires Javascript.
You are seeing this because something didn't load right. We suggest you, (a) try refreshing the page, (b) enabling javascript if it is disabled on your browser and, finally, (c) loading the non-javascript version of this page . We're sorry about the hassle.

2 solutions

Frank Petiprin
Mar 1, 2021 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 
import numpy as np
#Spheres are referenced as 1,2,3,4
#The following code can be expanded for any number of Spheres
#where the Target Box contains the smallest Sphere.
def Get_torf_FourSpheres(Ssize,SphN):
    global X,Y,Z,SphC,SphR,torfS1,torfS2,torfS3,torfS4,TorF
    ######### Find The x,y,z Intervals that define the smallest Sphere. 
    xl,xr=SphC[SphN][0]-SphR[SphN],SphC[SphN][0]+SphR[SphN]
    yl,yr=SphC[SphN][1]-SphR[SphN],SphC[SphN][1]+SphR[SphN]
    zl,zr=SphC[SphN][2]-SphR[SphN],SphC[SphN][2]+SphR[SphN]
    #Found Ssize Random Numbers For X Y Z in itervals (xl,xr),(yl,yr),(zl,zr)#
    X,Y=np.random.uniform(xl,xr,Ssize),np.random.uniform(yl,yr,Ssize)
    Z=np.random.uniform(zl,zr,Ssize)
    ##Code determines if the Random (x,y,z) in the Target Box lies in the intersection
    ## of the four spheres.
    torfS1=(X-SphC[0][0])**2+(Y-SphC[0][1])**2+(Z-SphC[0][2])**2<=(SphR[0]**2)
    torfS2=(X-SphC[1][0])**2+(Y-SphC[1][1])**2+(Z-SphC[1][2])**2<=(SphR[1]**2)
    torfS3=(X-SphC[2][0])**2+(Y-SphC[2][1])**2+(Z-SphC[2][2])**2<=(SphR[2]**2)
    torfS4=(X-SphC[3][0])**2+(Y-SphC[3][1])**2+(Z-SphC[3][2])**2<=(SphR[3]**2)
    TorF=torfS1 & torfS2 & torfS3 & torfS4 
############### Main Progran#############   
global X,Y,Z,SphC,SphR,torfS1,torfS2,torfS3,torfS4,TorF
cntTrue,cntFalse=0,0
#Ssize Is the size of X,Y,Z Random Arrays 
a,b,c,Ssize=2-np.sqrt(3),2-(2/3)*np.sqrt(3),np.sqrt(8/3)+2,200000000
#Arrays Radius and Centers For The Four Spheres 1,2,3,4
SphR,SphC=[2,2,2,2],[[1,a,2],[2,2,2],[3,a,2],[2,b,c]]
SphNum=int(input("Enter Smallest Sphere In Radius: "))
#Target Box Contains The Smallest Sphere
TarBoxVol=(2*SphR[SphNum-1])**3
Get_torf_FourSpheres(Ssize,SphNum-1)
#Find the count of Trues To Determine Hits In The Intersecting Volumes
cntTrue,cntFalse=np.sum(TorF),len(TorF)-cntTrue
approx=((cntTrue/Ssize)*TarBoxVol)
print('Total Number Of Attempts',Ssize)
print('Approximate Intersecting Volumne For The Four Spheres',approx)
######### Sample Runs ###############Enter Smallest Sphere In Radius: 1

#Total Number Of Attempts 50000000
#Approximate Intersecting Volumne For The Four Spheres 3.37761536 

#Total Number Of Attempts 100000000
#Approximate Intersecting Volumne For The Four Spheres 3.37795328

#Total Number Of Attempts 100000000
#Approximate Intersecting Volumne For The Four Spheres 3.37871552

#Total Number Of Attempts 200000000
#Approximate Intersecting Volumne For The Four Spheres 3.37692544

1 Helpful 2 Interesting 2 Brilliant 0 Confused
Steven Chase
Mar 1, 2021

I am generally a fan of Monte Carlo, but I opted for a more conventional approach here. I did a numerical triple integral in spherical coordinates over the entire volume of the first sphere, dividing the volume into N 3 N^3 partitions. Each infinitesimal sub-volume is added to a running volume sum if it also simultaneously lies within spheres 2, 3, and 4. There is reasonably good convergence of the result as the spatial resolution increases. With N = 100 N = 100 (one million sub-volumes), I get V = 3.47 V = 3.47 . With N = 500 N = 500 (125 million sub-volumes), I get V = 3.37 V = 3.37 . With N = 1000 N = 1000 (1 billion sub-volumes), I get V = 3.39 V = 3.39 .

Code attached: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
 
import math

Num = 1000

R = 2.0

a = 2.0 - math.sqrt(3.0)
b = 2.0 - (2.0/3.0)*math.sqrt(3.0)
c = math.sqrt(8.0/3.0) + 2.0

x1 = 1.0
y1 = a
z1 = 2.0

x2 = 2.0
y2 = 2.0
z2 = 2.0

x3 = 3.0
y3 = a
z3 = 2.0

x4 = 2.0
y4 = b
z4 = c

###################################################

dr = R/Num
dtheta = 2.0*math.pi/Num
dphi = math.pi/Num

###################################################

V = 0.0

r = 0.0

while r <= R:

    theta = 0.0

    while theta <= 2.0*math.pi:

        phi = 0.0

        while phi <= math.pi:

            x = x1 + r*math.cos(theta)*math.sin(phi)
            y = y1 + r*math.sin(theta)*math.sin(phi)
            z = z1 + r*math.cos(phi)

            dV = (r**2.0)*math.sin(phi) * dr*dtheta*dphi

            D2sq = (x-x2)**2.0 + (y-y2)**2.0 + (z-z2)**2.0
            D3sq = (x-x3)**2.0 + (y-y3)**2.0 + (z-z3)**2.0
            D4sq = (x-x4)**2.0 + (y-y4)**2.0 + (z-z4)**2.0

            if (D2sq < R**2.0) and (D3sq < R**2.0) and (D4sq < R**2.0):
                V = V + dV

            phi = phi + dphi

        theta = theta + dtheta

    r = r + dr

###################################################

print Num
print V

#>>> 
#100
#3.47049534078
#>>> ================================ RESTART ================================
#>>> 
#200
#3.42377483278
#>>> ================================ RESTART ================================
#>>> 
#400
#3.41535186496
#>>> ================================ RESTART ================================
#>>> 
#500
#3.37129464636
#>>>
#>>> 
#1000
#3.38655401639
#>>>

0 Helpful 0 Interesting 0 Brilliant 0 Confused

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...