zhukov000

HyperShpere

Oct 15th, 2021
745
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import pylab
  3. import time
  4. from math import pi
  5.  
  6. RSp = 1     # радиус сферы
  7. MaxD = 16   # максимальная размерность
  8.  
  9. # реальные объемы с графика
  10. Vreal = [1, 2, pi, 4*pi/3, pi**2/2, 8*pi**2/15, pi**3/6, 16*pi**3/105, pi**4/24, 32*pi**4/945, pi**5/120, 64*pi**5/10395, pi**6/720, 128*pi**6/135135, pi**7/5040, 15*pi**7/118771, 0.23533063035889312]
  11.  
  12. # x - количество точек, D - число измерений
  13. def calcVNumPy(x, D):
  14.     Xi = np.random.uniform(0, 1, x * D).reshape(x, D)    # набор случайных точек - вектора из [0,1]
  15.     if D > 1:
  16.       Dist = np.sum(Xi**2, axis=1)**0.5                    
  17.       PtsInSphere = np.sum(Dist < RSp)                     # число точек внутри сферы
  18.     elif D == 1:
  19.       PtsInSphere = x
  20.    
  21.     return 2**d * PtsInSphere / x                        # объем методом Монте Карло
  22.  
  23. data = []
  24.  
  25. for d in range(1, MaxD+1):           # d - число измерений
  26.     startTime = time.time()          # время начало расчета
  27.     n = 10**5
  28.     m = 10**5
  29.     Vn = calcVNumPy(n, d)
  30.     while True:
  31.           Vm = calcVNumPy(m, d)
  32.           Vnm = (n * Vn + m * Vm) / (n + m)
  33.  
  34.           V1 = Vnm
  35.           epsA = abs(V1 - Vreal[d])          # абсолютная погрешность
  36.           epsR = epsA / Vreal[d]             # относительная погрешность
  37.          
  38.           dV = abs(Vnm - Vn)
  39.           s = epsA + abs(Vnm) * epsR
  40.  
  41.           if (dV == 0 and s == 0) or (epsA < 0.001 and dV / s < 1.0):
  42.             break
  43.           n += m
  44.           Vn = Vnm
  45.     endTime = time.time()
  46.    
  47.     # число измерений, количество точек, абс.погр., относит.погр., время в мс
  48.     print("Время расчета для {:2d}-d t = {:.2f}с".format(d, endTime-startTime), end='. ')
  49.     print("Значение: {:.7f}. Расчет: {:.7f}. Погрешности абс={:.6f}, относ={:.6f}".format(Vreal[d], V1, epsA, epsR))
  50.     data.append([d, V1, n]) # сохраняем результаты
  51.    
  52. D, V, P = np.array(data).T
  53. # График 1: расчитанный объем шара в зависимости от размерности
  54. fig, axs = pylab.subplots(2, 1, figsize=(10,20))
  55. axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
  56. axs[0].plot(D, Vreal[1:MaxD+1], 'C0', marker='o')
  57. axs[0].plot(D, V, 'C2', marker='x')
  58. # График 2: количество точек в зависимости от размерности пространства
  59. axs[1].set_title('Количество точек', color='C1')
  60. axs[1].plot(D, P, marker='x')
  61.  
  62. pylab.show()
RAW Paste Data