问小白 wenxiaobai
资讯
历史
科技
环境与自然
成长
游戏
财经
文学与艺术
美食
健康
家居
文化
情感
汽车
三农
军事
旅行
运动
教育
生活
星座命理

计算圆周率的五种方法

创作时间:
作者:
@小白创作中心

计算圆周率的五种方法

引用
CSDN
1.
https://blog.csdn.net/zhang_dengyu/article/details/136756168

圆周率π是一个在数学、物理乃至工程领域都具有重要地位的常数。自古以来,人们就一直在寻找更精确、更高效的计算π的方法。本文将介绍五种计算圆周率的经典方法:割圆法、无穷级数法、蒙特卡洛法、梅钦法和拉马努金法。每种方法都包含了相关知识背景、编程要求和代码实现,帮助读者全面了解这些算法的原理和应用。

割圆法

相关知识

割圆术是古代中国数学家刘徽在公元263年提出的一种计算圆周率的方法。刘徽形容他的“割圆术”说:“割之弥细,所失弥少,割之又割,以至于不可割,则与圆合体,而无所失矣。”这种方法包含了求极限的思想。通过求圆内接正多边形的周长来近似求圆的周长,并使正多边形的周长无限接近圆的周长,进而求得较为精确的圆周率。

设圆半径为1,圆内接正6边形边长也为1,可以如下计算正12边形的边长:

  • OC = √(OB² - BC²) = √(1 - (AB/2)²)
  • AB = OD = OA = 1
  • CD = 1 - OC
  • 新的边长AD值为:AD = √(AC² + CD²) = √((AB/2)² + CD²)

编程实现

import math

def cutting_circle(n):  # n为分割次数
    side_length = 1  # 初始边长
    edges = 6  # 初始边数

    def length(x):  # x为AB
        h = math.sqrt(1 - (x / 2) ** 2)  # h为OC
        return math.sqrt((x / 2) ** 2 + (1 - h) ** 2)

    for i in range(n):
        side_length = length(side_length)
        edges *= 2
        pi = side_length * edges / 2
    return edges, pi

if __name__ == '__main__':
    times = int(input())  # 割圆次数
    edges, pi = cutting_circle(times)  # 调用函数返回值
    print(f'分割{times}次,边数为{edges},圆周率为{pi:.6f}')
    print(f'math库中的圆周率常量值为{math.pi:.6f}')

测试说明

平台会对你编写的代码进行测试:

  • 输入格式:输入为一行,一个正整数n,表示在内接正六边形的基础上割圆的次数。
  • 输出格式:输出为两行,第一行输出在内接正六边形基础上再次分隔n次,分隔得到的正多边形的边数m,以及割圆法计算得到的圆周率pi(保留六位小数);第二行输出math库里的pi值(保留六位小数)。具体输出格式见示例。

示例

输入:

4

输出:

分割4次,边数为96,圆周率为3.141032
math库中的圆周率常量值为3.141593

无穷级数法

π是个超越数,圆周率的超越性否定了化圆为方这种尺规作图精确求解问题的可能性。有趣的是,π可以用无穷级数表示:

1/1 - 1/3 + 1/5 - 1/7 + ⋯ = π/4

左边的展式是一个无穷级数,被称为莱布尼茨级数(Leibniz),这个级数收敛到π/4,它通常也被称为格雷戈里-莱布尼茨级数,用以纪念莱布尼茨同时代的天文学家兼数学家詹姆斯·格雷戈里。

编程实现

def leibniz_of_pi(error):
    """接收用户输入的浮点数阈值为参数,返回圆周率值"""
    f = 1
    i = 1
    result = 0  # 初始化
    while 1 / (2 * i - 1) >= error:
        result += f * 1 / (2 * i - 1)
        f = -f
        i += 1
    return result * 4

if __name__ == '__main__':
    threshold = float(input())
    print("{:.8f}".format(leibniz_of_pi(threshold)))  # 保留小数点后八位

测试说明

平台会对你编写的代码进行测试:

  • 输入格式:输入为一行,一个正浮点数threshold,用以限定级数法求解pi值时,迭代累加只加绝对值大于threshold的项。
  • 输出格式:输出为一个浮点数,是程序使用级数法求解的pi值,要求保留小数点后八位。

示例

输入:

0.000002

输出:

3.14158865

输入:

1e-6

输出:

3.14159065

蒙特卡洛法

蒙特卡洛(Monte Carlo)方法是由数学家冯·诺伊曼提出的,诞生于上世纪40年代美国的“曼哈顿计划”。蒙特卡洛是一个地名,位于赌城摩纳哥,象征概率。蒙特卡洛方法的原理是通过大量随机样本,去了解一个系统,进而得到所要计算的值。

用蒙特卡洛方法计算圆周率π的原理如下:一个边长为2r的正方形内部相切一个半径为r的圆,圆的面积是πr²,正方形的面积为4r²,二者面积之比是π/4,因为比值与r大小无关,所以可以假设半径r的值为1。

在这个正方形内部,随机产生n个点,坐标为(x,y),当随机点较多时,可以认为这些点服从均匀分布的规律。计算每个点与中心点的距离是否小于圆的半径(x²+y²<r²),以此判断是否落在圆的内部。统计圆内的点数c,c与n的比值乘以4,就是π的值。理论上,n越大,计算的π值越准,但由于随机数不能保证完全均匀分布,所以蒙特卡洛法每次计算结果可能不同。

随机数种子
计算机中的随机数是通过算法产生,相同的随机数种子可以帮助我们获得相同的随机数序列,所以自动评测的题目都要求使用随机数种子。本题中只需要执行以下语句就可以。

sd = int(input())  # 读入随机数种子
random.seed(sd)  # 设置随机数种子

编程实现

import random

def monte_carlo_pi(num):
    """接收正整数为参数,表示随机点的数量,利用蒙特卡洛方法计算圆周率
    返回值为表示圆周率的浮点数"""
    c = 0
    for _ in range(num):
        x, y = random.uniform(-1, 1), random.uniform(-1, 1)
        if x * x + y * y <= 1:
            c += 1
    time = c * 4 / num
    return time

if __name__ == '__main__':
    sd = int(input())  # 读入随机数种子
    random.seed(sd)  # 设置随机数种子
    times = int(input())  # 输入正整数,表示产生点数量
    print(monte_carlo_pi(times))  # 输出圆周率值,浮点数

测试说明

平台会对你编写的代码进行测试:

示例

输入:

100
100000

输出:

3.15316

梅钦法

1706年,英国人约翰·梅钦(John Machin)发明了一个用于计算π/4值的公式:

4π = 4arctan(1/5) - arctan(1/239)

梅钦公式是格里高利/莱布尼茨计算的公式的变体,但是更实用,它的收敛速度显著增加,这使得它成为了更实用的计算的方法,虽然有若干种类梅钦(Machin-like)公式,但梅钦公式仍然是计算值的主要公式。

编程实现

from math import *

def machin_of_pi():
    """用梅钦级数计算圆周率,返回圆周率值"""
    pi = 4 * (4 * atan(1 / 5) - atan(1 / 239))
    return pi

if __name__ == '__main__':
    cal_pi = machin_of_pi()  # 调用判断类型的函数
    print(cal_pi)  # 输出函数运行结果

测试说明

平台会对你编写的代码进行测试:

输入:无

输出:

3.1415926535897936

拉马努金法

拉马努金曾经提出过很多关于求π的公式,这些公式都有以下几个特点:等号右边的构造超乎常人想象,收敛速度极快!

π = 9801 / (2 * sqrt(2)) * ∑(k=0 to ∞) [(k!)^4 * (1103 + 26390k)] / [(396^(4k)) * (4k)!]

输入一个正整数n,根据些公式计算累加n次时的圆周率值。

编程要求

根据提示,在右侧编辑器补充代码,调用math库中的阶乘函数计算并输出圆周率的大小。

代码实现

from math import *

def ramanujan_of_pi(n):
    """接收一个正整数n为参数,用拉马努金公式的前n项计算圆周率并返回。"""
    c = 0
    b = 2 * sqrt(2) / 9801
    for i in range(n):
        a = pow(factorial(i), 4) * pow(396, 4 * i)  # 分母
        c += (factorial(4 * i) * (1103 + 26390 * i)) / a

    pi = 1 / (b * c)
    return pi

if __name__ == '__main__':
    n = int(input())
    cal_pi = ramanujan_of_pi(n)
    print(cal_pi)  # 输出函数运行结果

测试说明

平台会对你编写的代码进行测试:

测试输入:

1

预期输出:

3.1415927300133055
© 2023 北京元石科技有限公司 ◎ 京公网安备 11010802042949号