Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import pylab
- import time
- from math import pi
- RSp = 1 # радиус сферы
- MaxD = 16 # максимальная размерность
- # реальные объемы с графика
- 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]
- # x - количество точек, D - число измерений
- def calcVNumPy(x, D):
- Xi = np.random.uniform(0, 1, x * D).reshape(x, D) # набор случайных точек - вектора из [0,1]
- if D > 1:
- Dist = np.sum(Xi**2, axis=1)**0.5
- PtsInSphere = np.sum(Dist < RSp) # число точек внутри сферы
- elif D == 1:
- PtsInSphere = x
- return 2**d * PtsInSphere / x # объем методом Монте Карло
- data = []
- for d in range(1, MaxD+1): # d - число измерений
- startTime = time.time() # время начало расчета
- n = 10**5
- m = 10**5
- Vn = calcVNumPy(n, d)
- while True:
- Vm = calcVNumPy(m, d)
- Vnm = (n * Vn + m * Vm) / (n + m)
- V1 = Vnm
- epsA = abs(V1 - Vreal[d]) # абсолютная погрешность
- epsR = epsA / Vreal[d] # относительная погрешность
- dV = abs(Vnm - Vn)
- s = epsA + abs(Vnm) * epsR
- if (dV == 0 and s == 0) or (epsA < 0.001 and dV / s < 1.0):
- break
- n += m
- Vn = Vnm
- endTime = time.time()
- # число измерений, количество точек, абс.погр., относит.погр., время в мс
- print("Время расчета для {:2d}-d t = {:.2f}с".format(d, endTime-startTime), end='. ')
- print("Значение: {:.7f}. Расчет: {:.7f}. Погрешности абс={:.6f}, относ={:.6f}".format(Vreal[d], V1, epsA, epsR))
- data.append([d, V1, n]) # сохраняем результаты
- D, V, P = np.array(data).T
- # График 1: расчитанный объем шара в зависимости от размерности
- fig, axs = pylab.subplots(2, 1, figsize=(10,20))
- axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
- axs[0].plot(D, Vreal[1:MaxD+1], 'C0', marker='o')
- axs[0].plot(D, V, 'C2', marker='x')
- # График 2: количество точек в зависимости от размерности пространства
- axs[1].set_title('Количество точек', color='C1')
- axs[1].plot(D, P, marker='x')
- pylab.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement