近似计算圆周率
正方形内部有一个相切的圆,它们的面积之比是π/4。
由圆面积和正方形面积公式可得:
S圆面积S正方形=πr2(2r)2=π4
现在向上面的正方形区域均匀的随机撒点,通过计算落在圆内的点数与总的点数来估计 π
# -*- coding:utf8 -*-
import numpy as np
import matplotlib.pyplot as plt
def pi_exp(R, N):
x = np.random.uniform(-R * 1.0, R * 1.0, (N,))
y = np.random.uniform(-R * 1.0, R * 1.0, (N,))
incircle = np.sum(np.power(x, 2) + np.power(y, 2) <= (R * 1.0)**2)
error = np.pi - (4.0 * incircle / N)
print error
######################################
theta = np.linspace(0, 2*np.pi, 3600)
cx = R * np.cos(theta)
cy = R * np.sin(theta)
plt.plot(cx,cy, color="r")
######################################
rx = np.asarray([-1, 1, 1, -1, -1] *R, dtype = np.float)
ry = np.asarray([1, 1, -1, -1, 1] *R, dtype = np.float)
plt.plot(rx,ry, color="b")
######################################
cl = ["r" if i**2 + j**2 <= R else "b" for i, j in zip(x,y)]
plt.scatter(x, y, c=cl)
plt.text(-1.2, 1.2, s = "R={R}, N={N}, error={error}, pi={pi}".format(R=R, N=N, error=error, pi=np.pi))
plt.show()
if __name__ == "__main__":
pi_exp(1, 10000)
运行结果如下图:
随着采样数的增大,预测的 π 与真实 π 之间的误差会越来越小。
参考:http://www.ruanyifeng.com/blog/2015/07/monte-carlo-method.html