多种方法计算Pi并且精确度比较

多种方法计算圆周率并比较精确度

【摘要】本文介绍了多种方法求圆周率的近似值并对各种方法进行精确度的比较得出具体情况选择的方法，且通过mathematica 编程模拟实验过程，得出各种方法的特点。 【关键字】圆周率 数值积分法 泰勒级数法 蒙特卡罗法 拉马努金公式法

0.引言

平面上圆的周长与直径之比是一个常数，称为圆周率，记作π。在很长的一段时期，计算π的值是数学上的一件重要的事情。有数学家甚至说：“历史上一个国家所得的圆周率的准确程度，可以作为衡量一个国家当时数学发展的一面旗帜。”足以见圆周率扮演的是角色是如此举足轻重。

π作为经常使用的数学常数，它的计算已经持续了2500多年了，到今天都依然在进行着，中间涌现出许多的计算方法，它们都各有千秋，在此，我们选择几种较典型的方法，包括数值积分法，泰勒级数法，蒙特卡罗法，韦达公式法，拉马努金公式法以及迭代法来和大家一起体验π的计算历程，同时利用mathematica 通过对各种方法作精确度的比较得出选择的优先顺序，为相关的理论研究提供一定的参考价值。

1.数值积分法

以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G 是一个扇形，由曲线y=

2

11

x

+(x ∈[0,1])及两条坐标轴围成，它的面积S=4π。算出了S 的近似值，它的4倍就是π的近似值。（

4111

2π

=+?dx x ） 用一组平行于y 轴的直线x=i x (1≤i ≤n-1，a=0x <1x <2x <...

地看作抛物线，就得到辛普森公式。

梯形公式：S ≈

]2

...[10121y y y y y n a

b n +++++-- 辛普森公式：S ≈

)]...(4)...(2)[(6212

3211210--+++++++++-n n n y y y y y y y y n a

b

Mathematica 程序如下：

n=1000;y[x_]:=4/(1+x^2);

s1=(Sum[y[k/n],{k,1,n-1}]+(y[0]+y[1])/2)/n;

s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/(6*n); Print[{N[s1,20],N[s2,20],N[Pi,30]}]

注：以上s1,s2分别是用梯形共识和辛普森公式计算出的π。最后一句中的N[s1,20]表示s1的前20为准确有效数字组成的近似值。N[Pi,30]是π的前30位有效数字组成的近似值。Print[]语句表示将方括号内的数显示出来。

再取n=5000,10000时分别算出π的近似，记录入表格如下：

梯形公式

辛普森公式

圆周率的准确值

n=1000 3.1415924869231265718 3.1415926535897932385 3.14159265358979323846264338328 n=5000 3.1415926469231265718 3.1415926535897932385 3.14159265358979323846264338328 n=10000

3.1415926519231265718

3.1415926535897932385

3.14159265358979323846264338328

注：下划线标记的数字为用对应方法所求得的精确到的位数。

观察发现，当n 的取值越大，即当[a,b]之间分的越细时，计算所得的π的精确度越高。n 从

1000到10000，梯形公式所求得的近似值精确度增加两位，而辛普森公式所求的的近似值有效位数没有变化，可见n 的变化对梯形公式影响较大。但总的而言利用辛普森公式所求得的近似值精度比用梯形公式所求得的近似值很多。当要求精确度较高时，两种方法中可优先选择辛普森公式法。

2.泰勒级数法

利用反正切函数的泰勒级数：Arctan x=x- (1)

2)1( (531)

2153+--+-+--k x x x k k ①X=1代入上式得到莱布尼兹级数：∑∞

=---=11

1

2)1(4k k k π

②令a=arctan 21,b=4π-a,则tan b=tan(4

π-a)=a a

tan 4tan 1tan 4tan

ππ

+-=

2

11121

1?

+-

=31 因此b=arctan 31,即4π-arctan 21=atctan 3

1

，从而得到

4π=arctan 21+arctan 31 即π=4(arctan 21+arctan 3

1）

③令a=arctan

51，tan a=51，则容易算出tan2a=125,tan4a=119

120

, Tan(4a-4π)=a a 4tan 4tan 14tan 4tan ππ+-=119

12011119120

+

-=2391

4a-4π=arctan 2391

4π=4a-arctan 2391=4arctan 51-arctan 2391 即π=16arctan 51-4arctan 239

1 利用mathematica 编写①②③程序如下：

T[x_,n_]:=Sum[(-1)^(k-1)x^(2k-1)/(2 k-1),{k,1,n}]

4N[T[1,1000],20]

N[4(T[1/2,1000]+T[1/3,1000]),50] N[16T[1/5,1000]-4 T[1/239,1000],50] N[Pi,50]

注：4N[T[1,1000],20]为用①的方法求解近似值，N[4(T[1/2,1000]+T[1/3,1000]),50]表示用②的方法求解近似值，N[16T[1/5,1000]-4 T[1/239,1000],50]表示用③的方法求解近似值。 输出结果为： ① 3.1405926538397929260 ② 3.1415926535897932385 ③

3.1415926535897932385 π

3.1415926535897932385

观察可见，①的方法精确度不高，当n=5000,10000,时用①方法所得结果为

{3.1413926535917932384,3.1414926535900432385}，即当n 的值增大较大幅度时近似值的精确度提高不大，有效数依然只有4位。而由上面表格中的结果可知当n 取1000时，②③两种方法所得的前20位有效数字完全与π的前20位数相同，对②③两种方法取150个有效数字时所得结果为： ② 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813 ③

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813 π

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813

观察可知，用②③两种方法所求近似值前150位有效数字均和π依然相同，可见用这两种方法求解π的近似值的精确度非常高！当计算要求精确度比较高时，可以选取用②③两种方法均可。

3.蒙特卡罗法

单位圆的

4

1

是一个扇形G,它是边长为1的单位正方形G1的一部分，如图，0.0

0.2

0.4

0.6

0.8

1.00.00.20.40.60.81.0

单位正方形G1的面积S1=1。只要能够求出扇形G 的面积S 在正方形G1的面积S1中所占的比例k ，即可以求得S ，乘与4即为π的值。

此处利用正方形中随机地投入很多点，使所投的每个点落在正方形中每个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m 与所投点的总数n 的比n

m

可以作为k 的近似值。在mathematica 中用Random[]函数产生[0,1]之间的随机数。 Mathematica 编写程序如下：

Module[{n=1000,p={},m,x,y,k},Do[m=0;Do[x=Random[];y=Random[];If[x^2+y^2≤1,m++],{k,1,n}];AppendTo[p,N[4m/n]],{t,1,10}]; Sum[p[[t]],{t,1,10}]/10]

由于程序中采用随机数思想，故而所得结果不唯一，运行三次程序所得结果为：3.158，3.1348，3.1504。

观察可知，用该方法所得的近似值精确度不是很高，只是得出了两位有效数字。但思想方法比较简单，比较直观。进一步投点5000,10000。

投点5000，三次运行结果为：3.1516，3.15056，3.14168 投点10000，三次运行结果为：3.15032，3.1412，3.13656

观察发现，随着投点数的增加，会增加结果的精确度，但是效果非常不明显。对比数值积分法和泰勒级数法，蒙特卡罗法的精确度确实比较低，但是当所求不是一个扇形面积而是其它比较特殊的图形时（如100个已知圆的公共部分G 的面积）时，利用蒙特卡罗法利用计算机实现是非常简单的事情。所以，蒙特卡罗法在很多场合，特别是在对精确度要求不太高的情况下是大有用武之地的。

4.拉马努金公式法

1914年，印度天才数学家拉马努金在他的论文里发表了一系列共14条圆周率的计算公式。

n

n n n n 404396)

26391103()!()!4(980221

+=∑∞=π 这个公式每计算一项可以得到8位的十进制精度。1985年Gosper 用这个公式计算到了圆周率的17,500,000位。

Mathematica 编写程序如下：

a[n_]:=Sum[(4k)!*(1103+26390k)/((k!)^4*369^(4k)),{k,0,n}] N[1/(2Sqrt[2]/9801*a[1000]),15] N[Pi,15]

输出结果为：3.14159262864498 ，3.14159265358979 可见用该公式可以精确到8位有效数字。

1989年，大卫·丘德诺夫斯基和格雷高里·丘德诺夫斯基兄弟将拉马努金公式改良，这个公式被称为丘德诺夫斯基公式，每计算一项可以得到15位的十进制精度。1994年丘德诺夫斯基兄弟利用这个公式计算到了4,044,000,000位。以下是丘德诺夫斯基公式的一个便于计算机计算的版本。

Mathematica 编写程序如下：

b[n_]:=Sum[((6k)!(13591409+545140134k))/((3k)!(k!)^3(-640320)^(3k)),{k,0,n}] N[426880*Sqrt[10005]/b[1000],150]

N[Pi,150]

输出结果：

近似值

3.1415926535897932384626433832795028841971693993751058209749445923078164

0628620899862803482534211706798214808651328230664709384460955058223172535940813

准确值

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813

由表格可知，当n 取1000时，该公式计算所得的近似值前150位与圆周率的精确值前150位完全相同，可见该公式计算所得的精确度已经非常高，同时这个公式也非常适合计算机编程，

是目前计算机使用较快的一个公式。

5.算法比较

由以上计算圆周率的方法和计算结果，我们可以得知，当要求计算的精确度非常高时，我们可以选择泰勒级数法中的公式

π=4(arctan 2

1+arctan 3

1）以及

π=16arctan 51-4arctan

239

1 和丘德诺夫斯基公式;当精确度要求较高时我们可以选择数值

积分法中的梯形公式，辛普森公式以及拉马努金公式n

n n n n 404396

)

26391103()!()!4(980221

+=∑∞=π;当精确度要求不高时，我们可以选择泰勒级数法中的莱布尼兹级数以及蒙特卡罗法，特别是当所求不是一个扇形面积而是其它比较特殊的图形时（如100个已知圆的公共部分G 的面积）时，可利用蒙特卡罗法快速简单的求出近似值。

6.结束语

随着计算机技术的飞速发展和计算方法的突破与创新, 计算的世界纪录正在迅速地被

刷新。1949 年计算机 ENIAC 用 70 分钟的时间把 P 计算到 2037位以来, 到 1999 年东京大学的 Kana-da 和Takahashi 用两台超级计算机花费了 73 小时把 P 计算到 68, 719, 470, 000 位, 这期间由于计算机运算速度的加快, 计算位数越来越多, 相对使用时间越来越短。 虽然这样高的精确度已经没有太多的实际意义,但这反映了现代数学科学的日新月异,反映了人类智慧向极限的挑战。就像东京大学的 Kanada 所说的，“我们进行

π运算的主要原因

是检验我们计算机的可靠性和对于运算、程序和算法的正确性,成为世界记录的保持者只是它的一个副产品而己”。

参考文献

[1]张晓贵.圆周率计算的四个时期［J ］.辽宁教育学院学报，2009，9(17)，32 －33. [2]何光.用蒙特卡罗方法计算圆周率的近似值.内江师范学院学报,第23卷第4期. [3]陶珊珊，孟凤娟.计算π的几种方法. 科技信息，2010年第35期

[4]强春晨,刘兴祥,岳育英.圆周率计算方法发展史.延安大学学报(自然科学版), 2012年6

月第31卷第2期

[5]王广超.圆周率的数学实验算法设计.博士·专家论坛

[6]李尚志，陈发来，张韵华，吴耀华著.数学实验第二版

[7]王维平.实验数学导引:数学原来可以这样做（2004）

Several ways to calculate pi and compare the accuracy

LUO Yan

(School of Mathematics and Applied Mathematics, Sun Yat-sen

university,GuangZhou,Guangdong province)

Abstract:This article describes several ways to seek an approximation of pi and

tells us when to use these ways bying comparing the accuracy of the approximation.Besides,by using mathematica program to make simulation experiments , we can get the features of these ways.

Key words:Pi,Numerical integration method,Taylor series method,Monte Carlo method,Ramanujan formula method