周率

南京师范大学泰州学院

毕业论文（设计）

（一三届）

题目：漫谈圆周率值的计算

院（系、部）：数学科学与应用学院

专业：数学与应用数学

姓名：

学号

指导教师：

南京师范大学泰州学院教务处制

南京师范大学泰州学院本科毕业论文

1

摘要：圆周率



在数学中是一个非常重要的常数，受到广泛关注。古今中外一代代的数

学家为计算



献出了自己的智慧和劳动，人类对



值的认识过程，反映了数学和计算技

术发展情形的一个侧面。本文首先研究圆周率计算的发展历程，然后给出了利用计算机

求解圆周率的三种算法的基本原理。借助Mathematica软件编程给出了每一种算法的计

算结果和误差分析，利用Matlab软件对所得数据进行了分析和比较，对三种算法的优

缺点进行了讨论，最后阐述了从计算圆周率的过程中得到的启示。

关键词：圆周率；计算；近似值；数学实验

Abstract:een

paper,wefirstdiscussthedevelopment

processofthecalculatingof,andthenwegivethebasicprinciplesofthree

algorithmstocalculateathematicaandMatlab,we

givetheresultsofthethreealgorithmsanddiscusstheadvantagesand

y,wegetsomeinspirationfromtheaboverearch.

Keywords:；calculation；approximation；mathematicalexperiment

南京师范大学泰州学院本科毕业论文

2

目录

1绪论.........................................................3

1.1研究意义..........................................................3

1.2国内外研究现状....................................................4

1.3本文的研究方法和主要解决问题......................................4

2圆周率简介和圆周率计算的四个时期.............................5

2.1圆周率的简史及其重要性............................................5

2.2圆周率计算的四个时期..............................................6

2.2.1无算法记录时期..............................少先队员手抄报内容 .................6

2.2.2几何推算时期.................................................7

2.2.3解析计算时期.................................................9

2.2.4计算机运算时期...........................相生相克五行 ...................10

3借助计算机求解圆周率的方法..................................12

3.1数值积分法.......................................................12

3.1.1算法原理...................................................12

3.1.2计算结果及误差分析.........................................14

3.2泰勒级数法.......................................................17

3.2.1算法原理...................................................17

3.2.2计算结果及误差分析.........................................18

3.3蒙特卡洛法.......................................................21

3.3.1正方形内投点法..............................................21

3.3.2蒲丰投针法..................................................25

3.3.3随机整数互素法.............................................27

4从圆周率计算中得到的启示....................................30

谢辞.........................................................31

参考文献......................................................32

附录.........................................................33

南京师范大学泰州学院本科毕业论文

3

1绪论

我们知道，平面上圆的周长与直径之比是一个常数，称为圆周率，记作



。在日常

生活中，人们经常与



打交道，



的计算伴随着人类的进步而发展，许多数学家在其计

算上发费了巨大的精力。有些数学家甚至说：“历史上一个国家所算得的圆周率的准确

程度，可以作为衡量一个国家当时数学发展的一面旗帜[1]”。我国伟大的科学家祖冲之（公

元429—500年）在前人的基础上深入地研究了圆周率，经过长期坚持不懈的努力，求

出了当时世界上最好的近似值，他利用割圆术，求出了精确到小数点后七位数字的圆周

率，并且明确指出了圆周率的取值在3.1415926和3.1415927之间，这一举世瞩目的成

就在世界上领先了一千多年。在很长一个时期里，计算的



值是数学上一件重要的事情。

本课题将从圆周率的简史和重要性开始，了解人们为什么这么执着的计算着圆周率。接

着分析圆周率计算的四个时期，更清楚的知道人们为探索圆周率花费的心血，使我们当

代大学生受益匪浅，启发我们不仅要学习前人的数学思维方式，更要学习他们孜孜不倦，

开拓进取，为科学奋斗终身的精神。最后，利用我们所学的数学知识，运用数学实验方

法，结合积分，迭代，随机试验等，在计算机的帮助下，介绍了三种计算圆周率的快捷

方法：数值积分法，泰勒级数法和蒙特卡罗法。借助Mathematica软件编程给出了每一

种算法的计算结果和误差，利用Matlab软件对所得数据进行了分析，对三种算法的优

缺点进行了讨论，

1.1研究意义

作为数学上的一个重要常数,



不仅用于圆的计算,而且也在很多的公式中出现,

就我们现在的中学数学教材来说,数学中的初等几何、高中的立体几何、代数中的三角

函数、统计学等等都要用到



,它是我们最熟悉的无理数；在物理学科中也有很多的公

式要用到



,比如单摆周期T的公式、库仑研究的两个带电质点的相互作用力的公式中

也有它的身影,还有其他的科学分支中也要用到



,在科学史上有重要的地位。同时从

其发展史可以看出在计算圆周率的过程中用到了极限的概念、微积分的思想、概率统计

的理论,我们还要有实数的理论,除了这些以外,还要靠数学和科学技术的发展,它展

示了数学思想、方法的发展历程,在寻求圆周率的计算过程中也发现了很多的问题,从

而推动了数学的发展,促使人们不断为提高计算速度而寻找新的计算方法、改进计算的

手段。正如前文中提到的“历史上一个国家所算得的圆周率的准确程度,可以作为衡量

这个国家当时数学发展的一个标志。”因此,圆周率的发展历史从一个侧面反映了数学

的发展历史尤其是算法的发展史,而且还是计算机科学的发展史,代表了当时的计算

机科学的水平[2]。

几千年来作为数学家们的奋斗目标，古今中外一代代的数学家为此献出了自己的智

慧和劳动。在中国有刘徽、祖冲之等，在国外有阿基米德、卡西等。近代科学家如华罗

庚、严士健等在其数论论文中也对圆周率问题进行了探讨。古往今来，从未有哪一个数

学常数能向圆周率那样吸引众多的学者。圆周率在各个时期的文明中都像一颗闪耀的明

珠，它往往能够在一定程度上折射出该文明数学发展的水平[3]。为求得圆周率的值，人

南京师范大学泰州学院本科毕业论文

4

类走过了漫长而曲折的道路，它的历史是饶有趣味的。通过几何、微积分、概率等广泛

的范围和渠道发现



，这充分显示了数学方法的奇异美。也使我们充分认识到数学的奥

秘，促使我们从多角度思考解决问题，不断的创新以满足数学科学的发展。

1.2国内外研究现状

由于对圆周率的认识过程在一定程度上反映了数学和计算技术发展情形的一个侧

面，因此圆周率的相关问题一直受到古今中外众多学者的关注，如263年刘徽在注释《九

章算术》时求得了



的近似值，南北朝时代祖冲之进一步得出精确到小数点后7位的



值。15世纪初卡西在求得圆周率17位精确小数值，打破祖冲之保持近千年的记录。此

后，在数学家的探索下，到1948年英国的弗格森和美国的伦奇共同发表了



的808位

小数值，成为人工计算圆周率值的最高记录。随着计算机的发展，数学算法的不断创新，



的计算更加精确。进入新世纪以来，国内许多学者开始借助计算机来探讨圆周率的计

算[3-9]，他们设计了许多的算法，取得了很好的结果，例如利用Excel软件产生随机数

来模拟撒芝麻的试验来估计圆周率的值；利用Mathematica软件，借助随机数的思想，

利用蒙特卡洛法研究圆周率的近似值等。

1.3本文的研究方法和主要解决问题

本课题将首先通过查阅文献法研究圆周率的发展历程，首先了解一下什么是圆周率

和圆周率计算经历的四个时期：经验型获得时期、几何推算时期、解析计算时期和计算

机运算时期。通过这些历史漫谈古今中外圆周率的计算史。接着，借助借助Mathematica

和Matlab软件，利用数学实验的方法研究利用计算机求解圆周率



的算法，主要研究

三种算法：数值积分法、泰勒级数法和蒙特卡罗法。针对每一种算法，我们利用

Mathematica软件编程给出每一种算法的计算结果和误差，利用Matlab软件对所得数据

进行了绘图分析，并对三种算法的优缺点进行比较、讨论，最后运用归纳总结法，阐述

从计算圆周率的过程中得到的启示。

南京师范大学泰州学院本科毕业论文

5

2圆周率简介和圆周率计算的四个时期

2.1圆周率的简史及其重要性

所谓圆周率,通俗地说,就是圆的周长与直径之比,它是一个常数。这个数不仅

是无理数,而且是超越数[3]。在距今天4000年前的巴比伦王国它已被发现,当时认为

圆周率的值是3或

8

1

3

。大约2600年前,就提出“化圆为方”问题即“作与圆相等面积

的正方形”，此问题成为世界三大难题之一。公元前3世纪，古希腊数学、物理学家阿

基米德（Archimdes,B.C.287-B.C.212）提出“将圆的半径作为高，将圆周的长度作为底

边的三角形的面积就等于圆的面积”，通过这种方法得到圆面积的计算公式为2r。

作为圆周率的符号，目前全世界都在使用



。



的语源是希腊语的第一

个字母。计算圆周率



的方法虽然很多，但归纳起来主要有4种：割圆术（又称为古典

算法）、分析法、“沙-波法”、椭圆积分法。这里只重点列举几位很有影响的数学家的结

果。

利用割圆术阿基米德科学而准确地首次确定

7

22

71

223

,取两位实用值为3.14或

7

22

,在理论上指出了利用割圆术可以求得任意准确度的



值,第一次在科学中提出误差

估计及其精确度和如何确定的问题。中国南北朝时期南朝的科学家、数学家祖冲之

(429—500)算出的密率为

113

335

,精确度是6位小数,化为循环小数时实际上循环节达到

112位。祖冲之的密率不但当时最准,而且领先了世界1000多年。335/113便于记忆,将

最小的奇数1,3,5各重复一次后“平均”斩为两段,再让大的“住楼上”小的“住楼下”

即可。有趣的是,它的分子和分母,都可以用完全平方数简单地表示出

来:

22

222

87

1597

113

335





。更有趣的是7,8,9是连续的自然数,而且1587。

分析法算



主要是将



展开为无穷幂级数来求



值,在分析法算



的大军中有著名

的数学家莱布尼茨、牛顿、欧拉等。“沙—波法”即相关二次算法,其代表人物是欧仁沙

拉明(EugeneSalamin)和理查德波伦特(RichardBrent)。椭圆积分法是建立在椭圆积分

变换的理论上,其代表人物是印度传奇的数学家拉马努金(1887-1920)。

“人工”算



经历了三个世纪,最高记录是由美国数学家列维史密斯(LeviSmith)

和雷恩奇在1949年6月算出的1121位



值。自从1946年有了第一台电子计算机，数

学家们便开始了用各种公式借助电子计算机算



的历程。据统计到2002年12月,利用

计算机已算出12411亿位



值。

1995年,美国的达维德贝利和加拿大的彼德波尔文、西蒙普洛菲(SimonPlouffe)

发表了一个以三人的姓氏命名的算法,简称BBP公式[4]:

南京师范大学泰州学院本科毕业论文

6











































68

1

58

1

48

2

18

4

16

1

0

iiii

i

i



它打破了传统的算



法,可以计算



的16进制数字任意第n位而不用算前面的1n位。



与“无穷”关系密切,其中



的无穷表达式主要指无限连分式(数)、无穷乘积、

无穷级数、反正切式等。



与“化圆为方”、超越数、希尔伯特第7问题、近似计算、

逼近理论、空隙、转圈悖论、伯努利难题、欧拉公式、黄金分割、弧度制、曲线长度、

曲线图形面积、旋转体体积等有关。

对



的深入研究,扩充、发展了数系理论；对



计算的方法和思路可以引发新的概

念和思想；能否计算出位数越来越多的



值,已成了许多专家用来检验计算机可靠性、

精确性、运算速度及计算容量的有力方法、手段和衡量计算进展的指标；



的研究成果,

在一定程度上反映出一个国家的数学水平。仅从这里就可以看出,



在自然科学中有着

多么重要的地位和作用。

2.2圆周率计算的四个时期

古往今来,从没有哪一个数学常数能象圆周率那样吸引众多的学者。关于圆周率确

切值的计算，几千年来它都是数学家的奋斗目标之一，一代又一代的数学家为此献出了

自己的智慧和劳动，回顾圆周率计算的发展历程，大致可分为四个时期。

2.2.1无算法记录时期

该阶段的特点是,



的获得并没有理论上的根据,而是从实际经验中得到的,一般

说来,精确度是不高的。

古埃及和巴比仑的



属于经验性获得阶段。在古埃及所留下的两批草纸之一的莱

登草纸上有一个例子：“有一块9凯特(即直径为9)的圆形土地,其面积多大?今取其直径

的九分之一,即1,则余8,作8乘以8,得64,这个大小就是面积。”[3]由此可见,他们认为

圆的面积等于一个边长为此圆直径的九分之八的正方形面积,通过简单的推算,就可得

出圆周长与其直径之比是

81

256

,大约是3.1605。在巴比仑,他们把圆的面积取为圆周平方

的十二分之一,由此似乎可以看出,他们认为圆周是直径的三倍即



取3。但在给出正六

边形及外接圆周长之比时,实际上又用了

8

25

即3.125作为



的值。以上的时间大约是公

元前2000年左右。

在希腊,公元五世纪以后,希腊人为了解决“化圆为方”问题,对计算圆面积和圆周

率作出了很大努力,其中有数学家Hippocrates，Anaxagoras，Antiphon,精舞门罗志祥 Eudoxus等人,他

们的工作为阿基米德的进一步的研究开辟了道路,欧几里得的《几何原本》虽集古希腊

几何之大成,但在圆周率上却无贡献。在阿基米德之前,希腊人的圆周率仍然是经验性

的。

公元前416年,罗马取代希腊成为地中海一带的统治者。罗马人的数学是不值一提

的。ViruviusPollio是罗马著名的建筑师,他应该熟知阿基米德的工作,可他却使用

南京师范大学泰州学院本科毕业论文

7

125.3,其误差比

7

22

大了十几倍,显然罗马的圆周率是经验性的。

在中国,成书大约在一世纪的《周髀算经》上记述了周公和商高的问答,在商高曰“数

之法出于圆方”下,有赵爽(公元220年)注（“周三而径一”）。东汉科学家张衡提出

10,而在西汉缉为定本的中国古典数学名著《九章算术》中仍沿用周三径一之说,

其精度比不上古埃及和巴比仑,这种状况一直延续到公元三世纪的魏晋时期,因为数学

家刘徽的出现而得以改变。

在古印度,宗教活动中的庙宇和祭坛等的建筑设计,需要用到数学知识,在梵文经典

《测绳的法规》中对此作了总结,所包含的内容可以上溯到公元前五世纪或更早的年代,

其中使用



的值往往用复杂的式子[3]表示如:

0044.3

15

2

14

2

















0883.3

12

3

1

1

4

2



















显然,这些近似值比3强不了多少,并且也都是经验性。印度这种使用经验性



值的

年代一直延续到公元六世纪数学家阿耶波多。

2.2.2几何推算时期

从几何上推算出圆周率近似值的应该首推古希腊的

Archimedes(287-212BC),他得到的圆周率是大于223/71而

小于22/7,他是用96边的圆内接正多边形和圆外切正多边形

来进行推算的。下面简单地介绍Archimedes的推算。

半径为1的圆,分别内截和外切正132n边形,设它们周

长的一半分别为

n

a和

n

b，如图2-1是当2n时的情形。当n递

增时可以得到递增的数列

12

,,,,

n

aaa,和递减的数列

12

,,,,

n

bbb,此二数列有相同

的极限



,显然有,

K

Ka

n



tan

，

K

Kb

n



sin

和

K

Ka

n2

tan2

1







，

K

Kb

n2

sin2

1







这里132nK，由上式可以推导出下式[2]:

1

211





nnn

aba

（2-1）

2

11nnn

abb



（2-2）

南京师范大学泰州学院本科毕业论文

8

33

3

tan3

1





a，

2

33

3

sin3

1





b

由（2-1）得到

2

a，再由（2-2）得到

2

b，又通过（2-1）得到

3

a，通过（2-2）得到

3

b，

如此一直得到

6

a和

6

b。而

66

ab。

需要说明的是,Archimedes并不是用我们这里的代数和三角符号而是用纯几何的方

法推导的,并且也没有使用我们现在使用的小数表示(小数的正式使用是在十六、十七世

纪的事),所以他从

1

a,

1

b出发推导出

6

a，

6

b是极为烦琐的,计算量是惊人的。

在公元前150年左右,希腊的天文学家Ptolemy制作了一份弦表,以半径的

60

1

作为

长度单位,每一单位分为60分,每一分又分为50秒。他算出了圆心角1度所对弦长为1

单位2分50秒,于是圆内接360边形的周长与直径之比是(1360单位2分50秒)120

单位,即

2

830377

33.1416

6060120



这该是自Archimedes以来的巨大进步。

印度在公元500—1000年间,出现了四、五个有名的数学家,印度数学由此而出现了

繁荣的景象。对圆周率得出最好近似值的是阿耶波多,他所得到的近似值是3.1416,但直

到十二世纪前后印度数学家始终没有使用过该值。在他的《阿耶波多书》里,他是这样

说的:100加4,乘以8,再加62000,结果是直径为20000的圆周的近似值,这就导致了圆

周率为3.1416,由于书中没有一处地方提示过证明的方法,所以我们无从得知他是如何

得出该结果的,但从其准确性上看,他应该是通过推算得出的。出生于1119年的巴士卡

拉给圆周率提出了两个值,“当圆的直径乘以3927再除以1250时,商就接近于圆周,要

么乘以22再除以7,便是大致的圆周了。”因此,我们自然会认为巴士卡拉能够推算出圆

周率的值是:

392722

3.14163.1429

12507



下面看中国,刘徽是三世纪中国著名的数学家,他是用割圆术来求圆周率的。作圆内

接正n和2n边形,设

n

L为内接n边形的周长，

n

S为内接n边形的面积,S表示圆的面积,L

表示圆的周长,圆的半径为r。如下的式子[2]成立:



2

2,

nnn

LrSSLLn

从而有2SLr

他算到192边形时得到314.1024100314.2704。刘徽用14.3

50

157

表示圆周率，被

称为“徽率”。刘徽所建立的一般公式

22

2

nnnn

SSSSS可以把圆周率计算到任意

的精度，它比阿基米德用内接和外切双方逼近的方法更为简洁。在刘徽之后两百年，南

北朝人祖冲之应用刘徽的割圆术，在刘徽的基础上继续推算，球出了精确的七位有效数

字的圆周率的值：3.14159263.1415927。在《中国科学技术史》中,李约瑟博士指：

南京师范大学泰州学院本科毕业论文

9

“在这个时期,中国人不久赶上了希腊人,并且在公元五世纪祖冲之和他的儿子祖暅的

计算中又出现了跃进,从而使他们领先了一千年。”

祖冲之所得圆周率的精度保持了记录达一千年,直到十五世纪中亚数学家al-Kashi

和十六世纪法国数学家Viete才计算出更精确的值,前者到第十四位,后者到第九位。到

欧洲文艺复兴之前,圆周率的最好结果是公元1600年VanCeulen所得的第35位。

2.2.3解析计算时期

欧洲的文艺复兴带来了一个崭新的数学世界



,数学公式的出现使圆周率的计算进

入了一个新的阶段,最早的公式之一是数学家Willis所得的:

2133557

224466







二最著名的公式是：

111

1

4357





该公式有时被归功于Leibniz(1646-1716),但首先发现该公式的似乎应该是数学家

JamesGregory(1638-1675)。

在当时这些公式是令人惊奇的,原因是等式的右边完全是算术的,而



是从几何中来

的,在以前的计算中一直使用的是几何的方法。但从计算



的角度看,这些公式其实并没

有太大的价值。例如,在Gregory的公式中,要使



的准确度达到小数点后4位,要求误差

小于0.00005=1/20000,这将需要计算到级数的第10000项。但是Gregory给出的一个一

般的公式却在计算



时非常有用,该公式[2]是：

1111

1

6335337333

3













该结果右边的级数收敛是很快的,其第10项是

9

1

1933





它小于0.00005,所以只要计

算到级数第9项,就可以得到至少小数点后4位准确值,利用下面的公式也可以计算出



的值,只需要把

2

1

和

3

1

分别代如前面的Gregory的一般公式就可以了。

11

11

tantan

423













如果能够找到一个形如11

11

tantan

4ab













具有较大的数a

和b,这样利用一般公式时就能得到较快收敛的级数,在1706年Machin发现了一个这样

的公式：11

11

4tantan

45239















在这里插入符号



的由来。在1647年的时

候,Oughtred曾使用符号



d

来表示圆的直径与周长之比,在1697年,DavidGregory曾用

r



表示圆的周长与半径之比,首次用



表示现在意义的是数学家WilliamJones在1706年

使用的,大数学家Euler在1737年采用了该符号,接着它便很快成为了标准的符号。

南京师范大学泰州学院本科毕业论文

10

由于计算公式己经有了,除了需要耐心外,



的计算已经没有什么困难了,而追求更

高精度的计算其实没有什么太大的意义了。但还是有一些人花费了巨大的时间和精力用

于对



的计算。一个相当典型的例子就是英国人Shanks,他用二十多年的时间算出了小

数点后707位,其结果于1837年发表,遗憾的是1945年,Ferguson发现只有前572位是

正确的。

在Shanks的结果发表之后不久,中国著名的科学家、数学家李善兰在用尖锥术求圆

面积时得到了如下式子:

1111131

44

232452467















这也是中国在这一

时期圆周率计算的代表之一。

在Shanks己知道



是一个无理数。因为早在1761年数学家Lambert就己证明该结

果。在Shanks计算



之后不久,Lindemann证明了



是一个超越数,这就意味着



不可能

是任何整系数多项式方程的解,这也说明了所谓的“化圆为方”是不可能。

通过纸笔计算去追求圆周率的更高精度,既没意义,也需要太多的精力,并且所得精

度也是有限度的,表2-1给出了此时期几个有名的成果[5],Shanks的527的精度保持了很

长时间的世界记录,一直到现代计算机的运用。

表2-1解析法计算圆周率的成果

年代计算人国籍位数

英国16

英国72

英国100

y法国127

奥地利140

德国205

英国527

on英国620

美国808

2.2.4计算机运算时期

前面长达3000多年的时代,数学家们都是用手计算



值,有的甚至花毕生的精力来

做这些繁杂的计算工作,多少代数学家的努力,在1948年两个美国人将记录推至小数点

后808位,这是人工手算圆周率的新的记录。但是,随着1945年第一台电子计算机问世

后,



值的计算不断迈入新的阶段,记录不时被刷新,1949年,三位美国科学家利用计算

机将圆周率算至小数点后2037位,前后花去70个小时的上机时间（见表2-2）。随着数

学的发展,数学家们在进行数学研究时有意无意的发现了许多计算圆周率的公式,运用

数学分析和计算机技术使得



值越来越精确,2002年12月日本东京大学的金田康正教授

宣布,耗费601小时56分更新了圆周率计算位数的全球记录,最新的为12411亿位。

在使用计算机计算的时代,圆周率的计算公式和计算方法也不断更新,其中基于

1914年印度数学家拉玛奴扬(Ramanujan,S.1887-1920年)的文章“模方程和



的逼近”

中提出的当时计算



值最快的公式,建立了椭圆积分变换理论上的计算方法在日本由金

田教授和日立共同开发的名为分解有理数化法(DRM法)的计算方法。

南京师范大学泰州学院本科毕业论文

11

表2-2借助计算机计算圆周率部分结果

年代所用计算机计算者国籍计算位数

1949ENIAC美国2035

1955NORC美国3089

1957PEGASUS英国7480

1958IBM704法国10000

1959IBM704法国16167

1959IBM7090英国20000

1961IBM7090美国100265

1966IBM7030法国250000

1967CDC6600法国500000

1986CLAY-2美国29360000

1987SX-2日本133326000

1988HITACS-820/80日本201326000

事实上,我们在实际运算中往往只取



的前几位数就可以了,但是人们为什么仍然

对



的精确推算乐此不疲呢?德国的一位数学家曾经说过：“历史上一个国家所算得的圆

周率的准确程度,可以作为衡量这个国家当时数学发展的一个标志。”纵观



的计算发展

史,可见此话确有一番道理.在计算机技术高度发达的今天,计算



值又被认为对测试计

算机的性能具有科学价值,如上述提到的日立公司认为通过计算圆周率,可以进一步提

高编译器、数值计算和节点间通信的程序库、磁记录设备的输入输出性能调节以及长时

间高速稳定运行技术。

南京师范大学泰州学院本科毕业论文

12

3借助计算机求解圆周率的方法

在1000多年前祖冲之就把圆周率精确的计算到小数点后六位。作为一个现代的大

学生，我们所知道的数学知识比祖冲之所知多得多，并且有高性能的计算机作为辅助工

具,所以我们有理由相信自己一定能够找到求圆周率精度更高的方法。下面我们来介绍

一下三种计算圆周率的数学实验方法。

3.1数值积分法

3.1.1算法原理

半径为1的圆称为单位圆，它的面积等

于



。只要算出单位圆的面积，就算出了



。

以单位圆的圆心为原点建立直角坐标

系，则单位圆在第一象限内的部分G是一个

扇形，由曲线210,1yxx及两条坐

标轴围成，它的面积

4

S



。算出了S的近似值，

它的4倍就是



的近似值。

怎样计算扇形G的面积S的近似值？如

图3-1，作平行于y轴的直线将x轴上的区

间[0,1]（也就是扇形在x轴上的半径）分

成n等份（图3–1中21n）相应地将扇

形G分成n个同样宽度

n

1

的部分

k

G

（nk1）。每部分

k

G是一个曲边梯形：

它的左方、右方的边界是相互平行的直线段，类似于梯形的两底；上方边界为一段曲线，

因此称为曲边梯形。如果n很大，每个曲边梯形

k

G的上边界可以近似的看成直线段，从

而将

k

G近似的看成一个梯形来计算它的面积

k

S：梯形的高（也就是它的宽度）

n

h

1



，

两条底边的长分别是

2

1

1

1

















n

k

y

k

和

2

1















n

k

y

k

，于是这个梯形面积



n

yyT

kkk

1

2

1

1





可以作为曲边梯形面积

k

S的近似值。所有这些梯形面积的和T就可

以作为扇形面积S的近似值：

南京师范大学泰州学院本科毕业论文

13

11

0

121

4

1

2

knkn

n

n

SSSSTTT

yy

yyy

n

















n越大，计算出来的梯形面积之和T就越接近扇形面积S，而4T就越接近



的准确值。

图3-1中的扇形面积S实际上就是定积分

4

11

0

2



dxx

。与



有关的定积分很多，比

如

21

1

x

的定积分

4

1

11

0

2







dx

x

就比21x的定积分更容易计算，更适合于用来计算



。

一般的，要计算定积分dxxfS

b

a，也就是计算曲线xfy与直线

bxaxy,,0

所围成的曲边梯形G的面积S。为此，用一组平行于y轴的直线

bxxxxxanixx

nni



1210

,11将曲边梯形T分成n个小曲边梯形，

总面积S分成这些小曲边梯形的面积之和，如果取n很大，使每个小曲边梯形的宽度都

很小，可以将它上方的边界

ii

xxxxf

1

近似的看作直线段，将每个小曲边梯形近

似的当作梯形来求面积，就得到梯形公式。如果更准确些，将每个小曲边梯形的上边界

近似的看作抛物线段，就得到辛普森公式，具体公式如下：

梯形公式：设分点

11

,

n

xx将积分区间ba,分成n等份，即



n

abi

ax

i





，

ni0

。所有的曲边梯形的宽度都是

n

ab

h





。记

ii

xfy则第i个曲边梯形的面积

i

S

近似的等于梯形面积.

2

1

1

hyy

ii





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

0

121

()

2

n

n

yy

ba

Syyy

n







这就是梯形公式[3]。

辛普森公式：仍用分点



n

abi

ax

i





，

11ni

将区间ba,分成n等份，直线

11nixx

i

将曲边梯形分成n个小曲边梯形，在作每个小区间

ii

xx,

1

的中点

.

2

1

2

1n

ab

iax

i



















将第i个小曲边梯形的上边界

ii

xxxxfy

1

近似的看作经

过三点





















i

i

i

xxxxxfx,,

2

1,1

的抛物线段，则可求得

,4

6

2

11























i

i

ii

yyy

n

ab

S

南京师范大学泰州学院本科毕业论文

14

其中

.

2

1

2

1

















ii

xfy

于是得到

.42

6

2

1

2

3

2

11210









































n

nn

yyyyyyyy

n

ab

S

这就是辛普森公式[5]。

3.1.2计算结果及误差分析

我们利用Mathematica软件,用梯形公式和Simpson公式计算

1

2

0

4

1

dx

x





。可以

完成相应的程序编写。详细的计算程序见附录1，计算结果见表3-1

表3-1利用数值积分法计算圆周率结果及误差

1

梯形公式3.312643.312663.1275

误差-1.6666710-7-4.1666710-8-1.8518510-8

Simpson公式3.979053.978873.9798

误差-2.6645410-15-4.4408910-154.8849810-15

4

梯形公式3.312653.312633.01635

误差-1.0416710-8-6.6666710-9-4.6296310-9

Simpson公式3.980563.979053.97927

误差1.2434510-14-2.6645410-15-4.4408910-16

为方便观察，利用Matlab软件根据表3-1中的数据进行绘图（程序见附录3-2），

可得图像3-2至图3-5.通过观察发现，Simpson公式的计算结果要明显优于梯形公式，

观察图3-2和图3-3.可以看到利用梯形公式计算时，当迭代次数达到5000以上时，计

算结果的精度可达到7位数，而图像3-4和3-5显示，利用Simpson公式计算时，当迭

代次数达到5000时，计算结果的精度可达到11位数。

迭

代

次

数

算

法

南京师范大学泰州学院本科毕业论文

15

南京师范大学泰州学院本科毕业论文

16

南京师范大学泰州学院本科毕业论文

17

3.2泰勒级数法

3.2.1算法原理

利用反正切函数的泰勒级数











12

1

53

arctan

12

1

53

k

xxx

xx

k

k（3-1）

将1x代入(3-1)式得到：

.

12

1

1

5

1

3

1

11arctan

4

1





n

n



（3-2）

即1

111

4(1(1))

3521

n

n





。

n越大越精确。利用Mathematica软件编程（程

序见附录3-3）计算发现花费的时间很长，所得到的结果的准确度却很差。分析其原因

是由于当1x时得到的1arctan的展开式(3-2)收敛得太慢。

为使泰勒级数（3-1）收敛得快，应当使x的绝对值小于1，最好是远比1小，这样，

随着指数的增加，x的幂快速接近于0，泰勒级数就会快速收敛。比如，取

2

1

x

得到的

2

1

arctan

就收敛的快。例如：

12

1

3

2

1

12

1

1

2

1

3

1

2

1

2

1

tan

































n

n

n

ara

中取6312n得到的

2

1

arctan

的近似值

的误差就小于

652

1

，准确度已经非常非常高。

令

2

1

arctan，





4

，则

.

3

1

2

1

11

2

1

1

tan

4

tan1

tan

4

tan

4

tantan









































因此

,

3

1

arctan即

3

1

arctan

2

1

arctan

4





，从而得到

.

3

1

arctan

2

1

arctan

4





(3-3)

3

1

arctan比

2

1

arctan收敛的更快。利用泰勒级数计算出

2

1

arctan与

3

1

arctan的近似值再

相加，然后再乘以4，就得到



的近似值。

还可以考虑用

5

1

arctan来计算



。由

5

1

tan易算出

南京师范大学泰州学院本科毕业论文

18

.

239

1

arctan

4

4

.

239

1

119

120

1

1

119

120

4

tan4tan1

4

tan4tan

4

4tan

,

119

120

4tan,

12

5

2tan



















































从而得到

239

1

arctan4

5

1

arctan16（3-4）

利用xarctan的泰勒展开式求出

239

1

arctan,

5

1

arctan

的近似值，再代入公式（3-4）就

可以求出



的近似值。

3.2.2计算结果及误差分析

我们利用Mathematica软件,可以完成相应的程序编写（见附录3-3）。计算结果如

下表3-2和表3-3：

表3-2利用泰勒级数法计算圆周率结果及误差(4(arctan(1/2)+arctan(1/3))

展开项数123

计算结果3.3333333333.1172839513.145576132

误差0.191740680-0..003983478

展开项数456

计算结果3.1408505623.1417411973.141561588

误差-0..000148544-0.000031066

展开项数789

计算结果3.1415993413.1415911843.141592981

误差6.68710-6-1.46910-63.2810-7

表3-3利用泰勒级数法计算圆周率结果及误差(16arctan(1/5)-4arctan(1/239))

展开项数123

计算结果3.1712230223.1285568513.129580851

误差0.029630368-0.013035802-0.012011802

展开项数456

计算结果3.1295515943.1295525043.129552475

误差-0.012041059-0.012040149-0.012040179

展开项数789

计算结果3.1295524763.1295524763.129552476

误差-0.012040178-0.012040178-0.012040178

利用Matlab软件根据表3-2和3-3中的数据进行绘图（程序见附录3-4），可得图

南京师范大学泰州学院本科毕业论文

19

3-6、3-7、3-8和3-9。观察发现利用公式

.

3

1

arctan

2

1

arctan

4





计算圆周率的收敛速

度很快，并且当开展项数为8项时，计算精度可达7位数，精度较高。

南京师范大学泰州学院本科毕业论文

20

将两种展开方式所得实验结果进行比较（见图3-10），可以看到利用公式（3-3）

的收敛速度要远远高于公式（3-4）。

南京师范大学泰州学院本科毕业论文

21

3.3蒙特卡洛法

3.3.1正方形内投点法

3.3.1.1算法原理

利用单位圆与边长为1的正方形面积之比来计算



的近似值。具体思想如下:单位

圆的

4

1

是一个扇形

G

，它是边长为1的单位正方形

1

G的一部分，如图所示，单位正方

形

1

G面积.1

1

S只要能够求出扇形G的面积S在正方形

1

G的面积

1

S中多占的比例

1

S

S

k，就能立即得到S，从而得到



的值。

怎样求出扇形面积在正方形面积中所占的比例k？一个办法是在正方形中随机的

投入很多点，使所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点

落在扇形内。将落在扇形内的点的个数m与投点的总数n的比

n

m

可作为k的近似值[6]。

产生两个这样的随机数yx,，则以yx,为坐标的点就是单位正方形内的一点P，它

南京师范大学泰州学院本科毕业论文

22

落在正方形内每个位置的机会均等。P落在扇形内的充分必要条件是.122yx这样利

用随机数来解决问题的数学方法称为蒙特卡罗法。

3.3.1.2计算过程演示

通过Mathematica可以编写相应的程序（见附录3-7）来演示投点的过程（图3-11

—图3-13）。

南京师范大学泰州学院本科毕业论文

23

3.3.3.3计算结果及误差分析

通过观察,可以看到这种方法得到了一个



的近似值,虽然精确度不高,但是思想方

法简单,比较直观。观察发现,随着次数的增加会改善结果的精度，由于此方法采用随机

数的思想,故而结果不唯一，为尽可能的减小随机因素的影响，在利用程序计算圆周率

时，采用如下的流程计算：每投n个点得到一个



的近似值,将其存放在数组p中,同样

操作重复20次得到20个近似值,最后用Print语句显示全部近似值,并求出20个近似

值的平均值。这样可以减小随机因素的影响（程序见附录3-7）。计算结果见表3-4.

表3-4利用蒙特卡洛法计算圆周率计算结果级误差

投点数量100200300

模拟结果3.1100000003.1530000003.162666667

误差-0...021074013

投点数量400500600

模拟结果3.1265000003.1296000003.129333333

误差-0.015092654-0.011992654-0.012259320

投点数量700800900

模拟结果3.1231428573.1365000003.139555556

误差-0.018449796-0.005092654-0.002037098

投点数量1

模拟结果3.1422000003.1324000003.147733333

误差0.000607346-0.009192650.006140680

投点数量4

模拟结果3.1485000003.1473200003.143120000

误差0...001527346

投点数量200

模拟结果3.1411300003.1412466673.142050000

误差-0.000462654-0..000457346

利用Matlab软件根据表3-4中的数据进行绘图（程序见附录3-8），可得图3-14

和3-15。观察图像发现，随着投点数的增加，模拟的精度将会提高，即投点数的增加会

南京师范大学泰州学院本科毕业论文

24

改善结果的精度，当投点数达到10000次以上后，模拟的精度将达到2位数以上。

南京师范大学泰州学院本科毕业论文

25

3.3.2蒲丰投针法

3.3.2.1算法原理

(1)取白纸一张,在上面画许多间距为d的等距平行线;

(2)取一根长度为l(d1)的均匀直针,随机地向画有平行线的纸投去,共投n次(n

是一个很大的整数),观察针和直线相交的次数m;

(3)直线与针相交概率p的近似值可用

n

m

得到,进而可得



的近似值为

md

nl2

。

“针与直线相交”这一事件可以用图

3-16中的阴影部分表示。其中x表示针

落下后针的中点到最近一条平行线的距

离,表示针与平行线所成的夹角,阴影

部分的点满足:

.sin

2



l

x

.于是可得到

概率

d

l

p



2



，p的近似值可用

n

m

得到,

进而可以计算出



的近似值。

3.3.2.1计算结果及误差分析

通过Mathematica可以编写相应的程序（见附录3-9）,其中针长与间距分别取为2

和4,即有



1

p

，计算结果见表3-5.

表3-5利用蒲丰投针法计算圆周率计算结果级误差

投点数量100200300万的草书

模拟结果3..1062355553.200680525

误差-0.048057003-0..059087872

投点数量400500600

模拟结果3.2074640223.1684384723.184616969

误差0...043024316

投点数量700800900

模拟结果3.1164604973.1479673613.131756008

误差-0..006374708-0.009836646

投点数量1

模拟结果3.1822463743.1633047713.117298460

误差0..021712117-0.024294194

投点数量4

模拟结果3.1242825053.1242658893.126721964

误差-0.017310149-0.017326765-0.014870689

投点数量200

模拟结果3.1435165033.1352043433.138897469

误差0.001923849-0.006388311-0.002695185

利用Matlab软件根据表3-5中的数据进行绘图（程序见附录3-10），可得图3-17

和3-18.观察图像发现，随着投点数的增加，模拟的精度将会提高，即投点数的增加会

南京师范大学泰州学院本科毕业论文

26

改善结果的精度，但改善的结果并不明显，当投点数达到20000次以上后，模拟的精限购城市 度

将达到2位数以上。

南京师范大学泰州学院本科毕业论文

27

3.3.3随机整数互素法

3.3.3.1算法原理

利用随机整数互素的概率来得到



的近似值.具体过程如下。取一大整数N,在1到

N

之间随机地取一对整数

ba,

,找到它们的最大公约数ba,,做n次这样的实验,记录

1,ba的情况次数m,计算出

n

m

p

的值.理论分析,随机整数互素的概率为

,

6

21

1

222









p

因而

m

n6



。

3.3.3.2计算结果及误差分析

通过Mathematica可以编写相应的程序（见附录3-11）,计算结果见表3-6.

表3-6利用随机整数互素法计算圆周率计算结果级误差

投点数量100200300

模拟结果3.1227050763.1306873633.157066770

误差-0.018887578-0..015474116

投点数量400500600

模拟结果3.1329844413.1570345743.133481643

误差-0..015441920-0.008111011

投点数量700800900

模拟结果3.1622345343.1667099133.152577039

误差0...010984385

投点数量1

模拟结果3.1469378393.1416609333.146955890

误差0...005363236

投点数量4

模拟结果3.1486791223.1428631513.144321271

误差0...002728617

投点数量200

模拟结果3.1421845733.1414875433.142078386

误差0.000591920-0..000485732

利用Matlab软件根据表3-6中的数据进行绘图（程序见附录3-12），可得图3-19

和3-20.通过观察,发现随着次数的增加会提高结果的精度,同样不是很明显。观察投点

为10000时的结果:3.144321271。可以看到此时的精度比较稳定,基本达到2位以上的

有效数字。

南京师范大学泰州学院本科毕业论文

28

最后,通过表3-4，3-4和3-6及图3-14至图3-20，来观察三种计算



的近似值的

方法的精度情况。利用Matlab软件将三种方式的模拟结果放在一起进行比较，并将他

南京师范大学泰州学院本科毕业论文

29

们计算所得误差进行比较，可得图3-21和3-22.（程序见附录3-13）

南京师范大学泰州学院本科毕业论文

30

从图3-18和3-19可以看出,当投点次数不大,如1000和5000的时候,三种方法的

结果区别不明显;当次数足够大,如10000时,方法三(即利用随机整数互素的概率来得

到



的近似值)展示出一定的优势,特别是当随机整数范围取更大时,精度比其它两种方

法要好。

综上所述,三种计算



的近似值的方法,其特点在于用相关数学软件清晰明了地展

现出最后结果,通过编写程序能够模拟正真的实验过程,让我们不仅对蒙特卡罗方法的

思想--随机投点思想有了较清楚的认识,更对蒙特卡罗法的简洁实用的特点有进一步的

领悟。

4从圆周率计算中得到的启示

本文给出了圆周率



的简介和计算圆周率的四个时期，使人们更加熟悉了圆周率的

相关知识，并对计算圆周率的各种方法有初步了解。在此基础上，文章更加深入讲述了

运用数学分析和概率等的思想，结合数学实验的动手操作方法，利用计算机，采用现代

方法快捷的计算了



。



值的计算方法千变万化，各种计算方法锻炼着人们的思维，激

发想象力和创造力。通过漫谈圆周率的计算，让我们也学了一回祖冲之，探索了一下



的奥秘，同时感叹古人的智慧，对于当代人要站在巨人的肩膀上，利用所学，不断进取。

回顾了圆周率



从古至今的计算史后，不但认识到古今中外人类为计算它付出的辛

劳与智慧，还学习到了前人为追求科学而孜孜不倦的精神，并被



深深吸引。同时应用

所学知识，借助计算机，用数学实验方法试着多角度计算了圆周率的值，感受到了知识

就是力量，探索是无止境的。

通过这次对计算



的探索，增强了自己学习数学的兴趣和锻炼了动手能力。我认识

到在学习中可以把数学作为一门实验性的科学，从问题出发，借助计算机，亲自设计和

动手，体验解决问题的过程，从试验中去学习、探索和发现数学规律，这样很容易提升

学习数学的兴趣。学好数学很重要一点就是要是要自己动手去体验，追求内容的系统性、

完整性，激发动手探索的兴趣，培养动手操作能力。

南京师范大学泰州学院本科毕业论文

31

谢辞

论文得以完成，首先要感谢。。。，

同时要感谢四年来教导过我的各科老师，学院的各位领导，还有在我写论文过程中，

帮我一起搜集资料的朋友们。正是因为有你们，才使得这篇论文能完整的呈现在这里，

才能是自己完成了这个令人兴奋的任务。

任何一篇优秀的论文都离不开老师和朋友的参与、支持和帮助。而每一篇好的论文

又能为大家所分享和阅读，这真是一种善缘，愿我们在这样的关系中能成长和进步。

南京师范大学泰州学院本科毕业论文

32

参考文献

[1]王广超.圆周率的数学实验算法设计[J].科技信息,2009,(14):8-9.

[2]刘迪.圆周率



的发展简史[J].数学爱好者，2006，（7）：50-51.

[3]张晓贵.圆周率计算的四个时期[J].辽宁教育学适合初中生看的书 院院报,2000,(05):66-69.

[4]刘玉玲.圆周率



的简介[J].高等教学研究,2007,(02):62-65.

[5]孙宏安.圆周率计算简史[J].中学数学教学参考，1988,（11）：48-49.

[6]李尚志.数学实验[M].北京：高等教育出版社,2004:28-29.

[7]王铂强,陈军.MonteCarlo方法计算圆周率[J].南通职业大学学

报,2005,(04):61-96.

[8]焦青霞,王俊芳.关于圆周率



值的计算[J].统计与咨询,2008,(06):46-47.

[9]何光.用蒙特卡洛方法计算圆周率的近似值[J].内江师范学院学

报,2007,(04):14-16.

南京师范大学泰州学院本科毕业论文

33

附录

附录3-1

使用软件Mathematica7.0

y[x_]:=4/(1+x^2);(*定义积分函数*)

trapezium[n_]:=((y[0]+y[1])/2+Sum[y[k/n],{k,1,n-1,1}])/n;(*梯形公式计算*)

Simpson[n_]:=(y[0]+y[1]+2Sum[y[k/n],{k,1,n-1,1}]+4Sum[y[(k-0.5)/n],{k,1,n,1}])/(6n)

(*Simpson公式计算*)

N[trapezium[1000]]

N[Simpson[1000]]

wctrapezium=N[trapezium[1000]]-N[Pi]（*计算误差*）

wcsim鹅蛋做法 pson=N[Simpson[1000]]-N[Pi]

表3-1中其他数据计算程序类似，不再一一写出

**************************************************************************************

****************************************************

附录3-2

使用软件Matlab7.0

d=0:6;

h1=line([06],[3.97933.9793])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[03.312643.312663.12753.31265

3.312633.01635];

plot(d,tx,'--d',d,tx)

title('图3-2利用梯形公式计算圆周率随迭代次数增加模拟值的变化图')

xlabel('迭代次数')

ylabel('模拟值')

legend('圆周率值','实验模拟值')

t(gca,'XTicklabel',['0000';'1000';'2000';'3000';'4000';'5000';'6000'])

axis([063.14159243.141592654])

*******************************************************************************d=0:6;

h1=line([06],[00])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[0-1.66667-0.416667-0.185185-0.104167-0.0666667-0.0462963];

stem(d,tx,'--d','fill')

title('图3-3利用梯形公式计算圆周率随迭代次数增加误差的变化趋势图')

xlabel('迭代次数')

ylabel('误差（单位：10的-7次方）')

legend('参考线','实验误差值')

t(gca,'XTicklabel',['0000';'1000';'2000';'3000';'4000';'5000';'6000'])

*******************************************************************************

南京师范大学泰州学院本科毕业论文

34

d=0:6;

h1=line([06],[3.97933.9793])

get(h1)

t(h1,'LineWidth',0.5)

holdon

tx=[03.979053.92312663.97983.98056

3.979053.97927];

plot(d,tx,'--d',d,tx)

title('图3-4利用simpson公式计算圆周率随迭代次数增加模拟值的变化图')

xlabel('迭代次数')

ylabel('模拟值')

legend('圆周率值','实验模拟值')

t(gca,'XTicklabel',['0000';'1000';'2000';'3000';'4000';'5000';'6000'])

axis([063.9793.9805])

*******************************************************************************d=0:8;

h1=line([08],[00])

get(h1)

t(h1,'LineWidth',0.5)

holdon

simp=[0-0.266454-0.4440890.4884980.124345-0.266454-0.04440890.0033540.0015457];

stem(d,simp,'--d','fill')

title('图3-5利用simpson公式计算圆周率随迭代次数增加误差的变化趋势图')

xlabel('迭代次数')

ylabel('误差（单位：10的-15次方）')

legend('参考线','实验误差值')

t(gca,'XTicklabel',['0000';'1000';'2000';'3000';'4000';'5000';'6000';'7000';'8000'])

**************************************************************************************

************************************************************************

附录3-3

使用软件Mathematica7.0

a[n_]:=N[Sum[((-1)^(i-1)/(2i-1))*(1/2)^(2i-1),{i,1,n}],20];

b[n_]:=N[Sum[((-1)^(i-1)/(2i-1))*(1/3)^(2i-1),{i,1,n}],20];

4*(a[1]+b[1])(*展开项数为1时，计算的近似值*)

wc=4*(a[1]+b[1])-N[Pi](*计算误差*)

表3-2中其他数据计算程序类似，不再一一写出

**************************************************************************************

************************************************************************

附录3-4

使用软件Matlab7.0

d=0:9;

h1=line([09],[3.97933.9793])

get(h1)

t(h1,'LineWidth',0.5)

holdon

南京师范大学泰州学院本科毕业论文

35

tx=[03.3333333333.1172839513.1455761323.1408505623.1417411973.1415615883.141599341

3.1415911843.141592981];

plot(d,tx,'--d',d,tx)

title('图3-5利用泰勒级数法计算圆周率随展开项数增加模拟值的变化图')

xlabel('展开项数')

ylabel('模拟圆周率值')

legend('圆周率真实值','实验模拟值')

%t(gca,'XTicklabel',['0000';'1';'2';'3';'4';'5';'6';'7';'8';'9'烤肉蘸料汁怎么调 ])

axis([093.13.333334])

*******************************************************************************

d=0:9;

h1=line([09],[00])

get(h1)

t(h1,'LineWidth',0.5)

holdon

simp=[00.191740680-0..003983478-0..000148544-0.000031066

0.000006687-0..000000328];

stem(d,simp,'--d')

title('图3-6利用taylor公式计算圆周率随迭代次数增加误差的变化趋势图')

xlabel('展开项数')

ylabel('误差')

legend('参考线','实验误差值')

%t(gca,'XTicklabel',['0000';'1000';'2000';'3000';'4000';'5000';'6000';'7000';'8000'])

**************************************************************************************

************************************************************************

附录3-5

使用软件Mathematica7.0

a[n_]:=N[Sum[((-1)^(i-1)/(2i-1))*(1/5)^(2i-1),{i,1,n}],20];

b[n_]:=N[Sum[((-1)^(i-1)/(2i-1))*(1/139)^(2i-1),{i,1,n}],20];

16*a[1]-4*[1](*展开项数为1时，计算的近似值*)

wc=4*(a[1]+b[1])-N[Pi](*计算误差*)

表3-3中其他数据计算程序类似，不再一一写出

**************************************************************************************

************************************************************************

附录3-6

使用软件Matlab7.0

d=0:9;

h1=line([09],[3.97933.9793])

get(h1)

t(h1,'LineWidth',0.5)

holdon

tx=[03.1712230223.1285568513.1295808513.1295515943.1295525043.1295524753.129552476

3.1295524763.129552476];

plot(d,tx,'--d',d,tx)

南京师范大学泰州学院本科毕业论文

36

title('图3-7利用泰勒级数法计算圆周率随展开项数增加模拟值的变化图')

xlabel('展开项数')

ylabel('模拟圆周率值')

legend('圆周率真实值','实验模拟值')

%t(gca,'XTicklabel',['0000';'1';'2';'3';'4';'5';'6';'7';'8';'9'])

axis([093.13.2])

*******************************************************************************d=0:9;

h1=line([09],[00])

get(h1)

t(h1,'LineWidth',0.5)

holdon

simp=[00.029630368-0.013035802-0.012011802-0.012041059-0.012040149-0.012040179

-0.012040178-0.012040178-0.012040178];

stem(d,simp,'--d')

title('图3-8利用taylor公式计算圆周率随展开项数增加误差的变化趋势图')

xlabel('展开项数')

ylabel('误差')

legend('参考线','实验误差值')

%t(gca,'XTicklabel',['0000';'1000';

%'2000';'3000';'4000';'5000';'6000';'7000';'8000'])

*******************************************************************************

d=0:9;

h1=line([09],[3.97933.9793])

get(h1)

t(h1,'LineWidth',0.5)

holdon

tx1=[03.3333333333.1172839513.1455761323.1408505623.1417411973.1415615883.141599341

3.1415911843.141592981];

plot(d,tx1,'-d')

tx2=[03.1712230223.1285568513.1295808513.1295515943.1295525043.1295524753.129552476

3.1295524763.129552476];

plot(d,tx2,':o')

title('图3-8(1)利用泰勒级数法计算圆周率两种展开方式收敛速度比较图')

xlabel('展开项数')

ylabel('模拟圆周率值')

legend('圆周率真实值','方法一实验模拟值','方法二实验模拟值')

%t(gca,'XTicklabel',['0000';'1';'2';'3';'4';'5';'6';'7';'8';'9'])

axis([093.13.333334])

**************************************************************************************

****************************************************

附录3-7

使用软件Mathematica7.0

pic1=Graphics[Line[{{0,1},{1,1},{1,0},{0,0},{0,1}}],AspectRatio1];

pic2=Plot[Sqrt[1-x^2],{x,0,1},AspectRatio1];

南京师范大学泰州学院本科毕业论文

37

pic3=Show[pic1,pic2]

n=2000过敏奶粉 0;t={};

m=0;

Do[x=Random[];y=Random[];

t=AppendTo[t,{x,y}];

If[x^2+y^21,m=m+1],{n}];

pic4=ListPlot[t,PlotStyle{{RGBColor[1,0,0]}}];

Show[pic3,pic4]

N[4m/n]

*********************************************************************

n=2000;p={}；

Do[m=0;Do[x=Random[];y=Random[];If[x^2+y^21,m++],{k,1,n}];

AppendTo[p,N[4m/n,10]],{t,1,20}];

Print[p];

Sum[p[[t]],{t,1,20}]/20

N[Pi,10]

wc=Sum[p[[t]],{t,1,20}]/20-N[Pi,10](*计算误差*)

表3-4中其他数据计算程序类似，不再一一写出

***********************************************************************

附录3-8

使用软件Matlab7.0

d=0:18;

h1=line([018],[3.97933.9793])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[03.1100000003.1530000003.1626666673.1265000003.1296000003.1293333333.123142857

3.1365000003.1395555563.1422000003.1324000003.1477333333.1485000003.147320000

3.1431200003.1411300003.1412466673.142050000];

plot(d,tx,'--d',d,tx)

title('图3-12利用投点法计算圆周率随投点数量增加模拟值的变化图')

xlabel('投点数量')

ylabel('模拟值')

legend('圆周率值','实验模拟值')

t(gca,'XTi眼睛化妆 cklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

axis([0183.13.17])

********************************************************************

d=0:18;

h1=line([018],[00])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[021074013-0.015092654-0.011992654-0.012259320-0.018449796

-0.005092654-0..005727346

南京师范大学泰州学院本科毕业论文

38

0.001527346-0.000462654-0..000457346];

stem(d,tx,'--d')

title('图3-13利用投点法计算圆周率随投点数量增加误差的变化图')

xlabel('投点数量')

ylabel('模拟误差')

legend('参考线','实验模拟误差值')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

%axis([0183.13.17])

**************************************************************************************

************************************************************************附录3-9

使用软件Mathematica7.0

n=100;p={};d=4;l=2;

Do[m=0;

Do[y=Random[Real,{0,d/2}];x=Random[Real,{0,Pi}];

If[y0&&ySin[x],m++],{k,1,n}];AppendTo[p,N[n/m]],{t,1,20}];

Print[p];

Sum[p[[t]],{t,1,20}]/20

wc=Sum[p[[t]],{t,1,20}]/20-N[Pi,10](*计算误差*)

表3-5中其他数据计算程序类似，不再一一写出

**************************************************************************************

************************************************************************

附录3-10

使用软件Matlab7.0

d=0:18;

h1=line([018],[3.97933.9793])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[03..1062355553.2006805253.2074640223.1684384723.1846169693.116460497

3.1479673613.1317560083.1822463743.1633047713.1172984603.1242825053.124265889

3.1267219643.1435165033.1352043433.138897469];

plot(d,tx,'*',d,tx)

title('图3-14利用蒲丰投针法计算圆周率随投针数量增加模拟值的变化图')

xlabel('投针数量')

ylabel('模拟值')

legend('圆周率值','实验模拟值')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

axis([0183.083.22])

*******************************************************************************d=

0:18;

h1=line([018],[00])

get(h1)

t(h1,'LineWidth',0.5)

holdon

南京师范大学泰州学院本科毕业论文

39

tx=[043024316-0.025132156

021712117-0.024294194-0.017310149-0.017326765

-0..001923849-0.006388311-0.002695185];

stem(d,tx,'*')

title('图3-15利用蒲丰投针法计算误差随投针数量增加模拟值的变化图')

xlabel('投针数量')

ylabel('误差值')

legend('参考线','实验模拟值')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

%axis([0183.083.22])

**************************************************************************************

************************************************************************附录3-11

使用软件Mathematica7.0

n=1000;p={};

Do[m=0;

Do[x=Random[Integer,{1,50000}];y=Random[Integer,{1,50000}];

If[GCD[x,y]==1,m++],{k,1,n}];

AppendTo[p,N[(6n/m)^(1/2)]],{t,1,20}];

Print[p];

Sum[p[[t]],{t,1,20}]/20

**************************************************************************************

************************************************************************

附录3-12

使用软件Matlab7.0

d=0:18;

h1=line([018],[3.97933.9793])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx=[03.1227050763.1306873633.1570667703.1329844413.1570345743.1334816433.162234534

3.1667099133.1525770393.1469378393.1416609333.1469558903.1486791223.142863151

3.1443212713.1421845733.1414875433.142078386];

plot(d,tx,'o',d,tx)

title('图3-19利用随机整数互素计算圆周率随取数数量增加模拟值的变化图')

xlabel('取数数量')

ylabel('模拟值')

legend('圆周率值','实验模拟值')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

axis([0183.083.22])

*******************************************************************************

d=0:18;

h1=line([018],[00])

get(h1)

t(h1,'LineWidth',0.5)

南京师范大学泰州学院本科毕业论文

40

holdon

tx=[0-0.018887578-0..015474116-0..015441920-0..020641880

0..000591920-0..000485732];

stem(d,tx,'o')

title('图3-20利用随机整数互素计算误差随取数数量增加模拟值的变化图')

xlabel('取数数量')

ylabel('误差值')

legend('参考线','实验模拟值')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

%axis([0183.083.22])

**************************************************************************************

************************************************************************

附录3-12

使用软件Matlab7.0

d=0:18;

h1=line([018],[3.97933.9793])

get(h1)

t(h1,'LineWidth',2.5)

holdon

tx1=[03.1100000003.1530000003.1626666673.1265000003.1296000003.1293333333.123142857

3.1365000003.1395555563.1422000003.1324000003.1477333333.1485000003.147320000

3.1431200003.1411300003.1412466673.142050000];

plot(d,tx1,'-.dr')

tx2=[03..1062355553.2006805253.2074640223.1684384723.1846169693.116460497

3.1479673613.1317560083.1822463743.1633047713.1172984603.1242825053.124265889

3.1267219643.1435165033.1352043433.138897469];

plot(d,tx2,':*')

tx3=[03.1227050763.1306873633.1570667703.1329844413.1570345743.1334816433.162234534

3.1667099133.1525770393.1469378393.1416609333.1469558903.1486791223.142863151

3.1443212713.1421845733.1414875433.142078386];

plot(d,tx3,'o',d,tx3)

title('图3-21MoteCarlo法三种方式模拟结果的比较图')

xlabel('投点数量')

ylabel('模拟值')

legend('圆周率值','单位圆投点法','蒲丰投针法','随机整数互素法')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

axis([0183.083.22])

*******************************************************************************

d=0:18;

h1=line([018],[00])

get(h1)

t(h1,'LineWidth',0.5)

holdon

南京师范大学泰州学院本科毕业论文

41

tx1=[021074013-0.015092654-0.011992654-0.012259320

-0.018449796-0.005092654-0..000607346-0.009192650..006907346

0..001527346-0.000462654-0..000457346];

stem(d,tx1,'--d')

tx2=[043024316-0.025132156

021712117-0.024294194-0.017310149-0.017326765

-0..001923849-0.006388311-0.002695185];

stem(d,tx2,'*')

tx3=[0-0.018887578-0..015474116-0..015441920-0..020641880

0..000591920-0..000485732];

stem(d,tx3,'o')

title('图3-22MoteCarlo法三种方式模拟误差的比较图')

xlabel('投点数量')

ylabel('误差值')

legend('参考线','单位圆投点法','蒲丰投针法','随机整数互素法')

t(gca,'XTicklabel',['00';'200';'400';'600';'800';'1000';'3000';'5000';'20000';'40000'])

%axis([0183.13.17])