圆周率Π的数值计算，附Matlab代码和详细步骤

坚果橙

于 2023-12-25 18:15:00 发布

阅读量1k

点赞数 7

文章标签：算法

本文链接：https://blog.csdn.net/m0_54534574/article/details/135186337

版权

文章介绍了历史上计算圆周率的传统方法，如无穷级数和割圆术，然后重点讲述了蒲丰投针实验这一统计学方法，以及如何利用MATLAB实现的随机取样（蒙特卡洛法）来近似圆周率。作者通过模拟演示了这一过程和结果的准确性。

摘要由CSDN通过智能技术生成

对于圆周率 $\pi$ 的具体数值的计算，在过去一直是一个非常困难的问题。一般的思路都是试图将 $\pi$ 无限分割来逐渐逼近具体数值。无论欧拉的无穷级数思想还是阿基米德和刘徽的割圆术都是这个想法。

后来看到一个很“清奇”的计算方法——蒲丰投针实验法。利用统计学思想，来计算 $\pi$ 值。细看蒲丰投针实验，个人感觉其实验的一个思维亮点在于选取了针与平行线的夹角 $\theta$ 来作为横坐标，这样就引入了 $\pi$ ，非常巧妙。

蒲丰投针试验具体内容如下，假设一张白纸上有很多条距离为 $d$ 的平行线，现在有一根长度为 $l$ 的针，将其随机投下。那么则有针与平行线相交的概率 $p$ ：

$p = \frac{2}{\pi }\int_0^{\frac{\pi }{2}} {\frac{{l\sin \theta }}{d}} d\theta = \frac{{2l}}{{\pi d}}[ - \cos \theta |_0^{\frac{\pi }{2}}] = \frac{{2l}}{{\pi d}}$

其中，针与平行线的夹角 $\theta$ ， ${l\sin \theta }$ 是针在纸上的高，上式是对可能的角度取平均。在实验中，需要保证 $l\sin \theta \le d$ 。

后来也有很多学者重复了这个实验，用这种算概率的方法计算 $\pi$ 的具体值，如下图所示：

显然，这个方法也有其局限性——看上去十分笨拙。尤其在那个没有计算机的时代，需要人工一次又一次的投针，这样的方法真的看上去就是对时间的一种变相浪费。

我们一般称这种方法为随机取样法，或者“洋气”地称之为“蒙特卡洛法”（它源于世界著名的赌城——摩纳哥的Monte Carlo（蒙特卡洛））。

出于好玩，我自己用MATLAB，按照这个数学思想，构建了如下模型：

在坐标轴中选取2*2的正方形，在其中画一个半径为1的圆（如下图），然后在该正方形内进行随机取点（N），计算点在圆内的次数（n），正方形的面积（S），其中S=4，而Π对应的是其中圆的面积。

然后用该方法求出 $\pi$ 。

${\pi = \frac{n}{N}S}$

自己编写代码运行，进行模拟，得出结果。

步骤如下：

（1）首先打开Matlab，新建实时脚本。

（2）输入如下代码，取采样次数为10^8。

clc,clear
format long
N=1e8;
x=unifrnd(-1,1,[1,N]);
y=unifrnd(-1,1,[1,N]);
pinshu=sum(y.^2+x.^2<=1);
area_pi=pinshu/N*4

（3）得到结果如下：

这个结果还算是比较令人满意的。如果将N取得再大一些，理论上能够获得更加准确的 $\pi$ 值。

坚果橙

关注

7
点赞
踩
18

收藏

觉得还不错? 一键收藏
0
评论
圆周率Π的数值计算，附Matlab代码和详细步骤

从蒲丰投针到现代计算机的仿真模拟（蒙特卡洛法）
复制链接

扫一扫