统计计算统计模拟的估计.docVIP

统计模拟作业 班级 学号 姓名 作业一：π的估计 1. 分析方法 从十七世纪中叶起，人们开始用更先进的分析方法来求∏的近似值，其中应用的主要工具是收敛的无穷乘积和无穷级数，下面列举一些用此类方法求∏近似值的实例。 (1)由公式 推出=4 编写程序： syms k x=symsum((-1)^k/(2*k+1),k,0,10) y=4*x 得出当k=10时，π≈3.5593 编写程序 syms k x=symsum((-1)^k/(2*k+1),k,0,20) y=4*x 得出当k=20时，π≈3.7595 依次，加大k的值 K=50，π≈3.7050 K=100，π≈3.0990 K=200，π≈3.2986 … K π 10 3.5593 50 3.7050 100 3.0990 200 3.2986 （2）沃里斯(Wallis)方法 在积分学中我们经常会遇到如下的沃利斯(Wallis)公式。沃利斯(Wallis)公式揭示了π与整数之间的一种很不寻常的关系。但在实际学习中很少注意到沃利斯(Wallis)公式,更不会关注它的应用。实际上,沃利斯(Wallis)公式有许多作用,经常有以下几方面的应用。1应用于极限计算中。由于沃利斯(Wallis)公式与极限有关,所以有些极限的计算可以通过沃利斯(Wallis)公式很容易计算出来。2.应用于积分计算中。对于一些用积分法不易求得原函数的积分，而用沃利斯(Wallis)公式却很容易解决问题。 编写程序： 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.3498 增加k的值 K=20，π≈3.9230 K=50，π≈3.5409 K=100，π≈3.8159 K=10000，π≈3.1864 K=，π≈3.1880 … k π 10 3.3498 20 3.9230 50 3.5409 100 3.8159 10000 3.1864 3.1880 （3）蒙特卡罗法 取一正方形A，以A的一个顶点为圆心，A的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B。已知A的面积，只要求出B的面积与A 的面积之比，就能得出，再由B的面积为圆面积的四分之一，利用公式即可求出的值。因此，我们的目的就是要找出的值。 可以把A和B看成是由无限多个点组成，而B内的所有点都在A内。随机产生个点，若落在B内的有个点（假定A的边长为1，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标、满足的点，就是落在B内的点），则可近似得出的值，即，由此就可以求出的值。 编写程序： 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 = 3. 泰勒级数法 函数的泰勒级数展开式为： 将代入上式有 . 利用这个式子就可以求出的值了。 编写程序： 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 = 3. 当取的值为10000时，就会花费很长时间，而且精度也不是很高。原因是时，的展开式收敛太慢。因此就需要找出一个使得收敛更快。 若取，则我们只有找出与的关系，才能求出的值。 令，，根据公式有，则有。所以可以用来计算的值。 （5）利用公式 推出π=4() 编写程序： i=1;n=1000;s=0;s1=0;s2=0; for i=1:n s1=s1+(-1)^(i-1)*(1/2)^(2*i-1)/(2*i-1); s2=s2+(-1)^(i-1)*(1/3)^(2*i-1)/(2*i-1); end s=s1+s2; p=vpa(4*s,30) 程序运行结果： p = 3. 显然，级数收敛越快，取同样的值可以得到更高的精度。以同样的方法，能得出，程序和上面的一样。这样的近似值可以精确到几百位。 麦琴公式 由 推出π=4()，这个公式由英国天文学教授约翰·麦琴于1706年发现。他利用这个公式计算到了100位的圆周率。麦琴公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。 编写程序： syms n; f1=(-1)^(n-1)*(1/5)^(2*n-1)/(2*n-1); f2=(-1)^(n-1)*(1/239)^