Skip to content

Instantly share code, notes, and snippets.

Created December 14, 2018 06:32
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
What would you like to do?
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
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
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))
[ax.legend() for ax in axs]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment