pi的计算 实验报告.docx
《pi的计算 实验报告.docx》由会员分享，可在线阅读，更多相关《pi的计算 实验报告.docx（16页珍藏版）》请在冰豆网上搜索。
pi的计算实验报告
数学应用软件设计
实验报告
班级：
02
姓名：
张海洋
学号:
9
圆周率（选做）
要求：
先查资料看看古人是怎样计算π的，再对π的各种计算方法进行研究和讨论（收敛速度等），并给出不同算法算出π的小数点后第10000位的数字是什么，你觉得该数字应该是多少
第一部分：
圆周率简介
圆周率是指平面上圆的周长与直径之比（ratioofthecircumferenceofacircletothediameter）。
用符号π（读音：
pài）表示。
中国古代有圆率、圆率、周等名称。
它是一个常数（约等于）。
它是一个无理数，即无限不循环小数。
在日常生活中，通常都用代表圆周率去进行近似计算。
而用十位小数便足以应付一般计算。
即使是工程师或物理学家要进行较精密的计算，充其量也只需取值至小数点后几百个位。
计算圆周率的方法
“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。
”历史上最马拉松式的计算，其一是德国的LudolphVanCeulen，他几乎耗尽了一生的时间，计算到圆的内接正262边形，于1609年得到了圆周率的35位精度值，以至于圆周率在德国被称为Ludolph数；其二是英国的WilliamShanks，他耗费了15年的光阴，在1874年算出了圆周率的小数点后707位。
可惜，后人发现，他从第528位开始就算错了。
把圆周率的数值算得这么精确，实际意义并不大。
现代科技领域使用的圆周率值，有十几位已经足够了。
如果用LudolphVanCeulen算出的35位精度的圆周率值，来计算一个能把太阳系包起来的一个圆的周长，误差还不到质子直径的百万分之一。
以前的人计算圆周率，是要探究圆周率是否循环小数。
自从1761年Lambert证明了圆周率是无理数，1882年Lindemann证明了圆周率是超越数后，圆周率的神秘面纱就被揭开了。
在中国，公元263年前后，刘徽提出著名的“割圆术”求出了比较精确的圆周率。
他发现：
当圆内接正多边形的边数不断增加后，多边形的周长会越来越逼近圆周长，而多边形的面积也会越来越逼近圆面积。
于是，刘徽利用正多边形面积和圆面积之间的关系，从正六边形开始，逐步把边数加倍：
正十二边形、正二十四边形，正四十八边形……，一直到正三○七二边形，算出圆周率等于三点一四一六，将圆周率的精度提高到小数点后第四位。
在刘徽研究的基础上，祖冲之进一步地发展，经过既漫长又烦琐的计算，一直算到圆内接正24576边形，而得到一个结论：
＜π＜同时得到π的两个近似分数：
约率为22／7；密率为355／113。
他算出的π的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。
以致于有数学史家提议将这一结果命名为“祖率”。
现在的人计算圆周率,多数是为了验证计算机的计算能力，还有，就是为了兴趣。
第二部分：
古人计算圆周率方法古人计算圆周率，一般是用割圆法。
即用圆的内接或外切正多边形来逼近圆的周长。
Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；LudolphVanCeulen用正262边形得到了35位精度。
这种基于几何的算法计算量大，速度慢，吃力不讨好。
随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。
下面挑选一些经典的常用公式加以介绍。
除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了
1.
Machin公式
这个公式由英国天文学教授JohnMachin于1706年发现。
他利用这个公式计算到了100位的圆周率。
Machin公式每计算一项可以得到位的十进制精度。
因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。
还有很多类似于Machin公式的反正切公式。
在所有这些公式中，Machin公式似乎是最快的了。
虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。
下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。
这些算法用程序实现起来比较复杂。
因为计算过程中涉及两个大数的乘除运算，要用FFT（FastFourierTransform）算法。
FFT可以将两个大数的乘除运算时间由O（n2）缩短为O（nlog（n））。
2、
Ramanujan公式
1914年，印度数学家SrinivasaRamanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。
这个公式每计算一项可以得到8位的十进制精度。
1985年Gosper用这个公式计算到了圆周率的17,500,000位。
3、
AGM（Arithmetic-GeometricMean）算法
Gauss-Legendre公式：
初值：
重复计算：
最后计算：
这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。
1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。
4、
Borwein四次迭代式：
初值：
重复计算：
最后计算：
这个公式由JonathanBorwein和PeterBorwein于1985年发表，它四次收敛于圆周率。
5、
Bailey-Borwein-Plouffe算法
这个公式简称BBP公式，由DavidBailey,PeterBorwein和SimonPlouffe于1995年共同发表。
它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。
这为圆周率的分布式计算提供了可行性。
1997年，FabriceBellard找到了一个比BBP快40％的公式：
第三部分：
对于π的几种计算的研究和讨论：
1、
数值积分法（I）利用积分公式
计算
n=10ans=;
n=20ans=;
n=50ans=;
n=100ans=;
n=200ans=;
n=500ans=;
n=1000ans=;
n=2000ans=;
半径为1的圆称为单位圆，它的面积等于
。
只要计算出单位圆的面积，就算出了
。
在坐标轴上画出以圆点为圆心，以1为半径的单位圆（如下图），则这个单位圆在第一象限的部分是一个扇形，而且面积是单位圆的1/4，于是，我们只要算出此扇形的面积，便可以计算出
。
但计算量较大。
2、
数值积分法（II）利用公式
n=10ans=;
n=20ans=;
n=50ans=;
n=100ans=;
n=200ans=;
n=500ans=;
n=1000ans=;
n=2000ans=;
设分点x1,x2,…xn-1将积分区间[0,1]分成n等分。
所有的曲边梯形的宽度都是h=1/n。
记yi=f（xi）.则第i个曲边梯形的面积A近似地等于梯形面积，即：
A=（y（i-1）+yi）h/2。
将所有这些梯形面积加起来就得到：
A≈2/n[2（y1+y2+…yn-1）+y0+yn]
3、
利用复化梯形算法求Pi的近似值.
=1/2（Tn+h
）
4、
泰勒级数法计算
：
利用反正切函数的泰勒级数
来计算
（I）x=1时
=4*
n=10ans=;
n=20ans=;
n=50ans=;
n=100ans=;
n=200ans=;
n=500ans=;
n=1000ans=;
n=2000ans=;
x=1时得到的
的展开式收敛太慢,逼近速度太慢，运算庞大，对速度造成了很大影响。
5、
泰勒级数法计算
：
利用反正切函数的泰勒级数
来计算
（II）进一步精细化
n=10ans=;
n=20ans=;
n=50ans=;
n=100ans=;
n=200ans=;
n=500ans=;
n=1000ans=;
n=2000ans=;
当x=1时得到的
的展开式收敛太慢。
要使泰勒级数收敛得快，容易想到，应当使x的绝对值小于1，最好是远比1小。
例如，因为
，所以我们可以计算出
的值，从而得到
的值。
这样，就使得收敛速度加快，逼近的速度大大增加.
6、
利用麦琴给出
，推出π=4（
）
n=10ans=;
n=20ans=;
n=50ans=;
n=100ans=;
n=200ans=;
n=500ans=;
n=1000ans=;
n=2000ans=;
对泰勒级数，随着｜x｜的减小，级数的收敛速度明显加快，这启示我们另外构造相关级数来逼近π。
7、
蒙特卡罗法计算
单位圆的1/4是一个扇形，它是边长为1的单位正方形的一部分，单位正方形的面积
。
只要能够求出扇形的面积
在正方形的面积中所占的比例
，就能立即得到
，从而得到
的值。
下面的问题归结为如何求
的值，这就用到了一种利用随机数来解决此种问题的蒙特卡罗法，其原理就是在正方形中随机的投入很多点，是所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。
降落在扇形内的点的个数
与所投店的总数
的比可以近似的作为
的近似值。
n=100ans=;
n=200ans=;
n=500ans=
n=1000ans=;
n=2000ans=;
n=5000ans=;
n=10000ans=;
n=20000ans=;
这种数据模拟算法收敛的速度很慢，从运行结果来看，蒙特卡罗法的计算结果为，虽然精确度不太高，但运行时间短，在很多场合下，特别是在对精确度要求不高的情况下很有用的。
除以上几种方法之，还有下列一些的求其他
的方法：
*利用高斯公式
=48
+32
—20
*、π=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））+……
*、e^（π*i）+1=0（欧拉公式，也称世界上最杰出的公式）
*、1+（1/2）^2+（1/3）^2+（1/4）^2+….（1/n）^2=π^2/6
*、1+（1/2）^4+（1/3）^4+（1/4）^4+….（1/n）^4=π^4/90
*、1+（1/2）^6+（1/3）^6+（1/4）^6+….（1/n）^6=π^6/945
*、1+（1/2）^8+（1/3）^8+（1/4）^8+….（1/n）^8=π^8/9450
*、1+（1/2）^10+（1/3）^10+（1/4）^10+….（1/n）^10=π^10/93555
第四部分：
求
的小数点后第位数字的几种方法：
1、
反正切函数的泰勒级数
得到
；
推出：
=4*
Mathematica程序
小数点后10000位数字为8
2、
沃里斯（Wallis）方法
Mathematica程序
小数点后10000位数字为8
3、
基于
的级数
麦琴给出：
；
推出π=4（
）
MATLAB程序
小数点后10000位数字为8
求取的方法很多，过程类似，不再累述。
验证
·
·
·
结论：
π的小数点后10000位数字为8
第五部分：
心得体会
对于π的值我们早已知晓，但是对于这个数值如何得到却不清楚。
通过这次学习，我们知道原来计算π的值有如此多的方法，并且知道了级数在一些计算中的作用。
通过本次实验，我们知道了运用各种方法计算π的近似值，知道了matlab和mathematica中关于级数的运算， 学习是一个枯燥的过程，有些东西是必须懂的。
也就是一些基础知识吧，包括基础的操作和相关的理论知识。
因此，有一本教材看一次，必定会有所收获。
看了，练习，思考。
MATLAB是一个实用性很强，操作相对容易，比较完善的工具软件。
使用起来比较方便，通过操作可以很快看到结果，能够清晰地感觉到成功与失败，编程问题最头疼的不是编程序，而是调程序，所以在你的程序编完之后，一定要进行验证其正确性，你要尽量多的设想你的问题的复杂性，当然，要一步一步复杂，这样才能保证你的程序的适用性很强。
以后虽然没有数学应用软件设计的课程了，但是对于一个学习数学的学生来说，以后用到MATLAB和mathematica的地方还有很多，我不会因为不用学而不去学，希望能在数学建模之前，把两个软件运用灵便，为以后更强大的数学世界打下夯实的基础。
升级会员 