计算圆周率
圆周率的计算,伴随了很长的历史时期。不同时期先后曾有各种大牛使用不同的方法,都计算过圆周率。具体历史我们 我们不在这里赘述,大家可以参考百度百科 或者维基百科(伟大的防火墙可能不想让你知道)。
今天我们当然是看看计算机如何计算圆周率,一种最简单的就是蒙特卡洛法,其实就是增加随机实验的次数,不断逼近 圆周率。中学内容肯定学过了,抛大头钉的实验,是不是有印象了。具体做法如下:
- 我们先画一个边长为1的正方形,(具体几厘米可以自选),然后画出内接圆。
- 往这个正方形中抛大头钉,分别数出正方形中大头钉的数目N,圆中大头钉的数目n。
圆的面积为S = pi/4 ,正方形的面积就是1。落入圆中的概率为n/N,这个值应该等于圆的面积与正方向的面积之比。 于是得到圆周率pi的计算公式:
$$ \pi = 4 \frac{n}{N}$$
为了编写程序简单,我们考虑四分之一圆,并且把这个四分之一圆放在第一象限,抛掷大头钉可以使用(0,1)的 随机数模拟。产生一对(0,1)的随机数。分别判断这个点是否在圆内。首先我们写判断圆内的函数。
import math
import random
def is_in_circle():
x, y = random.random(), random.random()
y_ = math.sqrt(1 - x ** 2)
stat = False
if y < y_:
stat = True
return stat
数学模块math
之前已经用过了,随机数模块random
的random
方法产生一个[0,1)的随机数。
初始的状态默认为stat=False
,如果y<y_
则表示在圆内,更新stat
的状态,否则就返回初始值False
。
每次运行这个函数都会模拟抛掷一次大头钉实验。
接着我们写统计函数,返回pi的计算值。
def count(num=10000):
# num 总的循环次数
in_circle = 0 # 圆内计数器
for i in range(1, num + 1): # 注意range函数终点值不会取得需要加1
stat = is_in_circle() # 模拟实验
if stat:
in_circle += 1
if i % 1000 == 0:
print(4 * in_circle / num) # 每隔1000次输出圆周率
return 4 * in_circle / num
count()
输出结果为:
0.3172
0.6372
0.9416
1.2588
1.5736
1.8888
2.2036
2.5128
2.8252
3.1444
每次输出结果都不一样,当循环次数较少时,结果差别很大,次数很大比如超过10万以后,最终结果较稳定, 读者可以自行完成程序,并且尝试改变实验的次数,观察结果的变化。 完整的程序如下:
import math
import random
def is_in_circle():
x, y = random.random(), random.random()
y_ = math.sqrt(1 - x ** 2)
stat = False
if y < y_:
stat = True
return stat
def count(num=10000):
"""
:param num: 总的循环次数
"""
in_circle = 0 # 圆内计数器
for i in range(1, num + 1): # 注意range函数终点值不会取得需要加1
stat = is_in_circle()
if stat:
in_circle += 1
if i % 1000 == 0:
print(4 * in_circle / num) # 每隔1000次输出圆周率
return 4 * in_circle / num
count(10000)