高等数学：圆周率的计算方法

本文将对圆周率 $\pi$的计算方法作简单整理， $\pi$的重要性不必多说，它的计算促进了诸多领域的发展，本文将通过割圆术、无穷级数(泰勒展开)、微积分、概率学、连分数等方法进行计算，同时给出了部分公式的证明，之后用Matlab对部分公式进行了验证，并在文末给出了所编写的函数。

圆周率$\pi$的小数点后前314位为：

$3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063$

1 割圆术

刘徽 (225-295)，祖冲之(429-500)
利用圆内接正多边形及正多边形每条边与圆所延伸出的矩形得到圆周率上界和下界，从而圆周率近似值。割圆过程如下图所示。当然，刘徽是从正六边形开始的，然后边长依次乘2，直至正3072边形，得到了圆周率上界为3.1416，下界为3.1415。

割圆术逼近示意图1
下面计算圆周率的上界和下界。设圆的半径为r，则图中橙色部分面积即圆内接正n边形的面积可计算为：

$$\begin{array}{rl}{S}_{\text{橙}}& =n\cdot \frac{1}{2}\cdot 2r\cos \frac{\pi -\frac{2\pi }{n}}{2}\cdot r\sin \frac{\pi -\frac{2\pi }{n}}{2}\\ & =\frac{n}{2}{r}^2\sin \left(\pi -\frac{2\pi }{n}\right)\\ & =\frac{n}{2}{r}^2\sin \frac{2\pi }{n}\end{array}$$

图中绿色部分面积可表示为：
$$\begin{array}{rl}{S}_{\text{绿}}& =S_{\text{橙}}+n\cdot 2r\cos \frac{\pi -\frac{2\pi }{n}}{2}\cdot \left(r-r\sin \frac{\pi -\frac{2\pi }{n}}{2}\right)\\ & =S_n+2n{r}^2\sin \frac{\pi }{n}-n{r}^2\sin \frac{2\pi }{n}\\ & =\frac{n}{2}{r}^2\sin \frac{2\pi }{n}+2n{r}^2\sin \frac{\pi }{n}-n{r}^2\sin \frac{2\pi }{n}\\ & =2n{r}^2\sin \frac{\pi }{n}-\frac{n}{2}{r}^2\sin \frac{2\pi }{n}\end{array}$$
显然，由于$S_{\text{橙}}<S_{\text{圆}}<S_{\text{绿}}$，即：
$$\frac{n}{2}{r}^2\sin \frac{2\pi }{n}<\pi r^2<2n{r}^2\sin \frac{\pi }{n}-\frac{n}{2}{r}^2\sin \frac{2\pi }{n}$$
化简得：
$$\frac{n}{2}\sin \frac{2\pi }{n}<\pi <2n\sin \frac{\pi }{n}-\frac{n}{2}\sin \frac{2\pi }{n}$$
下面将$n=3,6,12,\dots ,3072$边形对应的$\pi$的近似值绘制如下，显然，n越大，$\pi$的上下界近似值约精确。当$n=3072$时，有：
$$3.14159046322805<\pi <3.14159374877049$$
同时，上界和下界与圆周率$\pi$真实值的相对偏差变化情况如下，可见当n=3072时，相对偏差量级为\text{1e-7}。

割圆术对应的圆周率上下界和真实值对比

割圆术对应的圆周率上下界和真实值的相对偏差
另一方面，利用极限${\lim }_{x\to 0}\frac{\sin x}{x}=1$可得：
$$\begin{gathered} \underset{n\to \infty }{\lim }\frac{n}{2}\sin \frac{2\pi }{n}=\pi \frac{\sin \frac{2\pi }{n}}{\frac{2\pi }{n}}=\pi \\ \underset{n\to \infty }{\lim }\left(2n\sin \frac{\pi }{n}-\frac{n}{2}\sin \frac{2\pi }{n}\right)=\underset{n\to \infty }{\lim }\pi \left(2\frac{\sin \frac{\pi }{n}}{\frac{\pi }{n}}-\frac{\sin \frac{2\pi }{n}}{\frac{2\pi }{n}}\right)=\pi \end{gathered}$$
因此，$n\to \infty$时，即可得到圆周率。

小结：割圆术已经有点极限的意思了，除了利用上述形式的面积来近似，也可通过圆内接正多边形和外切正多边形的面积或周长来进行近似，下面画出利用圆内接正多边形和外切正多边形近似的示意图。

割圆术逼近示意图2
显然，对于面积有$S_{\text{扇形}OAB}<S_{△OCD}$因此：
KaTeX parse error: Undefined control sequence: \overparen at position 10: {A B} < \̲o̲v̲e̲r̲p̲a̲r̲e̲n̲{A B} < {C D} \…
那么三者周长L与面积S存在如下关系：
$$\begin{gathered} L_{\text{内接正多边形}}<L_{\text{圆}}<L_{\text{外切正多边形}} \\ {S}_{\text{内接正多边形}}<S_{\text{圆}}<S_{\text{外切正多边形}} \end{gathered}$$
于是，便可通过上式对圆周率进行逼近，数值逼近过程不再赘述，感兴趣的读者可自行推导并验证。

2 无穷级数

2.1 拉马努金(Ramanujan)圆周率公式

这里只举一个精度较高的计算公式，Ramanujan于1914年发现，其所发现的公式远不止这一个。该公式证明会涉及到椭圆积分、广义超几何函数、模函数等，笔者能力有限，就不详细推导了(手动狗头)。
$$\frac{1}{\pi }=\frac{2\sqrt{2}}{99^2}\sum _{k=0}^{\infty }\frac{(4k)!}{(k!{)}^4}\frac{26390k+1103}{396^{4k}}$$
取至第1项，有：$\pi =3.141592730013305660313996189025215518599\dots$；和真实值的相对偏差约为$\text{1e-8}$。

取至第2项，有：$\pi =3.141592653589793877998910597178909590809\dots$；和真实值的相对偏差约为$\text{1e-16}$。

取至第3项，有：$\pi =3.141592653589793238462653836575593749894\dots$；和真实值的相对偏差约为$\text{1e-24}$。

…

可见，每多一项，计算精度提升约8个数量级。

2.2 Chudnovsky圆周率公式

该公式由Chudnovsky兄弟于1988年发现，可认为是Ramanujan圆周率公式的变体，计算时每多一项，计算精度提升约14个数量级。
$$\frac{1}{\pi }=\frac{1}{53360\sqrt{640320}}\sum _{k=0}^{\infty }(-1{)}^k\frac{(6k)!}{(k!{)}^3(3k)!}\frac{13591409+545140134k}{640320^{3k}}$$

2.3 BBP公式

该公式由Baily，Borwein，Plouffe于1996发现，简称为BBP公式，且利用该公式证明了在十六进制下，不需要计算出$\pi$的小数点后前$n-1$位数，便可计算出$\pi$的小数点后第n位数，这也是该公式的创新之处。
$$\pi =\sum _{n=0}^{\infty }\left(\frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6}\right)\frac{1}{16^n}$$
下面给出该公式的一种证明。
$$\begin{array}{rl}I& =\sum _{n=0}^{\infty }\left(\frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6}\right)\frac{1}{16^n}\\ & =\sum _{n=0}^{\infty }\left({4\cdot \frac{1}{8n+1}{x}^{8n+1}|}_0^1-{2\cdot \frac{1}{8n+4}{x}^{8n+4}|}_0^1-{1\cdot \frac{1}{8n+5}{x}^{8n+5}|}_0^1-{1\cdot \frac{1}{8n+6}{x}^{8n+6}|}_0^1\right)\frac{1}{16^n}\\ & =\sum _{n=0}^{\infty }\left({\int }_0^{1}4x^{8n}\mathrm{d}x-{\int }_0^{1}2x^{8n+3}\mathrm{d}x-{\int }_0^1{x}^{8n+4}\mathrm{d}x-{\int }_0^1{x}^{8n+5}\mathrm{d}x\right)\frac{1}{16^n}\\ & ={\int }_0^1\sum _{n=0}^{\infty }{\left(\frac{x^8}{16}\right)}^n\left(4-2x^3-x^4-x^5\right)\mathrm{d}x\\ & ={\int }_0^1\frac{1}{1-\frac{x^8}{16}}\left(4-2x^3-x^4-x^5\right)\mathrm{d}x\\ & ={\int }_0^1\frac{16(1-x)}{\left(2-x^2\right)\left(x^2-2x+2\right)}\mathrm{d}x\\ & ={\int }_0^1\frac{4x}{x^2-2}+\frac{8-4x}{x^2-2x+2}\mathrm{d}x\\ & ={\int }_0^1\frac{4x}{x^2-2}+\frac{4-4x}{x^2-2x+2}+\frac{4}{x^2-2x+2}\mathrm{d}x\\ & =2\ln \left(2-x^2\right)-2\ln \left(x^2-2x+2\right)+{4\arctan (x-1)|}_0^1\\ & =\pi \end{array}$$

2.4 其他级数

下面将介绍一些比较常规的级数形式，并对部分作以证明。

(1) 自然数倒数偶次方和

下面仅给出次数为2、4、6的结果。
$$\begin{array}{rl}& \sum _{n=1}^{\infty }\frac{1}{n^2}=\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+\cdots =\frac{{\pi }^2}{6}\\ & \sum _{n=1}^{\infty }\frac{1}{n^4}=\frac{1}{1^4}+\frac{1}{2^4}+\frac{1}{3^4}+\cdots =\frac{{\pi }^4}{90}\\ & \sum _{n=1}^{\infty }\frac{1}{n^6}=\frac{1}{1^6}+\frac{1}{2^6}+\frac{1}{3^6}+\cdots =\frac{{\pi }^6}{945}\end{array}$$
因此，圆周率$\pi$可计算如下：
$$\pi =\sqrt{6\sum _{n=1}^{\infty }\frac{1}{n^2}}=\sqrt[4]{90\sum _{n=1}^{\infty }\frac{1}{n^4}}=\sqrt[6]{945\sum _{n=1}^{\infty }\frac{1}{n^6}}$$
下面给出第一个等式的一种证明，其余等式类似。

首先，写出函数$f(x)=x^2$在$[-\pi ,\pi ]$的Fourier级数，即：
$$f(x)=x^2=\frac{{\pi }^2}{3}+4\sum _{n=1}^{\infty }\frac{(-1{)}^n}{n^2}\cos (nx)$$
然后，将$x=\pi$带入可得：
$$f(\pi )={\pi }^2=\frac{{\pi }^2}{3}+4\sum _{n=1}^{\infty }\frac{(-1{)}^n}{n^2}\cos (n\pi )$$
由于$\cos (n\pi )=(-1{)}^n$，因此有：
$${\pi }^2=\frac{{\pi }^2}{3}+4\sum _{n=1}^{\infty }\frac{1}{n^2}$$
最后，求得：
$$\sum _{n=1}^{\infty }\frac{1}{n^2}=\frac{{\pi }^2}{6}$$
为验证其收敛效果，作出这三种求解方式随计算项数增加的结果，可见偶次方次数越高，收敛到\pi的速度越快。

自然数倒数偶次方和计算圆周率收敛效果

自然数倒数偶次方和计算圆周率相对偏差

(2) 泰勒展开

常见反三角函数的泰勒展开式如下所示：
$$\begin{array}{rl}& \arcsin x=\sum _{n=0}^{\infty }\frac{(2n)!}{4^n(n!{)}^2(2n+1)}{x}^{2n+1},|x|\le 1\\ & \arccos x=\frac{\pi }{2}-\sum _{n=0}^{\infty }\frac{(2n)!}{4^n(n!{)}^2(2n+1)}{x}^{2n+1},|x|\le 1\\ & \arctan x=\sum _{n=0}^{\infty }\frac{(-1{)}^n}{2n+1}{x}^{2n+1},|x|\le 1\end{array}$$
由于$\sin \frac{\pi }{2}=1$，$\cos \frac{\pi }{2}=0$，$\tan \frac{\pi }{4}=1$。因此，
$$\pi =2\arcsin 1=2\arccos 0=4\arctan 1$$
当然，通过该方法还可以得到更多的计算\pi的等式。如：
$$\begin{array}{rl}& \pi =4\left(\arctan \frac{1}{2}+\arctan \frac{1}{3}\right)\\ & \pi =2\left(\arcsin \sqrt{\frac{4}{5}}+\arctan \frac{1}{2}\right)\\ & \pi =4\left(4\arctan \frac{1}{5}-\arctan \frac{1}{239}\right)\end{array}$$

(3) 无穷乘积

首先给出计算公式如下：
$$\frac{\pi }{2}=\underset{n\to \infty }{\lim }\frac{1}{2n+1}{\left[\frac{(2n)!!}{(2n-1)!!}\right]}^2$$
证明：首先，给出Wallis公式，即：
$${\int }_0^{\frac{\pi }{2}}{\sin }^{n}x\mathrm{d}x={\int }_0^{\frac{\pi }{2}}{\cos }^{n}x\mathrm{d}x=\{\begin{array}{ll}\frac{(n-1)!!}{n!!}\frac{\pi }{2},& n\text{ 为偶数 }\\ \frac{(n-1)!!}{n!!},& n\text{ 为奇数 }\end{array}$$
其中：
$$n!!=\{\begin{array}{ll}n(n-2)\cdots 4\cdot 2,& n\text{ 为偶数 }\\ n(n-2)\cdots 3\cdot 1,& n\text{ 为奇数 }\end{array}$$
由于当$x\in \left[0,\frac{\pi }{2}\right]$时，有$\sin x\in \left[0,1\right]$，因此：
$${\sin }^{2n+1}x\le {\sin }^{2n}x\le {\sin }^{2n-1}x$$
根据Wallis公式得：
$$\frac{(2n)!!}{(2n+1)!!}\le \frac{(2n-1)!!}{(2n)!!}\frac{\pi }{2}\le \frac{(2n-2)!!}{(2n-1)!!}$$
移项有：
$$1\le \frac{(2n+1)!!(2n-1)!!}{(2n)!!(2n)!!}\frac{\pi }{2}\le \frac{(2n+1)!!(2n-2)!!}{(2n)!!(2n-1)!!}$$
化简得：
$$1\le (2n+1){\left[\frac{(2n-1)!!}{(2n)!!}\right]}^2\frac{\pi }{2}\le \frac{2n+1}{2n}$$
作简单变形：
$$1\le \frac{\frac{\pi }{2}}{(2n+1){\left[\frac{(2n-1)!!}{(2n)!!}\right]}^2}\le \frac{2n+1}{2n}$$
对上式取极限得：
$$1\le \underset{n\to \infty }{\lim }\frac{\frac{\pi }{2}}{(2n+1){\left[\frac{(2n-1)!!}{(2n)!!}\right]}^2}\le \underset{n\to \infty }{\lim }\frac{2n+1}{2n}=1$$
因此，根据夹逼准则有：
$$\underset{n\to \infty }{\lim }\frac{\frac{\pi }{2}}{(2n+1){\left[\frac{(2n-1)!!}{(2n)!!}\right]}^2}=1$$
即：
$$\frac{\pi }{2}=\underset{n\to \infty }{\lim }\frac{1}{(2n+1)}{\left[\frac{(2n)!!}{(2n-1)!!}\right]}^2$$
Wallis公式证明如下：

记：
$$\begin{array}{rl}{I}_n& ={\int }_0^{\frac{\pi }{2}}{\sin }^{n}x\mathrm{d}x=-{\int }_0^{\frac{\pi }{2}}{\sin }^{n-1}x\mathrm{d}\cos x\\ & =-{{\sin }^{n-1}x\cos x|}_0^{\frac{\pi }{2}}+(n-1){\int }_0^{\frac{\pi }{2}}{\sin }^{n-2}x{\cos }^{2}x\mathrm{d}x\\ & =(n-1){\int }_0^{\frac{\pi }{2}}{\sin }^{n-2}x\left(1-{\sin }^{2}x\right)\mathrm{d}x\\ & =(n-1)I_{n-2}-(n-1)I_n\end{array}$$
移项得：
$$I_n=\frac{n-1}{n}{I}_{n-2}$$
因此可得：
$$\begin{array}{rl}{I}_n& =\frac{n-1}{n}{I}_{n-2}=\frac{n-1}{n}\frac{n-3}{n-2}{I}_{n-4}=\cdots \\ & =\{\begin{array}{l}\frac{n-1}{n}\frac{n-3}{n-2}\cdots \frac{3}{4}\frac{1}{2}{I}_0,n\text{ 为偶数 }\\ \frac{n-1}{n}\frac{n-3}{n-2}\cdots \frac{4}{5}\frac{2}{3}{I}_1,n\text{ 为奇数 }\end{array}\\ & =\{\begin{array}{ll}\frac{(n-1)!!}{n!!}{I}_0,& n\text{ 为偶数 }\\ \frac{(n-1)!!}{n!!}{I}_1,& n\text{ 为奇数 }\end{array}\end{array}$$
其中，
$$\begin{array}{rl}& I_0={\int }_0^{\frac{\pi }{2}}{\sin }^{0}x\mathrm{d}x={\int }_0^{\frac{\pi }{2}}1\mathrm{d}x=\frac{\pi }{2}\\ & I_1={\int }_0^{\frac{\pi }{2}}{\sin }^{1}x\mathrm{d}x={\int }_0^{\frac{\pi }{2}}\sin x\mathrm{d}x=1\end{array}$$
那么：
$$I_n=\{\begin{array}{ll}\frac{(n-1)!!}{n!!}{I}_0=\frac{(n-1)!!}{n!!}\frac{\pi }{2},& n\text{ 为偶数 }\\ \frac{(n-1)!!}{n!!}{I}_1=\frac{(n-1)!!}{n!!},& n\text{ 为奇数 }\end{array}$$
另一方面，
$$\begin{array}{rl}{I}_n& ={\int }_0^{\frac{\pi }{2}}{\cos }^{n}x\mathrm{d}x={\int }_0^{\frac{\pi }{2}}{\sin }^n\left(x-\frac{\pi }{2}\right)\mathrm{d}x=(-1{)}^n{\int }_{-\frac{\pi }{2}}^0{\sin }^{n}t\mathrm{d}\left(t+\frac{\pi }{2}\right)\\ & =(-1{)}^n{\int }_{\frac{\pi }{2}}^0{\sin }^n(-m)\mathrm{d}(-m)=(-1{)}^n{\int }_{\frac{\pi }{2}}^0(-\sin (m){)}^n\mathrm{d}(-m)\\ & =(-1{)}^n(-1{)}^n{\int }_{\frac{\pi }{2}}^0{\sin }^{n}m\mathrm{d}(-m)=-{\int }_{\frac{\pi }{2}}^0{\sin }^{n}m\mathrm{d}m\\ & ={\int }_0^{\frac{\pi }{2}}{\sin }^{n}x\mathrm{d}x\end{array}$$
至此完成了对Wallis公式的证明。

关于无穷乘积还有其他公式，如：
$$\frac{2}{\pi }=\underset{n\to \infty }{\lim }\cos \frac{\pi }{4}\cos \frac{\pi }{8}\cdots \cos \frac{\pi }{2^n}$$
在此不做详细展开了。

3 微积分

有如下定积分成立：
$${\int }_{-1}^1\sqrt{1-x^2}\mathrm{d}x=\frac{\pi }{2}$$
根据其几何意义，该积分为半径为1的圆的面积的一半，因此，积分值为$\frac{\pi }{2}$，当然，也可直接求其原函数：
$${\int }_{-1}^1\sqrt{1-x^2}\mathrm{d}x={\frac{1}{2}\left(x\sqrt{1-x^2}+\arcsin x\right)|}_{-1}^1=\frac{\pi }{2}$$
于是有：
$$\pi =2{\int }_{-1}^1\sqrt{1-x^2}\mathrm{d}x$$
根据微积分定义：
$${\int }_a^{b}f(x)\mathrm{d}x=\underset{n\to \infty }{\lim }\sum _{i=1}^{n}f\left(a+\frac{b-a}{n}i\right)\frac{b-a}{n}$$
于是，$\pi$可通过下式进行计算：
$$\begin{array}{rl}\pi & =2{\int }_{-1}^1\sqrt{1-x^2}\mathrm{d}x\\ & =2\underset{n\to \infty }{\lim }\sum _{i=1}^n\sqrt{1-{\left(-1+\frac{2}{n}i\right)}^2}\frac{2}{n}\\ & =8\underset{n\to \infty }{\lim }\sum _{i=1}^n\sqrt{\frac{i}{n}\left(1-\frac{i}{n}\right)\frac{1}{n}}\end{array}$$
数值计算结果及相对偏差如下图：

定积分计算圆周率

定积分计算圆周率相对偏差
可见，当$n=2^{19}$时，该公式计算出的圆周率与真实值相对偏差约为$\text{1e-9}$量级。当然，构造其他积分也可进行计算，比如常见的正态分布概率密度函数。

4 概率学(Monte Carlo)

这里提供一种直观易懂的几何概型，即：如下图，

边长为1的正方形及其内切圆
在边长为1的正方形内有一内切圆，则内切圆的面积可计算为：
$$S_{\mathrm{O}}=\frac{\pi }{4}$$
由几何概型，内切圆面积还可近似计算如下：在正方形内随机均匀的生成N个点，统计落在圆内的点M，则内切圆的面积近似为：
$$S_{\mathrm{O}}=\frac{M}{N}$$
综上，圆周率可表示为：
$$\pi =\frac{4M}{N}$$
对总点数$N$从$10^1$到$10^9$进行试验，所计算出的圆周率近似值如下图：

几何概型计算圆周率

几何概型计算圆周率相对偏差
可见，随着点数N的增加，$\pi$的近似计算值越来越接近真实值，当点数N达到10^9量级时，相对偏差约为$\text{1e-5}\量级。

几何概型计算$\pi$的著名实验还有Buffon投针试验，即：平面上有一组等距平行线，相邻平行线之间间隔为d，向该平面上随机投任意长为l的针，则针与平行线相交的概率为：
$$P=\frac{2l}{\pi d}$$
那么便可通过该式计算圆周率$\pi$，其证明过程中会用到微积分相关思想，在此不作赘述。

5 连分数

不难证明，有欧拉连分数公式：
$$\begin{array}{rl}& A=a_0+a_0{a}_1+a_0{a}_1{a}_2+\cdots +a_0{a}_1{a}_2\cdots a_n\\ & =\frac{a_0}{1-\frac{a_1}{1+a_1-\frac{a_2}{1+a_2+\frac{\ddots }{1+a_{n-1}-\frac{a_n}{1+a_n}}}}}\end{array}$$
因此，可由函数的泰勒展开得到其连分数表示形式，如：
$$\arctan x=\frac{x}{1+\frac{x^2}{3-x^2+\frac{(3x{)}^2}{5-3x^2+\frac{(5x{)}^2}{7-5x^2+\ddots }}}}=\frac{x}{1+\frac{x^2}{3+\frac{(2x{)}^2}{5+\frac{(3x{)}^2}{7+\ddots }}}}$$
由于$\arctan 1=\frac{\pi }{4}$，因此：
$$\frac{\pi }{4}=\arctan 1=\frac{1}{1+\frac{1^2}{2+\frac{3^2}{2+\frac{5^2}{2+\ddots }}}}=\frac{1}{1+\frac{1^2}{3+\frac{2^2}{5+\frac{3^2}{7+\ddots }}}}$$
那么，
$$\pi =\frac{4}{1+\frac{1^2}{2+\frac{3^2}{2+\frac{5^2}{2+\ddots }}}}=\frac{4}{1+\frac{1^2}{3+\frac{2^2}{5+\frac{3^2}{7+\ddots }}}}$$
经过编程计算，两种连分式表示的$\pi$计算数值及相对偏差如下。其中，右侧连分式在计算至$n=21$时和$\pi$的真实值相对偏差已经达到$\text{1e-16}$量级，精度非常高。有关$\pi$的连分式表示还有很多，感兴趣的读者可自行了解。

6 总结

本文通过割圆术、无穷级数(泰勒展开)、微积分、概率学、连分数等多种不同的思路对圆周率\pi进行了计算，均得到了不错的计算精度。对圆周率的计算，一方面可以提升算法的计算效率，验证计算机的计算性能，另一方面也可以促进某些学科的发展。还有其他方法也可以计算出\pi，等有时间了，再加进来。

最后加一句：探索永无止境！