本科论文：圆周率额数值计算.doc

乐山师范学院毕业论文（设计）本科生毕业论文（设计）系（院） 专业 论文题目：圆周率的数值计算学生姓名 指导教师 （姓名及职称）班 级 学 号 完成日期：二一四年三月2 圆周率的数值计算姓名学院 专业 学号【摘要】：圆周率在数学中是一个非常重要的常数，受到广泛关注。古今中外一代代的数学家为的计算献出了一生的智慧和辛劳，人类对值的探知过程，反映了数学和计算技术发展的进程。本文先简单介绍了圆周率计算的发展史，继后介绍了利用计算机求圆周率数值解的两种方法的原理。借助Mathematica软件，给出结果、分析误差和比较数据，并分析其优缺点，最后总结了自己所受到的启发和感想。【关键词】：圆周率 计算 数值积分 1 引言1.1 圆周率在数学上的地位在还没有圆周率这个名词的时候，人们还不知道圆的周长和面积该怎么来计算，也许也是可以估算的，像历史上就有很多运用很多种方法来试图计算圆的周长和面积的。比较有名的当数切割法了，于是数学家们渐渐的接近了真相，圆周率问世了，让人们知道了圆周率其实也就是一个比值，是个常数！有了圆周率之后，可以说在数学上以及在历史上都是一个重大的发现，科学家们不用再纠结圆的周长和面积该怎么计算了，它就是一个数学上的里程碑。1.2 研究圆周率数值计算的意义我们都知道在上小学的时候就接触了圆周率，而在今后的更高层的学习中还会接触到圆周率。 比如大学里面的微积分，还有从事航天事业及研究工作的，都需要用到圆周率，而圆周率我们都知道它是一个无理数，也就是说我们无法知道它的具体数值，只能知道它的近似值。所谓近似值，那么它就存在着一定的误差，而有些计算中误差是很重要的，尽量小的误差才是安全的，才是人们所需要的。所以说圆周率数值的研究是非常必要及重要的。1.3 圆周率数值的计算方法在研究圆周率上祖先们花费了很多的时间和心血，从最开始的不知道怎么求值到求近似值，慢慢延伸到接近的数值。经过了一代代数学家呕心沥血的研究，从最初的无法计算到能计算，从能计算到尽量接近真实值，从手工计算到运用计算机计算，这些都是我们宝贵的财富。当然这一切过程都是漫长的，在刘徽的“割圆术”思想中我们看到了很多，也学到了很多，数学以及其他方面的研究，我们的思想要开阔，不要禁锢在已有的思想中。后来的数学家们发明了很多种方法来计算圆周率，为圆周率的计算做了很大的贡献，但我觉得最值得一说的事刘徽。继刘徽之后的数学家们多多少少都受到了“割圆术”的启示，所以刘徽在圆周率计算上算是一大功臣。2 圆周率简介和圆周率计算的发展史2.1 圆周率的简介及发现2.1.1 圆周率的简介 一般地说,圆周率就是圆的周长与直径之比，我们在中学时知道是无理数。 早在4000年前，也就是在公元前2000年左右的巴比伦王国时被发现,它就隐藏在自然数中。不过当时认为圆周率的值是3或。 大约2600年前,“化圆为方”问题的提出，成为了世界三大难题之一。公元前3世纪，古希腊数学家阿基米德（Archimdes,B.C.287-B.C.212）“将圆的半径作为高，将圆周的长度作为底边的三角形的面积就等于圆的面积”的思想，有了圆面积的计算公式为。作为圆周率的符号，它的语源是希腊语的第一个字母。而最早将作为符号使用的人是.乔治（1675-1749），这时是周围的意思。随后哥德巴赫、欧拉都使用作为圆周率的符号。但又据说早在乔治之前的1647年奥托莱多就使用作为圆周率了。所以究竟是谁在什么地方首先使用的，还不知道。计算圆周率的方法有很多，如：割圆术（又称为古典算法）、级数法、弧度法、蒙特卡罗法、数值积分法。2.1.2 圆周率的发现其实圆周率在很早以前就被人们所知道了，从人类诞生在地球起，经历了几千万年的时间。我们看到的太阳几乎总是圆的，我们还发现了很多如圆般的形状。随着人类文明的发展，人类开始研究一些他们疑惑的东西，当然他们对圆就感到比较奇妙，于是发现了无论圆的大小如何，它的周长与直径之比是个定值，这个常数就是圆周率，记作。2.2 圆周率的计算及时期2.2.1 圆周率实验时期发现了圆周率这样一个新名词，人们肯定迫切地想知道它的数值到底是多少，但并没有理论的根据，而是从实验中得到的，按理说得到的数值，精确度不高。公元前2000年的巴比伦人以为圆周率的数值是3或，而稍后时代的埃及人把圆周率的数值定为（已修改成了现在的表示形式）。此公式原式表示在有名的林德纸莎草纸上，写成小数就是3.16049.说到圆周率的计算， 那是经过无数的数学家用了无数的心血得来的成果。那是经历了漫长的岁月的沉淀科学财富！2.2.2 圆周率几何法时期经历了圆周率的实验时期，有了更多的人对圆周率有了各种各样的思考，计算方法也多种多样的，在此时期有一个有趣的故事。被称为“圆面积问题”。在大约2600年以前，希腊的印欧学派有位叫阿那克梅内斯的天文学家。他的弟子中，有位优秀青年叫阿那克萨哥拉，他有着聪明的头脑，总是能正确地判断、解释天文学上的许多现象。那时候的科学不发达，认为天文学上的事都是上天的旨意，也没有人进行科学的测量。而阿那克萨哥拉却对太阳的运行，昼夜的变化，月圆月缺等进行了一系列的调查研究，试图对这些现象的解释说明。但他的努力研究，在别人的眼中无疑是对神明的亵渎，因而受到当时人们的惩罚，被逮捕入狱，并没收全部书。面临着这些阻力，阿那克萨哥拉没有放弃对科学的研究，他思考了一个问题，也就是著名的“圆面积问题”了，他的这种研究精神真是值得敬佩！圆面积问题是世界三大难题之一，“作与圆相等面积的正方形”在几何作图上存在困难。而还有两大难题是“将给定的三角形三等分”、“作出给定立方体的2倍体积的立方体” 之后还有很多数学家及学者研究着几何法计算的问题，下面谈谈我国对圆周率数值计算上的成就。说到我国对圆周率的研究，首先要从刘徽的“割圆术”谈起，刘徽从圆内接正六边形开始算，然后边数一倍倍的增加，当正多边形边数足够多的时候，这时候的正多边形就接近圆形了。用足够多的正多边形的面积数值来逐步逼近圆周率，这样就能无限精确地逼近圆周率，但每一项都比圆周率小，这就是刘徽的“割圆术”的本质。刘徽的“割圆术”用数学语言写出来，就是有一个半径是1的圆O,作内接正六边形（如图1）。正六边形面积是的面积的六倍。由于： 所以，六边形的面积是:刘徽的“割圆术”有两个可贵之处，其一在于：怎样用已知的、可求的来逼近未知的、要求的；另一点在于：他把圆看做边数无穷的正多边形，用有限来逼近无限。这种想法，在后期都起到了重要的作用，以致将可能永远起着重要的作用。继刘徽的“割圆术”之后，有很多数学家就是用这种思想来计算圆周率的。由正多边形得到的的近似值的数学家有阿基米德（B.C.287-B.C.212)、皮沙诺（1175-？）、比埃塔（1540-1603）、罗马奴斯（1561-1615）、卢多夫（1540-1610）等，他们研究的正多边形的边数分别是、，正确值分别是3.14、3.141、小数点后10位、小数点后15位、小数点后35位。2.2.3 圆周率的分析法时期对于的值，发现了很多的无穷级数和连分数，随着的近似值的小数点后的位数的快速增长，用“割圆术”思想计算值的竞争停止了。但真正的竞争并没有结束，数学家们在寻找其他可以快速计算圆周率数值的方法。其中一种方法叫做连分数法，连分数的创始者是数学家盖托蒂。布朗开罗的连分数源于下式：在有限的范围内计算上式，可先用繁分数形式，又由得但计算中取的项数是有限的，导致的值没超过3，存在很大的误差。后来又出现了一种计算圆周率的方法，叫做弧度法。至于是谁发明的弧度法，无法得知，在圆周上取与半径长度相同的弧，把该弧所对的中心角的大小，作为角度单位。如中心角为弧度时，弧长和扇形面积就比较好求了，即：有了弧度法的公式问世，计算圆周率渐渐用上了很多式子，用展开式计算圆周率就是其中一种方法。比较有名的莱布尼茨、欧拉、马庭等展开式，只是莱布尼茨展开式收敛较慢，而欧拉、马庭的展开式，收敛性较其他展开式又好些。还有很多著名的数学家的计算公式和展开式，如：高斯公式：克拉森发表的公式：达泽的公式：莱布尼茨展开式：欧拉展开式：2.2.4 计算机计算时期数学家们有3000多年的时间,是在亲力亲为地研究计算圆周率的，有的甚至花掉了毕生的心血。但毕竟人类的力量是有限的，人工计算圆周率的精度也是有限的，在1948年圆周率算到小数点后808位,达到了手算圆周率的一个新高点。1945年第一台电子计算机问世后,值的计算已经不会再显得那么困难，人们一直在超越知识，超越自己；1949年人们才借助计算机的力量,将圆周率算到2037位。随着数学的发展,科技的进步，又出现了很多公式来计算圆周率，这些都是数学家们的成果，这也给计算机计算圆周率提供了更好的平台。 在用计算机计算圆周率的时代,不断更新的公式、方法让圆周率的数值越发的精准。1958年，达到1万位，仅隔一年就达到16167位，可见计算机的参与，让的位数成指数增长（稍有夸大）。1988年，用巨型计算机就可使的近似值达到以亿计位数，据说1989年，达到了4亿8千万。总之，记录在不断被刷新！我们在实际运算中其实只需取的前几位数就行了,但为什么还一直在进行延长位数的竞争?曾有数学家说过：“历史上一个国家所算得的圆周率的准确程度,可作为衡量这个国家当时数学发展的一个标志。”听着也不无道理。就说如今，人们也以计算机技术是否高明来判断国家的科技水平，而计算的近似值只是一个托而已，这场竞争将一直持续下去，不是说有人的地方就有竞争吗？3 用计算机求圆周率数值的方法十九世纪时，圆周率的计算已经处于炙热状态，而且圆周率小数点后的位数也越来越多，到了二十世纪，随着计算机的发明和科技的发展，圆周率的计算已经不是什么难题，这时候圆周率小数点后的位数更是快速增加，而且近似值也比较准确。下面主要介绍两种用计算机计算圆周率的方法：数值积分法和蒙特卡罗法。3.1 数值积分法3.1.1 算法原理如果一个圆的半径为1，那么我们知道这个圆的面积为，只要知道圆的面积为多少，就可以知道圆周率为多少。问题就在于要怎样算出或者估算出圆的面积，我们学习了积分，知道积分可以估算出一些不规则的形状的面积。我们假如以这个单位圆的圆心为原点，建立直角坐标系。则单位圆在第一象限内的部分G是一个扇形，由曲线及两条坐标轴围成，它的面积。算出了S的近似值，它的4倍便是的近似值。图3.1 扇形G如图3.1，作平行于y轴且将x轴上的区间0,1分成n等份（）的直线，将扇形G分成n个宽度为的部分（）。为曲边梯形：左右的边界是相互平行的直线段，类似梯形的两底，上边界为弧线，称为曲边梯形。若无限大，曲边梯形的上边界就近似成直线段，的面积就近似成一个梯形的面积。高为，两条底边的长分别是和，于是曲边梯形面积的近似值为。扇形面积S的近似值可以看作所有梯形面积之和T： 越大，T就越接近S，4T就越接近值。定积分，就是图中扇形的面积。是的定积分，比便于计算，计算也更合适。 一般地，要计算定积分，也就是计算曲线和直线所围成的面积S。用分点将扇形的半径分成等份，若足够大，这时可将曲边梯形的上边界看作直线段，；曲边梯形宽度都为，记为是第个曲边梯形的面积，则有，当时，有，这就是梯形公式。 同理，把上边界看作是抛物线，就是辛普森公式。3.1.2 结果和误差探讨由于用数值积分法本就是一种估算的方式，是将图形分割成很多分，再用规则图形的计算方式加以计算，所以只要分成的份数足够多，就更接近正确答案。也就是它的误差是取决于所取份数的多少。 我们用Mathematica软件,用梯形公式和Simpson公式计算。可以完成相应的程序编写。详细的计算程序如下：yx_:=4/(1+x2);(*定义积分函数*)trapeziumn_:=(y0+y1)/2+Sumyk/n,k,1,n-1,1)/n; (*梯形公式计算*) Simpsonn_:=(y0+y1+2Sumyk/n,k,1,n-1,1+4Sumy(k-0.5)/n,k,1,n,1 )/(6n) (*Simpson公式计算*)Ntrapezium1000NSimpson1000wctrapezium = Ntrapezium1000-NPi （*计算误差*）wcsimpson = NSimpson1000 -NPi 计算结果见表3.1为方便观察，利用Matlab软件根据表3.1中的数据进行绘图，程序：d = 0:6;h1 = line(0 6,3.141592653589793 3.141592653589793)get(h1)set(h1, LineWidth, 2.5) hold ontx = 0 3.1415924869231264 3.1415926119231266 3.141592635071275 3.1415926431731265 3.1415926469231263 3.1415926489601635;plot(d,tx,-d,d,tx)title(图3.2 )xlabel(迭代次数)ylabel(模拟值)legend(圆周率值,实验模拟值)set(gca,XTicklabel,0000;1000; 2000;3000;4000;5000;6000)axis(0 6 3.1415924 3.141592654)d = 0:6;h1 = line(0 6,0 0)get(h1)set(h1, LineWidth, 2.5) hold ontx = 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(参考线,实验误差值)set(gca,XTicklabel,0000;1000; 2000;3000;4000;5000;6000)*d = 0:6;h1 = line(0 6,3.141592653589793 3.141592653589793)get(h1)set(h1, LineWidth, 0.5) hold ontx = 0 3.1415926535897905 3.141592653589231266 3.141592653589798 3.1415926535898056 3.1415926535897905 3.1415926535897927;plot(d,tx,-d,d,tx)title(图3.4 )xlabel(同上)ylabel(模拟值)legend(圆周率值,实验模拟值)set(gca,XTicklabel,0000;1000; 2000;3000;4000;5000;6000)axis(0 6 3.14159265358979 3.141592653589805)* d = 0:8;h1 = line(0 8,0 0)get(h1)set(h1, LineWidth, 0.5) hold onsimp = 0 -0.266454 -0.444089 0.488498 0.124345 -0.266454 -0.0444089 0.003354 0.0015457;stem(d,simp,-d,fill)title(图3.5 )xlabel(同上)ylabel(误差（单位：10的-15次方）)legend(参考线,实验误差值)set(gca,XTicklabel,0000;1000; 2000;3000;4000;5000;6000;7000;8000)表3.1 利用数值积分法计算圆周率结果及误差迭代次数算法100020003000梯形公式3.14159248692312643.14159261192312663.141592635071275误差-1.6666710-7-4.1666710-8-1.8518510-8Simpson公式3.14159265358979053.14159265358978873.141592653589798误差-2.6645410-15-4.4408910-154.8849810-15400050006000梯形公式3.14159264317312653.14159264692312633.1415926489601635误差-1.0416710-8-6.6666710-9-4.6296310-9Simpson公式3.14159265358980563.14159265358979053.1415926535897927误差1.2434510-14-2.6645410-15-4.4408910-16可得图像3.2至图3.5.通过观察发现，Simpson公式的计算结果要明显优于梯形公式，观察图3.2和图3.3.可以看到利用梯形公式计算时，当迭代次数达到5000以上时，计算结果的精度可达到7位数，而图像3.4和3.5显示，利用Simpson公式计算时，当迭代次数达到5000时，计算结果的精度可达到11位数。图3.2利用梯形公式计算圆周率随迭代次数增加模拟值的变化图图3.3利用梯形公式计算圆周率随迭代次数增加误差的变化趋势图图3.4利用simpson公式计算圆周率随迭代次数增加模拟值的变化图图3.5利用simpson公式计算圆周率随迭代次数增加误差的变化趋势图 _ .-() o( )_-_ 3.2 蒙特卡罗法蒙特卡罗方法是用概率模型来进行近似计算的方法，又称为统计实验法。随着计算机技术的不断进步，蒙特卡罗方法的应用范围也越来越广泛。它能解决很多不同类型的数学和物理问题，在数学研究上有着很高的地位，如定积分的近似计算，在原子能技术、武装装备论证等问题中，也有着很高的应用价值。众多数学软件的相继出现，蒙特卡罗法在数学理论研究方面也越发有了发展空间，不仅能够对已经解决的问题进行强有力的佐证，而且为新的结论提供更好的发展平台。3.2.1 正方形内投点法3.2.1.1 算法原理 通过单位圆与边长为1的正方形面积的比来计算的近似值。具体如下:单位圆的是一个扇形，它是边长为1的正方形的一部分。如图所示，正方形面积为只要能够求出扇形的面积在正方形的面积占的比例，就能立即得到，从而得到的值。 如何求出扇形面积在正方形面积中所占的比例呢？我们可以在正方形中随机的投入很多点，让所投的每个点都是随机的，也就是说点落在正方形中每一个位置的机会是均等的，再看其中有多少个点落在扇形内。将落在扇形内的点的个数与投点的总数的比作为的近似值。 这样就产生两个随机数，则以为坐标的点就是正方形内的一点，它落在正方形内每个位置的机率是均等的。落在扇形内的充要条件是这种利用随机数来解决问题的方法称为蒙特卡罗法。3.2.1.2 计算过程演示通过Mathematica可以编写相应的程序：pic1=GraphicsLine0,1,1,1,1,0,0,0,0,1,AspectRatio1;pic2=PlotSqrt1-x2,x,0,1,AspectRatio1;pic3=Showpic1,pic2n=20000;t=;m=0;Dox=Random;y=Random; t=AppendTot,x,y; Ifx2+y21,m=m+1,n;pic4=ListPlott,PlotStyleRGBColor1,0,0;Showpic3,pic4N4 m/n*n=2000;p=；Dom=0;Dox=Random;y=Random;Ifx2+y21,m+,k,1,n; AppendTop,N4m/n,10,t,1,20;Printp;Sumpt,t,1,20/20NPi,10wc = Sumpt,t,1,20/20- NPi,10 (*计算误差*)演示投点的过程（图3.6图3.8）。图3.6利用投点法计算圆周率时的投点图像（投点数100,500）图3.7利用投点法计算圆周率时的投点图像（投点数1000,5000）图3.8利用投点法计算圆周率时的投点图像（投点数10000，20000）3.2.3.3 计算结果及误差分析通过观察,能看到这种方法获得了一个的近似值,虽精确度不高,但是方法简单,且直观。观察发现,随着所投次数的增加，精确度会变高，由于此方法采用随机数的思想,所以结果不唯一，为尽可能的减小随机因素对结果的影响，在利用程序计算圆周率时，采用如下的流程计算：每投个点得到一个的近似值,将结果放在数组中,同样操作重复20次得到20个近似值,最后用Print显示全部近似值,并求出20个近似值的平均值。我们可以发现，随着投点数的增加，模拟的精度将会提高，即投点数的增加会改善结果的精度，当投点数达到10000次以上后，模拟的精度将达到2位数以上。这样可以减小随机因素的影响。4 结论本文对圆周率进行了简单的介绍，还介绍了计算圆周率的四个时期和一些方法，让人们对计算圆周率的各种方法有初步了解。着重介绍了数值积分法计算圆周率，利用计算机，采用现代方法计算。给大家展现了数学的美妙，激发想象力和创造力，提升学习数学的兴趣。学好数学很重要一点就是要是要自己动手去体验，追求内容的系统性、完整性。在写这篇文章时，我自己也去了解了的发展，并被深深吸引。我在敬佩数学家们的研究精神的同时，在想我们人类的伟大，往往很多东西都是从无到有，我们后辈也应当具备这样不懈的研究精神，这样科技才会发展，社会才会发展，人类才会进步。很多人会说研究是数学家的事，我们只是普通人，其实科研研究不只是数学家们事，我们同样可以研究，哪怕我们的研究是失败的，是渺小的，但至少我们曾经努力过，具有科研精神，有不放弃的精神，只要具备这种精神，那么一切都变得具有意义，知识就是力量，探索是无止境的。我认为这就是好的，值得鼓励的。【参考文献】1 张晓贵.圆周率计算的四个时期J.辽宁教育学院院报,2000年.第一版2 孙宏安.圆周率计算简史J.中学数学教学参考，1988年.第一版3 何光.用蒙特卡洛方法计算圆周率的近似值J.内江