《数学软件作业圆周率的计算》由会员分享，可在线阅读，更多相关《数学软件作业圆周率的计算（6页珍藏版）》请在金锄头文库上搜索。
1、2007 数学应用软件设计内容怎样计算圆周率 的值实验报告姓名：。 。 学号：；实验目的：利用所学的数学应用软件，解决下列问题： 1. 用反正切函数的幂级数展开式结合有关公式求 ，若要精确到以 40 位、50 位数字，试比较简单公式和 Machin 公式所用的项数；2. 用数值积分计算 ，分别用梯形法和 Simpson 法精确到 10 位数字，用 Simpson 法精确到 15 位数字；3. 用 Monte Carlo 法计算 ，除了加大随机数，在随机数一定时可重复算若干次后求平均值，看能否求得 5 位精确数字？4. 设计方案用计算机模拟 Buffon 实验；5.利用学习过的知识(或查阅资料
2、)，提出其他计算 的方法。实验的基本理论与方法：1.利用反正切函数的幂级数表达式：利用积分公式两边积分得：在上式中令 x=1 得：用 Mathematica 计算In1 k=1000;S1=N4*Sum(-1)(n-1)/(2n-1),n,1,k,18Out2 3.14059265383979293In3 k=10000;Out4 3.14149265359004324In5 k=15000;Out6 3.14152598692320065In7 k=20000Out8 3.14154265358982449LL21422)(1nxxx LL12)(53arctnnxxx12)(5134n由上
3、述计算可以知道，用反正切函数的幂级数展开来计算 的值，运算速度比较慢，结果精度也不高,因此此方法一般不选取。下面对精确到相同精度时简单公式和 Machin 公式所用的项数进行比较：简单公式： Machin 公式： 再用 Mathematica 计算：当精确到 40 位时：对于简单公式：Cleark,n,SIn9 k =62;s2 = N4*Sum(-1)(n - 1)*(1/2)(2n - 1)/(2n - 1) + (-1)(n - 1)*(1/3)(2n - 1)/(2n - 1), n, 1, k, 40Out10 3.1415926535897932384626433832795028
4、84197对于 Machin 公式 ：In11 k = 28;s3 = N4*Sum4*(-1)(n - 1)*(1/5)(2n - 1)/(2n - 1) - (-1)(n - 1)*(1/239)(2n - 1)/(2n - 1), n, 1, k, 40Out12 3.141592653589793238462643383279502884197当精确到 50 位时：对于简单公式：In13 k = 78;s2 = N4*Sum(-1)(n - 1)*(1/2)(2n - 1)/(2n - 1) + (-1)(n - 1)*(1/3)(2n - 1)/(2n - 1), n, 1, k,
5、40Out14 3.1415926535897932384626433832795028841971693993750对于 Machin 公式 ：In15 k = 35;s3 = N4*Sum4*(-1)(n - 1)*(1/5)(2n - 1)/(2n - 1) - (-1)(n - 1)*(1/239)(2n - 1)/(2n - 1), n, 1, k, 50Out16 3.1415926535897932384626433832795028841971693993751431arctn2rtLL1253 )()(124 n )31(3532Ln4291arctn5rt4由上述比较可以看
6、出：在精确到相同精度时，用 Machin 公式计算的结果比较准确，并且所需的项数也比较少，当精确到 40 位时简单公式和 Machin 公式分别需要 62 和 28 位，而当精确到 50 位时，分别为 78 和 35 位。所以用 Machin 公式计算圆周率速度较快，结果也较精确，并且程序也较简单，是一种比较好的方法。2. 用数值积分计算 ：将0,1区间 n 等分，取 xk=k/n, yk= 1/ (1+xk2)用 Mathematica 计算：梯形法（精确到 10 位）：In17 yx_:=4/(1+x2);n=20000;S4=N1/(2*n)*(Sum2*yk/n,k,1,n-1+y0+
7、y1),20Out18 3.1415926531731265718Simpson 法（精确到 10 位）：In19 yx_ := 4/(1 + x2); m = 11;S5= N1/(3*2 m)*(Sum2*yk/(2 m), k, 1, 2m - 1 + Sum2*yk/(2 m), k, 1, 2m - 1, 2 + y0 + y1), 20Out20 3.1415926532398091222Simpson 法（精确到 15 位）：In21 yx_ := 4/(1 + x2); m = 80;S6 = N1/(3*2 m)*(Sum2*yk/(2 m), k, 1, 2m - 1 +
8、Sum2*yk/(2 m), k, 1, 2m - 1, 2 + y0 + y1), 20Out22 3.1415926535897908732梯形法精确到 10 位需要 20000 次的迭代，而 Simpson 法只需要 11 次就可以达到，并且精确到 15 位时也仅需 80 次计算，这也体现了 Simpson 法的优越性运算速度快，精dxA1024 21)(xy设 )(2011nnyAL梯 形 法 )()(324220 mmyySimpson法 )(131yL度高。3. 用 Monte Carlo 法计算 ：这个实验方法的操作很简单：找一根粗细均匀，长度为 d 的细针，并在一张白纸上画上一
9、组间距为 l 的平行线（方便起见，常取 l = d/2），然后一次又一次地将小针任意投掷在白纸上。这样反复地投多次，数数针与任意平行线相交的次数，于是就可以得到 的近似值。因为有科学家证明了针与任意平行线相交的概率为 p = 2l/d 。利用这一公式，可以用概率方法得到圆周率的近似值。曾经有人试验并得出一结果，选取 l = d/2 ，然后投针 2212 次，其中针与平行线相交 704 次，这样求得圆周率的近似值为 2212/704 = 3.142。当实验中投的次数相当多时，就可以得到 的更精确的值。调查资料显示，1850 年，一位叫沃尔夫的人在投掷 5000 多次后，得到 的近似值为 3.15
10、96。目前宣称用这种方法得到最好结果的是意大利人拉兹瑞尼，在 1901 年，他重复这项实验，作了 3408 次投针，求得 的近似值为 3.1415929，这个结果是如此准确。由此我们可以看出在随机数一定时重复算若干次后求平均值，是可以求到得 5 位的精确数字。4. 计算机模拟 Buffon 实验：设计方案： 在正方形 0 0, i-, m = m + IfRandom2 + Random2 = 1, 1, 0; N4*m/n, 10Out24 3.2016In25 n=50000;Out26 3.1384In27 n=100000;Out28 3.15516但是这种方法受外界的影响比较大，存在
11、许多偶然的因素，误差范围也比较大，所以结果有时候也很不精确，但就目前前人做过的实验结果来看，也有精确到了 7 位数字的。5其他的计算方法：（1）前人求解圆周率的方法：阿基米德求圆周率的更精确近似值的方法，体现在他的一篇论文圆的测定这一书中，阿基米德第一次创用上、下界来确定 的近似值，他用几何方法证明了“圆周长与圆直径之比小于 3+(1/7) 而大于 3 + (10/71) ”，即：圆周率 属于3.140853.14286 的范围内，他还提供了误差的估计。重要的是，这种方法从理论上而言，能够求得圆周率的更准确的值。AGM(Arithmetic-Geometric Mean) 算法 Gauss-L
12、egendre 公式：初值：重复计算： 最后计算： 这个公式每迭代一次将得到双倍的十进制精度，比如要计算 100 万位，迭代 20 次就够了。1999 年 9 月 Takahashi 和 Kanada 用这个算法计算到了圆周率的 206,158,430,000 位，创出新的世界纪录。此外还有：Ramanujan 公式 1914 年，印度数学家 Srinivasa Ramanujan 在他的论文里发表了一系列共 14 条圆周率的计算公式，这是其中之一。这个公式每计算一项可以得到 8 位的十进制精度。1985 年 Gosper 用这个公式计算到了圆周率的 17,500,000 位。1989 年，
13、David & Gregory Chudnovsky 兄弟将 Ramanujan 公式改良成为：这个公式被称为 Chudnovsky 公式，每计算一项可以得到 15 位的十进制精度。1994年 Chudnovsky 兄弟利用这个公式计算到了 4,044,000,000 位。Chudnovsky 公式的另一个更方便于计算机编程的形式是：（2）针对 1593 年，韦达给出的这一公式，若取前 10 项计算，则有：2/1.414213562*1.847759065*1.961570561*1.990369453*1.9975909121:GTRn*1.999397637*1.999849404*1.9
14、99962351*1.999990588*1.999997647/1024iPc/651.8989024/1024 0.636620022Pw8 2/0.6366200223.14159 lQ4由此可见，收敛速度还是挺快的，这也不失为一种好方法。（3）推了一个公式 nsin(180/n)= (当 n)设把圆无限等分，那么如图 AB=2r/n(r 为半径) OB=OA=rOBA 为等腰三角形作 OGABAG=GB=AB/2=（2r/n）/2=r/nSinAOG=AG/AOAOG=(AOB)/2=(360/n)/2=180/nSin(180/n)=AG/AO又AO=r Sin(180/n)=AG/rAG=r/nSin(180/n)=(r/n)/rSin(180/n)=/n= Sin(180/n)* n实验小结：根据题目的要求，一一作出了解答，但是有的方法不太理想 ，还需要进一步的改进。计算圆周率的方法有很多，本文采用了很多前人计算圆周率的方法，对它们进行了综合比较，得出了几种比较精确的计算方法。
