本文写给那些python初学者与对蒙特卡洛算法感兴趣,但却不知该如何理解或应用的人。

(虽然我发现这个貌似有许多人做过了,但是程序都相对复杂,不便于理解,于是我就自己编写了一段程序,海龟的可视化请看下集如果还有的话。)

头歌python计算圆周率公式法 python 圆周率的计算_python中圆周率咋打

蒙特卡罗方法,也称统计模拟方法,是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。

蒙特卡罗方法可以将所求解的问题可以转化为某种随机分布的特征数,比如随机事件出现的概率,然后通过随机抽样的方法,以随机事件出现的频率估计其概率,并将其作为问题的解。

假设我们要计算一个不规则图形的面积,那么图形的不规则程度和分析性计算(比如,积分)的复杂程度是成正比的。蒙特卡罗方法基于这样的想法:假設你有一袋豆子,把豆子均匀地朝这个图形上撒,然后数这个图形之中有多少颗豆子,这个豆子的数目就是图形的面积。当你的豆子越小,撒的越多的时候,结果就越精确。借助计算机程序可以生成大量均匀分布坐标点,然后统计出图形内的点数,通过它们占总点数的比例和坐标点生成范围的面积就可以求出图形面积。

——维基百科

蒙特卡洛算法到底是啥?听上去有些复杂,但如果结合实例来看就会简单许多。

在“计算π”这个问题中,可以先把它转化为一个“随机撒豆子”问题。同上面wiki说的,撒豆子的地方的面积可以看成是豆子的数量,于是问题就变为了面积之比的问题。

所以,想要撒豆子,先需要把问题转化为图形的面积之比。我们直接使用一个”框架“:正方形与圆

头歌python计算圆周率公式法 python 圆周率的计算_随机数_02

圆与外切方(自处留一个思考题,还有什么图形可以用来描述圆周率)

为了求π,我们需要一个圆,还需要一个与圆有关的图形,这个图形的部分元素需要使用圆的部分元素,这个正方形刚好符合,因为正方形的边长=圆的直径。

这个圆的面积为πr²,方的面积为d²,也就是(2r)²。

所以用圆的面积比上方的面积为:

πr²÷(2r)²

=πr²÷4r²

=π÷4

所以用圆面积比上方面积再乘4就是π了。

下一步,需要用坐标描述撒出去的豆子的位置,用来判断豆子撒在了圆里还是方里。这里我们使用直角坐标系。(其实只要判断是否撒在圆里就行,因为我们是往方形里撒豆子。)

头歌python计算圆周率公式法 python 圆周率的计算_python中圆周率咋打_03

直角坐标系(此处还有思考题,还可以用什么样的坐标系呢?)

python:
import random                     #为了随机撒豆,引入随机库
r=5                                  #设定一个半径
fang=0                                 #统计豆子落在方形中的数量
yuan-0                                 #统计豆子落在圆形中的数量
#从此往下所有的程序是有关撒一个豆子的,需要重复循环
while True:                           #重复循环
x = random.uniform(-r, r) #设定豆子的x轴位置,random.uniform生成随机小数,后面的两个值表示生成的范围
y = random.uniform(-r, r)    #设定豆子的y轴位置,同上

该如何判断是不是撒在圆里面呢?

可以用豆子与原点的距离,也就是半径来判断。如果半径大于圆的半径,那么撒在圆外面;如果半径小于圆的半径,那么撒在圆的里面。

直角坐标系中点a(x1,y1)与点b(x2,y2)的距离有如下公式:

头歌python计算圆周率公式法 python 圆周率的计算_python中圆周率咋打_04

请使用勾股定理

如果将点a设为原点的话,那么x1=0,y1=0

可以用以下程序表示这个公式。

python:

a = (0-x) ** 2            #在python中 ** 表示幂运算,例如 a² 就是 a**2
b = (0-y) ** 2
c = (a+b) ** 0.5        #当 ** 后面是0.5时,表示求平方根
#接下来统计豆子位置!
if c>r:                      #当落在圆外时
fang=fang+1        #方里豆子数+1(其实也可以写成fang+=1)
elif c
yuan=yuan+1      #圆里豆子数+1
fang=fang+1        #圆里面也算是方的面积,所以 方里豆子数+1
else:                        #都不是的情况时(撒到圆上面了,超极小概率)
continue               #跳过
if fang==0:                #圆周率计算要用除法,除数是零会报错,当除数为零时
continue              #跳过
if yuan==0:
continue
pai=yuan/fang*4       #计算圆周率!
print(pai)                   #直接打印

运行起来应该会看到字符的疯狂滚动,你还可以把循环加个条件,让它运行一会就停下。

理论上只要撒无数个豆子,就可以无限接近圆周率。(我觉得这可能要看python的小数精度,但你也可以把半径设一个巨大的值,再把输出的π乘一个很大的数,把

部分小数部分变成整数。

程序截图:

头歌python计算圆周率公式法 python 圆周率的计算_python_05

海龟可视化与储存数据下次再写(如果还有下次)