使用程序计算近似Π值

一、前言

现在大多数语言，只需要调用一下Math.PI就可以知道Π值了。但是你有没有想过这个PI是怎么来的，是直接存储吗？还是计算来的。虽然不知道具体是怎么实现的，但是我们可以使用一些简单的数学知识，来计算出近似的Π值。

二、实现原理

我们小学就学过圆的面积公式，只不过那个时候我们直接使用3.14作为Π。圆的面积公式如下：

S = \pi R^2

这个公式大家都熟悉，如果我们能知道R和S就能求出Π。那么除了上面的方法，还有什么方法可以根据R计算S呢，有一种可以参考的方法就是使用微积分的思想，即把圆拆分成无数个小矩形，不过在计算机中我们只能拆分出有限个小矩形。如下图：

图片2.png

这里先考虑1/4个圆。我们把1/4个圆拆分成n个小矩形，计算n个矩形的面积并相加。当我们的n足够大时，n个矩形的面积和就会趋近圆的面积。我们假设n个矩形面积和为A，那么近似的Π计算公式如下：

\pi = \frac{4A}{R^2}

现在我们知道了pi的具体公式，接下来我们看看A的计算。

A的计算主要是矩形面积的累加，假如我们把圆拆成n个矩形，每个矩形的宽都为1/n，那么第i个矩形的面积应该为：

dS = dx * y

因为是半圆，因此y的计算公式如下：

y = \sqrt{R^2 - (x-R)^2}

其中x应该为i/n*R（因为是把R分成n份，第i个矩形的x=i/n*R）。最后，n个矩形相加的公式为：

A = \sum_{i=1}^n\frac{\sqrt{R^2 - (\frac{i}{n}R-R)^2}}{n}

下面我们就可以根据公式用程序求出Π的近似值。

三、代码实现

我们以单位圆为例，将圆（1/4个）拆分成n份，拆分n份我们可以用循环来实现。具体代码如下：

n = 1000
A = 0	# 矩形面积
r = 1	# 半径
for i in range(1, n+1):
    dx = 1 / n	# 拆成n份，每一份x为1/n
    y = pow(pow(r, 2) - pow(i*r/n-r, 2), 0.5)	# 使用公式计算y，其中x=i/n
    dS = dx * y	# 单个小矩形的面积
    A += dS	# 对矩形面积进行累加
pi = 4 * A / pow(r, 2) # 计算Π
print(pi)

上面程序运行结果如下：

3.1410880051481085

结果的很接近圆周率。当我们把n增大时，比如设置为100000，结果为：

3.1415876489413512

当我们不断增大n时，结果会越来越趋近Π。