Pi的计算实验
数学实验之——
，无穷的神秘气息探秘
薛东明5030309891张晗雨5030309928
圆周率是一个极其驰名的数。从有文字记载的历史开始，这个数就引进了外行人和学者们的兴趣。作为一个非常重要的常数，圆周率最早是出于解决有关圆的计算问题。仅凭这一点，求出它的尽量准确的近似值，就是一个极其迫切的问题了。事实也是如此，几千年来作为数学家们的奋斗目标，古今中外一代一代的数学家为此献出了自己的智慧和劳动。回顾历史，人类对π 的认识过程，反映了数学和计算技术发展情形的一个侧面。π 的研究，在一定程度上反映这个地区或时代的数学水平。德国数学史家康托说：“历史上一个国家所算得的圆周率的准确程度，可以作为衡量这个国家当时数学发展水平的指标。”直到19世纪初，求圆周率的值应该说是数学中的头号难题。之所以如此，看看下面的方程就会明白。这只是一小部分，后面有更多。下面我和我的同伴就计算π的历史过程讨论一下π值的计算方法。
几何学：
若圆的半径为r，圆周为C = 2 πr
若圆的半径为r，其面积为A = πr2
若球的半径为r，其体积为V = (4/3) πr3
若球的半径为r，其表面积为r: A = 4 πr2
分析数学：
数论：
任意两个自然数，互质的概率是6/π/π。
概率论
取一枚长为l的针，再取一张白纸在上面画上一些距离为2的平行线。
把针从一定高度释放，让其自由落体到纸面上。针与平行线相交的概率
是圆周率的倒数。（蒲丰投针）
物理学
(海森堡测不准原理)
（爱因斯坦相对论场方程）
实验时期计算π
通过实验对π 值进行估算，这是计算π 的的第一阶段。这种对π 值的估算基本上都是以观察或实验为根据，是基于对一个圆的周长和直径的实际测量而得出的。在古代世界，实际上长期使用π ＝3这个数值。最早见于文字记载的有基督教《圣经》中的章节，其上取圆周率为3。这一段描述的事大约发生在公元前950
年前后。其他如巴比伦、印度、中国等也长期使用3这个粗略而简单实用的数值。在我国刘徽之前“圆径一而周三”曾广泛流传。我国第一部《周髀算经》中，就记载有圆“周三径一”这一结论。在我国，木工师傅有两句从古流传下来的口诀：叫做：“周三径一，方五斜七”，意思是说，直径为1的圆，周长大约是3，边长为5的正方形，对角线之长约为7。这正反映了早期人们对圆周率π 和√2 这两个无理数的粗略估计。东汉时期官方还明文规定圆周率取3为计算面积的标准。后人称之为“古率”。
早期的人们还使用了其它的粗糙方法。如古埃及、古希腊人曾用谷粒摆在圆形上，以数粒数与方形对比的方法取得数值。或用匀重木板锯成圆形和方形以秤量对比取值……由此，得到圆周率的稍好些的值。如古埃及人应用了约四千年的 4 (8/9)2 = 3.1605。在印度，公元前六世纪，曾取π= √10 = 3.162。在我国东、西汉之交，新朝王莽令刘歆制造量的容器――律嘉量斛。刘歆在制造标准容器的过程中就需要用到圆周率的值。为此，他大约也是通过做实验，得到一些关于圆周率的并不划一的近似值。现在根据铭文推算，其计算值分别取为3.1547，3.1992，3.1498，3.2031比径一周三的古率已有所进步。人类的这种探索的结果，当主要估计圆田面积时，对生产没有太大影响，但以此来制造器皿或其它计算就不合适了
几何法计算π
凭直观推测或实物度量，来计算π 值的实验方法所得到的结果是相当粗略的。真正使圆周率计算建立在科学的基础上，首先应归功于阿基米德。他是科学地研究这一常数的第一个人，是他首先提出了一种能够借助数学过程而不是通过测量的、能够把π 的值精确到任意精度的方法。由此，开创了圆周率计算的第二阶段。
圆周长大于内接正四边形而小于外切正四边形，因此2√2 ＜π ＜ 4 。当然，这是一个差
劲透顶的例子。据说阿基米德用到了正96边形才算出他的值域。
阿基米德求圆周率的更精确近似值的方法，体现在他的一篇论文《圆的测定》之中。在这一书中，阿基米德第一次创用上、下界来确定π 的近似值，他用几何方法证明了“圆周长与圆直径之比小于3+(1/7) 而大于3 +
(10/71) ”，他还提供了误差的估计。重要的是，这种方法从理论上而言，能够求得圆周率的更准确的值。到公元150年左右，希腊天文学家托勒密得出π ＝3.1416，取得了自阿基米德以来的巨大进步。
割圆术。不断地利用勾股定理，来计算正N边形的边长。
在我国，首先是由数学家刘徽得出较精确的圆周率。公元263年前后，刘徽提出著名的割圆术，得出π ＝3.14，通常称为“徽率”，他指出这是不足近似值。虽然他提出割圆术的时间比阿基米德晚一些，但其方法确有着较阿基米德方法更美妙之处。割圆术仅用内接正多边形就确定出了圆周率的上、下界，比阿基米德用内接同时又用外切正多边形简捷得多。另外，有人认为在割圆术中刘徽提供了一种绝妙的精加工办法，以致于他将割到192边形的几个粗糙的近似值通过简单的加权平均，竟然获得具有4位有效数字的圆周率π ＝
3927/1250 ＝3.1416。而这一结果，正如刘徽本人指出的，如果通过割圆计算得出这个结果，需要割到3072
边形。这种精加工方法的效果是奇妙的。这一神奇的精加工技术是割圆术中最为精彩的部分，令人遗憾的是，由于人们对它缺乏理解而被长期埋没了。
恐怕大家更加熟悉的是祖冲之所做出的贡献吧。对此，《隋书·律历志》有如下记载：“宋末，南徐州从事祖冲之更开密法。以圆径一亿为丈，圆周盈数三丈一尺四寸一分五厘九毫二秒七忽，朒数三丈一尺四寸一分五厘九毫二秒六忽，正数在盈朒二限之间。密率：圆径一百一十三，圆周三百五十五。约率，圆径七，周二十二。”
这一记录指出，祖冲之关于圆周率的两大贡献。其一是求得圆周率
3.1415926 ＜π ＜3.1415927
其二是，得到π 的两个近似分数即：约率为22／7；密率为355／113。
他算出的π 的8位可靠数字，不但在当时是最精密的圆周率，而且保持世界记录九百多年。以致于有数学史家提议将这一结果命名为“祖率”。
这一结果是如何获得的呢？追根溯源，正是基于对刘徽割圆术的继承与发展，祖冲之才能得到这一非凡的成果。因而当我们称颂祖冲之的功绩时，不要忘记他的成就的取得是因为他站在数学伟人刘徽的肩膀上的缘故。后人曾推算若要单纯地通过计算圆内接多边形边长的话，得到这一结果，需要算到圆内接正12288边形，才能得到这样精确度的值。祖冲之是否还使用了其它的巧妙办法来简化计算呢？这已经不得而知，因为记载其研究成果的著作《缀术》早已失传了。这在中国数学发展史上是一件极令人痛惜的事。
割圆术的Mathmatica直观表示：
For[k = 0, k≤5 ,k++, n= 6*2^k;
v1 = Graphics[{RGBColor[1, 0,0], Circle[{0,0} , 1]}];
v2 = Graphics[{RGBColor[0, 1,0], Line[Table[{Cos[2Pi*i/n],Sin[2Pi*i/n]}, {i, 0, n}]]}];
v3 = Graphics[{RGBColor[0,0,1], Line[Table[{Cos[2Pi*i/n]/Cos[pi/n], Sin[2Pi /n]/Cos[Pi/n]}, {i,0, n}]]}];
Show[v3, v1, v2,AspectRatio → Automatic]]
图像：
下面我们用多边形割圆术计算π
Print[" a1 " , " a2 ", " r", " n"]; For[k = 0, k ≤ 15, k++, n = 6 * 2^k; a1 = N[n/2 * Sin[2 Pi/n ], 20]; a2 = N[n * Tan[Pi/n], 20];r= a2- a1; Print[" ",a1, " ", a2, " ",r, " ",n]]
结果
a1 a2 r n
2.5980762113533159403
3.4641016151377545871 0.8660254037844386468 6
3.0000000000000000000 3.2153903091734724777 0.2153903091734724777 12
3.1058285412302491482 3.1596599420975004833 0.0538314008672513351 24
3.1326286132812381972 3.1460862151314349711 0.0134576018501967739 48
3.1393502030468672071 3.1427145996453682982 0.0033643965985010910 96
3.1410319508905096381 3.1418730499798238717 0.0008410990893142336 192
3.1414524722854620755 3.1416627470568485262 0.0002102747713864508 384
3.1415576079118576455 3.1416101766046895388 0.0000525686928318932 768
3.1415838921483184087 3.1415970343215261520 0.0000131421732077433 1536
3.1415904632280500957 3.1415937487713520280 3.2855433019322 106
3.1415921059992715505 3.14159292738509703358.213858254830 107
3.1415925166921574476 3.1415927220386138183 2.0534******** 107
3.1415926193653839552 3.1415926707019980479 5.133******** 108
3.1415926450336908967 3.1415926578678444198 1.28341535232 108
3.1415926514507676517 3.1415926546593060325 3.2085383808 109
3.1415926530550368417 3.14159265385717143698.021345952 1010
此种方法可以夹出Pi值
Pi = 3.141592653
数值积分方法计算π：
由于计算数学分析方法很多我们只就一些方法讨论：
此种方法已作出VC程序。
（韦达给出的1593年）
数值积分方法： 014 1
梯形法
a = 0; b= 1;y[x_]:= 4/(1 + x^2)
n = 1000
pis1 = N[(b - a )/n * (Sum[y[a + i * (b-a)/n], {i, 1, n-1}] + (y[a] + y[b]) / 2), 100]
输出：
3.1415924869231265717979608435969622546877898925997631081083727313902881097905625 85631933321397748830
两个端点和中点所成的二次曲线的面积和（辛普森公式）：
a = 0; b= 1;y[x_]:= 4/(1 + x^2)
n = 1000;
pis2 = N[(b-a)/6/n*((y[a] + y[b]) + 2 * Sum[y[a + i*(b - a)/n], {i, 1, n-1}] + 4*Sum[y[a + (i - 1/2) * (b
-a )/n],{i, 1, n}]), 100]
输出：
3.1415926535897932384620233435969635160810190557944194802104971343099602108249154
46285805237414852023
幂级数展开： k 0 1 k x 2k 1 当x = 1时，
4 1 k
For[i = 0, i ≤ 5000, i = i+500; d = N[4Sum[(-1)^j/(2j + 1 ), {j, 0, i }], 5]; r = N[ 4/ ( 2i + 3), 5]; Print["i = ", i, " Pi = ", d, " ri = ", r] ]
输出
i = 500 Pi = 3.1436 ri = 0.0039880
i = 1000 Pi = 3.1426 ri = 0.0019970
i = 1500 Pi = 3.1423 ri = 0.0013320
i = 2000 Pi = 3.1421 ri = 0.00099925
i = 2500 Pi = 3.1420 ri = 0.00079952
i = 3000 Pi = 3.1419 ri = 0.00066633
i = 3500 Pi = 3.1419 ri = 0.00057118
i = 4000 Pi = 3.1418 ri = 0.00049981
i = 4500 Pi = 3.1418 ri = 0.00044430
i = 5000 Pi = 3.1418 ri = 0.00039988
i = 5500 Pi = 3.1418 ri = 0.00036354
可以看出收敛的速度极慢，不适合计算π值。 其他级数法：
其中级数法已做出了VC 程序。
泰勒级数展开法
n = 10000;
pimaqin = N[16 * Sum[(-1)^(k - 1) * (1/5)^ (2k -1)/(2k -1), {k ,1, n}] - 4*Sum[(-1)^ (k -1) * (1/239)^(2k -1)/(2k-1),{k, 1 ,n}] , 1000]
输出：
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378
387528865875332083814206171776691473035982534904287554687311595628638823537875937519577 81857780532171226806613001927876611195909216420199
另一种展开法：
n = 100;
pimaqin = N[4* Sum[(-1)^(k - 1) * (1/2)^ (2k -1)/(2k -1), {k ,1, n}] +4*Sum[(-1)^ (k -1) * (1/5)^(2k -1)/(2k-1),{k, 1 ,n}] + 4* Sum[(-1)^(k - 1) * (1/8)^ (2k -1)/(2k -1), {k ,1, n}] , 1000]
输出:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280 348253421170679821480865132823066470938446095505822317253594081284811174502841027019385 211055596446229489549303819644288109756659334461284756482337867831652712019091456485669 234603486104543266482133936072602491412737245870066063155881748815209209628292540917153 643678925903600113305305488204665213841469519415116094330572703657595919530921861173819 326117931051185480744623799627495673518857527248912279381830119491298336733624406566430 860213949463952247371907021798609437027705392171762931767523846748184676694051320005681 271452635608277857713427577896091736371787214684409012249534301465495853710507922796892 589235420199561121290219608640344181598136297747713099605187072113499999983729780499510 597317328160963185950244594553469083026425223082533446850352619311881710100031378387528 865875332083814206171776691473035982534904287554687311595628638823537875937519577818577 80532171226806613001927876611195909216420199
以上两种方法比第一个收敛得快得多。
这里还有一种更快的方法：
π/4 = 12tan(1/38) + 20tan(1/57) + 7tan(1/239)+24tan(1/268)
程序：
n = 1000;
pimaqin = N[48 * Sum[(-1)^(k - 1) * (1/38)^ (2k -1)/(2k -1), {k ,1, n}] + 80*Sum[(-1)^ (k -1) * (1/57)^(2k -1)/(2k-1), {k ,1, n}]+28* Sum[(-1)^(k - 1) * (1/239)^ (2k -1)/(2k -1), {k , 1, n}] + 96*Sum[(-1)^ (k -1) * (1/268)^(2k -1)/(2k-1),{k, 1 ,n}] , 1000]
输出：
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089 986280348253421170679821480865132823066470938446095505822317253594081284811174502841027 019385211055596446229489549303819644288109756659334461284756482337867831652712019091456 485669234603486104543266482133936072602491412737245870066063155881748815209209628292540 917153643678925903600113305305488204665213841469519415116094330572703657595919530921861 173819326117931051185480744623799627495673518857527248912279381830119491298336733624406 566430860213949463952247371907021798609437027705392171762931767523846748184676694051320 005681271452635608277857713427577896091736371787214684409012249534301465495853710507922 796892589235420199561121290219608640344181598136297747713099605187072113499999983729780 499510597317328160963185950244594553469083026425223082533446850352619311881710100031378 387528865875332083814206171776691473035982534904287554687311595628638823537875937519577 81857780532171226806613001927876611195909216420199
此种方法收敛性极好，当n = 10时，输出
3.1415926535897932384626433832795013586099566818871243928967600620507409280583844 38167912417331072689
完全正确。可以正确的计算出100位以上。
其实，级数法的展开式arctg法可以有很多表示方法。即是每次用tan的二倍角公式解二元二次方程。用左右两端夹逼出来。
π有非常美观的连分数表示方法：
其中以经选择了第二种方法做出了VC程序
蒙特卡罗法：
圆中投点法
n =1000;
pimontecario1 = Block[{i , m = 0}, For[i = n, i >0, i--, m = m + If[Random[]^2 + Random[]^2 1, 1,0]]; N[4*m/n, 10]]
输出：
3.138232000
其实，不一定要圆中投点，任何和圆有关的都可以计算。
如立体的圆球投点法，球面投点法。
π的其它计算方法：
蒲丰投针法
在1777年出小版的《或然性算术实验》一书中，蒲丰提出了用实验方法计算π 。这个实验方法的操作很简单：找一根粗细均匀，长度为 d 的细针，并在一张白纸上画上一组间距为l 的平行线（方便起见，常取l = d/2），然后一次又一次地将针任意投掷在白纸上。这样反复地投多次，数数针与任意平行线相交的次数。
于是就可以得到π 的近似值。因为蒲丰本人证明了针与任意平行线相交的概率为p = 2l/πd 。
不过，蒲丰实验的重要性并非是为了求得比其它方法更精确的π 值。蒲丰投针问题的重要性在于它是第一个用几何形式表达概率问题的例子。计算π 的这一方法，不但因其新颖，奇妙而让人叫绝，而且它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导。
随机整数互质法：
在用概率方法计算π 值中还要提到的是：R·查特在1904年发现，两个随意写出的数中，互素的概率为6／π2。
此种方法已经作了VC程序。
自然现象：
1995年4月英国《自然》杂志刊登文章，介绍英国伯明翰市阿斯顿大学计算机科学与应用数学系的罗伯特·马修斯，如何利用夜空中亮星的分布来计算圆周率。马修斯从100颗最亮的星星中随意选取一对又一对进行分析，计算它们位置之间的角距。他检查了100万对因子，据此求得π 的值约为3.12772。这个值与真值相对误差不超过5％。
计算机出现后的计算机方法：
1946年，世界第一台计算机ENIAC制造成功，标志着人类历史迈入了电脑时代。电脑的出现导致了计算方面的根本革命。1949年，ENIAC根据梅钦公式计算到2035（一说是2037）位小数，包括准备和整理时间在内仅用了70小时。计算机的发展一日千里，其记录也就被频频打破。
ENIAC：一个时代的开始
1973年，有人就把圆周率算到了小数点后100万位，并将结果印成一本二百页厚的书，可谓世界上最枯燥无味的书了。1989年突破10亿大关，1995年10月超过64亿位。1999年9月30日，《文摘报》报道，日本东京大学教授金田康正已求到2061.5843亿位的小数值。如果将这些数字打印在A4大小的复印纸上，令每页印2万位数字，那么，这些纸摞起来将高达五六百米。来自最新的报道：金田康正利用一台超级计算机，计算出圆周率小数点后一兆二千四百一十一亿位数，改写了他本人两年前创造的纪录。据悉，金田教授与日立制作所的员工合作，利用目前计算能力居世界第二十六位的超级计算机，使用新的计算方法，耗时四百多个小时，才计算出新的数位，比他一九九九年九月计算出的小数点后二千六百一十一位提高了六倍。圆周率小数点后第一兆位数是二，第一兆二千四百一十一亿位数为五。如果一秒钟读一位数，大约四万年后才能读完。
现在科学家计算π值的动力：
现在打破记录，不管推进到多少位，也不会令人感到特别的惊奇了。实际上，把π 的数值算得过分精确，应用意义并不大。现代科技领域使用的π 值，有十几位已经足够。如果用鲁道夫的35位小数的π 值计算一个能把太阳系包围起来的圆的周长，误差还不到质子直径的百万分之一。我们还可以引美国天文学家西蒙·纽克姆的话来说明这种计算的实用价值：
“十位小数就足以使地球周界准确到一英寸以内，三十位小数便能使整个可见宇宙的四周准确到连最强大的显微镜都不能分辨的一个量。”
那么为什么数学家们还象登山运动员那样，奋力向上攀登，一直求下去而不是停止对π 的探索呢？为什么其小数值有如此的魅力呢？
这其中大概免不了有人类的好奇心与领先于人的心态作怪，但除此之外，还有许多其它原因。
奔腾与圆周率之间的奇妙关系……
1、它现在可以被人们用来测试或检验超级计算机的各项性能，特别是运算速度与计算过程的稳定性。这
对计算机本身的改进至关重要。就在几年前，当Intel公司推出奔腾（Pentium）时，发现它有一点小问题，这问题正是通过运行π 的计算而找到的。这正是超高精度的π 计算直到今天仍然有重要意义的原因之一。
2、计算的方法和思路可以引发新的概念和思想。虽然计算机的计算速度超出任何人的想象，但毕竟还需
要由数学家去编制程序，指导计算机正确运算。实际上，确切地说，当我们把π 的计算历史划分出一个电子计算机时期时，这并非意味着计算方法上的改进，而只是计算工具有了一个大飞跃而已。因而如何改进计算技术，研究出更好的计算公式，使公式收敛得更快、能极快地达到较大的精确度仍是数学家们面对的一个重要课题。在这方面，本世纪印度天才数学家拉马努扬得出了一些很好的结果。
他发现了许多能够迅速而精确地计算π 近似值的公式。他的见解开通了更有效地计算π 近似值的思路。现在计算机计算π 值的公式就是由他得到的。至于这位极富传奇色彩的数学家的故事，在这本小书中我们不想多做介绍了。不过，我希望大家能够明白π 的故事讲述的是人类的胜利，而不是机器的胜利。
3、还有一个关于π 的计算的问题是：我们能否无限地继续算下去？答案是：不行！根据朱达偌夫斯基
的估计，我们最多算1077位。虽然，现在我们离这一极限还相差很远很远，但这毕竟是一个界限。为了不受这一界限的约束，就需要从计算理论上有新的突破。前面我们所提到的计算，不管用什么公式
都必须从头算起，一旦前面的某一位出错，后面的数值完全没有意义。还记得令人遗憾的谢克斯吗？
他就是历史上最惨痛的教训。
4、于是，有人想能否计算时不从头开始，而是从半截开始呢？这一根本性的想法就是寻找并行算法公式。
1996年，圆周率的并行算法公式终于找到，但这是一个16进位的公式，这样很容易得出的1000亿位的数值，只不过是16进位的。是否有10进位的并行计算公式，仍是未来数学的一大难题。
5、作为一个无穷数列，数学家感兴趣的把π 展开到上亿位，能够提供充足的数据来验证人们所提出的
某些理论问题，可以发现许多迷人的性质。如，在π 的十进展开中，10个数字，哪些比较稀，哪些比较密？π 的数字展开中某些数字出现的频率会比另一些高吗？或许它们并非完全随意？这样的想法并非是无聊之举。只有那些思想敏锐的人才会问这种貌似简单，许多人司空见惯但却不屑发问的问题。
6、数学家弗格森最早有过这种猜想：在π 的数值式中各数码出现的概率相同。正是他的这个猜想为发
现和纠正向克斯计算π 值的错误立下了汗马功劳。然而，猜想并不等于现实。弗格森想验证它，却无能为力。后人也想验证它，也是苦于已知的π 值的位数太少。甚至当位数太少时，人们有理由对猜想的正确性做出怀疑。如，数字0的出现机会在开始时就非常少。前50位中只有1个0，第一次出现在32位上。可是，这种现象随着数据的增多，很快就改变了：100位以内有8个0；200位以内有19个0；……1000万位以内有999，440个0；……60亿位以内有599，963，005个0，几乎占1／10。
其他数字又如何呢？结果显示，每一个都差不多是1／10，有的多一点，有的少一点。虽然有些偏差，但都在1／10000之内。
7、人们还想知道：π 的数字展开真的没有一定的模式吗？我们希望能够在十进制展开式中通过研究数
字的统计分布，寻找任何可能的模型――如果存在这种模型的话，迄今为止尚未发现有这种模型。同时我们还想了解：π 的展开式中含有无穷的样式变化吗？或者说，是否任何形式的数字排列都会出现呢？著名数学家希尔伯特在没有发表的笔记本中曾提出下面的问题：π 的十进展开中是否有10个9连在一起？以现在算到的60亿位数字来看，已经出现：连续6个9连在一起。希尔伯特的问题答案似乎应该是肯定的，看来任何数字的排列都应该出现，只是什么时候出现而已。但这还需要更多π 的数位的计算才能提供切实的证据。
8、在这方面，还有如下的统计结果：在60亿数字中已出现连在一起的8个8；9个7；10个6；小数点
后第710150位与3204765位开始，均连续出现了七个3；小数点52638位起连续出现了14142135这八个数字，这恰是的前八位；小数点后第2747956位起，出现了有趣的数列876543210，遗憾的是前面缺个9；还有更有趣的数列123456789也出现了。
如果继续算下去，看来各种类型的数字列组合可能都会出现。通过几何、微积分、概率等广泛的范围和渠道发现π ，这充分显示了数学方法的奇异美。π 竟然与这么些表面看来风马牛不相及的试验，沟通在一起，这的确使人惊讶不已。难怪有这么多人探索他的奥妙。
计算方法上机实验报告
《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日
前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求 在 附近的数值解 ，并使其满足 . 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x
的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把 ) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果： 实验二：
PI值计算方法1
计算圆周率的一些公式 -|waruqi 发表于 2005-12-8 9:24:00 Machin公式 这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 还有很多类似于Machin公式的反正切公式: pi/4=arctg(1/2)+arctg(1/5)+ arctg(1/8) 1844.达塞利 = arctg(1/2)+ arctg(1/3) =2 arctg(1/3)+ arctg(1/7) =12 arctg(1/18)+8 arctg(1/57)-5 arctg(1/239) 在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Four ier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nl og(n))。 （FFT算法不在此文讲诉）
Ramanujan公式 1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共1 4条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。 1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为： 这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky 公式的另一个更方便于计算机编程的形式是： AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式： 初值： 重复计算：
数学实验报告
《数学实验》报告 题目：根据数值积分计算方法计 算山东省面积 学生姓名： 学号： 专业班级：机械工程17-1班
2019年4月15日
一、问题背景与提出 图1是从百度地图中截取的山东省地图，试根据前面数值积分计 算方法，计算山东省面积。 图 1 二、实验目的 1、 学会运用matlab 解决一些简单的数学应用问题。 2、 学会运用matlab 建立数学模型。 3、 学会运用一些常见的数值积分计算方法结算实际问题，并 了解其实际意义，建立积分模型。 三、实验原理与数学模型 将积分区间 [a , b] n 等分，每个区间宽度均为h = (b - a) / n , h 称 为积分步长。记 a = x 0 < x 1 < … < x k … < x n = b , 在小区间上用小矩形面积近似小曲边梯形的面积，若分别取左端点和右端点的函数值为小矩形的高，则分别得到两个曲边梯形的面积的近似公式： Ln = h ∑f (x k )n=1k=0 , h = b?a ?
R n =?∑f (x k )n k=1 , h = b?a ? 如果将二者求平均值，则每个小区间上的小矩形变为小梯形，整 个区间上的值变为： Tn =?∑f (X k )n=1 k=1+?2[f (x 0)+f (x n )] 将山东省边界上的点反映在坐标化，运用梯形公式积分计算得山 东省的面积。 四、实验内容（要点） 1、将山东省的地图区域在matlab 中画出 。 2、在坐标系上运用积分方法将所求区域的面积求出。 3、通过比例尺将山东省的实际面积求出。 五、实验过程记录（含基本步骤、主要程序清单及异常情况记录等） 1、 在百度地图中标识出山东省的区域范围，标明对应的比例： 图 2 2、 取出所截取图片中山东的边界的坐标，即将边界坐标化： （1） 运用imread 函数和imshow 函数导入山东省的区域 图片。
计算pi
一、实验目的 探索精确计算π值的方法，并且比较不同方法之间的不同之处和优缺点。掌握数值积分的辛普森公式。 二、问题描述 1. 任务1 1) 用反正切函数的幂级数展开式结合有关公式求π，若要精确到40位、50位数 字，试比较简单公式和Machin 公式所用的项数。 2) 验证公式 111=arctan arctan arctan 4 258π ++ 试试此公式右端做幂级数展开完成任务1所需要的步数。 2. 任务2 用数值积分计算π，分别用梯形法和Simpson 法精确到10位数字，用Simpson 法精确到15位数字。 3. 任务3 用Monte Carlo 法计算π，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得5位精确数字？ 设计方案用计算机模拟Buffon 实验 4. 任务4 利用积分 2 0(1)!!sin !!2 n n xdx n π π-=? ，n 为奇数 推导公式 224422213352121 n n n n π=-+ ……… 用此公式计算π的近似值，效果如何？ 5. 任务5 利用学过的知识（或查阅资料），提出其他计算π的方法（先用你学过的知识证明），然后实践这种方法。 对你在实验中应用的计算π的方法进行比较分析。 6. 任务6 e 是一个重要的超越数 1e lim 1)n n n →∞=+（ 1111...2!!(1)! e e n n θ =++++++ 试用上述公式或其他方法近似计算e 。
三、问题解法 1. 任务1 1) 根据幂级数展开的相关知识，易知： 24122211(1)1n n x x x x --=-+-+-++……… 因为2 1(arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 35 211arctan (1)3521 n n x x x x x n --=-+-+-+-……… 当x=1时， -11111--(-1)4352-1 n n π=+??++? 当叠加了十万次以后得到结果π=3.141582654…只有五位有效数字，可见其精度与效率极低。如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式—— 简单公式： 11=arctan arctan 4 23 π + Machin 公式： 11=4arctan arctan 45239 π- 对以上两式进行arctan 的幂级数展开可以非常快速的求得π比较精确的数值。下面比较π精确到40位和50位数字时两个公式各需要计算多少项。 用简单公式得到40位有效数字需要叠加62项： 3.141592653589793238462643383279502884197 用简单公式得到50位有效数字需要叠加79项： 3.1415926535897932384626433832795028841971693993751 用Machin 公式得到40位有效数字需要叠加27项： 3.141592653589793238462643383279502884197 用Machin 公式得到50位有效数字需要叠加35项： 3.1415926535897932384626433832795028841971693993751 从上面简单的对比可以看出Machin 公式要优于简单公式，简单公式要优于不用公式的arctan 幂级数展开。在得到相同精度的条件下，Machin 公式所需要的叠加步数要显著少于简单公式，并且在计算精度越高的情况下，优势越明显。道理很简单，因为 Machin 公式计算的收敛速度要显著快于普通公式。简单公式决定收敛速度的是12n ?? ??? ，而Machin 公式决定收敛速度的是15n ?? ???，因为15n ?? ???的收敛速度快于12n ?? ??? ，故Machin 公式计算pi 的时候收敛速度要快于普通公式。所以Machin 公式比普通公式更加精确，并且在计算高精度的时候有更大优势。 2) 根据三角函数公式有：
PI的计算算法
上机7：文件操作 一、问题： 古人计算pi （3.1415926535 89793 23846……），一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。大约公元前1200年，中国人计算到小数点后1位；Archimedes 在公元前250年用正96边形得到pi 小数点后3位的精度；公元263年，刘徽用正3072边形计算到小数点后5位 ； 公元480年，祖冲之计算到小数点后7位；Ludolph Van Ceulen 在1609年用正262边形得到了35位 精度；1706年，John Machin 计算到小数点后100位。1874年，William Shanks 计算到小数点后707位（527位正确）。这种基于几何的算法计算量大，速度慢，吃力不讨好。 1973年Guilloud & Bouyer 用 CDC 7600计算机计算到1,001,250位；1989年Kanada & Tamurar 用 HITACHI S-820/80计算机计算到1,073,741,799位；1999年Takahashi & Kanada 用 HITACHI SR8000计算机计算到206,158,430,000位；pi 的最新计算纪录由两位日本人Daisuke Takahashi 和Yasumasa Kanada 所创造。他们在日本东京大学的IT 中心，以Gauss-Legendre 算法编写程序，利用一台每秒可执行一万亿次浮点运算的超级计算机，从日本时间1999年9月18日19:00:52起，计算了37小时21分04秒，得到了pi 的 206,158,430,208(3*236)位十进制精度，之后和他们于1999年6月27日以Borwein 四次迭代式计算 了46小时得到的结果相比，发现最后45位小数有差异，因此他们取小数点后206,158,430,000位的p 值为本次计算结果。这一结果打破了他们于1999年4月创造的68,719,470,000位的世界纪录。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算pi 的公式。 Machin 公式 1 15239164arg arctg tg π=- 3 5 7 211()(1)35721 n n x x x x a rctg x x n --=-+-++-- 这个公式由英国天文学教授John Machin 于1706年发现。他利用这个公式计算到了100位的pi 值。Machin 公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 Bailey-Borwein-Plouffe 算法 0142 1 1 ()1681848586n n n n n n π∞ ==---++++∑ 这个公式简称BBP 公式，由David Bailey, Peter Borwein 和Simon Plouffe 于1995年共同发表。它打破了传统的p 的算法，可以计算p 的任意第n 位，而不用计算前面的n-1位。这为p 的分布式计算提供了可行性。 因此使用这个公式计算，与计算常数e （2.718282……）的公式非常类似，编程很容易实现。 二、要求： 编写程序计算常数π到小数点后任意位，并将数据写入文件中。 参考程序： #include
数学计算方法实验报告
数学计算方法实验报告 习题二 2.估计用二分法求方程f(x)=x3+4x2-10=0在区间[1,2]内根的近似值，为使方程不超过10时所需的二分次数。f(x k) 程序过程： function two (tolerance) a=1;b=2;counter=0; while (abs(b-a)>tolerance) c=(a+b)/2; fa=a^3+4*a^2-10;
fb=b^3+4*b^2-10; fc=c^3+4*c^2-10; if ((fa==0|fb==0)) disp(counter); elseif (fa*fc<0) b=c;counter=counter+1; elseif (fb*fc<0) a=c;counter=counter+1; elseif (fb==0) disp(counter); end end solution=(a+b)/2; disp(solution); disp(counter); 实验结果： 6.取x0=1.5,用牛顿迭代法求第三中的方程根.f(x)=x3+4x2-10=0的近似值（精确到||x k+1-x k|≦10-5,并将迭代次数与3题比较。 程序过程： function six (g) a=1.5; fa=a^3+4*a^2-10;
ga=3*a^2+8*a; b=a-fa/ga; k=1; while(abs(b-a)>g) a=b; fa=a^3+4*a^2-10; ga=3*a^2+8*a; b=a-fa/ga; k=k+1; end format long; disp(a); disp(k); 实验结果：程序结果计算结果 8.用弦割法求方程f(x)=x3-3x2-x+9=0在区间[-2,-1]内的一个实根近似值x k,|f(x k)|≦10-5. 程序过程： function eight (t) a=-2; b=-1; fa=a^3-3*a^2-a+9; fb=b^3-3*b^2-b+9; c=b-fb*(b-a)/(fb-fa); k=1; while(abs(c-b)>t) a=b; b=c; fa=a^3-3*a^2-a+9; fb=b^3-3*b^2-b+9; c=b-fb*(b-a)/(fb-fa); k=k+1; end
《数学软件》实验报告-符号计算基础与符号微积分
实验报告 课程名称：数学软件姓名： 学院： 专业： 年级： 学号： 指导教师： 职称： 年月日
实验项目列表
附件三： 实验报告（二） 系：专业：年级：姓名学号：实验课程： 实验室号：_ 实验设备号：实验时间： 指导教师签字：成绩： 1. 实验项目名称：符号计算基础与符号微积分 2. 实验目的和要求 1.掌握定义符号对象的方法 2.掌握符号表达式的运算法则以及符号矩阵运算 3.掌握求符号函数极限及其导数的方法 4.掌握求符号函数定积分和不定积分的方法 3. 实验使用的主要仪器设备和软件 方正商祺N260微机；MATLAB7. 0或以上版本 4. 实验的基本理论和方法 （1）符号函数;sym(x)；syms a b …… （2）平方根：sqrt(x) （3）分解因式：factor（s） （4）符号表达式化简：simplify（s） （5）逆矩阵：inv(x) （6）下三角矩阵：tril（x） （7）矩阵行列式的值:det(x)
（8）符号函数求极限：limit （f ，x ，a ）；limit （f ，x ，a ，‘right ’） （9）符号函数求导：diff （f ，v ，n ） （10）符号函数求不定积分：int （f ，v ） （11）符号函数求定积分：int （f ，v ，a ，b ） 5. 实验内容与步骤 （描述实验中应该做什么事情，如何做等，实验过程中记录发生的现象、中间结果、最终得到的结果，并进行分析说明） （包括：题目，写过程、答案） 题目： 1. 已知x=6，y=5，利用符号表达式求 y x x z -++= 31。 提示：定义符号常数)'5(')'6('sym y sym x ==，。 >> x=sym('6'); >> y=sym('5'); >> z=(x+1)/(sqrt(3+x)-sqrt(y)) z = 7/(3-5^(1/2)) 2. 分解因式：44y x - >> syms x y; >> A=x^4-y^4; >> factor(A) ans = (x-y)*(x+y)*(x^2+y^2) 3. 化简表达式 (1)2121sin cos cos sin ββββ- (2) 123842+++x x x (1) >> syms x y; >> f1=sin(x)*cos(y)-cos(x)*sin(y);
数学实验：怎样计算圆周率
怎样计算 姓名： 学号 班级：数学与应用数学4班
实验报告 实验目的：自己尝试利用Mathematica软件计算的近似值，并学会计算的近似值的方法。 实验环境：Mathematica软件 实验基本理论和方法： 方法一：数值积分法（单位圆的面积是，只要计算出单位圆的面积也就计算出了的值） 其具体内容是：以单位圆的圆心为原点建立直角坐标系，则单位圆在第一象限内的部分G是一个扇形， 由曲线（）及坐标轴围成，它的面积是，算出了S的近似值，它的4倍就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图
图一 扇形G中，作平行于y轴的直线将x轴上的区间[0，1]（也就是扇形在x轴上的半径）分成n等份（n=20），相应的将扇形G分成n个同样宽度1/n的部分（）。每部分是一个曲边梯形：它的左方、右方的边界是相互平行的直线段，类似于梯形的两底；上方边界是一段曲线，因此称为曲边梯形。如果n很大，每个曲边梯形的上边界可以近似的看成直线段，从而将近似的看成一个梯形来计算它的面积；梯形的高（也就是它的宽度）h=1/n，两条底边的长分别是和，于是这个梯形面积可以作为曲边梯形面积的近似值。所有这些梯形面积的和T就可以作为扇形面积S的近似值： n越大，计算出来的梯形面积之和T就越接近扇形面积S，而4T就越接近的准确值。 方法二：泰勒级数法 其具体内容是：利用反正切函数的泰勒级数 计算。 方法三：蒙特卡罗法
其具体内容是：单位正方形的面积=1，只要能够求出扇形G 的面积S在正方形的面积中所占的比例，就能立即得到S，从而得到的值。而求扇形面积在正方形面积中所占的比例k的值，方法是在正方形中随机地投入很多点，使所投的每个点落在正方形中每一个位置的机会均等，看其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投的点的总数n的比可以作为k的近似值。能够产生在区间[0，1]内均匀分布的随机数，在Mathematica中语句是 Random[ ] 产生两个这样的随机数x，y，则以（x，y）为坐标的点就是单位正方形内的一点P，它落在正方形内每一个位置的机会均等。P落在扇形内的充分必要条件是。这样利用随机数来解决数学问题的方法叫蒙特卡罗法。 实验内容、步骤及其结果分析： 问题1：在方法一中，取n=1000，通过计算图一中扇形面积计算的的近似值。 分析：图一中的扇形面积S实际上就是定积分。 与有关的定积分很多，比如的定积分
计算pi
实验四你会用几种方法计算PI（圆周率）的值 一、问题分析 若想计算π的值，就要将跟π有关联的联系在一起，找到与π近似等价的式子，利用计算其值来得到π的值，还有对于含有π的面积、体积等关系式，可以尽量使用较规则的图形来代替进行面积、体积的求解。 二、模型建立 2.1数值积分法 找一个积分值等于π的定积分,则只要利用定积分计算出的值，就可以得到π的近似值。2.2幂级数法 利用arctanx的泰勒展开式，计算π的近似值。 当x=1时，arctan1= 2.3迭代法 1976年的迭代算法： 2.4 随机模拟法（蒙特卡罗方法） 用随机模拟求单位圆面积 向单位正方形随机投n块小石头，n很大时小石头大致均匀第分布在正方形中，如果有k块落在单位圆内，单位圆面积的近似值 三、解决问题所需的基本理论和方法 (1)对于定积分,则只要计算出的值，就可以得到π的近似值,也就是计算出与直线 y=0,x=0,x=1所围成的曲边梯形，而对于此类计算往往采用数值积分梯形公式计算。 梯形公式：将积分区间n等分将所有梯形面积加起来得到 Trapz(x)：输出数组x，输出按梯形公式x的积分(单位步长) Trapz(x,y):计算y对x的梯形积分，其中x、y定义函数关系y=f(x) (2)利用arctanx的泰勒展开式，计算π的近似值。 函数taylor用于实现Taylor级数 r=taylor(f,n,v),指定自变量v和阶数n r= taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数 (3)级数法
由于利用arctanx的幂级数展开法的收敛较慢，可采用公式 的计算来求pi值。 (4)特殊公式（BBP） 四、设计算法、编程求解 4.1数值积分法 梯形公式Matlab代码： format long x=0:0.1:1; % x=0:0.01:1; x=0:0.02:1; x=0:0.001:1; x=0:0.0001:1; y=sqrt(1-x.^2); pi=4*trapz(x,y) MAtlab结果： 数值积分法（梯形公式）Pi值 n=10 3.104518326248318 n=100 3.140417031779045 n=500 3.141487477002141 n=1000 3.141555466911028 n=10000 3.141591477611324 4.2幂级数法Matlab代码： (1) format long syms x f=atan(x); t= taylor(f,10,x,0); % t= taylor(f,100,x,0); t= taylor(f,500,x,0); t= taylor(f,1000,x,0); t= taylor(f,10000,x,0); x=1; pi=4*eval(t) 结果： 级数法（arctanx）Pi值 n=10 3.339682539682540 n=100 3.121594652591011 n=500 3.137592669589472 n=1000 3.139592655589785
如何计算π的值
1、蒙特卡罗（Monte Carlo ）法 思想： 取一正方形A ，以A 的一个顶点为圆心，A 的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B 。已知A 的面积，只要求出B 的面积与A 的面积之比B A S k S =，就能得出B S ，再由B 的面积为圆面积的四分之一，利用公式2=S R π圆即可求出π的值。因此，我 们的目的就是要找出k 的值。 可以把A 和B 看成是由无限多个点组成，而B 内的所有点都在A 内。随机产生n 个点，若落在B 内的有m 个点（假定A 的边长为1，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标x 、y 满足221x y +≤的点，就是落在B 内的点），则可近似得出k 的值，即m k n =，由此就可以求出π的值。 程序（1）： i=1;m=0;n=1000; for i=1:n a=rand(1,2); if a(1)^2+a(2)^2<=1 m=m+1; end end p=vpa(4*m/n,30) 程序运行结果：
p = 2、泰勒级数法 思想： 反正切函数arctan x 的泰勒级数展开式为： 35 211arctan (1)3521 k k x x x x x k --=-+-???+-+???- 将1x =代入上式有 1 111arctan11(1)43521 n n π-==-+-???+--. 利用这个式子就可以求出π的值了。 程序（2）： i=1;n=1000;s=0; for i=1:n s=s+(-1)^(i-1)/(2*i-1); end p=vpa(4*s,30) 程序运行结果： p = 当取n 的值为10000时，就会花费很长时间，而且精度也不是很高。原因是1x =时，arctan1的展开式收敛太慢。因此就需要找出一个x 使得arctan x 收敛更快。
离散数学实验报告()
《离散数学》实验报告 专业网络工程 班级 姓名 学号 授课教师 二 O 一六年十二月
目录 实验一联结词的运算 实验二根据矩阵的乘法求复合关系 实验三利用warshall算法求关系的传递闭包实验四图的可达矩阵实现
实验一联结词的运算 一．实验目的 通过上机实验操作，将命题连接词运算融入到C语言的程序编写中，一方面加强对命题连接词运算的理解，另一方面通过编程实现命题连接词运算，帮助学生复习和锻炼C语言知识，将理论知识与实际操作结合，让学生更加容易理解和记忆命题连接词运算。二．实验原理 (1) 非运算, 符号: ,当P=T时，P为F, 当P=F时，P为T 。 (2) 合取, 符号: ∧ , 当且仅当P和Q的真值同为真，命题P∧Q的真值才为真；否则，P∧Q的真值为假。 (3) 析取, 符号: ∨ , 当且仅当P和Q的真值同为假，命题P∨Q的真值才为假；否则，P∨Q的真值为真。 (4) 异或, 符号: ▽ , 当且仅当P和Q的真值不同时，命题P▽Q的真值才为真；否则，P▽Q的真值为真。 (5) 蕴涵, 符号: →, 当且仅当P为T,Q为F时，命题P→Q的真值才为假；否则，P→Q 的真值为真。 (6) 等价, 符号: ?, 当且仅当P,Q的真值不同时，命题P?Q的真值才为假；否则，P→Q的真值为真。 三．实验内容 编写一个程序实现非运算、合取运算、析取运算、异或运算、蕴涵运算、等价运算。四．算法程序 #include
pi的计算 数值分析论文
古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发现了许多计算圆周率的公式。下面挑选一些经典的常用公式加以介绍。除了这些经典公式外，还有很多其它公式和由这些经典公式衍生出来的公式，就不一一列举了。 1、 Machin公式 [这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 Machin.c 源程序 还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT(Fast Fourier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为O(nlog(n))。 2、 Ramanujan公式 1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。 1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为： 这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky公式的另一个更方便于计算机编程的形式是： 3、AGM(Arithmetic-Geometric Mean)算法 Gauss-Legendre公式： 这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。 4、Borwein四次迭代式： 这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆周率。
pi的无穷级数算法
3π的计算模型—π的无穷级数算法
思想：利用一些特殊函数的幂级数展开式计算圆周率。 根据幂级数展开的相关知识，易知： 241222 1 1(1)1n n x x x x --=-+-+-++……… 因为2 1 (arctan )'1x x =+，故可以求得arctan x 的幂级数展开式为： 3521 1arctan (1)3521n n x x x x x n --=-+-+-+-……… 当x=1时， -11111--(-1)4352-1 n n π =+??++? 用Matlab 计算： 创建m 文件calpi1.m ，内容如下： function y=calpi1(k) for n=1:k a(n)=(-1).^(n-1)./(2*n-1); end ; 4*sum(a) 在命令窗口中输入如下命令：
它是一个与π有关的无穷级数，实际计算时，我们只能使用有限项。如果取级数前n 项之和作为π的近似值，其误差为 ｜｜≤ 为了保证误差不超过 ，就要取级数的前20000项进行计算，计算 量之大可想而知。现在看来，计算π的级数有明显的缺点：级数收敛太慢，计算量过大。其原因是｜x ｜偏大。如果想要精确计算π的数值的话，非常有必要寻找改进以后的方法，这就引出了两个能够提高计算效率的公式—— （1）欧拉公式：11=arctan arctan 4 23 π + 现取，令α=，显见0＜α＜， 记β=，而 = =， 所以β= ，就是11=arctan arctan 4 2 3 π + （2）马庭公式：11=4arctan arctan 4 5239 π - 可取，令，则 ， ，故4α≈，再令 ， 即，而，就有11=4arctan arctan 4 5 239 π - 将欧拉公式和马庭公式与arctanx 的泰勒级数相结合，会加快该级数的收敛速度，具有很强的实用性。
计算方法实验报告 插值
实验名称：插值计算 1引言 在生产和科研中出现的函数是多种多样的。常常会遇到这样的情况：在某个实际问题中，虽然可以断定所考虑的函数f（x）在区间[a,b]上存在且连续，但却难以找到它的解析表达式，只能通过实验和观测得到在有限个点上的函数值。用这张函数表来直接求出其他点的函数值是非常困难的，在有些情况下，虽然可以写出f(x)的解析表达式，但由于结构十分复杂，使用起来很不方便。面对这些情况，构造函数P(x)作为f(x)的近似，插值法是解决此类问题比较古老却目前常用的方法，不仅直接广泛地应用与生产实际和科学研究中，而且是进一步学习数值计算方法的基础。 设函数y=f(x)在区间[a,b]上连续，且在n+1个不同的点a≤x0,x1……,xn≤b上分别取值y0,y1……,yn. 插值的目的就是要在一个性质优良、便于计算的函数φ中，求一简单函数P(x),使P(xi)=yi(i=0,1…,n)而在其他点x≠xi上，作为f(x)的近似。 通常，称区间[a,b]为插值区间，称点x0,x1,…,xn为插值节点，上式为插值条件，称函数类φ为插值函数类，称P(x)为函数f(x)在节点x0,x1,…,xn处的插值函数，求插值函数P(x)的方法称为插值法。 2实验目的和要求 用matlab定义分段线性插值函数、分段二次插值函数、拉格朗日插值函数，输入所给函 数表，并利用计算机选择在插值计算中所需的节点，计算f(0.15),f(0.31),f(0.47)的近似值。
3算法描述 1.分段线性插值流程图
2.分段二次插值流程图
3.拉格朗日插值流程图
4程序代码及注释 1.分段线性插值
圆周率的计算数学实验报告
数学实验报告（二） 一、实验题目：圆周率的计算 二、实验目的： 1.用多种方法计算圆周率的值； 2.通过实验来说明各种方法的优劣； 3.尝试提出新的计算方法。 三、实验内容和方法： 1.古典方法： 用圆内接正多边形和圆外切正多边形来逼近 以阿基米德的圆内接96边形和圆外切96边形逼近为例 已知：sin < y=4*x 得出当k=10时，π≈3.232315809405593 编写程序 syms k x=symsum((-1)^k/(2*k+1),k,0,20) y=4*x 得出当k=20时，π≈3.189184782277595 依次，加大k 的值 K=50，π≈3.161198612987050 K=100，π≈3.151493401070990 K=200，π≈3.146567747182986e+159 … (2).沃里斯(Wallis)方法 ∏∞=??? ??+?-?=??? ??????? ??? ???? ????=1122122276565 43432122k k k k k π 编写程序： format long x=1; for k=1:10 x=x*(2*k/(2*k-1)*2*k/(2*k+1)); end y=2*x 得k=10时，π≈3.067703806643498 增加k 的值 K=20，π≈3.103516961539230 K=50，π≈3.126078********* K=100，π≈3.133787********* K=10000，π≈3.141514118681864 K=1000000，π≈3.141591868191880 … (3).利用公式3 1arctan 21arctan 1arctan 4+==π 推出π=4(31 arctan 21 arctan +) 编写程序： syms n; f1=(-1)^(n-1)*(1/2)^(2*n-1)/(2*n-1); f2=(-1)^(n-1)*(1/3)^(2*n-1)/(2*n-1); ans1=symsum(f1,n,1,79); Pi值计算(积分法) 程序清单： #include 差值法实验日志 实验题目：插值法 实验目的： 1．掌握拉格朗日插值、牛顿插值、分段低次插值和样条插值的方法。 2．对四种插值结果进行初步分析。 实验要求： （1）写出算法设计思想； （2）程序清单； （3）运行的结果； （4）所得图形； （5）四种插值的比较； （6）对运行情况所作的分析以及本次调试程序所取的经验。如果程序未通过，应分析其原因。 实验主要步骤： 1．已知函数) f满足： (x x0.0 0.1 0.195 0.3 0.401 0.5 f(0.39894 0.39695 0.39142 0.38138 0.36812 x ） 0.35206 （1）用分段线性插值； 打开MATLAB，按以下程序输入： x0=-5:5; y0=1./(1+x0.^2); x=-5:0.1:5; y=1./(1+x.^2); y1=lagr(x0,y0,x); y2=interp1(x0,y0,x); y3=spline(x0,y0,x); for k=1:11 xx(k)=x(46+5*k); yy(k)=y(46+5*k); yy1(k)=y1(46+5*k); yy2(k)=y2(46+5*k); yy3(k)=y3(46+5*k); end [xx;yy;yy2;yy3]' z=0*x; plot(x,z,x,y,'k--',x,y2,'r') plot(x,z,x,y,'k--',x,y1,'r') pause plot(x,z,x,y,'k--',x,y3,'r') 回车得以下图形: （2） 拉格朗日插值。 创建M 文件，建立lagr 函数： function y=lagr1(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 新建一个M 文件，输入： x0=[0.0 0.1 0.195 0.3 0.401 0.5]; y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; x=0.0:0.01:0.5; y1=lagr1(x0,y0,x); 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 怎样计算 姓名: 学号 班级:数学与应用数学4班 实验报告 实验目的:自己尝试利用Mathematica软件计算的近似值,并学会计算的近似值的方法。 实验环境:Mathematica软件 实验基本理论与方法: 方法一:数值积分法(单位圆的面积就是,只要计算出单位圆的面积也就计算出了的值) 其具体内容就是:以单位圆的圆心为原点建立直角坐标系,则单位圆在第一象限内的部分G就是一个扇 形, 由曲线()及坐标轴围成,它的面积就是,算出了S的近似值,它的4倍就就是的近似值。而怎样计算扇形G的面积S的近似值呢?如图 图一 扇形G中,作平行于y轴的直线将x轴上的区间[0,1](也就就是扇形在x轴上的半径)分成n等份(n=20),相应的将扇形G分成n个同样宽度1/n的部分()。每部分就是一个曲边梯形:它的左方、右方的边界就是相互平行的直线段,类似于梯形的两底;上方边界就是一段曲线,因此称为曲边梯形。如果n很大,每个曲边梯形的上边界可以近似的瞧成直线段,从而将近似的瞧成一个梯形来计算它的面积;梯形的高(也就就是它的宽度)h=1/n,两条底边的长分别就 是与,于就是这个梯形面积 可以作为曲边梯形面积的近似值。所有这些梯形面积的与T就可以作为扇形面积S的近似值: n越大,计算出来的梯形面积之与T就越接近扇形面积S,而4T就越接近的准确值。 方法二:泰勒级数法 其具体内容就是:利用反正切函数的泰勒级数 计算。 方法三:蒙特卡罗法 其具体内容就是:单位正方形的面积=1,只要能够求出扇形G 的面积S在正方形的面积中所占的比例,就能立即得到S,从而得到的值。而求扇形面积在正方形面积中所占的比例k的值,方法就是在正方形中随机地投入很多点,使所投的每个点落在正方形中每一个位置的机会均等,瞧其中有多少个点落在扇形内。将落在扇形内的点的个数m与所投的点的总数n的比可以作为k 的近似值。能够产生在区间[0,1]内均匀分布的随机数,在Mathematica 中语句就是 Random[ ] 产生两个这样的随机数x,y,则以(x,y)为坐标的点就就是单位正方形内的一点P,它落在正方形内每一个位置的机会均等。P落在扇形内的充分必要条件就是。这样利用随机数来解决数学问题的方法叫蒙特卡罗法。 实验内容、步骤及其结果分析: 问题1:在方法一中,取n=1000,通过计算图一中扇形面积计算的 的近似值。积分法计算pi值
数值计算方法实验报告
数学实验怎样计算圆周率