Skip to content

Instantly share code, notes, and snippets.

@dgerosa dgerosa/nsphere.py
Created Dec 14, 2018

Embed
What would you like to do?
nsphere
'''
Calculate the volume of the n-dimensional spehere a function of the dimension using a simple hit-or-miss Monte Carlo.
'''
from __future__ import print_function, division
import numpy as np
import scipy.special
import pylab as plt
import time
Dmax = 30 # Largest dimension
N = int(1e7) # Monte Carlo samples
data = []
for D in np.arange(2,Dmax+1): # Loop
t0=time.time()
all = np.random.uniform(-1,1, D*N).reshape(N,D) # Generate points in cube
distance = np.sum(all**2,axis=1)**0.5 # Distance from cube center
Ninside = np.sum(distance<=1) # Number of points within sphere
fractioninside = Ninside/N # Fraction of points within sphere
cube = 2**D # Volume of the cube (from -1 to 1)
sphere = fractioninside * cube # Volume of the sphere
solution = np.pi**(D/2) / scipy.special.gamma( D/2 + 1) # True value
diff= np.abs(sphere-solution)/solution # Error
t=time.time()-t0
print("D=%i\tN=%i\tN_in=%06d\tV=%.5f\t\tSol=%.5f\t\tdiff=%.5f\tt=%.2fs" %(D,N,Ninside,sphere,solution,diff,t))
data.append([D,sphere,solution,diff]) # Store data
# Plots...
dims, my,true,errors = np.array(data).T
fig, axs= plt.subplots(2, 1,figsize=(5,10))
axs[0].plot(dims,my,label='my',ls='-',marker='x')
axs[0].plot(dims,true,label='true',ls='-')
axs[1].plot(dims,errors,label='errors')
[ax.legend() for ax in axs]
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.