2-9π的计算

实验九　 的近似计算 【实验目的】 1． 了解圆周率 的计算历程。 2． 了解计算 的割圆术、韦达、级数法、拉马努金公式、迭代法。 3． 学习掌握MATLAB软件有关的命令。 【实验】 利用韦达(VieTa)公式 计算 的近似值。 【实验准备】 是人们经常使用的常数，对 的研究已经持续了2500多年。今天，这种探索还在继续中。 1．割圆术 汉代著名数学家刘徽在《九章算术注》中创造了割圆术。刘徽注意到圆内接正多边形的面积小于圆面积，边数加倍时正多边形的面积随之增加，圆内接正多边形的周长也逐渐逼近圆的周长。这充分体现了朴素极限的思想。利用这种，刘徽计算了圆内接96变形的边长与面积，得到 。 割圆术从单位圆开始，首先做圆的内接正6边形，然后边数加倍，正12边形，正24边形，正48边形，正69边形，等等。利用勾股定理可以建立边数和面积的递推公式，从而得到 的近似值。 图9.1 割圆术 设圆内接正 边形的边长和面积分别为 和 ，根据勾股定理，有如下的递推关系： 。 圆面积 与多边形的面积 之间有如下关系： 。 著名数学家祖冲之（公元429年），计算了圆内接正24576边形的边长，得到了 的取值范围为 。 作单位圆的 边内接和外切正多边形，则单位圆的面积介于内外两个多边形的面积之间，由此也可以建立夹逼方法近似计算 。德国数学家鲁道夫（1610年）用了几乎半生的时间，使用这种方法计算出 的35位小数。为了纪念他，德国人民把 称为鲁道夫数。 图9.2 多变形对圆面积的夹挤 2．韦达(VieTa)公式 1593年，韦达首次给出了计算 的精确表达式： 为达公式看起来有些神秘，其实他的推导过程都是用简单的数学方法： （1） 从 开始 所以，对任意 ，有 故有 令 ，有 取 ，得到 （2） 从 起 用数学归纳法可以证明， 其中有 重根号。结合（1）、（2）两式，可知韦达公式成立。 3． 利用级数计算 （1） 莱布尼茨级数（1674年发现） 1844年，数学家达布在没有计算机的情况下，利用此式计算出 的前200位。使用误差估计公式 计算一下要精确到 的200位小数需要取级数的多少项？由此可以看到达布的工作是多么艰巨。 （2） 欧拉的两个级数（1748年发现） 这两个级数收敛也非常缓慢，计算是实用价值不大。 （3） 基于 的级数 ， 取 时可得， ，即为莱布尼茨级数，直接使用时收敛速度极慢，必须考虑加速算法。 观测级数可知， 得知越接近于0，级数收敛越快。由此可以考虑令 则有 因此， 实际上， 所以 在下面的练习中可以看到，取前25项的部分和就可以使 精确到37位小数，加速效果非常明显。根据这一原理，还可以得到： 高斯公式： 斯托梅公式： 及下面一些类似的公式： ， ， ， ， ， ， ， . 4． 拉马努金（Ranmaunujan）公式 目前，计算 的一个极其有效的公式为 ， 这个级数收敛的非常快，级数每增加一项，可提高大约8位小数的精度。 1985年，数学家比尔.高斯帕依使用这个公式在计算机上算出了 的1750万位小数。这个神奇的公式归功于年轻的传奇数学家拉马努金（Ranmaunujan，1887-1920）。 另一个经过改进的计算公式为 。 级数每增加一项，可提高大约14位小数的精度。 5． 迭代方法 迭代公式1 1989年，Borwein发现了下列收敛于 的迭代公式： 迭代误差可以下式估计 在后面的练习中可以看到，迭代4次可精确到693位小数，迭代8次可精确到178814位小数。 迭代公式2 1996年，Baily发现了另一个收敛于 的迭代公式： 迭代误差可以下式估计 6． 的五百位近似值 这里，我们给出 的两百位近似值： 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491 随着计算机技术的飞速发展和计算方法的突破与创新，计算 的世界纪录正在迅速地被刷新。目前， 的数值已计算到小数点后2061.5843亿位。这一纪录是日本东京大学教授金田康正和他的助手于1999年9月创造的。计算用了37小时32分，检验用了46小时7分。虽然这样高的精度已经没有太多的实际意义，但这反映了现代数学科学的日新月异，反映了人类智慧向极限的挑战。 【实验方法与步骤】 练习1 用刘徽的迭代公式 计算 的近似值。 相应的MATLAB代码为： >>clear; >>x=1; >>for i=1:30 >>x=vpa(sqrt(2-sqrt(4-x^2)),15) %计算精度为15位有效数字 >>S=vpa(3*2^i*x,10) >>end 计算可得 练习２用韦达公式 计算 的近似值。 相应的MATLAB代码为： >>clear; >>x=1; >>for i=1:20 >>i >>x=vpa(x*cos(pi/2^(i+1)),30); %计算精度为30位有效数字 >>pai=vpa(2/x,22) >>error=vpa(pi-2/x,22) >>end 计算可得 1 pai=2.828427124746190097603 error =.313165528843603140860 2 pai=3.061467458920718232298 error =.80125194669075006165e-1 3 pai=3.121445152258052404216 error =.20147501331740834247e-1 4 pai=3.136548490545939249125 error =.5044163043853989338e-2 5 pai =3.140331156954752858963 error =.1261496635040379500e-2 6 pai =3.141277250932772720892 error =.315402657020517571e-3 7 pai =3.141513801144301048318 error =.78852445492190145e-4 8 pai =3.141572940367091461583 error =.19713222701776880e-4 9 pai =3.141587725277159716287 error =.4928312633522176e-5 10 pai =3.141591421511200086037 error =.1232078593152426e-5 11 pai =3.141592345570117830926 error =.308019675407537e-6 12 pai =3.141592576584872625535 error =.77004920612928e-7 13 pai =3.141592634338562821891 error =.19251230416572e-7 14 pai =3.141592648776985437337 error =.4812807801126e-8 15 pai =3.141592652386591008149 error =.1203202230314e-8 16 pai =3.141592653288992575505 error =.300800662958e-9 17 pai =3.141592653514592792966 error =.75200445497e-10 18 pai =3.141592653570993021726 error =.18800216737e-10 19 pai=3.141592653585093078916 error =.4700159547e-11 20 pai =3.141592653588617918820 error =.1175319643e-11 从上面结果可以看出，每增加两项，可以提高1位数的精确度。 练习3用莱布尼茨级数 计算 的近似值。 相应的MATLAB代码为： >>clear; >>x=0; >>for i=1:2000 >>x=vpa(x+(-1)^(i+1)/(2*i-1),30); >>pai(i)=vpa(4*x,22); >>error(i)=vpa(pi-4*x,22); >>end >>for i=1:10 >>i*200 >>pai(i*200) >>error(i*200) >>end 计算可得 200 3.136592684838816750415 .4999968750976488048e-2 400 3.139092657496012721466 .2499996093780516997e-2 600 3.139925988080529960462 .1666665509263278001e-2 800 3.140342654078073534793 .1249999511719703670e-2 1000 3.140592653839792925964 .999999750000312499e-3 1200 3.140759320401135705469 .833333188657532994e-3 1400 3.140878367966615337793 .714285623177900670e-3 1600 3.140967653650828364910 .624999938964873553e-3 1800 3.141037098077104607384 .555555512688631079e-3 2000 3.141092653621043228697 .499999968750009766e-3 上面的结果中，第一列是级数的项数，第二列是 的近似值，第三列式计算误差。从结果可以看到，是用前1000项可以精确到百分位，前2000项可以精确到千分位。 练习4利用级数加速后的公式 计算 的近似值。 相应的MATLAB代码为： >>clear; >>digits(160); %以160位相对精度进行计算 >>syms x >>x=sym(0); >>for k=1:15 >>k >>x=x+sym(-1)^(sym(k)+sym(1))/(sym(2)*sym(k)-sym(1))*(sym(16)/sym(5)^(sym(2)*sym(k)-sym(1))-sym(4)/sym(239)^(sym(2)*sym(k)-sym(1))); >>vpa(x,40) >>vpa(vpa(pi,60)-x,5) >>end 计算可得 1 3.183263598326359832635983263598326359833 -.41671e-1 2 3.140597029326060314304531106579228898150 .99562e-3 3 3.141621029325034425046832517116408069706 -.28376e-4 4 3.141591772182177295018212291112329795027 .88141e-6 5 3.141592682404399517240259836073575860490 -.28815e-7 6 3.141592652615308608149350747666502755367 .97448e-9 7 3.141592653623554761995504593820311845927 -.33762e-10 8 3.141592653588602228662171260486978513156 .11910e-11 9 3.141592653589835847485700672251684395509 -.42609e-13 10 3.141592653589791696917279619620105448141 .15415e-14 11 3.141592653589793294747374857715343543379 -.56285e-16 12 3.141592653589793236391840944671865282509 .20708e-17 13 3.141592653589793238539324592671865282509 -.76681e-19 14 3.141592653589793238459788161264457875102 .28552e-20 15 3.141592653589793238462750207675492357860 -.10682e-21 从结果可以看到，取前15项部分和就可以精确第21位小数。 练习5利用拉马努金公式 计算 的近似值。 相应的MATLAB代码为： >>clear; >>digits(1000); >>syms('x','a1','a2') >>x=sym(2)*sym(2)^sym(0.5)/sym(9801)*sym(1103); >>for i=1:10 >>i >>a1=sym(1); a2=sym(1); >>for j=1:i >>a1=a1*sym(j); >>end >>for j=1:4*i >>a2=a2*sym(j); >>end >>x=x+sym(2)*sym(2)^sym(0.5)/sym(9801)*a2/a1^sym(4)*(sym(1103)+sym(26390)*sym(i))/(sym(396)^(sym(4)*sym(i))); >>vpa((sym(pi)-sym(1/x)),100) >>end 计算可得 1 -.63953e-15 2 -.5682e-23 3 -.5238e-31 4 -.4944e-39 5 -.4741e-47 6 -.4598e-55 7 -.4499e-63 8 -.4433e-71 9 -.4391e-79 10 -.4369e-87 可以看出，级数每增加一项，可提高大约8位小数的精度。 练习6利用迭代公式一 计算 的近似值。 相应的MATLAB代码为： >>clear; >>digits(1000); >>syms('y','z','a') >>y=sym(2)^sym(0.5)-sym(1); a=sym(6)-sym(4)*sym(2)^sym(0.5); >>for k=1:4 >>k >>z=(sym(1)-sym(y)^sym(4))^sym(0.25); >>y=(sym(1)-z)/(sym(1)+z); >>a=(sym(1)+y)^sym(4)*a-sym(2)^(sym(2)*sym(k)+sym(1))*y*(sym(1)+y+y^sym(2)); >>vpa((sym(pi)-sym(1/a)),700) >>end 计算可得 1 .7376e-8 2 .5472e-40 3 .2308e-170 4 .1111e-693 从计算结果可以看到，迭代4次就可以精确到693位小数。根据误差估计式 相应的MATLAB代码为： >>clear; >>digits(1000); >>for n=1:10 >>n >>vpa(sym(164)^sym(n)*sym(exp(1))^(-sym(2)*sym(pi)*sym(4)^sym(n)),10) >>end 计算结果为 1 .1994495307e-8 2 .5883616903e-39 3 .1010085797e-167 4 .1989217109e-689 5 .6783390976e-2783 6 .2079569911e-11163 7 .4164280709e-44692 8 .1518191156e-178813 9 .6081128537e-715306 10 .3541017675e-2861282 说明迭代１０次就可以保证精确到小数２８６１２８２位！ 【练习与思考】 1. 利用勾股定理推导在割圆术一种给出的公式 ， 。 2. 利用韦达公式，构造出一种迭代算法来计算 的近似值，并进行实际计算，评价算法效果。 3. 利用欧拉公式 计算出 的前4位小数。 4．使用高斯公式 和斯托梅公式 计算出 的前60位小数。 5．用改进后的公式 ， 的前 项计算 的值，并列出计算误差。 6．使用Bailey的迭代算法 迭代计算4次，计算 的近似值，观测近似效果。并用迭代误差估计公式 计算出迭代次数和迭代误差界的数据表（ ）。 PAGE 1 _1130237245.unknown _1130247918.unknown _1130248806.unknown _1130250154.unknown _1130250219.unknown _1130389692.unknown _1130389696.unknown _1130267410.unknown _1130388775.unknown _1130389280.unknown _1130389499.unknown _1130389137.unknown _1130308235.unknown _1130388774.unknown _1130250526.unknown _1130250563.unknown _1130249155.unknown _1130249412.unknown _1130250096.unknown _1130249323.unknown _1130248895.unknown _1130248077.unknown _1130248174.unknown _1130248420.unknown _1130248496.unknown _1130248117.unknown _1130247978.unknown _1130248045.unknown _1130247919.unknown _1130237962.unknown _1130238450.unknown _1130238616.unknown _1130247803.unknown _1130238533.unknown _1130238145.unknown _1130238359.unknown _1130238132.unknown _1130237707.unknown _1130237876.unknown _1130237917.unknown _1130237778.unknown _1130237634.unknown _1130237649.unknown _1130237578.unknown _1130221980.unknown _1130236496.unknown _1130236844.unknown _1130237200.unknown _1130236978.unknown _1130237041.unknown _1130236768.unknown _1130236824.unknown _1130236543.unknown _1130236229.unknown _1130236332.unknown _1130236365.unknown _1130236295.unknown _1130236092.unknown _1130236118.unknown _1130236031.unknown _1130221087.unknown _1130221459.unknown _1130221499.unknown _1130221569.unknown _1130221297.unknown _1130221324.unknown _1130221345.unknown _1130221223.unknown _1130220774.unknown _1130221014.unknown _1130221038.unknown _1130219344.unknown _1130219376.unknown _1130219767.unknown _1130219244.unknown