pi计算实验报告（18页）

实验报竹

班级：02

姓名：张海洋

学号：9

圆周率（选做）

要求：先查资料看看古人是怎样计算 n的，再对n的各种计算 方法进行研究和讨论（收敛速度等），并给出不同算法算出 n

的小数点后第10000位的数字是什么，你觉得该数字应该是多少

第一部分:

圆周率简介

圆周率是指平面上圆的周长与直径之比 （ratio of the circumfere nee of a circle to the

diameter）。用符号n （读音：p a）i表示。中国古代有圆率、圆率、周等名称。它是一 个常数（约等于）。它是一个无理数，即无限不循环小数。在日常生活中，通常都用代 表圆周率去进行近似计算。而用十位小数便足以应付一般计算。即使是工程师或物理学 家要进行较精密的计算，充其量也只需取值至小数点后几百个位。计算圆周率的方法

“历史上一个国家所算得的圆周率的准确程度， 可以作为衡量这个国家当时数学发展

水平的指标。”

历史上最马拉松式的计算，其一是德国的 Ludolph Van Ceulen他几乎耗尽了一生的

时间，计算到圆的内接正262边形，于1609年得到了圆周率的35位精度值，以至于圆 周率在德国被称为Ludolph数；其二是英国的 William Shanks,他耗费了 15年的光阴，在 1874年算出了圆周率的小数点后707位。可惜，后人发现，他从第528位开始就算错了。

　把圆周率的数值算得这么精确，实际意义并不大。现代科技领域使用的圆周率值，有十 几位已经足够了。如果用Ludolph Van Ceulen算出的35位精度的圆周率值，来计算一个 能把太阳系包起来的一个圆的周长，误差还不到质子直径的百万分之一。以前的人计算 圆周率，是要探究圆周率是否循环小数。自从 1761年Lambert证明了圆周率是无理数，

1882年Lin dema nn证明了圆周率是超越数后，圆周率的神秘面纱就被揭开了。

在中国，公元263年前后，刘徽提出著名的 “割圆术”求出了比较精确的圆周率。

他发现：当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而 多边形的面积也会越来越逼近圆面积。于是，刘徽利用正多边形面积和圆面积之间的关 系，从正六边形开始，逐步把边数加倍：正十二边形、正二十四边形，正四十八边形……， 一直到正三。七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点 后第四位。在刘徽研究的基础上，祖冲之进一步地发展，经过既漫长又烦琐的计算，一 直算到圆内接正24576边形，而得到一个结论： V n V 同时得到n 的两个近似分

2

2、

数：约率为22/7;密率为355/113。他算出的n的8位可靠数字，不但在当时是最 精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为

“祖率”。 现在的人计算圆周率，多数是为了验证计算机的计算能力，还

有，就是为了兴趣。

第二部分：

古人计算圆周率方法

古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆 的周长。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072 边形得到5位精度；Ludolph Van Ceulen用正262边形得到了 35位精度。这种基 于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进 行数学研究时有意无意地发现了许多计算圆周率的公式。 下面挑选一些经典的常

116 arctg 亏14arg tg 239

1

16 arctg 亏

1

4arg tg 239

Machin公式

2 n 1

3 5 7

arctg(X)X； X5 X7 L

(1) n 1 x

2n

这个公式由英国天文学教授 John Machin于1706年发现。他利用这个公式计算到了 100位的圆周率。

Machin公式每计算一项可以得到位的十进制精度。 因为它的计算过程中被乘数和被除数都不大于长整

数，所以可以很容易地在计算机上编程实现。

还有很多类似于 Machin公式的反正切公式。在所有这些公式中， Machin公式似乎是最快的了。虽然

如此，如果要计算更多的位数，比如几千万位， Machin公式就力不从心了。下面介绍的算法，在 PC

机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因 为计算过程中涉及两个大数的乘除运算，要用 FFT(Fast Fourier Transform)算法。FFT可以将两个大数的

乘除运算时间由 0(n2)缩短为0(nlog(n))。

Ramanujan 公式

Ramanujan 公式

9801

1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共 14条圆周率的计算公式， 这是

其中之一。这个公式每计算一项可以得到 8位的十进制精度。1985年Gosper用这个公式计算到了圆周

率的 17,500,000 位。

3、

AGM(Arithmetic-Geometric Mea n)算法

Gauss-Legendre 公式：

初值：

a+6重复计算：

a+6

STWHHE HFW

这个公式每迭代一次将得到双倍的十进制精度，比如要计算 100万位，迭代20次就够了'

1999年9月Takahash和Kanada用这个算法计算到了圆周率的 206,158,430,000位，创出 新的世界纪录。

4、

Borwe in

四次迭代式：

初值：

'- : 二_ ： 1 重复计算:

sr泸 1

飞?厂1+（卜铲）"4

(1}' 1 -

最后计算： ■!，

这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。

5、

Bailey-Borwein-Plouffe 算法

丄（丄o

丄（丄

o16n 8n 1

2 1 1

8n 4 8n 5 8n 6

这个公式简称 BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995年共同发 表。它打破了传统的圆周率的算法，可以计算圆周率的任意第 n位，而不用计算前面的n-

1位。这为圆周率的分布式计算提供了可行性。 1997年，Fabrice Bellard找到了一个比BBP

32 G 二f快

32 G 二f

滋麗一一璽—八 葆一 或I 10?+1 ~ ib/i+S " 10n+o 10n+7 10m+9

第三部分：

对于n的几种计算的研究和讨论:

1、

Comma nd window

? n=lCOOO;

y=sq.rt (1 一貨."2);

h=l/n;

a(ri3=vps r aps (y.) ,11)

ans =

3,5

n=10

ans

n=20

ans

n=50

ans

n=100

ans

n=200

ans

n=500

ans

n=1000

ans

n=2000

ans

半径为1的圆称为单位圆，它的面积等于 。只要计算出单位圆的面积，就算出了

在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第 一象限的部分是一个扇形，而且面积是单位圆的 1/4,于是，我们只要算出此扇形的面

积，便可以计算出 。但计算量较大。

2、

数值积分法（II）利用公式

1 1

4 dx

01 x2

Command Window

Command Window

X

X

Command Window

Command Window

X

X

? n=50:

i=a； i/r； i； s=0;

Tou k=i; 1 enetli ti)-1 s=s+(l/(l-K((i(i)+i(k-kl)V2) ?nd

arts =

3, 14162536692300

n=10

ans

n=20

ans

n=50

ans

n=100

ans

n=200

ans

n=500

ans

n=1000

ans

n=2000 ans =;

设分点Xi,X2「Xn-1将积分区间[0,1]分成n等分。

所有的曲边梯形的宽度都是h=1/n。记yi=f(xi).则第i个曲边梯形的面积A近似地等于梯形 面积，即：

A=(y(i-1)+yi)h/2

将所有这些梯形面积加起来就得到：

A~ 2n[ 2(yi+y2+…yn-i)+y0+yn]

3、

利用复化梯形算法求Pi的近似值.

1 40 1 X

1 4

0 1 X2

=1/2 (Tn+h

n

h \

f a

i h o

i 1

2丿

Com Window

yy clear；a-0

* 1= (b-aj /2* (￡ (a) +f (bD ;

cr=l;n=l ;

^uhile er>l. Oe-6 h^(b-a)/n.

5=0；

for 1 s=s-H(a+i*h-h/2> ;

end t2=(tHh*s)/?: eu=abs (t2~t 1);

fprin+f t=K. 6f j r=9t. 6￡\n a 12S er); n=2tn; tl=t2,

end

end

t-3. IOCOOO, r=0.100000

t=5. 131170J r=C< C31]?G

t=3.. 138988, r=0. 007912

：i=3.140942, r=0. 001955

t=3, 141430, r=0. 00048&

t=3. 141 啤 r=0. 00C122

t=3.141B82, r=0.000031

t=3. 141590, r=C. 000008

t=3.141592, r=0.000002

t=3. 141592, r=0. OCOOOO

?

4、

泰勒级数法计算 ：

3 5

2k 1

利用反正切函数的泰勒级数arcta nx x x x

(1)k1 x

来计算

3 5

2k 1

(l)x=1 时一

1丄

1 丄

4

3

5 7

n i k 1

—2k 1

=4*

k 1

Command Window

Command Window

Command Window

Command Window

? n=lQOO;

? 5=0； for k=1：n

e=s44? 1) * frlrl) / (2*k-jj);

End

ypa(Ej 20)

ans -

J. 1405926538397941350

n = 10 ans=;

n =

=20

ans

n =

=50

ans

n =

=100

ans

n =

=200

ans

n =

=500

ans

n =

=1000

ans

n =

=2000

ans

x=1时得到的arctan1的展开式收敛太慢，逼近速度太慢，运算庞大，对速度造成了很大 影响。

5、

泰勒级数法计算 ：

3 x

5 x

(1)k1

2k 1 x

3

5

2k 1

arcta n

1

arcta n -

2

3

来计算利用反正切函数的泰勒级数

来计算

arcta n x x

(II)进一步精细化 arcta n1

? eyms n;

* (n-1) *(1/2) (2*n.-1V (2*n-l);

f2=C-lV ;

ans 1 - symsujn (f 1; 1, 79),

anssynsujnCf2, g I, 79);

ari￡=vp a (4# (^51+^52), 20)

ans -

3. H15926535897932385

n = 10 ans=; 当x=1时得到的arctan1的展开式收敛太慢。要使泰勒级数收敛得快，容易想到， 1应当使x

当x=1时得到的arctan1的展开式收敛太慢。要使泰勒级数收敛得快，容易想到，

1

应当使x的绝对值小于1，最好是远比1小。例如，因为arctan1 arctan- arctan-，

3

所以我们可以计算出arctan丄,arctan1的值，从而得到arctan1的值。这样，就使得收敛

2 3

速度加快，逼近的速度大大增加?

6、

11、 1 1

利用麦琴给出一 4arctan— arctan ，推出 n =4(4arctan— arctan )

5 239 5 239

n :

=20

ans

n =

:50

ans

n :

=100

ans

n :

=200

ans

n :

=500

ans

n :

=1000

ans

n =

:2000

ans

Command Window

、\ syiris FL；

tl- (-1) * Cn-1)4 (1/5) ;

f2=(-l) Cn-D* (1/233) * (2*n-l)/C2*n-n : ansl=s7nsujn (f lf n, 1, TS);

(f2,n, 1^ 79);

ans=vpa(4* (■4*angl-ans2)J 2Q)

=

3-141592663&857932386

n =

10 ans=;

n = 20

ans =;

n = 50

ans =;

n = 100

ans =;

n = 200

ans =;

n = 500

ans =;

n = 1000

ans =;

n = 2000

ans =;

对泰勒级数，

随着I

xI的减小，级数的收敛速度明显加快，这启示我们另

外构造相关级数来逼近n。

7、

蒙特卡罗法计算

单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积

S 1。只要能够求出扇形的面积S在正方形的面积中所占的比例k S/S，就能立即 得到S，从而得到 的值。下面的问题归结为如何求k的值，这就用到了一种利用随机 数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投 的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。降落 在扇形内的点的个数m与所投店的总数n的比可以近似的作为k的近似值。

1

1、

| Command Window

? k=2000,

? JlirD；

for rt= 1 ik

if 匚andU) " 2+rand (1) '2^.= 1 m=ni+l;

enA

end

vpa(4*VlJ)

ans =

3. 108DOQOOOOOOODDOOOOOOOOOODOOOOO

n = 100 ans=;

n

=200

ans =

n = 500

ans

n

=1000

ans :

n =

2000

ans =;

n =

5000

ans =;

n =

:10000

ans =

n =

20000

ans =;

这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为, 虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情 况下很有用的。

除以上几种方法之，还有下列一些的求其他 的方法:

*利用高斯公式

1 1 1

=48arctan — +32arcta n ——一20arcta n

、n

、n =2+1/3*(2+2/5*(2+3/7*(2+

(2+k/(2k+1)*(2+ )))……)))当….(

k=2799时可精确到800位）

、(6=1/2+1/2*1/(3*2八 3)+((1*3)/(2*4))*(1/(5*2八5))+

、eA( ni)+仁0 (欧拉公式，也称世界上最杰出的公式)

、1+(1/2)A2+(1/3)A2+(1/4)A2+ ….(1/rf)2/6= n

、1+(1/2)A4+(1/3)A4+(1/4)A4+ ….(1/rf)44=) n

、1+(1/2)A6+(1/3)A6+(1/4)A6+ ….(1/rf)6/945 n

、1+(1/2)八8+(1/3)八8+(1/4)八8+ ….(1/rf)8?450 n

、1+(1/2)A10+(1/3)A10+(1/4)A10+ ….(1/n)W=355n

第四部分：求的小数点后第位数字的几种方法:

第四部分：

求的小数点后第

位数字的几种方法:

反正切函泰勒级3arcta n x x32k 1(1*得到 14(1)kk 02k 1

反正切函

泰勒级

3

arcta n x x

3

2k 1

(1*

得到 1

4

(1)k

k 0

2k 1

推出：■ =4 *

1)k2k 1

Mathematica 程序

riB' s n = In finity;

5 = 4 ?

J 2 k + 丄

Z上￡ 4丄

Ort】。］-3> 14139265358979323&462C433a32795025841971693?9J\

75105520974^4459230751€4052&6205^962603402&343

小数点后10000位数字为8

2、

沃里斯(Wallis)方法

2 2 2

1 3

4

4

6 6

3

5

5 7

2 -2k

k 1 2k

2k

1

2k 1

Mathematica 程序

r~rn

S3

n s Infinity;

5 = 2 P[ (2 t / C2

t，1) t 2 1c / (2 J： + 1))；

N[s, 10 0021

3lI； ：＞ 3,14L552?535S9"532324626433S3279502Be4L971?93993\

75L05E2C 974 344S923C731€46622 620 2 9^3^250 3 4 3Z5S42--

回 S3

35yy5U545UTirawr2UgTTU53B7UlJT7^45tUUI^7T742BZr

9S922 5?9S9Eaai59205600101fi5525fi 375e7B?

小数点后10000位数字为8

3、

基于arctanx的级数

1

4 arcta n

5

麦琴给出：一

4

、 1 推出n =4( 4 arcta n— arcta n

5

arctan^^ ;

239 _1_

239 )

MATLAB程序

小数点后10000位数子为8

求取的方法很多，过程类似，不再累述

验证

lr- < B[Plr 10 002j

3.141 592€535B5_5323E462d433B32" 95023341371^353^3-.

?51C5i2CS74?4 45a23C7：l￡40^a^20：^5e￡25034525342-

11^0 67 92214SO￡€5132E23(j6^47Q938^460955C5:223172--

1 1*TdAn?*d1

?耒琴名丄'

■—

? 口

回

—T

Jb3S453C3^40704-oS12US*140^327QD1264560UL623'_422 = -

0210927fi457?310€57922&55249eB7275e4€L012￡4e2￡99\

9￡3225^9596?￡15^55256375€^36

*

结论：n的小数点后10000位数字为8

第五部分：

心得体会

对于n的值我们早已知晓，但是对于这个数值如何得到却不清楚。 通过这次

学习，我们知道原来计算n的值有如此多的方法，并且知道了级数在一些计 算中的作用。通过本次实验，我们知道了运用各种方法计算n的近似值，知

道了 matlab和mathematica中关于级数的运算， 学习是一个枯燥的过程， 有些东西是必须懂的。也就是一些基础知识吧， 包括基础的操作和相关的理

论知识。因此，有一本教材看一次，必定会有所收获。看了，练习，思考。

　MATLAB是一个实用性很强，操作相对容易，比较完善的工具软件。使用起 来比较方便，通过操作可以很快看到结果，能够清晰地感觉到成功与失败， 编程问题最头疼的不是编程序，而是调程序，所以在你的程序编完之后，一 定要进行验证其正确性，你要尽量多的设想你的问题的复杂性，当然，要一 步一步复杂，这样才能保证你的程序的适用性很强。以后虽然没有数学应用 软件设计的课程了，但是对于一个学习数学的学生来说，以后用到 MATLAB

和mathematica的地方还有很多，我不会因为不用学而不去学，希望能在数 学建模之前，把两个软件运用灵便， 为以后更强大的数学世界打下夯实的基 础。

1.

y y eyns ii；

f 1=(-1) * * tn-l)*(l/5)' (2*n-l)/(24n-l);

* Cn-l)?(l/239)^(2*n-l)/(2*n-l):

arx51=sjTns\m (f lj ns 1』ir^f),

ms2= syikEun (f2, n, 1』inf);

ans=vpa(4* (4*an^ 1 -ai\s2)} 10002)

STlE

3.1415926535597932384026433832705028&419?169399

1592055001016552553^56706-