描述
1.割圆法
国魏晋时期的数学家刘徽于公元263年撰写《九章算术注》,其中就有数学史上著名的“割圆术”。刘徽形容他的“割圆术”说:割之弥细,所失弥少,割之又割,以至于不可割,则与圆合体,而无所失矣。这包含了求极限的思想。通过求圆内接正多边形的周长来近似求圆的周长,并使正多边形的周长无限接近圆的周长,进而求得较为精确的圆周率。
设圆半径为1,圆内接正6边形边长也为1,可如下计算正12边形的边长:
2.莱布尼茨级数法
π是个超越数,圆周率的超越性否定了化圆为方这种尺规作图精确求解问题的可能性。有趣的是,π可以用无穷级数表示:
左边的展式是一个无穷级数,被称为莱布尼茨级数(Leibniz),这个级数收敛到π/4,它通常也被称为格雷戈里-莱布尼茨级数,用以纪念莱布尼茨同时代的天文学家兼数学家詹姆斯·格雷戈里。
编程用这个公式计算π值,输入一个小数作为阈值,当最后一项的绝对值小于给定阈值时停止计算并输出得到的π值。
- 蒙特卡洛法
蒙特卡洛(Monte Carlo)方法是由数学家冯·诺伊曼提出的,诞生于上世纪40年代美国的“曼哈顿计划”。蒙特卡洛是一个地名,位于赌城摩纳哥,象征概率。蒙特卡洛方法的原理是通过大量随机样本,去了解一个系统,进而得到所要计算的值。
用蒙特卡洛方法计算圆周率π的原理如下:一个边长为2r的正方形内部相切一个半径为r的圆,圆的面积是πr2,正方形的面积为4r2,二者面积之比是π/4,因为比值与r大小无关,所以可以假设半径 r的值为1。
在这个正方形内部,随机产生n个点,坐标为(x,y),当随机点较多时,可以认为这些点服从均匀分布的规律。计算每个点与中心点的距离是否大于圆的半径(x2+y2>r2),以此判断是否落在圆的内部。统计圆内的点数c,c与n的比值乘以4,就是π的值。理论上,n越大,计算的π值越准,但由于随机数不能保证完全均匀分布,所以蒙特卡洛法每次计算结果可能不同。
编程实现用蒙特卡洛方法计算π值,为了自动测评的需要,请先读入一个正整数sd作为随机数种子,并要求使用 x,y = random.uniform(-1,1) , random.uniform(-1,1) 语句来生成随机点的坐标值。
4. 梅钦法
梅钦公式是格里高利/莱布尼茨计算的公式的变体,但是更实用,它的收敛速度显著增加,这使得它成为了更实用的计算的方法,虽然有若干种类梅钦(Machin-like)公式,但梅钦公式至今仍然是计算值的主要公式,根据此公式计算圆周率值。
5. 拉马努金法
拉马努金曾经提出过很多关于求π的公式,这些公式都有以下几个特点: 等号右边的构造超乎常人想象,收敛速度极快!
比如这个拉马努金在1914年发布的以他自己名字命名著名公式。
输入一个正整数n,根据此公式计算累加n次时的圆周率值。
示例 1
输入:
割圆法
4
输出: 3.14103195089051
示例 2
输入:
无穷级数法
0.000002
输出: 3.141588653589781
示例 3
输入:
蒙特卡洛法
100000
100
输出: 3.15316
示例 4
输入:
梅钦法
输出: 3.1415926535897936
示例 5
输入:
拉马努金法
1
输出: 3.1415927300133055
# -------- ------- --------
# @File : 计算圆周率.py
# @Author : 赵广辉
# @Contact: vasp@qq.com
# @Company: 武汉理工大学
# @Version: 1.0
# @Modify : 2021/10/30 23:18
# Python程序设计基础,高等教育出版社
# -------- ------- --------
import math
import random
import matplotlib.pyplot as plt
def type_judge(pi_type):
"""接收一个字符串为参数,根据参数调用相应函数计算圆周率。 """
if pi_type == '割圆法':
times = int(input()) # 输入一个表示边数量的正整数
return cutting_circle(times) # 调用函数计算圆周率
elif pi_type == '无穷级数法':
threshold = float(input()) # 输入转为浮点数
return leibniz_of_pi(threshold) # 调用函数计算圆周率
elif pi_type == '蒙特卡洛法':
num = int(input()) # 输入转为整数
s = int(input()) # 输入随机数种子
return monte_carlo_pi(num, s) # 调用函数计算圆周率
elif pi_type == '梅钦法':
return machin_of_pi() # 调用函数计算圆周率
elif pi_type == '拉马努金法':
num = int(input()) # 输入转为整数
return ramanujan_of_pi(num)
else:
return f'未找到{pi_type}计算方法'
def cutting_circle(times):
"""参数times为分割次数
π = 周长/(2*圆的半径)得到π的近似值。
# 半径为1的圆内接正6边形边长也是1
# 边长 side_length
# 半径 radius
# 圆周率 pi
# 三角形的高 height
"""
side_length = 1 # 初始边长
edges = 6 # 初始边数
for i in range(times):
height = 1 - math.sqrt(1 - (side_length / 2) ** 2)
side_length = math.sqrt(height ** 2 + (side_length / 2) ** 2)
edges = edges * 2 # 每割一次,边数量加倍
pi = side_length * edges / 2
return pi
def leibniz_of_pi(error):
"""接收用户输入的浮点数阈值为参数,用格雷戈里-莱布尼茨级数计算圆周率,返回圆周率值"""
n, quarter_of_pi, sign = 1, 0, 1
while 1 / (2 * n - 1) > error:
quarter_of_pi = quarter_of_pi + sign * 1 / (2 * n - 1)
n = n + 1
sign = -sign
pi = 4 * quarter_of_pi
return pi
def monte_carlo_pi(num, s):
"""接收一个表示随机次数的整数和一个整数做随机数种子,用蒙特卡洛法计算圆周率,返回一个浮点数"""
random.seed(s)
hits = 0 # 落在圆内的计数器初值设为 0
for i in range(1, num + 1):
x, y = random.uniform(-1, 1), random.uniform(-1, 1), # 生成两个随机数模拟一个点的坐标
pos = (x ** 2 + y ** 2) ** 0.5 # 计算坐标(x,y)到原点的距离
if pos <= 1.0: # 如果距离小于等于1,点在圆内
hits = hits + 1
pi = 4 * (hits / num)
return pi
def machin_of_pi():
"""用梅钦级数计算圆周率,返回圆周率值"""
quarter_of_pi = 4 * math.atan(1/5) - math.atan(1/239)
pi = 4 * quarter_of_pi
return pi
def ramanujan_of_pi(n):
"""接收一个正整数n为参数,用拉马努金公式的前n项计算圆周率并返回。"""
result = 0
for i in range(n):
result = result + ((2 * 2 ** 0.5)) * (math.factorial(4 * i) * (1103 + 26390 * i)) / (math.pow(math.factorial(i), 4) * math.pow(396, 4 * i))/9801
pi = 1 / result
return pi
# 取消以下代码行前面的注释可以调用其中的函数对割圆法和蒙特卡洛法进行演示
# def draw_circle(r, side_num):
# """创建图形和轴,用折线图绘制正多边形。
# @参数 r:圆的半径
# @参数 side_num:正多边形的边数
# """
# fig, ax = plt.subplots() # 创建图形和轴
# plt.subplots_adjust(left=0.1, bottom=0.25)
# x, y = xy_of_polygon(r, side_num)
# plt.plot(x, y, lw=2, color='red') # 设置线宽和颜色
# ax.set_aspect('equal') # 设置坐标轴纵横比相等
# ax.set_xlim(-r-5, r+5) # 设置x轴刻度起止值
# ax.set_ylim(-r-5, r+5)
# plt.draw() # 重新绘制多边形
# plt.show() # 显示图形
#
#
# def cutting_circle(times):
# """
# 接收表示分割次数的整数n 为参数,计算分割n 次时正多边形的边数和圆周率值,返回边数和圆周率值。
# @参数 times:分割次数
# π = 周长/(2*圆的半径)得到π的近似值。
# # 半径为1的圆内接正6边形边长也是1
# # 边长 side_length
# # 半径 radius
# # 圆周率 pi
# # 三角形的高 height
# >>> cutting_circle(4)
# 3.14103195089051
# """
# side_length = 1 # 初始边长
# edges = 6 # 初始边数
# for i in range(times):
# height = 1 - math.sqrt(1 - (side_length / 2) ** 2)
# side_length = math.sqrt(height ** 2 + (side_length / 2) ** 2)
# edges = edges * 2 # 每割一次,边数量加倍
# draw_circle(300, edges) # 调用此函数可演示割圆效果,每次循环绘制一个近似圆,关闭后看下一个
# pi = side_length * edges / 2
# return pi
#
#
# def monte_carlo_show(in_circle_lst, out_circle_lst):
# """
# @参数 in_circle_lst:落在圆内的点的坐标列表
# @参数 out_circle_lst:落在圆外的点的坐标列表
#
# """
# plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
# plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
# plt.figure(figsize=(11, 11)) # 设置画布长宽
# x_in_circle = [x[0] for x in in_circle_lst] # 落在圆内的点的列表x
# y_in_circle = [x[1] for x in in_circle_lst] # 落在圆内的点的列表y
# x_out_circle = [x[0] for x in out_circle_lst] # 落在圆外的点的列表x
# y_out_circle = [x[1] for x in out_circle_lst] # 落在圆外的点的列表y
# plt.scatter(x_out_circle, y_out_circle, s=10, facecolors='blue') # 绘制散点,落在圆外在点颜色用蓝色
# plt.scatter(x_in_circle, y_in_circle, s=5, facecolors='red') # 绘制散点,落在圆内在点颜色用红色
# plt.title('蒙特卡洛法计算圆周率演示') # 图的标题
# plt.show() # 显示绘制结果
#
#
# # 给出turtle模板,学生填monte_carlo代码,帮助学生理解算法,了解turtle绘图
# def monte_carlo_pi_turtle(n, s):
# """用turtle绘图模拟蒙特卡洛n次结果"""
# random.seed(s)
# turtle.tracer(1000)
# turtle.pensize(4)
# turtle.penup()
# turtle.goto(-300, -300)
# turtle.pendown()
# for i in range(4):
# turtle.forward(600)
# turtle.left(90)
# turtle.forward(300)
# turtle.pencolor('green')
# turtle.circle(300)
# hits = 0 # 落在圆内的计数器初值设为 0
# for i in range(1, n + 1):
# x, y = random.uniform(-1, 1), random.uniform(-1, 1) # 生成两个随机数模拟一个点的坐标
# pos = (x ** 2 + y ** 2) ** 0.5 # 计算坐标(x,y)到原点的距离
# if pos <= 1.0: # 如果距离小于等于1,点在圆内
# hits = hits + 1
# turtle.pencolor('red') # 落在圆内在点颜色用红色
# else: # 如果距离大于1,点在圆外
# turtle.pencolor('blue') # 落在圆内在点颜色用蓝色
# turtle.penup()
# turtle.goto(x * 300, y * 300) # 画笔抬起并移动到数值放大400倍的x,y处
# turtle.pendown()
# turtle.dot(3) # 画一个半径为 3 的圆点
# if i % 10000 == 0: # 实验为10000的倍数次时
# turtle.pencolor('black')
# pi = 4 * (hits / i) # 根据落在圆内外的点数量估算PI值
# turtle.penup() # 画笔抬起
# turtle.goto(320, 150 - i // 1000 * 30) # 移动到区域外记录当前PI值
# turtle.pendown() # 画笔抬起
# turtle.write("{}次时PI的值是{:.4f}".format(i, pi), font=("宋体", 18, "normal"))
# turtle.hideturtle() # 隐藏光标
# turtle.update() # 刷新
# turtle.done() # 结束绘制
if __name__ == '__main__':
type_of_pi = input() # 接收用户输入的字符串
cal_pi = type_judge(type_of_pi) # 调用判断类型的函数
print(cal_pi) # 输出函数运行结果