无理数π的计算.doc

无理数的计算谈起数字，人们常常为文人以数字入诗词的才华叫绝。初唐的骆宾王就被人称作“算博士”。骆诗帝京篇有“秦塞重关一百二，汉家离宫三十六”，写得颇有气势。苏东坡与辽邦使臣的掌故众人皆知，苏东坡以“四诗风雅颂”回应“三光日月星”终让辽使称臣，传下一则数字对的佳话。在科学工作者的心目中，准确的数字是他们的终身追求。为提高数字的计算精度，古往今来的许许多多人们为此付出了艰苦的努力。，或者其近似值3.14，从小学的算术课本中，我们就认识了这个数字。而在距今天4000年前，也就是在公元前2000年左右的巴比伦王国这个数字就已经被发现。对这样一个古老的数字，人们已经是再熟悉不过了。现在，在计算机的帮助下，截止到1995年，人们已经将这个古老的常数精确地计算到了六十亿位以上。譬如，利用Matlab的命令，我们可以十分容易得到的小数点后200位的值如下：vpa (pi, 201) ans =3. 1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 59230781640628620899 8628034825 3421170679 8214808651 3282306647 0938446095 50582231725359408128 4811174502 8410270193 8521105559 6446229489 5493038196但是，在数学史上，圆周率的计算，却使得许多的数学家付出了难以形容的艰苦劳动。在计算机技术发展到今天的情况下，我们来介绍一点的计算方法，并自己动手试着计算的值，看看会得到什么样的结论。一古典方法1.阿基米德的方法通过使正多边形外切或内接于圆，将其边数逐步增多来计算圆周长度的方法，很早就用于圆周率的计算了。显然，若取圆的直径D = 1，则圆周长度就是圆周率！希腊数学家阿基米德就用这种方法从计算圆内接正6边形的周长的开始，最后算到圆内接正96边形，发现圆周率大于。进而他又继续计算圆外切正96边形的周长，发现圆周率小于，就是说， 圆周率 ；用近似小数表示就是 3.140845 nSAOBC = S2n S n(SAOBC + ) = S2n + (S2n Sn)从而 S2n S 2S2n Sn 上式就是刘徽用于计算的圆面积不等式。 下面我们用刘徽割圆术讨论计算的具体过程。为简单计，取圆半径等于1,这时圆面积S =，这样式中S的上、下界就是的上、下界。记an = |AB|，有 hn = |OF| = dn = |EF| = 1hn = 1 从图1可见，多边形由n边变为2n之后，其面积在原来的基础上增加了n/2个矩形ABED的面积，有 S2n = Sn + 新的正2n边形的边长为 a2n = = 上述的式就是我们计算的数学模型。从某个圆内接正多边形开始计算，就可以逐步逼近的真值。譬如，从正6边形开始计算是很方便的，这时的边长和面积为 a6 = 1， S6 = 依次使用式便可得到正12边形的面积S12和边长a12 h6 = d6 = 1h6 S12 = S6 a12 = 得到a12后，再次使用式便可得到正24边形的面积和边长，这一迭代过程反复进行，便可以得到边数为48,96,192,的内接正多边形的面积，它们不断地逼近值 ，所求的值上界和下界一定满足不等式。 根据式计算值上、下界的Matlab程序如下n=6;a=1;s=1.5*sqrt(3);sc=;for k=1:25 d=1-sqrt(1-a*a/4); s2=s+n*a*d*0.5; n=n+n;s3=2*s2-s; sc=sc;k,n,s2,s3; a=sqrt(2*d);s=s2;endsc表2是上述程序计算值上、下界的迭代过程记录。 表2 刘徽割圆术估计值的迭代过程次数 边数 下界 上界01 12 3.000000000000000 3.401923788646684 02 24 3.105828541230249 3.211657082460498 03 48 3.132628613281238 3.159428685332228 04 96 3.139350203046867 3.146071792812496 05 192 3.141031950890509 3.142713698734152 06 384 3.141452472285462 3.141872993680414 07 768 3.141557607911857 3.141662743538253 08 1536 3.141583892148318 3.141610176384779 09 3072 3.141590463228050 3.141597034307782 10 6144 3.141592105999271 3.141593748770493 11 12288 3.141592516692157 3.141592927385042 12 24576 3.141592619365383 3.141592722038610 13 49152 3.141592645033690 3.141592670701997 14 98304 3.141592651450767 3.141592657867844 15 196608 3.141592653055036 3.141592654659305 16 393216 3.141592653456104 3.141592653857171 23 50331648 3.141592653589786 3.141592653589810 24 100663296 3.141592653589792 3.141592653589798 25 201326592 3.141592653589794 3.141592653589795 26 402653184 3.141592653589794 3.141592653589795可以看到,随着迭代次数的增加，的有效位数也在增加，大约每三步增加两位，增加的速度很慢。据介绍，没有人能用几何方法手算出40位有效数字以上的值，这是可能的原因。荷兰人鲁道夫为计算35位值，几乎将精力拖垮。 我们在此介绍了计算的早期几何方法，特别详细介绍了如何用刘徽的割圆术来估计的上、下限，显然读者也可以根据“圆面积介于同边数的圆内接与外切正多边形之间”的几何思路来确定值的上、下限。二数值积分法我们知道，还有，所以， ，或者，于是，我们可以通过计算上述定积分的近似值来得到的近似值。 在许多情况下，我们可以用牛顿-莱布尼兹公式计算定积分的值，但在另外一些情况下，譬如被积函数f (x)的原函数不能用初等函数表示，或者如上式中的值为未知的情况下，我们就只能用数值积分法来计算定积分的近似值。 考虑到被积函数f (x)的可积性，用分点a = x0 x1 x21e-5 s=(-1)(n-1)/(2*n-1)+s; n=n+1;ends=4*s,n 为使计算结果精确到小数点后第五位，循环次数n达25000次。 现在看来，计算的级数有明显的缺点：级数收敛太慢，计算量过大。其原因是|x|偏大。现x=，令，显见，记，而，所以，就是 这就是著名的欧拉公式，由欧拉(L.Euler)在1737年提出。又，可取x =，令，则，1,故 ，再令，即，而，就有 这个公式由马庭(J.Machin)于1706年发现，故称为马庭公式。将欧拉公式和马庭公式与arctanx的泰勒级数相结合，会加快该级数的收敛速度，具有很强的实用性。据资料记载1873年英国人W.Shanks用马庭公式与级数结合将算到小数点以后707位，在很长时间内人们认为这是最高记录。二战以后，人们用计算机验算，发现直到小数点以后第527都是正确的。用笔将算到这个程度真是非常不容易，想必他一定付出了极大的努力。下面的计算值公式具有更快的优点。一个是1962年算出10万位值的Shanks等人所使用的公式 类似的公式还有 前者来自数学家高斯(1863)的著作，后者曾被18世纪著名数学家欧拉用来算，据说他只用了一个小时便手算出20位值。四更快的计算值公式如果要计算的值位数不很多，那么沿用前面给出的级数方法一般不会有什么问题。但当要求计算的位数很多时(例如几十万位以上),仅仅只靠计算机硬件的进步，也难在合理的时间内完成作业，因此，我们需要收敛速度更快的计算公式。在历史的长河中，人们用于计算值公式几乎都是线性收敛的，下面我们要介绍当今所获得的2阶及4阶计算值的新公式。1. 算术几何均值(AGM)公式设a0，b0是两个正数，今定义算术均值数列ak、几何均值数列bk 再定义数列ck 若ak和bk当时极限存在并相等，称该极限为(a0, b0)的算术几何均值，记为AGM(a0, b0)。理论推导可知，对于初值 ， 有的计算式 = 且式是2阶收敛的。在实际计算中，可用Sk表示式分母的前k项之和，并用ak和bk分别作为AGM(a0, b0)的估计值就可以算出的上、下限和(显然有)。建议读者自行编程计算值上、下限，并输出计算过程，观察需用几步可以达到理想的效果。2. 波尔文(Borwein)高阶公式在值的高阶算法研究中，最好的结果来自两个都叫波尔文的数学家。他们在1984年发表了一个2阶收敛公式： ， （21）式中。1986年，他们又发现了1/的4阶收敛公式： ， （22） （23）式中。用4阶Borwein公式计算值的Matlab程序如下：a=6-4*sqrt(2);b=sqrt(2)-1;ppi=;for k=1:2b4=(1-b4)(0.25); b=(1-b4)/(1+b4); a=a*(1+b)4-2(2*k+1)*b*(1+b+b2); ppi=ppi;k 1/a; endppi程序中的矩阵变量ppi中的元素为计算的中间过程，输出结果为 1 3.14159264621355 2 3.14159265358979从显示的情况可见第2步的结果已精确到小数点后14位，收敛速度确是很快。五Matlab的数值积分命令一元函数的数值积分命令Matlab提供了在有限区间内，计算函数积分近似值的命令，常用的有下面三个quad，quadl具体的调用格式如下：q=quad (fun, a, b, tol) 采用辛普森公式计算积分。q=quadl (fun , a, b, tol) 采用Lobatto方法求数值积分，计算精度高。说明 被积函数fun，可以是字符串、内联函数，或是M函数文件名。 a,b分别是积分的下限和上限，都是确定的数值。 前三个参数是必要的，后面两个参数可以缺省。 tol表示绝对误差，缺省值是10-6 。例 计算积分，其精确值为0.746 824 . 符号变量法syms x;ss=int(exp(-x2),x,0,1)ss = 1/2*erf(1)*pi(1/2)vpa(ss) %输出32位积分结果ans = .74682413281242702539946743613185 积分指令quad或quadl求积 F=inline(exp(-x.2),x);QQ1=quad(F,0,1),QQ2=quadl(F,0,1)QQ1 = 0.74682418072642QQ2 = 0.74682413398845例 计算曲线 x = cos2t , y = sint , z = t . t0 ,2 的长度L. 曲线弧长的计算公式： L = 作一个名为len.m的m文件，内容如下： function f=len(t) f=sqrt(1+4*(sin(2*t).2+(cos(t).2);在命令窗口键入以下命令g=quad8(len,0,2*pi)运算结果是：g =11.4609二重积分的数值积分命令1.积分限为常数的二重积分命令Q=dblquad(fun, inmin, inmax, outmin, outmax, tol, method)说明 对被积函数fun的要求同指令quad和quadl。 inmin和inmax分别是内变量的下限和上限，outmin和outmax分别是外变量的下限和上限，它们只能是常数。 tol的含义与quad指令中情况相同