圆周率π是一个无理数,没有任何一个精确公式能够计算π值,π的计算只能采用近似算法。
国际公认的π值计算采用蒙特卡洛方法。
蒙特卡洛(Monte Carlo)方法,又称随机抽样或统计试验方法。当所求解的问题是某种事件出现的概率,或某随机变量的期望值时,可以通过某种“试验”方法求解。
简单说,蒙特卡洛是利用随机试验求解问题的方法。
构造一个单位正方形和一个单位圆的1/4,往整个区域内随机投入点,根据点到原点的距离判断点是落在1/4的圆内还是在圆外,从而根据落在两个不同区域的点的数目,求出两个区域的比值。如此一来,就可以求出1/4单位圆的面积,从而求出圆周率π。
下面是用Python的实现:
# pi.py
from random import random
from math import sqrt
from time import clock
DARTS = 12000 # 总的实验次数
hits = 0
clock()
for i in range(1, DARTS):
x, y = random(), random();
dist = sqrt(x**2 + y**2)
if dist <= 1.0:
hits = hits + 1
pi = 4 * (hits/DARTS)
print('Value of PI is %s' % pi)
print('Total runtime : %-5.5ss' % clock())
代码中用到了random和math库的random()函数和sqrt()函数,为了统计时间,还用到了time库的clock()函数。
投入的点越多,计算值越精确。
另一个 Python 程序:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# 投点次数
n = 10000
# 圆的信息
r = 1.0 # 半径
a, b = (0., 0.) # 圆心
# 正方形区域边界
x_min, x_max = a-r, a+r
y_min, y_max = b-r, b+r
# 在正方形区域内随机投点
x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)
# 计算 点到圆心的距离
d = np.sqrt((x-a)**2 + (y-b)**2)
# 统计 落在圆内的点的数目
res = sum(np.where(d < r, 1, 0))
# 计算 pi 的近似值(Monte Carlo方法的精髓:用统计值去近似真实值)
pi = 4 * res / n
print('pi: ', pi)
# 画个图看看
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') # 防止图像变形
circle = Circle(xy=(a,b), radius=r, alpha=0.5)
axes.add_patch(circle)
plt.show()
某次运行截图:
蒙特卡洛方法提供了一个利用计算机中随机数和随机试验来解决现实中无法用公式求解问题的思路,广泛应用在金融工程学、宏观经济学、计算物理学等领域。
利用蒙特卡洛(Monte Carlo)方法计算π值[ 转载]
原文:https://www.cnblogs.com/juking/p/9498884.html