Volume d'une hyperboule

###################################################################
# Calcul du volume d'une hypersphère en dimension n
# Rescassol - Mai 2019
###################################################################

###################################################################
# Importations
###################################################################

import numpy as np
from scipy.integrate import nquad

###################################################################
# Fonctions
###################################################################

# MC comme Monte-Carlo

def VolMC1(n): # Version lente
    N=1000000
    cpt=0
    for k in range(N):
        T = np.random.rand(n)
        cpt += np.sum(T**2) < 1
    return 2**n*cpt/N

def VolMC2(n): # Version un peu plus rapide
    N=1000000
    T=np.random.random((n,N))  
    return 2**n*np.mean([np.sum(T[:,k]**2) < 1 for k in range(N)])

# Formule exacte du volume d'une n-sphère

def Volume(n): 
    if n%2==0:
        f=np.prod(range(1,n//2+1))
        return np.pi**(n//2)/f
    else:
        f=np.prod(range(1,n+1,2))
        return 2**((n+1)//2)*np.pi**((n-1)//2)/f
      
# Versions avec calcul approché d'intégrales
      
def vol2(): # Cas n = 2 - Aire d'un cercle
    def f(x0,x1):
        return 1
    
    def r1(x1):
        return [0,np.sqrt(1-x1**2)]
        
    return 4*nquad(f,[r1,[0,1]])[0]
    
def vol3(): # Cas n = 3 - Volume d'une sphère de R^3
    def f(x0,x1,x2):
        return 1
    
    def r1(x1):
        return [0,np.sqrt(1-x1**2)]
        
    def r2(x1,x2):
        return [0,np.sqrt(1-x1**2-x2**2)]
        
    return 8*nquad(f,[r2,r1,[0,1]])[0]
    
def vol4(): # Cas n = 4 - Volume d'une 4-sphère
    def f(x0,x1,x2,x3):
        return 1
    
    def r1(x1):
        return [0,np.sqrt(1-x1**2)]
        
    def r2(x1,x2):
        return [0,np.sqrt(1-x1**2-x2**2)]
        
    def r3(x1,x2,x3):
        return [0,np.sqrt(1-x1**2-x2**2-x3**2)]
        
    return 16*nquad(f,[r3,r2,r1,[0,1]])[0]
    
def voln(n): # Cas général - Volume d'une n-sphère
    def f(*args):
        return 1
    
    def rn(*args):
        y=np.sum(np.array(args)**2)
        return [0,np.sqrt(1-y)]
    
    r=[rn]*(n-1)+[[0,1]]
    
    return 2**n*nquad(f,r)[0]

###################################################################
# Script principal de test
###################################################################

print('Versions Monte-Carlo')
print('La première (la plus lente)')
print((VolMC1(2),Volume(2)))
print((VolMC1(3),Volume(3)))
print((VolMC1(4),Volume(4)))
print()

print('La deuxième (un peu plus rapide)')
print((VolMC2(2),Volume(2)))
print((VolMC2(3),Volume(3)))
print((VolMC2(4),Volume(4)))
print()

print('Versions avec nquad:')
print('Première version (3 fonctions vol):')
print((vol2(),Volume(2)))
print((vol3(),Volume(3)))
print((vol4(),Volume(4)))
print()

print('Seconde version (1 fonction voln unique):')
print((voln(2),Volume(2)))
print((voln(3),Volume(3)))
print((voln(4),Volume(4)))